On August 18, 2017, prompt observations identified gravitational wave emission (GW170817; [1,2]), shortly followed by a gamma-ray burst (short-GRB GRB170817A [3]).
Extensive followup observations identified a long-duration optical/near-infrared counterpart, AT2017gfo, later identified as a "kilonova" [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21].
A kilonova [22] is characterized by thermal emission from rapidly expanding, radioactively-powered, heavyelement material ejected from the associated progenitor merger. The detection of the joint gravitationaland electromagnetic-wave emission from GW170817 and AT2017gfo has initiated an era of precision kilonova observations.
Most interpretations of kilonova observations have relied on broadband photometry, in part owing to the relative sparsity of available spectra for AT2017gfo (and lack of spectral observations of other kilonovae; although see e.g. [23][24][25]) [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43]. A comprehensive review of kilonova broadband photometry has recently been compiled and presented in Ref. [44]. Many studies interpreted their observations of AT2017gfo shortly after detection principally by comparison to simplified models for kilonovae [4-10, 12, 14-21, 45, 46] consisting of one or more groups of non-accelerating (homologous) expanding material. Motivated both by binary merger simulations and the inability to fit observations with one component [8], at least two components are customarily employed [10,12,[15][16][17]20], with properties loosely associated with two expected features of merger simulations: promptly ejected material (the "dynamical" ejecta), associated with tidal tails or shocked material at contact; and material driven out on longer timescales by properties of the remnant system (the "wind" ejecta) [47].
Radiative transfer models of two-component kilonovae, whilst more physically accurate and informative, come with a significantly higher computational cost compared to semi-analytical or one-component models. As a result of this cost, many groups have resorted to surrogate models, or emulators, for the kilonova outflow, to reduce the computational impact associated with inference using these more complex models [48][49][50][51][52][53][54][55].
In this paper, we use a neural network emulator for light-curve interpolation and parameter inference. We train on a previously-generated library of ∼ 400 twocomponent kilonova light-curve simulations. Our method can be easily applied to any modestly-sized archive of adaptively-learned astrophysical transient light-curve simulations.
This paper is organized as follows. In Sec. II we discuss our simulation training library, the interpolation model and its architecture, and the associated light-curve interpolation methodology. In Sec. III we compare our emulator's performance with others employed in the literature and report inference results for observations of AT2017gfo. In Sec. IV, we summarize our findings. Our two-component kilonova model consists of a lanthanide-rich equatorial dynamical ejecta component and a lanthanide-poor axial wind ejecta component as described in [56,57] and motivated by numerical simulations [47,58]. Each component is homologously expanding and parameterized by a mass and velocity such that M d , v d and M w , v w describe the dynamical and wind components' masses and averaged velocities, respectively. The morphology for the dynamical component is an equatorially-centered torus, whereas the wind component is represented by an axially-centered peanut component; Figure 1 of [56] displays the torus-peanut, or "TP," schematic corresponding to the morphologies employed in this work [see 57, for detailed definition]. The lanthanide-rich dynamical ejecta is a result of the r-process nucleosynthesis from a neutron-rich material with a low electron fraction (Y e ≡ n p /(n p + n n )) of Y e = 0.04 with elements reaching the third r-process peak (A ∼ 195), while the wind ejecta originates from higher Y e = 0.27 which encapsulates elements between the first (A ∼ 80) and second (A ∼ 130) r-process peaks. The detailed breakdown of the elements in each component can be found in Table 2 of Ref. [56].
We use SuperNu [59], a Monte Carlo code for simulation of time-dependent radiation transport with matter in local thermodynamic equilibrium, to create simulated kilonova spectra F λ,sim assuming the aforementioned two-component model. Both components are assumed to have fixed composition and morphology for the duration of each simulation. SuperNu uses radioactive power sources calculated from decaying the r-process composition from the WinNet nuclear reaction network [60][61][62][63]. These radioactive heating contributions are also weighted by thermalization efficiencies introduced in Ref. [64] (see Ref. [65] for a detailed description of the adopted nuclear heating). We use detailed opacity calculations via the tabulated, binned opacities generated with the Los Alamos suite of atomic physics codes [66][67][68]. In the database that we use, the tabulated, binned opacities are not calculated for all elements; therefore, we produce opacities for representative proxy elements by combining pure-element opacities of nuclei with similar atomic properties [67]. Specifics of the representative elements for our composition are given in Ref. [56].
The SuperNu outputs are observing-angle-dependent, simulated spectra F λ,sim , post-processed to a source distance of 10 pc, in units of erg s -1 cm -2 Å-1 . The spectra are binned into 1024 equally log-spaced wavelength bins spanning 0.1 ≤ λ ≤ 12.8 microns.
For the purposes of this work, we consider the light curves for the 2MASS grizy and Rubin Observatory JHK broadband filters. As we only consider anisotropic simulations in this study, unless otherwise noted, we extract simulated light curves using 54 angular bins, uni-formly spaced in cos θ over the range -1 ≤ cos θ ≤ 1, where the angle θ is taken between the line of sight and the symmetry axis as defined in Equation 2. Specifically, simulations in our database cover all observing angles with a resolution ranging from ∆θ = arcsin (2/54) ≃ 2.1 • at the equator, to ∆θ = arccos (1 -2/54) ≃ 15.6 • near the axes. This limiting angular resolution near the axes is comparable to the angles inferred from long-term radio observations of the gamma ray burst jet afterglow [10,16,17,[69][70][71].
SuperNu Monte Carlo radiative transport results have modest but nonzero Monte Carlo error. The impact of this Monte Carlo error can be estimated both by resolution studies as well as by simple smoothness diagnostics (e.g., versus angle and time). For example, in a previous study [50] we performed Gaussian process interpolation over this same training set, inferring both an estimate of the light curve and a conservative estimate of its variance. The Monte Carlo error inherent in the underlying SuperNu simulations can be seen for example in Figure 6 of Ref. [50] as short-angular-scale roughness on top of the overall smooth trend, with at most ∼ 3% deviation at the 1σ level. Later, in Section III C, we describe a targeted resolution study. Based on these investigations, we expect that the statistical error of the underlying SuperNu simulations is smaller than the systematic errors introduced by interpolating between these simulations as described below.
Below we describe the approach taken to generate the simulation library in Ref. [50], hereafter R22. Our training library of 22248 kilonova light-curve simulations was constructed using iterative simulation placement guided by Gaussian process variance minimization. New simulations were placed with parameter combinations that were identified as having the largest bolometric luminosity variance by our Gaussian process regression approach. In other words, we placed new simulations in regions of parameter space where our bolometric luminosity interpolation root-mean-square uncertainty was largest. The Gaussian process variance s(⃗ x) 2 is defined as
where ⃗ x is the vector of input parameters, ⃗ x a is the training data vector, the function k(⃗ x, ⃗ x ′ ) is the kernel of the Gaussian process, and the indices a, a ′ are used to calculate the covariance between inputs ⃗ x and training data ⃗ x a , ⃗ x a ′ such that if a = a ′ , the variance is 0. In building the simulation library, we only considered the four-dimensional space of ejecta parameters
Each ejecta parameter combination yields simulations calculated for 54 equally-spaced viewing angles; as such, our training set of 22248 light curves corresponds to a core set of 412 unique ejecta parameter combinations.
For this work, we use the aforementioned light curves in the original simulation library as our training set. The light curves used in this work have the same parameters as those used for our light-curve interpolation approach in R22. No additional simulations were produced for the purposes of this work; all training data came from the simulation library presented in R22.
The original training data library consists of 22248 total light-curve simulations calculated at 264 times and 54 angular bins each. We do not perform any coordinate transformations, but rather interpolate directly in our ejecta parameter space and angle.
The entirety of our 22248 simulations is represented by a tensor, M abc , containing the AB magnitudes across the bands described in Section II A, which correspond to a set of input parameters ⃗ x. M abc has dimensions of 412×264×54, corresponding to 412 simulations evaluated at 264 log-spaced times between 0.125 and 37.24 days for 54 viewing angles equally spaced in cos θ for θ ranging from 0 to 180 • . We do not perform any normalization of our inputs or outputs, with ejecta parameters ranging from -3 ≤ log m/M ⊙ ≤ -1 and 0.05 ≤ v/c ≤ 0.3 and light-curves ranging from -18 to 8 AB magnitudes.
We split our four-dimensional ejecta parameter vector ⃗ x into training, validation, and test sets. We use 60% of the data for the training set, which contains information that the neural network uses to learn. Of the remaining 40%, 20% is used for the validation check, which tracks how well the network generalizes to off-sample inputs during training, and 20% goes into the test set, which is used to evaluate the network's predictions compared to known simulation data. Neither the validation set nor the test set data is used by the neural network for learning; therefore, we only use ∼ 247 simulations for training, while the rest are used in various steps of verifying generalization (i.e. avoiding overfitting to the training data).
After splitting the data into training, validation, and test sets, we incorporate viewing angle as a fifth input parameter. As mentioned above, the viewing angles in our simulations are equally-spaced in cos θ space across 54 angular bins, as presented in the following equation:
Temporarily ignoring training, validation, or test sets, our training library consists of a total of 22248×5 inputs, as illustrated in the following schematic matrix:
1 θ 3 . . . . . . . . . . . . . . . m d,1 v d,1 m w,1 v w,1 θ 54 m d,2 v d,2 m w,2 v w,2 θ 1 m d,2 v d,2 m w,2 v w,2 θ 2 . . . . . . . . . . . . . . . m d,412 v d,412 m w,412 v w,412 θ 54 D. Neural Network Architecture and Training
We use a standard feed-forward neural network called a multi-layer perceptron (MLP). Figure 1 shows our MLP architecture. The input data of dimension 5 (blue block) is propagated through the hidden layers of the MLP (pink blocks). These hidden layers apply a sequence of linear and non-linear transformations (black arrows) to progressively map the input to a higher-dimensional space. The network has six fully-connected layers (pink blocks) of dimension 128, 256, 64, 256, 128, and 264, respectively, which are followed by Rectified Linear Unit (ReLU) activation functions, except for the middle two layers. The final layer reduces the dimension to a 264 × 1 vector, which matches the length of our light curves.
For each observing band, we train a separate neural network and compare its predictions with the simulation results. During training, the neural network's predictions for the inputs in the validation set are compared to the simulation data for those same inputs. The residual, or difference, between the two is evaluated by the mean squared error (MSE) loss function, which we use to measure the average squared difference between simulation data and the neural network prediction. We calculate the MSE according to
where y i is the corresponding simulation data at time i and ŷi is the predicted value from the MLP model at the same time. We show the evolution of the training and validation losses in Figure 2 for the g-band network. We train a separate neural network for each of the broadband filters described in Section II A. Each network is trained for 1000 epochs, as this is enough time for the validation loss to convincingly stabilize at a floor value without beginning to increase, indicating overfitting. FIG. 1. A visual representation of the neural network architecture. The blue block represents our five-dimensional inputs ⃗ x θ . The pink blocks represent hidden layers, with labels under each block representing the input and output dimensions of each. Unlabeled right arrows indicate a linear mapping between layers, while those labeled "ReLU" have a Rectified Linear Unit activation function applied to their outputs. The orange block represents the neural network prediction in the form of a 264 × 1 vector matching the length of our broadband light curves.
in Figure 3 shows no indication of overfitting to the training data. For inference applications, our neural networks can generate the outcomes corresponding to five million five-dimensional inputs in about one minute. Training each neural network takes 20 minutes on a 2022-edition Macbook Pro with an M2 chip using the CPU. Since all the bands are independent and can be trained simultaneously in parallel, training all the emulators can be completed in this same 20 minute interval. Our training time is half the value reported in Ref. [49], though it is unclear whether their reported time assumes training in parallel or in serial. Training time of networks in parallel is expected to decrease appreciably on a highperformance computing cluster.
The top left plot in Figure 3 shows a histogram of the MSE values when evaluating the simulation library parameters using the neural network. We evaluate only the 412 unique ejecta parameter combinations, fixing the viewing angle to θ = 0 • in each case. In containing simulations from the training, validation, and test sets, this histogram represents the neural network's on-and offsample fidelity. The light-curve plots in Figure 3 We note that, as with the emulators presented in Ref. [50], predictions for inputs with low-mass (log(M ) ∼ -3) components or viewing angles θ ∼ 90 • may deviate substantially from expectations. In Monte Carlo radiative transfer simulations of a kilonova, the representa- 5 1 2 4 8 16 32 64 t (days) 18 16 14 12 10 8 6 AB Mag g y r J i H z K 0.125 0.5 1 2 4 8 16 32 64 t (days) 16 14 12 10 8 6 4 AB Mag g y r J i H z K 0.125 0.5 1 2 4 8 16 32 64 t (days) 16 14 12 10 8 6 4 2 0 AB Mag g y r J i H z K FIG. 3. Top left: A histogram of MSE values, averaged across all bands by the number of observations, which characterize the deviation of the MLP's predictions from the true simulation library light curves. For simplicity, we assume a fixed viewing angle of θ = 0 • for each light curve and evaluate the MSE only for the 412 unique ejecta parameter combinations. Top right: True (points) and predicted (lines) light curves for a randomly drawn set of parameters with M SE < 0.01. Bottom left: Same as top right, but for 0.01 ≤ M SE ≤ 0.1. Bottom right : Same as top right, but for M SE > 0.1.
tion, a finite number of particles are employed to represent photons escaping from the system, forming "packets" of energy. The system's luminosity is intrinsically linked with, among other physical parameters, the ejecta mass. Consequently, when dealing with extremely lowmass components or viewing angles that look into the high-opacity dynamical ejecta, the simulations become particularly sensitive to Poisson noise due to reduced photon count. This effect, particularly with respect to low-mass ejecta, likely arises from our SuperNu simulations preferably sampling photon packets from higherenergy regions of the ejecta. In future studies, we hope to enhance the simulation interpretation under these conditions by way of an increase in photon packet count.
As in R22, we infer the parameters of the kilonova AT2017gfo using our interpolated light curves and the AT2017gfo photometric data. The AT2017gfo data is originally presented in [4-18, 20, 21]. We use the RIFT framework [72] to adaptively perform the Monte Carlo integral and generate samples using a reduced χ 2 statistic. The parameter priors are the same as in R22, with uniform ejecta parameter priors of -3 ≤ log m/M ⊙ ≤ -1 and 0.05 ≤ v/c ≤ 0.3 and a Gaussian angle prior with µ = 20 and σ = 5 degrees. Unlike before, in this work we employ an adaptive volume Monte Carlo integrator, following closely the approach outlined in Ref. [73]. The adaptive volume integrator allows for more efficient sam-pling given the higher-dimensional space being explored in this work.
Each sample ⃗ x θ is evaluated by the MLP to produce a light-curve prediction ŷ for every one of the grizy JHK broadband filters. We calculate the residual between the MLP prediction ŷ and the AT2017gfo observed data d for every band B by way of the reduced-χ 2 statistic
In our χ 2 residual calculation, we include observational uncertainties from the AT2017gfo data σ d , as well as systematic uncertainties σ sys which we use as a catch-all term to encompass all uncertainties, quantifiable or otherwise, associated with the neural network interpolation process. As outlined above and as discussed in greater quantitative detail in Section III C, we adopt a systematic modeling uncertainty, σ sys , of 0.5 magnitudes for our inference analysis.
For inference, we adopt a purely Gaussian likelihood based on these residuals
where N is the number of observations. As our inference is performed via adaptive Monte Carlo integration, the reliability of our posterior can be expressed in terms of a number of effective samples n eff . [Several different conventions exist for this number; see the appendix of [74] for discussion.] For this study, we terminate our analyses when n eff ≃ 10 3 . To validate our inference strategy, we constructed random synthetic sources, with kilonova model parameters drawn from our prior (albeit adopting a uniform rather than gaussian angular prior) and using observational times and uncertainties precisely matching the AT2017gfo cadence and instruments. For our synthetic sources, the expected light curve is generated using our neural network. As demonstrated with one example in Figure 4, our inferences are always consistent with the injected kilonova parameters. To demonstrate that our implementation retains statistical purity, we also used 100 random synthetic sources to perform a conventional probability-probability (PP) test, available in an Appendix.
While we employ the fixed σ sys = 0.5 for most of our studies, in order to validate our results we also perform a few selected analyses with different choices, on the one hand adopting different discrete choices and on the other treating σ sys as a continuous unknown model parameter.
The analysis in Figure 4 demonstrates that, if the underlying model is correct, comparison with AT2017gfolike observations should very tightly constrain each of this model's parameters. This fiducial result thus has qualitatively different behavior than our and others' prior analyses of AT2017gfo, where posterior inferences arrive
logmd = 1.986 +0.028 0.027 0 .1 8 0 0 .1 9 5 0 .2 1 0 0 .2 2 5 0 .2 4 0 vd [c] vd = 0.210 +0.006 0.007 1 .9 5 1 .9 0 1 .8 5 1 .8 0 logmw [M ] logmw = 1.865 +0.018 0.024 0 .1 3 0 0 .1 3 5 0 .1 4 0 0 .1 4 5 0 .1 5 0 vw [c] vw = 0.141 +0.002 0.002 2 .1 6 2 .1 0 2 .0 4 1 .9 8 1 .9 2 logmd [M ] 8 4 8 7 9 0 9 3 9 6 [deg] 0 .1 8 0 0 .1 9 5 0 .2 1 0 0 .2 2 5 0 .2 4 0 vd [c] 1 .9 5 1 .9 0 1 .8 5 1 .8 0 logmw [M ] 0 .1 3 0 0 .1 3 5 0 .1 4 0 0 .1 4 5 0 .1 5 0 vw [c] 8 4 8 7 9 0 9 3 9 6 [deg] = 87.528 +1.868 1.389
FIG. 4. Posteriors derived from a single randomly-generated synthetic kilonova source are consistent with its assumed model parameters. For this analysis, we have adopted the "zero noise" realization, where the kilonova light curve is precisely equal to its expected value.
at much broader posterior intervals, as discussed below. That said, the posterior shown above is consistent with the standard Fisher matrix estimate of the inverse covariance matrix Γ = Σ -1 , derived for example by taking (the expected value over noise realizations of) a secondorder Taylor series expansion of the log-likelihood as ln L = ln L o -2 -1 Γ ab (x -x * ) a (x -x * ) b where x * are the true synthetic parameters:
In this expression, only first-order derivatives appear because we assume that the model has no systematic bias such that ⟨d⟩ = ŷ; however, this simple estimate also arises inevitably using the large-amplitude "linear signal approximation" [75]. This Fisher matrix can be estimated to order of magnitude by replacing the derivatives ∂ ŷB /∂x a by the ratio ∆y B /∆x a , which for the mass parameters we approximate as 2/2 = 1, so the Fisher matrix is approximately Γ ≃ N/σ 2 sys and the posterior in each mass hyperparameter should have a one-standarddeviation scale of order 1/ √ Γ ≃ σ sys / √ N ≃ 0.5/ √ 333 ≃ 0.027.
The posterior distributions for our input parameters ⃗ x θ are plotted in Figure 5. The black posteriors represent 0.125 0.5 1 2 4 8 16 32 64 t (days) 17.5 20.0 22.5 25.0 27.5 30.0 AB Mag g + 7 y + 3 r + 6 J + 2 i + 5 H + 1 z + 4 K + 0 FIG. 6. Light curves generated by the MLP for the median parameters presented in Figure 5 (lines with shaded regions) with AT2017gfo observational data overplotted (scatter points). Despite the different values recovered between this analysis and that of R22, especially for the wind ejecta parameters and θ, the MLP prediction is able to replicate the observations surprisingly well.
the parameters identified in this study. As a direct comparison to the results of R22, we overplot the posterior distributions from that study in red. Despite using the same training data and comparing to the same observational data, the Gaussian process emulator and the MLP emulator recover posterior distributions with substan-tially different median values and distribution widths.
The dynamical ejecta parameters log m d and v d are similar between the two emulators, but the wind ejecta parameters log m w and v d , as well as the viewing angle θ, are quite different, with two of the three parameters inferred by the MLP residing outside of the GP inference 1σ limits. We verify the fidelity of the MLP inference results by generating light curves corresponding to the parameter values identified at the top of each column in Figure 5. These light curves are shown in Figure 6 and indicate that, assuming a 0.5 magnitude systematic uncertainty, the inferred MLP parameters do indeed replicate the AT2017gfo data to a reasonable degree of accuracy. While the majority of data is well replicated, the earlytime g and r bands and the late-time J and H bands deviate slightly outside of our uncertainty bands.
The recovery of θ ≈ 6 • is surprising for several reasons. First, the recovered angle was modestly offset from the Gaussian adopted as our inclination prior. Although different studies find a variety of viewing angles associated with AT2017gfo [10,16,17,[69][70][71]76], none indicate that the viewing angle is as low as our inference suggests. Second, and more importantly, for angular binning described by Equation 2, the first angular bin encompasses all emitted photons for viewing angles ∼ 0 -16 • . Therefore, by recovering a narrowly-peaked posterior around θ ≈ 6 • , the MLP seems to indicate that it can identify angular variations within a single angular bin at a resolution much finer than what is provided by training data. In other words, these inference results are either overly constrained, or the MLP is able to identify fine angular variations in the light curves when trained on radiative transfer simulations using a coarser angular grid.
One conceivable explanation for the narrow posterior distribution seen in Figure 5 is an underestimate of the underlying systematic error. To investigate this possibility, Figure 7 shows the results of inferences performed when adopting dfferent choices for the white-noise systematic error parameter σ sys , adopting both discrete and continuously-distributed choices for this parameter. In all scenarios, we infer similar parameters for AT2017gfo, even though we allow for several magnitudes of potential systematic uncertainty. Conversely, this direct comparison between our models and the data directly infers a value for our systematic uncertainty parameter consistent with our fiducial choice.
The extremely narrow posterior distribution in kilonova parameters and angle motivates a focused investigation of our training simulations and MLP model in the neighborhood of that posterior. As a first step, we performed followup SuperNu simulations at the inferred parameters using a higher angular resolution. Specifi-
logmd = 1.711 +0.130 0.136 0 .0 6 0 .1 2 0 .1 8 0 .2 4 0 .3 0 vd [c] vd = 0.132 +0.017 0.014 2 .2 2 .0 1 .8 1 .6 1 .4 logmw [M ] logmw = 1.690 +0.062 0.069 0 .0 8 0 .1 2 0 .1 6 0 .2 0 vw [c] vw = 0.095
+0.009 0.008 2 .8 2 .4 2 .0 1 .6 1 .2 logmd [M ] 8 1 6 2 4 3 2 [deg] 0 .0 6 0 .1 2 0 .1 8 0 .2 4 0 .3 0 vd [c] 2 .2 2 .0 1 .8 1 .6 1 .4 logmw [M ] 0 .0 8 0 .1 2 0 .1 6 0 .2 0 vw [c] cally, we increased the resolution by a factor of four to get a total of 216 angular bins. By reducing the number of photon packets in each angular bin by a factor of four, we also increase statistical noise by a factor of two. To ensure that our finer angular resolution analysis is not affected by this increase in statistical noise, we compare three separate simulations in Figure 8. The
8 10 12 14 16 M AB, y 54 n 1/4 216 0.125 0.5 1 2 4 8 16 32 64 Time (days) 0.1 0.0 0.1 M AB, y n 1/4 54 216 54
FIG. 8. Plots of SuperNu y-band light curves for a simulation like the ones described in Section II A (θ54), a simulation like θ54, but with one-quarter as many photon packets (n 1/4 ), and a simulation like θ54, but with four times greater angular resolution (θ216). The top panel indicates that, on a macroscopic scale, the simulations are identical. The bottom plot indicates that deviations do exist, likely attributed to statistical noise from the increase in angular bins, or matching reduction in packet count. The ejecta parameters used to create these simulations are those presented in Figure 5.
blue curve, labeled θ 54 , shows a SuperNu simulation using the parameters from Figure 5, hereafter x MLP and a standard 54-bin angular grid, as used in the training data simulations. The orange curve, labeled n 1/4 , shows a SuperNu simulation with the same exact parameters as θ 54 , except it uses one-quarter as many photon packets in the simulation. If the statistical noise described above were significant, the noise in n 1/4 should be much more pronounced than in θ 54 . Finally, the green curve, labeled θ 216 , shows a SuperNu simulation using the same parameters as θ 54 , but with a 216-bin angular grid, representing a factor of four increase in resolution. All three light curves show AB magnitude in the y-band as a function of time in days. As seen in Figure 8, our followup simulations agree with one another, with small Monte Carlo error comparable to our initial estimate and small compared to our adopted systematic uncertainty (σ sys = 0.5).
We then compare our high angular resolution simulation θ 216 to the predictions of the MLP, computing the difference between model and prediction at all times and simulation angles. Figure 9 shows the residual ∆M AB in the y-band when we take the absolute difference between ŷMLP and y sim , capped at a maximum difference of 1 magnitude. The residual values ∆M AB are initially evaluated for the 216 discrete angular bins; for visual clarity and diagnostic power, we linearly interpolate ∆M AB across time t and angle θ using the RegularGridInterpolator function from the scipy.interpolate library [77]. Figure 9 exhibits both expected and unexpected behaviors. The large mismatch at early times when photon count is FIG. 9. Deviations in y-band AB magnitude (∆MAB) between a SuperNu simulation with four times the usual angular resolution (ysim) and the MLP prediction for the parameters reported in Figure 5 (yMLP) as a function of time t and angle θ. Notably, this comparison between our MLP (trained at low angular resolution) and the followup simulation (performed at high angular resolution) does not exhibit strong small-scale variation along the θ direction (i.e., between adjacent angular bins in the original training set and outside the original training resolution), suggesting that the MLP correctly interpolates to smaller angular scales. low, particularly as θ approaches 90 • where the higher opacity dynamical ejecta further reduces photon count, is to be expected. The MLP is not fitting the light curves in this region well as there are too few photons available in the training data. However, the low mismatch (i.e. good fitting) dark blue regions in the plot indicate some sort of structure in the MLP's underlying ability to reproduce light curves at different times and angles. While we only include the plot of ∆M AB for the y-band in Figure 9, it is worth noting that the same behavior can be observed in these plots for all bands; predictable, low photon count systematics are identified in expected regions, but other, unexpected regions of increased systematic error manifest in different regions of the parameter space. Overall, however, the systematic uncertainties illustrated here are consistent with our expectations from the reported validation loss: a conservative systematic error of √ M SE val ≃ 0.5 reflects our overall uncertainty. For this reason, we adopted this systematic uncertainty in our parameter inferences above. This systematic uncertainty is substantially more conservative than the uncertainties adopted in R22.
The investigations in Section III C have introduced a surprise. On the one hand, our neural network reliably reproduces its training and validation data, including fol-lowup off-sample simulations performed at higher resolution. Though not shown here, we have also confirmed that the neural network agrees well with the surrogate provided in R22, using a small sample of randomlyselected ejecta parameters. On the other hand, the inferences obtained in Section III C by comparing all AT2017gfo kilonova observations to our MLP produce strikingly different results than R22. However, as noted in R22 and other works, most investigations have some tension between their models and the data, particularly in the bluer bands. In this section, motivated by this discrepancy, we also examine the effects on our AT2017gfo parameter inference when we use only specific subsets of the observational data. We perform two additional parameter inference runs using two categories of broadband data subsets: blue bands represented by the griz data and red bands represented by the yJHK data.
The posteriors in Figure 10 show how these bandlimited results compare with each other, the results of R22, and the all-band analysis presented in Section III C. The most apparent result is that the angle prior is recovered in both sets of posteriors, and both cases match the R22 results well. The other interesting feature is that the red yJHK posterior matches the R22 results much more closely than the blue posterior. Even when using a smaller subset of the observational data, the blue posteriors remain narrowly peaked in the parameter space, while the red posteriors become broader. The narrowness of the blue posteriors indicates that the blue broadband data determines the overall shape of the posteriors in the full broadband data inference.
The seemingly disproportionate effect of the blue broadband data on the posteriors could be attributed to the rapid evolution of the bluer bands compared to the red bands. As can be seen in Figure 6, the evolution of the griz light curves is more rapid than the yJHK bands, with the griz light curves dimming by ∼ 4 magnitudes compared to the yJHK bands dimming by two magnitudes during the first 10 days of observations. As such, the griz light curves will be more restrictive regarding which model parameters fit the data, which is evident from the Figure 10 posteriors.
To summarize, following R22 we fit the same simulations and performed comparable inference of AT2017gfo. After assessing our training data and fit systematics, we adopted more conservative systematic uncertainties than R22. We nonetheless find dramatically narrower posteriors than R22, with inferred light curves consistent with observational predictions.
Several possible reasons for the qualitative discrepancies arise, primarily pertaining to our systematic uncertainty estimate. In R22, the Gaussian process interpolation provided a pointwise and parameter-dependent error estimate, which we qualitatively verified across the pa- rameter space. Also, the fitting strategy adopted in R22 independently fit each timestep. As a result, we expect that R22's models are unlikely to have correlated systematics in time, angle, and wavelength. For example, the R22 light curves occasionally have small but notable random discontinuities, consistent with their reported fitting uncertainty. By contrast, in this work, we fit all times and angles together within the same training set to generate a vector prediction. Our approach does not presently provide a pointwise error estimate. Therefore, for the method adopted in this work, we expect correlated errors in time and angle, but lack a method to characterize them versus those parameters or even the intrinsic kilonova parameters. Our MLP's vector lightcurve prediction necessitates fitting the entire light curve for a certain band given a single set of parameters. In addition, as mentioned in Section III D, the blue bands are particularly constraining due to their more significant evolution over the observation period. Combining both of these effects results in only an extremely narrow region of the parameter space fitting all broadband data consistently. We did not encounter such narrow posteriors in R22 as each prediction was made for a specific time point; as such, many samples could reasonably fit an observation at any given time, resulting in broader posteriors when all times were stitched together to form the light curve.
We use the validation loss value to roughly estimate the systematic fitting error associated with our neural network outcomes: the validation curve in Figure 2 suggests that differences of order √ M SE ≃ √ 0.2 magnitudes should occur in our predictions. In practice, as illustrated in Figure 7 below, we find that the average squared systematic error suggests a larger value than the validation MSE. We therefore anticipate that our naive estimate of σ sys = 0.5, though well motivated by our detailed followup study, may still understate the systematic uncertainty inherent in our fitting approach, resulting in narrow posterior distributions. We also emphasize that the differences are not simply a matter of scale: the investigations performed in Figure 7 suggest that larger white-noise systematic error cannot reconcile differences between our current analysis and previous results.
We note that we have experienced similar systematic uncertainty associated with observations in blue bands in R22 and an associated inference using simulations of spectra [54]. The systematics in those works were related to our inability to reproduce the observed blue flux at times past ∼ 2 days using our best-fit simulations. The systematics in this work, though similar in their connection to blue observations, introduce slightly different effects in our resultant inference. As we solve the bigger problem of matching our simulations to late-time blue observations, we anticipate that a more sophisticated treatment of our emulator systematics will allow us to better understand the effects of blue-band data on our inferences.
A thorough investigation of suitable fit systematics for this neural network is well beyond the scope of our study. In the meantime, the neural network is suitable for investigations such as the one presented in Figure 10, where we can examine our models' ability to fit certain subsets of the data. In the griz case, we see that our models require over 0.1 M ⊙ of slow-moving dynamical ejecta to fit the blue data. But, we expect dynamical ejecta to be less massive and faster, thus potentially suggesting a missing modeling component.
We present a neural network architecture that is useful for the interpolation of kilonova light curves. We report on the network's training and validation loss as a metric of successful training, as well as present examples of off-sample light-curve recovery. We use the neural network to infer the parameters of the AT2017gfo kilonova and compare to previous inference performed in Ref. [50]. We find that the inference results are quite different from those previously obtained, but the light curves generated by the recovered parameters align well with the observational data. In particular, we investigate the neural network's ability to seemingly infer narrow regions of the angle space despite being trained on light-curve data that should not allow for such specific inference. Given a detailed analysis of the mismatch between the neu-ral network's predictions and a simulation with higherresolution angular data, we find that the network's pointwise systematic errors are consistent with our error estimate. However, our investigations also suggest that the systematic errors are correlated, not independent, in time and angle, in a way that is not captured by our model for systematic uncertainties. In other words, we have discovered that the neural network's goodness-offit varies appreciably across the time-angle space. While some of these variations are expected, others form interesting features that we cannot readily explain. We leave the analysis of the interpretability of these features for a future investigation.
We also show that the systematic uncertainty may be more complex than assumed in our simple uncorrelated (white noise) error model. This was not the case in R22 due to the interpolation uncertainty, which naturally stemmed from the Gaussian process methodology. In recovering different parameters for AT2017gfo using two emulators trained on the same library of simulations, we highlight the importance of quantifiable uncertainty analysis in using emulators for robust inference. As we do not present a way to handle correlated uncertainties in this work, a detailed uncertainty analysis, along with the resultant effects on parameter inference, will be necessary in future work.
To validate the statistical purity of our brute-force inference technique and our understanding of the noise model, we constructed a standard probability-probability (PP) plot test [78,79]. Our description follows the notation and narrative used in [74]. For each source k, with true parameters λ k , we calculate the fraction of its posterior distribution with parameter λ α below the true source value λ k,α [ Pk,α (< λ k,α )]. After reindexing the sources so that Pk,α (λ k,α ) increases with k for some fixed α, a plot of k/N versus Pk (λ k,α ) can be compared with the expected diagonal result (P (< p) = p) and binomial uncertainty interval. Figure 11 shows the PP plot derived using kilonova light curves generated with our neural network interpolator. In these analyses, we adopt precisely the same observation cadence and uncertainties as AT2017gfo. As in our fiducial analysis of AT2017gfo, we adopt σ sys = 0.5. Each synthetic observation incorporates both observational and (white noise) systematic uncertainty, added in quadrature consistent with our assumed likelihood.
The PP plot in Figure 11, being sufficiently consistent with the binomial credible interval, suggests that the brute-force Monte Carlo inference strategy adopted in this work suffices for our purposes: in short, that the qualitative extent and character of the posteriors shown in our figures are reasonably correct, such that the considerable tension between our analysis and previous inferences accurately reflects the posterior. We have specifically chosen to present results from a brute-force inference technique to circumvent debates about our choice of implementation or our method of assessing convergence. With AT2017gfo, we have confirmed that the choice of integrator also doesn't qualitatively change our answer: alternative brute-force Monte Carlo integrator implementations produce similar results. That said, the PP plot above is clearly not as diagonal as would expected for a well-developed and calibrated Bayesian inference algorithm applied to this problem: its S-shape features suggests either modest overdispersion in our synthetic error model or modest underdispersion in our posterior distributions.