Nebular Emission from Lanthanide-rich Ejecta of Neutron Star Merger
Kenta Hotokezaka, Masaomi Tanaka, Daiji Kato, Gediminas Gaigalas
MMNRAS , 1–15 (2020) Preprint 17 February 2021 Compiled using MNRAS L A TEX style file v3.0
Nebular Emission from Lanthanide-rich Ejecta of Neutron StarMerger
Kenta Hotokezaka, , ★ Masaomi Tanaka, Daiji Kato, , Gediminas Gaigalas Research Center for the Early Universe, Graduate School of Science, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japanu University, Aoba, Sendai 980-8578, Japan Astronomical Institute, Tohoku University, Aoba, Sendai 980-8578, Japan National Institute for Fusion Science, 322-6 Oroshi-cho, Toki 509-5292, Japan Department of Advanced Energy Engineering Science, Kyushu University, Kasuga, Fukuoka 816-8580, Japan Institute of Theoretical Physics and Astronomy, Vilnius University, Saul˙etekio Ave. 3, Vilnius, Lithuania
17 February 2021
ABSTRACT
The nebular phase of lanthanide-rich ejecta of a neutron star merger (NSM) is studied byusing a one-zone model, in which the atomic properties are represented by a single species,neodymium (Nd). Under the assumption that 𝛽 -decay of 𝑟 -process nuclei is the heat andionization source, we solve the ionization and thermal balance of the ejecta under non-localthermodynamic equilibrium. The atomic data including energy levels, radiative transitionrates, collision strengths, and recombination rate coefficients, are obtained by using atomicstructure codes, GRASP2K and
HULLAC . We find that both permitted and forbidden lines roughlyequally contribute to the cooling rate of Nd II and Nd III at the nebular temperatures. We showthat the kinetic temperature and ionization degree increase with time in the early stage ofthe nebular phase while these quantities become approximately independent of time after thethermalization break of the heating rate because the processes relevant to the ionization andthermalization balance are attributed to two-body collision between electrons and ions at latertimes. As a result, in spite of the rapid decline of the luminosity, the shape of the emergentspectrum does not change significantly with time after the break. We show that the emission-line nebular spectrum of the pure Nd ejecta consists of a broad structure from 0 . 𝜇 m to 20 𝜇 mwith two distinct peaks around 1 𝜇 m and 10 𝜇 m. Key words: transients: neutron star mergers
Neutron star mergers (NSMs) have been considered as the sitesof 𝑟 -process nucleosynthesis (Lattimer & Schramm 1974). In Au-gust 2017, the LIGO/Virgo Collaboration (LVC) discovered the firstNSM, GW170817, which was accompanied by radiation across theentire electromagnetic spectrum (Abbott et al. 2017; Nakar 2020;Margutti & Chornock 2020). In particular, the spectrum and lightcurve of the uv-optical-infrared counterpart referred to as ‘kilonova’or ‘macronova’ indicate that a copious amount of 𝑟 -process elementsis produced in this event (Andreoni et al. 2017; Arcavi et al. 2017;Coulter et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017;Evans et al. 2017; Kasliwal et al. 2017; Pian et al. 2017; Smarttet al. 2017; Tanvir et al. 2017; Utsumi et al. 2017). The amount ofthe produced 𝑟 -process elements and the event rate estimated fromGW170817 suggest that NSMs could provide all the 𝑟 -process el- ★ E-mail: [email protected] ements in the Galaxy (e.g. Hotokezaka et al. 2018; Rosswog et al.2018).Lanthanide ions have unique optical properties, which enhancethe opacity of the NSM ejecta material (Barnes & Kasen 2013;Kasen et al. 2013; Tanaka & Hotokezaka 2013; Wollaeger et al.2018; Bulla 2019; Barnes et al. 2020). Thus, the existence of lan-thanide ions imprints observable signatures in kilonova light curvesand spectra. In fact, the late-time spectra of GW170817 peakingaround the near infrared (nIR) band implies that lanthanides existin the GW170817 ejecta (Kasen et al. 2017; Tanaka et al. 2017b).At the same time, the light curve rises on a short time scale of ∼ . ∼ © a r X i v : . [ a s t r o - ph . H E ] F e b K. Hotokezaka et al. day consist of a single temperature blackbody with structure pro-duced by atomic transitions. They found that the main structure ofthe spectra is consistent with the P Cygni profiles produced by Sr IIdoublet 4078 , , , 𝛾 -rays canbe one of the most robust identifications of radioactive isotopes(Hotokezaka et al. 2016; Li 2019; Wu et al. 2019b; Korobkin et al.2020). However, such measurements are very challenging and onlyweak upper limits were put by NuSTAR in GW170817 (Evans et al.2017).Here we consider the nebular phase of kilonovae, where theemergent spectrum is dominated by emission lines, and hence, spec-troscopic observations may enable to identify the elements producedin NSMs. Since the slower ejecta component can be observed in thenebular phase than the earlier phases one can expect that the Dopplerbroadening of lines is weaker so that the spectral structure arisingfrom individual lines may be more pronounced. In GW170817, the
Spitzer
Space Telescope detected the late-time nebular emission at4 . 𝜇 m and put upper limits at 3 . 𝜇 m (Kasliwal et al. 2019; Vil-lar et al. 2018), suggesting that a fraction of the luminosity of thenebula is radiated in infrared with a peculiar spectral shape.The primary goal of this paper is to address the evolutionof thermodynamic quantities and the emerging spectral shape oflanthanide-rich NSM nebulae. The early works on the nebular mod-elings of kilonovae assume local thermodynamic equilibrium (LTE)for ionization and level population (Waxman et al. 2018; Gillanderset al. 2021). However, the non-LTE effects are crucial for the late-time nebular modelings. Here we develop a NSM nebula modelunder non-LTE by following the studies of supernova (SN) neb-ular emission (Axelrod 1980; Fransson & Chevalier 1989; Ruiz-Lapuente & Lucy 1992; Mazzali et al. 2006; Maeda et al. 2006;Botyánszki et al. 2018). The paper is organized as follows. In §2,we describe the heating and ionization rates due to 𝛽 -decay of 𝑟 -process nuclei. In §3, we describe the equations and several ap-proximations used in the modeling. In §4, we show the atomic dataobtained by using the atomic codes. In §5, we apply our model toa lanthanide-rich NSM nebula and show the time evolution of tem-perature, ion abundances, and emission spectra. We conclude anddiscuss our study in §6. Calculating the nebular emission generally requires radiative trans-fer computations under non-LTE. However, such computations forNSM nebulae demand a lot of effort. As a first step, we focus hereon the nebular phase where the following conditions are satisfied:(i) the ejecta is optically thin and (ii) the recombination and coolingtimes are shorter than the dynamical time. The former allows us tosimplify the treatment of radiative transfer and the latter allows touse the steady-state approximation. Time since merger [day] ˙ Q t h [ e r g s − g − ] . M (cid:12) , . c . M (cid:12) , . c . M (cid:12) , . c SNIa , Ni , Co Time since merger [day] − − − − − Γ / n [ e r g c m s − ] . M (cid:12) , . c . M (cid:12) , . c . M (cid:12) , . c SNIa , Ni , Co Figure 1.
Specific heating rates and normalized heating rates of 𝛽 -decay of 𝑟 -process nuclei. The final composition is assumed to be the solar 𝑟 -processabundance pattern in a mass range of 130 (cid:54) 𝐴 (cid:54) Ni → Co → Fedenoted as SN Ia. For calculating the thermalization in SNe Ia efficiencies,we assume 𝑀 Ni = . 𝑀 (cid:12) , 𝑀 ej = . 𝑀 (cid:12) , and 𝑣 = − . Herethe heating rates are computed by using an open code (Hotokezaka & Nakar2020). The optical depth of the NSM ejecta is estimated by 𝜏 ≈ 𝜅𝑀 ej 𝜋𝑣 𝑡 , (1)where 𝑀 ej is the ejecta mass and 𝜅 is the ejecta opacity. Here,we assume a homologous expansion of the ejecta, i.e., 𝜌 ( 𝑡, 𝑣 ) ∝ 𝑡 − 𝑣 − 𝛼 between 𝑣 and 𝑣 , where 𝜌 is the ejecta density, 𝑡 is timesince merger, 𝑣 and 𝑣 are the minimum and maximum expansionvelocities, and 𝛼 describes the velocity profile of the ejecta. In thispaper, we consider spherical symmetric ejecta for simplicity.In the case that the ejecta material is mainly composed of 𝑟 -process elements, 𝜅 is typically ≈
10 cm g − for lanthanide-richmaterial and ≈ g − for lanthanide-free material (Barnes &Kasen 2013; Kasen et al. 2013; Tanaka & Hotokezaka 2013; Wol-laeger et al. 2018; Tanaka et al. 2020). The time when a kilonovaenters the NSM nebular phase is roughly estimated by 𝜏 ( 𝑡 thin ) ≈ 𝑡 thin ≈ √︄ 𝜅𝑀 ej 𝜋𝑣 , (2) ≈
35 day (cid:18) 𝜅
10 cm g − (cid:19) / (cid:18) 𝑀 ej . 𝑀 (cid:12) (cid:19) / (cid:16) 𝑣 . 𝑐 (cid:17) − , (3)where 𝑐 is the speed of light. Note that this time scale is significantlyshorter for lanthanide-free material and fast expanding material. Forinstance, 𝑡 thin is ≈ 𝑀 ej = . 𝑀 (cid:12) , 𝜅 = . g − , and 𝑣 = . 𝑐 .We assume the radioactivity of 𝑟 -process nuclei as the source MNRAS , 1–15 (2020) eutron Star Merger Nebula of heat and ionization. Note that the light curve of the GW170817kilonova is consistent with the picture that 𝛽 -decay of 𝑟 -processnuclei predominantly heats the ejecta material over the time scalesfrom 0 . . In the nebular phase,the ejecta is optically thin for 𝛾 -rays, and hence, we consider only 𝛽 -decay electrons. We use the heating rate with 130 (cid:54) 𝐴 (cid:54) 𝐴 is atomic mass number, provided by Hotokezaka & Nakar(2020). In the following, we describe the characteristic features ofthe 𝛽 -decay heating rate relevant to the NSM nebular modelings.At early times, the heating rate per unit mass approximatelyfollows (Metzger et al. 2010): (cid:164) 𝑄 th ( 𝑡 ) ∝ 𝑡 − . (4)This power law is valid as long as thermalization of 𝛽 -decay elec-trons occurs on a time scale much shorter than a dynamical time.The heating rate starts to deviate from equation (4) around the ther-malization time, 𝑡 th , estimated by 𝜏 eff , e = 𝜅 eff , e 𝜌𝑐𝑡 ∼ 𝑡 th ≈ (cid:32) 𝐶 𝜌 𝑐𝜅 𝑒, eff 𝑀 ej 𝑣 (cid:33) / , ≈
55 day (cid:18) 𝐶 𝜌 . (cid:19) / (cid:18) 𝑀 ej . 𝑀 (cid:12) (cid:19) / (cid:16) 𝑣 . 𝑐 (cid:17) − / × (cid:18) 𝜅 eff ,𝑒 . g − (cid:19) / (cid:18) 𝐸 𝑒 .
25 MeV (cid:19) − / , (5)where 𝐸 𝑒 is the initial energy of 𝛽 -decay electrons and 𝜅 eff ,𝑒 is aneffective opacity of the interaction of 𝛽 -decay electrons with theejecta material. Hereafter, we use the one-zone approximation, inwhich the density at a given time is represented by the mass weightedmean, 𝜌 𝑚 ( 𝑡 ) = 𝐶 𝜌 𝑀 ej 𝑡 − 𝑣 − , where 𝐶 𝜌 is a normalization constant(Hotokezaka & Nakar 2020).For 𝑡 (cid:38) 𝑡 th , the specific heating rate declines as (Kasen &Barnes 2019; Waxman et al. 2019; Hotokezaka & Nakar 2020) (cid:164) 𝑄 th ( 𝑡 ) ∝ 𝑡 − . ( for 𝑡 (cid:29) 𝑡 th ) . (6)This break in the 𝛽 -decay heating rate from ∝ 𝑡 − . to ∝ 𝑡 − . isreferred to as the thermalization break, which typically occurs inthe nebular phase (see equations 3 and 5).It is useful to introduce the normalized heating rate: Γ 𝑛 = 𝜌 𝑚 (cid:164) 𝑄 th 𝑛 ∝ (cid:26) 𝑡 . ( 𝑡 (cid:28) 𝑡 th ) ,𝑡 . ( 𝑡 (cid:29) 𝑡 th ) . (7)where Γ is the heating rate per unit volume and 𝑛 is the mass-weighted mean atomic number density 𝑛 ≈ · cm − (cid:18) (cid:104) 𝐴 (cid:105) (cid:19) − (cid:18) 𝐶 𝜌 . (cid:19) × (cid:16) 𝑣 . 𝑐 (cid:17) − (cid:18) 𝑀 ej . 𝑀 (cid:12) (cid:19) (cid:16) 𝑡
35 d (cid:17) − , (8)where (cid:104) 𝐴 (cid:105) is the mean atomic mass of the ejecta material. As wewill see later, the evolution of kinetic temperature and ionizationdegree roughly follows the evolution of the normalized heating rate.Figure 1 shows the specific heating rates, (cid:164) 𝑄 th , and normalizedheating rates, Γ / 𝑛 , with three different combinations of the ejectamass and velocity. For more massive and slower ejecta, the normal-ized heating rate at a given time is smaller, corresponding to that the Spontaneous fission and 𝛼 -decay of heavy nuclei can also be the energysource (Zhu et al. 2018; Wanajo 2018; Wu et al. 2019a) efficiency of ionization and heating is lower. As expected from equa-tion (7), the slope of the normalized heating rates becomes almostflat at later times. For comparison, figure 1 also shows the heatingrate of the decay chain powering SNe Ia, Ni → Co → Fe, with 𝑀 Ni = . 𝑀 (cid:12) , 𝑀 ej = . 𝑀 (cid:12) and 𝑣 = − . Unlike the 𝑟 -process cases the normalized heating rate of this decay chain turnsto decrease around the half-life of Co. Note that the normalizedheating rate of NSMs is much larger than that of SNe Ia because ofthe difference in the expansion velocity, suggesting that ionizationin the NSM ejecta is more efficient.The ionization rate of an 𝑖 -th ionized ion, 𝑋 𝑖 + , by 𝛽 -decayelectrons is characterized by the work per ion pair 𝑤 𝑖 (see AppendixA). With this quantity, the ionization rate per unit volume is givenby Υ 𝑖 = Γ 𝑤 𝑖 . (9)For the NSM nebulae, we estimate 𝑤 𝑖 / 𝐼 𝑖, ∼
30, where 𝐼 𝑖, isthe first ionization potential of 𝑋 𝑖 + . This value indicates that thesignificant fraction of 𝛽 -electrons’ energy is deposited to the thermalenergy and only ∼ In the nebular phase, the ejecta material is in non-LTE, i.e., onlyfree electrons are distributed according to Maxwell’ law and atomsare not in equilibrium. Thus, one must solve the ionization andthermal balance to obtain the kinetic temperature, 𝑇 𝑒 , and ionizationfractions. Here we use the nebular modeling developed by Axelrod(1980) with some modifications. In this section, we briefly describethe equations used and discuss some generic features of the NSMnebular emission that arise from the characteristic properties of the 𝑟 -process heating rate (equation 7) without specifying the details ofthe atomic structure.We consider the NSM nebular phase where the recombinationand cooling time scales are shorter than a dynamical time. Thiscondition allows us to use the steady state approximation. As wewill show later, it holds 𝑡 (cid:46)
100 day after merger. We assume thatthe ejecta is composed of a single atomic species, 𝑋 , for simplicity .Under these conditions, the equation for ionization balance is − Υ 𝑖 𝑓 𝑖 − ∑︁ 𝑗>𝑖 𝑃 𝑖 𝑗 𝛼 𝑗 + 𝜒 𝑓 𝑗 + 𝑛 + ( − 𝑃 𝑖𝑖 ) 𝛼 𝑖 + 𝜒 𝑓 𝑖 + 𝑛 ≈ (cid:54) 𝑖 (cid:54) 𝑁 − , (10)where 𝑓 𝑖 is the number fraction of 𝑋 𝑖 + , 𝑃 𝑖 𝑗 is the probability thatthe photons created by the recombination of an ion 𝑋 𝑗 + ionize anion 𝑋 𝑖 + , 𝛼 𝑖 is the recombination rate coefficient for 𝑋 ( 𝑖 + )+ → 𝑋 𝑖 + ,and 𝜒 is the free electron fraction. This equation can be rewrittenin the form (Axelrod 1980) − Γ 𝑛 𝑓 𝑖 𝑤 𝑖 − ∑︁ 𝑗>𝑖 𝑃 𝑖 𝑗 𝛼 𝑗 + 𝜒 𝑓 𝑗 + + ( − 𝑃 𝑖𝑖 ) 𝛼 𝑖 + 𝜒 𝑓 𝑖 + ≈ (cid:54) 𝑖 (cid:54) 𝑁 − , (11)The ion fractions 𝑓 𝑖 and free electron fraction 𝜒 are obtained bysolving equation (11) together with the conditions of (cid:205) 𝑍𝑖 = 𝑓 𝑖 = We consider a single atomic species only for ionization and thermalbalances. However we consider that the 𝛽 -decay heat is produced by manydifferent isotopes.MNRAS000
100 day after merger. We assume thatthe ejecta is composed of a single atomic species, 𝑋 , for simplicity .Under these conditions, the equation for ionization balance is − Υ 𝑖 𝑓 𝑖 − ∑︁ 𝑗>𝑖 𝑃 𝑖 𝑗 𝛼 𝑗 + 𝜒 𝑓 𝑗 + 𝑛 + ( − 𝑃 𝑖𝑖 ) 𝛼 𝑖 + 𝜒 𝑓 𝑖 + 𝑛 ≈ (cid:54) 𝑖 (cid:54) 𝑁 − , (10)where 𝑓 𝑖 is the number fraction of 𝑋 𝑖 + , 𝑃 𝑖 𝑗 is the probability thatthe photons created by the recombination of an ion 𝑋 𝑗 + ionize anion 𝑋 𝑖 + , 𝛼 𝑖 is the recombination rate coefficient for 𝑋 ( 𝑖 + )+ → 𝑋 𝑖 + ,and 𝜒 is the free electron fraction. This equation can be rewrittenin the form (Axelrod 1980) − Γ 𝑛 𝑓 𝑖 𝑤 𝑖 − ∑︁ 𝑗>𝑖 𝑃 𝑖 𝑗 𝛼 𝑗 + 𝜒 𝑓 𝑗 + + ( − 𝑃 𝑖𝑖 ) 𝛼 𝑖 + 𝜒 𝑓 𝑖 + ≈ (cid:54) 𝑖 (cid:54) 𝑁 − , (11)The ion fractions 𝑓 𝑖 and free electron fraction 𝜒 are obtained bysolving equation (11) together with the conditions of (cid:205) 𝑍𝑖 = 𝑓 𝑖 = We consider a single atomic species only for ionization and thermalbalances. However we consider that the 𝛽 -decay heat is produced by manydifferent isotopes.MNRAS000 , 1–15 (2020) K. Hotokezaka et al. and 𝜒 = (cid:205) 𝑍𝑖 = 𝑖 𝑓 𝑖 for a given 𝑛 and 𝑇 𝑒 . The details of the reprocessof recombination radiation, 𝑃 𝑖 𝑗 , are described in Appendix B.The thermal balance determines the temperature 𝑇 𝑒 and levelpopulation: Γ − ∑︁ 𝑖 Λ 𝑖 ≈ , (12)where Λ 𝑖 is the cooling rate of 𝑋 𝑖 + per unit volume. This equationcan be rewritten as Γ 𝑛 − 𝜒 ∑︁ 𝑖 (cid:18) Λ 𝑖 𝑛 𝑖 𝑛 𝑒 (cid:19) 𝑓 𝑖 ≈ , (13)where 𝑛 𝑖 = 𝑓 𝑖 𝑛 and 𝑛 𝑒 = 𝜒𝑛 . Note that the cooling rate due to free-free emission is much smaller than the atomic cooling rate in thetemperature range of the nebular phase and therefore we consideronly the atomic cooling in the following.The number density of 𝑋 𝑖 + in a level 𝑗 , 𝑛 𝑖, 𝑗 , is determined bya given kinetic temperature, density, electron fraction, and radiationfield. We use the escape probability approximation to solve the levelpopulation (e.g., Chapter 19 of Draine 2011): 𝑑𝑛 𝑖 𝑑𝑡 = ∑︁ 𝑗<𝑖 (cid:2) 𝑛 𝑗 𝑛 𝑒 𝑘 𝑗𝑖 − 𝑛 𝑖 𝑛 𝑒 𝑘 𝑖 𝑗 − (cid:104) 𝛽 𝑖 𝑗 (cid:105) 𝑛 𝑖 𝐴 𝑖 𝑗 (cid:3) + ∑︁ 𝑗>𝑖 (cid:2) 𝑛 𝑗 𝑛 𝑒 𝑘 𝑗𝑖 − 𝑛 𝑖 𝑛 𝑒 𝑘 𝑖 𝑗 + (cid:104) 𝛽 𝑖 𝑗 (cid:105) 𝑛 𝑗 𝐴 𝑗𝑖 (cid:3) ≈ , (14)where 𝑘 𝑖 𝑗 is the collisional rate coefficient for 𝑖 → 𝑗 , 𝐴 𝑖 𝑗 is theradiative transition rate for 𝑖 → 𝑗 , and (cid:104) 𝛽 𝑖 𝑗 (cid:105) is the escape probabilityof photons created by the transition 𝑖 → 𝑗 . Note that we omittedthe suffix that denotes the ionizing state in equation (14). Herewe use the Sobolev optical depth to evaluate (cid:104) 𝛽 𝑖 𝑗 (cid:105) (see AppendixC). Because this description includes only self-absorption of lines,the cooling function and spectrum of each ion can be computedseparately. However, this approximation is not valid around thefrequencies where the effect of the line overlapping is important.Such a situation can occur in the optical region for lanthanide-richnebulae as will be discussed in §4.The atomic cooling rate of 𝑋 𝑖 + is calculated by Λ 𝑖 = ∑︁ 𝑗> 𝑛 𝑖, 𝑗 ∑︁ 𝑘< 𝑗 𝐸 𝑗𝑘 (cid:104) 𝛽 𝑗𝑘 (cid:105) 𝐴 𝑗𝑘 , (15)where 𝐸 𝑗𝑘 is the energy difference between levels 𝑗 and 𝑘 . One canshow that the normalized cooling rate ( Λ 𝑖 / 𝑛 𝑒 𝑛 𝑖 ) depends only on 𝑇 𝑒 at sufficiently low densities. For NSM nebulae, as we will show later, ( Λ 𝑖 / 𝑛 𝑒 𝑛 𝑖 ) is almost independent of the density at 𝑛 (cid:46) cm − ,corresponding to 𝑡 (cid:38)
40 day for 𝑀 ej = . 𝑀 (cid:12) and 𝑣 = . 𝑐 .The ionization degree, 𝑓 𝑖 , free electron fraction, 𝜒 , and kinetictemperature, 𝑇 𝑒 , at each time are obtained by solving equations (11)and (13) iteratively for given Γ ( 𝑡 ) and 𝑛 ( 𝑡 ) . Roughly speaking, equa-tions (11) and (13) explicitly depend on time only through Γ / 𝑛 .Therefore the thermodynamic quantities evolve with time accord-ing to Γ / 𝑛 . This fact and equation (7) suggest that the ionizationdegree, electron fraction, and kinetic temperature roughly increaseas ∝ 𝑡 . for 𝑡 (cid:46) 𝑡 th and increase very slowly as ∝ 𝑡 . for 𝑡 (cid:38) 𝑡 th (figure 1). This property is somewhat naturally expected from thefact that almost all the processes relevant to the ionization and ther-mal balance after the thermalization breaks are two-body collisionbetween electrons and ions. NSM ejecta are composed of atoms with a wide range of atomicnumbers, 𝑍 (cid:38)
30, in reality. The experimental atomic data of these heavy elements are largely unavailable. To derive the atomicdata necessary for our purpose we use atomic structure codes,General Relativistic Atomic Structure Package (
GRASP2K ; Jönssonet al. 2013) and Hebrew University Lawrence Livermore Atomic(
HULLAC ; Bar-Shalom et al. 2001) codes.
HULLAC is an integratedcode for calculating atomic structures and cross sections for themodelings of atomic processes in plasmas and emission spec-tra, which employs a parametric potential method for calculationsof bound- and free-electron wavefunctions. The
GRASP2K codeprovides more rigorous bound-electron wavefunctions based onthe multiconfiguration Dirac–Hartree–Fock method, which enablesmore ab-initio calculations of atomic structures and bound-boundradiative transition probabilities, and therefore, we use
GRASP2K toderive the level spectra and radiative transition rates (see Gaigalaset al. 2019 for details) and compute recombination rate coefficientsby using
HULLAC .Lanthanide ions enhance the opacity of NSM ejecta (Kasenet al. 2013; Tanaka & Hotokezaka 2013; Tanaka et al. 2020; Barneset al. 2020), and thus, they are naturally expected to be strong emit-ters in the nebular phase. In addition, one may be able to capture,at least qualitatively, some important features of the nebular emis-sion of lanthanide-rich ejecta by using a single element because ofthe similarity in the spectral structure between lanthanide elements.Motivated by these, we focus on neodymium (Nd), in order to qual-itatively understand the nebular emission of lanthanide-rich ejectain this and following sections.
Before proceeding the details of the atomic data, here, we brieflysummarize some spectral properties of lanthanide ions (Gold-schmidt 1978). Lanthanide elements, Ln, are a group of elementswith atomic numbers 57 – 71 (La – Lu). Their ions are characterizedby the number of electrons in 4 𝑓 -shell, 𝑁 or 𝑁 +
1, where we usea number 𝑁 = 𝑍 −
57, e.g., 𝑁 = 𝑁 =
14 for Lu.Their configurations lying at low energies often have one to threeelectrons in the outer shells, 5 𝑑 and 6 𝑠 , which means that the energyscales of 4 𝑓 , 5 𝑑 , and 6 𝑠 are similar so that several different config-urations with the same parity consists of a group. Figure 2 shows acharacteristic spectral structure of first to third lanthanide ions (LnII-IV). For two groups connected by arrows, there are permittedtransitions between them. Forbidden transitions between differentorbital angular momenta, 𝐿 , as well as those associated with the finestructure are also important for the cooling rate and the emergentspectra. The characteristic spectra of lanthanides are summarizedas follows.(i) Permitted (E1) transitions of Ln II and Ln III exist in the nIRand optical bands. These lines lead to the enhancement of absorptionand can also be the source of nIR-optical emission in the nebularphase.(ii) Dipole forbidden (M1) transitions of Ln II - Ln IV betweendifferent configurations or between different total orbital angularmomenta produce emission lines in the nIR and optical bands.(iii) Transitions between fine stricture levels produce mid-IRlines ( 𝜆 ∼ 𝜇 m).Note that, among Ln II, Nd II has resonance lines at the lowesttransition energy and more E1 transitions at longer wavelengths.In fact, Tanaka et al. (2020) show that the abundance of Nd atomshas the most significant effect on the opacity in kilonovae. In thefollowing, we focus on the atomic data of Nd. MNRAS , 1–15 (2020) eutron Star Merger Nebula Ln II Ln III Ln IV N (5d+6s) N+1
6s 4f
N+1
N-1 (5d+6s) N N N+2 N+1 nIRoptical optical N+1 N
6s 4f N
5d 4f N N-1
6s 4f
N-1 optical nIR UV Figure 2.
Characteristic structure of the transitions from the lowest configurations for second to fourth lanthanide ions.
Ln II : The ground configurationof Ln II is one of the configurations in the two groups 4 𝑓 𝑁 ( 𝑑 + 𝑠 ) = 𝑓 𝑁 𝑑 , 𝑓 𝑁 𝑑 𝑠, 𝑓 𝑁 𝑠 and 4 𝑓 𝑁 + 𝑑, 𝑠 . The energy levels of theseconfigurations overlap and there exist permitted transitions between the two groups. This feature causes dense spectra in the nIR band. The transition linesfrom these lowest energy groups to the second lowest ones are typically in the optical region. Ln III and IV : The ground configurations of Ln III and Ln IV are4 𝑓 𝑁 + and 4 𝑓 𝑁 , respectively, except for La III, Gd III, and Lu III. The typical energy separations between the ground and the next lowest configurations arein the nIR to optical region for Ln III and in the ultraviolet region for Ln IV. Because only excited configurations overlap their spectra are less dense than thoseof Ln II. − E u [eV] o f li n e s p e r l og b i n Nd II
M1E1 − E u [eV] Nd III
M1E1 − E u [eV] Nd IV M1 T e [K] − − − − − − Λ / ( n e n i ) [ e r g c m s − ] Nd II , n i = 10 cm − totalE1M1 T e [K] Nd III , n i = 10 cm − totalE1M1 T e [K] Nd IV , n i = 10 cm − totalM1 Figure 3.
Number of transition lines per logarithmic interval of energy of excited states 𝐸 𝑢 computed by using GRASP2K ( top panels ). Red and blue histogramsshow M1 and E1 transitions, respectively. Cooling function of Nd II, Nd III, and Nd IV at an atomic number density of 10 cm − ( bottom panels ). This numberdensity corresponds to ∼
40 day after merger in the case of an ejecta mass of ∼ . 𝑀 (cid:12) and expansion velocity of ∼ . 𝑐 . Dotted and dashed curves show thecooling due to M1 and E1 transitions, respectively. For Nd II and Nd III, the contribution of E1 transitions is significant at 𝑇 𝑒 (cid:38)000
40 day after merger in the case of an ejecta mass of ∼ . 𝑀 (cid:12) and expansion velocity of ∼ . 𝑐 . Dotted and dashed curves show thecooling due to M1 and E1 transitions, respectively. For Nd II and Nd III, the contribution of E1 transitions is significant at 𝑇 𝑒 (cid:38)000 , 1–15 (2020) K. Hotokezaka et al. λ [ ˚ A ] − . − . − . . . . l og ( g f ) Nd II : even → odd GRASP2KDen Hartog 03 λ [ ˚ A ] − . − . − . − . − . . l og ( g f ) Nd II
GRASP2KHasselquist 16 λ [ ˚ A ] − . − . − . − . − . − . − . . l og ( g f ) Nd III : 4 f → f d GRASP2KRyabchikova 06
Figure 4.
Comparison of line strengths of Nd II and Nd III computed byusing
GRASP2K (Gaigalas et al. 2019) with the experimental results (DenHartog et al. 2003) and the APOGEE line list (Hasselquist et al. 2016) forNd II and the Nd III line list based on the stellar spectrum of a stronglymagnetic Ap star (Ryabchikova et al. 2006).
We include the excited levels of Nd ions up to ≈ ∼ than E1 transitions. We also examined E2 transitionsand found that their contribution to the cooling function is rather minor in the relevant temperature range and therefore we decide notto include E2 transitions.Note that there are excited states that can decay through E1transitions down to ∼ . GRASP2K with those from Den Hartog et al. (2003) andRyabchikova et al. (2006) in the optical region and that from Has-selquist et al. (2016) in the nIR region. Den Hartog et al. (2003)experimentally measured the wavelengths and oscillator strengthsof over 700 lines of Nd II. Here we focus on the intensive lines withlog 𝑔 𝑓 > − . 𝐸 𝑢 < − , and 𝐽 (cid:62) / < 𝜆 < ∼ GRASP2K statistically agreeswith the laboratory-based one. The
GRASP2K line distribution isalso roughly in agreement with the nIR lines of Nd II identifiedfrom the Apache Point Observatory Galactic Evolution Experiment(APOGEE) H-band spectra (Hasselquist et al. 2016).Because the line spectrum of Nd III is poorly known exper-imentally, here we compare the
GRASP2K result with the line listprovided by Ryabchikova et al. (2006), in which they propose theline classification for Nd III based on stellar spectra and a theoret-ical calculation of atomic structure. In figure 4, we show 23 linesassociated with the transitions between 4f and 4f
5d in the rangeof 4500 < 𝜆 < ∼ % level. We consider that the GRASP2K line list issufficiently accurate at least for E1 transitions to capture the spectralstructure of the NSM nebular emission.
We derive the collisional rate coefficients for the
GRASP2K atomicdata with the procedure described in Appendix D. Here we discussthe typical critical densities for Nd ions and implications to theevolution of the cooling functions and spectra. The critical densityfor a given upper level 𝑢 is estimated as 𝑛 crit ,𝑢 ≡ (cid:205) 𝑙<𝑢 𝐴 𝑢𝑙 (cid:205) 𝑙<𝑢 𝑘 𝑢𝑙 , (16) ∼ (cid:26) cm − (E1 transition) , cm − (M1 transition) , (17)where we used the typical value of 𝑘 𝑢𝑙 and 𝐴 𝑢𝑙 . These criticaldensities correspond to the critical times: 𝑡 crit ∼ (cid:26) ,
40 day (M1 transition) . (18)When the NSM ejecta becomes optically thin, the time scale ofE1 radiative deexcitation is much faster than that of excitation, i.e., 𝑡 crit ( E1 ) (cid:28) 𝑡 thin . Therefore, the level populations in the nebularphase are always far from those in collisional equilibrium, i.e., theLTE values. For 𝑡 > 𝑡 crit ( M1 ) , excited levels predominantly decaythrough radiative transition. Such a state is referred to as coronaequilibrium. In this case, the cooling rate is proportional to 𝑛 𝑒 𝑛 𝑖 ,i.e., the cooling function, Λ 𝑖 / 𝑛 𝑒 𝑛 𝑖 , is independent of the density,and therefore, the kinetic temperature is expected to evolve veryslowly with time after 𝑡 crit ( M1 ) because of Γ / 𝑛 ∝ 𝑡 . . MNRAS , 1–15 (2020) eutron Star Merger Nebula The bottom panels of figure 3 show the cooling functions of NdII, III and IV ions at an ion density of 10 cm − , corresponding to ∼
40 days after merger for 𝑀 ej = . 𝑀 (cid:12) and 𝑣 = . 𝑐 . We findthe overall trend of the cooling functions, Λ ( IV ) < Λ ( III ) < Λ ( II ) ,which can be understood from the fact that Nd II has more lines in theIR to optical region. M1 transitions dominate the cooling functionsof Nd II and Nd III for 𝑇 𝑒 < < left ) and thedensity effect ( right ) on the cooling rates. The cooling functionswithout the trapping effect are calculated with the assumption of (cid:104) 𝛽 𝑖 𝑗 (cid:105) =
1. We note that the absorption effect reduces the coolingfunction of Nd II by (cid:38)
50% at (cid:38) K, where E1 transitionsdominate the cooling rate. This effect is weaker for Nd III andabsent for Nd IV. Note that, however, we likely overestimate theescape probability of lines with 𝜆 (cid:46) 𝜇 m because these lines maybe absorbed by nearby permitted lines such as resonance lines (seemore details in §5). The density effect is quite small at the densitiesof lanthanide-rich NSM nebulae. Thus, for 𝑛 (cid:46) cm − , the ejectais in corona equilibrium and the cooling function can be consideredto be independent of the density. For the results presented in thefollowing section, the trapping and density effects are accountedfor. The cooling time scale is estimated as 𝑡 cool ≈ 𝑘𝑇 𝑒 Λ / 𝑛 ∼ s (cid:18) Λ / 𝑛 𝑒 𝑛 − erg cm s − (cid:19) − (cid:18) 𝑇 𝑒 𝐾 (cid:19) × (cid:18) 𝑀 ej . 𝑀 (cid:12) (cid:19) − (cid:18) (cid:104) 𝐴 (cid:105) (cid:19) (cid:16) 𝑣 . 𝑐 (cid:17) (cid:18) 𝑡 (cid:19) , (19)where Λ is the total cooling function. This time scale is much shorterthan a dynamical time until ∼
10 years after merger, and thus, thesteady-state approximation for thermal balance is valid on the timescale, 𝑡 (cid:46)
100 day, focused in this work.
Dielectronic recombination dominates over radiative recombinationfor lower ionized Nd ions. At nebular temperatures ( ∼ K), au-toionizing states lying slightly above the ionization threshold con-tribute to the dielectronic capture process so that resolving finestructure is important here. For this purpose, we use the level modeof
HULLAC to obtain the energy levels, radiative transition rates, andautoionization rates. With these quantities, we calculate the ratecoefficients by following the prescription of Nussbaumer & Storey(1983) (see also Appendix E). We include the following autoioniz-ing states: • Nd II: 4f 𝑛𝑙 ( 𝑛 (cid:54) , 𝑙 (cid:54) 𝑛𝑙 ( 𝑛 (cid:54) , 𝑙 (cid:54) 𝑛𝑙 ( 𝑛 (cid:54) , 𝑙 (cid:54) • Nd III: 4f 𝑛𝑙 ( 𝑛 (cid:54) , 𝑙 (cid:54) 𝑛𝑙 ( 𝑛 (cid:54) , 𝑙 (cid:54) 𝑛𝑙 ( 𝑛 (cid:54) , 𝑙 (cid:54) • Nd IV: 4f 𝑛𝑙 ( 𝑛 (cid:54) , 𝑙 (cid:54) 𝑛𝑙 ( 𝑛 (cid:54) , 𝑙 (cid:54) 𝑛𝑙 ( 𝑛 (cid:54) , 𝑙 (cid:54) 𝑛𝑙 ( 𝑛 (cid:54) , 𝑙 (cid:54) 𝛾𝑛𝑙 , where 𝛾 denotes thestate of the core electrons, 𝑛 and 𝑙 denote the principal and orbitalangular momentum quantum numbers of the captured electron. Wenote that the contribution of each configuration with higher 𝑛 and 𝑙 that is not included is less than ∼
1% for 𝑇 𝑒 (cid:46) K. Table 1.
Model parameters.model 𝑀 ej [ 𝑀 (cid:12) ] 𝑣 [ 𝑐 ] wind (fiducial) 0 .
05 0 . .
02 0 . .
05 0 . Figure 6 shows the recombination rate coefficients for Nd II- IV. The contribution of radiative recombination to the total ratecoefficient is less than 10% for all the cases. Note that, for Nd I, weassume that the rate coefficient of dielectronic recombination is 1 / 𝑡 rec ∼ 𝛼𝑛 𝑒 ∼ (cid:18) 𝛼 − cm s − (cid:19) − (cid:18) 𝑀 ej . 𝑀 (cid:12) (cid:19) − × (cid:18) (cid:104) 𝐴 (cid:105) (cid:19) (cid:16) 𝑣 . 𝑐 (cid:17) (cid:18) 𝑡
40 day (cid:19) , (20)where 𝛼 is the total recombination rate coefficient. The ionizationtime scale is estimated from the heating rate (see figure 1): 𝑡 ion ∼ Υ Nd III ∼ (cid:18) 𝑡
40 day (cid:19) . for 𝑡 (cid:38) 𝑡 th , (21)where we have used 𝑤 𝑖 ∼
60 eV for Nd III. For the fiducial model, 𝑀 ej = . 𝑀 (cid:12) and 𝑣 = . 𝑐 , these two time scales become com-parable to a dynamical time at ∼
100 day. Thus, we consider thenebular phase at (cid:46)
100 day after merger, where the steady-stateapproximation is valid.
By solving the equations described in §3 with the atomic data ofNd ions shown in §4, we obtain the evolution of the thermodynamicquantities and emergent spectrum in the NSM nebular phase (seeAppendix G for an application of our method to SN Ia nebulae).Table 1 shows the three cases studied here and we choose the windmodel, ( 𝑀 ej , 𝑣 ) = ( . 𝑀 (cid:12) , . 𝑐 ) , as the fiducial model.Figure 7 shows the evolution of the fractional ion abundancesand the kinetic temperature in the fiducial case. The temperatureslowly increases with time from ∼ 𝑡 th ≈
50 day, where the normalizedheating function changes its slope from ∝ 𝑡 . to 𝑡 . . The frac-tional ion abundances also very slowly change with time after thethermalization break.Figure 8 shows the temperature evolution for the dynamicalejecta and slow wind models. The characteristic temperatures forthe dynamical ejecta and slow wind models are ≈ K and 3 · K, respectively. The ionization degrees of dynamical ejecta andslow wind models are higher and lower than the fiducial model,respectively.The individual spectra of Nd II – IV at 𝑛 𝑖 = 𝑛 𝑒 = . · cm − and 𝑇 = MNRAS000
50 day, where the normalizedheating function changes its slope from ∝ 𝑡 . to 𝑡 . . The frac-tional ion abundances also very slowly change with time after thethermalization break.Figure 8 shows the temperature evolution for the dynamicalejecta and slow wind models. The characteristic temperatures forthe dynamical ejecta and slow wind models are ≈ K and 3 · K, respectively. The ionization degrees of dynamical ejecta andslow wind models are higher and lower than the fiducial model,respectively.The individual spectra of Nd II – IV at 𝑛 𝑖 = 𝑛 𝑒 = . · cm − and 𝑇 = MNRAS000 , 1–15 (2020)
K. Hotokezaka et al. T e [K] − − − − − − Λ / ( n e n i ) [ e r g c m s − ] Nd II , n e = 10 cm − total(trap)total(notrap)E1(trap)E1(notrap) T e [K] − − − − − − Λ / ( n e n i ) [ e r g c m s − ] Nd II cm − cm − cm − Figure 5.
Cooling function of Nd II with and without the radiation trapping effect ( left ) and at different densities ( right ). T e [K] − − − − − − α [ c m s − ] Nd III → Nd II
DR (total)DR : 4f nl DR : 4f nl DR : 4f nl RR T e [K] Nd IV → Nd III
DR (total)DR : 4f nl DR : 4f nl DR : 5p nl RR T e [K] Nd V → Nd IV
DR(total)DR : 5p nl DR : 5p nl DR : 4f nl DR : 4f5d nl RR Figure 6.
Rate coefficients for dielectronic recombination (DR) and radiative recombination (RR). The rate coefficients of dielectronic recombination areobtained by using
HULLAC . For dielectronic recombination, each line shows the contribution of a specific configuration of the autoionizing state 𝛾𝑛𝑙 , where 𝛾 denotes the core configuration, 𝑛 and 𝑙 denote the principal and orbital angular momentum quantum numbers of the captured electron. The range of 𝑛 and 𝑙 ofeach autoionizing state included in the calculation are described in the text.
20 40 60 80 100
Time since merger [day] − − F r a c t i o n a li o n a bund a n ce wind : 0 . M (cid:12) , . c IIIIIIIV
10 20 30 40 50 60 70 80 90 100
Time since merger [day] T e [ K ] wind : 0 . M (cid:12) , . c Figure 7.
Evolution of fractional ion abundances, Nd I – Nd IV ( left ) and kinetic temperature ( right ) for the fiducial model. The ionization degree and kinetictemperature increase until the thermalization time 𝑡 th and then become roughly constant with time. these ions have two distinct peaks, one around 5–10 𝜇 m producedby fine structure transitions and another around optical-nIR region.Nd II has among the richest spectral structure and its luminosity peratom is the brightest. The dense emission line distribution and theDoppler broadening result in a continuum-like spectrum with somestructures. We find that the following transitions predominatelyproduce the Nd II spectrum: 4f → → → → → → → and 4f → for Nd III and 4f → for Nd IV. Notethat individual M1 lines are more pronounced at 𝜆 (cid:46) 𝜇 m becausethe line population in this wavelength region is less dense.There are more E1 transition lines at 𝜆 (cid:46) 𝜇 m for Nd II andNd III (see figure 3). This implies that these E1 lines may absorb MNRAS , 1–15 (2020) eutron Star Merger Nebula
20 40 60 80 100
Time since merger [day] T e [ K ] dynamical ejectawind
50 100 150 200 250 300
Time since merger [day] slow windwind
Figure 8.
Kinetic temperature evolution for the dynamical ejecta model ( left : 𝑀 ej = . 𝑀 (cid:12) and 𝑣 = . 𝑐 ) and slow wind model ( right : 𝑀 ej = . 𝑀 (cid:12) and 𝑣 = . 𝑐 ). The fiducial model (wind) is also shown as a dash-dotted curve for comparison. The time scales on which the ejecta enters the nebular phase forthe dynamical ejecta and slow models are ≈
10 day and 70 day, respectively. ν [Hz] ν L ν / f i [ e r g s − ] totalE1M1 ν [Hz] totalE1M1 ν [Hz] totalM1 λ [ µ m] Nd II , n e = 1 . · cm − , T e = 4500 K λ [ µ m] Nd III , n e = 1 . · cm − T e = 4500 K λ [ µ m] Nd IV , n e = 1 . · cm − T e = 4500 K Figure 9.
Normalized spectra for Nd II, Nd III, and Nd IV. Here we use a kinetic temperature of 𝑇 𝑒 = 𝑛 𝑒 = . · cm − , andelectron fraction of 𝜒 =
1. These values roughly correspond to those around 40 day after merger in the fiducial model. Solid, dashed, and dash-dotted curvesdepict the total spectrum, the contribution of E1 transitions, and the contribution of M1 transitions, respectively. Also shown as vertical lines are individual E1(blue) and M1 (red) lines. The Doppler broadening of each line at a frequency 𝜈 𝑖 is incorporated by using a Gaussian distribution with a standard deviation of ≈ ( 𝑣 / 𝑐 ) 𝜈 𝑖 = . 𝜈 𝑖 . ν [Hz] ν L ν [ e r g s − ] ν [Hz] totalNd IINd IIINd IV λ [ µ m] t = 40 day 10 λ [ µ m] t = 80 day Figure 10.
Spectra for the fiducial model at 40 day ( left ) and 80 day ( right ) after merger. The contributions of Nd II, Nd III, Nd IV are also shown. Filledcircle and triangle are the detection at 4 . 𝜇 m and 3 𝜎 upper limit at 3 . 𝜇 m obtained by Spitzer telescope at 43 day ( left ) and 74 day ( right ) after GW170817(Kasliwal et al. 2019).MNRAS000
Spectra for the fiducial model at 40 day ( left ) and 80 day ( right ) after merger. The contributions of Nd II, Nd III, Nd IV are also shown. Filledcircle and triangle are the detection at 4 . 𝜇 m and 3 𝜎 upper limit at 3 . 𝜇 m obtained by Spitzer telescope at 43 day ( left ) and 74 day ( right ) after GW170817(Kasliwal et al. 2019).MNRAS000 , 1–15 (2020) K. Hotokezaka et al. ν [Hz] ν L ν [ e r g s − ] totalNd IIINd IV ν [Hz] totalNd IINd III λ [ µ m] dynamical ejecta, t = 20 day 10 λ [ µ m] slow wind, t = 80 day Figure 11.
Same as figure 10 but for the dynamical ejecta model at 20 day ( left ) and the slow wind model at 80 day ( right ). other emission lines and reduce the emission at 𝜆 (cid:46) 𝜇 m. In fact,Nd II and Nd III respectively have ∼
50 and ∼
15 resonance lines inthe range of 0 . (cid:46) 𝜆 (cid:46) 𝜇 m and 0 . (cid:46) 𝜆 (cid:46) . 𝜇 m. This radiationtransfer effect is not accounted for in our modeling, and thus, ourmodeling likely overpredicts the optical emission.Figure 10 shows the total spectra at 40 and 80 day for thefiducial model with the fractional ion abundances shown in figure7. The Nd II lines dominate the total spectrum particularly in thenIR band. The spectral shape does not change significantly from 40to 80 day while the amplitude decreases by a factor of ∼
10. Thisfreeze-out of the nebular spectrum is a characteristic feature of theNSM nebular emission.The spectra of the dynamical ejecta and slow wind models areshown in figure 11. For dynamical ejecta, each line is significantlybroaden because of the fast expansion velocity, 0 . 𝑐 . As a result,the structures are completely smeared out. Nevertheless, there aretwo distinct peaks around the optical and IR bands. On the contrary,for the slow wind model, more lines can be seen in the IR region(1 (cid:46) 𝜆 (cid:46) 𝜇 m) and the optical emission is very weak. The spectralshape does not evolve significantly during the nebular phase in theboth models.We show the detectability of the structure of the nebular spec-trum by the James Webb Space Telescope ( JWST ) for a future kilo-nova event in figure 12. The
JWST is promising to resolve thespectral structure of the nebular emission around 40 day for eventsout to ∼
100 Mpc.
The emission-line nebular phase of the NSM ejecta is studied byusing a one-zone nebula model under non-LTE, in which the ejectais considered to be composed of one of lanthanide elements, Nd.The atomic data necessary for the modeling are calculated by usingthe atomic structure codes,
GRASP2K and
HULLAC . We find that thekinetic temperature and ionization fraction are nearly constant withtime after the thermalization break of the beta-decay heating rate.Consequently, the spectral shape of the emergent emission is alsoexpected to be frozen after the break. For the ejecta parameters of 𝑀 ej = . 𝑀 (cid:12) and 𝑣 = . 𝑐 , we show that Nd II and Nd IIIare the most abundant ions and the kinetic temperature approaches ≈ 𝛽 -decay heating rate re-sults in a deviation in the ionization state from LTE. In particular, λ [ µ m] − − F l u x d e n s i t y a t M p c [ µ J y ] JWST NIRSpec JWST MIRI
Figure 12.
Spectrum at 40 day for the fiducial model at a distance of100 Mpc. Also depicted are the 1 𝜎 sensitivity curves of NIRSpec FS andMIRI LRS with 10 s integration. we find that the neutral fraction is significantly suppressed in thenebular phase. Although we do not account for the velocity distri-bution in this work, we speculate that this deviation can occur evenat the earlier times, e.g., (cid:46) ∼ . 𝜇 m to 20 𝜇 m with two dis-tinct peaks around ∼ 𝜇 m and ∼ 𝜇 m. Fine-structure transitionsproduce the mid-IR peak. This spectral structure may be an uniquefeature of lanthanide-rich nebulae. It is worth emphasizing that indi-vidual M1 lines are more pronounced at 𝜆 (cid:46) 𝜇 m because the linepopulation in this wavelength region is less dense. Importantly, the JWST will be able to resolve such structure in the nIR and midIRregions for events at ∼
100 Mpc. Note, however, that this struc-ture may be suppressed once more elements are included. Anothercaveat of our modelling is that we have neglected the absorptiondue to line overlapping, which may lead to an overestimate of theoptical-nIR emission ( 𝜆 (cid:46) 𝜇 m), where Nd II and Nd III have anumber of permitted lines.We use a crude approximation for the collisional strength offorbidden lines, i.e., Ω 𝐹 =
1, in the case of the
GRASP2K calculation.
MNRAS , 1–15 (2020) eutron Star Merger Nebula While this approximation is statistically consistent with the colli-sional strengths derived with
HULLAC and can reasonably reproducethe cooling rates, the predicted line intensity ratios are by no meansaccurate. Thus, we need more accurate collisional strengths for thefuture studies.
ACKNOWLEDGMENTS
We thank B. T. Draine, M. M. Kasliwal, K. Kawaguchi, and E.Nakar for useful discussion and M. Busquet for the generous supporton the
HULLAC code. K. H. was supported by Japan Society for thePromotion of Science (JSPS) Early-Career Scientists Grant Number20K14513.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Abbott B. P., et al., 2017, ApJ, 848, L12Andreoni I., et al., 2017, Publ. Astron. Soc. Australia, 34, e069Arcavi I., et al., 2017, Nature, 551, 64Arnett W. D., Fryer C., Matheson T., 2017, ApJ, 846, 33Axelrod T. S., 1980, PhD thesis, California Univ., Santa Cruz.Bar-Shalom A., Klapisch M., Oreg J., 2001, J. Quant. Spectrosc. RadiativeTransfer, 71, 169Barnes J., Kasen D., 2013, ApJ, 775, 18Barnes J., Zhu Y. L., Lund K. A., Sprouse T. M., Vassh N., McLaugh-lin G. C., Mumpower M. R., Surman R., 2020, arXiv e-prints, p.arXiv:2010.11182Beigman I. L., Chichkov B. N., 1980, Journal of Physics B Atomic MolecularPhysics, 13, 565Bohr N., 1913, The London, Edinburgh, and Dublin Philosophical Magazineand Journal of Science, 25, 10Botyánszki J., Kasen D., Plewa T., 2018, ApJ, 852, L6Bulla M., 2019, MNRAS, 489, 5037Coulter D. A., et al., 2017, Science, 358, 1556Cowperthwaite P. S., et al., 2017, Astrophys. J., 848, L17Den Hartog E. A., Lawler J. E., Sneden C., Cowan J. J., 2003, ApJS, 148,543Draine B. T., 2011, Physics of the Interstellar and Intergalactic MediumDrout M. R., et al., 2017, Science, 358, 1570Evans P. A., et al., 2017, Science, 358, 1565Fransson C., Chevalier R. A., 1989, ApJ, 343, 323Gaigalas G., Kato D., Rynkun P., Radži¯ut˙e L., Tanaka M., 2019, ApJS, 240,29Gillanders J. H., McCann M., Smartt S. A. S. S. J., Ballance C. P., 2021,arXiv e-prints, p. arXiv:2101.08271Goldschmidt Z. B., 1978, in Handbook on the Physics andChemistry of Rare Earths, Vol. 1, Metals. Elsevier, pp 1– 171, doi:https://doi.org/10.1016/S0168-1273(78)01005-3,
Hasselquist S., et al., 2016, ApJ, 833, 81Hoffknecht A., et al., 1998, Journal of Physics B: Atomic, Molecular andOptical Physics, 31, 2415Hotokezaka K., Nakar E., 2020, ApJ, 891, 152Hotokezaka K., Wanajo S., Tanaka M., Bamba A., Terada Y., Piran T., 2016,MNRAS, 459, 35Hotokezaka K., Beniamini P., Piran T., 2018, International Journal of Mod-ern Physics D, 27, 1842005 Jönsson P., Gaigalas G., Bieroń J., Fischer C. F., Grant I. P., 2013, ComputerPhysics Communications, 184, 2197Kasen D., Barnes J., 2019, ApJ, 876, 128Kasen D., Badnell N. R., Barnes J., 2013, ApJ, 774, 25Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, Na-ture, 551, 80Kasliwal M. M., et al., 2017, Science, 358, 1559Kasliwal M. M., et al., 2019, MNRAS, p. L14Kawaguchi K., Shibata M., Tanaka M., 2018, ApJ, 865, L21Kawaguchi K., Fujibayashi S., Shibata M., Tanaka M., Wanajo S., 2020,arXiv e-prints, p. arXiv:2012.14711Korobkin O., et al., 2020, ApJ, 889, 168Kozma C., Fransson C., 1992, ApJ, 390, 602Lattimer J. M., Schramm D. N., 1974, ApJ, 192, L145Li L.-X., 2019, ApJ, 872, 19Maeda K., Nomoto K., Mazzali P. A., Deng J., 2006, ApJ, 640, 854Margutti R., Chornock R., 2020, arXiv e-prints, p. arXiv:2012.04810Mazzali P. A., et al., 2006, ApJ, 645, 1323Mazzali P. A., et al., 2015, MNRAS, 450, 2631Metzger B. D., et al., 2010, MNRAS, 406, 2650Nahar S. N., 1996, Phys. Rev. A, 53, 2417Nahar S. N., 1997, Phys. Rev. A, 55, 1980Nahar S. N., Bautista M. A., Pradhan A. K., 1997, ApJ, 479, 497Nakar E., 2020, Phys. Rep., 886, 1Nussbaumer H., Storey P. J., 1983, A&A, 126, 75Perego A., Radice D., Bernuzzi S., 2017, ApJ, 850, L37Perego A., et al., 2020, arXiv e-prints, p. arXiv:2009.08988Pian E., et al., 2017, Nature, 551, 67Rosswog S., Sollerman J., Feindt U., Goobar A., Korobkin O., WollaegerR., Fremling C., Kasliwal M. M., 2018, A&A, 615, A132Ruiz-Lapuente P., Lucy L. B., 1992, ApJ, 400, 127Ryabchikova T., Ryabtsev A., Kochukhov O., Bagnulo S., 2006, A&A, 456,329Schippers S., et al., 2011, Phys. Rev. A, 83, 012711Shibata M., Fujibayashi S., Hotokezaka K., Kiuchi K., Kyutoku K.,Sekiguchi Y., Tanaka M., 2017, Phys. Rev. D, 96, 123012Smartt S. J., et al., 2017, Nature, 551, 75Spencer L. V., Fano U., 1954, Physical Review, 93, 1172Spruck K., et al., 2014, Phys. Rev. A, 90, 032715Storey P. J., 1981, MNRAS, 195, 27PTanaka M., Hotokezaka K., 2013, ApJ, 775, 113Tanaka M., et al., 2017a, Publ. Astron. Soc. Jap.Tanaka M., et al., 2017b, PASJ, 69, 102Tanaka M., Kato D., Gaigalas G., Kawaguchi K., 2020, MNRAS, 496, 1369Tanvir N. R., et al., 2017, ApJ, 848, L27Utsumi Y., et al., 2017, PASJ, 69, 101Villar V. A., et al., 2017, ApJ, 851, L21Villar V. A., et al., 2018, ApJ, 862, L11Wanajo S., 2018, ApJ, 868, 65Watson D., et al., 2019, Nature, 574, 497Waxman E., Ofek E. O., Kushnir D., Gal-Yam A., 2018, MNRAS, 481, 3423Waxman E., Ofek E. O., Kushnir D., 2019, ApJ, 878, 93Wollaeger R. T., et al., 2018, MNRAS, 478, 3298Wu M.-R., Barnes J., Martínez-Pinedo G., Metzger B. D., 2019a, Phys.Rev. Lett., 122, 062701Wu M.-R., et al., 2019b, ApJ, 880, 23Yagi S., Nagata T., 2001, J. Phys. Soc. Jpn., 9, 2559Zhu Y., et al., 2018, ApJ, 863, L23van Regemorter H., 1962, ApJ, 136, 906
APPENDIX A: WORK PER ION PAIR
The ionization efficiency of fast electrons for a stopping plasma is inprinciple obtained by solving the Boltzmann equation under someapproximations (Spencer & Fano 1954; Kozma & Fransson 1992).Here we take a simple approach employed by Axelrod (1980), which
MNRAS000
MNRAS000 , 1–15 (2020) K. Hotokezaka et al. describes the radioactive ionization rate in terms of work per ionpair. The work per ion pair of 𝑋 𝑖 + is defined by 𝑤 𝑖 = 𝑓 𝑖 𝐸 dis 𝑁 𝑖 , (A1)where 𝑓 𝑖 is the number fraction 𝑋 𝑖 + , 𝐸 dis is the total dissipatedenergy of injected fast elections and 𝑁 𝑖 is the total number of ionpairs (ion-electron pairs) of 𝑋 ( 𝑖 + )+ produced through by the fastelectrons. The value of 𝑤 𝑖 simply represents the amount of energythat is dissipated in each ion-electron pair production.Let us consider first the work per ion pair for a primary electronwith an initial kinetic energy of 𝐸 𝑝 injected in a stopping plasma.The number of ion pairs of 𝑋 ( 𝑖 + )+ through the thermalization ofthe primary is given by 𝑁 𝑖, 𝑝 = ∫ 𝑠 th 𝑛 𝑖 𝜎 𝑖 ( 𝐸 ( 𝑠 )) 𝑑𝑠, (A2)where 𝑛 𝑖 is the number density of 𝑋 𝑖 + , 𝑠 is the travel distance of theelectron, and 𝜎 𝑖 is the ionization cross section. The energy loss perdistance interval is 𝑑𝐸𝑑𝑠 ≈ − 𝑛 ( 𝐿 ( 𝐸 ) + 𝐿 th ( 𝐸, 𝜒 )) , (A3)where 𝐿 and 𝐿 th are the stopping cross sections due to collisionalionization and excitation, and due to the Coulomb collision withthermal electrons, respectively. Equation (A2) is rewritten as 𝑁 𝑖, 𝑝 = 𝑓 𝑖 ∫ 𝐸 𝑝 𝑘𝑇 𝑒 𝜎 𝑖 ( 𝐸 ) 𝐿 ( 𝐸 ) + 𝐿 th ( 𝐸, 𝜒 ) 𝑑𝐸. (A4)The work per ion pair of the primary electron is then 𝑤 𝑝𝑖 = 𝐸 𝑝 ∫ 𝐸 𝑝 𝑑𝐸𝜎 𝑖 /( 𝐿 ( 𝐸 ) + 𝐿 th ( 𝐸, 𝜒 )) . (A5)To evaluate equation (A5), we use the total ionization crosssection of 𝑋 𝑖 + by electron-ion collision given by Axelrod (1980) 𝜎 𝑖 ≈ 𝜋𝑒 𝑚𝑣 𝑁 ∑︁ 𝑗 = 𝑞 𝑗 𝑃 𝑗 (cid:20) ln (cid:18) 𝑚 𝑒 𝑣 𝑃 𝑗 (cid:19) − ln ( − 𝛽 ) − 𝛽 (cid:21) , (A6)where 𝑣 is the electron’s velocity, 𝛽 = 𝑣 / 𝑐 , 𝑞 𝑗 and 𝑃 𝑗 are thenumber of electrons and the ionization potential of a subshell 𝑗 . Wenote that this formula (A6) agrees with the experimental data (Yagi& Nagata 2001). The stopping power for electrons is given by theBethe formula: 𝐿 ( 𝐸 ) ≈ 𝜋𝑍𝑒 𝑚 𝑒 𝑣 (cid:34) ln (cid:32) √︁ 𝑚 𝑒 𝑣 𝑇 / (cid:104) 𝐼 (cid:105)( − 𝛽 ) / (cid:33) (A7) −( √︃ − 𝛽 − − 𝛽 ) ln 2 + − 𝛽 + ( − √︃ − 𝛽 ) (cid:21) , where 𝑍 is the charge of the target ion, 𝑇 is the kinetic energy ofthe electron, and (cid:104) 𝐼 (cid:105) is the mean ionization energy of the stoppingmaterial. The value of (cid:104) 𝐼 (cid:105) is taken from the ESTAR database . Thestopping power of thermal plasma for with thermal velocity 𝑣 th (cid:28) 𝑣 is given by (Bohr 1913) 𝐿 th ( 𝐸, 𝜒 ) ≈ 𝜋𝑒 𝜒𝑚 𝑒 𝑣 ln (cid:32) . 𝑚 𝑒 𝑣 𝑒 𝜔 𝑝 (cid:33) , (A8)where 𝜔 𝑝 is the plasma frequency, 𝑛 𝑒 is the electron number density,and 𝜒 is the electron fraction 𝑛 𝑒 / 𝑛 . https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html For comparison between different ions, it is useful to definework per ion pair normalized by the first ionization potential, 𝐼 𝑖, .For instance, in the case of 𝐸 𝑝 =
250 keV and 𝜒 =
2, we find 𝑤 𝑝𝑖 / 𝐼 𝑖, ∼
45, 45, 40, and 35 for Nd I, Nd II, Nd III, and Nd IV,respectively. In addition to ionization by primary electrons, sec-ondary electrons may cause further ionization. This means thatthe total number of ion pairs 𝑁 𝑖 in equation (A1) is larger than 𝑁 𝑖, 𝑝 . Secondaries are ejected with recoil energy typically aroundthe binding energy of the target electron. For a weakly ionized Ndplasma, the stopping power of thermal electrons dominates overthe ionization energy loss at electron energies (cid:46) (cid:104) 𝐼 (cid:105) ≈ . ∼ 𝑤 𝑖 / 𝐼 𝑖, ∼
30. In this paper, we use 𝑤 𝑖 / 𝐼 𝑖, = APPENDIX B: PHOTOIONIZATION
The recombination processes emit photons that may be reprocessedby photoelectric absorption. This reprocess reduces the recombi-nation rate. The recombination of 𝑋 𝑗 + may emit ionizing photonsfor 𝑋 ( 𝑖< 𝑗 )+ . The number of ion pairs of 𝑋 ( 𝑗 + )+ and 𝑒 − producedby photoionization due to the photons emitted in a recombinationprocess of 𝑋 ( 𝑗 + )+ → 𝑋 𝑗 + is estimated by 𝑃 𝑖 𝑗 = ∫ (cid:0) − 𝑒 − 𝜏 (cid:1) 𝑓 𝑖 𝜎 ph 𝑖 (cid:205) 𝑘 𝑓 𝑘 𝜎 ph 𝑘 (cid:18) 𝑑𝑁 ph 𝑑𝜈 (cid:19) 𝑗 𝑑𝜈, (B1)where ( 𝑑𝑁 ph / 𝑑𝜈 ) 𝑗 is the number of photons per frequency intervalemitted in recombination of 𝑋 ( 𝑗 + )+ . Here the optical depth forphotons with frequency 𝜈 is given by 𝜏 ( 𝜈 ) ≈ ∑︁ 𝑘 = 𝑓 𝑘 𝜎 ph 𝑘 ( 𝜈 ) 𝑛𝑅, (B2)where 𝑓 𝑘 and 𝜎 ph 𝑘 are the number fraction and the photoionizationcross section of 𝑋 𝑘 + , and 𝑅 is the radius of the ejecta. Because 𝜎 ph is O( − cm ) , the optical depth is ∼ ( 𝑡 /
40 day ) − andtherefore the ejecta is optically thick for recombination photons inthe nebular phase.As shown in figure 6, Nd ions recombine predominantlythrough dielectronic recombination, where recombination photonsare produced through the radiative cascade of auto-ionization statesto the ground state. Therefore, it is not straightforward to deter-mine the recombination photon spectrum ( 𝑑𝑁 ph / 𝑑𝜈 ) 𝑗 . For auto-ionization states that have a large radiative transition rate, eachauto-ionization state contributes substantially to the rate coefficienteven though the number of such states is relatively small. In thischannel, auto-ionization states are typically stabilized through theemission of a photon with energy close to the first ionization poten-tial of the recombined ion and therefore this cascade produces oneionizing photon and several low energy photons. At the same time,there are many auto-ionization states that are stabilized through theemission of photons with energy sufficiently lower than the firstionization potential but high enough to ionize ions in lower ionizedstates. As a result, the recombination photons are likely to have asomewhat flat spectrum per logarithmic frequency interval. Thus,we assume that the number of recombination photons is constantat each energy scale below the sum of the first ionization potential MNRAS , 1–15 (2020) eutron Star Merger Nebula and the thermal energy of free electrons, i.e., ( 𝑑𝑁 ph / 𝑑𝜈 ) 𝑗 ∝ / 𝜈 for ℎ𝜈 < 𝐼 𝑗, + 𝑘𝑇 𝑒 and its normalization is set such that the totalenergy of recombination photons is 𝐼 𝑗, + 𝑘𝑇 𝑒 . APPENDIX C: SELF-ABSORPTION OF STRONG LINES
The absorption due to strong lines may have significant impacts onthe cooling functions and emergent spectra. In general, absorptionoccurs non-locally so that one must solve radiation transfer, whichis beyond the framework of our one-zone modeling. Here we use theescape probability approximation, which allows to include the ef-fects of self-absorption of lines in one-zone modelings (e.g, Chapter19 of Draine 2011).In homologously expanding ejecta, the escape probability isapproximated by (cid:104) 𝛽 𝑖 𝑗 (cid:105) = − 𝑒 − 𝜏 𝑠,𝑖𝑗 𝜏 𝑠,𝑖 𝑗 , (C1)where 𝜏 𝑠,𝑖 𝑗 is the Sobolev optical depth: 𝜏 𝑠,𝑖 𝑗 = 𝑔 𝑖 𝐴 𝑖 𝑗 𝜋 𝜆 𝑖 𝑗 (cid:18) 𝑛 𝑗 𝑔 𝑗 − 𝑛 𝑖 𝑔 𝑖 (cid:19) 𝑡 ( 𝑖 > 𝑗 ) . (C2)For resonance lines, the optical depth is estimated by 𝜏 𝑠, res ∼ (cid:18) 𝑔 𝑢 𝐴 𝑢 s − (cid:19) (cid:18) 𝜆 𝑢 . 𝜇 m (cid:19) (cid:18) 𝑀 ej . 𝑀 (cid:12) (cid:19) × (cid:16) 𝑣 . 𝑐 (cid:17) − (cid:18) 𝑡
30 day (cid:19) − . (C3)This suggests that the resonance lines are trapped in the ejecta ontime scales focused in this paper, (cid:46)
100 day.
APPENDIX D: COLLISIONAL EXCITATION ANDDEEXCITATION
With the usual convention, the velocity averaged rate coefficient forcollisional deexitation from an upper level 𝑖 to a lower level 𝑗 isgiven by 𝑘 𝑖 𝑗 = . · − Ω 𝑖 𝑗 ( 𝑇 𝑒 ) 𝑔 𝑖 𝑇 / 𝑒 cm s − , (D1)where Ω 𝑖 𝑗 is the velocity averaged collision strength connectinglevels 𝑖 and 𝑗 . The collisional excitation rate coefficient is given by 𝑘 𝑗𝑖 = 𝑔 𝑖 𝑔 𝑗 𝑘 𝑖 𝑗 𝑒 − 𝐸 𝑖𝑗 / 𝑘𝑇 𝑒 , (D2)where 𝑔 𝑖 is the level degeneracy, 𝐸 𝑖 𝑗 is the energy-level difference.The collisional strengths are currently not available for the GRASP2K atomic data. Therefore, we use the following approxi-mations for the collisional strengths Ω 𝑖 𝑗 for the GRASP2K atomicdata. For E1 transitions, we calculate Ω 𝑖 𝑗 by using the approximateformula (van Regemorter 1962): Ω 𝑖 𝑗 ≈ . 𝑃 (cid:18) 𝐸 𝑖 𝑗 𝑘𝑇 𝑒 (cid:19) (cid:18) 𝜆 𝜇 m (cid:19) (cid:18) 𝑔 𝑖 𝐴 𝑖 𝑗 𝑠 − (cid:19) , (D3)where 𝑃 ( 𝑥 ) is the Gaunt factor integrated over the electron velocitydistribution and 𝐴 𝑖 𝑗 is computed with GRASP2K . Here we use 𝑃 ( 𝑥 ) ≈ .
2, which is a good approximation for 𝑥 (cid:46) Ω 𝑖 𝑗 = Ω 𝐹 , where Ω 𝐹 is a constant value. Figure C1 shows the collisional strengths for M1 transitions at 𝑇 𝑒 = HULLAC . We findthat the averaged values around 𝐸 𝑖 𝑗 (cid:46) Ω 𝐹 ≈ GRASP2K with that of
HULLAC .The cooling functions due to M1 transitions derived with the twocodes are in a good agreement. This fact justifies our choice of Ω 𝐹 ≈
1. However, the E1 transition cooling of
GRASP2K is higherthan that of
HULLAC , suggesting that the van Regemorter formulaslightly overestimates the collisional strengths. Figure C3 shows thespectrum of each ion at 𝑛 = . · cm − and 𝑇 𝑒 = HULLAC . We note that the spectralstructures computed with the two codes are qualitatively similarbut the Nd II spectrum of
HULLAC has significant emission around3 𝜇 m. APPENDIX E: DIELECTRONIC RECOMBINATION
Dielectronic recombination occurs via the following process: 𝑋 ( 𝑖 + )+ 𝑝 + 𝑒 − → 𝑋 𝑖 + 𝑎 → 𝑋 𝑖 + 𝑏 + ℎ𝜈, (E1)where 𝑎 and 𝑏 denote an autoionizing state of 𝑋 𝑖 + and a bound stateof 𝑋 𝑖 + , respectively. The bound state, 𝑋 𝑖 + 𝑏 , is stabilized by radiativedecays. At the nebular temperature, radiative decays of both thecore and captured electrons contribute to the stabilization of 𝑋 𝑖 + 𝑎 (Beigman & Chichkov 1980; Storey 1981).The dielectronic recombination rate coefficient of the captureprocess (E1) is calculated by 𝛼 di ( 𝑝, 𝑎 ; 𝑇 𝑒 ) = (cid:32) 𝑁 𝑆 ( 𝑋 𝑖 + 𝑎 ) 𝑁 𝑒 𝑁 𝑆 ( 𝑋 ( 𝑖 + )+ 𝑝 ) (cid:33) (cid:205) 𝑗 𝐴 𝑎 𝑗 (cid:205) 𝑐 Γ 𝑎𝑐 (cid:205) 𝑐 Γ 𝑎𝑐 + (cid:205) 𝑘 𝐴 𝑎𝑘 , (E2)where 𝑁 𝑆 ( 𝑋 𝑖 + 𝑎 ) 𝑁 𝑒 𝑁 𝑆 ( 𝑋 ( 𝑖 + )+ 𝑝 ) = 𝑔 𝑎 𝑔 𝑝 (cid:18) ℎ 𝜋𝑚 𝑒 𝑘𝑇 𝑒 (cid:19) / 𝑒 − 𝐸 𝑎 / 𝑘𝑇 𝑒 , (E3)where 𝐸 𝑎 is the energy of a state 𝑎 relative to 𝑋 ( 𝑖 + )+ 𝑝 , and 𝑔 𝑝 is thestatistical weight of the state 𝑝 , Γ 𝑎𝑐 is the autoionization rate. Herethe sum for 𝑗 is over the levels that are stable against autoionizationand the sum for 𝑘 is over all the lower levels. We assume that therecombining ion 𝑋 ( 𝑖 + )+ 𝑝 is in the ground state. Then the total ratecoefficient is 𝛼 totdi ( 𝑇 𝑒 ) = ∑︁ 𝑎 𝛼 di ( 𝑝, 𝑎 ; 𝑇 𝑒 ) . (E4)This capture is a resonant process such that 𝐸 ( 𝑒 − ) = 𝐸 ( 𝑋 𝑖 + 𝑎 ) − 𝐸 ( 𝑋 ( 𝑖 + )+ 𝑝 ) must be satisfied and the autoinizing states that areaccessible via collision with thermal electrons contribute to thecapture rate. This indicates that ions with denser autoionizing statessuch as open f-shell ions have larger recombination rate coefficients.In fact, the measured values of the dielectronic recombination ratecoefficient of Au + , W + , and W + , nearly half open f-shell ions,are larger than the radiative recombination rate coefficient by twoto three orders of magnitude at nebular temperatures (Hoffknechtet al. 1998; Schippers et al. 2011; Spruck et al. 2014).For the nebular temperatures ( (cid:46) K), the kinetic energy ofthermal electrons is typically much smaller than the first ionizationpotential of ions, and therefore, autoionizing states only slightlyabove the ionization threshold contribute to the recombination pro-cess. Thus, levels in a small energy range from 𝐼 to ∼ 𝐼 + HULLAC thatresolves the fine structure.
MNRAS000
MNRAS000 , 1–15 (2020) K. Hotokezaka et al. E ij [eV] − − Ω i j ( M ) Nd II mean E ij [eV] Nd III mean E ij [eV] Nd IV mean
Figure C1.
Velocity averaged collision strengths for M1 transitions at 𝑇 𝑒 = HULLAC . Theaveraged value is shown by a blue line. T e [K] − − − − − − Λ / ( n e n i ) [ e r g c m s − ] Nd II , n i = 10 cm − GRASP2KHULLAC : fullHULLAC : approx . T e [K] Nd III , n i = 10 cm − GRASP2KHULLAC : fullHULLAC : approx . T e [K] Nd IV , n i = 10 cm − GRASP2KHULLAC : fullHULLAC : approx . Figure C2.
Cooling functions for different atomic data.
Solid curve shows the cooling function for the
GRASP2K atomic data with the formula (D3) and Ω 𝐹 =
1. Also shown are the cooling functions derived the
HULLAC atomic data, where collisional strengths are computed ( dashed curve ). Dotted curve showsthe cooling function with the
HULLAC energy levels and 𝐴 𝑖 𝑗 but collisional strengths are approximated by the formula (D3) and Ω 𝐹 = ν [Hz] ν L ν / f i [ e r g s − ] totalE1M1 ν [Hz] totalE1M1 ν [Hz] totalM1 λ [ µ m] Nd II , n e = 1 . · cm − , T e = 6000 K λ [ µ m] Nd III λ [ µ m] Nd IV
Figure C3.
Same as figure 9 but the energy levels, radiative transition rates, and collisional strengths are computed with
HULLAC . Here we set 𝑇 𝑒 = GRASP2K due to the low cooling functionof Nd II of
HULLAC . APPENDIX F: RADIATIVE RECOMBINATION
Radiative (direct) recombination occurs via 𝑋 ( 𝑖 + )+ 𝑝 + 𝑒 − → 𝑋 𝑖 + 𝑏 + ℎ𝜈. (F1)A photon produced by the recombination of 𝑋 𝑖 + directly to theground state is most likely absorbed by 𝑋 𝑖 + . This rate coefficient of this process is denoted customary as 𝛼 𝐴 and that of the recombina-tion to the other states is denoted 𝛼 𝐵 . Axelrod (1980) provides 𝛼 𝐴 ( 𝑇 ) = − 𝑖 (cid:18) 𝑇 K (cid:19) − / cm s − , (F2) MNRAS , 1–15 (2020) eutron Star Merger Nebula
100 200 300 400 500
Time since explosion [day] − − F r a c t i o n a li o n a bund a n ce SN Ia
IIIIIIIV
Figure F1.
Evolution of fractional ion abundances of Fe I – IV. Here weuse 𝑀 Ni = . 𝑀 (cid:12) and 𝑣 = − . λ [ ˚ A ] . . . . F λ [ e r g s − c m − ˚ A − ] × Fe IIFe II Fe IIIFe III modelSN2011fe , λ [ ˚ A ] . . . . . . . F λ [ e r g s − c m − ˚ A − ] × Fe IIFe IFe II , IIIFe IIFe II , III modelSN2011fe , Figure F2.
Spectra of a pure Fe nebula. Here we use 𝑀 Ni = . 𝑀 (cid:12) and 𝑣 = − . Ω 𝐹 = . HULLAC and the simple approximated method. Also shown is the ob-served nebular spectrum of SN2011fe (Mazzali et al. 2015). Because ourmodel includes only Fe ions some lines are absent in the synthetic spectra,e.g., Co III lines around 6000. and 𝛼 𝐵 ( 𝑇 ) = · − 𝑖 (cid:34)(cid:18) 𝑇 K (cid:19) − / − (cid:18) 𝑇 K (cid:19) − / (cid:35) cm s − . (F3)This form is provided for iron but it is not significantly differentfor heavy elements. We include the case B recombination (equationF3) in our modeling. APPENDIX G: NEBULAR SPECTRA OF SNE IA
Our nebula modeling is by no means accurate because we use anumber of approximations and assumptions. In order to show theability of our simple modeling, here we apply our method to the neb-ular emission of SNe Ia for comparison. Here we consider the decaychain of Ni → Co → Fe as the heat and ionization sourceand the heating rate computed by a code developed by Hotokezaka& Nakar (2020). As we did for NSM nebulae, we assume that theatomic properties are represented by a single atomic species, Fe. Thework per ion pair for Fe ions for primary electrons is 𝜔 𝑝 / 𝐼 ≈ 𝜔 / 𝐼 ≈
25. Because the properties of the transition linesof Fe ions relevant to the SN Ia nebula modelings are experimen-tally known, we use the NIST line list instead of preparing themwith the atomic codes. The collisional strengths are computed inthe prescription shown in Appendix D. Here we use Ω 𝐹 = . HULLAC in the relevant temperature range.Note that our one-zone modeling is fully characterized by only twoparameters: the total Ni mass, 𝑀 Ni , and the ejecta velocity, 𝑣 .We discussed in §4.5 that dielectronic recombination domi-nates over radiative recombination for Nd ions. Likewise, dielec-tronic recombination is more important for lower ionized Fe ions(see Nahar et al. 1997; Nahar 1997, 1996 for the results of theR-matrix method). We obtain the rate coefficients of dielectronicrecombination of Fe ions by using HULLAC . We find that our ratecoefficients are higher by a factor of ∼ 𝑀 Ni = . 𝑀 (cid:12) and 𝑣 = − . Figure F2 showsthe spectrum of pure Fe emission at ∼
200 day and ∼
360 day. Alsodepicted is the observed spectrum of a typical SN Ia, SN 2011fe(Mazzali et al. 2015). Our simple one-zone model reproduces thecharacteristic Fe-line structure (see Mazzali et al. 2015 and Botyán-szki et al. 2018 for more detailed modelings). Note that our choiceof 𝑀 Ni = . 𝑀 (cid:12) agrees with the mass estimate by using thepre-nebular light curve of SN 2011fe (Arnett et al. 2017). At ∼ HULLAC , are slightly larger than those used inthe literature.
This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000