1Introduction

The detection of the binary neutron star (BNS) merger GW170817 [4,5] on August 17, 2017, by the Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) [6] and Advanced Virgo [7] marked a watershed moment in multimessenger astronomy. This event, accompanied by the short gamma-ray burst (GRB) GRB170817A detected by NASA’s Fermi Gamma-ray Space Telescope [8] and the European Space Agency’s INTEGRAL [9], and followed by the kilonova AT2017gfo [10,11] observed in the galaxy NGC4993, provided an unprecedented opportunity to study the Universe through multiple messengers: gravitational waves (GW) and electromagnetic (EM) waves. Occurring approximately 40Mpc away, GW170817 signature chirp lasted nearly 100 seconds, far longer than typical black hole merger signals, indicated the collision of two neutron stars with masses of about \(1.1\) and \(1.6\) solar masses, confirmed by the near-instantaneous gamma-ray flash just \(1.74\) seconds later. This joint detection, the first of its kind, alerted astronomers worldwide, leading to a follow-up campaign that captured the fading glow of AT2017gfo, interpreted as a kilonova—an EM transient powered by the radioactive decay of r-process elements synthesized during the merger.

The significance of GW170817 extends across multiple fields, enhancing our understanding of cosmology and Hubble tension [12–21], nuclear physics [22–25,4,26,27,19,28], modified gravity theories [29–31], and the chemical evolution of the Universe [32–37,36]. While GW170817 remains the only confirmed multimessenger BNS merger, observational evidence suggests that GW190425 [38], GRB211211 [39–41], and GRB230307A [42] may also originate from similar mergers, hinting at an increasing likelihood of future detections as observational facilities improve.

The interpretation of GW170817 and future multimessenger events relies on comparing observational data with robust theoretical models. For GWs, systematic uncertainties are well-controlled for current detectors [43–47], though next-generation observatories like the Einstein Telescope [48–50] and Cosmic Explorer [51] will demand even higher precision. However, with the EM counterpart, the challenges are more complex. Kilonovae, spanning optical, ultraviolet, and infrared wavelengths, are driven by the radioactive decay of unstable heavy elements formed via r-process nucleosynthesis [52,53]. Despite progress from simplified semi-analytical models [54,11,55,37,53,56–60] to sophisticated 3D radiative transfer simulations [61–66] uncertainties persist in ejecta properties (mass, velocity, composition) [67–71], nuclear heating rates, thermalization efficiencies, and opacities of r-process elements [72–79]. These uncertainties, if unaddressed, risk biasing parameter estimates even if the data quality improves.

This thesis investigates how these uncertainties evolve over time and per filter, and how they impact the inferred ejecta properties of kilonovae. By leveraging the Nuclear-Physics and Multi-Messenger Astrophysics (NMMA) framework and high-fidelity radiative transfer simulations from POSSIS  [64], we aim to develop a more robust methodology for incorporating time- and filter-dependent uncertainties into kilonova analyses [1]. Our study focuses on AT2017gfo, the best-observed kilonova to date, as a case study to assess the effects of evolving systematic uncertainties on light curve modeling and parameter estimation. The insights gained from this work will contribute to refining kilonova models for future multimessenger detections.

2Anatomy of a Binary Neutron Star Merger

2.1Neutron Stars

Neutron stars emerge from the remnants of massive stars after a supernova explosion expels the outer layers, leaving a core that collapses under gravity. When the core’s mass exceeds the Chandrasekhar limit (approximately \(1.4\ M_\odot\)), electron capture and neutronization convert protons and electrons into neutrons. If the core mass exceeds a critical value (around \(2-3\ M_\odot\)), the collapse continues until neutron degeneracy pressure, supported by strong nuclear interactions at extreme densities, halts further collapse and forms a neutron star. Fig. 2.1 illustrates the evolution of a BNS system, highlighting the potential outcomes based on the masses and rotational energy of the remnant.

2.1.1Tolman-Oppenheimer-Volkoff Equation

The Tolman-Oppenheimer-Volkoff equations [81,82] describe the structure of a static, spherically symmetric neutron star in general relativity, extending Newtonian hydrostatic equilibrium to account for relativistic effects. They emerge from solving the Einstein Field Equations [83] for a fluid ball under these conditions, providing the pressure, density, and mass profiles as functions of radius.

In Newtonian mechanics, the balance between gravitational and pressure forces in a spherical body is governed by the equation of hydrostatic equilibrium. However, as the density increases to extreme levels within compact objects like neutron stars, relativistic corrections become essential. The TOV equations are derived from the Einstein Field Equations for a perfect fluid, under the assumptions of staticity (no motion) and spherical symmetry. The general relativistic metric for a static, spherically symmetric star is given by

where \(\Phi(r)\) and \(\lambda(r)\) are functions of the radial coordinate \(r\), representing the gravitational potential and the space-time curvature, respectively. The perfect fluid model is described by the stress-energy tensor \(T_{\mu\nu}\), which for a static, spherically symmetric star takes the form,

where \(\rho\) is the energy density, \(p\) is the pressure, and \(u_\mu\) is the four-velocity of the flluid.

The TOV equations are derived by solving the Einstein Field Equations (EFE) with the metric in Eq. 2.1 and the stress-energy tensor in Eq. 2.2, yielding three first-order differential equations

where \(M(r)\) is the mass enclosed within radius \(r\), and \(P(r)\) is the pressure at that radius.

Eq. 2.3a expresses the mass conservation within a spherical shell of radius \(r\), Eq. 2.3b describes the pressure gradient in the star, balancing the gravitational force with the pressure gradient, and Eq. 2.3c relates the gradient of the gravitational potential to the pressure gradient, essentially describing how the spacetime metric varies throughout the star. The TOV equations are typically solved numerically to obtain the mass, radius, and density profiles of neutron stars, providing insights into their structure and properties.

2.1.2Equation of State

The Tolman-Oppenheimer-Volkoff (TOV) equations are fundamentally dependent on the equation of state (EOS) of neutron star material, which defines the relationship between pressure \(P\) and density \(\rho\). The EOS encapsulates the microphysics of ultra-dense matter, governing how neutron stars resist gravitational collapse under extreme conditions. However, the EOS at supranuclear densities (\(\gtrsim10^{14} \ \rm g/cm^3\)) remains poorly constrained due to limitations in terrestrial experiments and theoretical uncertainties in quantum chromodynamics and nuclear interactions.

For instance, a soft EOS, characterized by lower pressures at high densities, results in greater compressibility, yielding smaller and less massive neutron stars with radii typically around \(10 - 11 \ \rm km\) and maximum masses closer to \(1.4 - 2 \ M_\odot\) [84]. Conversely, a stiff EOS, with higher pressures resisting compression, supports more massive neutron stars—potentially exceeding \(2.5 \ M_\odot\)with larger radii, up to \(13\ \rm km\) [85,27]. These variations arise because the TOV pressure gradient balances gravity more effectively with a stiff EOS, allowing greater mass accumulation before collapse to a black hole occurs, as supported by recent observational constraints from NICER and GW170817 [86,27].

2.2Kilonova

Kilonovae, transient electromagnetic events are triggered by the merger of neutron stars or a neutron star and a black hole. The name kilonova [87] stems from their luminosity, which sits roughly 1000 times above that of a classical nova yet below a typical supernova [52]. Kilonovae are characterized by their multi-wavelength emission, spanning optical, ultraviolet, and infrared bands, and exhibit a distinctive light curve evolution over time. Few of the first possible kilonova candidates were observed as the aftermath of the short gamma-ray bursts include GRB060614 [88,89], GRB070809 [90], GRB080503 [91] GRB130603B [92], followed by the detection of AT2017gfo associated with the binary neutron star merger GW170817 [10,11].¹Even though a number of GRBs have been associated with neutron star mergers (these associations arise due to the lack of a supernova like component in the spectral features of these GRBs), AT2017gfo/GRB170817A/GW170817 is the only confirmed kilonova candidate so far.

The kilonova’s optical and infrared emission, peaking within a few days to weeks post-merger, is powered by the radioactive decay of r-process elements, primarily lanthanides and actinides, synthesized during the merger. The decay chain of these heavy elements releases a large amount of energy, heating the ejecta and driving the kilonova’s luminosity. The kilonova’s spectral evolution is characterized by a blue to red color transition, reflecting the changing composition and temperature of the expanding ejecta.

The early blue emission, peaking in the ultraviolet and optical bands, arises from the radioactive decay of light r-process elements, while the later red emission, peaking in the infrared, is powered by the decay of heavier elements. The kilonova’s light curve exhibits a rapid rise to peak luminosity, followed by a gradual decline over weeks to months as the ejecta cools and expands. The kilonova’s peak luminosity, color, and duration depend on the ejecta mass, velocity, composition, and thermalization efficiency, making them valuable probes of the merger dynamics and the properties of neutron star matter. Theoretical models of kilonovae have evolved from simple analytical prescriptions [52,87] to sophisticated radiative transfer simulations [61–66], capturing the complex interaction between nuclear heating, radiative transport, and opacities of r-process elements.

2.2.1Ejecta Properties

The material ejected during a neutron star merger forms the basis of a kilonova, driving its observable signatures through radioactive decay and thermal emission. Ejecta originate from multiple mechanisms, each contributing distinct mass, velocity, and compositional profiles that shape the resulting light curves and spectra.

2.2.1.1Dynamical Ejecta

Dynamical ejecta from NS-NS mergers consist of matter expelled on the characteristic dynamical timescale (\(\sim 10 \rm ms\)) of the system, with typical masses on the order of \(10^{-4}-10^{-2}\,M_\odot\) and velocities of \(0.1-0.3c\) [93–96]. The total mass of this ejecta depends strongly on the binary’s total mass, its mass ratio, and the neutron star EOS. Eccentric NS-NS mergers can yield higher ejecta masses [97,98], but such mergers likely occur less frequently in nature due to the high stellar densities required for significant dynamical capture events [99]. Simulations also indicate that very high neutron star spin can enhance the amount of dynamically ejected matter [100–102]. When a neutron star merges with a black hole (NSBH), the tidal disruption of the star by the black hole can unbind up to \(\sim 0.1\,M_\odot\) of material, typically at similar velocities to those in NS-NS collisions [103–105].

In NS-NS systems, two hydrodynamical processes dominate the dynamical ejection. First, the contact interface between merging stars experiences strong compression and shock heating, which drives material outward as the remnant undergoes quasi-radial pulsations [106,107,93]. This “shock-driven” component of the ejecta emerges in a broad range of angular directions. Second, tidal interactions form spiral arms in the equatorial plane, and angular momentum redistribution within these arms propels additional matter outward. The relative importance of these mechanisms depends on the binary mass ratio, strongly asymmetric mergers (\(q \ll 1\)) tend to eject more mass [107,108]. Ejecta masses also decrease if the remnant collapses promptly into a black hole, since prompt collapse swallows the shocked interface region before much matter can be expelled.

In NSBH mergers, the dominant ejection mechanism is a tidal disruption of the neutron star by the black hole, which also focuses the outflow in the equatorial plane [109,104]. The azimuthal extent of this tidal ejecta can be limited, implying that kilonova emission may depend more sensitively on viewing angle in NSBH systems than in NS-NS mergers.

Besides mass and velocity, the electron fraction \(Y_e\) (ratio of protons to total nucleons) of the dynamical ejecta is crucial for determining its role in \(r\)-process nucleosynthesis and kilonova emission. Early simulations that neglected weak interactions found that merger outflows are highly neutron-rich (\(Y_e \lesssim 0.1\)), suggesting robust production of heavy elements with \(A \gtrsim 130\) through fission recycling [110,111,107,112]. More recent calculations that include electron captures and neutrino irradiation reveal a broader range of electron fractions \(Y_e \sim 0.1-0.4\), thereby allowing for the synthesis of lighter \(r\)-process elements as well [113,94,114]. In NS-NS mergers, the higher \(Y_e\) ejecta tends to be more isotropic, whereas the lower \(Y_e\) (and hence lanthanide-rich) material often remains concentrated in the equatorial plane, reflecting its origin in the tidal arms. This results in lanthanide-rich tidal ejecta (\(Y_e < 0.25\)) and lanthanide-poor polar ejecta (\(Y_e > 0.25\)), influencing the kilonova’s spectral evolution and color.

2.2.1.2Wind Ejecta

Wind ejecta are launched after the merger from the accretion disk or the central compact object, in contrast to the dynamical ejecta expelled promptly during the coalescence phase. The formation of a massive accretion disk (on the order of \(0.01$--$0.3 \,M_\odot\)) is well established in numerical simulations and depends sensitively on the total mass of the binary, the mass ratio, the spins of the merging bodies, and the neutron star equation of state [115,116,53]. Prompt collapse of the merger remnant into a black hole typically leads to lower disk masses because there is little time for the remnant to redistribute its mass and angular momentum prior to collapse [117].

Once the disk forms, a combination of neutrino heating, magnetic stresses, and nuclear recombination drives sustained mass loss on timescales ranging from milliseconds to seconds. During the earliest, high-accretion phase, strong neutrino irradiation from the disk (and, if present, a hypermassive neutron star) can unbind material from the disk surface [118]. If the central neutron star survives for longer than \(\sim 50\) ms, additional neutrino luminosity from the star itself ejects non-negligible amounts of mass [119–121], potentially enhanced further by spiral density waves excited in the disk [122]. Over time, as the disk expands viscously and its temperature drops, it transitions from an initially neutrino-cooled configuration to a radiatively inefficient, geometrically thicker state [123,116,124,125]. In this later phase, weak interactions “freeze out,” allowing the outflows to remain neutron-rich [126,123].

The wind outflows typically attain velocities of \(0.1$--$0.3c\) and can have electron fractions \(Y_e\) spanning a broad range, from below 0.25 to above 0.3, depending on the level and duration of neutrino irradiation [127,117,128]. Outflows with lower \(Y_e\) remain rich in lanthanides and produce red or infrared components of the kilonova, while higher-\(Y_e\) material can generate bluer emission [129,130].

Numerical simulations find that the fraction of the disk ultimately ejected through these winds can be substantial, ranging from a few percent for low black hole spin to as high as 30–40% or more for rapidly spinning black holes [127,117,131,132]. Even higher fractions, up to \(\sim 90\%\) of the disk, may be unbound if the central remnant remains a long-lived neutron star instead of collapsing promptly to a black hole [133,134]. The wind outflows typically attain velocities of \(0.1$--$0.3c\) and can have electron fractions \(Y_e\) spanning a broad range, from below 0.25 to above 0.3, depending on the level and duration of neutrino irradiation [127,117,128].

Magnetic fields also play a crucial role in shaping the wind dynamics. Turbulence generated by the magneto-rotational instability transports angular momentum outward and boosts mass loss. If the neutron star remnant is backed by a strong, large-scale magnetic field (\(\sim10^{14}$--$10^{15}\,\mathrm{G}\)), additional outflows can be driven centrifugally along open field lines [135,136]. Although the exact origin and geometry of these large-scale fields remains under investigation, they may account for particularly fast ejecta components observed in kilonovae such as GW170817 [137].

Overall, wind ejecta from the disk and central remnant can rival or exceed the mass of the prompt dynamical ejecta and thus represent a significant channel for both mass loss and the synthesis of heavy elements in neutron star mergers [138,128,131]. Their observed features, especially the broad range of velocities and electron fractions indicate that wind-driven outflows are key contributors to the multi-wavelength kilonova emission, often shaping the luminosity and color evolution on timescales of hours to days after the merger.

2.2.2Chemical Composition

The chemical composition of kilonova ejecta, determined by the electron fraction \(Y_e\), governs nucleosynthesis pathways and radiative properties of the emitted light. \(Y_e\) governs the neutron-to-proton ratio during the r-process, which occurs in the neutron-rich conditions immediately following ejection. In neutron-rich conditions (\(Y_e < 0.25\)) typical of tidal dynamical ejecta, enables the rapid neutron capture process onto the seed nuclei, producing isotopes far beyond the stability line that decay into heavy r-process elements, including lanthanides and actinides [139]. These elements, with their dense forest of atomic transitions, yield opacities of \(\kappa \sim 10-100 \ \rm{cm^2 g^{-1}}\), reddening emission and delaying photon escape [53].

On the other hand, the lanthanide-free ejecta (\(Y_e > 0.25\)) prevalent in disk wind ejecta and polar dynamical ejecta exposed to neutrino irradiation, limits neutron availability, yielding lighter r-process elements such as Strontium [36] and Yttrium [140]. These elements exhibit simpler atomic structures, resulting in opacities of \(\kappa \sim 1-10 \ \rm{cm^2 g^{-1}}\), permitting greater photon escape in the optical and ultraviolet bands. The transition between these regimes is not sharp; realistic ejecta likely feature \(Y_e\) gradient, blending lanthanide-rich and lanthanide-free regions within a single outflow [11]. This gradient complicates radiative transfer calculations but explains the multi-peaked light curves observed in events like AT2017gfo.

2.2.3Heating

The primary energy source of kilonovae lies in the radioactive decay of r-process nuclei within the ejecta, a process that heats the material and sustains its thermal emission. This decay process releases energy through beta decay, alpha decay, and fission, producing gamma rays, electrons, and neutrinos that deposit heat into the surrounding material [87]. Following Refs. [53,111], this heating rate can be expressed as

where \(X_{r,v}\) is the r-process mass fraction in mass layer \(M_v\) and \(e_r\) is the specific heating rate. For \(Y_e \lesssim 0.25\) i.e. neutron-rich ejecta \(\dot{e}_r\) can be approximated as

where \(t_0=1.3\)s and \(\sigma=0.11\)s are constants, and \(\epsilon_{\mathrm{th}, \mathrm{m}}\) is the thermalization efficiency. While Eq. 2.5 describes a constant heating rate during the first second of the kilonova, the later times can be well approximated by a power-law decay as reflecting the cumulative contribution of multiple r-process isotopes decay,

This power-law behavior characterizes the long-term energy deposition that sustains kilonova emission over days to weeks.

2.2.4Emission

The emission from kilonovae emerges as thermal radiation from the heated ejecta, evolving as photons diffuse through the expanding material. Initially, the ejecta are optically thick, meaning radiation remains trapped until the diffusion timescale is reached, approximated by,

where \(\kappa\) is the opacity, \(M\) is the mass of radius \(R\), asumming an uniform average density. Until the expansion timescale \(t_{\rm diff}= t\) (Arnett’s Law [141]) is reached, the diffusion time decreases with \(t_{\rm diff} \propto t^{-1}\) as the ejecta expands and eventually the radiation escapes the ejecta. The peak luminosity of the kilonova occurs when the diffusion time matches the expansion time, marking the transition from the optically thick to the optically thin regime, given by,

where \(\beta \approx 3\) depends on the density profile of the ejecta [107]. By varying the opacity \(\kappa\) over a broad range from lanthanide-free to lanthanide-rich matter (\(0.5-30 \ \rm{cm^2 g^{-1}}\) [142]), Eq. 2.8 yeilds characteristic kilonova durations from 1 day to 1 week.

2.3Toy Model

From Eq. 2.8, it is evident that kilonova observables are sensitive to the ejecta mass, velocity, opacity and the sources contributing to the ejecta heating \(\dot{Q}(t)\). As discussed in Sec. 2.2.1 and Sec. 2.2.2, the different ejecta components and their characteristics can be combined to model the overall kilonova emission. Following from Sec. 2.2.3, we construct a one-dimensional analytical toy model from Ref. [53], under the assumption of homologous expansion [143]. In this framework, the velocity of the ejecta is proportional to its distance from the merger center and can be expressed as,

where \(r_i\) is the radius of the material at time \(t - t_{\rm merger}\). In simpler terms, this implies that material located farther from the center moves at a higher velocity. We first describe the mass distribution of the ejecta, to this end we define it as a power-law function of the velocity,

where \(M\) is the total mass of the ejecta, \(v_0\) is the minimum velocity of the ejecta, and \(\beta\) is the power-law index and takes the value \(\approx 3\). Given the simplicity of this model and complexity of numerical simulations, a single power-law index cannot truly capture the diversity of ejecta properties, however this framework is flexible enough in order to use multiple power-law indices to describe the ejecta mass distribution. Similar to Eq. 2.7, the radiation escapes from a given mass layer \(M_v\) on a diffusion timescale \(t_{d,v}\), which can be expressed as,

where \(\kappa_v\) is the opacity of the mass layer \(M_v\) and \(R_v = vt\) reflects the homologously expanding ejecta. Subsitiuting the power-law distribution Eq. 2.10 into the diffusion timescale, we obtain,

Setting \(t_{d,v} = t\) defines the mass depth from which radiation peaks at time \(t\)

where \(t_{\text{peak}}\) is the time when diffusion occurs through the entire ejecta mass, typically evaluated at \(v = v_0\). Outer layers (lower \(M_v\)) emit first, while the innermost shell (\(M_v \sim M\)) peaks at \(t_{\text{peak}}\), often dominating the total light curve unless heating rates or opacities vary significantly with depth (e.g., due to free neutrons in outer layers). As the ejecta expands, the radius of each mass layer evolves simply as: \[\frac{dR_v}{dt} = v,\] reflecting free expansion. The thermal energy \(\delta E_v\) of a layer with mass \(\delta M_v\) evolves according to

where the first term represents adiabatic losses due to expansion, \(L_v\) is the luminosity, and \(\dot{Q}\) accounts for heating sources. The luminosity is given by,

where \(t_{lc,v} = R_v / c\) (the light-crossing time) ensures energy loss does not exceed physically realistic limits, particularly when the layer becomes optically thin. Heating arises from multiple mechanisms

including radioactive decay (\(\dot{Q}_{r,v}\)), a millisecond magnetar (\(\dot{Q}_{\text{mag}}\)), or fallback accretion (\(\dot{Q}_{\text{fb}}\)). For radioactive heating, the energy input is typically smaller than the initial kinetic energy, allowing the initial power-law distribution to hold. However, significant energy contribution from a central engine restricts the free expansion. In such cases, the innermost shell’s velocity \(v_0\) (mass \(M\), energy \(E_{v_0}\)) evolves as,

where \(R_0\) is the inner shell’s radius, and the right-hand side balances adiabatic losses, neglecting relativistic effects (at \(M_\rm{ej}\) \(\lesssim 10^{-2} \ M_\odot\)) or secondary shocks in outer layers for simplicity. Assuming blackbody emission, the effective temperature is

where \(L_{\text{tot}} = \int L_v \frac{dM_v}{dv} dv\) sums the luminosity across all layers, and \(R_{\text{ph}}\) is the photosphere radius, defined where optical depth \(\tau_v = 1\). The spectral flux density at frequency \(\nu\) is then given by,

where \(D\) is the distance to the source, and cosmological corrections are omitted for simplicity.

3Radiative Transfer

While an analytical toy model from Sec. 2.3 offer valuable insights into the basic physical principles governing kilonovae, they rely on simplifying assumptions that limit their applicability to real events. In reality, modeling kilonovae poses formidable challenges due to the intricate physics of their ejecta. The material expelled during a neutron star merger is neither uniform nor spherically symmetric, often comprising regions with diverse compositions—some rich in lanthanides and actinides, others relatively free of these heavy elements. These compositional variations lead to significant differences in opacity, influencing the radiation. Furthermore, the ejecta expand rapidly, necessitating models that account for time-dependent changes in density, temperature, and chemical makeup. To study the observable signatures of kilonovae—such as spectra, light curves, and polarization—radiative transfer simulations are important. These simulations trace the journey of photons through the ejecta, bridging the gap between theoretical predictions and observed data collected.

At the core of radiative transfer simulations is the radiative transfer equation (RTE), which describes how photons interact with matter as they traverse,

where \(I_\nu\) is the specific intensity of radiation at frequency \(\nu\), \(j_\nu\) is the emission coefficient, \(\alpha_\nu\) is the absorption coefficient, and \(s\) is the distance traveled by the photons. Solving the RTE analytically is often impractical, particularly in optically thick media or complex geometries, so numerical methods are employed. One approach is the Monte Carlo method, which simulates the trajectories of many individual photons or photon packets. It tracks their interactions—absorption, scattering, or transmission—with the medium and uses a large sample to build a statistical picture of the radiation field. This method works well for three-dimensional setups and frequency-dependent opacities but requires significant computational resources.

In this work, we primarily use the 3-dimensional Monte Carlo Radiative Transfer (MCRT) code, POSSIS [64] and 1-dimensional models from SEDONA [35].

3.1POSSIS

POSSIS is a radiative transfer code which leverages the Monte Carlo method as its computational backbone, a probabilistic approach ideally suited to study the radiation transport in kilonova ejecta. In this technique, the photon packets are characterized by frequency, direction, and position, as they traverse through the expanding medium. These packets originate from the thermal energy pool, driven primarily by radioactive decay, undergoing scattering, absorption, or re-emission based on local physical conditions. By tracking individual photon histories, POSSIS preserves the full complexity of the system, enabling detailed predictions of how emission varies with viewing angles, time, and frequency.

The physical foundation of POSSIS begins with the ejecta model typically a 3-dimensional grid from Numerical Relativity simulations, which describes the structure and evolution of the material expelled during the neutron star merger. This geometry is setup in such that, we have a lanthanide rich component in the equatorial region and lanthanide free component around in the polar region.

A key assumption is homologous expansion, where the velocity of each ejecta element scales linearly with its radial distance from the explosion’s center, expressed as \(v \propto r\). This approximation reflects the rapid, free expansion following the merger’s dynamical phase. As the ejecta expands, the density scales as,

where \(\rho_{ij}\) is the density at grid cell \(i\) at time \(t_j\), \(\rho_{i,0}\) is the initial density, \(t_0\) is the initial time, and \(t_j\) is the time step. The factor of 3 arises from the volumetric expansion of the ejecta, directly influencing the optical depth and the timescales over which radiation escapes. As already established in Sec. 2.2, the ejecta’s composition is crucial, with lanthanides and actinides playing a significant role in shaping the kilonova’s emission. POSSIS accounts for this by incorporating several opacity sources, including bound-bound (lanthanide rich: \(Y_e \leq 0.25\), lanthanide free: \(Y_e >0.25\)), electron scattering, bound-free, and free-free opacities. Bound-bound opacities, \(\kappa_{\rm bb}\), govern the interaction of radiation with spectral lines. Electron scattering opacity, \(\kappa_{\rm es}\), is critical for polarization calculations, while bound-free \(\kappa_{\rm bf}\) and free-free \(\kappa_{\rm ff}\) contribute to continuum absorption. Opacities may be specified at each time step or interpolated from a reference time using an appropriate functional form (based on atomic calculations from Ref. [142]), facilitating efficient modeling of time-dependent spectra.

At each simulation time step, a predetermined number of photon packets is generated, each characterized by an initial position, propagation direction, energy, frequency, and polarization state. The frequency of a packet is assigned by sampling the thermal emissivity function,

where \(B(\nu, T)\) is the Planck function and \(\kappa_{\rm tot}(\nu)\) is the total opacity. The photon packets propagate through the ejecta, interacting with the medium via absorption, re-emission, and scattering processes. These interactions are computed in the comoving frame, with appropriate transformations applied to ensure proper treatment of relativistic effects.

Photon packets travel through the ejecta until they reach the computational boundary, where they are either absorbed, scattered, or allowed to escape. The propagation path of each photon is determined probabilistically based on the mean free path associated with the local opacity. In the case of electron scattering, the photon direction is altered while its frequency remains unchanged. If absorption occurs, the packet is re-emitted with a new frequency determined by the thermal emissivity. This iterative process continues until all photon packets either escape the system or are fully absorbed.

3.1.1Synthetic Observables

To extract synthetic observables, POSSIS employs two complementary techniques: the direct counting technique (DCT) and the event-based technique (EBT). In the DCT approach, escaping photon packets are binned according to their final propagation directions to generate spectra and light curves. The EBT, in contrast, assigns probabilities to virtual packets at each interaction site, weighting them according to their likelihood of reaching the observer. This method significantly reduces Monte Carlo noise and improves the accuracy of polarization predictions.

The synthetic observables generated by POSSIS include spectral energy distributions (SEDs), broadband light curves, and polarization spectra. The SEDs provide wavelength-dependent flux predictions at different epochs and viewing angles, revealing the interaction between the ejecta’s composition and its radiative properties. Broadband light curves, computed in various photometric filters, enable direct comparisons with observational data from transient surveys. Polarization spectra further constrain ejecta asymmetries by quantifying the anisotropic scattering of photons.

Broadband light curves exhibit significant viewing-angle dependence, reflecting the asymmetric nature of the ejecta. When viewed along the polar axis, the light curves appear brighter and decay more rapidly, whereas equatorial views lead to fainter, longer-lasting emission due to increased opacity from lanthanide-rich material. This anisotropy is particularly pronounced in kilonovae, where distinct ejecta components contribute differently across various wavelengths. The analysis of broadband light curves allows for the inference of key parameters such as the total ejecta mass, lanthanide fraction, and inclination angle relative to the observer.

3.2SEDONA

Although SEDONA is a 3D code, in this thesis, we focus on a grid of one-dimensional spherically symmetric models presented in Ref. [35]. This model is parameterized by the ejecta total mass \(M_{\rm ej}\), average velocity \(v_{\rm ej}\), and lanthanide fraction \(X_\rm{lan}\). The bulk of the freely expanding ejecta is determined by the ejecta mass \(M_{\rm ej}\). The density profile of the ejecta is described using a broken power-law that transits from the gradually declining interior with \(v_{\rm ej}t / r\) to the steeply dropping outer layer with \((v_{\rm ej}t / r)^{10}\). Finally, the lanthanide fraction \(X_\rm{lan}\) influences the opacity and color evolution of the kilonova, where larger lanthanide fractions result in increased opacity and longer-duration emissions shifted towards the infrared.

The energy source powering the kilonova emission is the radioactive decay of r-process nuclei synthesized during the merger. The heating rate follows a power-law which is very similar to the on defined in Eq. 2.6, \(\dot{\epsilon}_{\rm nuc} \propto t^{-1.3}\), due to the wide distribution of half-lives of these freshly formed isotopes. Consequently, the radiation escapes once the diffusion timescale becomes comparable to the expansion time, and the time-dependent luminosity is primarily governed by the product \(\kappa M / v\), where \(\kappa\) is the opacity. This opacity is highly sensitive to the composition—lanthanide-rich ejecta have significantly higher opacities (up to \(\sim\)10 cm\(^2\) g\(^{-1}\)) than lanthanide-poor material (\(<\)1 cm\(^2\) g\(^{-1}\)), due to the dense forest of line transitions from f-shell electrons in heavy elements.

As a result, low-\(X_{\rm lan}\) (\(\leq 10^{-4}\)) material emits briefly and in the optical (a “blue” kilonova), while high-\(X_{\rm lan}\) (\(\geq 10^{-2}\)) ejecta lead to longer, redder, infrared-dominated emission (a “red” kilonova). In AT2017gfo, both components are observed, the early-time optical light is attributed to fast (\(\sim0.3c\)), lanthanide-poor ejecta likely from polar shock-driven outflows; the longer infrared tail arises from slower (\(\sim0.1-0.2c\)), lanthanide-rich material, possibly originating from disk winds or tidal ejecta.

4Statistical Methods

Repeated observations of astronomical phenomena are often characterized by complex, noisy data that require sophisticated statistical techniques for their analysis. In this chapter, we introduce the fundamental concepts of probability theory and Bayesian statistics, which form the basis of our approach to parameter estimation and model comparison in multi-messenger astronomy. We also discuss the computational techniques required for practical Bayesian inference in complex astrophysical problems.

4.1Bayesian Statistics

4.1.1Probability Theory

Probability theory forms the mathematical backbone of statistical inference. A probability distribution describes how probability mass or density is assigned to possible values of a random variable. For a discrete random variable \(X\), the probability of an event \(A\) occurring is given by,

where \(A_i\) are mutually exclusive outcomes. For a continuous variable \(\theta\), the probability density function (PDF) is given by,

One of the fundamental properties of probability theory is the law of total probability, which states that for any partition \(B_i\) of the sample space, the probability of an event \(A\) is given by,

Another key principle is conditional probability, which quantifies the probability of an event \(A\) given that another event \(B\) has occurred, given by

In the context of astrophysical inference, we are often interested in estimating the parameters \(\theta\) of a model, given some observed data \(D\). The interpretation of probability in such problems can be broadly categorized into two perspectives:

• Frequentist Interpretation: Probability is defined as the long-run relative frequency of an event in repeated trials. This approach underlies classical statistical methods, including confidence intervals and hypothesis testing, but does not incorporate prior knowledge in parameter estimation.

• Bayesian interpretation: Probability represents a degree of belief in an event, updated as new evidence is obtained. This approach allows for the incorporation of prior information through Bayes’ theorem, making it particularly suited for problems where data are limited or where prior knowledge is available.

In astrophysical applications, particularly in transient event analysis, the frequentist approach is often inadequate because no two astronomical events are strictly identical. Consequently, Bayesian inference, which systematically updates prior knowledge with new observations, provides a more robust framework for parameter estimation in such scenarios.

4.1.2Bayes’ Theorem

Bayes’ theorem is derived directly from the definition of conditional probability. For two events \(A\) and \(B\), Eq. 4.4 can be rewritten as

which can be rearranged to yield Bayes’ theorem as follows,

where each term carries a specific statistical interpretation. The prior probability, \(P(A)\), quantifies the initial probability of event \(A\) based on pre-existing information or assumptions, established prior to observing the specific event. The likelihood, \(P(B | A)\) quantifies the probability of observing data \(B\) given that the hypothesis \(A\) is true. The posterior probability, \(P(B|A)\) is the updated probability of \(A\) after accounting for the data \(B\). The evidence, or marginal likelihood, \(P(B)\), acts as a normalization factor, ensuring that the posterior is a valid probability distribution. For discrete hypothesis, it is calculated as

summing over all possible states \(A_i\) and for continuous hypothesis, it is given by

where \(\theta\) denotes the parameter vector and the integral spans the prior distribution.

This is the fundamental form of Bayes’ theorem, which allows the inversion of conditional probabilities based on prior knowledge. The real power of Bayes’ theorem lies in its systematic way to update our beliefs about a parameter as more data become available. This makes it particularly powerful in astrophysical applications, where observations are often sequentially obtained over time, such as in the case of transient events.

4.1.3Parameter Estimation

To perform parameter estimation in a Bayesian framework, we seek to infer the posterior distribution of the model parameters \(\vec\theta\) given the observed data \(d\) under the assumption of some hypothesis \(\mathcal{H}\). This is achieved by combining the likelihood function, which quantifies the probability of observing the data given the model parameters, with the prior distribution, which encodes our prior knowledge about the parameters. To this end, Bayes’ theorem takes the form

or more familiarly as

and the evidence is given by

where integration is performed over the entire parameter space \(\Omega_{\vec\theta}\).

In this context, the hypothesis \(\mathcal{H}\) represents the underlying model or theoretical framework that we assume to be true when performing our analysis. It encapsulates our fundamental assumptions about the system under study and defines the space of possible parameter values \(\vec\theta\) that can describe the data. For instance, in astronomical observations, \(\mathcal{H}\) might represent the hypothesis that Newton’s law of gravitation correctly describes the motion of celestial bodies, or alternatively, that Einstein’s theory of general relativity is the correct framework. The explicit conditioning on \(\mathcal{H}\) throughout the equations reminds us that all our inferences about parameters \(\vec\theta\) are conditional on the validity of the underlying model. The posterior distribution \(P(\vec\theta | d,\mathcal{H})\) provides the probability of different parameter values given both the observed data and our assumed model. Similarly, the prior \(P(\vec\theta | \mathcal{H})\) encodes our beliefs about parameter values within the context of hypothesis \(\mathcal{H}\) before observing any data. The evidence, cf. Eq. 4.11, also known as the marginal likelihood, quantifies how well the hypothesis \(\mathcal{H}\) explains the observed data, averaged over all possible parameter values weighted by their prior probabilities.

4.1.4Hypothesis Testing

In Bayesian hypothesis testing, we extend the framework of parameter estimation to evaluate the relative plausibility of competing hypotheses given observed data. Unlike parameter estimation, where we infer the distribution of parameters within a single model, hypothesis testing involves comparing multiple models or hypotheses to determine which one provides a better explanation for the observed phenomena.

Consider two competing hypotheses, \(\mathcal{H}_1\) and \(\mathcal{H}_2\), each representing different models or theoretical frameworks. Eq. 4.9 can be transformed to represent the posterior probability of hypothesis \(\mathcal{H}_i\) given the data \(d\) as follows

where \(P(d | \mathcal{H}_i)\) is the evidence for hypothesis \(\mathcal{H}_i\), \(P(\mathcal{H}_i)\) is the prior probability assigned to hypothesis \(\mathcal{H}_i\) before observing the data, and \(P(d)\) is the total probability of observing the data across all possible hypotheses.

The ratio of posterior probabilities for two competing hypotheses, known as the posterior odds or odds ratio, provides a direct measure of their relative plausibility

The first term on the right-hand side is the Bayes factor, denoted as \(\mathcal{B}^{1}_{2}\), which quantifies the relative evidence provided by the data in favor of \(\mathcal{H}_1\) over \(\mathcal{H}_2\)

where \(\mathcal{Z}_i\) represents the evidence for hypothesis \(\mathcal{H}_i\) as defined in Eq. 4.11. The second term in Eq. 4.13 is the prior odds, which represents the relative plausibility of the hypotheses before observing the data,

With Eq. 4.14 and Eq. 4.15, the odds ratio can be expressed as,

In practical applications, the prior odds are often set to unity under the assumption that both hypotheses are equally likely, hence the posterior odds are same as the Bayes factor.

The Bayes factor \(\mathcal{B}^{1}_{2}\) provides an objective measure of the strength of evidence in the data, independent of prior beliefs about the hypotheses.

When the Bayes factor \(\mathcal{B}^{1}_{2} > 1\), the data provides evidence in favor of \(\mathcal{H}_1\) over \(\mathcal{H}_2\), with larger values indicating stronger evidence. Conversely, when \(\mathcal{B}^{1}_{2} < 1\), the data supports \(\mathcal{H}_2\) over \(\mathcal{H}_1\). The magnitude of the Bayes factor can be interpreted according to established scales, with values above 10 typically considered strong evidence and values above 100 decisive evidence.

For nested models, where hypothesis \(\mathcal{H}_1\) represents a more general model that includes \(\mathcal{H}_2\) as a special case, the Bayes factor automatically implements Occam’s razor by penalizing unnecessary complexity. If the additional parameters in \(\mathcal{H}_1\) do not significantly improve the fit to the data, the evidence will favor the simpler model \(\mathcal{H}_2\), as expressed by

When dealing with multiple hypotheses \(\{\mathcal{H}_1, \mathcal{H}_2, \ldots, \mathcal{H}_n\}\), we can compute the posterior probability of each hypothesis as

This formulation allows for a coherent assessment of the relative plausibility of all competing hypotheses simultaneously, rather than through pairwise comparisons alone.

4.2Computational Methods

The Bayesian approach to parameter estimation and hypothesis testing, while theoretically elegant, presents significant computational challenges in practical applications. In particular, evaluating the posterior distribution, cf. Eq. 4.10, requires computing the multidimensional integral Eq. 4.11 for the evidence \(\mathcal{Z}\). For all but the simplest problems with analytic solutions, this integral is intractable using traditional numerical integration methods, especially as the dimensionality of the parameter space increases—a phenomenon known as the "curse of dimensionality".

In astrophysical applications, we typically encounter complex likelihood functions that may be multimodal, exhibit strong parameter degeneracies, or contain sharp peaks in high-dimensional spaces. For instance, in multi-messenger astronomy, the joint analysis of gravitational wave and electromagnetic signals involves parameter spaces with dimensions often exceeding 10-15 parameters. Direct numerical integration in such spaces would require an exponentially increasing number of function evaluations that quickly becomes computationally expensive, more often than not, prohibitive.

Consequently, stochastic sampling methods have emerged as the preferred approach for Bayesian computation for astrophysical problems. These methods approximate the posterior distribution by drawing samples from it, without requiring the explicit calculation of the normalizing evidence. The generated samples can then be used to estimate parameter values, construct credible intervals, and compute quantities of interest. Additionally, some sampling methods provide estimates of the evidence \(\mathcal{Z}\), enabling rigorous model comparison through Bayes factors as described in Eq. 4.14.

In the following texts, we discuss two principal sampling techniques that have proven particularly effective for astronomical applications: Markov Chain Monte Carlo (MCMC) methods, which excel at efficiently exploring posterior distributions for parameter estimation, and nested sampling, which was specifically designed to compute the evidence while simultaneously providing posterior samples.

4.2.1Markov Chain Monte Carlo

MCMC methods constitute a class of algorithms designed to draw samples from a probability distribution by constructing a Markov chain whose equilibrium distribution converges to the target posterior distribution. The fundamental principle underlying MCMC is that, after an initial "burn-in" period, the states of the chain represent correlated samples from the posterior distribution, which can be used to estimate any quantity of interest.

The Metropolis-Hastings algorithm [144,145], a canonical MCMC method, generates samples from the posterior distribution, \(P(\vec\theta | d, \mathcal{H})\), through the following procedure:

For symmetric proposal distributions, where \(q(\vec\theta' | \vec\theta_t) = q(\vec\theta_t | \vec\theta')\), this simplifies to the Metropolis algorithm with acceptance ratio,

A crucial aspect of the Metropolis-Hastings algorithm is that it only requires evaluating the posterior up to a normalizing constant, as the ratio in Eq. 4.20 cancels out the evidence term \(\mathcal{Z}\)

While conceptually straightforward, the practical performance of MCMC methods depends critically on the choice of proposal distribution. Proposals that are too conservative (small step sizes) result in high acceptance rates but slow exploration of the parameter space, whereas proposals that are too ambitious (large step sizes) lead to many rejections and inefficient sampling. This challenge is especially pronounced in high-dimensional spaces with strong parameter correlations, which are common in astrophysical applications.

A critical aspect of MCMC implementation is the handling of burn-in samples. Since the Markov chain typically starts from an arbitrary position in parameter space, the initial samples do not accurately represent the target posterior distribution. These early iterations, known as the burn-in phase, must be discarded to ensure that only samples from the converged stationary distribution are used for inference. Formally, if \(B\) represents the number of burn-in iterations, the samples used for analysis are \({\vec\theta_{B+1}, \vec\theta_{B+2}, \ldots, \vec\theta_N}\).

The primary output of MCMC methods, after discarding burn-in samples, is a set of samples that approximates the posterior distribution. From these samples, we can estimate moments, compute credible intervals, and construct marginal posterior distributions. The posterior mean for each parameter \(\theta_i\) is estimated as,

Similarly, the posterior covariance between parameters \(\theta_i\) and \(\theta_j\) is estimated as,

4.2.2Nested Sampling

While MCMC methods are well-suited for exploring the posterior distribution for parameter estimation, they are less efficient at computing the evidence \(\mathcal{Z}\), which is essential for model comparison. nested sampling, introduced by Skilling [146], is a computational method specifically designed to efficiently calculate the evidence \(\mathcal{Z}\) while simultaneously generating samples from the posterior distribution. Unlike MCMC methods, which focus primarily on exploring the posterior, nested sampling deliberately targets the evidence integral in Equation [eq:evidence]. The key idea of nested sampling is to transform the multidimensional evidence integral over parameter space, \(\Omega_{\vec\theta}\), into a one-dimensional integral over the prior volume. We define the prior volume \(X(\lambda)\) as the fraction of prior volume enclosed by the iso-likelihood contour \(\mathcal{L}(\vec\theta) = \lambda\),

As \(\lambda\) increases from 0 to the maximum likelihood value, \(X(\lambda)\) decreases monotonically from 1 to 0. This allows us to rewrite the evidence integral as a one-dimensional integral over the prior volume,

where \(\mathcal{L}(X)\) is the inverse function of \(X(\lambda)\), representing the likelihood value that corresponds to prior volume \(X\).

Stopping Criterion

Since nested sampling proceeds by iteratively removing the point with the lowest likelihood and replacing it with a new sample, the algorithm must determine when to terminate. A common stopping criterion is based on the remaining prior volume contributing negligibly to the evidence estimate. Specifically, the algorithm terminates when the change in evidence between successive iterations falls below a specified threshold. At a given iteration \(i\) we compute an estimate \(\hat{Z_i}\) of the accumulated evidence and compare it with the remaining evidence \(\Delta \hat{Z_i}\). However, this remaining evidence is unknown, but one can put strict upper bounds on it by considering the maximum likelihood ever encountered in the entire parameter space and the remaining prior volume \(X_i\). Under this assumption, the remaining evidence can be bounded as,

With this bound, the stopping criterion can be expressed as the log ratio of evidences,

Evidence and Posterior

Once the nested sampling algorithm terminates, we can estimate the evidence \(\mathcal{Z}\) by integrating the likelihood values of the removed dead points weighted by the change in prior volume,

Since, the entire discussion is based on how unfeasible it is to compute the integral, we can approximate the integral using the trapezoidal sum as follows,

The posterior distribution can be estimated from the dead points by weighing each sample by the importance weight \(p_i\),

5Nuclear-Physics and Multi-Messenger Astrophysics

Building upon the foundations established in the preceding chapters, we have developed a comprehensive understanding of the physical mechanisms governing BNS mergers, the EM signatures they produce, and the statistical frameworks required to extract meaningful constraints from observational data.

All of these components converge in the NMMA framework, which serves as an advanced Bayesian inference tool for the simultaneous analysis of gravitational wave signals, kilonova light curves, and other multi-messenger data. NMMA synthesizes information from different observational channels while incorporating nuclear-physics constraints and gravitational wave data to provide a self-consistent framework for understanding BNS mergers.

5.1Overview

The NMMA framework is designed to address several fundamental problems in neutron star astrophysics. The equation of state (EOS) of neutron star matter remains one of the central unresolved questions, as it governs the internal structure and dynamics of these compact objects. Constraints on the EOS can be obtained from gravitational-wave observations of the tidal deformations experienced by neutron stars during inspiral, from X-ray and radio measurements of pulsars, and from nuclear theory calculations. Similarly, the properties of BNS and neutron star-black hole (NSBH) mergers, including their ejecta masses, velocities, and composition can be constrained by analyzing the EM signals emitted following the merger. Additionally, multi-messenger observations of neutron star mergers provide an independent means of measuring the Hubble constant, thus offering a novel probe of cosmology. The NMMA framework integrates all of these constraints into a unified inference pipeline, ensuring that information from different messengers is combined in a statistically rigorous manner.

At its core, NMMA is based on Bayesian inference (through BILBY [148,149]), allowing for a consistent incorporation of prior knowledge, observational data and physical measurements to estimate posterior probability distributions of physical parameters. The framework primarily uses nested sampling to efficiently explore the high-dimensional parameter space of BNS mergers. The likelihood function within NMMA is constructed as a product of the likelihoods associated with different observational channels: gravitational waves, kilonova light curves, and, when available, GRB afterglows. This multi-component likelihood function enables NMMA to leverage all available information to obtain more precise constraints on neutron star properties and merger physics.

5.1.1Nuclear Physics Constraints

NMMA incorporates nuclear physics constraints into the inference process to ensure that the derived properties of neutron stars and their mergers are physically consistent with our understanding of dense nuclear matter. The EOS of neutron star matter, is a critical component that influences both gravitational wave signals through tidal deformabilities and kilonova emission through the properties of merger ejecta. The framework employs a two-region approach to model the EOS. At low densities (up to approximately \(2n_{\rm sat}\), where \(n_{\rm sat} = 0.16\ \rm{fm}^{-3}\) is the nuclear saturation density), NMMA utilizes Quantum Monte Carlo calculations based on Chiral Effective Field Theory (EFT), providing a systematic description of nuclear interactions with quantifiable uncertainties. For the higher densities reached in neutron star cores, where Chiral EFT becomes unreliable, NMMA implements a model independent speed-of-sound extension scheme bounded only by causality (\(c_s \leq c\)) and mechanical stability (\(c_s \geq 0\)).¹Refer to and of Ref. [150] and of Ref. [151] for a detailed discussion.

This methodology allows NMMA to generate EOS realizations that are consistent with known nuclear physics while adequately exploring uncertainties at high densities. For each sampled EOS, the framework solves the TOV equations to determine neutron star properties such as the mass-radius relationship and tidal deformabilities, which directly influence observable gravitational wave and EM signatures.

5.1.2Gravitational Waves

Gravitational-wave data play a fundamental role in NMMA’s inference pipeline. The inspiral phase of a BNS merger encodes crucial information about the component masses and tidal deformabilities, which in turn provide constraints on the EOS. To extract this information, NMMA utilizes various waveform models from LALSuite [152]. The likelihood function for gravitational-wave data is computed by cross-correlating the observed detector strain with theoretical waveform templates, weighted by the detector’s noise power spectral density. The likelihood \(\mathcal{L}_{\rm GW}\) for gravitational wave data, given a set of parameters \(\theta\), is defined as,

where \(d\) is the detector strain data, \(h(\theta)\) is the GW signal generated from parameters \(\theta\), and \(\langle \cdot \mid \cdot \rangle\) is the inner product defined as,

where \(\tilde{a}(f)\) and \(\tilde{b}(f)\) are the Fourier transforms of \(a(t)\) and \(b(t)\), respectively, \(S_{n}(f)\) is the one-sided noise power spectral density of the noise, and \(*\) is the complex conjugate. The integration limits \(f_{\mathrm{min}}\) and \(f_{\mathrm{max}}\) are chosen to match the frequency range of the signal.

5.1.3Kilonovae and Gamma-Ray Bursts

Beyond gravitational waves, NMMA incorporates EM data from kilonovae and GRB afterglows to further constrain the properties of neutron star mergers. To model kilonovae, NMMA primarily employs POSSIS. Additionally, NMMA models the GRB afterglows using afterglowpy [153,154], a semi-analytic code that simulates the synchrotron emission from relativistic jets produced in neutron star mergers. The likelihood \(\mathcal{L}_\rm{EM}\) for kilonova and GRB can be combined as a single EM likelihood, given a set of parameters \(\theta\) as

where \(m^j_i\) is the AB magnitude in filter \(j\) at time \(i\) with the corresponding measurement error \(\sigma^j_i \equiv \sigma^j(t_i)\), \(m_i^{j, \text { est }}(\vec{\theta})\) is the estimated AB magnitude for the source parameters \(\vec{\theta}\) (e.g., ejecta masses, velocities) from the model and \(\sigma_{\mathrm{sys}}\) is the additional uncertainty to account for systematic errors in the kilonova and afterglow modelling. Computing likelihood in this form is equivalent to including an additional shift to the light curve by \(\Delta m\), and marginalizing it with a normal distribution with a mean of 0 and variance \(\sigma_{\mathrm{sys}}\).

Using Eq. 5.1 and Eq. 5.3, NMMA combines GW and EM to obtain a joint likelihood function,

in order to capture the full multi-messenger picture.

5.1.4Measurements from Isolated Neutron Stars

While the primary focus of NMMA is on BNS and NSBH mergers, the framework also incorporates information from isolated neutron stars to further constrain the EOS and priors on neutron star properties. This information includes measurements of neutron star masses and radii from NICER X-ray observations, as well as radio pulsar timing data.²Refer to and of Ref. [150], of Ref. [151], and Ref. [155] for a detailed discussion.

5.2Surrogate Modelling

Radiative transfer simulations, such as those performed by POSSIS and SEDONA (hereafter referred to interchangeably as Bu2019lm and Ka2017, respectively, in order to be coherent with their names within the NMMA community), are computationally expensive, making them impractical to run repeatedly for any arbitrary set of parameters during the Bayesian inference sampling. To address this challenge, NMMA employ surrogate models [156] that efficiently interpolates between pre-computed simulation grids in order to reduce the memory footprint and overall computational cost of the inference process. To this end, NMMA employs two step approach, first reducing the dimensionality of the data using singular value decomposition (SVD), and then training a surrogate model either using Gaussian Process Regression (GPR) or Neural Network (NN) to interpolate between the reduced data points.

5.2.1Dimensionality Reduction

While one can directly train a machine learning model on the observables, it is much efficient to first perform some kind of dimensionality reduction in order to get a more compact representation of the data. This is achieved using SVD. For a given set of light curves \(\mathbb{L}^{j}\), generated from a set of parameters \({\theta}^j\), the full light curve matrix is represented as \(\mathbb{L}^j_i \equiv [\mathbb{L}_i(\theta^j)]\), where \(i\) indexes the time steps of the light curve. The singular value decomposition of \(\mathbb{L}^j_i\) is performed first by normalizing the light curves using MinMax Normalisation and then decomposing the matrix as

where \(U\) are the left singular vectors, \(\Sigma\) is the diagonal matrix of singular values, and \(V\) are the right singular vectors. With these new basis vectors, the original light curve \(\mathbb{L}\) can be projected onto the left singular vectors

where \(s_k(\theta^j)\) are the weights of the principal components of the input light curve matrix \(\mathbb{L}^j_i\). In general, this allows for a more compact representation of the light curves, such that, in NMMA only first 10 basis vectors are enough to capture the majority of the information in the light curves.

5.2.2Interpolation Methods

Once the dimensionality of the data has been reduced, the next step is to train a surrogate model that can interpolate between the pre-computed simulation grids. However, in cases such as, the requested time interval is larger than the original data, a simpler linear interpolation is used.

5.2.2.1 Gaussian Process Regression

GPR provides a non-parametric method for estimating smooth functions while incorporating uncertainty quantification. A Gaussian process is defined by a mean function \(m(\theta)\) and a covariance function \(k(\theta, \theta')\), such that the function values at any set of points \(\theta\) are jointly Gaussian distributed such that the coefficients \(s_k\) can be represented as

Consequently, the coefficients \(s_k\) are then interpolated with a rational-quadratic kernel [157]

where \(\vec{\theta}\) and \(\vec{\theta}'\) are vectors of input parameters and hyperparameters \(\alpha\) and \(l\) are chosen by maximizing the evidence for the data under a zero-mean Gaussian process. Post interpolation, the full light curve can be reconstructed by projected the interpolated coefficients back onto the original basis vectors

NMMA uses scikit-learn [158] to implement the GPR interpolation.

5.2.2.2Neural Networks

Following Ref. [159], NMMA has implemented a NN based interpolation method to reduce memory consumption and accelerate inference compared to the GPR implementation. The primary advantage of this approach is its ability to efficiently handle high-dimensional kilonova models as opposed to GPR, which tends to breakdown with more than 4 model parameters.³Refer Ref. [159] for a detailed discussion on NN accelerated kilonova surrogate modeling.

Similar to GPR, the NN is trained on a reduced set of data points and then used to interpolate between them. The neural network follows a feed-forward architecture implemented using TensorFlow [160] with Adam optimizer [161], a mean squared error loss function and Rectified linear unit (ReLU) [162] activation function. Once trained, the NN can rapidly generate interpolated values across the parameter space, significantly reducing the computational cost associated with traditional interpolation methods. Similar to the GPR approach, the full light curve can be reconstructed by performing inverse transformations on the interpolated coefficients.

The specific architecture, hyperparameter choices, and further implementation details of the neural networks used in this work are discussed below.

POSSIS

The NN surrogate model is trained on 1596 parameter combinations and their corresponding light curves generated from the POSSIS simulation. The NN architecture consists of three fully connected hidden layers with 128, 256, and 128 neurons, respectively. The input layer has neurons corresponding the total number model parameters \((\Phi\), \(\iota\), \({\log_{10}M^{\rm{wind}}_{\rm{ej}}}\) and \({\log_{10}M^{\rm{dyn}}_{\rm{ej}}}\), cf.Sec. 3.1). The output layer consists of neurons corresponding to the number of principal components, which is set to 10 in our case. The model is trained for 100 epochs with a batch size of 128, using a learning rate of \(10^{-3}\). To ensure numerical stability and consistency in the output range, output is rescaled using a MinMax scaler, which normalizes them to a fixed range before applying the inverse transformation to reconstruct the light curves. Additionally, 20% of the dataset is set aside for validation to monitor the model’s generalization performance and prevent overfitting.

SEDONA

The total number of publicly available⁴https://github.com/dnkasen/Kasen_Kilonova_Models_2017 light curves for SEDONA is relatively small, comprising only 329 parameter combinations. Despite the limited dataset, the pre-training steps remain the same as those used for POSSIS. However, due to the reduced data volume, a different neural network architecture is adopted.

Here we used a single hidden layer of 2048 neurons, with input a layer of three neurons corresponding to the model parameters (\({\log_{10}v_{\rm{ej}}}\), \({\log_{10}M_{\rm{ej}}}\), and \({\log_{10}X_{\rm{lan}}},\) cf.Sec. 3.2). As in the case of POSSIS, the output layer consists of 10 neurons, corresponding to the number of principal components used to reduce the dimensionality of the light curves. To improve generalization and mitigate overfitting, a dropout layer with a rate of 0.6 is included during training.

5.2.3Neural Network Benchmarking

Since the two radiative transfer models, POSSIS and SEDONA, have fundamentally different geometries and are trained on distinct datasets, their corresponding NN surrogate models are also distinct. This difference implies that the performance of the two models should not be directly compared without careful consideration of their respective training datasets and methodologies.

To assess the accuracy of the NN surrogate models, we use the root mean square error (RMSE) as a metric to quantify the average deviation between the predicted and actual light curves. The RMSE is defined as

where \(y_i\) is the actual value obtained through the radiative-transfer simulation at time \(t_i\), \(\hat{y}_i\) denotes the prediction from the surrogate model, and \(N\) is the total number of time points.

We compute this metric across the entire parameter grid and all photometric filters. The median RMSE values are found to be \(0.103~\rm mag\) for Bu2019lm (computed over a 14-day period) and \(0.485~\rm mag\) for Ka2017 (computed over a 7-day period). The rationale behind selecting these specific time intervals will be discussed in the next sections.

6Systematic Uncertainties

Despite the robustness of NMMA’s inference methodology, one of the key challenges in kilonova modeling is the presence of uncertainties, which arise from factors such as nuclear heating rates, opacities, and ejecta morphology, cf. Ch. 2, all of which evolve over time. Historically, the NMMA [151] framework has addressed these uncertainties by applying a fixed systematic uncertainty of \(1\) magnitude [19,151,163,41], a choice validated by studies demonstrating consistency across different ejecta morphologies [67]. However, assuming a fixed uncertainty in the kilonova likelihood function (\(\sigma_{\rm sys}\) in Eq. 5.3) can lead to biases in parameter estimation.

More recent approaches have sought to refine this treatment by allowing systematic uncertainty to be sampled as a free parameter [3,164]¹Ref. [3] used the freely sampled systematic error treatment as implemented in commit 32aadf8., with some works incorporating filter-dependent priors [165]. However, systematic uncertainty in kilonova models is not solely filter-dependent but also evolves over time, reflecting changes in opacity descriptions at early epochs [142,74] and non-local thermal equilibrium effects that become significant after about one week [166,167]. To address this limitation, this thesis and the corresponding paper introduce a time-dependent error model, which accounts for variations in kilonova properties across different observational epochs. Additionally, the methodology allows for a combined treatment of filter- and time-dependent uncertainties, providing a more nuanced representation of systematic effects. Additionally, we also the treatment of systematic uncertainty as a sampling parameter. These modifications, already implemented in NMMA,²See PRs  383 and  387. is model-agnostic and can be applied to any EM model irrespective of the nature of the transient, e.g., also for gamma-ray-burst afterglows or supernovae.

6.1Implementation Details

6.1.1Freely Sampled Systematic Uncertainty

In this implementation, systematic uncertainty is treated as a free parameter, denoted by \({\sigma_{\text{sys}}}\), which is sampled from a predefined prior probability distribution rather than being fixed. This is expressed as

where \(P(\sigma_{\rm sys})\) represents the prior distribution. The choice of prior can vary depending on empirical constraints, with common choices including a uniform or log-uniform distribution

or

During inference, \({\sigma_{\text{sys}}}\) is sampled along with other parameters, and its effect is incorporated into the likelihood function as the quadrature sum of statistical and systematic uncertainties as illustrated in the denominator of Eq. 5.3. This allows for a more flexible treatment of systematic effects by marginalizing over \({\sigma_{\text{sys}}}\) in the posterior distribution,

where \(\vec\theta\) represents the set of model parameters, \(d\) is the data, and \(\mathcal{L}\) is the likelihood function.

By defining \({\sigma_{\text{sys}}}\) as a variable sampled from a prior probability distribution, we allow the inference process to adaptively account for systematic effects without fixing them a priori. This flexibility is critical when empirical constraints on systematic uncertainties are uncertain or context-dependent, as it permits the model to adaptively explore the plausible range of \({\sigma_{\text{sys}}}\) during sampling. As the posterior distribution, Eq. 6.4, integrates this nuisance parameter alongside the model parameters \(\vec\theta\), this approach also provides a mechanism to assess the goodness of fit by examining the marginalized behavior of \({\sigma_{\text{sys}}}\) [3].

Model selection

The inclusion of the systematic uncertainty parameter \({\sigma_{\text{sys}}}\) in the Bayesian evidence is crucial for obtaining an accurate and representative inference of model parameters, particularly in cases where systematic effects may significantly influence the data.

Incorporating \({\sigma_{\text{sys}}}\) into the Bayesian evidence ensures that all potential systematic effects are properly integrated into the analysis, thus reflecting a more complete description of the data. If \({\sigma_{\text{sys}}}\) were excluded, the evidence would effectively disregard the systematic uncertainties, leading to an incomplete or biased model of the data. By allowing \({\sigma_{\text{sys}}}\) to be sampled from a prior distribution, it is possible to marginalize over these uncertainties and evaluate the likelihood of the data under varying assumptions about the magnitude of systematic errors. This marginalization is particularly important in cases when the prior knowledge does not provide exact estimates for systematic uncertainties, but rather offers a range of plausible values.

Furthermore, \({\sigma_{\text{sys}}}\) influences the Bayesian model comparison by affecting the model evidence, which becomes important while determining the relative plausibility of competing models. As systematic uncertainties may significantly alter the interpretation of the data, their treatment in the evidence calculation directly impacts the model selection. Ignoring these uncertainties would undermine the validity of model comparisons and potentially lead to misleading conclusions about the best-fitting model. Hence, one may need to be cautious when interpreting the Bayesian evidence when using such a nuisance parameter.³Subtle foreshadowing regarding the analysis of AT2017gfo.

6.1.2Time-Dependent Systematic Uncertainty

Radiative transfer models have made significant progress in kilonova modeling, yet discrepancies persist, particularly at very early times (<1 day) and late times (>1 week). At early times, the ejecta is extremely hot, dense, and expanding rapidly, leading to highly ionized atomic states for which opacity data remains incomplete [168,142]. Additionally, aspherical geometries complicate spectral modeling, while inefficient thermalization of radioactive decay products [78,74] affects early light curve predictions. These challenges hinder accurate parameter estimation from observed early-time spectra, such as those from AT2017gfo.

At later times, as the ejecta expands and cools, it transitions into a nebular phase where local thermodynamic equilibrium (LTE) approximations break down. The spectrum becomes dominated by line emission, requiring a full non-LTE treatment, which is computationally demanding [166]. Furthermore, evolving opacities due to the recombination of lanthanides and actinides introduce additional uncertainties in the predicted light curves [169]. Multi-component ejecta structures and viewing-angle dependencies further complicate light curve interpretations [170]. While recent advances, such as improved atomic databases and multi-dimensional radiative transfer models, have helped mitigate these issues, systematic uncertainties remain a critical challenge for inferring merger properties from kilonova observations.

To address these challenges, we introduce a time-dependent systematic uncertainty model that accounts for the evolving nature of kilonova properties across different epochs. Specifically, we define \({\sigma_{\text{sys}}}\) as a function of time, \({\sigma_{\text{sys}}}(t)\), allowing systematic effects to vary over the course of the kilonova event. The systematic uncertainty is modeled using a piecewise linear interpolation scheme, where \({\sigma_{\text{sys}}}(t)\) is sampled from a prior distribution at fixed time intervals during the inference process and then linearly interpolated at measurement times. Mathematically, this is expressed as

where \(n = 0, 1, ..., N - 1\), \(N\) is the total number of sampling times, and \(\sigma_n\) follows a prior distribution \(P({\sigma_{\text{sys}}})\).

The likelihood function for parameter inference, incorporating both statistical and systematic uncertainties, is given by

where \(m^j_i\) is the observed AB magnitude in filter \(j\) at time \(t_i\) with corresponding measurement error \(\sigma^j_i \equiv \sigma^j(t_i)\), while \(m_i^{j, \text { est }}(\vec{\theta})\) is the estimated AB magnitude for the source parameters \(\vec{\theta}\) (e.g., ejecta masses, velocities) from the model. The term \({\sigma_{\text{sys},i}}\) represents the interpolated systematic error at time \(t_i\).

The described procedure is illustrated in Fig. 6.1, where black crosses represent the mock data and red line denotes the best-fit light curve.⁴Best-fit light curves are obtained by maximizing the likelihood function, also known as Maximum Likelihood Estimation (MLE), where the MLE estimate \(\hat\theta\) given by, \(\hat{\theta}=\underset{\theta \in \Theta}{\arg \max } \mathcal{L}(\theta)\). This \(\hat\theta\) is the set of parameters that maximizes the likelihood function \(\mathcal{L}(\theta)\). The light curves are then generated using the best-fit parameters \(\hat\theta\). The red-shaded region corresponds to the systematic uncertainty in the best-fit posterior. Additionally, the surrounding Gaussian distribution represents the denominator of Eq. 6.6, encapsulating both statistical and systematic uncertainties.

6.1.3Time- and Filter-Dependent Systematic Uncertainty

At early times, when emission peaks in the ultraviolet and blue bands, significant uncertainties arise due to incomplete atomic data, relativistic expansion effects, and deviations from spherical symmetry. In contrast, late-time infrared emission is affected by non-local thermodynamic equilibrium (non-LTE) conditions, evolving opacities of lanthanides and actinides, and uncertainties in atomic structure data. The distinct sources of error across different spectral regions necessitate a systematic treatment that accounts for their wavelength dependence to improve parameter inference from kilonova observations.

A single, time-independent uncertainty model is insufficient given the spectral evolution of kilonovae, as different opacity sources and thermalization efficiencies dominate at different wavelengths and epochs. The early-time blue emission is particularly sensitive to uncertainties in the opacity of light r-process elements, while late-time infrared emission depends on the complex atomic physics of heavy elements.

To address these challenges, we extend the time-dependent interpolation scheme to incorporate filter-specific systematic errors. This is achieved by defining \({\sigma_{\text{sys}}}\) as a function of time and filter, \({\sigma_{\text{sys}}}^j(t_i)\), where \(j\) denotes the filter index. The selection of filters requiring independent and joint treatment of systematic uncertainties is determined by evaluating the Mean Absolute Deviation (MAD) for each filter and comparing it to the overall MAD computed across all filters. The MAD serves as a statistical measure of dispersion, by determining the variability within individual filters relative to the entire dataset. The filter-specific MAD is computed as

where \(x_{i,j}\) is the AB magnitude in filter \(j\) at time \(t_i\), \(\mu_{j}\) is the mean of all AB magnitudes in filter \(j\), and \(n_j\) is the total number of data points in filter \(j\). Similarly, the total MAD is calculated as

where \(x_{i}\) is the AB magnitude at time \(t_i\), \(\mu\) is the mean of all AB magnitudes and \(n\) is the total number of data points across all available filters.

Given this, the likelihood in Eq. 6.6 is modified to include filter-specific systematic uncertainties, resulting in

where \({\sigma_{\text{sys},i}}^j\) is the interpolated systematic error at time \(t_i\) and filter \(j\).

6.2Validation

To assess the validity of our approach, we performed a series of injection-recovery tests on synthetic light curves generated using the Ka2017 and Bu2019lm models. These synthetic light curves were constructed with a uniform time step of 0.5 days over a duration of 20 days. To simulate real-world observational constraints, such as weather conditions, instrumental limitations, and scheduling gaps, we randomly selected 45% of the data points for analysis.

Since the synthetic light curves generated by these models provide single values at each time step without measurement (statistical) uncertainties, we introduced a random noise to each data point. This noise follows a Gaussian distribution with zero mean and a standard deviation of 0.1 mag, thereby mimicking realistic observational noise. The perturbed light curves then serve as injections for our parameter estimation pipeline.

Given the limitations of the underlying models, we restrict the data used for parameter estimation to the first seven days for Ka2017 and the first fourteen days for Bu2019lm. In particular, we limit the Ka2017 model to one week because some light curves exhibit unphysical features at later epochs, such as an artificial increase in luminosity, likely due to the low signal-to-noise ratio in the underlying simulations at late times (cf. App. C).

We then performed cross-matching injection-recovery tests, wherein the light curves generated from Ka2017 were injected into the Bu2019lm model and vice versa. A total of four injection-recovery tests were performed, incorporating different systematic uncertainty treatments: (1) a fixed \(1~\rm mag\) systematic uncertainty, (2) a freely sampled systematic uncertainty, (3) a time-dependent systematic uncertainty, and (4) a time- and filter-dependent systematic uncertainty.

The injection parameters used in this work are motivated by the observed properties of AT2017gfo and are summarized in Tab. 6.1, along with the respective priors. To facilitate a direct comparison between the posterior distributions of the two competing models, we focus primarily on their shared parameters—namely, the luminosity distance, \(D_L\), and the total ejecta mass, \(\log_{10} M_{\rm{ej}}\). In the case of the Bu2019lm model, the total ejecta mass is defined as the sum of the wind and dynamical ejecta components, \(M_{\textrm{ej}} = M_{\mathrm{ej,wind}} + M_{\mathrm{ej,dyn}}\).

6.2.1Fixed Systematic Uncertainty

Following Refs. [19,67,151], we first consider a fixed systematic uncertainty of \(1~\rm mag\) for all data points. From Fig. 6.2, it is evident that the assumed 1 mag uncertainty is sufficient for the recovered light curves to reliably approximate the injected light curves, even when different models are used for the injection and recovery. Regarding the recovery of the injected parameters, we find that when the same model is used for both the injection and recovery, the injected parameters are recovered reliably within the statistical uncertainties, without any noticeable bias. However, when a different model is used for recovery, we observe systematic biases in the recovered posteriors of the source parameters.

In general, the Bu2019lm model recovers the injected values more accurately for both models than the Ka2017 model. This difference can be attribute to the fact, that, the two-component ejecta in the Bu2019lm model provides greater flexibility during light curve fitting.

6.2.2Freely Sampled Systematic Uncertainty

The inclusion of systematic uncertainty as a freely sampled parameter provides a more realistic representation of the uncertainties encountered in real-world kilonova observations. In this test, the systematic uncertainty is allowed to vary, drawn from a uniform prior between 0 and 2, thereby enabling the inference framework to marginalize over potential deviations between the models and the observed data. The results of this analysis are shown in Fig. 6.3. While the recovered light curves are largely consistent with the injected light curves, the recovered posteriors exhibit significant deviations, particularly in the ejecta mass parameter.

6.2.2.1Self-Consistent Injection and Recovery

When performing injection and recovery with the same model, Bu2019lm for both injection and recovery, the light curve fits remain consistent with the injected data. However, despite this agreement in photometric evolution, the posterior distributions indicate that the ejecta mass parameter is not well recovered within the 90% credible interval (CI). The systematic uncertainty posterior is notably high, indicating that the inference framework is compensating for model discrepancies through systematic uncertainty rather than constraining the physical parameters. This behavior suggests that even within a self-consistent framework, freely sampling systematic uncertainty can obscure the true ejecta properties by absorbing differences between the model and the data.

Similarly, for the Ka2017 self-consistent case, the light curve fit remains accurate, and the systematic uncertainty posterior is relatively lower. However, despite the lower inferred systematic uncertainty, the ejecta mass parameter is still not recovered within the 90% CI, suggesting that, Ka2017, despite being self-consistent in this scenario, lacks the necessary complexity to fully constrain the ejecta mass.

6.2.2.2Cross-Model Injection and Recovery

When injecting with the Bu2019lm model and recovering with Ka2017, the light curve fits remain satisfactory, but the ejecta mass is not recovered reliably. The posterior distribution of systematic uncertainty remains low, suggesting that the Ka2017 model can adequately fit the light curves without requiring significant adjustments to systematic uncertainty. However, this comes at the cost of parameter accuracy, as the ejecta mass is poorly constrained. This indicates that Ka2017, being a simpler model compared to Bu2019lm, forces a fit to the data rather than fully capturing the underlying ejecta physics.

Conversely, when injecting with Ka2017 and recovering with Bu2019lm, the inference framework is able to recover both the ejecta mass and luminosity distance parameters. However, the posteriors exhibit increased uncertainties, as reflected in the broader light curve credible regions and parameter distributions. This suggests that Bu2019lm, having a more flexible two-component ejecta prescription, can better accommodate the Ka2017 injection while maintaining some level of parameter recovery. Nonetheless, this flexibility results in wider uncertainties, indicating that while Bu2019lm avoids strong biases, it does not necessarily lead to tighter constraints on source properties as compared to Ka2017 while also trying to compensate for the systematic uncertainties.

6.2.2.3Discussion

The results indicate a fundamental trade-off between light curve fitting and parameter recovery. In cases where systematic uncertainty is inferred to be low, the light curves are well-matched to the observed data, but the parameter recovery is poor. On the other hand, when systematic uncertainty is inferred to be high, the parameter recovery improves but at the expense of larger uncertainties in both the light curves and posteriors.

Furthermore, model mismatch between injection and recovery plays a significant role in parameter biases. Ka2017 recovery struggles to capture the ejecta mass of Bu2019lm injections, whereas Bu2019lm is more flexible and can accommodate Ka2017 injections, albeit with increased uncertainties. This suggests that more complex models, while providing better overall recovery, may introduce additional uncertainty in their predictions.

A key observation is that Ka2017 recovery fails to constrain the ejecta mass, even in self-consistent scenarios. This suggests that the issue is not merely due to systematic uncertainty but is intrinsic to the Ka2017 model’s ability to constrain the ejecta mass from light curve data. On the other hand, Bu2019lm exhibits greater flexibility but at the cost of increased uncertainty.

6.2.3Time-Dependent Systematic Uncertainties

Relaxing the assumption of a constant 1 mag systematic uncertainty and instead allowing for a time-dependent uncertainty significantly impacts the accuracy of the recovered light curves. This approach acknowledges that systematic uncertainties evolve over time, better reflecting the complexity of real observational data. By permitting systematic uncertainty to vary dynamically, we aim for a more realistic marginalization over model discrepancies and data uncertainties.

6.2.3.1Self-Consistent Injection and Recovery

Performing self-consistent injection and recovery represents an idealized, best-case scenario where the inference model fully encapsulates the physics governing the injected kilonova light curve. In such a case, we expect the estimated systematic uncertainty, \(\sigma_{\rm sys}(t)\), to remain minimal throughout the entire evolution of the light curve. This expectation is confirmed in our analysis and is evident in Fig. 6.4. The 90% posterior light curve bands remain tightly bound around the injected data, with the injected light curve consistently falling within these bands across all filters and times.

The inferred systematic uncertainty remains on the order of \(\mathcal{O}(0.1 \rm mag)\), which primarily stems from the observational uncertainties artificially introduced into the injection data. This improved accuracy of the recovered light curves translates directly into more precise parameter inference. Specifically, the posterior distributions of source parameters exhibit a faithful parameter recovery and reduced spread around the injected values, with the width of the posterior shrinking by up to a factor of two. As seen in Fig. 6.4, the posterior distributions for distance and ejecta mass are more tightly constrained around the true injected values compared to the previous scenarios.

6.2.3.2Cross-Model Injection and Recovery

The impact of time-dependent systematic uncertainty becomes more complex when performing cross-model injection and recovery. When injecting Ka2017 light curves, both the Ka2017 and Bu2019lm models are able to recover the injected parameters within reasonable accuracy. However, a significant difference emerges in the uncertainty of the recovered parameters. The Bu2019lm model’s posteriors exhibit considerably broader distributions—over five times larger than those of Ka2017. This suggests that Bu2019lm has the flexibility to fit the injected light curves well, but at the cost of increased uncertainty in the recovered parameters. The systematic uncertainty posterior further supports this interpretation, as it shows that Bu2019lm compensates for discrepancies between the models by inflating the inferred uncertainty.

Conversely, the situation is different when injecting Bu2019lm light curves and attempting to recover them using the Ka2017 model. In this case, the recovered ejecta mass deviates significantly from the injected value, indicating a substantial bias in the recovered parameters. The root cause of this discrepancy lies in the fundamental differences between the models. The Ka2017 model, being a spherically symmetric, single-component ejecta model, lacks the flexibility to fully capture the complex structure of Bu2019lm light curves. As a result, it struggles to accurately recover the source parameters, leading to biased posteriors and poor parameter inference.

6.2.3.3Discussion

These finding underscores the importance of selecting physically motivated models that adequately capture the multi-component nature of kilonovae. While time-dependent systematic uncertainty certainly improves light curve fitting, it cannot fully compensate for a model that lacks the necessary physical complexity. The failure of Ka2017 to recover Bu2019lm injections exhibits the limitations of simplistic ejecta prescriptions.

Nonetheless, time-dependent systematic uncertainty provides a more flexible and realistic approach to modeling observational uncertainties in kilonova light curves. By allowing \(\sigma_{\rm sys}(t)\) to evolve dynamically, we achieve improved light curve fits and reduced parameter uncertainty, particularly when the inference model closely matches the injected model.

The key takeaway is that while systematic uncertainty marginalization can mitigate some model discrepancies, it cannot fully correct for an inadequate physical description. The failure of Ka2017 to recover Bu2019lm injections demonstrates that kilonova models must capture multi-component ejecta structures to provide reliable parameter estimates.

6.2.4Time- and Filter-Dependent Uncertainties

In our final test analysis, we introduce a time- and filter-dependent uncertainty model, allowing systematic uncertainties to vary not only with time but also across different photometric bands. This approach is motivated by initial tests, particularly those analyzing AT2017gfo, which revealed the most significant differences in the ultraviolet (\(u\)-band) and infrared (\(K\)-band) filters. To account for this, we allow separate systematic uncertainties for the \(u\)- and \(K\)- bands while grouping all other bands (\(g\) to \(H\)) under a single uncertainty parameter.

This choice, while practical, is not unique. A fully generalized approach would involve assigning independent systematic uncertainties to each filter, which could provide even finer control over model discrepancies. However, increasing the number of free uncertainty parameters also introduces computational challenges. Compared to the purely time-dependent uncertainty model discussed in the previous section, our filter-dependent analysis requires three times as many free parameters, significantly increasing the runtime of the inference process. Despite this trade-off, incorporating filter-dependent uncertainties allows us to better assess the impact of model discrepancies in different spectral regions.

6.2.4.1Self-Consistent Injection and Recovery

As expected, when using the same model for both injection and recovery, the resulting light curve reconstruction remains highly accurate. The inferred systematic uncertainties remain small across all filters, reflecting the fact that the model is trying it’s best to encapsulate the underlying physics of the injected kilonova light curve. As a result, the recovered source parameters show minimal deviation from the injected values, with uncertainties comparable to those observed in the purely time-dependent uncertainty case.

This result reinforces our previous findings: when the inference model is well-matched to the true underlying kilonova physics, marginalizing over systematic uncertainty (whether time-dependent or filter-dependent) does not introduce significant bias. Instead, the recovered light curves and posteriors exhibit only minor broadening, reflecting the observational uncertainties rather than fundamental modeling limitations. The consistency of these results suggests that incorporating filter-dependent uncertainties does not fundamentally alter the conclusions drawn from the time-dependent uncertainty approach, at least in the case of self-consistent injection and recovery.

6.2.4.2Cross-Model Injection and Recovery

The limitations of model-dependent recovery become more pronounced when applying this filter-dependent approach to cross-model injection and recovery. When attempting to recover Bu2019lm injections using the Ka2017 model, we again observe a failure to accurately infer the ejecta mass. This is consistent with our previous findings, where the Ka2017 model—being a single-component ejecta model—lacks the flexibility to fully capture the complex, multi-component structure of Bu2019lm injections.

Despite the added degrees of freedom introduced by filter-dependent uncertainties, the underlying model discrepancies persist. While the inferred light curves may achieve better fits due to the additional flexibility, the posteriors remain biased. This suggests that marginalizing over systematic uncertainties—whether through a time-dependent or filter-dependent approach—cannot fully compensate for a fundamental mismatch between the injection and recovery models.

In contrast, when injecting Ka2017 light curves and recovering them with Bu2019lm, the results mirror those of the time-dependent case. The Bu2019lm model successfully fits the injected light curves, perhaps with broader posteriors, particularly for the ejecta mass. This suggests that while Bu2019lm has the flexibility to accommodate Ka2017-like light curves, it does so at the cost of increased parameter uncertainty. The systematic uncertainty posteriors further confirm this interpretation, as they remain elevated, indicating that Bu2019lm compensates for model discrepancies by inflating the inferred uncertainty.

6.2.4.3Discussion

Our results demonstrate that incorporating filter-dependent uncertainties does not fundamentally alter the conclusions drawn from time-dependent uncertainty models. While allowing systematic uncertainties to vary across different photometric bands improves light curve fits, it does not resolve the underlying issue of model mismatch. In particular, simplistic models like Ka2017 continue to struggle when applied to complex, multi-component ejecta structures, leading to biased parameter inference.

While, the systematic uncertainty marginalization, whether time-dependent or filter-dependent, is a valuable tool for addressing observational uncertainties it still cannot fully correct for deficiencies in the physical modeling of kilonovae. Hence, even with the added flexibility of time- and filter-dependent uncertainties, the uncertainties must be limited at the physical modeling stage.

7Analyzing AT2017gfo

In the previous section, we validated our uncertainty quantification method using synthetic injections. These tests confirmed the robustness of our approach but also highlighted the importance of employing more physically complete models, such as those incorporating non-spherical symmetry when interpreting real observational data. Building on these results, we now apply our uncertainty framework to the observed kilonova AT2017gfo. Given the relatively better performance of the Bu2019lm model during our validation tests, we focus exclusively on this model for further analysis. Tab. 7.1 reports the priors used for the analysis of AT2017gfo.

7.1Time-Dependent Systematic Uncertainties

An important aspect of this investigation is to assess the impact of the number of time nodes on the inferred posteriors. To this end, we perform parameter estimation using seven different time node configurations, with \(N \in \{4,6,8,10,12,14,16\}\). For comparison, we also analyze AT2017gfo using a constant \(1~\rm mag\) systematic uncertainty, as well as a freely sampled systematic error modeled as \({\sigma_{\text{sys}}} = \mathcal{U}(0,2)\). The results of these analyses are summarized in Tab. 7.2. The best-fit light curves for different time node configurations are shown in Fig. 7.1. The recovered posteriors are visualized in Fig. B.5 and further summarized in Tab. 7.2.

7.1.1Impact of Time Nodes on Posterior Distributions

Our results indicate that the number of time nodes has a relatively minor influence on the recovered source parameters. The inferred posteriors remain consistent across different configurations, with their uncertainties overlapping significantly. This suggests that our method remains relatively stable regardless of the specific choice of time nodes, as long as a sufficient number is used to capture the variability of the kilonova light curve.

However, an important trend emerges when examining the behavior of the systematic uncertainty as a function of time. As shown in the bottom panel of Fig. 7.1, the inferred systematic uncertainty is largest at early times (\(<1\) day) and at late times (\(>5\) days). This can be attributed to two key factors:

• Sparse Observational Coverage at Early Times: The scarcity of early-time data points introduces larger uncertainties in the inferred light curves (for e.g., until 2 days there is only one data point in \(u\) band and \(\leq 4\) data points in \(zyJHK\) bands). This highlights the critical need for rapid follow-up observations of kilonovae, as emphasized in recent observational studies (e.g., [171]).

• Limitations of the Bu2019lm Model: The systematic uncertainties at early (\(\lesssim 1\) day) and late (\(\gtrsim 5\) days) times likely reflect intrinsic limitations in the underlying assumptions of the Bu2019lm model. Specifically, early-time uncertainties arise due to the opacities implemented in POSSIS, which are computed only for low-ionization states  [142,74]. At late times, systematic uncertainties increase as the assumption of local thermodynamic equilibrium starts to break down [167,166].

7.1.2Selecting an Optimal Time Node Configuration

While all time node configurations yield broadly consistent posterior distributions, an important consideration is determining which setup is sufficiently flexible to represent the underlying systematic uncertainties without introducing excessive computational overhead. As mentioned in Sec. 4.1.4, model selection is often performed using Bayes factors. However, in our case, the systematic uncertainty parameters do not represent physical quantities but rather introduce artificial degrees of freedom. To avoid invoking Occam’s razor¹Occam’s razor is a philosophical principle stating that, among competing hypotheses that explain the same phenomenon, the simplest one is preferred. It suggests that unnecessary complexity should be avoided unless additional assumptions significantly improve explanatory power. inappropriately, we instead compare the maximum log-likelihood across different configurations, as reported in Tab. 7.2.

As expected, the most flexible configuration, with 16 time nodes, achieves the highest maximum log-likelihood. However, the improvement in log-likelihood is marginal compared to configurations with fewer time nodes, while the computational cost increases substantially. Beyond the eight-time-node configuration, the runtime scales significantly, making higher node configurations computationally expensive without offering proportional gains in model accuracy. Based on this trade-off, we adopt the eight-time-node configuration as our reference for further time- and filter-dependent analysis.

7.2Time- and Filter-Dependent Systematic Uncertainties

In addition to the previous analyses, we apply our method to AT2017gfo while incorporating time- and filter-dependent uncertainties.

The choice of filters to be sampled independently is motivated by the MAD, cf. Sec. 6.1.3. From Fig. 7.3, we observe that the \(u\) and \(K\) bands exhibit the highest positive deviation from the baseline MAD. This indicates that the systematic uncertainties in these bands biases the parameter estimation from those in the other bands, justifying the need to sample them independently.

The best-fit light curves obtained from this analysis are presented in Fig. 7.4, and the recovered source parameters are summarized in Tab. 7.2.²See App.  for full posteriors and light curves. From Fig. 7.4, it is evident that the systematic uncertainty in the \(u\) band remains relatively low, below 1 magnitude for the first six days, despite the presence of only a single detection point. However, after six days, the uncertainty increases as the model struggles to provide an accurate fit due to the sparsity of observational data. For the \(K\) band, while for first six days, the error is below 0.5 magnitudes, the uncertainty increases to as high as 2 magnitudes at later times. This is likely due to the fact that for \(K\) model is likely to breakdown due to reason discussed previously. The other bands, such as \(z\), \(y\), and \(J\), exhibit a more stable behavior, with uncertainties remaining below 1 magnitude throughout the observed time range.

Although, introducing the filter-dependent systematic uncertainty leads to a substantial increase and decrease in the maximum log-likelihood compared to the 8 time nodes configuration of filter-independent uncertainty approach, the recovered parameters themselves remain largely unchanged (cf. Tab. B.1). This suggests that model uncertainties are absorbed by the systematic error term, there by inflating the uncertainties in the light curves.

8Conclusion

In this thesis, we have developed and validated a comprehensive framework for quantifying systematic uncertainties in kilonova light curve modeling, with particular emphasis on their time and filter dependencies. Through a systematic investigation of different uncertainty parameterizations, we have demonstrated that properly accounting for these uncertainties is crucial for accurate parameter inference in kilonova observations.

Our validation tests using synthetic light curves from the Ka2017 and Bu2019lm models revealed several important findings. First, when using a fixed systematic uncertainty approach (1 mag), we observed that parameters could be recovered accurately within statistical uncertainties when the same model was used for both injection and recovery. However, significant biases emerged during cross-model analyses, highlighting the inherent limitations of simplistic models and uncertainty treatments.

The freely sampled systematic uncertainty approach, while providing more flexibility, illustrated a fundamental trade-off between light curve fitting and parameter recovery. Low inferred systematic uncertainties often corresponded to well-fitted light curves but poor parameter recovery, while high inferred uncertainties improved parameter recovery at the cost of larger uncertainties in light curves and posteriors.

A key insight from our work is that time-dependent systematic uncertainties significantly improved parameter recovery across all scenarios. By allowing uncertainties to evolve dynamically over time, we achieved more accurate light curve fits and reduced parameter uncertainty, particularly when the inference model closely matched the injected model. This approach proved especially valuable in self-consistent scenarios, where the estimated systematic uncertainty remained minimal throughout the light curve evolution.

Our most sophisticated approach, involving both time and filter dependence in systematic uncertainties, provided additional flexibility by allowing uncertainties to vary not only with time but also across different photometric bands. While this approach improved light curve fits, it did not fundamentally alter our conclusions regarding model mismatch. In particular, simplistic models like Ka2017 continued to struggle when applied to complex, multi-component ejecta structures, leading to biased parameter inference regardless of uncertainty treatment.

Application of our methodology to AT2017gfo revealed that systematic uncertainties are most prominent at early times (<1 day) and late times (>5 days). This pattern reflects both observational limitations (sparse early-time coverage) and inherent model limitations (opacity approximations at early times and breakdown of local thermodynamic equilibrium assumptions at late times). Our analysis of AT2017gfo using different time node configurations showed that while all configurations yielded broadly consistent posterior distributions, an eight-time-node configuration offered the optimal balance between flexibility and computational efficiency.

Perhaps most significantly, our work demonstrates that while systematic uncertainty marginalization can mitigate some model discrepancies, it cannot fully compensate for fundamental limitations in physical modeling. The persistent failure of Ka2017 to recover Bu2019lm injections underscores the critical importance of employing physically complete models that accurately capture the multi-component nature of kilonovae.

The framework developed in this thesis represents a significant advancement in uncertainty quantification for kilonova modeling and can be readily extended to other electromagnetic transient phenomena. By providing a more nuanced approach to systematic error treatment, our work contributes to improving the reliability of parameter estimation in future kilonova observations, which will be crucial for extracting meaningful astrophysical insights from upcoming gravitational wave events with electromagnetic counterparts.

As we enter an era of multi-messenger astronomy with increasingly sensitive gravitational wave detectors and rapid electromagnetic follow-up capabilities, the accurate characterization of kilonova properties will play an essential role in constraining neutron star equation of state, heavy element nucleosynthesis, and cosmological parameters. The methodologies developed in this thesis provide a robust foundation for these future endeavors by ensuring that parameter uncertainties accurately reflect both statistical and systematic sources of error.

AppendixAPOSSIS Systematics

As highlighted in Ch. 2 and Sec. 3.1, the luminosity of a kilonova is affected the opacity of the ejecta material, governed by how easily photons can escape from the expanding ejecta, thereby shaping both the peak brightness and the temporal evolution of the light curve.

To systematically investigate the role of opacity in shaping the observable light curves, we perform POSSIS simulations by scaling the opacity \(\kappa\) by a scaling factor \(\alpha \in \{0.5,~1,~2\}\) and compute the synthetic light curves. These scalings do not correspond to physically different compositions but rather mimic the uncertainty or variability in opacities due to, for example, incomplete atomic data or variations in ejecta ionization states.

Fig. A.1 shows the resultant broadband light curves for 11 viewing angles (uniformly distributed in cosine of observation angle), across three values of opacity scaling. With \(\alpha = 1\), the fiducial opacity \(\kappa\) is used without modification. The case of \(\alpha = 0.5\) corresponds to a reduced opacity, allowing photons to diffuse more readily through the ejecta, resulting in a brighter and earlier-peaking kilonova. In contrast, increasing the opacity by a factor of two (\(\alpha = 2\)) leads to a delayed and fainter emission profile due to longer photon diffusion timescales and increased trapping of radiation within the ejecta.

An additional observation from the opacity-scaling is that light curves computed with higher opacity values (i.e., \(\alpha = 2\)) exhibit more erratic and noisy behavior. This effect is most noticeable in the \(\theta = 90^\circ\) (equatorial) viewing direction, where the line-of-sight optical depth is already higher due to the geometry of the ejecta. This variability can be attributed to the fact that, with higher opacity, photons undergo more scatterings and absorptions before escaping the ejecta. This results in a more complex and stochastic propagation history for each photon packet. Even small differences in the location of emission or path taken can significantly affect the final observed flux, leading to increased variability in the light curve, especially when sampled over discrete time bins.

Moreover, at high opacity, photon packets can become “trapped” for longer periods, undergoing a large number of interactions. This increases the computational burden and decreases the statistical efficiency of the sampling. Unless the number of photon packets is dramatically increased (which becomes computationally expensive), this leads to larger shot noise and visible fluctuations in the synthetic light curves.

Relevance to Systematic Uncertainties in Kilonova Modeling

This opacity-scaling analysis provides a physically motivated illustration of how microphysical modeling uncertainties propagate into observable features of kilonova light curves. Variations in opacity encapsulate systematic errors arising from uncertain atomic physics, ionization conditions, and thermalization processes—factors that are particularly impactful at early and late times. Because opacity¹We encourage the readers to refer Refs. [142,74,76,172] for a detailed discussion on the opacities. directly influences photon diffusion, even small deviations can lead to observable differences in brightness and temporal structure, introducing degeneracies between ejecta properties and radiative transfer assumptions.

The increased irregularity observed at higher opacities further highlights the interplay between physical uncertainty and numerical limitations. These effects underscore the necessity of modeling frameworks that do not assume fixed or static uncertainties but instead allow for time- and filter-dependent treatments, as developed in this thesis.

AppendixBAT2017gfo Plots

Constant Systematic Uncertainty

Freely Sampled Systematic Uncertainty

Time- and Filter- Dependent Systematic Uncertainties

.

AppendixCSEDONA Light Curves

Here, we present eleven light curves from the SEDONA model grid, each exhibiting a characteristic drop in magnitude at approximately 8-10 days (gray shaded region). Given this behavior, we restricted the use of SEDONA light curves to the first seven days in the main text to ensure a more consistent and reliable analysis. By limiting the time range, we focus on the phase where the model remains well-constrained and avoids potential biases introduced by uncertainties in late-time evolution.

AppendixDCode and Documentation

D.1Code

D.1.1Parsing systematics from .yaml file

import inspect
import warnings
from pathlib import Path

import yaml
from bilby.core.prior import analytical

warnings.simplefilter("module", DeprecationWarning)


class ValidationError(ValueError):
    def __init__(self, key, message):
        super().__init__(f"Validation error for '{key}': {message}")


ALLOWED_FILTERS = [
    "2massh",
    "2massj",
    "2massks",
    "atlasc",
    "atlaso",
    "bessellb",
    "besselli",
    "bessellr",
    "bessellux",
    "bessellv",
    "ps1__g",
    "ps1__i",
    "ps1__r",
    "ps1__y",
    "ps1__z",
    "sdssu",
    "uvot__b",
    "uvot__u",
    "uvot__uvm2",
    "uvot__uvw1",
    "uvot__uvw2",
    "uvot__v",
    "uvot__white",
    "ztfg",
    "ztfi",
    "ztfr",
]

ALLOWED_DISTRIBUTIONS = dict(inspect.getmembers(analytical, inspect.isclass))


def get_positional_args(cls):
    init_method = cls.__init__

    signature = inspect.signature(init_method)
    params = [
        param.name
        for param in signature.parameters.values()
        if param.name != "self" and param.default == inspect.Parameter.empty
    ]

    return params


DISTRIBUTION_PARAMETERS = {k: get_positional_args(v) for k, v in ALLOWED_DISTRIBUTIONS.items()}


def load_yaml(file_path):
    return yaml.safe_load(Path(file_path).read_text())


def validate_only_one_true(yaml_dict):
    for key, values in yaml_dict["config"].items():
        if "value" not in values or not isinstance(values["value"], bool):
            raise ValidationError(key, "'value' key must be present and be a boolean")
    true_count = sum(value["value"] for value in yaml_dict["config"].values())
    if true_count > 1:
        raise ValidationError("config", "Only one configuration key can be set to True at a time")
    elif true_count == 0:
        raise ValidationError("config", "At least one configuration key must be set to True")


def validate_filters(filter_groups):
    used_filters = set()
    for filter_group in filter_groups:
        if isinstance(filter_group, list):
            filters_in_group = set()  # Keep track of filters within this group
            for filt in filter_group:
                if filt not in ALLOWED_FILTERS:
                    raise ValidationError(
                        "filters",
                        f"Invalid filter value '{filt}'. Allowed values are {', '.join([str(f) for f in ALLOWED_FILTERS])}",
                    )
                if filt in filters_in_group:
                    raise ValidationError(
                        "filters",
                        f"Duplicate filter value '{filt}' within the same group.",
                    )
                if filt in used_filters:
                    raise ValidationError(
                        "filters",
                        f"Duplicate filter value '{filt}'. A filter can only be used in one group.",
                    )
                used_filters.add(filt)
                filters_in_group.add(filt)
                # Add the filter to the set of used filters within this group
        elif filter_group is not None and filter_group not in ALLOWED_FILTERS:
            raise ValidationError(
                "filters",
                f"Invalid filter value '{filter_group}'. Allowed values are {', '.join([str(f) for f in ALLOWED_FILTERS])}",
            )
        elif filter_group in used_filters:
            raise ValidationError(
                "filters",
                f"Duplicate filter value '{filter_group}'. A filter can only be used in one group.",
            )
        else:
            used_filters.add(filter_group)


def validate_distribution(distribution):
    dist_type = distribution.get("type")
    if dist_type not in ALLOWED_DISTRIBUTIONS:
        raise ValidationError(
            "distribution type",
            f"Invalid distribution '{dist_type}'. Allowed values are {', '.join([str(f) for f in ALLOWED_DISTRIBUTIONS])}",
        )

    required_params = DISTRIBUTION_PARAMETERS[dist_type]

    missing_params = set(required_params) - set(distribution.keys())
    if missing_params:
        raise ValidationError(
            "distribution", f"Missing required parameters for {dist_type} distribution: {', '.join(missing_params)}"
        )


def create_prior_string(name, distribution):
    dist_type = distribution.pop("type")
    _ = distribution.pop("value", None)
    _ = distribution.pop("time_nodes", None)
    _ = distribution.pop("filters", None)
    prior_class = ALLOWED_DISTRIBUTIONS[dist_type]
    required_params = DISTRIBUTION_PARAMETERS[dist_type]
    params = distribution.copy()

    extra_params = set(params.keys()) - set(required_params)
    if extra_params:
        warnings.warn(f"Distribution parameters {extra_params} are not used by {dist_type} distribution and will be ignored")

    params = {k: params[k] for k in required_params if k in params}

    return f"{name} = {repr(prior_class(**params, name=name))}"


def handle_withTime(values):
    validate_distribution(values)
    filter_groups = values.get("filters", [])
    validate_filters(filter_groups)
    result = []
    time_nodes = values["time_nodes"]

    for filter_group in filter_groups:
        if isinstance(filter_group, list):
            filter_name = "___".join(filter_group)
        else:
            filter_name = filter_group if filter_group is not None else "all"

        for n in range(1, time_nodes + 1):
            prior_name = f"sys_err_{filter_name}{n}"
            result.append(create_prior_string(prior_name, values.copy()))

    return result


def handle_withoutTime(values):
    validate_distribution(values)
    return [create_prior_string("sys_err", values)]


config_handlers = {
    "withTime": handle_withTime,
    "withoutTime": handle_withoutTime,
}


def main(yaml_file_path):
    yaml_dict = load_yaml(yaml_file_path)
    validate_only_one_true(yaml_dict)
    results = []
    for key, values in yaml_dict["config"].items():
        if values["value"] and key in config_handlers:
            results.extend(config_handlers[key](values))
    return results

D.1.2Injecting systematic configuration into PriorDict

from bilby.core.prior import PriorDict as _PriorDict
from . import systematics

def from_list(self, systematics):
    """
    Similar to `from_file` but instead of file buffer, takes a list of Prior strings
    See `from_file` for more details
    """

    comments = ["#", "\n"]
    prior = dict()
    for line in systematics:
        if line[0] in comments:
            continue
        line.replace(" ", "")
        elements = line.split("=")
        key = elements[0].replace(" ", "")
        val = "=".join(elements[1:]).strip()
        prior[key] = val
    self.from_dictionary(prior)

setattr(_PriorDict, "from_list", from_list)

class ConditionalGaussianIotaGivenThetaCore
                (ConditionalTruncatedGaussian):
    ...
    ...
    ...
def create_prior_from_args(model_names, args):
    AnBa2022_intersect = list(set(['AnBa2022_linear', 'AnBa2022_log']) & set(model_names))

    if len(AnBa2022_intersect) > 0:
        ...
        ...
        ...
    else:
        priors = bilby.gw.prior.PriorDict(args.prior)

    if args.systematics_file is not None:
            systematics_priors = systematics.main(args.systematics_file)
            priors.from_list(systematics_priors)
    return priors

D.1.3Likelihood Calculation

from __future__ import division

import numpy as np
import scipy.stats
from scipy.interpolate import interp1d
from scipy.stats import truncnorm
from pathlib import Path
import yaml

from . import utils, systematics

class OpticalLightCurve(Likelihood):
    ...
    ...
    ...
        self.systematics_file = systematics_file
    ...
    ...
    ...
    def log_likelihood(self):
        lbol, mag_abs = self.light_curve_model.generate_lightcurve(
            self.sample_times, self.parameters
        )
        ...
        ...
        ...
        if self.systematics_file is not None:
            yaml_dict = yaml.safe_load(
                Path(self.systematics_file).read_text())
            systematics.validate_only_one_true(yaml_dict)

            if yaml_dict["config"]["withTime"]["value"]:
                n = yaml_dict["config"]["withTime"]
                ["time_nodes"]
                time_array = np.round(np.linspace(self.tmin, self.tmax, n), 2)
                yaml_filters = yaml_dict["config"]["withTime"]
                ["filters"]
                systematics.validate_filters(
                    yaml_filters)

                for filter_group in yaml_filters:
                    if isinstance(filter_group, list):
                        filt = "___".join(filter_group)
                    elif filter_group is None:
                        filt = "all"
                    else:
                        filt = filter_group

                    globals()[f"sys_err_{filt}_array"] = np.array([])

                    for i in range(1, n + 1):
                        value = self.parameters.get(f"sys_err_
                        {filt}{i}")
                        globals()[f"sys_err_{filt}_array"] = np.append(globals()[f"sys_err_
                        {filt}_array"], value)

                for filter_group in yaml_dict["config"]["withTime"]
                ["filters"]:
                    if isinstance(filter_group, list):
                        filt = "___".join(filter_group)
                    elif filter_group is None:
                        filt = "all"
                    else:
                        filt = filter_group
                    globals()[f"sys_err_{filt}_interped"] = interp1d(time_array, globals()[f"sys_err_{filt}_array"], fill_value="extrapolate", bounds_error=False)
        

        # compare the estimated light curve and the measured data
        minus_chisquare_total = 0.0
        gaussprob_total = 0.0
        for filt in self.filters:
            # decompose the data
            data_time = self.light_curve_data[filt][:, 0]
            data_mag = self.light_curve_data[filt][:, 1]
            data_sigma = self.light_curve_data[filt][:, 2]

            if self.systematics_file is not None:
                if yaml_dict["config"]["withTime"]
                ["value"]:
                    yaml_filters = yaml_dict["config"]["withTime"]
                    ["filters"]

                    filter_group_finite_idx_match = False

                    if yaml_filters is not None and yaml_filters != []:
                        for yaml_filt in yaml_filters:
                            if yaml_filt is not None and filt in yaml_filt:
                                if isinstance(yaml_filt, str):
                                    filters_to_use = yaml_filt
                                elif isinstance(yaml_filt, list):
                                    filters_to_use = "___".join(
                                        yaml_filt)
                                filter_group_finite_idx
                                _match = True
                                break
                    if not filter_group_finite_idx_match:
                        filters_to_use = "all"
                    data_sigma = np.sqrt(data_sigma**2 + (globals()[f"sys_err_{filters_
                    to_use}_interped"](data_time)) ** 2)
                    
                elif yaml_dict["config"]["withoutTime"]
                ["value"]:
                    data_sigma = np.sqrt(data_sigma**2 + self.parameters["sys_err"] ** 2)

            # include the error budget into calculation
            elif 'sys_err' in self.parameters:
                data_sigma = np.sqrt(data_sigma**2 + self.parameters['sys_err']**2)
            else:
                data_sigma = np.sqrt(data_sigma**2 + self.error_budget[filt]**2)

            # evaluate the light curve magnitude at the data points
            mag_est = mag_app_interp[filt](data_time)

            # seperate the data into bounds (inf err) and actual measurement
            infIdx = np.where(~np.isfinite(data_sigma))[0]
            finiteIdx = np.where(np.isfinite(data_sigma))[0]
            ...
            ...
            ...
            if len(infIdx) > 0:
                if self.systematics_file is not None:
                    if yaml_dict["config"]["withTime"]
                    ["value"]:
                        yaml_filters = yaml_dict["config"]["withTime"]
                        ["filters"]

                        filter_group_infinite_idx_match = False

                        if yaml_filters is not None and yaml_filters != []:
                            for yaml_filt in yaml_filters:
                                if yaml_filt is not None and filt in yaml_filt:
                                    if isinstance(yaml_
                                    filt, str):
                                        upperlim_sigma = globals()[f"sys
                                        _err_{yaml_filt}_
                                        interped"]
                                        (data_time)
                                        [infIdx]
                                    elif isinstance(yaml_
                                    filt, list):
                                        filters_to_use = "___".join(yaml
                                        _filt)
                                        upperlim_sigma = globals()[f"sys
                                        _err_{filters_to_
                                        use}_interped"]
                                        (data_time)
                                        [infIdx]
                                    filter_group_infinite
                                    _idx_match = True
                                    break
                        if not filter_group_infinite_idx
                        _match:
                            filters_to_use = "all"
                            upperlim_sigma = globals()[f"sys_err_
                            {filters_to_use}_interped"]
                            (data_time)[infIdx]

                        gausslogsf = scipy.stats.norm.logsf(
                            data_mag[infIdx], mag_est[infIdx], upperlim_sigma)

                    elif yaml_dict["config"]["withoutTime"]
                    ["value"]:
                        upperlim_sigma = self.parameters["sys_err"]
                        gausslogsf = scipy.stats.norm.logsf(
                            data_mag[infIdx], mag_est[infIdx], upperlim_sigma)
                elif 'sys_err' in self.parameters:
                    upperlim_sigma = self.parameters['sys_err']
                    gausslogsf = scipy.stats.norm.logsf(data_mag
                    [infIdx], mag_est[infIdx], upperlim_sigma)
                else:
                    ...
                    ...
                    ...

        return log_prob

D.2Documentation

NMMA currently uses –error-budget to specify the constant systematic uncertainties to be added to the likelihood quadrature.

However, it is now possible to use systematic error (\(\sigma_{sys}\)) prior in form of a freely sampled parameter, time dependent and/or filter dependent systematic error. This can done by specifying the file path using the –systematics-file flag in lightcurve-analysis command.

For more information on systematics error, refer to the main text or Ref. [1].

D.2.1Freely Sampled Systematic Uncertainty

In this case the systematic error is freely sampled and is not dependent on time or filter.

config:
    withTime:
        value: false
        filters:
            - null
        time_nodes: 4
        type: Uniform
        minimum: 0
        maximum: 2
    withoutTime:
        value: true
        type: Uniform
        minimum: 0
        maximum: 2

D.2.2Time Dependent Systematic Uncertainty

In this configuration, a time dependent single systematic error is applied across all filters.

config:
    withTime:
        value: true
        filters:
            - null
        time_nodes: 4
        type: Uniform
        minimum: 0
        maximum: 2
    withoutTime:
        value: false
        type: Uniform
        minimum: 0
        maximum: 2

D.2.3Time and Filter Dependent Systematic Uncertainty

In this configuration, the sdssu and ztfr filters are sampled together for systematic errors, while the 2massks filter is sampled independently. All other filters are grouped and sampled together.

config:
    withTime:
        value: true
        filters:
            - [sdssu, ztfr]
            - [2massks]
            - null
        time_nodes: 4
        type: Uniform
        minimum: 0
        maximum: 2
    withoutTime:
        value: false
        type: Uniform
        minimum: 0
        maximum: 2

D.2.4Distribution types

Distribution can be of any of the analytical prior from bilby. Refer to bilby documentation for more information on the available distribution type and their usage. Only positional arguments are required for any of the distrbutions.