# Designing an Optimal Kilonova Search using DECam for Gravitational Wave Events

Clecio R. Bom, J. Annis, A. Garcia, A. Palmese, N. Sherman, M.

Soares-Santos, L. Santana-Silva, R. Morgan, K. Bechtol, T. Davis, et al.

► To cite this version:

Clecio R. Bom, J. Annis, A. Garcia, A. Palmese, N. Sherman, et al.. Designing an Optimal Kilonova Search using DECam for Gravitational Wave Events. Astrophys.J., 2024, 960 (2), pp.122. 10.3847/1538-4357/ad0462. hal-04006449

HAL Id: hal-04006449

https://hal.science/hal-04006449

Submitted on 6 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Designing an Optimal Kilonova Search using DECam for Gravitational Wave Events

Clecio R. Bom, J. Annis, A. Garcia, A. Palmese, N. Sherman, M.

Soares-Santos, L. Santana-Silva, R. Morgan, K. Bechtol, T. Davis, et al.

► To cite this version:

Clecio R. Bom, J. Annis, A. Garcia, A. Palmese, N. Sherman, et al. Designing an Optimal Kilonova Search using DECam for Gravitational Wave Events. The Astrophysical Journal, 2024, 960 (2), pp.122. 10.3847/1538-4357/ad0462. hal-04006449

HAL Id: hal-04006449

https://hal.science/hal-04006449v1

Submitted on 6 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Designing an Optimal Kilonova Search Using DECam for Gravitational-wave Events

C. R. Bom<sup>1,2</sup>, J. Annis<sup>3</sup>, A. Garcia<sup>4</sup>, A. Palmese<sup>5,56</sup>, N. Sherman<sup>4</sup>, M. Soares-Santos<sup>4</sup>, L. Santana-Silva<sup>6</sup>, R. Morgan<sup>7,55</sup>, K. Bechtol<sup>7</sup>, T. Davis<sup>8</sup>, H. T. Diehl<sup>3</sup>, S. S. Allam<sup>3,57</sup>, T. G. Bachmann<sup>9</sup>, B. M. O. Fraga<sup>1</sup>, J. García-Bellido<sup>10</sup>, M. S. S. Gill<sup>11</sup>, K. Herner<sup>3</sup>, C. D. Kilpatrick<sup>12</sup>, M. Makler<sup>1,13</sup>, F. Olivares E.<sup>14</sup>, M. E. S. Pereira<sup>15</sup>, J. Pineda<sup>16,17</sup>, A. Santos<sup>1</sup>, D. L. Tucker<sup>3</sup>, M. P. Wiesner<sup>18</sup>, M. Aguena<sup>19</sup>, O. Alves<sup>4</sup>, D. Bacon<sup>20</sup>, P. H. Bernardinelli<sup>21</sup>, E. Bertin<sup>22,23</sup>, S. Bocquet<sup>24</sup>, D. Brooks<sup>25</sup>, M. Carrasco Kind<sup>26,27</sup>, J. Carretero<sup>28</sup>, C. Conselice<sup>29,30</sup>, M. Costanzi<sup>31,32,33</sup>, L. N. da Costa<sup>19</sup>, J. De Vicente<sup>34</sup>, S. Desai<sup>35</sup>, P. Doel<sup>25</sup>, S. Everett<sup>36</sup>, I. Ferrero<sup>37</sup>, J. Frieman<sup>3,38</sup>, M. Gatti<sup>39</sup>, D. W. Gerdes<sup>4,40</sup>, D. Gruen<sup>24</sup>, R. A. Gruendl<sup>26,27</sup>, G. Gutierrez<sup>3</sup>, S. R. Hinton<sup>41</sup>, D. L. Hollowood<sup>42</sup>, K. Honscheid<sup>43,44</sup>, D. J. James<sup>45</sup>, K. Kuehn<sup>46,47</sup>, N. Kuropatkin<sup>3</sup>, P. Melchior<sup>48</sup>, J. Mena-Fernández<sup>34</sup>, F. Menanteau<sup>26,27</sup>, A. Pieres<sup>19,49</sup>, A. A. Plazas Malagón<sup>48</sup>, M. Raveri<sup>50</sup>, M. Rodriguez-Monroy<sup>34</sup>, E. Sanchez<sup>34</sup>, B. Santiago<sup>19,51</sup>, I. Sevilla-Noarbe<sup>34</sup>, M. Smith<sup>52</sup>, E. Suchyta<sup>53</sup>, M. E. C. Swanson<sup>25,45</sup>, G. Tarle<sup>4</sup>, C. To<sup>43</sup>, and N. Weaverdyck<sup>4,54</sup>

<footnote n="1"><sup>1</sup>Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, CEP 22290-180, Rio de Janeiro, RJ, Brazil; debom@cbpf.br</footnote>

<footnote n="2"><sup>2</sup>Centro Federal de Educação Tecnológica Celso Suckow da Fonseca, Rodovia Mário Covas, lote J2, quadra J, CEP 23810-000, Itaguaí, RJ, Brazil</footnote>

<footnote n="3"><sup>3</sup>Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA</footnote>

<footnote n="4"><sup>4</sup>Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="5"><sup>5</sup>Department of Physics, University of California Berkeley, 366 LeConte Hall, MC 7300, Berkeley, CA 94720, USA</footnote>

<footnote n="6"><sup>6</sup>NAT-Universidade Cruzeiro do Sul/Universidade Cidade de São Paulo, Rua Galvão Bueno, 868, 01506-000 São Paulo, SP, Brazil</footnote>

<footnote n="7"><sup>7</sup>Physics Department, 2320 Chamberlin Hall, University of Wisconsin–Madison, 1150 University Avenue, Madison, WI 53706-1390, USA</footnote>

<footnote n="8"><sup>8</sup>School of Mathematics and Physics, The University of Queensland, Brisbane QLD 4072, Australia</footnote>

<footnote n="9"><sup>9</sup>Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="10"><sup>10</sup>Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="11"><sup>11</sup>SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA</footnote>

<footnote n="12"><sup>12</sup>Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA</footnote>

<footnote n="13"><sup>13</sup>International Center for Advanced Studies & Instituto de Ciencias Físicas, ECYT-UNSAM & CONICET, 1650 Buenos Aires, Argentina</footnote>

<footnote n="14"><sup>14</sup>Instituto de Astronomía y Ciencias Planetarias, Universidad de Atacama, Av. Copayapu 485, Copiapó, Chile</footnote>

<footnote n="15"><sup>15</sup>Hamburger Sternwarte, Universität Hamburg, Gothenburgsweg 112, D-21029 Hamburg, Germany</footnote>

<footnote n="16"><sup>16</sup>Departamento de Ciencias Físicas, Universidad Andres Bello, Avda. Republica 252, Santiago, Chile</footnote>

<footnote n="17"><sup>17</sup>Millennium Institute of Astrophysics (MAS), Nuncio Monseñor Sótero Sanz 100, Providencia, Santiago, Chile</footnote>

<footnote n="18"><sup>18</sup>Benedictine University, Department of Physics, 5700 College Road, Lisle, IL 60532, USA</footnote>

<footnote n="19"><sup>19</sup>Laboratório Interinstitucional de e-Astronomia—LIneA, Rua Gal. José Cristino 77, Rio de Janeiro RJ—20921-400, Brazil</footnote>

<footnote n="20"><sup>20</sup>Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK</footnote>

<footnote n="21"><sup>21</sup>Astronomy Department, University of Washington, Box 351580, Seattle, WA 98195, USA</footnote>

<footnote n="22"><sup>22</sup>CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014 Paris, France</footnote>

<footnote n="23"><sup>23</sup>Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014 Paris, France</footnote>

<footnote n="24"><sup>24</sup>University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, D-81679 Munich, Germany</footnote>

<footnote n="25"><sup>25</sup>Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="26"><sup>26</sup>Center for Astrophysical Surveys, National Center for Supercomputing Applications, 1205 West Clark Street, Urbana, IL 61801, USA</footnote>

<footnote n="27"><sup>27</sup>Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 West Green Street, Urbana, IL 61801, USA</footnote>

<footnote n="28"><sup>28</sup>Instituto de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona), Spain</footnote>

<footnote n="29"><sup>29</sup>Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK</footnote>

<footnote n="30"><sup>30</sup>University of Nottingham, School of Physics and Astronomy, Nottingham NG7 2RD, UK</footnote>

<footnote n="31"><sup>31</sup>Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo 11, I-34131 Trieste, Italy</footnote>

<footnote n="32"><sup>32</sup>INAF-Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="33"><sup>33</sup>Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy</footnote>

<footnote n="34"><sup>34</sup>Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain</footnote>

<footnote n="35"><sup>35</sup>Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India</footnote>

<footnote n="36"><sup>36</sup>Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="37"><sup>37</sup>Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029, Blindern, NO-0315 Oslo, Norway</footnote>

<footnote n="38"><sup>38</sup>Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="39"><sup>39</sup>Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA</footnote>

<footnote n="40"><sup>40</sup>Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="41"><sup>41</sup>School of Mathematics and Physics, University of Queensland, Brisbane QLD 4072, Australia</footnote>

<footnote n="42"><sup>42</sup>Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA</footnote>

<footnote n="43"><sup>43</sup>Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="44"><sup>44</sup>Department of Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="45"><sup>45</sup>Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="46"><sup>46</sup>Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia</footnote>

<footnote n="47"><sup>47</sup>Lowell Observatory, 1400 Mars Hill Road, Flagstaff, AZ 86001, USA</footnote>

<footnote n="48"><sup>48</sup>Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA</footnote>

<footnote n="49"><sup>49</sup>Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro RJ—20921-400, Brazil</footnote>

<footnote n="50"><sup>50</sup>Department of Physics, University of Genova and INFN, Via Dodecaneso 33, I-16146 Genova, Italy</footnote>

<footnote n="51"><sup>51</sup>Instituto de Física, UFRGS, Caixa Postal 15051, Porto Alegre RS—91501-970, Brazil</footnote>

<footnote n="52"><sup>52</sup>School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, UK</footnote>
<figure><img src="image_2.png" /></figure>

# Designing an Optimal Kilonova Search Using DECam for Gravitational-wave Events

C. R. Bom¹,² ID, J. Annis³ ID, A. Garcia⁴ ID, A. Palmese⁵,⁵⁶ ID, N. Sherman⁴ ID, M. Soares-Santos⁴ ID, L. Santana-Silva⁶ ID, R. Morgan⁷,⁵⁵ ID, K. Bechtol⁷ ID, T. Davis⁸ ID, H. T. Diehl³ ID, S. S. Allam³,⁵⁷ ID, T. G. Bachmann⁹ ID, B. M. O. Fraga¹, J. García-Bellido¹⁰ ID, M. S. S. Gill¹¹ ID, K. Herner³ ID, C. D. Kilpatrick¹², M. Makler¹,¹³ ID, F. Olivares E.¹⁴ ID, M. E. S. Pereira¹⁵ ID, J. Pineda¹⁶,¹⁷, A. Santos¹, D. L. Tucker³ ID, M. P. Wiesner¹⁸ ID, M. Aguena¹⁹ ID, O. Alves⁴ ID, D. Bacon²⁰ ID, P. H. Bernardinelli²¹ ID, E. Bertin²²,²³, S. Bocquet²⁴ ID, D. Brooks²⁵ ID, M. Carrasco Kind²⁶,²⁷ ID, J. Carretero²⁸ ID, C. Conselice²⁹,³⁰ ID, M. Costanzi³¹,³²,³³, L. N. da Costa¹⁹, J. De Vicente³⁴ ID, S. Desai³⁵ ID, P. Doel²⁵, S. Everett³⁶ ID, I. Ferrero³⁷ ID, J. Frieman³,³⁸ ID, M. Gatti³⁹, D. W. Gerdes⁴,⁴⁰ ID, D. Gruen²⁴ ID, R. A. Gruendl²⁶,²⁷ ID, G. Gutierrez³ ID, S. R. Hinton⁴¹ ID, D. L. Hollowood⁴² ID, K. Honscheid⁴³,⁴⁴ ID, D. J. James⁴⁵ ID, K. Kuehn⁴⁶,⁴⁷ ID, N. Kuropatkin³ ID, P. Melchior⁴⁸ ID, J. Mena-Fernández³⁴ ID, F. Menanteau²⁶,²⁷ ID, A. Pieres¹⁹,⁴⁹ ID, A. A. Plazas Malagón⁴⁸ ID, M. Raveri⁵⁰, M. Rodriguez-Monroy³⁴, E. Sanchez³⁴ ID, B. Santiago¹⁹,⁵¹, I. Sevilla-Noarbe³⁴ ID, M. Smith⁵² ID, E. Suchyta⁵³ ID, M. E. C. Swanson²⁵,⁴⁵, G. Tarle⁴ ID, C. To⁴³ ID, and N. Weaverdyck⁴,⁵⁴ ID

¹Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, CEP 22290-180, Rio de Janeiro, RJ, Brazil; debom@cbpf.br

²Centro Federal de Educação Tecnológica Celso Suckow da Fonseca, Rodovia Mário Covas, lote J2, quadra J, CEP 23810-000, Itaguaí, RJ, Brazil

³Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA

⁴Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

⁵Department of Physics, University of California Berkeley, 366 LeConte Hall, MC 7300, Berkeley, CA 94720, USA

⁶NAT-Universidade Cruzeiro do Sul/Universidade Cidade de São Paulo, Rua Galvão Bueno, 868, 01506-000 São Paulo, SP, Brazil

⁷Physics Department, 2320 Chamberlin Hall, University of Wisconsin-Madison, 1150 University Avenue, Madison, WI 53706-1390, USA

⁸School of Mathematics and Physics, The University of Queensland, Brisbane QLD 4072, Australia

⁹Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA

¹⁰Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain

¹¹SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

¹²Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA

¹³International Center for Advanced Studies & Instituto de Ciencias Físicas, ECYT-UNSAM & CONICET, 1650 Buenos Aires, Argentina

¹⁴Instituto de Astronomía y Ciencias Planetarias, Universidad de Atacama, Av. Copayapu 485, Copiapó, Chile

¹⁵Hamburger Sternwarte, Universität Hamburg, Gothenburgsweg 112, D-21029 Hamburg, Germany

¹⁶Departamento de Ciencias Físicas, Universidad Andres Bello, Avda. Republica 252, Santiago, Chile

¹⁷Millennium Institute of Astrophysics (MAS), Nuncio Monseñor Sótero Sanz 100, Providencia, Santiago, Chile

¹⁸Benedictine University, Department of Physics, 5700 College Road, Lisle, IL 60532, USA

¹⁹Laboratório Interinstitucional de e-Astronomia—LIneA, Rua Gal. José Cristino 77, Rio de Janeiro RJ—20921-400, Brazil

²⁰Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK

²¹Astronomy Department, University of Washington, Box 351580, Seattle, WA 98195, USA

²²CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014 Paris, France

²³Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014 Paris, France

²⁴University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, D-81679 Munich, Germany

²⁵Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK

²⁶Center for Astrophysical Surveys, National Center for Supercomputing Applications, 1205 West Clark Street, Urbana, IL 61801, USA

²⁷Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 West Green Street, Urbana, IL 61801, USA

²⁸Instituto de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona), Spain

²⁹Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK

³⁰University of Nottingham, School of Physics and Astronomy, Nottingham NG7 2RD, UK

³¹Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo 11, I-34131 Trieste, Italy

³²INAF-Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, I-34143 Trieste, Italy

³³Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy

³⁴Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain

³⁵Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India

³⁶Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

³⁷Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029, Blindern, NO-0315 Oslo, Norway

³⁸Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

³⁹Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

⁴⁰Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

⁴¹School of Mathematics and Physics, University of Queensland, Brisbane QLD 4072, Australia

⁴²Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA

⁴³Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA

⁴⁴Department of Physics, The Ohio State University, Columbus, OH 43210, USA

⁴⁵Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

⁴⁶Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia

⁴⁷Lowell Observatory, 1400 Mars Hill Road, Flagstaff, AZ 86001, USA

⁴⁸Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA

⁴⁹Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro RJ—20921-400, Brazil

⁵⁰Department of Physics, University of Genova and INFN, Via Dodecaneso 33, I-16146 Genova, Italy

⁵¹Instituto de Física, UFRGS, Caixa Postal 15051, Porto Alegre RS—91501-970, Brazil

⁵²School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, UK
⁵³ Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA

⁵⁴ Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA

Received 2023 February 18; revised 2023 August 31; accepted 2023 October 16; published 2024 January 8

## Abstract

We address the problem of optimally identifying all kilonovae detected via gravitational-wave emission in the upcoming LIGO/Virgo/KAGRA observing run, O4, which is expected to be sensitive to a factor of ~7 more binary neutron star (BNS) alerts than previously. Electromagnetic follow-up of all but the brightest of these new events will require >1 m telescopes, for which limited time is available. We present an optimized observing strategy for the DECam during O4. We base our study on simulations of gravitational-wave events expected for O4 and wide-prior kilonova simulations. We derive the detectabilities of events for realistic observing conditions. We optimize our strategy for confirming a kilonova while minimizing telescope time. For a wide range of kilonova parameters, corresponding to a fainter kilonova compared to GW170817/AT 2017gfo, we find that, with this optimal strategy, the discovery probability for electromagnetic counterparts with the DECam is ~80% at the nominal BNS gravitational-wave detection limit for O4 (190 Mpc), which corresponds to an ~30% improvement compared to the strategy adopted during the previous observing run. For more distant events (~330 Mpc), we reach an ~60% probability of detection, a factor of ~2 increase. For a brighter kilonova model dominated by the blue component that reproduces the observations of GW170817/AT 2017gfo, we find that we can reach ~90% probability of detection out to 330 Mpc, representing an increase of ~20%, while also reducing the total telescope time required to follow up events by ~20%.

*Unified Astronomy Thesaurus concepts: Gravitational wave astronomy (675); Astronomical methods (1043); Astrostatistics tools (1887); Transient detection (1957)*

## 1. Introduction

The discovery of an electromagnetic (EM) counterpart to gravitational-wave (GW) event GW170817 ushered in a new era of astrophysics and multimessenger astronomy (Abbott et al. 2017b; Coulter et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017; Nicholl et al. 2017; Soares-Santos et al. 2017; Utsumi et al. 2017). The discovery of the optical counterpart occurred ~12 hr (Abbott et al. 2017b; Arcavi et al. 2017; Coulter et al. 2017; Lipunov et al. 2017; Soares-Santos et al. 2017) after merger, allowing for extensive imaging and spectroscopic observations (e.g., Andreoni et al. 2017; Díaz et al. 2017; Balasubramanian et al. 2021; Kilpatrick et al. 2022). This enabled the community to not only confirm the existence of long-hypothesized kilonovae (KNe) and characterize their light curves but also to derive the first cosmological standard siren (Schutz 1986) constraint (Abbott et al. 2017b).

The EM counterpart to GW170817, AT 2017gfo, was observed in γ-rays 2 s after the merger signal (e.g., Abbott et al. 2017a; Goldstein et al. 2017; Savchenko et al. 2017), then a few hours later in optical bands as a fast-decaying blue object, with longer-lived emission in the infrared, X-ray, and radio bands. The astrophysics community is still learning from the vast data set compiled for this event. Early analyses showed that AT 2017gfo had small ejecta masses with relativistic outflow velocities (from 0.1 to 0.3 of the speed of light c) and was likely powered by a combination of shocked material and r-process radioactive decay, which was modulated by the presence of high-opacity lanthanides synthesized during the

explosion and a highly nonspherical ejecta and jet (e.g., Drout et al. 2017; Kasen et al. 2017; Thielemann et al. 2017; Gottlieb et al. 2018). Many questions about this event and the diversity of KNe remain open, including the dependency on the neutron star (NS) equation of state, the geometry of the ejecta, the nature of the early blue emission, the physics of the relativistic jet launch, the exact contribution of KNe to the r-process element production, and the difference between KNe produced by binary neutron star (BNS) and neutron star–black hole (NS–BH) mergers.

The gravitational waveform observed from binary compact object mergers allows first-principles distance measurements. The GW distance measurement and the optical redshift of GW170817 provided constraints on the Hubble constant from a single event (Schutz 1986; Del Pozzo 2012; Abbott et al. 2017b), and further constraints came from measurements of $H_0$ from multiple GW events using a low-precision statistical method (i.e., dark sirens; Soares-Santos et al. 2019; Palmese et al. 2020, 2021a; Finke et al. 2021; LIGO Scientific Collaboration & Virgo Collaboration 2021; Gray et al. 2023).

From a cosmological perspective, the detection of GW170817 provided the grounds for the application of a new approach to understand the well-known Hubble tension. Measurements of the Hubble constant $H_0$ from early-time Universe observations (e.g., Planck observations of the cosmic microwave background in Planck Collaboration et al. 2020) are in tension with measurements from observations made at late times, e.g., the Hubble Space Telescope observations of Cepheids in concert with supernovae as standard candles (Riess et al. 2021). These measurements have been made at increasing precision over the last few years, only adding to the open questions surrounding the tension, and thus have provoked a series of investigations of variations on the cosmological model (for survey papers, see, e.g., Knox & Millea 2020; Abdalla et al. 2022; Schöneberg et al. 2022). There is no obvious resolution (Valentino et al. 2021), and the Hubble tension remains unresolved.

<figure><img src="image_3.png" /></figure>

<footnote n="55">⁵⁵ Legacy Survey of Space and Time Corporation Data Science Fellowship Program.</footnote>

<footnote n="56">⁵⁶ NASA Einstein Fellow.</footnote>

<footnote n="57">⁵⁷ In memoriam.</footnote>
⁵³ Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA

⁵⁴ Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA

Received 2023 February 18; revised 2023 August 31; accepted 2023 October 16; published 2024 January 8

## Abstract

We address the problem of optimally identifying all kilonovae detected via gravitational-wave emission in the upcoming LIGO/Virgo/KAGRA observing run, O4, which is expected to be sensitive to a factor of ~7 more binary neutron star (BNS) alerts than previously. Electromagnetic follow-up of all but the brightest of these new events will require >1 m telescopes, for which limited time is available. We present an optimized observing strategy for the DECam during O4. We base our study on simulations of gravitational-wave events expected for O4 and wide-prior kilonova simulations. We derive the detectabilities of events for realistic observing conditions. We optimize our strategy for confirming a kilonova while minimizing telescope time. For a wide range of kilonova parameters, corresponding to a fainter kilonova compared to GW170817/AT 2017gfo, we find that, with this optimal strategy, the discovery probability for electromagnetic counterparts with the DECam is ~80% at the nominal BNS gravitational-wave detection limit for O4 (190 Mpc), which corresponds to an ~30% improvement compared to the strategy adopted during the previous observing run. For more distant events (~330 Mpc), we reach an ~60% probability of detection, a factor of ~2 increase. For a brighter kilonova model dominated by the blue component that reproduces the observations of GW170817/AT 2017gfo, we find that we can reach ~90% probability of detection out to 330 Mpc, representing an increase of ~20%, while also reducing the total telescope time required to follow up events by ~20%.

*Unified Astronomy Thesaurus concepts: Gravitational wave astronomy (675); Astronomical methods (1043); Astrostatistics tools (1887); Transient detection (1957)*

## 1. Introduction

The discovery of an electromagnetic (EM) counterpart to gravitational-wave (GW) event GW170817 ushered in a new era of astrophysics and multimessenger astronomy (Abbott et al. 2017b; Coulter et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017; Nicholl et al. 2017; Soares-Santos et al. 2017; Utsumi et al. 2017). The discovery of the optical counterpart occurred ~12 hr (Abbott et al. 2017b; Arcavi et al. 2017; Coulter et al. 2017; Lipunov et al. 2017; Soares-Santos et al. 2017) after merger, allowing for extensive imaging and spectroscopic observations (e.g., Andreoni et al. 2017; Díaz et al. 2017; Balasubramanian et al. 2021; Kilpatrick et al. 2022). This enabled the community to not only confirm the existence of long-hypothesized kilonovae (KNe) and characterize their light curves but also to derive the first cosmological standard siren (Schutz 1986) constraint (Abbott et al. 2017b).

The EM counterpart to GW170817, AT 2017gfo, was observed in γ-rays 2 s after the merger signal (e.g., Abbott et al. 2017a; Goldstein et al. 2017; Savchenko et al. 2017), then a few hours later in optical bands as a fast-decaying blue object, with longer-lived emission in the infrared, X-ray, and radio bands. The astrophysics community is still learning from the vast data set compiled for this event. Early analyses showed that AT 2017gfo had small ejecta masses with relativistic outflow velocities (from 0.1 to 0.3 of the speed of light c) and was likely powered by a combination of shocked material and r-process radioactive decay, which was modulated by the presence of high-opacity lanthanides synthesized during the

explosion and a highly nonspherical ejecta and jet (e.g., Drout et al. 2017; Kasen et al. 2017; Thielemann et al. 2017; Gottlieb et al. 2018). Many questions about this event and the diversity of KNe remain open, including the dependency on the neutron star (NS) equation of state, the geometry of the ejecta, the nature of the early blue emission, the physics of the relativistic jet launch, the exact contribution of KNe to the r-process element production, and the difference between KNe produced by binary neutron star (BNS) and neutron star–black hole (NS–BH) mergers.

The gravitational waveform observed from binary compact object mergers allows first-principles distance measurements. The GW distance measurement and the optical redshift of GW170817 provided constraints on the Hubble constant from a single event (Schutz 1986; Del Pozzo 2012; Abbott et al. 2017b), and further constraints came from measurements of $H_0$ from multiple GW events using a low-precision statistical method (i.e., dark sirens; Soares-Santos et al. 2019; Palmese et al. 2020, 2021a; Finke et al. 2021; LIGO Scientific Collaboration & Virgo Collaboration 2021; Gray et al. 2023).

From a cosmological perspective, the detection of GW170817 provided the grounds for the application of a new approach to understand the well-known Hubble tension. Measurements of the Hubble constant $H_0$ from early-time Universe observations (e.g., Planck observations of the cosmic microwave background in Planck Collaboration et al. 2020) are in tension with measurements from observations made at late times, e.g., the Hubble Space Telescope observations of Cepheids in concert with supernovae as standard candles (Riess et al. 2021). These measurements have been made at increasing precision over the last few years, only adding to the open questions surrounding the tension, and thus have provoked a series of investigations of variations on the cosmological model (for survey papers, see, e.g., Knox & Millea 2020; Abdalla et al. 2022; Schöneberg et al. 2022). There is no obvious resolution (Valentino et al. 2021), and the Hubble tension remains unresolved.

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="55">⁵⁵ Legacy Survey of Space and Time Corporation Data Science Fellowship Program.</footnote>

<footnote n="56">⁵⁶ NASA Einstein Fellow.</footnote>

<footnote n="57">⁵⁷ In memoriam.</footnote>
The GW distances measured from compact binary mergers do not rely on the distance ladder and therefore will have different systematics than the astronomical distance ladder late-time measurements of $H_0$ and possibly smaller systematics. With additional standard siren measurements, the precision on the Hubble constant could reach the percent level (Chen et al. 2018; Bom & Palmese 2023). This level of precision would be an important contribution to resolving the Hubble tension (see discussion later in this section).

The third LIGO/Virgo/KAGRA (LVK) collaboration observing run (O3) resulted in over 60 new events (see the third Gravitational-wave Transient Catalog, GWTC; LIGO Scientific Collaboration et al. 2023). One GW event originated from a BNS merger, and two high-confidence events originated from NS-BH coalescences (Abbott et al. 2021), but no EM counterparts were confirmed from any event despite extensive follow-up campaigns (e.g., Andreoni et al. 2019, 2020; Goldstein et al. 2019; Garcia et al. 2020; Morgan et al. 2020; Kilpatrick et al. 2021; Oates et al. 2021; Tucker et al. 2022). One study (Graham et al. 2020) proposed an AGN counterpart to the binary black hole merger GW190521, but the association with the GW event cannot be made with confidence (Ashton et al. 2021; Palmese et al. 2021b).

LVK's O4 campaign was scheduled to start in early 2023, surveying on a factor of twice the median luminosity distance and a factor of ~7 detections of BNS events than the previous O3 campaign (see Table 2 in Petrov et al. 2022). The expected rates of BNSs are uncertain but in the range of 9–88 yr⁻¹ (Abbott et al. 2020; Petrov et al. 2022). It is unlikely that each subsequent KN event will obtain the same amount of follow-up resources as in O3. As GW detectors become more sensitive and able to detect events at larger distances, the optical follow-up of BNS events will become more challenging, including for campaigns using the Dark Energy Camera (DECam; Flaugher et al. 2015), such as those coordinated by the DES Gravitational Wave (DESGW) follow-up group.

Here, we present improved and optimized strategies for discovery of KNe. The methodology presented here focuses on DECam resources, which the DESGW group plans to use to follow up LVK O4 events. However, the method is generic and could be easily adapted to any telescope. Our primary science goal is the construction of the standard siren Hubble diagram via maximizing the number of GW-detected events with known redshifts.

Our paper is structured as follows. In Section 2, we describe our GW and KN simulations. In Section 3, we discuss how we measure the success of each observing strategy. In Section 4, we detail the various optimization options and types of KN models. In Section 5, we detail the strategies that we consider to be most successful for ensuring that DECam detects KNe with our follow-up observations. We additionally describe how this code may be used in a real-time follow-up in Section 6. Finally, we discuss and summarize our results in Sections 7 and 8.

## 2. Simulation Data

### 2.1. Simulated O4 BNS Events

We start by producing a set of simulated BNS mergers that are expected to be detectable in the upcoming LVK O4 observing run. The procedure is similar to that in Petrov et al. (2022).

<figure><img src="image_4.png" /><figcaption>**Figure 1.** The distance distribution and 90% credible interval sky localization area of the O4 simulations of BNS mergers used in this paper. The simulated primary masses follow a $x = 1.5$, $\sigma = 1.1$, $x_{min} = 1.1$, $x_{max} = 3M_\odot$ Gaussian distribution, and the spins (see Equation (2) of first reference in LIGO Scientific Collaboration et al. 2023) of the events are distributed uniformly in the range $\pm 0.05$. Simulated detections are limited to events with a network S/N greater than 12. Throughout the rest of this paper, we will assume follow-up of events only with areas of <300 deg², which is 611 of the 860 events.</figcaption></figure>

GW events are simulated using the BAYESTAR software (Singer & Price 2016; Singer et al. 2016a, 2016b), which uses LALSuite (LIGO Scientific Collaboration 2018) tools. We assume sensitivity curves for Advanced LIGO, Virgo, and KAGRA as O4 sensitivities discussed in Abbott et al. (2020)<fnref n="58" /> though we assume a sensitivity for KAGRA of a BNS range of ~80 Mpc.<fnref n="59" /> All detectors have a duty cycle of 70%, which is consistent with LVK predictions (Abbott et al. 2018). Assuming a Planck Collaboration et al. (2020) cosmology, we create 10,000 BNS events of the type O4 could theoretically observe, following a uniform in comoving volume distribution, and then inject them into the GW search-and-discovery pipeline (Herner et al. 2020). TaylorF2 waveforms (Buonanno et al. 2009) are assumed for both injections and reconstructions. The primary mass distribution of our injections follows the NS mass function found in Abbott et al. (2023), normally distributed with mean $1.5 M_\odot$ and standard deviation $1.1 M_\odot$, truncated to be within $1.1 M_\odot \le M_{NS} \le 3M_\odot$ in order to stay consistent with the Kasen models' parameter space. The NS spin distribution was uniformly distributed between -0.05 and 0.05. After injecting the BNS mergers, a matched filter search retrieves the detected events. We consider as detections those events for which a single-detector signal-to-noise ratio (S/N) >4 is reached by at least two detectors, and the overall network S/N is >12, resulting in 860 detected events. The measured S/N is added with Gaussian noise. Finally, we produce BAYESTAR sky maps for the detected events. In Figure 1, we present the area (90% credible interval) and luminosity distance (integrated over the whole sky) for all of the simulations used in this analysis.

<footnote n="58">58 The curves are available at https://dcc.ligo.org/LIGO-T2000012/public and https://dcc.ligo.org/LIGO-T2000012-v1/public.</footnote>

<footnote n="59">59 LVK has updated the assumed sensitivity for KAGRA, as seen in v2 of the in v2 of the table in https://dcc.ligo.org/LIGO-T2000012-v1/public; this makes little difference to our analysis, which is concerned mostly with distant events.</footnote>
The GW distances measured from compact binary mergers do not rely on the distance ladder and therefore will have different systematics than the astronomical distance ladder late-time measurements of $H_0$ and possibly smaller systematics. With additional standard siren measurements, the precision on the Hubble constant could reach the percent level (Chen et al. 2018; Bom & Palmese 2023). This level of precision would be an important contribution to resolving the Hubble tension (see discussion later in this section).

The third LIGO/Virgo/KAGRA (LVK) collaboration observing run (O3) resulted in over 60 new events (see the third Gravitational-wave Transient Catalog, GWTC; LIGO Scientific Collaboration et al. 2023). One GW event originated from a BNS merger, and two high-confidence events originated from NS-BH coalescences (Abbott et al. 2021), but no EM counterparts were confirmed from any event despite extensive follow-up campaigns (e.g., Andreoni et al. 2019, 2020; Goldstein et al. 2019; Garcia et al. 2020; Morgan et al. 2020; Kilpatrick et al. 2021; Oates et al. 2021; Tucker et al. 2022). One study (Graham et al. 2020) proposed an AGN counterpart to the binary black hole merger GW190521, but the association with the GW event cannot be made with confidence (Ashton et al. 2021; Palmese et al. 2021b).

LVK's O4 campaign was scheduled to start in early 2023, surveying on a factor of twice the median luminosity distance and a factor of ~7 detections of BNS events than the previous O3 campaign (see Table 2 in Petrov et al. 2022). The expected rates of BNSs are uncertain but in the range of 9–88 yr⁻¹ (Abbott et al. 2020; Petrov et al. 2022). It is unlikely that each subsequent KN event will obtain the same amount of follow-up resources as in O3. As GW detectors become more sensitive and able to detect events at larger distances, the optical follow-up of BNS events will become more challenging, including for campaigns using the Dark Energy Camera (DECam; Flaugher et al. 2015), such as those coordinated by the DES Gravitational Wave (DESGW) follow-up group.

Here, we present improved and optimized strategies for discovery of KNe. The methodology presented here focuses on DECam resources, which the DESGW group plans to use to follow up LVK O4 events. However, the method is generic and could be easily adapted to any telescope. Our primary science goal is the construction of the standard siren Hubble diagram via maximizing the number of GW-detected events with known redshifts.

Our paper is structured as follows. In Section 2, we describe our GW and KN simulations. In Section 3, we discuss how we measure the success of each observing strategy. In Section 4, we detail the various optimization options and types of KN models. In Section 5, we detail the strategies that we consider to be most successful for ensuring that DECam detects KNe with our follow-up observations. We additionally describe how this code may be used in a real-time follow-up in Section 6. Finally, we discuss and summarize our results in Sections 7 and 8.

## 2. Simulation Data

### 2.1. Simulated O4 BNS Events

We start by producing a set of simulated BNS mergers that are expected to be detectable in the upcoming LVK O4 observing run. The procedure is similar to that in Petrov et al. (2022).

<figure><img src="image_4.png" /><figcaption>**Figure 1.** The distance distribution and 90% credible interval sky localization area of the O4 simulations of BNS mergers used in this paper. The simulated primary masses follow a $x = 1.5$, $\sigma = 1.1$, $x_{min} = 1.1$, $x_{max} = 3M_\odot$ Gaussian distribution, and the spins (see Equation (2) of first reference in LIGO Scientific Collaboration et al. 2023) of the events are distributed uniformly in the range $\pm 0.05$. Simulated detections are limited to events with a network S/N greater than 12. Throughout the rest of this paper, we will assume follow-up of events only with areas of <300 deg², which is 611 of the 860 events.</figcaption></figure>

GW events are simulated using the BAYESTAR software (Singer & Price 2016; Singer et al. 2016a, 2016b), which uses LALSuite (LIGO Scientific Collaboration 2018) tools. We assume sensitivity curves for Advanced LIGO, Virgo, and KAGRA as O4 sensitivities discussed in Abbott et al. (2020)<fnref n="58" /> though we assume a sensitivity for KAGRA of a BNS range of ~80 Mpc.<fnref n="59" /> All detectors have a duty cycle of 70%, which is consistent with LVK predictions (Abbott et al. 2018). Assuming a Planck Collaboration et al. (2020) cosmology, we create 10,000 BNS events of the type O4 could theoretically observe, following a uniform in comoving volume distribution, and then inject them into the GW search-and-discovery pipeline (Herner et al. 2020). TaylorF2 waveforms (Buonanno et al. 2009) are assumed for both injections and reconstructions. The primary mass distribution of our injections follows the NS mass function found in Abbott et al. (2023), normally distributed with mean $1.5 M_\odot$ and standard deviation $1.1 M_\odot$, truncated to be within $1.1 M_\odot \le M_{NS} \le 3M_\odot$ in order to stay consistent with the Kasen models' parameter space. The NS spin distribution was uniformly distributed between -0.05 and 0.05. After injecting the BNS mergers, a matched filter search retrieves the detected events. We consider as detections those events for which a single-detector signal-to-noise ratio (S/N) >4 is reached by at least two detectors, and the overall network S/N is >12, resulting in 860 detected events. The measured S/N is added with Gaussian noise. Finally, we produce BAYESTAR sky maps for the detected events. In Figure 1, we present the area (90% credible interval) and luminosity distance (integrated over the whole sky) for all of the simulations used in this analysis.

<footnote n="58">58 The curves are available at https://dcc.ligo.org/LIGO-T2000012/public and https://dcc.ligo.org/LIGO-T2000012-v1/public.</footnote>

<footnote n="59">59 LVK has updated the assumed sensitivity for KAGRA, as seen in v2 of the in v2 of the table in https://dcc.ligo.org/LIGO-T2000012-v1/public; this makes little difference to our analysis, which is concerned mostly with distant events.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** A sample of 10 randomly chosen sets of simulated GW localization maps used in this work. The inner (outer) lines represents 50% (90%) confidence level regions.</figcaption></figure>

Petrov et al. (2022) argues that the alerts produced by the Ligo-Virgo Collaboration (LVC) during O3 are better modeled by dropping the two-detector coincident detection requirement and using a minimum S/N for BNS events of >8. These criteria would have the effect of increasing the number of low-S/N, and therefore large sky area, events in our simulations. Figure 1 shows a median luminosity distance of ~150 Mpc and sky area of ~20 deg², whereas Petrov et al. (2022) find 352 ± 10 Mpc and 1820±170 deg², a difference from previous work that they attribute to the changing of the S/N requirements. There are 249 of the total of 860 events that have 90% sky area >300 deg². To change our sample to have a median sky area of ~350 deg², we would need an additional ≈360 events, all of which would need 90% sky area >350 deg². In our strategy definitions and proposed decision-making process to trigger limited target-of-opportunity (ToO) resources, we choose a sky area limit of <300 deg², so we choose not to consider these additional high sky area events here. We do not make any other selections in the events. In Figure 2, we present a set of 10 randomly chosen GW localization skymaps simulated by our pipeline.

## 2.2. KN Physical Models

We model KNe using the time-evolving theoretical spectral energy distributions (SEDs) of KN atmospheres from Kasen et al. (2017), hereafter referred to as the Kasen models. The Kasen models are parameterized by the mass ejected in the explosion, $M_{ej}$; the abundance of lanthanide elements in the ejecta, $X_{lan}$; and the velocity of the ejecta, $v_{ej}$. The models are a set of 329 time-dependent SEDs on a grid of discretized parameters: $0.001 M_{⊙} \le M_{ej} \le 0.1 M_{⊙}$, $0.03c \le v_{ej} \le 0.3c$, and $1 \times 10^{-9} \le X_{lan} \le 1 \times 10^{-1}$.

There are other KN models available in the literature (e.g., Bulla 2019; Darbha & Kasen 2020; Hotokezaka & Nakar 2020; Wollaeger et al. 2021; Gillanders et al. 2022). In comparison to these, the Kasen models do not build in the dependence on the geometry of the mergers and the viewing angle (see, e.g., Stewart et al. 2022 for a visualization of an asymmetric KN supported by accretion disk simulations) but instead provide atmosphere models to be used in building a geometric model.

## 2.3. KN Light-curve Simulations

The light-curve simulation pipeline is similar to that used by Morgan et al. (2020), which constrains the physical KN properties of GW190814. The SuperNova ANALysis (SNANA) software (Kessler et al. 2009, 2019) enables the simulation of light curves of KNe as they would be measured by DECam. SNANA produces light curves by simulating fluxes and uncertainties in observations by incorporating information about cadence, image zero-points, and noise levels in search and template images. The light curves are in absolute magnitudes and converted to observed magnitudes using a given cosmology. In particular, SNANA chooses a grid of 15 redshifts at $0.003 \le z \le 0.2$, which is used to transform the Kasen model SEDs without evolution. The Kasen model SED is redshifted to a z on our grid. SNANA also takes into consideration the reddening by the dust in the Milky Way (O'Donnell 1994; Kessler et al. 2019). Later, the SEDs are convolved with DECam transmission curves accounting for atmosphere, telescope, filters, and CCDs. Given a cadence, the light curves are calculated by sampling the magnitude grid.

The 329 Kasen models, $\alpha_j$, each with a range of time after burst, $\tau$ (0.00 days $\le \tau \le$ 16.7 days), at the SNANA grid of redshifts, $z_i$ (0.003 $\le z \le 0.2$), describe 829,080 simulated magnitudes, $m_{\lambda}(\lambda, \tau, \alpha_j, z_i)$, for each of the four filters that we use (g, r, i, z).

## 2.4. DECam Limiting Magnitudes

We define the limiting magnitude as the magnitude at which we can measure a point source, e.g., a star, with 0.1 mag error, which corresponds to S/N ~ 10. We then define the limiting magnitude $m_0$ for a total effective exposure time of 90 s, i.e., $t_{exp} \times t_{eff} = 90$ s, where $t_{eff}$ is a unitless quantity that scales the exposure time to the “effective” exposure time when taking into account sky conditions (higher $t_{eff}$ being better sky conditions). Then, to scale up for different exposures, we construct

$$m_{\text{lim}} = m_0 + 1.25 \log \left( t_{\text{eff}} \times \frac{t_{\text{exp}}}{90 \text{ s}} \right), \quad (1)$$

which reflects S/N going as the square root of time in a sky noise-limited observation. We use the $m_0$ measured in the DES data by
<figure><img src="image_5.png" /><figcaption>**Figure 2.** A sample of 10 randomly chosen sets of simulated GW localization maps used in this work. The inner (outer) lines represents 50% (90%) confidence level regions.</figcaption></figure>

Petrov et al. (2022) argues that the alerts produced by the Ligo-Virgo Collaboration (LVC) during O3 are better modeled by dropping the two-detector coincident detection requirement and using a minimum S/N for BNS events of >8. These criteria would have the effect of increasing the number of low-S/N, and therefore large sky area, events in our simulations. Figure 1 shows a median luminosity distance of ~150 Mpc and sky area of ~20 deg², whereas Petrov et al. (2022) find 352 ± 10 Mpc and 1820+190-170 deg², a difference from previous work that they attribute to the changing of the S/N requirements. There are 249 of the total of 860 events that have 90% sky area >300 deg². To change our sample to have a median sky area of ~350 deg², we would need an additional ≈360 events, all of which would need 90% sky area >350 deg². In our strategy definitions and proposed decision-making process to trigger limited target-of-opportunity (ToO) resources, we choose a sky area limit of <300 deg², so we choose not to consider these additional high sky area events here. We do not make any other selections in the events. In Figure 2, we present a set of 10 randomly chosen GW localization skymaps simulated by our pipeline.

## 2.2. KN Physical Models

We model KNe using the time-evolving theoretical spectral energy distributions (SEDs) of KN atmospheres from Kasen et al. (2017), hereafter referred to as the Kasen models. The Kasen models are parameterized by the mass ejected in the explosion, $M_{ej}$; the abundance of lanthanide elements in the ejecta, $X_{lan}$; and the velocity of the ejecta, $v_{ej}$. The models are a set of 329 time-dependent SEDs on a grid of discretized parameters: $0.001 M_{⊙} \le M_{ej} \le 0.1 M_{⊙}$, $0.03c \le v_{ej} \le 0.3c$, and $1 \times 10^{-9} \le X_{lan} \le 1 \times 10^{-1}$.

There are other KN models available in the literature (e.g., Bulla 2019; Darbha & Kasen 2020; Hotokezaka & Nakar 2020; Wollaeger et al. 2021; Gillanders et al. 2022). In comparison to these, the Kasen models do not build in the dependence on the geometry of the mergers and the viewing angle (see, e.g., Stewart et al. 2022 for a visualization of an asymmetric KN supported by accretion disk simulations) but instead provide atmosphere models to be used in building a geometric model.

## 2.3. KN Light-curve Simulations

The light-curve simulation pipeline is similar to that used by Morgan et al. (2020), which constrains the physical KN properties of GW190814. The SuperNova ANALysis (SNANA) software (Kessler et al. 2009, 2019) enables the simulation of light curves of KNe as they would be measured by DECam. SNANA produces light curves by simulating fluxes and uncertainties in observations by incorporating information about cadence, image zero-points, and noise levels in search and template images. The light curves are in absolute magnitudes and converted to observed magnitudes using a given cosmology. In particular, SNANA chooses a grid of 15 redshifts at $0.003 \le z \le 0.2$, which is used to transform the Kasen model SEDs without evolution. The Kasen model SED is redshifted to a z on our grid. SNANA also takes into consideration the reddening by the dust in the Milky Way (O'Donnell 1994; Kessler et al. 2019). Later, the SEDs are convolved with DECam transmission curves accounting for atmosphere, telescope, filters, and CCDs. Given a cadence, the light curves are calculated by sampling the magnitude grid.

The 329 Kasen models, $\alpha_j$, each with a range of time after burst, $\tau$ (0.00 days $\le \tau \le 16.7$ days), at the SNANA grid of redshifts, $z_i$ (0.003 $\le z \le 0.2$), describe 829,080 simulated magnitudes, $m_{\lambda}(\lambda, \tau, \alpha_j, z_i)$, for each of the four filters that we use (g, r, i, z).

## 2.4. DECam Limiting Magnitudes

We define the limiting magnitude as the magnitude at which we can measure a point source, e.g., a star, with 0.1 mag error, which corresponds to S/N ~ 10. We then define the limiting magnitude $m_0$ for a total effective exposure time of 90 s, i.e., $t_{exp} \times t_{eff} = 90$ s, where $t_{eff}$ is a unitless quantity that scales the exposure time to the “effective” exposure time when taking into account sky conditions (higher $t_{eff}$ being better sky conditions). Then, to scale up for different exposures, we construct

$$m_{\text{lim}} = m_0 + 1.25 \log \left( t_{\text{eff}} \times \frac{t_{\text{exp}}}{90 \text{ s}} \right), \quad (1)$$

which reflects S/N going as the square root of time in a sky noise-limited observation. We use the $m_0$ measured in the DES data by
<table><caption>**Table 1**
Observational Conditions (Photometric Filter, Limiting Magnitude, and Effective Exposure Time) Averaged from DECam Follow-ups of Previous Observations</caption><thead><tr><th>Filter</th><th>m<sub>0</sub></th><th>t<sub>eff</sub> (Bright Time)</th><th>t<sub>eff</sub> (Dark Time)</th></tr></thead><tbody><tr><td>g</td><td>23.4</td><td>0.05</td><td>0.7</td></tr><tr><td>r</td><td>23.1</td><td>0.15</td><td>0.8</td></tr><tr><td>i</td><td>22.5</td><td>0.45</td><td>0.7</td></tr><tr><td>z</td><td>21.8</td><td>0.6</td><td>0.6</td></tr></tbody></table>

**Note.** Note that the higher the $t_{eff}$ value (which ranges from 0 to 1), the better the observational conditions.

Neilsen et al. (2016), thus the 90 s normalization factor, and present those values in Table 1 and the derived limiting magnitudes in Table 2. The $t_{eff}$ is closely related to the observational conditions during the night. Therefore, we break observing into nights of bright time and dark time and use $t_{eff}$ from previous ToO programs in DECam, in particular the observations from past GW follow-up events (Garcia et al. 2020; Morgan et al. 2020). Although we split our tests into only dark and bright time, rather than on gray nights, which correspond to ~50% of the nights available in telescopes, we show in Section 5.2 that the difference in performance between dark and bright and therefore gray is negligible due to the adaptability of the presented method. The differences are in the particular configurations chosen, e.g., the filters and exposures times selected in different observational conditions. Throughout this paper, we focus on detectability on dark nights unless otherwise stated.

## 2.5. Cadences and Observational Parameters

We use SNANA to find effective search strategies that maximize candidate detection considering realistic conditions, including the maximum duration of the night, intervals between observations, and sky brightness, while also minimizing telescope time and enabling earlier discovery. This is at base a trade-off of exposure times versus sky area coverage but with additional complications of filter choice and the time since the event occurred. We limit our telescope time expenditures to 8 hr night⁻¹ and assume that Blanco/DECam has a telescope/readout slew time of 30 s between the exposures, which is true for small slews, such as the ones less than ~10°. In fact, long slews of ~100°, which might be necessary to cover disjoint GW maps, take on the order of 3 minutes. Thus, if there are two disjoint regions, which would add one or possibly two long slews depending on the observational conditions, it would add a negligible amount of total telescope time, relatively. We test four filters, *g*, *r*, *i*, and *z*, starting 12 hr after the trigger and going to 4 days postmerger in half-day increments with several exposure times. These are summarized in Table 3.

Note that the time required to respond to a GW alert depends upon several factors, including the human decision-making to trigger the ToO program, when the event is visible in the sky, observation planning, and time to ask for and obtain the ToO interrupt. The time required to get the telescope on-sky also depends on when the alert is given. It is possible that it may be significantly less than 12 hr until we are on-sky; however, for this study, we opted for a conservative option of 12 hr, by which time, for any significant and observable event, we should be able to observe.

<table><caption>**Table 2**
DECam $m_{lim}$(10σ), Dark Time</caption><thead><tr><th rowspan="2">Filter</th><th colspan="4">Exposure Time (s)</th></tr><tr><th>10</th><th>100</th><th>1200</th><th>3600</th></tr></thead><tbody><tr><td>g</td><td>22.0</td><td>23.3</td><td>24.6</td><td>25.2</td></tr><tr><td>r</td><td>21.8</td><td>23.0</td><td>24.4</td><td>25.0</td></tr><tr><td>i</td><td>21.1</td><td>22.4</td><td>23.7</td><td>24.3</td></tr><tr><td>z</td><td>20.3</td><td>21.6</td><td>22.9</td><td>23.5</td></tr></tbody></table>

**Note.** $\Delta m_{lim}$(5σ − 10σ) = 0.75 mag.

<table><caption>**Table 3**
Observational Parameter Space $\Theta$ Explored</caption><tbody><tr><td>Time after Burst (days)</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0</td></tr><tr><td colspan="2">Scenario 1</td></tr><tr><td>Filter</td><td>g, r, i, z</td></tr><tr><td>2D credible sky area covered</td><td>0.9, 0.85, 0.8, 0.75, 0.7</td></tr><tr><td>Exposure time (s)</td><td>60, 90, 120, 200, 300, 600, 1200, 3600</td></tr><tr><td colspan="2">Scenario 2 (with inner and outer region)</td></tr><tr><td>Filter</td><td>g, r, i, z</td></tr><tr><td>2D credible sky area covered</td><td></td></tr><tr><td>Outer: 0.9</td><td>Inner: 0.8, 0.7, 0.5, 0.3</td></tr><tr><td>Outer: 0.8</td><td>Inner: 0.7, 0.5, 0.3</td></tr><tr><td>Outer: 0.7</td><td>Inner: 0.5, 0.4, 0.3</td></tr><tr><td>Outer/inner exposure time pairs (s)</td><td>60/90; 90/120; 120/200; 200/300; 300/600; 300/1200; 600/1200; 600/2400; 1200/2400; 2400/3600; 3600/5400</td></tr></tbody></table>

We define two exposure time scenarios. In scenario 1, we cover the area of a given GW event with a single set of exposure times, and in scenario 2, we explore the use of two different exposure times for a single search. The latter is motivated by the need to cover the high-probability area sky with deep exposures while covering the larger low-probability localization area outskirts with shorter exposure images. We designate the central high-probability areas as the “inner region” and the rest of the area inside the localization region as the “outer region.” The last section of Table 3 presents the combinations considered for the inner region, ranging from 30% to 80% sky probability coverage, for three different values of the total (deep+shallow) sky map probability coverage, from 70% to 90%. For instance, a combination of 40% probability for the inner region and 70% total coverage means that the 40% highest-probability region is covered with a higher exposure time, and the 70% − 40% = 30% left over is covered with the shorter exposures. Each of these combinations is considered for all possible deep and shallow exposures presented in Table 3. For scenarios where we cover the sky area twice in a single night, we additionally take into consideration the KN variability a few hours after the first search.

## 2.6. Simulation Data Summary

We simulate a set of GW detections $S = S_1, ..., S_n$, where $n = 611$, and $S_i$ is the ith simulation with distance $d_i$. For each $S_i$, we evaluate each of the parameter sets $\Theta$ in the two scenarios. In scenario 1, there are eight passes since burst, four filters, eight exposure times, and five sky area probability
<table><caption>**Table 1**
Observational Conditions (Photometric Filter, Limiting Magnitude, and Effective Exposure Time) Averaged from DECam Follow-ups of Previous Observations</caption><thead><tr><th>Filter</th><th>m<sub>0</sub></th><th>t<sub>eff</sub> (Bright Time)</th><th>t<sub>eff</sub> (Dark Time)</th></tr></thead><tbody><tr><td>g</td><td>23.4</td><td>0.05</td><td>0.7</td></tr><tr><td>r</td><td>23.1</td><td>0.15</td><td>0.8</td></tr><tr><td>i</td><td>22.5</td><td>0.45</td><td>0.7</td></tr><tr><td>z</td><td>21.8</td><td>0.6</td><td>0.6</td></tr></tbody></table>

**Note.** Note that the higher the $t_{eff}$ value (which ranges from 0 to 1), the better the observational conditions.

Neilsen et al. (2016), thus the 90 s normalization factor, and present those values in Table 1 and the derived limiting magnitudes in Table 2. The $t_{eff}$ is closely related to the observational conditions during the night. Therefore, we break observing into nights of bright time and dark time and use $t_{eff}$ from previous ToO programs in DECam, in particular the observations from past GW follow-up events (Garcia et al. 2020; Morgan et al. 2020). Although we split our tests into only dark and bright time, rather than on gray nights, which correspond to ~50% of the nights available in telescopes, we show in Section 5.2 that the difference in performance between dark and bright and therefore gray is negligible due to the adaptability of the presented method. The differences are in the particular configurations chosen, e.g., the filters and exposures times selected in different observational conditions. Throughout this paper, we focus on detectability on dark nights unless otherwise stated.

## 2.5. Cadences and Observational Parameters

We use SNANA to find effective search strategies that maximize candidate detection considering realistic conditions, including the maximum duration of the night, intervals between observations, and sky brightness, while also minimizing telescope time and enabling earlier discovery. This is at base a trade-off of exposure times versus sky area coverage but with additional complications of filter choice and the time since the event occurred. We limit our telescope time expenditures to 8 hr night⁻¹ and assume that Blanco/DECam has a telescope/readout slew time of 30 s between the exposures, which is true for small slews, such as the ones less than ~10°. In fact, long slews of ~100°, which might be necessary to cover disjoint GW maps, take on the order of 3 minutes. Thus, if there are two disjoint regions, which would add one or possibly two long slews depending on the observational conditions, it would add a negligible amount of total telescope time, relatively. We test four filters, *g*, *r*, *i*, and *z*, starting 12 hr after the trigger and going to 4 days postmerger in half-day increments with several exposure times. These are summarized in Table 3.

Note that the time required to respond to a GW alert depends upon several factors, including the human decision-making to trigger the ToO program, when the event is visible in the sky, observation planning, and time to ask for and obtain the ToO interrupt. The time required to get the telescope on-sky also depends on when the alert is given. It is possible that it may be significantly less than 12 hr until we are on-sky; however, for this study, we opted for a conservative option of 12 hr, by which time, for any significant and observable event, we should be able to observe.

<table><caption>**Table 2**
DECam $m_{lim}$(10σ), Dark Time</caption><thead><tr><th rowspan="2">Filter</th><th colspan="4">Exposure Time (s)</th></tr><tr><th>10</th><th>100</th><th>1200</th><th>3600</th></tr></thead><tbody><tr><td>g</td><td>22.0</td><td>23.3</td><td>24.6</td><td>25.2</td></tr><tr><td>r</td><td>21.8</td><td>23.0</td><td>24.4</td><td>25.0</td></tr><tr><td>i</td><td>21.1</td><td>22.4</td><td>23.7</td><td>24.3</td></tr><tr><td>z</td><td>20.3</td><td>21.6</td><td>22.9</td><td>23.5</td></tr></tbody></table>

**Note.** $\Delta m_{lim}$(5σ − 10σ) = 0.75 mag.

<table><caption>**Table 3**
Observational Parameter Space $\Theta$ Explored</caption><tbody><tr><td>Time after Burst (days)</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0</td></tr><tr><td colspan="2">Scenario 1</td></tr><tr><td>Filter</td><td>g, r, i, z</td></tr><tr><td>2D credible sky area covered</td><td>0.9, 0.85, 0.8, 0.75, 0.7</td></tr><tr><td>Exposure time (s)</td><td>60, 90, 120, 200, 300, 600, 1200, 3600</td></tr><tr><td colspan="2">Scenario 2 (with inner and outer region)</td></tr><tr><td>Filter</td><td>g, r, i, z</td></tr><tr><td>2D credible sky area covered</td><td></td></tr><tr><td>Outer: 0.9</td><td>Inner: 0.8, 0.7, 0.5, 0.3</td></tr><tr><td>Outer: 0.8</td><td>Inner: 0.7, 0.5, 0.3</td></tr><tr><td>Outer: 0.7</td><td>Inner: 0.5, 0.4, 0.3</td></tr><tr><td>Outer/inner exposure time pairs (s)</td><td>60/90; 90/120; 120/200; 200/300; 300/600; 300/1200; 600/1200; 600/2400; 1200/2400; 2400/3600; 3600/5400</td></tr></tbody></table>

We define two exposure time scenarios. In scenario 1, we cover the area of a given GW event with a single set of exposure times, and in scenario 2, we explore the use of two different exposure times for a single search. The latter is motivated by the need to cover the high-probability area sky with deep exposures while covering the larger low-probability localization area outskirts with shorter exposure images. We designate the central high-probability areas as the “inner region” and the rest of the area inside the localization region as the “outer region.” The last section of Table 3 presents the combinations considered for the inner region, ranging from 30% to 80% sky probability coverage, for three different values of the total (deep+shallow) sky map probability coverage, from 70% to 90%. For instance, a combination of 40% probability for the inner region and 70% total coverage means that the 40% highest-probability region is covered with a higher exposure time, and the 70% − 40% = 30% left over is covered with the shorter exposures. Each of these combinations is considered for all possible deep and shallow exposures presented in Table 3. For scenarios where we cover the sky area twice in a single night, we additionally take into consideration the KN variability a few hours after the first search.

## 2.6. Simulation Data Summary

We simulate a set of GW detections $S = S_1, ..., S_n$, where $n = 611$, and $S_i$ is the ith simulation with distance $d_i$. For each $S_i$, we evaluate each of the parameter sets $\Theta$ in the two scenarios. In scenario 1, there are eight passes since burst, four filters, eight exposure times, and five sky area probability
<table><caption>**Table 4**
KN Gaussian Model Priors</caption><thead><tr><th rowspan="2"></th><th colspan="2">Blue and Bright</th><th colspan="2">Reddish and Slow</th><th colspan="2">Red and Faint</th></tr><tr><th>μ</th><th>σ</th><th>μ</th><th>σ</th><th>μ</th><th>σ</th></tr></thead><tbody><tr><td>M<sub>ej</sub>(M<sub>☉</sub>)</td><td>0.025</td><td>0.001</td><td>0.025</td><td>0.01</td><td>0.035</td><td>0.15</td></tr><tr><td>log(X<sub>lan</sub>)</td><td>-5.0</td><td>1.0</td><td>-5.0</td><td>10.0</td><td>-2.0</td><td>5.0</td></tr><tr><td>v<sub>ej</sub> (c)</td><td>0.25</td><td>0.01</td><td>0.25</td><td>0.10</td><td>0.15</td><td>0.30</td></tr><tr><td>Peak &lang;M<sub>i</sub>&rang;</td><td colspan="2">-16.3</td><td colspan="2">-15.7</td><td colspan="2">-14.6</td></tr></tbody></table>

Note. Table of means, μ, and standard deviations, σ, of the Gaussian priors on the Kasen model parameter ranges. The Kasen models have parameter ranges of 0.001 $M_{☉} \le M_{ej} \le 0.1 M_{☉}$, $-9 \le \log(X_{lan}) \le -1$, and $0.03c \le v_{ej} \le 0.3c$.

coverages, i.e., $8 \cdot 4 \cdot 8 \cdot 5 = 1280$ possible parameter sets. In scenario 2, there are $8 \cdot 4 \cdot 11 \cdot 10 = 3520$ possible parameter sets to be evaluated. Over both, 4800 observation models are evaluated. Each model is evaluated with the machinery described in Section 2, resulting in SNANA KN-measured magnitudes. We have done this for each of the 329 Kasen models.

# 3. Detection Methodology

## 3.1. Discovery Probability

We define the probability of detection, $p_{\alpha_j}$, of the KN model $\alpha_j$ for the jth combination of ($M_{ej}$, log($X_{lan}$), $v_{ej}$), with observed magnitude $m_{\lambda} = m(\lambda, \tau, \alpha_j, z_i)$ in a given filter $\lambda$ for a given exposure time $t_{exp}$ weighted by observing condition $t_{eff}$ for an event at mean redshift, $\bar{z}$, over the SNANA grid of redshifts, $z_i$:

$$p_{\alpha_j} \propto \sum_i \begin{cases} \text{pr}(\alpha_j) \cdot \text{pr}(z_i), & \text{if } m_{\lambda} < m_{\text{lim}} \\ 0, & \text{if } m_{\lambda} \ge m_{\text{lim}} \end{cases}, \quad (2)$$

where

$$\begin{align*} \alpha_j &\equiv \text{jth combination of } (M_{ej}, \log(X_{lan}), v_{ej}), \\ p_{\alpha_j} &\equiv p_{\alpha_j}(\text{detection}|\tau, t_{\text{eff}}, \Theta, \bar{z}), \\ m_{\lambda} &\equiv \text{mag}_{\lambda}(\lambda, \tau, \alpha_j, z_i), \\ m_{\text{lim}} &\equiv \text{mag}_{\text{lim}}(\lambda, t_{\text{eff}}, \Theta), \quad \text{and} \end{align*}$$

$$\text{pr}(z_i) = \exp\left(-\frac{(z_i - \bar{z})^2}{2\sigma_{\bar{z}}^2}\right),$$

where $\tau$ is the time after the GW detection; $\Theta$ contains the specific observation strategy characteristics, including $t_{exp}$, filter, and area coverage as described in Table 3; $m_{lim}$ is the limiting magnitude of the observation for a $10\sigma$ detection; and $\text{pr}(\alpha_j)$ are the model priors defined by each KN model described in Section 4.2 and Table 4. The summation is over all redshifts for a given model, j, at a given $\tau$. Equation (2) represents a Gaussian prior for choosing light-curve models from the grid of SNANA-defined redshifts given the GW distance. Explaining it differently, Equation (2) is, for a given Kasen model, examining whether the resulting apparent magnitude is less than the limiting model with a Gaussian prior on the redshift, using the GW event mean redshift, $\bar{z}$, and the grid of SNANA redshifts the KN could be at, weighted by the GW event variance in redshift, $\sigma_{\bar{z}}^2$.

The GW localization maps present the probability that the event is located at a given sky position, the luminosity distance at that position, and its uncertainty. Therefore, we define the total probability of detecting an event as

$$P_d = \frac{\int_{\hat{\Omega}} d\Omega d_L(\Omega) p(\Omega)}{\int_{\hat{\Omega}} d\Omega d_L(\Omega)} \times \sum_{j=1}^{329} p_{\alpha_j}, \quad (3)$$

where

$$P_d \equiv P(\text{discovery}|\hat{\Omega}, \tau, t_{\text{eff}}, \Theta, \bar{z}),$$

and $\hat{\Omega}$ is the entire sky area observed in the follow-up, $d\Omega$ is the voxel, $p(\Omega)$ is the probability in the voxel, and $d_L(d\Omega)$ is the luminosity distance to the voxel. The sum over $p_{\alpha_j}$ includes the priors and thus indicates the model used. The sky coverage and exposure times determine the total telescope time for a given KN detection. We note that we use weighted spatial probability, rather than 2D on-sky probabilities; therefore, we have the $d_L$ weighting in Equation (3) and the attendant implication that we can have a higher detection probability than on-sky credible area covered. We can now evaluate $P_d$ for a given set of $\Theta$.

## 3.2. Confirmation Probability

In order for an object to be confirmed as a KN candidate, we require it be detected twice (in two observing “epochs”). This requirement can be lifted if, for example, there are sufficient spectroscopic resources available to follow up all of the candidates found after a single epoch. Generally, this is not the case, however. The second detection eliminates spurious detections, including image artifacts, asteroids, and other possible contaminants (Morgan et al. 2020; Shandonay et al. 2022). Another reason to consider the detections independent is that we typically observe while working on postprocessing and making target selection for spectroscopy in the data from the previous epoch.

Given that we want to make two detections to positively identify KN candidates, we define the probability of confirming the transient with two independent detections as

$$P_c \equiv P(\text{confirmed}) = P(\text{discovery}|\hat{\Omega}, \tau, t_{\text{eff}}, \Theta, \bar{z}) \times P(\text{discovery}|\tau_2, t_{\text{eff}}, \Theta, \bar{z}), \quad (4)$$

where

$$P(\text{discovery}|\tau, t_{\text{eff}}, \Theta, \bar{z}) = \sum_{\alpha_j} p_{\alpha_j}.$$

We design the strategies along the paper optimizing for $P_c$ given a set of constraints. Due to implementation choices and to speed up the numerical optimization, we use $P_c^* \equiv P_{d,1} \cdot P_{d,2}$, which gives us equivalent results in terms of strategy,<fnref n="60" /> such as telescope time or duration of the night. We further present our main findings as a function of discovery probability, which is the most relevant outcome for the proposed strategies.

<footnote n="60">⁶⁰ In this form, $P^*(\text{confirmed})$ has an extra spatial/volume probability term multiplied by $V_p$, in which $P_c^* = V_p \cdot P_c$ is the same for both $P_{d1}$ and $P_{d2}$. Therefore, this is irrelevant from an optimization viewpoint, since arg $\max(P_c^*) = \arg \max(P_c)$.</footnote>
<table><caption>**Table 4**
KN Gaussian Model Priors</caption><thead><tr><th rowspan="2"></th><th colspan="2">Blue and Bright</th><th colspan="2">Reddish and Slow</th><th colspan="2">Red and Faint</th></tr><tr><th>μ</th><th>σ</th><th>μ</th><th>σ</th><th>μ</th><th>σ</th></tr></thead><tbody><tr><td>M<sub>ej</sub>(M<sub>☉</sub>)</td><td>0.025</td><td>0.001</td><td>0.025</td><td>0.01</td><td>0.035</td><td>0.15</td></tr><tr><td>log(X<sub>lan</sub>)</td><td>-5.0</td><td>1.0</td><td>-5.0</td><td>10.0</td><td>-2.0</td><td>5.0</td></tr><tr><td>v<sub>ej</sub> (c)</td><td>0.25</td><td>0.01</td><td>0.25</td><td>0.10</td><td>0.15</td><td>0.30</td></tr><tr><td>Peak &lang;M<sub>i</sub>&rang;</td><td colspan="2">-16.3</td><td colspan="2">-15.7</td><td colspan="2">-14.6</td></tr></tbody></table>

Note. Table of means, μ, and standard deviations, σ, of the Gaussian priors on the Kasen model parameter ranges. The Kasen models have parameter ranges of 0.001 $M_{☉} \le M_{ej} \le 0.1 M_{☉}$, $-9 \le \log(X_{lan}) \le -1$, and $0.03c \le v_{ej} \le 0.3c$.

coverages, i.e., $8 \cdot 4 \cdot 8 \cdot 5 = 1280$ possible parameter sets. In scenario 2, there are $8 \cdot 4 \cdot 11 \cdot 10 = 3520$ possible parameter sets to be evaluated. Over both, 4800 observation models are evaluated. Each model is evaluated with the machinery described in Section 2, resulting in SNANA KN-measured magnitudes. We have done this for each of the 329 Kasen models.

# 3. Detection Methodology

## 3.1. Discovery Probability

We define the probability of detection, $p_{\alpha_j}$, of the KN model $\alpha_j$ for the jth combination of ($M_{ej}$, log($X_{lan}$), $v_{ej}$), with observed magnitude $m_{\lambda} = m(\lambda, \tau, \alpha_j, z_i)$ in a given filter $\lambda$ for a given exposure time $t_{exp}$ weighted by observing condition $t_{eff}$ for an event at mean redshift, $\bar{z}$, over the SNANA grid of redshifts, $z_i$:

$$p_{\alpha_j} \propto \sum_i \begin{cases} \text{pr}(\alpha_j) \cdot \text{pr}(z_i), & \text{if } m_{\lambda} < m_{\text{lim}} \\ 0, & \text{if } m_{\lambda} \ge m_{\text{lim}} \end{cases}, \quad (2)$$

where

$$\begin{align*} \alpha_j &\equiv \text{jth combination of } (M_{ej}, \log(X_{lan}), v_{ej}), \\ p_{\alpha_j} &\equiv p_{\alpha_j}(\text{detection}|\tau, t_{\text{eff}}, \Theta, \bar{z}), \\ m_{\lambda} &\equiv \text{mag}_{\lambda}(\lambda, \tau, \alpha_j, z_i), \\ m_{\text{lim}} &\equiv \text{mag}_{\text{lim}}(\lambda, t_{\text{eff}}, \Theta), \quad \text{and} \end{align*}$$

$$\text{pr}(z_i) = \exp\left(-\frac{(z_i - \bar{z})^2}{2\sigma_{\bar{z}}^2}\right),$$

where $\tau$ is the time after the GW detection; $\Theta$ contains the specific observation strategy characteristics, including $t_{exp}$, filter, and area coverage as described in Table 3; $m_{lim}$ is the limiting magnitude of the observation for a $10\sigma$ detection; and $\text{pr}(\alpha_j)$ are the model priors defined by each KN model described in Section 4.2 and Table 4. The summation is over all redshifts for a given model, j, at a given $\tau$. Equation (2) represents a Gaussian prior for choosing light-curve models from the grid of SNANA-defined redshifts given the GW distance. Explaining it differently, Equation (2) is, for a given Kasen model, examining whether the resulting apparent magnitude is less than the limiting model with a Gaussian prior on the redshift, using the GW event mean redshift, $\bar{z}$, and the grid of SNANA redshifts the KN could be at, weighted by the GW event variance in redshift, $\sigma_{\bar{z}}^2$.

The GW localization maps present the probability that the event is located at a given sky position, the luminosity distance at that position, and its uncertainty. Therefore, we define the total probability of detecting an event as

$$P_d = \frac{\int_{\hat{\Omega}} d\Omega d_L(\Omega) p(\Omega)}{\int_{\hat{\Omega}} d\Omega d_L(\Omega)} \times \sum_{j=1}^{329} p_{\alpha_j}, \quad (3)$$

where

$$P_d \equiv P(\text{discovery}|\hat{\Omega}, \tau, t_{\text{eff}}, \Theta, \bar{z}),$$

and $\hat{\Omega}$ is the entire sky area observed in the follow-up, $d\Omega$ is the voxel, $p(\Omega)$ is the probability in the voxel, and $d_L(d\Omega)$ is the luminosity distance to the voxel. The sum over $p_{\alpha_j}$ includes the priors and thus indicates the model used. The sky coverage and exposure times determine the total telescope time for a given KN detection. We note that we use weighted spatial probability, rather than 2D on-sky probabilities; therefore, we have the $d_L$ weighting in Equation (3) and the attendant implication that we can have a higher detection probability than on-sky credible area covered. We can now evaluate $P_d$ for a given set of $\Theta$.

## 3.2. Confirmation Probability

In order for an object to be confirmed as a KN candidate, we require it be detected twice (in two observing “epochs”). This requirement can be lifted if, for example, there are sufficient spectroscopic resources available to follow up all of the candidates found after a single epoch. Generally, this is not the case, however. The second detection eliminates spurious detections, including image artifacts, asteroids, and other possible contaminants (Morgan et al. 2020; Shandonay et al. 2022). Another reason to consider the detections independent is that we typically observe while working on postprocessing and making target selection for spectroscopy in the data from the previous epoch.

Given that we want to make two detections to positively identify KN candidates, we define the probability of confirming the transient with two independent detections as

$$P_c \equiv P(\text{confirmed}) = P(\text{discovery}|\hat{\Omega}, \tau, t_{\text{eff}}, \Theta, \bar{z}) \times P(\text{discovery}|\tau_2, t_{\text{eff}}, \Theta, \bar{z}), \quad (4)$$

where

$$P(\text{discovery}|\tau, t_{\text{eff}}, \Theta, \bar{z}) = \sum_{\alpha_j} p_{\alpha_j}.$$

We design the strategies along the paper optimizing for $P_c$ given a set of constraints. Due to implementation choices and to speed up the numerical optimization, we use $P_c^* \equiv P_{d,1} \cdot P_{d,2}$, which gives us equivalent results in terms of strategy,<fnref n="60" /> such as telescope time or duration of the night. We further present our main findings as a function of discovery probability, which is the most relevant outcome for the proposed strategies.

<footnote n="60">⁶⁰ In this form, $P^*(\text{confirmed})$ has an extra spatial/volume probability term multiplied by $V_p$, in which $P_c^* = V_p \cdot P_c$ is the same for both $P_{d1}$ and $P_{d2}$. Therefore, this is irrelevant from an optimization viewpoint, since arg $\max(P_c^*) = \arg \max(P_c)$.</footnote>
## 4. Optimizations and KN Types

After calculating $P_c$ for each O4 event in two observing scenarios, each with a grid of observational parameters, we can evaluate what works best to optically find the KN. The answer to this depends on the science goal. Our primary science objective is standard siren cosmology, so we aim to identify optical counterparts to every KN progenitor detected by the LVK. To this end, we chose to focus on optimizations that require two detections at least 30 minutes apart in order to remove spurious detections due to asteroids. In this section, we discuss optimizing the strategy given the science goal and then discuss the detailed metric, which involves exploring the meaning of covering the space of Kasen models.

### 4.1. Optimizations

We derive $P_c$ for each of the 1280 parameter combinations of scenario 1 and 3520 combinations of scenario 2, resulting in 4800 total for each simulated merger detection event, $S_i$. Not all of these combinations have an appreciable $P_c$ much above zero, as most of the predicted magnitudes are below sky noise.

We use $P_c$ as the variable to optimize on. We choose the highest $P_c$ for each sim $S_i$, look up the set of observational parameters $\Theta_i$ for it, and define this as the top strategy. Choosing the highest $P_c$ is the simplest optimization, but for our evaluation, we need at least two more, reference and low telescope time (low-TT).

1. Top is the $\Theta_i$ producing the highest $P_c$ for each $S_i$, the observational parameters producing the highest probability of confirmation for every O4 simulation.

2. Low-TT is the $\Theta_i$ combination that uses the lowest telescope time given while retaining a $P_c$ within 10% of the highest confirmation probability strategy, by definition, top. For example, if the top strategy finds $P_c = 0.85$, then a low-TT strategy will have $P_c \ge 0.75$, usually with a much reduced telescope time. We will find it interesting to vary the threshold away from 10%.

3. Reference has the 90% probability sky area observed in the i and z bands with 90 s exposures on the first two nights after the merger. This strategy models previous DECam searches, in particular the extensive search of Morgan et al. (2020), and has been used as the DECam strategy for the predictions of Chase et al. (2022).

The top strategy uses as much telescope time as possible to explore the volume, given the parameter exploration presented in Table 3. All strategies work within the constraint of requiring two passes over two 8 hr nights.

### 4.2. Bayesian Average Models

Not every Kasen model atmosphere is equally likely to be a good model for a real KN light curve. Most models for GW170817 are 2+ component models, as in Kilpatrick et al. (2017), Villar et al. (2017), Coughlin et al. (2018), and Gillanders et al. (2022). If the Kasen models define a linear space of KN models, then the Kilpatrick models are in that space; if the Kasen models are eigenmodels of KNe, then the models in Kilpatrick et al. are defined by the eigenvalues multiplying the eigenmodels. Both the values of the nonzero eigenvalues and the number of nonzero eigenvalues are highly model-dependent. Dropping the eigenvector language, it is clear that in the current situation of very few well-studied KNe,

the number of components in models describing KN candidates is uncertain.

One of the most common ways to define detection efficiency in the literature is to set up a grid of KN models, for example, over a viewing angle, as done for a $\chi^2$ analysis, and then calculating the fraction of models detected given an observation. This makes the detection probability explicitly dependent on nonphysical choices of the grid breadth and spacing. In our case, the grid would be the 329 Kasen models, even though these models were meant to extend past the range of models likely to describe real KNe. In a Bayesian framework, each of these models would come with a prior describing our belief in their applicability.

We will employ the useful idea of a Bayesian model average. We evaluate the entire grid of Kasen models, but instead of a uniform weighting, we place a Bayesian prior, $pr(\alpha)$, on each model. The Bayesian average model detection probability is implicit in Equation (2) but can be thought of as

$$ \langle P_d \rangle = \sum_{j=1}^{329} P_d \cdot \text{pr}(\alpha_j), \qquad (5) $$

where $\sum \text{pr}(\alpha_j) = 1$. Here we will use Gaussian priors to produce three Bayesian average models—bright and blue, reddish and slow, and red and faint—as given in Table 4. It is useful to guide the intuition to form the Bayesian average model absolute magnitude,

$$ M_{\lambda}(\lambda, \tau) = \sum_{j=1}^{329} M_{\alpha_j}(\lambda, \tau, \alpha_j) \cdot \text{pr}(\alpha_j), \qquad (6) $$

of which the peak $M_i$ is also given in Table 4 and the light curves shown in Figure 3. This means that the quantity, while illuminating, is incomplete, as the Bayesian formalism is designed to make the uncertainties explicit; the curves are thought of as the median value of a band of light curves weighted by the prior.

The bright and blue model is defined as the means and uncertainties of the blue component model in Kilpatrick et al. (2017) interpreted as Gaussian priors. Since Kilpatrick et al. (2017) gives no uncertainties for the blue model, we assume for $M_{ej}$ a relatively narrow 0.001 and a $v_{ej}$ of 0.01. Where Kilpatrick et al. (2017) has the blue model lanthanide fraction evolving from $\log X_{lan} = -4$ to $-6$ as the opacity falls due to KN atmosphere expansion, we take $\log X_{lan} = -5$ with an uncertainty of 1.0. This model is, on average, blue and reaches a peak luminosity a half-day after the trigger in the g band. It is also, on average, the brightest of the three models and 0.8 mag brighter than GW170817's peak, $M_{g,r,i,z} \approx -15.5$.

The reddish and slow model is defined as in bright and blue, except that we take the $\sigma$ to be 10 times the uncertainties there, to reflect our ignorance of the KN population. This results in a prior allowing the entire $\log X_{lan}$ range of the Kasen models to contribute. Such wide priors make our model, on average, redder and slower to peak than GW170817, though with the same $\langle M_i \rangle$. Most of our results use this model.

The red and faint model is defined as the means and 10 times the uncertainties of the red component model in Kilpatrick et al. (2017) interpreted as Gaussian priors. The predominantly lanthanide-rich Kasen models contribute. This results in a model that is, on average, redder and fainter than GW170817, with a peak $\langle M_i \rangle$ fainter by 1 mag and a $\langle M_r \rangle$ fainter by 1.4 mag.
## 4. Optimizations and KN Types

After calculating $P_c$ for each O4 event in two observing scenarios, each with a grid of observational parameters, we can evaluate what works best to optically find the KN. The answer to this depends on the science goal. Our primary science objective is standard siren cosmology, so we aim to identify optical counterparts to every KN progenitor detected by the LVK. To this end, we chose to focus on optimizations that require two detections at least 30 minutes apart in order to remove spurious detections due to asteroids. In this section, we discuss optimizing the strategy given the science goal and then discuss the detailed metric, which involves exploring the meaning of covering the space of Kasen models.

### 4.1. Optimizations

We derive $P_c$ for each of the 1280 parameter combinations of scenario 1 and 3520 combinations of scenario 2, resulting in 4800 total for each simulated merger detection event, $S_i$. Not all of these combinations have an appreciable $P_c$ much above zero, as most of the predicted magnitudes are below sky noise.

We use $P_c$ as the variable to optimize on. We choose the highest $P_c$ for each sim $S_i$, look up the set of observational parameters $\Theta_i$ for it, and define this as the top strategy. Choosing the highest $P_c$ is the simplest optimization, but for our evaluation, we need at least two more, reference and low telescope time (low-TT).

1. Top is the $\Theta_i$ producing the highest $P_c$ for each $S_i$, the observational parameters producing the highest probability of confirmation for every O4 simulation.

2. Low-TT is the $\Theta_i$ combination that uses the lowest telescope time given while retaining a $P_c$ within 10% of the highest confirmation probability strategy, by definition, top. For example, if the top strategy finds $P_c = 0.85$, then a low-TT strategy will have $P_c \ge 0.75$, usually with a much reduced telescope time. We will find it interesting to vary the threshold away from 10%.

3. Reference has the 90% probability sky area observed in the i and z bands with 90 s exposures on the first two nights after the merger. This strategy models previous DECam searches, in particular the extensive search of Morgan et al. (2020), and has been used as the DECam strategy for the predictions of Chase et al. (2022).

The top strategy uses as much telescope time as possible to explore the volume, given the parameter exploration presented in Table 3. All strategies work within the constraint of requiring two passes over two 8 hr nights.

### 4.2. Bayesian Average Models

Not every Kasen model atmosphere is equally likely to be a good model for a real KN light curve. Most models for GW170817 are 2+ component models, as in Kilpatrick et al. (2017), Villar et al. (2017), Coughlin et al. (2018), and Gillanders et al. (2022). If the Kasen models define a linear space of KN models, then the Kilpatrick models are in that space; if the Kasen models are eigenmodels of KNe, then the models in Kilpatrick et al. are defined by the eigenvalues multiplying the eigenmodels. Both the values of the nonzero eigenvalues and the number of nonzero eigenvalues are highly model-dependent. Dropping the eigenvector language, it is clear that in the current situation of very few well-studied KNe,

the number of components in models describing KN candidates is uncertain.

One of the most common ways to define detection efficiency in the literature is to set up a grid of KN models, for example, over a viewing angle, as done for a $\chi^2$ analysis, and then calculating the fraction of models detected given an observation. This makes the detection probability explicitly dependent on nonphysical choices of the grid breadth and spacing. In our case, the grid would be the 329 Kasen models, even though these models were meant to extend past the range of models likely to describe real KNe. In a Bayesian framework, each of these models would come with a prior describing our belief in their applicability.

We will employ the useful idea of a Bayesian model average. We evaluate the entire grid of Kasen models, but instead of a uniform weighting, we place a Bayesian prior, $pr(\alpha)$, on each model. The Bayesian average model detection probability is implicit in Equation (2) but can be thought of as

$$ \langle P_d \rangle = \sum_{j=1}^{329} P_d \cdot \text{pr}(\alpha_j), \qquad (5) $$

where $\sum \text{pr}(\alpha_j) = 1$. Here we will use Gaussian priors to produce three Bayesian average models—bright and blue, reddish and slow, and red and faint—as given in Table 4. It is useful to guide the intuition to form the Bayesian average model absolute magnitude,

$$ M_{\lambda}(\lambda, \tau) = \sum_{j=1}^{329} M_{\alpha_j}(\lambda, \tau, \alpha_j) \cdot \text{pr}(\alpha_j), \qquad (6) $$

of which the peak $M_i$ is also given in Table 4 and the light curves shown in Figure 3. This means that the quantity, while illuminating, is incomplete, as the Bayesian formalism is designed to make the uncertainties explicit; the curves are thought of as the median value of a band of light curves weighted by the prior.

The bright and blue model is defined as the means and uncertainties of the blue component model in Kilpatrick et al. (2017) interpreted as Gaussian priors. Since Kilpatrick et al. (2017) gives no uncertainties for the blue model, we assume for $M_{ej}$ a relatively narrow 0.001 and a $v_{ej}$ of 0.01. Where Kilpatrick et al. (2017) has the blue model lanthanide fraction evolving from $\log X_{lan} = -4$ to $-6$ as the opacity falls due to KN atmosphere expansion, we take $\log X_{lan} = -5$ with an uncertainty of 1.0. This model is, on average, blue and reaches a peak luminosity a half-day after the trigger in the g band. It is also, on average, the brightest of the three models and 0.8 mag brighter than GW170817's peak, $M_{g,r,i,z} \approx -15.5$.

The reddish and slow model is defined as in bright and blue, except that we take the $\sigma$ to be 10 times the uncertainties there, to reflect our ignorance of the KN population. This results in a prior allowing the entire $\log X_{lan}$ range of the Kasen models to contribute. Such wide priors make our model, on average, redder and slower to peak than GW170817, though with the same $\langle M_i \rangle$. Most of our results use this model.

The red and faint model is defined as the means and 10 times the uncertainties of the red component model in Kilpatrick et al. (2017) interpreted as Gaussian priors. The predominantly lanthanide-rich Kasen models contribute. This results in a model that is, on average, redder and fainter than GW170817, with a peak $\langle M_i \rangle$ fainter by 1 mag and a $\langle M_r \rangle$ fainter by 1.4 mag.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The absolute magnitudes *x* days after the trigger for the weighted average of three sets of priors in the KN parameter space models considered.</figcaption></figure>

This model aligns with what is expected by most models for viewing KNe in the orbital plane, where the lanthanide-poor material is hidden from view.

We will use our three models in various ways to evaluate the optimized strategies. The bright and blue model will be the easiest to detect, as the light curves peak at brighter than −16 in *g*, *r*, *i*, and *z*. The red and faint model will be the hardest, as the absolute magnitude only peaks at brighter than −13 in *g*, barely −14 for *r*, and approaching −15 for the observationally more difficult *i* and *z*. The reddish and slow model is intermediate. In this paper, we often use bright and blue to inform the reader's intuition to low inclination angle KNe. We use red and faint to show the effect of a faint KN on the strategies. It can be argued that red and faint represents the most likely KN from an NS–BH merger, as those events are expected to be redder and with low absolute luminosity compared to NS–NS mergers (Anand et al. 2021). We often use reddish and slow as an intermediate case that shows the behavior of choices in strategy optimization well.

There remains considerable uncertainty in the KN population statistics. While our three models are comparable to those in, e.g., Coughlin et al. (2020), Sagués Carracedo et al. (2021), Petrov et al. (2022), and Zhu et al. (2023), the population study of Setzer et al. (2023) has a two-peaked absolute magnitude distribution at *i*-band absolute magnitudes of −15 and −12. In that study, GW170817 is a 95th percentile event in luminosity, and the 50th percentile luminosity corresponds to our $\langle M_i \rangle$ for red and faint. The Setzer et al. (2023) study is for a random distribution of inclination angles and thus corresponds to the intrinsic KN population, whereas it is known that GW observatories predominantly select merging compact objects with inclination angles near ≈30° (see, e.g., Finn & Chernoff 1993; Nissanke et al. 2010; Schutz 2011). The observed KN population from GW event follow-up will mostly be from the first peak of Setzer et al. (2023). While a KN sample with $M \approx -12$ would be challenging even for a 4 m telescope, the first peak $M \approx -15$ is inside our model ranges.

We take our priors from Kilpatrick et al. (2017) because we wish to emphasize the uncertainties in the models in this study. The uncertainties in the models reported by Villar et al. (2017), for example, are much smaller. For observations in O4, the models and uncertainties in Coughlin et al. (2018) are likely more accurate and constraining for GW170717. However, we prefer to keep a wide variety of possible KNe. We suggest using the two-component light curve red and blue models as the counterparts to our bright and blue and red and faint models, interpreted as low and large inclination angle models, and using the one-component model with $10\sigma$ as the counterpart to our reddish and slow, interpreted as a maximally uncertain KN population model.

### 4.3. Strategies and KN Models

We begin our study of the strategies. As the community routinely uses the probability of detecting a KN once, calling this discovery, we will be showing $P_d \equiv P_{d,1}$ in most of the succeeding plots. Note that the strategies are all optimized on $P_c$, and in our language, “chance of discovery” is to be interpreted as “fraction of models detected given our priors on the space of models.”

Reducing our GW population statistics to only distance, we can describe the NS–NS merger population in Figure 1 as having a median $d_L$ of 150 Mpc, a 75% tile of 225 Mpc, and only 1% at $d_L > 300$ Mpc. We take 200 Mpc as a characteristic distance.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The absolute magnitudes *x* days after the trigger for the weighted average of three sets of priors in the KN parameter space models considered.</figcaption></figure>

This model aligns with what is expected by most models for viewing KNe in the orbital plane, where the lanthanide-poor material is hidden from view.

We will use our three models in various ways to evaluate the optimized strategies. The bright and blue model will be the easiest to detect, as the light curves peak at brighter than −16 in *g*, *r*, *i*, and *z*. The red and faint model will be the hardest, as the absolute magnitude only peaks at brighter than −13 in *g*, barely −14 for *r*, and approaching −15 for the observationally more difficult *i* and *z*. The reddish and slow model is intermediate. In this paper, we often use bright and blue to inform the reader's intuition to low inclination angle KNe. We use red and faint to show the effect of a faint KN on the strategies. It can be argued that red and faint represents the most likely KN from an NS–BH merger, as those events are expected to be redder and with low absolute luminosity compared to NS-NS mergers (Anand et al. 2021). We often use reddish and slow as an intermediate case that shows the behavior of choices in strategy optimization well.

There remains considerable uncertainty in the KN population statistics. While our three models are comparable to those in, e.g., Coughlin et al. (2020), Sagués Carracedo et al. (2021), Petrov et al. (2022), and Zhu et al. (2023), the population study of Setzer et al. (2023) has a two-peaked absolute magnitude distribution at *i*-band absolute magnitudes of −15 and −12. In that study, GW170817 is a 95th percentile event in luminosity, and the 50th percentile luminosity corresponds to our $\langle M_i \rangle$ for red and faint. The Setzer et al. (2023) study is for a random distribution of inclination angles and thus corresponds to the intrinsic KN population, whereas it is known that GW observatories predominantly select merging compact objects with inclination angles near ≈30° (see, e.g., Finn & Chernoff 1993; Nissanke et al. 2010; Schutz 2011). The observed KN population from GW event follow-up will mostly be from the first peak of Setzer et al. (2023). While a KN sample with $M \approx -12$ would be challenging even for a 4 m telescope, the first peak $M \approx -15$ is inside our model ranges.

We take our priors from Kilpatrick et al. (2017) because we wish to emphasize the uncertainties in the models in this study. The uncertainties in the models reported by Villar et al. (2017), for example, are much smaller. For observations in O4, the models and uncertainties in Coughlin et al. (2018) are likely more accurate and constraining for GW170717. However, we prefer to keep a wide variety of possible KNe. We suggest using the two-component light curve red and blue models as the counterparts to our bright and blue and red and faint models, interpreted as low and large inclination angle models, and using the one-component model with $10\sigma$ as the counterpart to our reddish and slow, interpreted as a maximally uncertain KN population model.

## 4.3. Strategies and KN Models

We begin our study of the strategies. As the community routinely uses the probability of detecting a KN once, calling this discovery, we will be showing $P_d \equiv P_{d,1}$ in most of the succeeding plots. Note that the strategies are all optimized on $P_c$, and in our language, “chance of discovery” is to be interpreted as “fraction of models detected given our priors on the space of models.”

Reducing our GW population statistics to only distance, we can describe the NS-NS merger population in Figure 1 as having a median $d_L$ of 150 Mpc, a 75% tile of 225 Mpc, and only 1% at $d_L > 300$ Mpc. We take 200 Mpc as a characteristic distance.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Discovery probability vs. luminosity distance using the low-TT observing strategy compared to the reference strategy for the 611 events with a 90% probability area of <300. Each plot describes a different KN model (top left, bright and blue; top right, reddish and slow; bottom left, red and faint). The shaded region around the median line represents the 68% confidence interval of the scatter among the simulations. We also mark the distance to the upper limit of the LVK O4 projected range, as well as noted GW events for reference. Bottom right: distribution of simulations for the low-TT strategy, binned by discovery probability in our wide prior model, reddish and slow. The marginal histograms are discovery-weighted.</figcaption></figure>

We show in Figure 4 the low-TT strategy discovery probability as a function of distance and compare it to the reference strategy. Starting with the easy case model, bright and blue, we see that the low-TT detection probability has a ceiling at ≈90%. The reference strategy $P_d$ falls after 200 Mpc, whereas the low-TT $P_d$ remains high out to 330 Mpc. The red and faint model is difficult for the reference strategy after about 100 Mpc, whereas the low-TT has $P_d > 50\%$ out to 330 Mpc. The reddish and slow model is intermediate, and the low-TT performs well. The luminosity distance versus area distribution shows the probability-weighted histograms on the margins. Not surprisingly, small spatial localizations are easiest to make identifications for, but less obvious is that the most likely distance for detection is flat, between 125 and 175 Mpc.

The exposure time must be balanced by the sky area to be covered; in our maximum 8 hr night⁻¹, DECam using 1 hr exposures can survey 24 deg² night⁻¹. Figure 1 shows that the number of events with a sky area of less than 24 deg² drops rapidly at distances greater than 200 Mpc.

Turning to the distribution of telescope time required per event in Figure 4, we will explore the motivation and design of the low-TT strategy. For the best (by construction) $P_c$ strategy, top, the mode of the time required is ≈13–15 hr for reddish and slow and red and faint but saturates at 3–5 hr for bright and blue. The top strategy likes to use all the time available over two nights to maximize detection. The low-TT strategy is the lowest telescope time within 10% of top's $P_c$, but clearly one

can tune how much loss in $P_c$ one is willing to accept. We have compared choosing 5%, 10%, and 15% (see Table 5), and not surprisingly, the best choice depends on the KN model chosen. For the bright and blue model, a threshold of 5% gives a strategy that outperforms the reference strategy in $P_d$, in particular after 200 Mpc, while using less total telescope time considering following up 90%–100% of the best-localized events after our initial cut of 300 deg². This is, in fact, one of the key advantages of our proposed method; we scale up the exposure times with the distance, while the reference strategy keeps the same exposure times of 90 s; thus, we keep high performance in farther-away events and use only the necessary time for close events. Another critical point to make optimized use of time is the fact that our proposed method allows the use of shallower exposures in the outer regions that contain less probability while focusing more deep exposures in the core probability region. This configuration, named scenario 2 in Table 3, was preferred, as discussed later in Section 5.1, and is the one presented in the plots unless otherwise stated. This is particularly important in the low-TT strategy, since it selects, most of the time, 60 s exposures in the bright and blue model in the outer regions and therefore saves time compared to the reference strategy. For the reddish and slow and red and faint models, we prefer to set the threshold at 10%. This produces, in the reddish and slow case, a strategy that uses a factor of 3 less telescope time than top and, in the red and faint case, a factor of 2 less, for a loss in $P_c$ of <10%. The gain in $P_d$ over the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Discovery probability vs. luminosity distance using the low-TT observing strategy compared to the reference strategy for the 611 events with a 90% probability area of <300. Each plot describes a different KN model (top left, bright and blue; top right, reddish and slow; bottom left, red and faint). The shaded region around the median line represents the 68% confidence interval of the scatter among the simulations. We also mark the distance to the upper limit of the LVK O4 projected range, as well as noted GW events for reference. Bottom right: distribution of simulations for the low-TT strategy, binned by discovery probability in our wide prior model, reddish and slow. The marginal histograms are discovery-weighted.</figcaption></figure>

We show in Figure 4 the low-TT strategy discovery probability as a function of distance and compare it to the reference strategy. Starting with the easy case model, bright and blue, we see that the low-TT detection probability has a ceiling at ≈90%. The reference strategy $P_d$ falls after 200 Mpc, whereas the low-TT $P_d$ remains high out to 330 Mpc. The red and faint model is difficult for the reference strategy after about 100 Mpc, whereas the low-TT has $P_d > 50\%$ out to 330 Mpc. The reddish and slow model is intermediate, and the low-TT performs well. The luminosity distance versus area distribution shows the probability-weighted histograms on the margins. Not surprisingly, small spatial localizations are easiest to make identifications for, but less obvious is that the most likely distance for detection is flat, between 125 and 175 Mpc.

The exposure time must be balanced by the sky area to be covered; in our maximum 8 hr night$^{-1}$, DECam using 1 hr exposures can survey $24 \text{ deg}^2 \text{ night}^{-1}$. Figure 1 shows that the number of events with a sky area of less than $24 \text{ deg}^2$ drops rapidly at distances greater than 200 Mpc.

Turning to the distribution of telescope time required per event in Figure 4, we will explore the motivation and design of the low-TT strategy. For the best (by construction) $P_c$ strategy, top, the mode of the time required is ≈13–15 hr for reddish and slow and red and faint but saturates at 3–5 hr for bright and blue. The top strategy likes to use all the time available over two nights to maximize detection. The low-TT strategy is the lowest telescope time within 10% of top's $P_c$, but clearly one

can tune how much loss in $P_c$ one is willing to accept. We have compared choosing 5%, 10%, and 15% (see Table 5), and not surprisingly, the best choice depends on the KN model chosen. For the bright and blue model, a threshold of 5% gives a strategy that outperforms the reference strategy in $P_d$, in particular after 200 Mpc, while using less total telescope time considering following up 90%–100% of the best-localized events after our initial cut of 300 deg$^2$. This is, in fact, one of the key advantages of our proposed method; we scale up the exposure times with the distance, while the reference strategy keeps the same exposure times of 90 s; thus, we keep high performance in farther-away events and use only the necessary time for close events. Another critical point to make optimized use of time is the fact that our proposed method allows the use of shallower exposures in the outer regions that contain less probability while focusing more deep exposures in the core probability region. This configuration, named scenario 2 in Table 3, was preferred, as discussed later in Section 5.1, and is the one presented in the plots unless otherwise stated. This is particularly important in the low-TT strategy, since it selects, most of the time, 60 s exposures in the bright and blue model in the outer regions and therefore saves time compared to the reference strategy. For the reddish and slow and red and faint models, we prefer to set the threshold at 10%. This produces, in the reddish and slow case, a strategy that uses a factor of 3 less telescope time than top and, in the red and faint case, a factor of 2 less, for a loss in $P_c$ of <10%. The gain in $P_d$ over the
<table><caption>**Table 5**
Average Telescope Time per Event in Hours Required for Two Detections, Discovery, and Confirmation</caption><thead><tr><th>Strategy</th><th colspan="3">Telescope Time (hr)</th></tr></thead><tbody><tr><td>Bright and blue</td><td>50%</td><td>90%</td><td>100%</td></tr><tr><td>Reference</td><td>0.4</td><td>1.1</td><td>1.6</td></tr><tr><td>Top</td><td>2.3</td><td>3.0</td><td>3.2</td></tr><tr><td>Low-TT (5%)</td><td>0.4</td><td>0.9</td><td>1.3</td></tr><tr><td>Low-TT (10%)</td><td>0.4</td><td>0.9</td><td>1.2</td></tr><tr><td>Low-TT (15%)</td><td>0.2</td><td>0.7</td><td>0.9</td></tr><tr><td colspan="4">Wide prior, reddish and slow</td></tr><tr><td>Reference</td><td>0.4</td><td>1.1</td><td>1.5</td></tr><tr><td>Top</td><td>10.5</td><td>12.1</td><td>12.2</td></tr><tr><td>Low-TT (5%)</td><td>5.2</td><td>6.8</td><td>7.0</td></tr><tr><td>Low-TT (10%)</td><td>3.2</td><td>4.3</td><td>4.5</td></tr><tr><td>Low-TT (15%)</td><td>2.2</td><td>3.0</td><td>3.2</td></tr><tr><td colspan="4">Red and faint</td></tr><tr><td>Reference</td><td>0.4</td><td>1.1</td><td>1.6</td></tr><tr><td>Top</td><td>10.4</td><td>12.3</td><td>12.6</td></tr><tr><td>Low-TT (5%)</td><td>7.4</td><td>9.1</td><td>9.2</td></tr><tr><td>Low-TT (10%)</td><td>5.4</td><td>6.3</td><td>6.4</td></tr><tr><td>Low-TT (15%)</td><td>4.1</td><td>4.8</td><td>4.8</td></tr></tbody></table>

**Note.** The events used were all that had a 90% probability area of <300 deg². Of 860 events, 611 are retained after this cut. The 50%, 90%, and 100% columns give average times for events ordered by the statistics of the size of the 90% localization area, low to high.

reference strategy is particularly dramatic for red and faint at higher distances.

Figure 5 shows how the low-TT strategy selects its optimized choice as compared to the top strategy. All 4800 $P_{d,1}$ for the parameter combinations in Table 3, assuming scenario 2, for a single event simulation are shown, color coded by the required telescope time. The highest $P_{d,1}$ tend to use the highest telescope times, but there are many high $P_{d,1}$ parameter combinations $\Theta_i$ that use much less telescope time. The bottom panel has the top and low-TT strategies marked. Generically, our strategies prefer balanced $P_{d,1}$ and $P_{d,2}$ to maximize $P_c$. For the event presented in Figure 5, the low-TT configuration chosen for the first detection was a core exposure of 200 s (120 s) using the *i* band, and the second detection was done with 300 s (200 s) in the core (outer) region using *z*. The core (outer) region is 80%(90%) during day 1 after the burst. The top configuration prefers to explore the core (outer) region with 600 s (300 s) exposures in the first detection and 300 s (200 s) in the second detection, with both cases using filter *z*. The core (outer) area represents 70% (90%) regions. It is worth noticing that this reddish and slow model peaks in the *i* band around 1 day after the burst and in 1.2 days in the *z* band; however, this model presents a slow decay in *z*, allowing a significant probability of discovery 2 days after the burst.

There is the question of “when is it good enough?,” or diminishing returns. We can adopt the reddish and slow model compared to the reference strategy, for example, and work with the best 90% of events. Then low-TT at a cost of 4× more telescope time detects 20% more KNe, and top at a cost of ~12× more telescope time detects only ~30% more KNe. Whether one is willing to accept the cost depends on the science case. For the standard siren cosmology case, which wants to maximize the number of KNe detected, one would prefer the maximum return of top, but might be willing to accept the rate of low-TT. For the

<figure><img src="image_8.png" /></figure>

**Figure 5.** Discovery probabilities of different $\Theta_i$ sets of configurations on a single event in scenario 2. The event is located at 160 Mpc with a 90% sky area of 168 deg² in a reddish and slow KN model peaking at 0.9 day after the burst in *i* and 1.2 days in *z*. Top: discovery probability vs. time to completion of first pass denoted for each configuration $\Theta_i$. Total telescope time in hours is color coded. Bottom: The best 10% of the confirmation probability from the top panel, displayed in a $P_{d,1}$ vs. $P_{d,2}$ plot with total telescope time color coded. The $\Theta_i$ corresponding to the top and low-TT strategies are marked. The histogram gives the distribution of telescope times for the $\Theta_i$ in the bottom panel.

science case of studying the astrophysics of KNe, where one wants to select good objects for detailed astrophysical study, gathering the next 10 expected in O4 might well be worth expending what top requires for the right events.

Returning to the performance of the strategies, we delve deeper in Figure 6 for the low-TT strategy. In this figure, we present detection probability as a function of sky area. Here the distance weight in Equation (3) becomes important. In each event, the detection probability calculated for each voxel in the sky map is weighted by the distance to that voxel to form a distance-weighted detection probability for the whole event. This accurately weights the probability for the often smaller distances in the lower-probability areas of the sky map (see, e.g., Singer et al. 2016a). Thus, the detection probability for bright and blue at $d_L = 50$ Mpc can be 99% when we only cover the 90% sky probability region.

The easiest way to understand Figure 6 is to start with the right column, especially the bright and blue model. We have seen from Figure 4 that low-TT is very effective at discovery
<table><caption>**Table 5**
Average Telescope Time per Event in Hours Required for Two Detections, Discovery, and Confirmation</caption><thead><tr><th>Strategy</th><th colspan="3">Telescope Time (hr)</th></tr></thead><tbody><tr><td>Bright and blue</td><td>50%</td><td>90%</td><td>100%</td></tr><tr><td>Reference</td><td>0.4</td><td>1.1</td><td>1.6</td></tr><tr><td>Top</td><td>2.3</td><td>3.0</td><td>3.2</td></tr><tr><td>Low-TT (5%)</td><td>0.4</td><td>0.9</td><td>1.3</td></tr><tr><td>Low-TT (10%)</td><td>0.4</td><td>0.9</td><td>1.2</td></tr><tr><td>Low-TT (15%)</td><td>0.2</td><td>0.7</td><td>0.9</td></tr><tr><td colspan="4">Wide prior, reddish and slow</td></tr><tr><td>Reference</td><td>0.4</td><td>1.1</td><td>1.5</td></tr><tr><td>Top</td><td>10.5</td><td>12.1</td><td>12.2</td></tr><tr><td>Low-TT (5%)</td><td>5.2</td><td>6.8</td><td>7.0</td></tr><tr><td>Low-TT (10%)</td><td>3.2</td><td>4.3</td><td>4.5</td></tr><tr><td>Low-TT (15%)</td><td>2.2</td><td>3.0</td><td>3.2</td></tr><tr><td colspan="4">Red and faint</td></tr><tr><td>Reference</td><td>0.4</td><td>1.1</td><td>1.6</td></tr><tr><td>Top</td><td>10.4</td><td>12.3</td><td>12.6</td></tr><tr><td>Low-TT (5%)</td><td>7.4</td><td>9.1</td><td>9.2</td></tr><tr><td>Low-TT (10%)</td><td>5.4</td><td>6.3</td><td>6.4</td></tr><tr><td>Low-TT (15%)</td><td>4.1</td><td>4.8</td><td>4.8</td></tr></tbody></table>

**Note.** The events used were all that had a 90% probability area of <300 deg². Of 860 events, 611 are retained after this cut. The 50%, 90%, and 100% columns give average times for events ordered by the statistics of the size of the 90% localization area, low to high.

reference strategy is particularly dramatic for red and faint at higher distances.

Figure 5 shows how the low-TT strategy selects its optimized choice as compared to the top strategy. All 4800 $P_{d,1}$ for the parameter combinations in Table 3, assuming scenario 2, for a single event simulation are shown, color coded by the required telescope time. The highest $P_{d,1}$ tend to use the highest telescope times, but there are many high $P_{d,1}$ parameter combinations $\Theta_i$ that use much less telescope time. The bottom panel has the top and low-TT strategies marked. Generically, our strategies prefer balanced $P_{d,1}$ and $P_{d,2}$ to maximize $P_c$. For the event presented in Figure 5, the low-TT configuration chosen for the first detection was a core exposure of 200 s (120 s) using the *i* band, and the second detection was done with 300 s (200 s) in the core (outer) region using *z*. The core (outer) region is 80%(90%) during day 1 after the burst. The top configuration prefers to explore the core (outer) region with 600 s (300 s) exposures in the first detection and 300 s (200 s) in the second detection, with both cases using filter *z*. The core (outer) area represents 70% (90%) regions. It is worth noticing that this reddish and slow model peaks in the *i* band around 1 day after the burst and in 1.2 days in the *z* band; however, this model presents a slow decay in *z*, allowing a significant probability of discovery 2 days after the burst.

There is the question of “when is it good enough?,” or diminishing returns. We can adopt the reddish and slow model compared to the reference strategy, for example, and work with the best 90% of events. Then low-TT at a cost of 4× more telescope time detects 20% more KNe, and top at a cost of ~12× more telescope time detects only ~30% more KNe. Whether one is willing to accept the cost depends on the science case. For the standard siren cosmology case, which wants to maximize the number of KNe detected, one would prefer the maximum return of top, but might be willing to accept the rate of low-TT. For the

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Discovery probabilities of different $\Theta_i$ sets of configurations on a single event in scenario 2. The event is located at 160 Mpc with a 90% sky area of 168 deg² in a reddish and slow KN model peaking at 0.9 day after the burst in *i* and 1.2 days in *z*. Top: discovery probability vs. time to completion of first pass denoted for each configuration $\Theta_i$. Total telescope time in hours is color coded. Bottom: The best 10% of the confirmation probability from the top panel, displayed in a $P_{d,1}$ vs. $P_{d,2}$ plot with total telescope time color coded. The $\Theta_i$ corresponding to the top and low-TT strategies are marked. The histogram gives the distribution of telescope times for the $\Theta_i$ in the bottom panel.</figcaption></figure>

science case of studying the astrophysics of KNe, where one wants to select good objects for detailed astrophysical study, gathering the next 10 expected in O4 might well be worth expending what top requires for the right events.

Returning to the performance of the strategies, we delve deeper in Figure 6 for the low-TT strategy. In this figure, we present detection probability as a function of sky area. Here the distance weight in Equation (3) becomes important. In each event, the detection probability calculated for each voxel in the sky map is weighted by the distance to that voxel to form a distance-weighted detection probability for the whole event. This accurately weights the probability for the often smaller distances in the lower-probability areas of the sky map (see, e.g., Singer et al. 2016a). Thus, the detection probability for bright and blue at $d_L = 50$ Mpc can be 99% when we only cover the 90% sky probability region.

The easiest way to understand Figure 6 is to start with the right column, especially the bright and blue model. We have seen from Figure 4 that low-TT is very effective at discovery
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Using the low-TT case. Left column: detection probability for a given luminosity distance and sky area for bright and blue (top), reddish and slow (middle), and red and faint (bottom). Right column: average required telescope time per event area as a function of luminosity distance, ordered the same as in the left column.</figcaption></figure>

for this model, and the two-dimensional projection shown in Figure 6 shows little variation about that high efficiency. What variation is present is the expected loss of efficiency with increasing distance and area. In the right column, for the bright and blue case, it is clear that the telescope time required increases with increasing distance and area. The increase with area is very nearly the ratio of the sky areas.

The red and faint model efficiency plot shows that the dominant variation is a decrease in efficiency with increasing distance as the faint objects fall below the limiting magnitude for the maximum exposure time. Table 3 shows a maximum exposure of 1 hr, so in our 8 hr maximum night, DECam can cover 24 deg². The variation in sky area is not as dramatic as with distance and shows the success of the inner/outer split in scenario 2, when deep exposures over the high-probability area are combined with shallower

exposures over the lower-probability area. (Note that Equation (2) was evaluated separately in each area.) The corresponding telescope time plot behaves as expected, with increasing telescope time with increasing distance and area.

The intermediate reddish and slow model behaves as expected in efficiency. It is more *m*<sub>lim</sub>-dominated than bright and blue and does not require an *m*<sub>lim</sub> as deep as red and faint. The smallest-area bin is likely using 1 hr to maximize area coverage. The corresponding telescope time plot cell shows 5.7 hr for the average event, likely less than 8 hr because of the events with less than 24 deg² sky area. The telescope time plot as a whole shows a surprise. There is a peak in the time at intermediate distances and area, and then the time falls with further distance. One can see some of the same behavior in the top row of the red and faint model telescope time plot. The explanation is that there are strategies (here meaning
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Using the low-TT case. Left column: detection probability for a given luminosity distance and sky area for bright and blue (top), reddish and slow (middle), and red and faint (bottom). Right column: average required telescope time per event area as a function of luminosity distance, ordered the same as in the left column.</figcaption></figure>

for this model, and the two-dimensional projection shown in Figure 6 shows little variation about that high efficiency. What variation is present is the expected loss of efficiency with increasing distance and area. In the right column, for the bright and blue case, it is clear that the telescope time required increases with increasing distance and area. The increase with area is very nearly the ratio of the sky areas.

The red and faint model efficiency plot shows that the dominant variation is a decrease in efficiency with increasing distance as the faint objects fall below the limiting magnitude for the maximum exposure time. Table 3 shows a maximum exposure of 1 hr, so in our 8 hr maximum night, DECam can cover 24 deg². The variation in sky area is not as dramatic as with distance and shows the success of the inner/outer split in scenario 2, when deep exposures over the high-probability area are combined with shallower

exposures over the lower-probability area. (Note that Equation (2) was evaluated separately in each area.) The corresponding telescope time plot behaves as expected, with increasing telescope time with increasing distance and area.

The intermediate reddish and slow model behaves as expected in efficiency. It is more $m_{lim}$-dominated than bright and blue and does not require an $m_{lim}$ as deep as red and faint. The smallest-area bin is likely using 1 hr to maximize area coverage. The corresponding telescope time plot cell shows 5.7 hr for the average event, likely less than 8 hr because of the events with less than 24 deg² sky area. The telescope time plot as a whole shows a surprise. There is a peak in the time at intermediate distances and area, and then the time falls with further distance. One can see some of the same behavior in the top row of the red and faint model telescope time plot. The explanation is that there are strategies (here meaning
$\Theta_t$) that use high telescope time to maximize $P_{d,1}$, but these take so long to cover the area that $P_{d,2}$ is compromised by the fading of the KN, and thus $P_c$ is lowered. A higher $P_c$ is obtained by using a shorter time to cover the area before the object fades, thus maximizing the product of $P_{d,1} \cdot P_{d,2}$.

## 5. Exploring Parameter Space

We have constructed 4800 parameter-set evaluations of DECam KN magnitudes produced by SNANA for 611 Bayestar simulation-detected O4 BNS merger events. We have developed a methodology that uses discovery and confirmation probabilities and a set of optimization rules to produce strategies optimized for the discovery of KNe under certain constraints. We have discussed maximized discovery (top strategy) and minimal telescope time (low-TT) already and will discuss several more strategies in Section 5.2. Here we will show the detailed behavior of the strategies.

We use the reddish and slow (GW1701817-blue, $10 \times \sigma$) KN model for the results in this section unless otherwise stated. Likewise, we will use the low-TT strategy for the results in this section unless otherwise noted. To recap this strategy, we go through all combinations of inner and outer region sizes and exposure times as listed in Table 3 and present the combinations that give the lowest telescope time within 10% of the highest confirmation probability strategy.

### 5.1. Exposure Times

One of the dominant features of the observing parameters described in Table 3 is the splitting of scenarios 1 and 2. In scenario 1, the sky is covered with uniform exposure times. In scenario 2, we allow the splitting of the sky localization area into an outer region and a longer-exposure inner region. Scenario 2 is best thought of as a homogeneous pass with a deeper exposure in the high-probability region. Putting aside the distance weighting, covering a sky localization area sets a ceiling on the discovery probability to the probability contained in that area, and thereafter it is maximizing the limiting magnitude in the sky area. Our optimal selections almost always prefer scenario 2 over scenario 1 as, implicitly or explicitly, the total telescope time is constrained.

Figure 7 shows the percentage of simulations that preferred each exposure time broken into inner and outer areas. The mode of the 2D distribution for the inner region is 20 minute exposures, with 31% of events using that exposure time and 84% of events using 20 minutes or less. For the outer region, the mode is 5 minutes, with 33% of events using that exposure time and 74% of events using 5 minutes or less. This is model-dependent. The bright and blue model uses 90 s 99% of the time in the inner region and 60 s 99% of the time in the outer region. Most often, 87% of the time, this is a very small inner region of 0.3, and 93% of the simulations chose the 0.9 outer probability area. This is the lowest exposure time combination for scenario 2 in Table 3, and it is likely that the strategies would have used a shorter exposure time if available, although it is notable that the strategy did not prefer the available 60 s homogeneous pass in scenario 1. The red and faint model has a more complicated exposure time pattern for the inner area. In order of use, 2400, 600, 3600, and 1200 s exposures are used in 26%, 18%, 16%, and 15% of the simulations, respectively. The outer area exposure time has a mode at 1200 s of 25%, with 13% and 14% for 2400 and 600 s, respectively, and 17% of simulations use 300 s.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Top: discovery probability vs. distance for low-TT comparing scenarios 1 and 2, as well as reference, which only uses scenario 1. The middle and bottom panels show the exposure times that were most common when using scenario 2, which breaks up the area into a shallower outer region and a deeper inner region of the spatial sky area.</figcaption></figure>

For distances out to 125 Mpc, the difference between scenarios 1 and 2 is minimal, but at around 300 Mpc, using scenario 2 gains ~5%–10% in discovery probability. In other words, the slope in probability versus distance is shallower for the two-zone scenario than it is for the one-zone scenario. Therefore, we believe that using the deeper exposures in the core region will in general be more successful in the more distant events.
$\Theta_t$) that use high telescope time to maximize $P_{d,1}$, but these take so long to cover the area that $P_{d,2}$ is compromised by the fading of the KN, and thus $P_c$ is lowered. A higher $P_c$ is obtained by using a shorter time to cover the area before the object fades, thus maximizing the product of $P_{d,1} \cdot P_{d,2}$.

## 5. Exploring Parameter Space

We have constructed 4800 parameter-set evaluations of DECam KN magnitudes produced by SNANA for 611 Bayestar simulation-detected O4 BNS merger events. We have developed a methodology that uses discovery and confirmation probabilities and a set of optimization rules to produce strategies optimized for the discovery of KNe under certain constraints. We have discussed maximized discovery (top strategy) and minimal telescope time (low-TT) already and will discuss several more strategies in Section 5.2. Here we will show the detailed behavior of the strategies.

We use the reddish and slow (GW1701817-blue, $10 \times \sigma$) KN model for the results in this section unless otherwise stated. Likewise, we will use the low-TT strategy for the results in this section unless otherwise noted. To recap this strategy, we go through all combinations of inner and outer region sizes and exposure times as listed in Table 3 and present the combinations that give the lowest telescope time within 10% of the highest confirmation probability strategy.

### 5.1. Exposure Times

One of the dominant features of the observing parameters described in Table 3 is the splitting of scenarios 1 and 2. In scenario 1, the sky is covered with uniform exposure times. In scenario 2, we allow the splitting of the sky localization area into an outer region and a longer-exposure inner region. Scenario 2 is best thought of as a homogeneous pass with a deeper exposure in the high-probability region. Putting aside the distance weighting, covering a sky localization area sets a ceiling on the discovery probability to the probability contained in that area, and thereafter it is maximizing the limiting magnitude in the sky area. Our optimal selections almost always prefer scenario 2 over scenario 1 as, implicitly or explicitly, the total telescope time is constrained.

Figure 7 shows the percentage of simulations that preferred each exposure time broken into inner and outer areas. The mode of the 2D distribution for the inner region is 20 minute exposures, with 31% of events using that exposure time and 84% of events using 20 minutes or less. For the outer region, the mode is 5 minutes, with 33% of events using that exposure time and 74% of events using 5 minutes or less. This is model-dependent. The bright and blue model uses 90 s 99% of the time in the inner region and 60 s 99% of the time in the outer region. Most often, 87% of the time, this is a very small inner region of 0.3, and 93% of the simulations chose the 0.9 outer probability area. This is the lowest exposure time combination for scenario 2 in Table 3, and it is likely that the strategies would have used a shorter exposure time if available, although it is notable that the strategy did not prefer the available 60 s homogeneous pass in scenario 1. The red and faint model has a more complicated exposure time pattern for the inner area. In order of use, 2400, 600, 3600, and 1200 s exposures are used in 26%, 18%, 16%, and 15% of the simulations, respectively. The outer area exposure time has a mode at 1200 s of 25%, with 13% and 14% for 2400 and 600 s, respectively, and 17% of simulations use 300 s.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Top: discovery probability vs. distance for low-TT comparing scenarios 1 and 2, as well as reference, which only uses scenario 1. The middle and bottom panels show the exposure times that were most common when using scenario 2, which breaks up the area into a shallower outer region and a deeper inner region of the spatial sky area.</figcaption></figure>

For distances out to 125 Mpc, the difference between scenarios 1 and 2 is minimal, but at around 300 Mpc, using scenario 2 gains ~5%–10% in discovery probability. In other words, the slope in probability versus distance is shallower for the two-zone scenario than it is for the one-zone scenario. Therefore, we believe that using the deeper exposures in the core region will in general be more successful in the more distant events.
## 5.2. Filter Choice and Bright/Dark Nights

In following up an LVK event, there is a similar chance of observing in dark conditions as in bright conditions. In Figure 8, we show the effect of dark versus bright time. The bright/dark distinction is handled in our methodology via a change in $t_{\text{eff}}$, as seen in Table 1. Bright time lowers the discovery probability by 5%–10%, with the loss being mitigated by filter choice and exposure times. The filters used in the two passes are, in dark time, rr (33%), zz (50%), and $ii+iz$ (17%). In bright time, the filters are zz (71%) and $ii+iz$ (25%). The filter choices are strategy-dependent, and, for comparison, the top strategy used in dark time rr (12%), ri (7%), rz (27%), iz (29%), and zz (17%) and in bright time iz (45%) and zz (41%). The filter choices are also model-dependent, driven by the color evolution of the reddish and slow model as seen in Figure 3. The g filter is never going to be favored in this model, and *i* will be picked the first night and *z* on the second, except that it is easier to go deeper in *r* than in *i*. We can predict that the filter selection for the red and faint model will be nearly the same but that the bright and blue model will predominantly use *g* and *r*. Notably, it is not straightforward for us to say which filters we use in our best strategy.

## 5.3. Impact of the Two-detection Requirement

Detecting a counterpart twice in DECam images is an important step in our experience, as it allows one to distinguish extragalactic transients from asteroids or other objects within a short period of time. This requirement is common; for example, Zhu et al. (2021), Sagués Carracedo et al. (2021), and Petrov et al. (2022) all demand two detections for a confirmation of a KN transient. Also important is covering the sky area with multiple filters to distinguish the KN from other transients.

The DECam search does not happen in isolation, and we emphasize the importance of the broader GW community during these times. While our programs’ main line is to proceed to confirmation via a second image, smaller telescopes may chose to follow up the first image preliminary candidates in order to efficiently reduce the candidate list that will be sent to expensive spectroscopic efforts. There are reasons to be interested in both $P_{d,1}$ and $P_c$.

In Figure 9, we show $P$(discovery) in the first detection, $P_{d,1}$, compared to the second detection, $P_{d,2}$. In the low-TT strategy, the detection probability remains nearly constant between pass 1 and pass 2. This is not the case for the reference strategy. For the reddish and slow model, there is a nearly constant offset between $P_{d,1}$ and $P_{d,2}$ favoring the second pass within a few percent. It is worth noticing that this is our model that peaks at later times; therefore, the second pass might happen closest to the peak. However, the difference is inside the 68% confidence interval as presented in Figure 4. For the red and faint model, $P_{d,1}$ is higher at $d < 100$ Mpc, and $P_{d,2}$ is higher at $d > 100$ Mpc. For the bright and blue model, we see the reference strategy become less efficient than the low-TT strategy at $d > 220$ Mpc, as both $P_{d,1}$ and $P_{d,2}$ drop with increasing distance. We infer that the low-TT strategy strongly prefers to balance $P_{d,1}$ and $P_{d,2}$.

## 5.4. Other Strategy Options

Our choice of optimization has flexibility. We might place a high priority on the earliest possible discovery, or we may have lost several nights due to weather conditions and need to find an optimal approach for the first clear night several days after trigger. Let us, therefore, explore three other optimizations.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Top: discovery probability for the different strategies adopted in bright and dark time. Middle: heat map showing the percentage of simulations using each filter in the first and second passes in dark time for the low-TT 10% strategy and the reddish and slow model. Bottom: heat map showing the percentage of simulations using each filter in the first and second passes in bright time for the low-TT 10% strategy and the reddish and slow model.</figcaption></figure>

1. Early discovery (ED) is the $\Theta_i$ that produces the earliest confirmation limited by the $P_c$ from low-TT (5%) for each $S_i$.

2. Late discovery (LD) is the $\Theta_i$ that produces the latest confirmation limited by the $P_c$ from low-TT (5%) for each $S_i$. This family is intended to find a competitive
## 5.2. Filter Choice and Bright/Dark Nights

In following up an LVK event, there is a similar chance of observing in dark conditions as in bright conditions. In Figure 8, we show the effect of dark versus bright time. The bright/dark distinction is handled in our methodology via a change in $t_{\text{eff}}$, as seen in Table 1. Bright time lowers the discovery probability by 5%–10%, with the loss being mitigated by filter choice and exposure times. The filters used in the two passes are, in dark time, rr (33%), zz (50%), and $ii + iz$ (17%). In bright time, the filters are zz (71%) and $ii + iz$ (25%). The filter choices are strategy-dependent, and, for comparison, the top strategy used in dark time rr (12%), ri (7%), rz (27%), iz (29%), and zz (17%) and in bright time iz (45%) and zz (41%). The filter choices are also model-dependent, driven by the color evolution of the reddish and slow model as seen in Figure 3. The g filter is never going to be favored in this model, and *i* will be picked the first night and *z* on the second, except that it is easier to go deeper in *r* than in *i*. We can predict that the filter selection for the red and faint model will be nearly the same but that the bright and blue model will predominantly use *g* and *r*. Notably, it is not straightforward for us to say which filters we use in our best strategy.

## 5.3. Impact of the Two-detection Requirement

Detecting a counterpart twice in DECam images is an important step in our experience, as it allows one to distinguish extragalactic transients from asteroids or other objects within a short period of time. This requirement is common; for example, Zhu et al. (2021), Sagués Carracedo et al. (2021), and Petrov et al. (2022) all demand two detections for a confirmation of a KN transient. Also important is covering the sky area with multiple filters to distinguish the KN from other transients.

The DECam search does not happen in isolation, and we emphasize the importance of the broader GW community during these times. While our programs’ main line is to proceed to confirmation via a second image, smaller telescopes may choose to follow up the first image preliminary candidates in order to efficiently reduce the candidate list that will be sent to expensive spectroscopic efforts. There are reasons to be interested in both $P_{d,1}$ and $P_c$.

In Figure 9, we show $P$(discovery) in the first detection, $P_{d,1}$, compared to the second detection, $P_{d,2}$. In the low-TT strategy, the detection probability remains nearly constant between pass 1 and pass 2. This is not the case for the reference strategy. For the reddish and slow model, there is a nearly constant offset between $P_{d,1}$ and $P_{d,2}$ favoring the second pass within a few percent. It is worth noticing that this is our model that peaks at later times; therefore, the second pass might happen closest to the peak. However, the difference is inside the 68% confidence interval as presented in Figure 4. For the red and faint model, $P_{d,1}$ is higher at $d < 100$ Mpc, and $P_{d,2}$ is higher at $d > 100$ Mpc. For the bright and blue model, we see the reference strategy become less efficient than the low-TT strategy at $d > 220$ Mpc, as both $P_{d,1}$ and $P_{d,2}$ drop with increasing distance. We infer that the low-TT strategy strongly prefers to balance $P_{d,1}$ and $P_{d,2}$.

## 5.4. Other Strategy Options

Our choice of optimization has flexibility. We might place a high priority on the earliest possible discovery, or we may have lost several nights due to weather conditions and need to find an optimal approach for the first clear night several days after trigger. Let us, therefore, explore three other optimizations.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Top: discovery probability for the different strategies adopted in bright and dark time. Middle: heat map showing the percentage of simulations using each filter in the first and second passes in dark time for the low-TT 10% strategy and the reddish and slow model. Bottom: heat map showing the percentage of simulations using each filter in the first and second passes in bright time for the low-TT 10% strategy and the reddish and slow model.</figcaption></figure>

1. Early discovery (ED) is the $\Theta_i$ that produces the earliest confirmation limited by the $P_c$ from low-TT (5%) for each $S_i$.

2. Late discovery (LD) is the $\Theta_i$ that produces the latest confirmation limited by the $P_c$ from low-TT (5%) for each $S_i$. This family is intended to find a competitive
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Discovery probability vs. distance to event for the low-TT observing strategy in the first and second passes compared to the reference strategy for KNe described by models bright and blue (top), reddish and slow with wide priors (middle), and red and faint (bottom).</figcaption></figure>

strategy when one cannot observe during some early/intermediate nights or in case the event was not confirmed in the first days.

3. Half-nights (HN) are intended to find a competitive strategy when one cannot observe half the night—if the object rises or sets, for example, or the telescope is only allocated for half-nights. Thus, from the subset limited by $P_c$ from low-TT (5%) for each $S_i$ configuration, we constrained the strategy to have both passes in less than 4 hr if they are in the same night, or each of the passes takes 4 hr individually if they are in different nights. In this exercise, we consider that the observation starts in the first half of the night.

All of these configurations are restricted in telescope time. In Figure 10, we show the distributions of the time of second pass completion relative to the merger. It is in the second pass that we achieve a confirmation. The ED scenario, not surprisingly, has earlier times than the top strategy, and, by design, earlier times than LD. The LD strategy has a different optimization and thus a different use. If the first night or two is not useful for observations, then the LD strategy is useful for pursuing the discovery at late times. It peaks around 2.5 days after the merger. Note that this strategy is optimized on $P_c$, so it does not describe the case where the event is unobservable due to clouds for a night or two but rather is working the scenario where the object is detected but cannot be confirmed for several days. Figure 10 suggests that since the first pass $P_d$ is higher on the first night, the most likely night to capture the confirmation pass is the third night. LD is representative of the kinds of strategies that would be necessary to deal with weather. The HN strategy enforces a limitation on the amount of time spent in pass 1, and the resulting performance is similar to low-TT in terms of both $P_c$ and telescope time expenditures, although without the guarantee of being within 10% of the best strategy: top.

Figure 11 shows the discovery probability versus distance for the strategies discussed here. For the reddish and slow model, all of our strategies have roughly equal performance (but see the top strategy outperforming low-TT by the amount it is allowed to), although with different telescope time cost. This is true to $d = 180$ Mpc, but after that, the LD strategy becomes less efficient. For the red and faint model, we see, interestingly, that the low-TT and top strategies have equal performance. Finally, the HN strategy performs very well until $d = 180$ Mpc, then it becomes infeasible.

For the reddish and slow model, the ED strategy for exposure time distribution is more or less equally split among the possible outer/inner exposure time pairs in Table 3 up to the outer region exposure time of 600 s, with the most likely outer exposure time of 300 s. There are no outer exposure times greater than 600 s or inner exposure times greater than 2400 s. The LD strategy's most likely exposure time is 1200 s, and the strategy tries to cover the largest inner region possible. The outer region has the most likely exposure time of 300 s and is never longer than 600 s. The strategy is to cover the largest inner region possible with relatively shallow outer region exposures. The second pass exposure time is weighted deeper than the first pass, but the most likely exposure time remains 1200 s. The HN strategy for exposure time distribution is 31% of simulations using 1200 s and 17% using 2400, most often in a 0.7 inner core region. The outer region uses 300 and 600 s exposures over half the time.

For the red and faint model, the ED strategy for exposure time distribution is deeper than for the same in reddish and slow, preferring 600, 1200, and 2400 s instead of 600 and
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Discovery probability vs. distance to event for the low-TT observing strategy in the first and second passes compared to the reference strategy for KNe described by models bright and blue (top), reddish and slow with wide priors (middle), and red and faint (bottom).</figcaption></figure>

strategy when one cannot observe during some early/intermediate nights or in case the event was not confirmed in the first days.

3. Half-nights (HN) are intended to find a competitive strategy when one cannot observe half the night—if the object rises or sets, for example, or the telescope is only allocated for half-nights. Thus, from the subset limited by $P_c$ from low-TT (5%) for each $S_i$ configuration, we constrained the strategy to have both passes in less than 4 hr if they are in the same night, or each of the passes takes 4 hr individually if they are in different nights. In this exercise, we consider that the observation starts in the first half of the night.

All of these configurations are restricted in telescope time. In Figure 10, we show the distributions of the time of second pass completion relative to the merger. It is in the second pass that we achieve a confirmation. The ED scenario, not surprisingly, has earlier times than the top strategy, and, by design, earlier times than LD. The LD strategy has a different optimization and thus a different use. If the first night or two is not useful for observations, then the LD strategy is useful for pursuing the discovery at late times. It peaks around 2.5 days after the merger. Note that this strategy is optimized on $P_c$, so it does not describe the case where the event is unobservable due to clouds for a night or two but rather is working the scenario where the object is detected but cannot be confirmed for several days. Figure 10 suggests that since the first pass $P_d$ is higher on the first night, the most likely night to capture the confirmation pass is the third night. LD is representative of the kinds of strategies that would be necessary to deal with weather. The HN strategy enforces a limitation on the amount of time spent in pass 1, and the resulting performance is similar to low-TT in terms of both $P_c$ and telescope time expenditures, although without the guarantee of being within 10% of the best strategy: top.

Figure 11 shows the discovery probability versus distance for the strategies discussed here. For the reddish and slow model, all of our strategies have roughly equal performance (but see the top strategy outperforming low-TT by the amount it is allowed to), although with different telescope time cost. This is true to $d = 180$ Mpc, but after that, the LD strategy becomes less efficient. For the red and faint model, we see, interestingly, that the low-TT and top strategies have equal performance. Finally, the HN strategy performs very well until $d = 180$ Mpc, then it becomes infeasible.

For the reddish and slow model, the ED strategy for exposure time distribution is more or less equally split among the possible outer/inner exposure time pairs in Table 3 up to the outer region exposure time of 600 s, with the most likely outer exposure time of 300 s. There are no outer exposure times greater than 600 s or inner exposure times greater than 2400 s. The LD strategy's most likely exposure time is 1200 s, and the strategy tries to cover the largest inner region possible. The outer region has the most likely exposure time of 300 s and is never longer than 600 s. The strategy is to cover the largest inner region possible with relatively shallow outer region exposures. The second pass exposure time is weighted deeper than the first pass, but the most likely exposure time remains 1200 s. The HN strategy for exposure time distribution is 31% of simulations using 1200 s and 17% using 2400, most often in a 0.7 inner core region. The outer region uses 300 and 600 s exposures over half the time.

For the red and faint model, the ED strategy for exposure time distribution is deeper than for the same in reddish and slow, preferring 600, 1200, and 2400 s instead of 600 and
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Distribution of confirmation days using the reddish and slow model. The top panel depicts the time of confirmation probability (i.e., how many days it took to observe the area twice) for our ED strategy, and the bottom panel depicts the LD strategy. Most simulated events are confirmed by the first day. For reference, we also show the distribution using the top strategy, which has no restrictions on when to perform the follow-up.</figcaption></figure>

1200 s, never using 90 s exposure times for the inner core like 15% of the reddish and slow simulations do. The LD strategy simulations 57% of the time use 600, 1200, or 2400 s exposures but a wide range of inner areas. The outer region uses all exposure times from 60 to 2400 s but is most heavily weighted toward 300, 600, and especially 1200 s. The HN strategy exposure time distribution is complicated, using a wide variety of inner areas and exposure times skewing deep at 1200 and 2400 s exposures. The outer region is similar, although here, as most often throughout the strategies, the outer region coverage of 0.9 is preferred. Table 6 summarizes the average telescope time required on each strategy and KN model.

## 6. Real Observations

We have presented a variety of strategies optimized for a variety of purposes. Here we describe how to use them when an observing team receives an LVK alert.

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 11.** We illustrate the flexibility of the strategy families. We present the curves constrained with 5% of the top strategy, e.g., low-TT (5%) and derived strategies. Upper: reddish and slow. Lower: red and faint. HN are not always available. For reddish and slow, 395 HN strategies produced detections, while for red and faint, 137 HN strategies produced detections. As the exposure times skew deep, the success rate is likely anticorrelated with the sky area.</figcaption></figure>

We can compute $P(\text{confirmed})$ for the grid of $\Theta$ in Table 3 (both scenarios) using estimations of $t_{\text{eff}}$ for the upcoming night. Given the $\Theta_i$, we can choose a strategy to follow, placing the top, low-TT, ED, or HN optimizations as appropriate for the event. The optimization gives us the single $\Theta_i$ for each of the strategy families.

The strategy computation takes about 30 minutes on a single core. This includes the four strategy families (top, low-TT, ED, and HN), as well as considerations about bright/dark time. The time to complete the computation can be brought down to <1 s by using the simulations as an approximation for the real event. Here our approach would be to choose a $\Theta_i$ by a nearest-neighbor search or to build a simple neural net on the simulation parameters and $P(\text{confirmed})$ to choose $\Theta_i$.

The choice of the KN model is important. For example, if we are very early on-sky, with great observational conditions for just a couple of hours, or some limitation on telescope access is imposed upon us, we might also consider a fast detection of a blue flash (bright and blue model). This strategy could be
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Distribution of confirmation days using the reddish and slow model. The top panel depicts the time of confirmation probability (i.e., how many days it took to observe the area twice) for our ED strategy, and the bottom panel depicts the LD strategy. Most simulated events are confirmed by the first day. For reference, we also show the distribution using the top strategy, which has no restrictions on when to perform the follow-up.</figcaption></figure>

1200 s, never using 90 s exposure times for the inner core like 15% of the reddish and slow simulations do. The LD strategy simulations 57% of the time use 600, 1200, or 2400 s exposures but a wide range of inner areas. The outer region uses all exposure times from 60 to 2400 s but is most heavily weighted toward 300, 600, and especially 1200 s. The HN strategy exposure time distribution is complicated, using a wide variety of inner areas and exposure times skewing deep at 1200 and 2400 s exposures. The outer region is similar, although here, as most often throughout the strategies, the outer region coverage of 0.9 is preferred. Table 6 summarizes the average telescope time required on each strategy and KN model.

## 6. Real Observations

We have presented a variety of strategies optimized for a variety of purposes. Here we describe how to use them when an observing team receives an LVK alert.

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 11.** We illustrate the flexibility of the strategy families. We present the curves constrained with 5% of the top strategy, e.g., low-TT (5%) and derived strategies. Upper: reddish and slow. Lower: red and faint. HN are not always available. For reddish and slow, 395 HN strategies produced detections, while for red and faint, 137 HN strategies produced detections. As the exposure times skew deep, the success rate is likely anticorrelated with the sky area.</figcaption></figure>

We can compute $P(\text{confirmed})$ for the grid of $\Theta$ in Table 3 (both scenarios) using estimations of $t_{\text{eff}}$ for the upcoming night. Given the $\Theta_i$, we can choose a strategy to follow, placing the top, low-TT, ED, or HN optimizations as appropriate for the event. The optimization gives us the single $\Theta_i$ for each of the strategy families.

The strategy computation takes about 30 minutes on a single core. This includes the four strategy families (top, low-TT, ED, and HN), as well as considerations about bright/dark time. The time to complete the computation can be brought down to <1 s by using the simulations as an approximation for the real event. Here our approach would be to choose a $\Theta_i$ by a nearest-neighbor search or to build a simple neural net on the simulation parameters and $P(\text{confirmed})$ to choose $\Theta_i$.

The choice of the KN model is important. For example, if we are very early on-sky, with great observational conditions for just a couple of hours, or some limitation on telescope access is imposed upon us, we might also consider a fast detection of a blue flash (bright and blue model). This strategy could be
<table><caption>**Table 6**
Average Telescope Time per Event in Hours Required for Two Detections, Discovery, and Confirmation, as in Table 5</caption><thead><tr><th>Strategy</th><th colspan="3">Telescope Time (hr)</th></tr></thead><tbody><tr><td>Bright and blue</td><td>50%</td><td>90%</td><td>100%</td></tr><tr><td>ED</td><td>0.4</td><td>1.0</td><td>1.4</td></tr><tr><td>LD<sup>a</sup></td><td>0.4</td><td>1.0</td><td>1.4</td></tr><tr><td>HN<sup>a</sup></td><td>0.4</td><td>1.0</td><td>1.3</td></tr><tr><td>Low-TT—bright night</td><td>0.5</td><td>1.2</td><td>1.6</td></tr><tr><td>Reddish and slow</td><td></td><td></td><td></td></tr><tr><td>ED</td><td>3.3</td><td>4.3</td><td>4.4</td></tr><tr><td>LD<sup>a</sup></td><td>3.1</td><td>3.9</td><td>4.0</td></tr><tr><td>HN<sup>a</sup></td><td>2.9</td><td>4.0</td><td>4.0</td></tr><tr><td>Low-TT—bright night</td><td>3.2</td><td>4.7</td><td>5.3</td></tr><tr><td>Red and faint</td><td></td><td></td><td></td></tr><tr><td>ED</td><td>6.7</td><td>7.7</td><td>7.8</td></tr><tr><td>LD<sup>a</sup></td><td>6.6</td><td>9.1</td><td>9.1</td></tr><tr><td>HN<sup>a</sup></td><td>2.6</td><td>3.9</td><td>4.0</td></tr><tr><td>Low-TT—bright night</td><td>7.6</td><td>9.6</td><td>9.9</td></tr></tbody></table>

**Note.**

<sup>a</sup> Due to the strong constraints in HN and LD, they are not always available for a given GW event simulation. For the bright and blue, reddish and slow, and red and faint models, these strategies are only defined for 602, 496, and 131 simulations, respectively, in the HN scenario, while LD is defined for 532, 511, and 251.

interesting, in particular for low-TT, due to the cost of deep exposures over a wide area.

There is also a flow of decision-making external to what we have described. Whether or not there is a short gamma-ray burst (GRB), we choose distinct approaches. In the case where there is, we are likely looking for a GW170817-like event (bright and blue). If not, then a conservative model is indicated (reddish and slow), as we might not be looking for an event with a high inclination angle. If the alert indicates it is an NS-BH merger, then the KN is likely more consistent with the red and faint model. After this, if there is a possible HN strategy, it is low-budget and high-performing by definition. It also gives enough time for us to use spectroscopy to quickly confirm. In an observing run where we expect one BNS merger every month, we could plan on allocating telescope time through the run considering the amount of telescope time remaining versus where the particular event lies in an S/N distribution (Chen & Holz 2014) or among a population of simulated events. This would lead to spending more toward the end of the run if considerable time was left or spending more if it was a particularly good event compared to the simulations.

## 6.1. Plan for Usage

To demonstrate how using this code might work during a real GW follow-up campaign, we randomly select one of the simulations used in this analysis. For this exercise, we are assuming that the LVK trigger announcement was received toward the beginning of the night, meaning we will aim to be on-sky within a few hours of merger. The general logic flow for determining the observing plan goes as follows.

1. Assess our external factors.

(a) Find out how many hours we are allotted for the night and about telescope availability.

(b) Find out the sky conditions.

(c) Check if there was a GRB reported in the area at the same time (as was the case with GW170817). This is important for choosing a KN model to use.

2. Using the external factors as input, run the code.

3. Assess output csv; see Table 7. Here we would decide if this is an event that we would like to follow up, considering the discovery probability.

4. Choose filters, exposure times, and area coverage from strategy and compile an observing plan.

While most external factors are easy to identify, choosing which KN model to follow may not be obvious and will depend on the particular science goal of the project. For instance, if one's science goal is to have the most complete set of KNe, then using the faintest KN model is useful; thus, red and faint is the favored model. For simplicity, in this example, we will assume that a GRB was reported within the LVK sky area, thus favoring bright and blue. Once the code has generated an output csv file, we can compare strategies.

The configuration with the highest detection probability for each strategy is displayed in Table 7. Here we can see that for this event, the top strategy yields the highest probability of detection, as expected. In this example, however, the top strategy has a factor of almost 5 more in the amount of telescope time needed for any other strategy for only an ~7% greater probability of detection.

## 7. Discussion

In this paper, our science goal is to maximize the number of KNe with GW-measured $d_L$ and securely identified redshifts. We put a high weight on the completeness of detection of the KNe, given the considerable uncertainty in the KN population. Our strategies go deep. The strategy chosen by our optimization considers a constraint in telescope time, i.e., in low-TT mode for a reddish and slow KN reach $m_{lim}(10\sigma)$ of $r \ge 24.4$ and $z \ge 22.9$ for the $\ge 1200$ s exposures (see Table 2) for the inner regions of 47% (as defined in Table 3) from our set of 611 simulations, regardless of distance (see Figure 7).

We can gain insight into our results by comparing how the literature handles a set of three questions:

1. the modeling of GW merger event distance and sky area distribution;

2. the range of KN model physical parameters, including inclination angle; and

3. the telescope, search cadence, detectability versus distance and $m_{lim}$.

*The GW event properties.* We simulate a merging NS-NS population expected in O4 given the expected LIGO sensitivities, drawing from the NS population, simulating the GW waveform, projecting onto a GW observatory network, and converting the GW observations to sky maps. Our approach to drawing, simulating, projecting, and forming sky maps has been performed before (Chen et al. 2021; Petrov et al. 2022), though it is sometimes done without sky maps (Zhu et al. 2021). There are simplified approaches. One can just assume a KN population to observe (Sagués Carracedo et al. 2021; Chase et al. 2022; Setzer et al. 2023). Or, one can use toy GW models, such as a low-significance event sky map at 200 Mpc (Coughlin et al. 2020). All of these approaches assume that the KN properties are independent of the GW event properties, other than distance. Colombo et al. (2022)
<table><caption>**Table 6**
Average Telescope Time per Event in Hours Required for Two Detections, Discovery, and Confirmation, as in Table 5</caption><thead><tr><th>Strategy</th><th colspan="3">Telescope Time (hr)</th></tr></thead><tbody><tr><td>Bright and blue</td><td>50%</td><td>90%</td><td>100%</td></tr><tr><td>ED</td><td>0.4</td><td>1.0</td><td>1.4</td></tr><tr><td>LD<sup>a</sup></td><td>0.4</td><td>1.0</td><td>1.4</td></tr><tr><td>HN<sup>a</sup></td><td>0.4</td><td>1.0</td><td>1.3</td></tr><tr><td>Low-TT—bright night</td><td>0.5</td><td>1.2</td><td>1.6</td></tr><tr><td>Reddish and slow</td><td></td><td></td><td></td></tr><tr><td>ED</td><td>3.3</td><td>4.3</td><td>4.4</td></tr><tr><td>LD<sup>a</sup></td><td>3.1</td><td>3.9</td><td>4.0</td></tr><tr><td>HN<sup>a</sup></td><td>2.9</td><td>4.0</td><td>4.0</td></tr><tr><td>Low-TT—bright night</td><td>3.2</td><td>4.7</td><td>5.3</td></tr><tr><td>Red and faint</td><td></td><td></td><td></td></tr><tr><td>ED</td><td>6.7</td><td>7.7</td><td>7.8</td></tr><tr><td>LD<sup>a</sup></td><td>6.6</td><td>9.1</td><td>9.1</td></tr><tr><td>HN<sup>a</sup></td><td>2.6</td><td>3.9</td><td>4.0</td></tr><tr><td>Low-TT—bright night</td><td>7.6</td><td>9.6</td><td>9.9</td></tr></tbody></table>

**Note.**

<sup>a</sup> Due to the strong constraints in HN and LD, they are not always available for a given GW event simulation. For the bright and blue, reddish and slow, and red and faint models, these strategies are only defined for 602, 496, and 131 simulations, respectively, in the HN scenario, while LD is defined for 532, 511, and 251.

interesting, in particular for low-TT, due to the cost of deep exposures over a wide area.

There is also a flow of decision-making external to what we have described. Whether or not there is a short gamma-ray burst (GRB), we choose distinct approaches. In the case where there is, we are likely looking for a GW170817-like event (bright and blue). If not, then a conservative model is indicated (reddish and slow), as we might not be looking for an event with a high inclination angle. If the alert indicates it is an NS–BH merger, then the KN is likely more consistent with the red and faint model. After this, if there is a possible HN strategy, it is low-budget and high-performing by definition. It also gives enough time for us to use spectroscopy to quickly confirm. In an observing run where we expect one BNS merger every month, we could plan on allocating telescope time through the run considering the amount of telescope time remaining versus where the particular event lies in an S/N distribution (Chen & Holz 2014) or among a population of simulated events. This would lead to spending more toward the end of the run if considerable time was left or spending more if it was a particularly good event compared to the simulations.

### 6.1. Plan for Usage

To demonstrate how using this code might work during a real GW follow-up campaign, we randomly select one of the simulations used in this analysis. For this exercise, we are assuming that the LVK trigger announcement was received toward the beginning of the night, meaning we will aim to be on-sky within a few hours of merger. The general logic flow for determining the observing plan goes as follows.

1. Assess our external factors.

(a) Find out how many hours we are allotted for the night and about telescope availability.

(b) Find out the sky conditions.

(c) Check if there was a GRB reported in the area at the same time (as was the case with GW170817). This is important for choosing a KN model to use.

2. Using the external factors as input, run the code.

3. Assess output csv; see Table 7. Here we would decide if this is an event that we would like to follow up, considering the discovery probability.

4. Choose filters, exposure times, and area coverage from strategy and compile an observing plan.

While most external factors are easy to identify, choosing which KN model to follow may not be obvious and will depend on the particular science goal of the project. For instance, if one’s science goal is to have the most complete set of KNe, then using the faintest KN model is useful; thus, red and faint is the favored model. For simplicity, in this example, we will assume that a GRB was reported within the LVK sky area, thus favoring bright and blue. Once the code has generated an output csv file, we can compare strategies.

The configuration with the highest detection probability for each strategy is displayed in Table 7. Here we can see that for this event, the top strategy yields the highest probability of detection, as expected. In this example, however, the top strategy has a factor of almost 5 more in the amount of telescope time needed for any other strategy for only an ~7% greater probability of detection.

## 7. Discussion

In this paper, our science goal is to maximize the number of KNe with GW-measured $d_L$ and securely identified redshifts. We put a high weight on the completeness of detection of the KNe, given the considerable uncertainty in the KN population. Our strategies go deep. The strategy chosen by our optimization considers a constraint in telescope time, i.e., in low-TT mode for a reddish and slow KN reach $m_{lim}(10\sigma)$ of $r \ge 24.4$ and $z \ge 22.9$ for the $\ge 1200$ s exposures (see Table 2) for the inner regions of 47% (as defined in Table 3) from our set of 611 simulations, regardless of distance (see Figure 7).

We can gain insight into our results by comparing how the literature handles a set of three questions:

1. the modeling of GW merger event distance and sky area distribution;

2. the range of KN model physical parameters, including inclination angle; and

3. the telescope, search cadence, detectability versus distance and $m_{lim}$.

*The GW event properties.* We simulate a merging NS–NS population expected in O4 given the expected LIGO sensitivities, drawing from the NS population, simulating the GW waveform, projecting onto a GW observatory network, and converting the GW observations to sky maps. Our approach to drawing, simulating, projecting, and forming sky maps has been performed before (Chen et al. 2021; Petrov et al. 2022), though it is sometimes done without sky maps (Zhu et al. 2021). There are simplified approaches. One can just assume a KN population to observe (Sagués Carracedo et al. 2021; Chase et al. 2022; Setzer et al. 2023). Or, one can use toy GW models, such as a low-significance event sky map at 200 Mpc (Coughlin et al. 2020). All of these approaches assume that the KN properties are independent of the GW event properties, other than distance. Colombo et al. (2022)
<table><caption>**Table 7**

Excerpt of Output csv for Example Simulated Event</caption><thead><tr><th>Strategy</th><th>Detection Prob. (%)</th><th>Filter</th><th>Exposure Time Outer, Inner (s)</th><th>Integrated Prob. Area Outer, Inner</th><th>Total Telescope Time (Two Detections; hr)</th><th>Confirmation Day (Days after Alert)</th></tr></thead><tbody><tr><td>Top</td><td>75</td><td>r</td><td>3600, 5400</td><td>0.9, 0.5</td><td>14.7</td><td>1.5</td></tr><tr><td>Low-TT</td><td>68</td><td>z</td><td>300, 1200</td><td>0.9, 0.7</td><td>3.2</td><td>1.0</td></tr><tr><td>HN</td><td>69</td><td>z</td><td>600, 1200</td><td>0.9, 0.7</td><td>3.5</td><td>1.0</td></tr><tr><td>ED</td><td>65</td><td>i</td><td>600, 1200</td><td>0.9, 0.8</td><td>2.9</td><td>0.5</td></tr><tr><td>LD</td><td>65</td><td>i</td><td>600, 1200</td><td>0.9, 0.7</td><td>3.5</td><td>2.0</td></tr><tr><td>Reference</td><td>48</td><td>i</td><td>90</td><td>0.9</td><td>0.5</td><td>1.5</td></tr></tbody></table>

**Note.** The table displays the highest-probability configuration for each strategy using scenario 2 (using deeper exposures in a higher-probability sky area) for the bright and blue model. The “Exposure Time Outer, Inner” column represents the time used for the shallow and deep areas used out to the area listed in the “Integrated Prob. Area Outer, Inner” column. This summary information will be used in the decision tree of which strategy to use for an actual LVC event trigger. Dual exposure times and probability areas describe the exposure time and area covered for the “deep” inner and “shallow” outer region exposures of Table 3.

performs the more sophisticated analysis of connecting the KN properties to the GW event properties by going through the chirp mass and mass ratio, which informs $M_{ej}$ and $v_{ej}$. Of the analyses that go through the draw, simulate, and project methodology, all assume the $S/N > 8$ or network $S/N > 12$ in selection, except for Petrov et al. (2022), who derive the effective $S/N$ threshold from the published GWTC events. They find that there will likely be more events with 90% credible sky areas >300 deg$^2$ than we simulate but do not include in our follow-up identification analysis. For the analyses that extend to jet production (e.g., Zhu et al. 2021; Colombo et al. 2022), about 10% of the events are sufficiently pole-on for jets to be observed. It is worth noticing, however, that the simulations used are based on assumptions about the sensitivity of the LVK detectors. The predicted sensitivities are roughly equivalent to those anticipated around mid-2022 for the upcoming observing cycle O4. The actual operational conditions diverge from this expected design sensitivity. The Virgo detector, for instance, has not joined the O4 from its beginning, and it remains uncertain whether it will achieve the proposed sensitivity for this cycle.

Moreover, the KAGRA detector has participated in the run for a limited period but at a reduced BNS inspiral range from that estimated here (80 Mpc). The adjustments in observational capabilities were publicly disclosed only when this work was nearing completion. As a result, the predictions offered in this study may be optimistic when evaluated in the context of the present O4 run. Nonetheless, these forecasts retain their significance in terms of methodology for guiding observations, even if they do not reflect the current O4 schedule.

**KN physical properties and inclination angles.** We have only one well-studied KN in the literature, so there is considerable uncertainty in the KN population. The literature has three ways of modeling KNE: using model atmospheres as building blocks (e.g., Kasen et al. 2017), full physical models of KN atmospheres (e.g., Bulla 2019; Wollaeger et al. 2021), and using scaling relations and fitting functions from relativistic numerical simulations (e.g., Dietrich & Ujevic 2017; Coughlin et al. 2019). The first two approaches are parameterized by at least three variables (for us, $M_{ej}$, $v_{ej}$, $\log X_{\text{lan}}$). The Bulla models, which explicitly include inclination angle, are widely used in the relevant literature (e.g., Coughlin et al. 2020; Sagués Carracedo et al. 2021; Zhu et al. 2021; Petrov et al. 2022). Often, analyses in the literature will set up a grid of KN parameters and inclination angles and proceed to simulate detection of each entry in the grid, either for a fixed distance or for an NS merger population. Typically, the fraction of models

detected is termed the detection efficiency. This makes the detection efficiency depend on the model space in unfortunate ways. Consider the case of inclination angle-dependent KN properties. One can set up a grid of inclination angles. Better would be to use the probability distribution function of inclination angles for a random isotropic inclination sample, $PDF(i) \propto \sin(i)$. This weights edge-on, i.e., red and faint, KNe more than a grid is likely to. The inclination angles sample-selected by a GW detector network search are not isotropic. The amplitude of the strain detected is generically described by $PDF(i) = 0.076076(1 + 6 \cos^2(i) + \cos^4(i))^{3/2} \sin(i)$; that is, GW observatories prefer inclination angles near 30° (Schutz 2011). For the Bulla models, this could be accounted for by using Bayesian priors on the models. It would be interesting to have a version of the Setzer et al. (2023) KN population absolute magnitude distribution weighted by the expected GW network inclination angle distribution. In our work, we use Bayesian model averaging on our calculated detection probabilities. We do not model polar versus equatorial directly, but one could map the appropriate blue and red Bayesian model average into those. Our red and faint average absolute magnitude is about 0.4 mag fainter than the brighter peak of Setzer et al. (2023).

**Detection efficiency.** Our study is for DECam on the Blanco 4 m with its 3 deg$^2$ field of view, covering the 90% credible sky area using real cadences in two passes to ensure a confirmation. A complete search-and-discovery simulation over full LVK simulated sky maps has been done by Coughlin et al. (2020) and Petrov et al. (2022). As the 1 m class telescopes have very large fields of view, up to the ZTF 47 deg$^2$, much more common in the literature is to assume a KN model, usually a GW170817 analog, and ask what exposure time or limiting magnitude is necessary to detect it. If instead of Bayesian model averaging of the detection probability, we had used Bayesian model average absolute magnitudes, our study would be very different. The average absolute magnitudes for our bright and blue, reddish and slow, and red and faint models for the r band are -16.3, -15.2, and -14.2, respectively. In 100 s, DECam reaches $r = 23.0$, sufficient to detect $M_r = -15.2$ to 500 Mpc. In 1200 s, DECam reaches $r = 24.4$, sufficient to detect $M_r = -14.2$ to 500 Mpc and GW170817's $M_r = -15.5$ to 1 Gpc. Why, then, do events at $d_L = 200$ Mpc have in our study detection probabilities $P_d$ of 90%, 73%, and 60% for the bright and blue, reddish and slow, and red and faint models, respectively, using the low-TT strategy and routinely require 60–90 s exposures for bright and blue and 300–1200 s exposures for the last two models? The effect of using
<table><caption>**Table 7**

Excerpt of Output csv for Example Simulated Event</caption><thead><tr><th>Strategy</th><th>Detection Prob. (%)</th><th>Filter</th><th>Exposure Time Outer, Inner (s)</th><th>Integrated Prob. Area Outer, Inner</th><th>Total Telescope Time (Two Detections; hr)</th><th>Confirmation Day (Days after Alert)</th></tr></thead><tbody><tr><td>Top</td><td>75</td><td>r</td><td>3600, 5400</td><td>0.9, 0.5</td><td>14.7</td><td>1.5</td></tr><tr><td>Low-TT</td><td>68</td><td>z</td><td>300, 1200</td><td>0.9, 0.7</td><td>3.2</td><td>1.0</td></tr><tr><td>HN</td><td>69</td><td>z</td><td>600, 1200</td><td>0.9, 0.7</td><td>3.5</td><td>1.0</td></tr><tr><td>ED</td><td>65</td><td>i</td><td>600, 1200</td><td>0.9, 0.8</td><td>2.9</td><td>0.5</td></tr><tr><td>LD</td><td>65</td><td>i</td><td>600, 1200</td><td>0.9, 0.7</td><td>3.5</td><td>2.0</td></tr><tr><td>Reference</td><td>48</td><td>i</td><td>90</td><td>0.9</td><td>0.5</td><td>1.5</td></tr></tbody></table>

**Note.** The table displays the highest-probability configuration for each strategy using scenario 2 (using deeper exposures in a higher-probability sky area) for the bright and blue model. The “Exposure Time Outer, Inner” column represents the time used for the shallow and deep areas used out to the area listed in the “Integrated Prob. Area Outer, Inner” column. This summary information will be used in the decision tree of which strategy to use for an actual LVC event trigger. Dual exposure times and probability areas describe the exposure time and area covered for the “deep” inner and “shallow” outer region exposures of Table 3.

performs the more sophisticated analysis of connecting the KN properties to the GW event properties by going through the chirp mass and mass ratio, which informs $M_{ej}$ and $v_{ej}$. Of the analyses that go through the draw, simulate, and project methodology, all assume the $S/N > 8$ or network $S/N > 12$ in selection, except for Petrov et al. (2022), who derive the effective $S/N$ threshold from the published GWTC events. They find that there will likely be more events with 90% credible sky areas >300 deg$^2$ than we simulate but do not include in our follow-up identification analysis. For the analyses that extend to jet production (e.g., Zhu et al. 2021; Colombo et al. 2022), about 10% of the events are sufficiently pole-on for jets to be observed. It is worth noticing, however, that the simulations used are based on assumptions about the sensitivity of the LVK detectors. The predicted sensitivities are roughly equivalent to those anticipated around mid-2022 for the upcoming observing cycle O4. The actual operational conditions diverge from this expected design sensitivity. The Virgo detector, for instance, has not joined the O4 from its beginning, and it remains uncertain whether it will achieve the proposed sensitivity for this cycle.

Moreover, the KAGRA detector has participated in the run for a limited period but at a reduced BNS inspiral range from that estimated here (80 Mpc). The adjustments in observational capabilities were publicly disclosed only when this work was nearing completion. As a result, the predictions offered in this study may be optimistic when evaluated in the context of the present O4 run. Nonetheless, these forecasts retain their significance in terms of methodology for guiding observations, even if they do not reflect the current O4 schedule.

**KN physical properties and inclination angles.** We have only one well-studied KN in the literature, so there is considerable uncertainty in the KN population. The literature has three ways of modeling KNE: using model atmospheres as building blocks (e.g., Kasen et al. 2017), full physical models of KN atmospheres (e.g., Bulla 2019; Wollaeger et al. 2021), and using scaling relations and fitting functions from relativistic numerical simulations (e.g., Dietrich & Ujevic 2017; Coughlin et al. 2019). The first two approaches are parameterized by at least three variables (for us, $M_{ej}$, $v_{ej}$, $\log X_{\text{lan}}$). The Bulla models, which explicitly include inclination angle, are widely used in the relevant literature (e.g., Coughlin et al. 2020; Sagués Carracedo et al. 2021; Zhu et al. 2021; Petrov et al. 2022). Often, analyses in the literature will set up a grid of KN parameters and inclination angles and proceed to simulate detection of each entry in the grid, either for a fixed distance or for an NS merger population. Typically, the fraction of models

detected is termed the detection efficiency. This makes the detection efficiency depend on the model space in unfortunate ways. Consider the case of inclination angle-dependent KN properties. One can set up a grid of inclination angles. Better would be to use the probability distribution function of inclination angles for a random isotropic inclination sample, $PDF(i) \propto \sin(i)$. This weights edge-on, i.e., red and faint, KNe more than a grid is likely to. The inclination angles sample-selected by a GW detector network search are not isotropic. The amplitude of the strain detected is generically described by $PDF(i) = 0.076076(1 + 6 \cos^2(i) + \cos^4(i))^{3/2} \sin(i)$; that is, GW observatories prefer inclination angles near 30° (Schutz 2011). For the Bulla models, this could be accounted for by using Bayesian priors on the models. It would be interesting to have a version of the Setzer et al. (2023) KN population absolute magnitude distribution weighted by the expected GW network inclination angle distribution. In our work, we use Bayesian model averaging on our calculated detection probabilities. We do not model polar versus equatorial directly, but one could map the appropriate blue and red Bayesian model average into those. Our red and faint average absolute magnitude is about 0.4 mag fainter than the brighter peak of Setzer et al. (2023).

**Detection efficiency.** Our study is for DECam on the Blanco 4 m with its 3 deg$^2$ field of view, covering the 90% credible sky area using real cadences in two passes to ensure a confirmation. A complete search-and-discovery simulation over full LVK simulated sky maps has been done by Coughlin et al. (2020) and Petrov et al. (2022). As the 1 m class telescopes have very large fields of view, up to the ZTF 47 deg$^2$, much more common in the literature is to assume a KN model, usually a GW170817 analog, and ask what exposure time or limiting magnitude is necessary to detect it. If instead of Bayesian model averaging of the detection probability, we had used Bayesian model average absolute magnitudes, our study would be very different. The average absolute magnitudes for our bright and blue, reddish and slow, and red and faint models for the r band are -16.3, -15.2, and -14.2, respectively. In 100 s, DECam reaches $r = 23.0$, sufficient to detect $M_r = -15.2$ to 500 Mpc. In 1200 s, DECam reaches $r = 24.4$, sufficient to detect $M_r = -14.2$ to 500 Mpc and GW170817's $M_r = -15.5$ to 1 Gpc. Why, then, do events at $d_L = 200$ Mpc have in our study detection probabilities $P_d$ of 90%, 73%, and 60% for the bright and blue, reddish and slow, and red and faint models, respectively, using the low-TT strategy and routinely require 60–90 s exposures for bright and blue and 300–1200 s exposures for the last two models? The effect of using
Equation (2) is to extend the search to lower $M_{ej}$ and $\log X_{\text{lan}}$ and thus lower luminosities. Our method of accounting for uncertainty in the KN population is driving our results. If we adopt the Coughlin et al. (2018) model parameters and uncertainties, our detection efficiencies will increase and exposure times decrease. Using the GW170817 absolute magnitude for detection is likely overly optimistic, as suggested by Colombo et al. (2022) placing GW170817 at the 75th percentile bright and Setzer et al. (2023) placing it at the 95th percentile bright. The studies using 1 m class telescopes are most likely to assume GW170817 analogs, though Petrov et al. (2022) uses both that and a lower-luminosity red model, and Sagués Carracedo et al. (2021) does a careful analysis of viewing angle-dependent models. The studies assuming the Rubin Observatory (Cowperthwaite et al. 2019; Chen et al. 2021) or a variety/network of telescopes (Coughlin et al. 2020; Chase et al. 2022) tend to analyze detection probabilities for lower-luminosity events. In summary, one will have to be careful comparing our detection probabilities with others in the literature, which often use a single luminous model or evaluate and average detection efficiencies over a grid of models using a uniform prior. In our language, “chance of finding it” in Figure 4 is to be interpreted as “fraction of models detected given our priors on the space of models.”

## 7.1. Applicability to NS-BH and Mass Gap Events

For standard siren studies, NS-BH mergers are just as valuable as BNS mergers, as long as they produce EM counterparts. NS-BH merger events have higher distances for similar S/N than BNS mergers, as can be seen in Table 8.

The dominant factor for use in standard siren cosmology is the probability of a KN given an NS-BH merger. No counterpart to an NS-BH merger has been detected (e.g., Kawaguchi et al. 2020a; Morgan et al. 2020; Anand et al. 2021). Zhu et al. (2021) argue that no detectable KN counterpart was expected for NS-BH mergers in O3. Kawaguchi et al. (2020b) and Darbha et al. (2021) study the brightness of KNe from NS-BH mergers, and Drozda et al. (2020) does the same for mass gap objects. The summary is that only a fraction of NS-BH events will produce KNe, primarily those mergers with low mass ratios and high spin.

The Kawaguchi et al. (2020b) models have absolute magnitudes that peak for *r* and *i* at −14.5 and −15.0, respectively. There is a spread of about 1 mag fainter in the *i*-band absolute magnitude, going fainter as the binary mass ratio increases and the effective spin gets smaller. Our red and faint model has an *i*-band absolute magnitude of −14.5, midway through the range of Kawaguchi et al. (2020b). Petrov et al. (2022) adopt the Bulla (2019) models, broken into BNS and NS-BH models both optimistic ($M_{ej} = 0.05, 0.08$) and conservative ($M_{ej} = 0.01, 0.01$). The optimistic BNS model has an absolute magnitude in the *r* band of −16.0, and the conservative NS-BH model has −14.8; these correspond well to our bright and blue and red and faint models, respectively. Taking our red and faint model as appropriate for dynamical ejecta-dominated NS-BH mergers, our low-TT strategy has a 50% detection probability out to 330 Mpc (see Figure 4). Our strategies are sufficient to obtain the majority of NS-BH events that have EM counterparts if they are at distances ≪330 Mpc.

<table><caption>**Table 8**
LVC GWTC Events (O1–O3) Containing NSs and with $M_2 < 5 M_\odot$ and $S/N > 8$</caption><thead><tr><th>ID</th><th><i>M</i><sub>1</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>2</sub><br>(<i>M</i><sub>☉</sub>)</th><th>Distance<br>(Mpc)</th><th>S/N</th><th>Class</th></tr></thead><tbody><tr><td>GW170817</td><td>1.5</td><td>1.3</td><td>40</td><td>33</td><td>BNS</td></tr><tr><td>GW190425</td><td>2.0</td><td>1.4</td><td>160</td><td>13</td><td>BNS</td></tr><tr><td>GW190814</td><td>23.2</td><td>2.6</td><td>240</td><td>22</td><td>BH-NS</td></tr><tr><td>GW200105_162426*</td><td>9.0</td><td>1.9</td><td>270</td><td>14</td><td>BH-NS</td></tr><tr><td>GW200115_042309</td><td>5.9</td><td>1.4</td><td>290</td><td>11</td><td>BH-NS</td></tr><tr><td>GW190426_152155</td><td>5.7</td><td>1.5</td><td>370</td><td>10</td><td>BH-NS</td></tr><tr><td>GW191219_163120</td><td>31.1</td><td>1.2</td><td>550</td><td>9</td><td>BH-NS</td></tr><tr><td>GW190917_114630</td><td>9.3</td><td>2.1</td><td>720</td><td>10</td><td>BH-NS</td></tr><tr><td>GW200210_092254</td><td>24.1</td><td>2.8</td><td>940</td><td>8</td><td>BH-NS</td></tr></tbody></table>

Note. S/N is the matched filter S/N. GW200105_162426* has $p_{astro} = 0.3$ and is thus considered a marginal candidate. Only two events have a 90% confidence sky area of <30 deg$^2$, two at <300 deg$^2$, and two at <400 deg$^2$. The median sky area is 1700 deg$^2$. Data are from https://www.gw-openscience.org/eventapi/html/allevents/.

## 7.2. Blanco/DECam and Rubin LSST

It is of interest to compare the strategies defined here with the program outlined in Chen et al. (2021). They assume inclination independence and a GW170817-like KN and argue for two filter observations. The program conservatively assumes 30% of A+ events by dedicating 7 hr of Rubin Observatory time in 30 s exposures, capturing $12 \text{ yr}^{-1}$. Our expectation is that the Rubin ToO program will use 3% of the available LSST time, so on order of 100 hr, which can pursue all BNS events in LVK O4 assuming $1 \text{ month}^{-1}$ and 8 hr event$^{-1}$, so pursuing light curves. Alternatively and more likely, Rubin will choose to observe the 50% best events by sky area in both the BNS and the NS-BH categories. In this scenario, a good use of DECam/Blanco would be to follow up the others that have sky area <300 deg$^2$. Table 5 suggests that this would be viable.

In fact, if the results of Petrov et al. (2022) hold, then there will be many merger events containing NSs that have sky areas greater than 300 deg$^2$; our simulations would have to be extended by another ~360 events, all with sky area >300 deg$^2$, to match their statistics. For the bright siren cosmology, every NS event is important. We demonstrate here that the Blanco/DECam, especially in combination with the Zwicky observatory and its counterparts PS1, OAJ, and LS4, is capable of following up the sources with sky area <300 deg$^2$. The optimal use of the Vera C. Rubin Observatory, with its immense endowment, is to follow up the LVK sources with >300 deg$^2$. The combination of sky coverage and depth is unmatched. Petrov et al. (2022) predicts that the median sky coverage for BNS events in O4 is $1820_{-170}^{+190}$ deg$^2$, the median luminosity distance is $352 \pm 10$ deg$^2$, and the NS-BH median distance is further away. For the Rubin field of view of 9.6 deg$^2$, the number of exposures to cover the sky area once is ≈200, which at 100 s exposures can be done in less than 6 hr, assuming a reddish and slow model and 1.2 mag deeper $m_0$ for Equation (1). Likely one could build a two-visit strategy that would take 10 hr per event, allowing Rubin to follow up 10 additional events per year without light curves. The Rubin time-domain ecosystem of data, brokers, and routine spectroscopic follow-up is likely to minimize positives, though perhaps not until after O4.