XKN: a Semi-analytic Framework for the Modelling of Kilonovae

After GW170817, kilonovae have become of great interest for the astronomical, astrophysics and nuclear physics communities, due to their potential in revealing key information on the compact binary merger from which they emerge, such as the fate of the central remnant or the composition of the expelled material. Therefore, the landscape of models employed for their analysis is rapidly evolving, with multiple approaches being used for different purposes. In this paper, we present xkn, a semi-analytic framework which predicts and interprets the bolometric luminosity and the broadband light curves of such transients. xkn models the merger ejecta structure accounting for different ejecta components and non-spherical geometries. In addition to light curve models from the literature based on time scale and random-walk arguments, it implements a new model, xkn-diff, which is grounded on a solution of the radiative transfer equation for homologously expanding material. In order to characterize the variety of the ejecta conditions, it employs time and composition dependent heating rates, thermalization efficiencies and opacities. We compare xkn light curves with reference radiative transfer calculations, and we find that xkn-diff significantly improves over previous semi-analytic prescriptions. We view xkn as an ideal tool for extensive parameter estimation data analysis applications.


INTRODUCTION
The detection of electromagnetic counterparts of gravitational wave signals represents one of the key aspects of gravitational wave astrophysics and, more in general, of multimessenger astronomy.While the gravitational wave signal produced by a coalescing compact binary encodes many properties and information about the merging system (e.g. the chirp mass, the masses of the two compact objects or their tidal deformation, if at least one of the two is not a black hole), the electromagnetic signal can provide complementary information, including for example the amount of matter expelled during the merger and its chemical composition.Other aspects of the merger, such as the nature of the coalescing objects or of the remnant that forms after the merger, could affect both the gravitational and the electromagnetic emission.In this case, the presence of more than one signal can provide tighter constraints and help discriminating between ambiguous or degenerate situations (see e.g.Radice & Dai 2019;Hinderer et al. 2019;Barbieri et al. 2019bBarbieri et al. , 2021;;Dietrich et al. 2020;Raaĳmakers et al. 2021).
The potential of multimessenger astrophysics was recently re- ★ Contact e-mail: giacomo.ricigliano@gmail.comvealed by GW170817 (Abbott et al. 2017a(Abbott et al. ,b, 2019)).Just from the analysis of the gravitational wave signal, it was impossible to exclude that the coalescing objects were two black holes, since the posterior of the binary tidal deformability was extending down to 0 (Abbott et al. 2017a(Abbott et al. , 2019)), the value expected for a binary black hole.The identification of the system as a binary neutron star was mostly based on the values of the masses of the merging bodies and, more importantly, on the detection of two electromagnetic counterparts, namely a short gamma-ray burst and a kilonova (Chornock et al. 2017;Cowperthwaite et al. 2017;Coulter et al. 2017;Drout et al. 2017;Evans et al. 2017;Goldstein et al. 2017;Hallinan et al. 2017;Kasliwal et al. 2017;Murguia-Berthier et al. 2017;Nicholl et al. 2017;Smartt et al. 2017;Soares-Santos et al. 2017;Tanvir et al. 2017;Troja et al. 2017;Villar et al. 2017;Waxman et al. 2018;Ghirlanda et al. 2019).Such a combination of signals was indeed very useful in providing constraints on the equation of state of nuclear matter or in shedding light on the central engine of gamma-ray bursts (see e.g.De et al. 2018;Abbott et al. 2018;Coughlin et al. 2019).On the other hand, in the case of the subsequent binary neutron star merger, GW190425 (Abbott et al. 2020), or in the first observed black hole-neutron star systems (Abbott et al. 2021a,b), no electromagnetic counterparts were observed (see e.g.Coughlin et al. 2020).In these cases, the nature of the binary was deduced only from the inferred masses, while the gravitational signal alone was not informative on the tidal deformability of the system, due to the lower expected values and to the not sufficiently high signal-to-noise ratio.Among the different electromagnetic counterparts, the kilonova is one of the most peculiar transients associated to compact binary mergers involving at least one neutron star (Li & Paczynski 1998;Metzger et al. 2010).It arises when the merger and its remnant expel a non-negligible amount of neutron-rich matter, which undergoes -process nucleosynthesis (Lattimer & Schramm 1974;Eichler et al. 1989;Freiburghaus et al. 1999); see also Cowan et al. (2021); Perego et al. (2021) for recent reviews.The decay of the freshly produced radioactive elements moving from the neutron-rich side towards the bottom of the nuclear valley of stability releases nuclear energy that, despite the fast expansion, keeps the expanding ejecta hot.Due to expansion the matter opacity to the electromagnetic radiation decreases until photons can eventually diffuse out, producing a kilonova (see Metzger 2020, for a recent review and references therein).Depending on the mass of the ejecta, on their expansion velocity and composition, the peak of the kilonova emission is expected to occur between a few hours and several days after the merger (see, for example, Arnett 1982;Pinto & Eastman 2000;Barnes & Kasen 2013;Tanaka & Hotokezaka 2013;Radice et al. 2018;Kawaguchi et al. 2020).At the same time, the emission is expected to evolve, moving from bluer to redder frequencies as a consequence of the photospheric expansion, of the decrease in the nuclear energy input and in the opacity of matter, as well as of the viewing angle (see e.g.Metzger & Fernández 2014;Perego et al. 2014;Martin et al. 2015;Rosswog et al. 2017;Wollaeger et al. 2018;Kasen & Barnes 2019;Korobkin et al. 2021).As long as the opacity of the innermost ejecta is large enough, the kilonova is characterized by the presence of a photosphere and the resulting emission can be described, in good approximation, as quasi-thermal.Non-thermal and non-local thermodynamics equilibrium effects become more and more relevant as time increases, until the transient enters its nebular phase.In the case of AT2017gfo (the kilonova associated to GW170817) the transition from a full photospheric regime to the nebular phase happened between a few days to a week after merger (see, for example, Smartt et al. 2017;Watson et al. 2019;Wu et al. 2019;Pognan et al. 2022a,b;Gillanders et al. 2023).
The modelling of kilonovae is extremely challenging.It requires the solution of a radiative transfer (RT) problem in a fast expanding, radioactive and radiation dominated medium.Not only the composition of matter changes with time due to nuclear reactions and decays, but due to the expansion and to the interaction between matter and radiation, atoms inside the ejecta (which are initially fully ionized due to the large matter temperature) span different ionization levels, following the progressive electron recombination.The presence of heavy elements, and in particular, of lanthanides and actinides, largely increase the photon opacity due to bound-bound and boundfree transitions involving the  electron shells (Roberts et al. 2011;Kasen et al. 2013;Tanaka & Hotokezaka 2013).For most of the heavy elements, opacities due to ionized species and for matter in the thermodynamics regime relevant for kilonovae are experimentally unknown and their values are usually provided by non-trivial atomic structure calculations (see, for example, Fontes et al. 2020;Tanaka et al. 2020;Banerjee et al. 2022Banerjee et al. , 2023)).Furthermore, large uncertainties still affect the calculation of the detailed nuclear energy released by -process elements (Rosswog et al. 2018;Barnes et al. 2020;Zhu et al. 2020Zhu et al. , 2022;;Lund et al. 2023), as well as the estimation of the fraction of the energy that the expanding matter is able to thermalize (Hotokezaka et al. 2016;Barnes et al. 2016;Kasen & Barnes 2019).In addition to the difficulties related with the problem of transporting photons inside an expanding, radioactive medium, an additional challenge is represented by the fact that a binary neutron star merger or a black hole-neutron star merger can expel matter with different properties and, possibly, with a high degree on anisotropy (see, e.g., Wanajo et al. 2014;Just et al. 2015;Sekiguchi et al. 2016;Foucart et al. 2016;Perego et al. 2017).This implies that the medium inside which the photons are produced, diffused and emitted can have a non-trivial stratification, as well as angular distribution.
It is not surprising that, given the complexity of the kilonova scenario and the variety of aspects involved, so far the problem of predicting or producing kilonova light curves and spectra has been tackled by a large variety of models, employing very different levels of sophistications and approximations.Some models solve the photon transport problem in an expanding medium considering wavelenght-and composition-dependent opacities, computed consistently and coupled to the calculation of the abundances of the different ion species, assuming local thermodynamics equilibrium (see, e.g., Tanaka & Hotokezaka 2013;Kasen et al. 2017;Wollaeger et al. 2018;Shingles et al. 2023).These models are the most sophisticated and reliable ones, but they necessarily require large computational resources, especially in three dimensions.Other examples of kilonova models include TARDIS (Kerzendorf & Sim 2014), which solves the 1D photon transport problem in the optically thin atmosphere above a predefined photosphere, POSSIS (Bulla 2019(Bulla , 2023)), which uses pre-computed wavelenght-and time-dependent opacities on a 3D Cartesian grid, or SNEC (Morozova et al. 2015;Wu et al. 2022), which solves radiation hydrodynamics in spherical symmetry through a gray flux-limited diffusion approach.These more approximated approaches clearly reduce the computational effort, especially if some symmetry is invoked.
At the opposite extreme, kilonova light curves have also been computed by using simplified kilonova models that avoid the direct solution of the RT problem, since they are often based on the solution of the energy conservation equation inside the ejecta or on time scale arguments mimicking the mean features of the photon diffusion problem (see, e.g., Grossman et al. 2014;Martin et al. 2015;Hotokezaka & Nakar 2020).They usually employ gray, constant opacities and can reproduce some of the most relevant features of the kilonova emission, at least at a qualitative level.The extremely reduced computational costs of these models allows their usage in multi-dimensional parameter estimate analysis, which requires the evaluation of millions, if not billions of kilonova light curves (as done, e.g., by Breschi et al. 2021).
With the increase of the number and sensitivity of gravitational wave detectors, the amount of multimessenger signals, and in particular, of kilonova counterparts of gravitational wave events, is expected to significantly grow in the years to come.For example, during the fourth observational campaign of LIGO, Virgo and KAGRA, the number of detected binary neutron star merger is expected to be a few tens per year (Petrov et al. 2022;Colombo et al. 2022).Additionally, the careful (re)analysis of the afterglow signals of close-by gamma-ray bursts can reveal signatures of kilonova emission, as in the case of the exceptionally bright GRB211211A (Troja et al. 2022), or more recently also in the case of the long GRB230307A (Levan et al. 2023).Given the present scenario, characterized by a growing number of kilonovae, which could significantly differ in terms of intrinsic properties, as well as in the quality and quantity of the data, it is still imperative to improve on the accuracy of fast and approximated kilonova models.The latter can be complementary to more sophisticated models, since they can be used for extensive parameter estimations, and to provide a robust and reliable framework to analyze coherently kilonova emissions coming from very different events.They can also be coupled to gravitational wave data analysis in the quest for coherent and genuine multimessenger analysis.
In this paper, we present xkn, a semi-analytic framework to perform analysis of kilonova emission, both in terms of bolometric luminosity and broadband light curves.xkn inherits the possibility of adding several ejecta components and of prescribing non-trivial ejecta geometries from previous implementations (Martin et al. 2015;Perego et al. 2017).However, with respect to the latter, it aims to improve on the accuracy and on the reliability of the resulting light curves by replacing the kilonova model grounded on time scale arguments with a different model, xkn-diff, based on a semi-analytic solution of the diffusion equation for homologously expanding ejecta.Such a solution was presented in Wollaeger et al. (2018), and based on works reported in Pinto & Eastman (2000).Here, we expand the class of solutions, by considering time dependent heating rates, thermalization efficiencies and opacities.Moreover, we improve the physical input, by including composition dependent heating rates and opacities.
The paper is structured as follows: in Sec. 2 we present in detail the spherically symmetric, kilonova emission model xkn-diff, distinguishing between the optically thick (Sec.2.1) and the optically thin (Sec.2.2) part.The general multi-component, anisotropic framework of xkn to compute light curves is presented in Sec. 3, while in Sec. 4 we detail the input physics entering the model, listing the heating rates (Sec.4.1) and the opacity (Sec.4.2) prescriptions.In Sec. 5 we compare the results of the various xkn models with the ones obtained using a RT code (Tanaka & Hotokezaka 2013;Kawaguchi et al. 2018Kawaguchi et al. , 2021)), taken as reference, in order to address the degree of accuracy and the limitations of our approach.We provide a summary and the conclusions of our work in Sec. 6.

SEMI-ANALYTIC 1D KILONOVA MODEL
Our kilonova model is based on a one-dimensional model for the diffusion and emission of photons from homologously expanding, radioactive matter.More specifically, the kilonova emission is calculated as the combination of two contributions, one emitted at the ejecta photosphere, i.e. the surface delimiting the optically thick bulk of the ejecta and from which photons can escape and move inside the atmosphere, and one coming from the optically thin layers above it.In the following, we separately present these two contributions.

Optically thick ejecta treatment
The contribution to the luminosity arising from the photosphere is computed starting from the semi-analytic formula originally proposed by Pinto & Eastman (2000) with the intent to treat the ejecta from Type Ia supernovae (SNe Ia), and later adapted by Wollaeger et al. (2018) to model kilonovae.In spite of its simplifying assumptions, this formula can qualitatively reproduce the thermal evolution of the ejecta.In the following, we report in broad lines its derivation.
The ejecta fluid is assumed to be optically thick throughout its entire depth and the radiation field properties are evolved on the basis of the time-dependent equation of RT.In particular, we consider the first two frequency-integrated moments of such equation in the comoving frame, calculated to order  (/): 1 where  is the fluid velocity field,  the radial coordinate,   the extinction coefficient and   the volume emissivity, while the subscript  indicates the frequency dependence.Here ,  and  are the energy density, flux and radiation pressure of the radiation field, respectively./ indicates the comoving (Lagrangian) derivative.The solution of this set of equations is found by adopting a series of hypotheses.
(i) We assume a homologous expansion of the fluid, i.e. each fluid element expands with constant radial speed.Under this assumption, the ejecta maintain their proportions while expanding with an external radius () ≃  max , where  max is the maximum outflow velocity.The homologous expansion hypothesis is consistent as long as the energy heating the fluid does not affect the fluid motion in a significant way.
(ii) We resort to the Eddington approximation, wherein the radiation field is isotropic and the relation  = 3 holds.The latter is valid since the outflow is optically thick, as one can expect at least in the early stages of its evolution.When the fluid becomes transparent at later times, even if this approximation breaks down, the error has still little effect on the bolometric luminosity.
(iii) Regarding the energy balance, we assume that the gas internal energy, that is its thermal kinetic energy as well as the ionization energy, is subdominant with respect to the radiation field energy, and therefore we ignore the former (radiation dominated conditions).
(iv) Additionally, we assume that the absorbed heat is immediately re-radiated as thermal emission, and thus: where  heat is the energy deposition rate per unit volume.
In light of the above considerations, we act on Eq. (1) and Eq. ( 2) with the aim to simplify them.Eq. ( 1) is an equation for the energy density field,  (, ), and in order to ensure the correct radiation energy balance we retain all terms to the order  (/), as  (, ) changes considerably on the hydrodynamic timescale.Eq. ( 2) is instead an equation for the radiation momentum  (, ) and we solve it at lower order by discarding all time and velocity-dependent terms.This choice is appropriate on the fluid-flow timescale as we assumed that  (, ) is not relevant for the outflow dynamics.Hence by inverting Eq. ( 2) we obtain: with  a properly defined frequency-averaged inverse mean free path.
The above expression for the flux can be inserted in Eq. (1), resulting in: where  denotes the derivative with respect to time.We now introduce  as an absorption opacity, homogeneous in space (but not necessarily in time), such that  = , with  the fluid density.Moreover, we express  heat as  heat =   th , where  is the energy release rate per unit mass and  th a thermalization efficiency coefficient, and both are function of time.From the assumption of homologous expansion, we recall that: where  ∈ [0, 1] is the dimensionless radius coordinate.Moreover, in the radiation transport equation we adopt the following single-zone homologous solution for the expansion profile: with  0 the initial time of the expansion, and  0 the density at  0 .We approximate the latter as: where  ej is the ejecta mass.Furthermore, from the hypothesis of radiation dominated gas, we employ the polytropic equation of state in the Eddington approximation to obtain  ∝  −4 .If we assume that the residual dependences of  can be separated into a spatial profile () and a temporal profile (), we can write: where  0 is treated as a free parameter.In particular, assuming that the radiation field has a black-body spectrum, we relate it to an initial black-body temperature: with  = 4 SB / = 7.5657×10 −15 erg cm −3 K −4 being the radiation constant and  SB the Stefan-Boltzmann constant.Using Eq. ( 6), Eq. ( 7) and Eq. ( 9), the transport equation Eq. ( 5) becomes: with the prime superscript on a function indicating the derivative with respect to its variable and where The latter is a comprehensive factor comparable to the diffusion time scale that carries a possible dependence on  through .The homogeneous form of Eq. ( 11) can be solved by means of variable separation, according to which the resulting two functions in  and  must be identically equal to the same separation constant, : resulting in two ordinary differential equations to be solved.The equation for the spatial profile can be expressed as an eigenvalue equation for the operator : The eigenfunctions of Eq. ( 14) can be determined by imposing suitable boundary conditions to the problem.While for the temporal profile of the energy density  (, ) we naturally assume ( 0 ) = 1 and (∞) = 0, for the spatial part it is reasonable to consider a reflection symmetry at  = 0 and a radiative-zero condition at  = 1, being the outflow optically thick: Eq. ( 15) can be directly translated into identical conditions for the eigenfunctions of Eq. ( 14).If we impose the normalization requirement: as we adopt the notation: the resulting homogeneous spatial eigenfunctions are: with  =  2  2 and  any positive integer.
Eq. ( 18) represents a complete orthonormal basis on the interval of interest, and therefore we can expand the general solution of Eq. ( 11) on the latter: Here the coefficients   () retain a generic time dependence, and we can conveniently redefine them in the form: thus obtaining: Using Eq. ( 21) in Eq. ( 11), we can exploit the orthonormality of   () integrating over  ∈ [0, 1] to find: Finally, the bolometric luminosity is found by employing Eq. ( 4) and Eq. ( 21) to compute the flux at the surface of the ejecta: where   () are the solutions of Eq. ( 22), that can be obtained once the time-dependence of ,  th and  has been specified.If the formal solution of Eq. ( 22) is complex, when we compute the luminosity through Eq. ( 23) we take only the real part of it.Since this model is valid in the limit of optically thick matter, the outcome of Eq. ( 23) is rescaled by a factor  thick / ej , where  thick is the mass of the optically thick portion of ejecta, defined as the region enclosed by the photosphere: Differently from the single-zone approximation adopted in the solution of the RT equation, here we choose a more accurate spacedependent density profile such as the self-similar homologous solution (Wollaeger et al. 2018): The photospheric radius evolution  ph () can be found analytically by imposing the condition with   the optical depth of the material: and (, ) the density of Eq. ( 25).The determination of  ph () implies the solution of a seventh order polynomial equation.However, regardless of the ejecta parameters, the temporal evolution of  ph resembles a parabolic behaviour.Then, we approximate it as a parabolic arc with extremes fixed by the condition  ph = 0 applied to Eq. ( 26), i.e.: and curvature fixed by a third point,  3 , taken in the proximity of  1 , where Eq. ( 26) can be solved by assuming By adopting this approximate solution, the error on the photosphere position with respect to that provided by the exact solution of Eq. ( 26) is contained within 8%.
We characterize the emission at the photosphere by assuming a Plankian black-body spectrum, and thus we compute the associated photospheric temperature  ph () by means of the Stefan-Boltzmann law: with  thick () the luminosity of the thick part of the ejecta.A temperature floor  floor is applied in order to approximately account for the electron-ion recombination during the ejecta expansion (see e.g.Barnes & Kasen 2013;Villar et al. 2017).When  ph () reaches the temperature floor,  ph () is thus recomputed solving the implicit equation obtained from the Stefan-Boltzmann law.The value of  floor , generally treated as model parameter, is in fact dependent on the ejecta opacity and therefore closely linked to its composition.
For this reason, we also include the possibility to interpolate the floor temperature between two model parameters  Ni and  La , corresponding to characteristic recombination temperatures in a Lanthanidespoor (  ≳ 0.3) and a Lanthanides-rich (  ≲ 0.2) environment, respectively.
In the following, we discuss a few cases for which the solutions of Eq. ( 22) can be obtained analytically, based on temporal dependence of ,  and  th .

Constant opacity, constant thermalization efficiency and power-law heating rate
We first consider the case in which  is not only uniform in space, but also constant in time, i,.e. =  0 .In this case, the quantity  becomes a constant,  0 = 3 0  0 ( max  0 ) 2 /.Additionally, we consider  th to be a constant,  th,0 , while In this case, the solutions of Eq. ( 22) take the explicit form: where Γ(, ) is the upper incomplete gamma function, defined as: while   are constants: and   are integration constants fixed by the boundary conditions.In order to find the latter, we exploit the assumptions over the temporal profile of the energy density.While (∞) = 0 is automatically satisfied by the form of   (), one needs to translate ( 0 ) = 1 into a condition on   ().We choose to assign   ( 0 ) =  ,1 to ensure the convergence of Eq. ( 23):

Constant opacity, power-law thermalization efficiency and heating rate
The previous solution can be trivially generalized to the case in which both the specific heating rate and the thermalization efficiency follow a power-law evolution: In this case, Eq. ( 31)-Eq.( 33) are still a solution of Eq. ( 22), once  has been replaced by  ′ and  ′ ≡  + .

Constant opacity, constant thermalization efficiency and power-law heating rate with exponential terms
We then consider the case in which the opacity and the thermalization efficiency are constant in time, while the specific heating rate can be written as The solutions of Eq. ( 22) becomes: where erfi() ≡ − erf () is the imaginary error function and the error function is defined as The   ,   and   coefficients read This solution can be trivially generalized to the case in which more than one exponential term is added to the power law term in Eq. (36).

Power-law opacity, thermalization efficiency and heating rate
Finally, we consider the case in which the opacity, the thermalization efficiency and the specific heating rate have a temporal power-law dependence: In this case,  can be expressed as  =  0 (/ 0 ) − so that Eq. ( 22) becomes: In this case, the solution of Eq. ( 43) becomes: where   and   are defined as and

Optically thin ejecta treatment
At the times relevant for the kilonova, we expect the r-process material outside of the ejecta photosphere to provide a non-negligible contribution to the heating powering the emission, especially when a proper photosphere will eventually not be identifiable anymore.However this contribution is expected to be considerably different with respect to the one provided by the same portion of the ejecta if the latter were optically thick, i.e. if we assumed  ph () =  max ().
Therefore, in addition to the radiation emitted at the photosphere, we approximate the bolometric luminosity  thin () from the thin region outside of it, following the prescription of Grossman et al. (2014); Martin et al. (2015); Wu et al. (2019).We thus divide this region into  thin layers of equal mass d  assuming local thermodynamics equilibrium within each layer, and we express such contribution as: where the sum runs over the discrete thin shells, while  () and  th,i () are the specific radioactive heating rate and the spacedependent binned thermalization efficiency, respectively, as described in Sec.4.1.The total bolometric luminosity of the ejecta is simply obtained by summing the contributions from the two separate regions: Specific energy rate [erg s Bolometric luminosity per unit of ejecta mass and its two contributions from the optically thick and thin ejecta as a function of time, for a spherically symmetric model with  ej = 0.01  ⊙ ,  ej = 0.2  and  = 5 cm 2 g −1 .Also shown in red is the specific radioactive heating rate powering the emission and its thermalized fraction.The heating curve was obtained following the procedure described in Sec.4.1 using an entropy of  = 10  B baryon −1 and an expansion timescale of  exp = 5 ms.
In Fig. 1, the total luminosity is displayed together with its components for a simple spherically symmetric model.As visible, the early light curve is dominated by the thick ejecta, which constitutes the majority of the total mass.In this phase most of the energy provided by the radioactive decays is trapped within the ejecta due to the high optical depths.After a few days, the ejecta density has decreased enough for this energy to escape, enhancing the emission up to be instantaneously greater than the thermalized heating rate.Meanwhile, a second contribution to the luminosity steps in, as a relevant portion of optically thin mass emits radiation as well.Finally, after several days, the thick bulk of the ejecta disappears, and the luminosity is completely determined by the optically thin matter.Since the latter is transparent to thermal radiation, the thermalized decay energy escapes without further processing, and the luminosity is equal to the thermalized heating rate.However, the latter is now only a small fraction of the total decay energy rate, since the lower densities make the thermalization process inefficient.Thus, the thermal emission will eventually fade away.
Along with the temperature of the photosphere  ph , we want to characterize also the temperature in the layers outside of it.For this purpose, we assign to each bin a temperature   () on the basis of the radial profile proposed by Wollaeger et al. (2018) and derived for a radiation dominated ideal gas using Eq. ( 25): For each time, we fix the factor  0 () by requiring the continuity of the profile with the photospheric temperature  ph ().Therefore, for every bin  we obtain: where   is the position of the bin and  ph () =  ph ()/ max () is the position of the photosphere.
Recently, Pognan et al. (2022a) showed that a synthetic non-local thermodynamics equilibrium evolution of the temperature in the late expanding ejecta features a re-increase from a minimum reached around a few tens of days post merger.This result was obtained by taking into account the material excitation and ionization states in a more careful calculation of the ejecta heating and cooling processes, using the spectral synthesis code SUMO.In light of the above computations, we expect the temperature to remain roughly constant at the late times still relevant for the first kilonova phase.Therefore, in order to describe the thermal emission, here we find sufficient to set a unique time-independent minimum temperature value  floor for both the thick and the thin part of the ejecta.

MULTI-COMPONENT ANISOTROPIC SEMI-ANALYTIC KILONOVA MODEL
Ejecta from compact binary mergers are expected to occur in different components, characterized by different properties.Moreover, the ejection mechanisms can result in a anisotropic structure of the ejecta.Motivated by this, we set up a multicomponent, anisotropic kilonova framework.In particular, we closely follow the set-up first proposed by Perego et al. (2017), originally based on Martin et al. (2015) and reprised by Barbieri et al. (2019aBarbieri et al. ( , 2020)); Camilletti et al. (2022).
The framework assumes axial symmetry around the rotation axis of the binary (denoted as ), as well as reflection symmetry about the  = 0 plane.The polar angle  is discretised in a series of   bins which can be equally spaced either in the angle itself or in cos .A kilonova model is specified once the polar distribution of all the relevant quantities (i.e.mass, velocity, opacity or electron fraction, entropy, expansion time scale) are given for each of the ejecta components.Inside each angular bin and for each component, the radial kilonova model described in Sec. 2 (or alternatively the model from Grossman et al. (2014)) can be employed to compute the contribution to the luminosity emerging from that angular bin.Being an extensive quantity, the mass inside the bin needs to be scaled by the factor 4/ΔΩ, where ΔΩ is the bin solid angle.All the other input quantities are otherwise intensive and do not need any rescaling.Once computed, the isotropic luminosity resulting from the 1D model is scaled again based on the actual emitting solid angle, i.e. it is multiplied by ΔΩ/4.Within the same angular bin and in the presence of more than one ejecta component, the corresponding luminosity contributions are summed together, assuming that photons emitted from the innermost components irradiate the outermost ones and are subsequently re-processed and re-emitted on a time scale smaller than the expansion one.Moreover, at each time we locate the photosphere of the overall ejecta at the position of the larger individual photosphere.This approach assumes that the different components are nested and that they do not cross each other significantly during the kilonova emission.We expect these hypotheses to be appoximately verified once the homologous expansion phase has been reached and if the late time ejecta are systematically slower than the first expelled ones.
The present implementation includes characteristic analytic functional forms for the angular dependences: uniform distributions, step functions, sin , sin 2 , cos  and cos 2  dependences.Despite their simplicity, some of these distributions were demonstrated to be remarkably valid in broadly reproducing the outcomes of generalrelativistic hydrodynamical simulations, accounting for the preferential equatorial direction of the dynamical component, as well as the excursion in the electron fraction caused by high-latitudes neutrino irradiation.Additionally, the code can interpolate on its angular grid arbitrary distributions, such as azimuthally-averaged angular profiles extracted from numerical simulations (see, for example, Camilletti et al. 2022).
Typical kilonova models employed in the past used up to three different components (Perego et al. 2017;Breschi et al. 2021, see, e.g.,).In the case of two components, the fastest one usually refers to the dynamical ejecta, while the slowest one to the disc wind ejecta of viscous origin.A third, intermediate component is sometimes used, possibly originated by magnetic-(see e.g. ) or neutrino-driven wind components (e.g.Perego et al. 2014), as well as from spiral wave wind ejecta (e.g.Nedora et al. 2019Nedora et al. , 2021)).
In the source frame, the emission is assumed to be thermal and the spectral fluxes are described by a Planckian spectral distribution   (), i.e.: with  B the Boltzmann constant and ℎ the Planck constant, both at the photosphere as well as within each thin external layer.In the former case,  is the photospheric temperature, while in the latter it is the temperature inside each mass shell.
If the source is located at a luminosity distance   , corresponding to a redshift , for an observer on Earth characterized by a viewing angle  iew , the radiant flux at frequency  and time  (measured in the observer frame) will be the sum over the angular bins of the contributions from the thick ejecta  thick , and the thin ejecta  thin

𝜈,𝑘
(computed in the source frame), once the redshift correction has been applied to the time, frequency and luminosity: with: and: where  thick,k is the photospheric luminosity of the bin , characterized by a photospheric temperature  ph, .The factors   ( iew ) in Eq. ( 52) account for the effective emission area as seen by the observer (Martin et al. 2015), and are calculated using the formula: where ì ( iew ) is the unit vector in the observer direction, while ì   is the unit vector pointing radially outwards from the surface of the bin .Finally, we compute the AB magnitude at a photon frequency  as:  AB, () = −2.5 log 10 (   ()) − 48.6 . (56)

Heating rates
The heating rate powering the kilonova originates from the many decays of heavy elements produced in the r-process nucleosynthesis, and as such it can be computed by employing a nuclear reaction network.The latter calculates the time evolution of the nuclides abundances while keeping track of the energy released in the process.Results obtained by nuclear network calculations retain a strong dependence on the properties of the ejecta, and in particular on the entropy, electron fraction and expansion timescale at the freeze-out from nuclear statistical equilibrium (NSE) (see, e.g., Hoffman et al. 1997;Lippuner & Roberts 2015).Furthermore, nuclear network calculations also depend on the nuclear physics employed, e.g. on the choice of the theoretical nuclear mass model, the reaction rates or the fission fragment distribution.This sensitivity is particularly strong at low electron fractions and the nuclear physics uncertainties can lead to changes in the predicted heating rates of about one order of magnitude (Mendoza-Temis et al. 2015;Rosswog et al. 2017;Zhu et al. 2020).
With the purpose to provide our kilonova model with a heating rate valid for arbitrary initial conditions, we consider the results of the broad nucleosynthesis calculations reported in Perego et al. (2022).In that work, the nuclear composition evolution of a set of Lagrangian fluid elements is computed using the nuclear reaction network SkyNet (Lippuner & Roberts 2017) with the finite-range droplet macroscopic nuclear mass model (FRDM) (Möller et al. 2016).Each SkyNet run is initialized at a temperature of 6 GK in NSE, and identified by the values of the initial electron fraction   , entropy , and expansion timescale  exp .The latter are considered as initial parameters and later evolved consistently by the network.More details about these nucleosynthesis calculations can be found in Perego et al. (2022).In particular, the heating rates we employ are computed over a comprehensive grid of ∼ 11700 distinct trajectories with 0.01 ≤   ≤ 0.48 linearly spaced, 1.5  B baryon −1 ≤  ≤ 200  B baryon −1 and 0.5 ms ≤  exp ≤ 200 ms logarithmically spaced.These intervals are expected to bracket the properties of the ejecta from BNS and NSBH mergers.We fit the heating rate trajectories obtained with SkyNet over the time interval 0.1 days ≤  ≤ 50 days after merger, using the following power-law dependence: where  1d and  are fit parameters, with typical values  ∼ 1.3 and  1d ∼ 10 10 erg s −1 g −1 .Such a temporal dependence in the heating rate is expected from the decay of large sample of unstable nuclei (Metzger et al. 2010;Korobkin et al. 2012).Moreover, it is equivalent to Eq. ( 30), one of the functional forms used in the optically thick kilonova model described in Sec. 2, provided a conversion factor between the fit reference time (1 day) and  0 .The quality of each single fit is evaluated using a mean fractional log error as employed in Lippuner & Roberts (2015), defined as: where   () is the original SkyNet heating rate trajectory, while the mean is performed over the fit time window without weighing over the time steps, in order to account for the original SkyNet resolution.
For most trajectories we find the average relative errors to be smaller than ∼ 1%.The largest errors are found at the boundary of the SkyNet grid, where the relative error can be as large as ∼ 5%.In Fig. 2, the values of the fitting coefficients are plotted against   ,  and  exp for representative sections of the SkyNet grid.As shown in the left column, for a fixed   the fit parameters are generally smooth functions of the two other thermodynamic variables, and in particular the value of  remains roughly constant (for   ≲ 0.2 it hardly deviates from ∼ 1.3, as already found in Korobkin et al. (2012)), while the value of  0 varies within a factor of a few.On the other hand, the variability of the fit parameters increases as the electron fraction is left free to vary.This strong and non-trivial dependence of the heating rate on the electron fraction is more evident for high   values, where the radioactive heating can be dominated by the decay of individual nuclear species, depending on the specific ejecta conditions.However, we find that the continuity in the fit parameters endures at least in the region which is more relevant to our study, i.e. for   ≤ 0.36,  ≤ 90  B baryon −1 and  exp ≤ 30 ms.We therefore adopt a trilinear interpolation of the fitting coefficients as functions of   , , and  exp in that region, while isolated points or boundary areas for which the continuity of the fitting coefficients is poor are treated by using a nearest-point interpolation.
In order to account for the efficiency with which decay products thermalize in the ejecta, we apply a thermalization efficiency factor to the heating rate as follows.For the thick core of the ejecta, we consider both a constant thermalization efficiency (compatible with all the analytic solutions presented in Sec. 2) and a thermalization efficiency with a time evolution  th =  th,0 (/ 0 ) − , as described in Sec.2.1.2,and Sec.2.1.4.The latter formula approximately mimics the decreasing in the thermalization behaviour expected during the first day in the optically thick ejecta.The values of  th,0 and  th can be fixed by imposing, for example, a thermalization efficiency of ∼ 0.7 and 0.4 at 0.1 days and 1 day, respectively.For the thin layers of the ejecta instead, we model a thermalization efficiency profile starting from the analytic formula proposed in Barnes et al. (2016) and fitted on the properties of the ejecta: where ,  and  are the fit parameters.In that work, this expression was obtained by assuming Eq. ( 7), and  (, ) = .Here instead, we adopt Eq. ( 25), and  (, ) =  1 −  2 −1 .We interpolate the fit parameters in Eq. ( 59) on the tabulated grid reported in Barnes et al. (2016), which spans the intervals 1 × 10 −3  ⊙ <  ej < 5 × 10 −2  ⊙ for the total ejecta mass and 0.1  <  ej < 0.3  for its average velocity.This combination of different efficiencies is motivated by the fact that, on one hand, we expect the decay energy in the thick bulk to thermalize in a similar way as long as the density is sufficiently high.In particular, roughly ∼ 35% of the energy escapes in the form of neutrinos, ∼ 45% is constituted by -rays which efficiently heat the material only within the first day post-merger, and the remaining ∼ 20% is carried by -particles, -particles and fission yields (Barnes et al. 2016).On the other hand, the thermalization efficiency drops significantly in the outer layers of the ejecta, where the lower density makes it harder for the decay products to deposit their energy through thermal processes.Fig. 3 shows the modeled thermalization efficiency profile for different times after merger.By construction, the efficiency in the thin ejecta rapidly declines to values < 20% after a few days post merger.Concurrently, the photosphere radius receeds inward until it disappears.Despite being an artifact, the discontinuity in the efficiency profile at the photosphere radius is not inconsistent with our photosphere model, which assumes a sharp difference in the properties of matter between the thick and the thin ejecta regions.However, we acknowledge the crudeness of the overall thermalization treatment, which does not rigorously account for the dependency on the ejecta conditions of the specific deposition processes involved.Therefore, we leave for a future investigation the    25) is shown as reference.
impact of more detailed thermalization descriptions on the resulting kilonova light curves.

Opacities
In our framework, we can consider the opacity for the r-process material in the ejecta as a free parameter of the model.Alternatively, we can also provide composition-dependent opacity values following the work of Tanaka et al. (2020), in which systematic atomic structure calculations on each element between Fe ( = 26) and Ra ( = 88) are performed using the integrated code HULLAC ( Bar-Shalom et al. 2001).That study mainly focuses on the ejecta conditions around 1 day after the merger, where the temperature is low enough ( ≲ 20000 K) to find the heavy elements ionization stages typically between I-IV.At this time, the density is assumed to be  = 1 × 10 −13 g cm −3 (which is a typical value for an ejecta with mass  ej ∼ 0.01  ⊙ and velocity  ej ∼ 0.1 ), and from here on the opacity in the IR, optical and UV is dominated by bound-bound transitions (Kasen et al. 2013).Bound-bound opacities on a fixed wavelength grid for the homologously expanding material are computed using the widely employed expansion opacity formalism: where the sum runs over all the transitions in the RT simulations within the bin Δ.Planck mean opacities are then computed for a representative ejecta model with different mixtures of heavy elements, characterized by the value of the initial electron fraction   (see Sec. 4.1).2020) for temperatures 5000 K <  < 10000 K and densities  ∼ 10 −13 g cm −3 and for different ejecta compositions characterized by a specific value of electron fraction.The mass fraction of Lanthanides and Actinides is reported for every considered composition.
In Fig. 4 we report the grey opacity values derived by Tanaka et al. (2020) for ejecta temperatures of 5000 K <  < 10000 K, whereas a stronger temperature dependence is found for  < 5000 K.We note that in more recent works (see, e.g., Banerjee et al. 2023Banerjee et al. , 2022Banerjee et al. , 2020) ) such opacity calculations are extended to the ionization stages V-XI of the elements up to Ra, which are expected to be present for ejecta temperatures up to ∼ 10 5 K at times shorter than 1 day post-merger.We therefore leave the corresponding suggested grey opacities as a possible alternative to the Tanaka et al. (2020) dataset.In general, around 1 to a few days, if the electron fraction is low enough (  ≲ 0.25), the grey opacity is dominated by lanthanides and actinides, with values  ≳ 10 cm 2 g −1 .Instead, an increase in the electron fraction between 0.25 ≲   ≲ 0.35 causes a general decrease of the opacity to values  ∼ 1 − 10 cm 2 g −1 , as the fraction of  -valence shell elements present in the ejecta decreases, leaving room for the -shell atoms to provide the leading contribution.Finally, at even higher electron fractions   ≳ 0.4, the contribution from Fe-like elements dominates the opacity, which reaches values  ∼ 0.1 − 1 cm 2 g −1 .In this instance, we interpolate the values in Fig. 4 to uniquely determine the ejecta opacity on the basis of the input   .

Radiative transfer code
In order to assess the level of reliability of our model, we set up a comparison between the light curves obtained by our semi-analytical model and the ones obtained by a RT kilonova simulation.For the latter, we refer to Kawaguchi et al. (2018Kawaguchi et al. ( , 2021)), who employ the wavelength-dependent Monte Carlo RT code originally presented in Tanaka & Hotokezaka (2013).For a given density structure and abundance distribution, the code computes the time evolution of the photon spectrum in the UVOIR wavelength range, together with multicolor light curves.Differently from the first 3D version, Kawaguchi et al. (2018Kawaguchi et al. ( , 2021) ) assume the ejecta to be axisymmetric.This allows for an increase in the simulation spatial grid resolution, and for the inclusion of special-relativistic effects in the photon transport.Photon-matter interaction is described by considering elastic scattering off electrons, together with free-free, bound-free and bound-bound transitions.The contribution to the opacity from the latter is computed using the expansion opacity formalism described in Sec.4.2, while the atomic transition line list employed in the code is the one already used in Tanaka et al. (2017); Tanaka et al. (2020).Since these atomic data concern the ionization stages I-III, the code is used only for temperatures up to ∼ 10000 K, below which further ionization stages are subdominant.Nuclear heating rates and elemental abundances are directly imported from the nucleosynthesis calculations of Wanajo et al. (2014), based on the post-processing of Lagrangian tracer particles obtained by a fully general relativistic simulation of a BNS merger with approximate neutrino transport.Each reaction network calculation starts from a representative thermodynamic trajectory with an initial electron fraction in the range   = 0.09 − 0.44.The fraction of thermalized energy is computed using the analytic formulae reported in Barnes et al. (2016) for the different decay products.These formulae depend on the mass and velocity of the ejecta in a similar fashion to Eq. ( 59).
In particular, while the velocity parameter in the thermalization formulae is fixed to  ej = 0.3 , the mass parameter is set starting from the local density and considering a uniform sphere of radius  ej .

Comparison setup
We prepare our comparison by setting the same ejecta properties in both codes.We consider two ejecta configurations, namely a lighter anisotropic dynamical component and a more massive spherically symmetric secular component.This choice is motivated by the general necessity of modelling multiple components of matter ejection, which are required in order to reproduce the color bands of observed kilonovae, as in the case of AT2017gfo (Cowperthwaite et al. 2017;Tanvir et al. 2017;Tanaka et al. 2017).
Regarding the secular component, we assume a total mass of  sec = 2.64 × 10 −2  ⊙ , an average velocity of  rms = 0.06  and constant values for the electron fraction and specific entropy, i.e.   = 0.2 and  = 10 k B baryon −1 .We compute the associated expansion time scale as  exp = / rms ≈ 17 ms.These values are representative of the outcomes of simulations that investigate the evolution of disks emerging as remnants of compact binary mergers and accreting onto the central object.In these simulations, a fraction between ∼ 20 − 40% of the disk mass is expelled during the secular evolution, with the initial disk mass  disk ∼ 10 −4 − 10 −1  ⊙ (see, e.g., Fahlman & Fernández 2022;Fujibayashi et al. 2020;Fernández et al. 2019;Hotokezaka et al. 2013).Instead, for the dynamical component, we use the properties of the dynamical ejecta extracted from one GRHD simulation of a BNS merger with M0 neutrino transport approximation, chosen among those performed by Perego et al. (2019) and compatible with the GW170817 event.Despite the simulations considered in that work include different EOSs, they all lead to similar ejecta angular distributions, and therefore we arbitrarily select the simulation employing the HS(DD2) EOS (Hempel & Schaffner-Bielich 2010).The dynamical ejecta are identified with the matter unbound within the end of the simulation according to the geodesic criterion, i.e. the matter for which |  | ≥ , with   the time-component of the four-velocity.For an equal-mass binary with masses  1,2 = 1.364  ⊙ , a total dynamical ejecta mass  dyn = 2.7 × 10 −3  ⊙ was found.The properties of this component Table 1.Fit parameters boundaries considered in the fitting procedure.The wider intervals adopted for the parameters associated to the secular ejecta reflect the major variability in composition of such component, based on the outcome of different numerical simulations.25), while a radially constant electron fraction is assumed.Angular distributions are extracted from the GRHD simulation of an equal-mass BNS system with masses  1,2 = 1.364  ⊙ using the DD2 EOS (Perego et al. 2019).

Secular wind NR dynamical ejecta
are recorded as matter crosses a extraction spherical surface characterized by a coordinate radius  E = 294 km, and are then reduced to an axisymmetric configuration by averaging over the azimuthal angle.In particular, xkn is informed with the angular distributions of the ejecta mass, average electron fraction and entropy, and average velocity at infinity, calculated as  ∞ rms =  √︁ 1 − (/  ) 2 .We choose the profile given by Eq. ( 25) to describe the radial density structure of each ejecta component both in the RT simulation and in xkn.Moreover, we assume a radially constant electron fraction in order to fix the composition.The resulting configuration of the dynamical ejecta as depicted in Fig. 5 reflects the general characteristics of this component as obtained in many merger simulations: neutron-rich matter is expelled preferentially across the equatorial plane partially through tidal forces, while shock-heated material subject to stronger neutrino irradiation and thus less neutron-rich escapes also at small polar angles.
The input profiles described above uniquely determine all the components of both the models, including energy deposition rates, elemental abundances and opacities, with the only exception of one remaining free parameter in xkn, that is the photospheric floor temperature  floor .We remark that the employed radioactive heating rates, as well as the prescription for the thermalization efficiency, are not coincident between the two models, although derived from the same initial conditions.However, the final energy deposition rates agree within a factor of a few, and we account for this discrepancy as being part of the general difference between the kilonova models.Furthermore, we acknowledge that the opacity treatment in our xkn model is significantly approximated: in addition to the adoption of grey values, we assume the opacity to be constant in time or at most to evolve according to a power-law, when characterizing the ejecta through their entire evolution and depth.In reality we expect the Planck mean opacity to vary by at least one order of magnitude between different regions and epochs.Therefore, the adoption of the opacity values derived by Tanaka et al. (2020) and described in Sec.4.2 is not more physically motivated than treating the opacity as a free parameter, and, for this reason, in the comparison we consider both possibilities.
The RT data employed in the comparison consist of the bolometric luminosity,  RT bol (), and of the AB magnitudes,  RT AB,,  () at different wavelengths , observed from multiple viewing angles  iew ∈ [0 • , 90 • ].We thus fit our free parameters to both sets of data separately, considering a logarithmically spaced time mesh, from 0.5 to 15 days.Within this time frame, we assume that the assumptions of our model are better verified, and that the RT calculations are more reliable, whereas temperatures throughout the ejecta are well below 10000 K, justifying the employed atomic data.We define two error functions in order to establish the fit procedure.For the bolometric luminosity, we compute the absolute logarithmic error between our model luminosity,  M bol , and the RT result,  RT bol , averaged over all   data points in the considered time frame: For the AB magnitudes, we consider three representative broadband filters, namely the  ( = 2157 nm),  ( = 972 nm) and  ( = 475 nm) filters.The light curves are calculated assuming a source luminosity distance of   = 40 Mpc, corresponding to the estimated distance for the merger associated to the AT2017gfo signal.Since our kilonova model is better suited to reproduce the light curve behaviour around the emission peak, only data points such that  RT AB, < 30 are considered, in order to avoid having the fits influenced by too dim values.In a similar fashion to the bolometric luminosity, we compute the absolute error between the magnitudes across the three different wavebands and two different viewing angles, i.e. 0 • and 90 We perform two sets of runs, one for each ejecta configuration, i.e. one for the secular isotropic ejecta and one for the dynamical anisotropic ejecta.Furthermore, for each set, we take into account two possibilities.In one case we allow the opacity in our model to vary freely, and in particular for the anisotropic setup we assume it to follow a step function, i.e. we adopt a higher value,  low lat , at low latitudes ( ≳ 45 • ) in correspondence of a neutron-rich environment with   ≲ 0.25, and a lower value,  high lat , at high latitudes ( ≲ 45 • ) where   ≳ 0.25.In the other case instead we compute the opacity using the   parametrization from Tanaka et al. (2020), leaving us with only the photospheric floor temperature to be fitted.
To be consistent with the opacity prescription, for the anisotropic ejecta setup we consider a   -dependent floor temperature parameterized by the two values  Ni and  LA .In Table 1 we report the adopted ranges for the parameters included in the fit procedure.Finally, each calculation is repeated using the semi-analytic kilonova model presented in Perego et al. (2017) for comparison purposes.The latter shares the same multicomponent, anisotropic framework as the model presented in this work.However, the underlying kilonova model is not based on the solution of the diffusion equation, but it is a phenomenological description based on timescale arguments, presented in Grossman et al. (2014) and Martin et al. (2015).Due to this distinction, we name the previous model as xkn-ts, as opposed to our new xkn-diff model.

Comparison results
The results for all the different models and configurations considered are summarized in Table 2.As visible, the fit procedure returns reasonable fit parameters values falling in the prior intervals, with the exception of a minor number of cases useful to let the modelling limits emerge.In general, both the fits on the bolometric luminosities and the magnitudes derived from the RT simulation show an overall improved fit quality when using xkn-diff with respect to the previous xkn-ts model.
In particular, when fitting on the bolometric luminosity, xkn-ts is limited by having only the degree of freedom associated to the ejecta opacity (when the latter is left free to vary), since the temperature floor does not enter the luminosity calculation, as opposed to the xkn-diff case.This difference arises because the floor temperature affects the late time photospheric radius in both models, but while in xkn-ts the latter is used only for the magnitudes computation through the Stefan-Boltzmann law and does not modify the volume of the radiative zone, in the xkn-diff case the photosphere position has a feedback on the allocation of mass to the optically thick and optically thin regimes, thus altering the bolometric luminosity as well.As a result, for the case in which the opacity is prescribed using the value of the electron fraction, the bolometric luminosity in xkn-ts is completely fixed for both the secular and the dynamical component configurations, and the correspondent fit errors are the worst in the set.On the other hand, in the xkn-diff model the floor temperature is adjusted in such a way to increase the late time agreement.Indeed such parameter can ultimately affect light curves only when temperatures in the ejecta have decreased enough, as it is commonly assumed to be the case around a few days post-merger.Specifically, in the secular component configuration this dependency drives the floor temperature to almost rail against the upper boundary, in order to accelerate the photosphere recession and maximize the amount of thin ejecta contributing to the late time emission.One can note that in all cases, and more rapidly for the faster dynamical component, both semi-analytic models converge to the same curve, since in the xkn-diff model the treatment of the thin ejecta, which eventually constitutes the totality of the outflow, is analytically equivalent to the one used for the entire ejecta in xkn-ts.Furthermore, since this treatment does not properly model the material opacity outside of the photosphere, varying the latter cannot affect the computed luminosity around 10 days, thus not improving the late time matching.As visible in Fig. 6, in all the investigated configurations xkn-ts is sistematically underestimating the early time luminosity with respect to both the RT simulation and xkn-diff of a factor from a few to even multiple orders of magnitude depending on the specific case.This evidence establishes qualitatively the error hierarchy in using an approximate scheme based on the calculation of the diffusion timescale of photons, versus an analytic solution of the simplified RT problem, with respect to a full RT simulation.Therefore, once the opacity is left free to vary, the physiologic behaviour of xkn-ts is to compensate this systematic by lowering the latter to very small values in order to increase the emission brightness especially before ∼ 1 day, even incurring in the opacity boundaries in the dynamical component case.Also the xkn-diff model is partially subject to the same mechanism, as visible specifically in the secular component case.This result suggests a general limitation on using the bolometric luminosity to fit parameters related to local features of the ejecta configuration.However, we also note that, especially for the dynamical component, the magnitude and shape of the bolometric luminosity are in good agreement with the ones derived from the RT calculations.
With respect to the fit on the bolometric luminosity, when the same procedure is applied to the AB magnitudes, the overall qualitative results remain roughly unaltered.However, in such a case we include by construction more information, coming from different wavebands and viewing angles.In addition, the temperature floor has a direct role in determining the color bands, since they strongly depend on the photosphere effective temperature, and thus this parameter influences the fit outcome regardless of the model employed.Therefore, we obtain best fit values not necessarily close to the ones found in the previous case.Both in the secular and the dynamical component configurations, we find the same error hierarchy as in the bolometric luminosity fits, with xkn-ts tipically underestimating the overall brightness up to a few days in all bands with respect to xkn-diff.Furthermore, Fig. 7 shows how magnitudes confirm the trend already evident for the bolometric luminosity, by which xkn-diff tends to underestimate the emission brightness at early times ≲ 1 day, indicating a model limitation.The retrieved opacities are now generally higher for both models, with values which can also be closer to the ones fixed by the atomic calculations.The fact that the opacity values differ significantly from the previous fits is not surprising, since in this case they are informed with the light curves in multiple filters as seen in edge-on and face-on configuration: especially in the dynamical ejecta configuration, the latter is valuable information in determining the opacity angular distribution with a better accuracy.In addition, we recall that color bands in the model are derived by composing a sprectrum mainly based on pure black-body emission, which is therefore not able to reproduce the black-body deviations found in the RT calculations.In particular, as pointed out in Gillanders et al. (2022), we note that Table 2. Best fit parameters obtained in the case of the bolometric luminosity (top) and magnitudes (bottom) with corresponding fit errors for the dynamical (right) and the secular (left) component configurations, as obtained by using xkn-diff and xkn-ts models, with and without a fixed opacity.Cases for which the fit rails against the chosen boundaries are highlighted in gray.(*) is used to indicate that the parameter value is fixed, either manually for the secular ejecta or from the NR simulation for the dynamical ejecta.(-) is placed when the value does not affect the fit outcome.m AB (θ view = 45

Bolometric luminosity
• ) Figure 7. AB magnitudes in the ,  and  filters as seen from a viewing angle  iew = 45 • , calculated for the best fit parameters in the dynamical and the secular component configurations, using xkn-diff and xkn-ts models, with and without a fixed opacity.Curves are shown in comparison to the data on which the fits were performed, derived from the RT calculations obtained with the code developed in Tanaka & Hotokezaka (2013); Kawaguchi et al. (2018Kawaguchi et al. ( , 2021)).
realistically part of the UV radiation is reprocessed by the heavy elements into the optical and NIR bands, thus shifting the emitted energy distribution significantly, without an heavy alteration in the bolometric luminosity.As a consequence, the more sensitive fits on the magnitudes retrieve opacity values which are increased in order to compensate for the lack of such feature in the model.The floor temperatures derived from the magnitudes are substantially different from the ones derived from the bolometric luminosity.This is due a combination of effects, whereby the floor temperature is not trivially connected to the final magnitudes and its value is subject to stronger variability.
On one side, higher floor temperatures are associated to stronger radiation fluxes at late times and, for a given energy emission rate, to smaller photospheric radii, with a net increase in the late broadband magnitudes.Being this the only effect in xkn-ts, the magnitude fit finds the best temperature floor parameter value up to ∼ 3400 K, in order to compensate for the systematic underestimation of the model, partially relieving the opacity parameter from such burden.This behaviour has also to be ascribed to our fit procedure, which tries to minimize quantitatively the separation between different curves, rather than trying to reproduce the same shape.For this reason, the detailed values of floor temperature that we obtain are not meant to be reliable, but they can nevertheless highlight the internal structure of the model.On the other side, on top of the above effect, as already pointed out for the bolometric luminosity fits, in the xkn-diff model higher floor temperatures also cause a faster decrease in the amount of optically thick ejecta at late times.In particular, in this case the temperature floor recovery which results from such interplay cures the drift towards the upper boundary that is found in the bolometric luminosity fit for the secular ejecta configuration with fixed opacity.As a general consequence, for xkn-diff, values are almost systematically and significantly lower than both xkn-ts and their counterparts in the previous fits, being in some cases down to only a few hundreds Kelvin degrees, and indicating a tendency to decrease the overall radiation fluxes in order to match color bands after a few days.

CONCLUSIONS
In this work we presented the new framework xkn for the computation of the kilonova emission from compact binary mergers, starting from the characterization of the merger ejecta.The framework allows for non-trivial ejecta structures as can be inferred from numerical relativity merger simulations, and it employs the results of recent efforts in nuclear astrophysics and atomic physics in terms of inputs for the kilonova model, such as radioactive heating rates from nuclear reaction network calculations (Perego et al. 2022) and grey opacities from systematic atomic structure calculations (Tanaka et al. 2020).With respect to previous iterations, xkn includes the model xkn-diff, which encapsulate different possible semi-analytic solutions of the diffusion equation for the radiation energy density field, derived from the RT problem under the assumption of homologously expanding material (Wollaeger et al. 2018).xkn-diff constitutes an improvement with respect to previous semi-analytic models based on simpler laws of energy conservations and approximate radiation diffusion timescale estimations.In addition, the model tracks the position of the ejecta photosphere in time in order to distinguish between the optically thick internal bulk and the optically thin external layers.The latter are treated with a simplified shell model, which approximately accounts for the non-negligible contribution to the total luminosity coming from this region from a few days post-merger on.We tested xkn models by comparing their results with the ones obtained from two-dimensional RT simulations obtained with the code developed in Tanaka & Hotokezaka (2013); Kawaguchi et al. (2018Kawaguchi et al. ( , 2021)), which are based on the same ejecta configurations.In particular, we considered two representative scenarios, i.e. a lighter anisotropic dynamical component and a more massive spherical secular component, and we fit the free parameters of the model to the bolometric luminosity and the AB magnitudes in different filters, as seen from multiple viewing angles.We found that xkn-diff is able to reproduce the overall behaviour of the light curves obtained from the RT simulations, with a better agreement with respect to the previous semi-analytic model, despite the simplified treatment of the decay energy thermalization process and of the ejecta opacity.However, as highlighted by the fit procedure, the latter still constitutes a limitation to this modelling approach and it will be the subject of future improvements.In particular, the average constant grey opacity values that the model employs are a crude approximation of the real effective opacity inside the ejecta, which significantly varies of more than one order of magnitude with time and across the different regions of the outflow, depending on the local temperature, density and composition.As a result, the emission brightness at early times, i.e. around a few hours post-merger, predicted by the model in the fit procedure, can be systematically lower with respect to the RT calculation, of a factor of a few in the bolometric luminosity and of up to 2 magnitudes in the color bands.We also note that the temperature floor, a secondary parameter in xkn which often appears in other semi-analytic models, is not easily constrained, since it is not trivially connected to the final magnitudes.We conclude that xkn constitutes a valid tool to model the kilonova emission from compact binary mergers, with the main strength being its computational efficiency, which allows for extensive explorations of the ejecta parameter space in a reasonable time frame.This is particularly useful in the context of the now thriving multi-messenger astronomy, whereas the kilonova is only one of the possible electromagnetic counterparts of the merger event.Coupling this model with information from other sources, such as the GRB afterglow or the GW signal, in a statistical framework, can sinergically help to constrain the properties of the original binary, the central remnant or the merger ejecta, and thus to shed light on the nature of the detected event itself.

Figure 2 .
Figure 2. Heating rate fit parameters as functions of initial electron fraction   and entropy , for fixed values of the expansion timescale  exp .The dashed line separates the region where we interpolate the parameters linearly (left) from the region where the continuity of the parameters is poor and we use a nearest-point interpolation (right).

Figure 3 .
Figure 3. Constructed thermalization efficiency radial profile at different days post-merger.The ejecta density radial profile from Eq. (25) is shown as reference.

Figure 4 .
Figure 4. Grey opacity values derived byTanaka et al. (2020) for temperatures 5000 K <  < 10000 K and densities  ∼ 10 −13 g cm −3 and for different ejecta compositions characterized by a specific value of electron fraction.The mass fraction of Lanthanides and Actinides is reported for every considered composition.

Figure 5 .
Figure 5. Density and electron fraction distributions the velocity space for the dynamical ejecta at  = 100 s post-merger.The density radial profile follows the analytic prescription of Eq. (25), while a radially constant electron fraction is assumed.Angular distributions are extracted from the GRHD simulation of an equal-mass BNS system with masses  1,2 = 1.364  ⊙ using the DD2 EOS(Perego et al. 2019).

Figure 6 .
Figure6.Bolometric luminosity correspondent to the best fit parameters in the dynamical and the secular component configurations, as obtained by using xkn-diff and xkn-ts models, with and without a fixed opacity.Curves are shown in comparison to the data on which the fits are performed, derived from the RT calculations obtained with the code developed inTanaka & Hotokezaka (2013);Kawaguchi et al. (2018Kawaguchi et al. ( , 2021)).