© The Authors 2024

1. Introduction

Neutron stars in binary systems are subject to accretion when matter is transferred to them by their companion star. This highly luminous phenomenon is observed especially well in the X-ray band of the electromagnetic spectrum. Accretion can occur through stellar winds or through an accretion disk if the companion star has evolved into a Roche Lobe donor. In the case of accretion by Roche Lobe overflow, accreted material crashes onto the neutron star surface (spinning up the star in the process) when a loss of angular momentum is triggered in the accretion disk.

The composition and the equation of state of the crust of an accreting neutron star is different from that of the crust of an isolated (catalyzed) neutron star. As material originating from the companion star crashes onto the neutron star envelope, thermonuclear burning occurs (Meisel et al. 2018). As a result, ashes with a nucleon number of A = 50 − 110 are deposited on the crust surface and progressively pushed deeper in the star. Contrary to a catalyzed crust that is at global equilibrium, for an accreted crust, the element that falls onto the surface is compressed and subject to non-equilibrium and exothermic reactions. The energy released by electron captures and pycnonuclear fusions deposit heat in the crust: this process is called deep crustal heating (Brown et al. 1998). The heat deposited is transported in the star and radiated through the atmosphere. Then, the resulting luminosity exhibited by the star can be observed after the active accretion process has stopped, that is to say, during quiescence. Consequently, models of accretion related heating of the crust can be compared to the thermal relaxation of X-ray transient sources.

The most common approximation to modeling the crust of an accreting neutron star is the fully accreted crust approximation, where the original crust (i.e., the crust present before any accretion has taken place) is neglected. In this approximation, the crust is solely made of accreted material and the heat sources are always located the same in the crust. Partially accreted crusts are considered when the amount of accreted matter is small enough that the original crust must be studied while it is being compressed under accreted material (more details on this topic are given in Suleiman et al. 2022 and references therein). The usual approach to the crust equation of state and composition is the single nucleus model (Haensel & Zdunik 2008; Fantina et al. 2018; Potekhin & Chabrier 2021). The role of neutron diffusion was investigated recently (Shchechilin et al. 2022), but it was found that it does not significantly affect the heat release of the outer crust (Potekhin et al. 2023).

Electron captures on an even-even parent nucleus (A, Z), with A as the nucleon number and Z as the proton number, occur according to the following pair of reactions:

$$\begin{array}{cc}& (A,Z)+e^{-}\to (A,Z-1)+\nu \text{slow},\hfill \end{array}$$(1)

$$\begin{array}{cc}& (A,Z-1)+e^{-}\to (A,Z-2)+\nu \text{fast},\hfill \end{array}$$(2)

with ν denoting the electronic neutrinos, (A, Z − 1) the daughter nucleus, and (A, Z − 2) the grand-daughter nucleus. The first electron capture is slow, whereas the second electron capture can be considered instantaneous due to the nuclear energy drop when passing from an odd-odd nucleus (A, Z − 1) to an even-even (A, Z − 2) one. Therefore, the kinetics of the pair of reactions is dominated by the first electron capture. The prevailing approach to the kinetics of electron captures in deep crustal heating is that all nuclei in an infinitesimally small piece of accreted matter are instantaneously changed from parent to daughter nuclei once the pressure threshold of the reaction has been reached Sato (1979), Haensel & Zdunik (1990, 2003, 2008), Fantina et al. (2018), Shchechilin & Chugunov (2019), Chugunov & Shchechilin (2020), Shchechilin et al. (2022); in the following, this approximation is referred to as the “instantaneous approach”.

Full reaction networks beyond the one-component approach were discussed in Gupta et al. (2007, 2008), Steiner (2012), Lau et al. (2018), Schatz et al. (2022). Shchechilin et al. (2021) found that complete reaction networks lead to a total heat deposited similar to the heat deposited in the “instantaneous approach”. The impact of time-dependent accretion on the heat release subsequent to the pycnonuclear fusion of ³⁴Ne in the inner crust while considering a finite reaction rate and transient accretion model was studied in Yakovlev et al. (2006). However, the physical processes leading to the reactions are completely different from those considered in this paper, as we study electron captures. Finite reaction rates of electron captures in the crust of accreting neutron stars were discussed in Bildsten & Cumming (1998), Ushomirsky et al. (2000): their study which includes temperature effects, focused solely on stationary accretion and does not include nucleus dependence of the half life in electron capture rates. In this paper, we discuss how including the reaction rate of electron captures with a complete nuclei dependence informed by nuclear experiments affects the composition and the heat release in a fully accreted outer crust constructed within the single nucleus model at zero temperature and subject to time-dependent and transient accretion.

In Section 2, we present the simple approximation taken to model the behavior of a piece of matter with a mixture of parent and grand-daughter nuclei in the outer crust of a neutron star in an accreting binary system. The equation of state followed by the mixture of parent and grand-daughter nuclei is derived in Sect. 2.1. Transitions between the different shells of the outer crust in this approach are established in Sect. 2.2. Section 2.3 details the derivation of electron capture reaction rates and the equation to solve the parent nuclei abundance is presented in Sect. 2.4. Astrophysically motivated cycles of accretion alternating active and quiescent phases, as well as values of accretion rates are discussed in Sect. 2.5. The heat release and neutrino loss related to the parent nuclei abundance is derived in Sect. 2.6 and Sect. 2.7. Results for the parent nucleus abundance are presented in Sect. 3 for the shells of an ⁵⁶Fe ashes accreted outer crust for various scenarios of the accretion cycle. In addition, the results for the heat released as well as the neutrino loss are presented in this section and compared to the results in the instantaneous approach.

2. Methods

2.1. The accreted outer crust equation of state

A simple approach is taken to establish the zero-temperature relation between the pressure and the density (equation of state) in the outer crust of an accreted neutron star: a lattice allowing for the mixture of two nuclei is permeated by a gas of electrons. The pressure and energy density of a degenerate gas of electrons are given by

$$\begin{array}{cc}\hfill P_e& =\frac{{(m_e{c}^2)}^4}{{(ħc)}^3}\varphi (\xi ),\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {\mathcal{E}}_e& =\frac{{(m_e{c}^2)}^4}{{(ħc)}^3}\chi \left(\xi \right),\hfill \end{array}$$(4)

with ℏ as the reduced Planck constant and c as the light velocity such that ℏc = 197.33MeVfm. The electron Fermi momentum is p_Fe = n_e^1/3(3π²)^1/3ℏ and the relativity parameter is ξ = p_Fe/m_ec with n_e and m_e, respectively, the number density and mass of the electron (see e.g., Shapiro & Teukolsky 1986). The expression of the electron’s pressure and energy density involves the dimensionless functions:

$$\begin{array}{cc}\hfill \varphi (\xi )& =\frac{\xi {(1+{\xi }^2)}^{1/2}(\frac{2}{3}{\xi }^2-1)+ln(\xi +{(1+{\xi }^2)}^{1/2})}{8{\pi }^2},\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill \chi (\xi )& =\frac{\xi {(1+{\xi }^2)}^{1/2}(1+2{\xi }^2)-ln(\xi +{(1+{\xi }^2)}^{1/2})}{8{\pi }^2}.\hfill \end{array}$$(6)

In the outer crust, below the liquid ocean, electrons are relativistic: ξ ≫ 1. In the considered reaction layers, the relativity parameter ξ is between 8 (for the first layer), up to ∼40 (for the deepest one). With an accuracy better than ∼1%, one can use ultra-relativistic approximation ϕ(ξ)≃ξ⁴/(12π²).

An ion lattice correction to the pressure, denoted P_lat, is added to the rigid electron background pressure by considering nuclei of number density n_I. The volume per nucleus is 1/n_I and the volume of a sphere centered at a nucleus (the ion sphere) has the radius a_I = (4π/3 n_I)^−1/3. The number density of nuclei n_I is related to the electron density, n_e, by

$$\begin{array}{c}\hfill n_e=n_I(X{Z}_0+(1-X)Z_2)\end{array}$$(7)

for the capture of two electrons, where the parent nucleus proton number is denoted Z₀ and the grand-daughter nucleus proton number is denoted Z₂ = Z₀ − 2. The quantity X = N₀/(N₀ + N₂) designates the nucleus (A, Z₀) abundance, with N₀ the number of parent nuclei (A, Z₀) and N₂ the number of grand-daughter nuclei (A, Z₂).

The lattice pressure P_lat for a mixture of two nuclei is given by

$$\begin{array}{c}\hfill P_{\mathrm{lat}}=-0.3(\frac{4\pi }{3}{)}^{1/3}{n}_e^{4/3}{e}^2\mathcal{F}(X),\end{array}$$(8)

with e the elementary charge and ℱ a function determined by the linear mixing rule (for details, see Section 2.4.7 of Haensel et al. 2007), which for a mixture of parent and grand-daughter ions is given by the formula

$$\begin{array}{c}\hfill \mathcal{F}(X)=\frac{{(Z_0-2)}^{5/3}+X[Z_0^{5/3}-{(Z_0-2)}^{5/3}]}{Z_0+2(X-1)}.\end{array}$$(9)

The relation between the baryon density denoted n and the density of electrons in the linear mixing rule approximation is

$$\begin{array}{c}\hfill n(P,X)=\frac{n_e(P,X)}{\gamma (X)},\end{array}$$(10)

with

$$\begin{array}{c}\hfill \gamma (X)=\frac{X{Z}_0+(1-X)(Z_0-2)}{A}.\end{array}$$(11)

The relation between the pressure and the baryon density in this approximation is given by

$$\begin{array}{c}\hfill P=\frac{{(m_e{c}^2)}^4}{{(ħc)}^3}\varphi \left(\alpha (X)n^{1/3}\right)+\beta (X)n^{4/3},\end{array}$$(12)

with

$$\begin{array}{cc}\hfill \alpha (X)& =\frac{{(3{\pi }^2)}^{1/3}ħc}{m_e{c}^2}(\gamma (X){)}^{1/3},\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill \beta (X)& =-0.3(\frac{4\pi }{3}{)}^{1/3}{e}^2\mathcal{F}(X)(\gamma (X){)}^{4/3}.\hfill \end{array}$$(14)

The chemical potential of a degenerate gas of electrons in a lattice is given by

$$\begin{array}{c}\hfill {\mu }_e=\sqrt{n^{2/3}{[\alpha (X)m_e{c}^2]}^2+{(m_e{c}^2)}^2}+4\frac{\beta (X)}{\gamma (X)}{n}^{1/3}.\end{array}$$(15)

The energy of the nucleus (A, Z) in the heavy nucleus approximation is given by its nuclear mass

$$\begin{array}{c}\hfill \mathcal{M}(A,Z)c^2={\mathcal{M}}_{\mathrm{at}}(A,Z)c^2-Z{m}_e{c}^2+B_{\mathrm{el}}(Z)+E_{\mathrm{exc}}(A,Z),\end{array}$$(16)

with ℳ_at(A, Z) as the atomic mass of the ground state nucleus (A, Z) and B_el(Z) the binding energy of Z electrons, as reported in Lunney et al. (2003). The total energy density is given by the formula

$$\begin{array}{c}\hfill \mathcal{E}=X\frac{\mathcal{M}(A,Z)c^2}{A}n+(1-X)\frac{\mathcal{M}(A,Z-2)c^2}{A}n+{\mathcal{E}}_{\mathrm{lat}}+{\mathcal{E}}_e,\end{array}$$(17)

where the lattice energy is ℰ_lat = 3P_lat.

2.2. Composition of the outer crust

In this paper, the composition of the outer crust is established starting from the fully accreted crust approximation. The crust of the neutron star is entirely composed of accreted matter. The boundaries of shells are defined by the pressure threshold of reactions. We consider the five shells of an accreted outer crust made originally of ⁵⁶Fe ashes from thermonuclear flashes. Nuclei in the entire outer crust have the nucleon number A = 56 because pycnonuclear fusions are not allowed there. The proton numbers from the neutron star crust surface to the outer and inner crust transition are Z = 26,  24,  22,  20,  18, namely: there are four pairs of electron captures that occur during the accretion process onto a fully accreted outer crust.

The composition of the outer crust is determined either by the table of experimentally measured nuclei mass table AME2020 presented in Wang et al. (2021), or by the theoretical model HFB-21 presented in Goriely et al. (2010) if masses of nuclei are not available from experimental data (it is the case for Z = 19, 18). The nuclear masses that we use here, as well as their β-decay schemes, are in vacuum. This stems from the fact that the electron gas does not influence the structure of the nuclei involved, including their excited states.

In the reactions studied in this paper, we only consider the excited state of the daughter nuclei of the first electron capture, with an excitation energy of E_exc(A, Z − 1). Actually, it is only the first electron capture in the first shell in the crust that involves an excited state of the nucleus ⁵⁶Mn, with an excited energy of 110keV, see Appendix B. The pressure threshold and consequently the boundaries of the five shells of the outer crust are defined by the energy threshold, W₁ and W₂, for the first and second electron capture, respectively, and is given by:

$$\begin{array}{cc}\hfill W_1& ={\mathcal{M}}_{\mathrm{at}}(A,Z-1)c^2-{\mathcal{M}}_{\mathrm{at}}(A,Z)c^2+m_e{c}^2\hfill \\ \hfill & +B_{\mathrm{el}}(Z-1)-B_{\mathrm{el}}(Z)+E_{\mathrm{exc}}(A,Z-1),\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill W_2& ={\mathcal{M}}_{\mathrm{at}}(A,Z-2)c^2-{\mathcal{M}}_{\mathrm{at}}(A,Z-1)c^2+m_e{c}^2\hfill \\ \hfill & +B_{\mathrm{el}}(Z-2)-B_{\mathrm{el}}(Z-1).\hfill \end{array}$$(19)

The atomic mass ℳ_at(A, Z) is extracted from experimental data for nuclei Z = [26 − 20] and calculated, using a theoretical approach for Z = 19, 18. The pressure and density threshold denoted P_th and n_th, respectively, defining the transition between each shell is determined by solving μ_e(X = 1, n_th) = W. In Table 1, we present the transition pressures and densities, as well as the energy threshold of the first and second electron captures, respectively, denoted W₁ and W₂. Our approach is different than the single-nucleus model considered in Fantina et al. (2018) in two respects. On the one hand, experimentally measured nuclei masses used in Fantina et al. (2018) are taken from the AME2016 Wang et al. (2017) data table and we used the updated AME2020 data table in our calculations; this is the main source of discrepancy (for the first transition pressure P_th(Z = 26), including the excitation energy plays a significant role). On the other hand, we do not require a simultaneous transition of all nuclei (A, Z)→(A, Z − 1) and, as a result, the threshold pressure in our model is slightly different (for different treatment of threshold pressure for two models, see Chamel et al. (2015), Chamel & Fantina (2016)).

2.3. Rates of β-decays and electron captures

The reaction rate of β-decays and electron captures can be established from the Fermi golden rule (originally presented and applied by Dirac 1927). For details on the derivation of the reaction rates, see Appendix A. For the β-decay rate (A, Z − 1)→(A, Z) denoted ℛ_β the calculation is done in vacuum and exclusively for allowed type of decays, resulting in the following formula:

$$\begin{array}{c}\hfill {\mathcal{R}}_{\beta }=\frac{2J_Z+1}{2J_{Z-1}+1}{|{\mathcal{M}}_{\beta }|}^2\frac{4m_e^5{c}^4}{{(2\pi )}^3{ħ}^7}f,\end{array}$$(20)

with J_Z as the spin of nucleus (A, Z) and here the factor f is the available momentum space factor; notably, the matrix element, ℳ_β, does not depend on the nuclei spins.

The electron capture rate on nuclei immersed in the Fermi sea of electrons (A, Z)→(A, Z − 1) is denoted by ℛ_ec and can then be expressed as:

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathrm{ec}}=\frac{2J_{Z-1}+1}{2(2J_Z+1)}{|{\mathcal{M}}_{\mathrm{ec}}|}^2\frac{4m_e^5{c}^4}{{(2\pi )}^3{ħ}^7}\mathcal{G},\end{array}$$(21)

where 𝒢 is a dimensionless available momentum space factor, given by

$$\begin{array}{cc}\hfill \mathcal{G}({\overline{E}}_F,\overline{W})& ={\int }_{\overline{W}}^{{\overline{E}}_F}\overline{E}\sqrt{{\overline{E}}^2-1}{(\overline{E}-\overline{W})}^2\mathrm{d}\overline{E},\hfill \\ \hfill & =[{\mathcal{G}}_b({\overline{E}}_F,\overline{W})-{\mathcal{G}}_b(\overline{W},\overline{W})],\hfill \end{array}$$(22)

with ${\overline{E}}_F=E_F/(m_e{c}^2)$, $\overline{W}=W/(m_e{c}^2)$ and E_F the Fermi energy of electrons. For the β-decay reactions, the factor f is $\propto {\mathcal{G}}_b({\overline{W}}_1,{\overline{W}}_1)-{\mathcal{G}}_b(1,{\overline{W}}_1)$ and often includes corrections for the Coulomb interaction. The generalized function, 𝒢_b, can be obtained via an analytical formula, first presented in Frank-Kamenetskii (1962)¹ as:

$$\begin{array}{cc}\hfill {\mathcal{G}}_b(x,w)& =\frac{\sqrt{x^2-1}}{60}(x^3(12x-30w)-4x^2(1-5w^2)\hfill \\ & +15wx-20w^2-8)+\frac{1}{4}wlog(\sqrt{x^2-1}+x).\hfill \end{array}$$(23)

In the ultra-relativistic approach (${\overline{E}}_F\ge \overline{W}\gg 1$), the expansion up to the first order in ${\overline{W}}^{-2}$ yields the approximation 𝒢 ≈ 𝒢_ultra:

$$\begin{array}{cc}\hfill {\mathcal{G}}_{\mathrm{ultra}}(E_F,W)& =\frac{{\overline{W}}^5}{3}(\frac{E_F}{W}-1{)}^3[1-\frac{1}{2{\overline{W}}^2}+\frac{3}{2}(\frac{E_F}{W}-1)\hfill \\ \hfill & +\frac{3}{5}(\frac{E_F}{W}-1{)}^2].\hfill \end{array}$$(24)

For an electron Fermi energy close to the reaction threshold energy (E_F − W ≪ W), the leading order expression is:

$$\begin{array}{c}\hfill {\mathcal{G}}_{\mathrm{ultra}}(E_F,W)=\frac{{\overline{W}}^5}{3}(\frac{E_F}{W}-1{)}^3,\end{array}$$(25)

here, we also refer to Eq. (4) and the zero temperature limit of Eq. (6) in Bildsten & Cumming (1998). We note that the ultra-relativistic approximation is used in this paper only to make analytical estimates, the full numerical solution makes use of Eq. (23) to formulate the reaction rate.

Using charge conjugation symmetry as well as time reversal symmetry, we can show that the matrix element for the electron capture and its corresponding β-decay in the above mentioned approximations are equal. Therefore, we can make use of the experimental data available for β-decay reactions to compute the reaction rate of the first electron capture per parent nucleus:

$$\begin{array}{cc}\hfill {\mathcal{R}}_{\mathrm{ec}}& =\frac{ln(2)}{f{t}_{1/2}}\frac{2J_{\mathrm{ec};\mathrm{Z}-1}+1}{2(2J_{\mathrm{ec};\mathrm{Z}}+1)}\frac{2J_{\beta ;Z-1}+1}{2J_{\beta ;Z}+1}\mathcal{G}({\overline{E}}_F,{\overline{W}}_1),\hfill \end{array}$$(26)

$$\begin{array}{cc}& \equiv \frac{1}{{\tau }_{\mathrm{ec}}}\mathcal{G}({\overline{E}}_F,{\overline{W}}_1),\hfill \end{array}$$(27)

where J_ec; Z and J_β; Z are nuclear spins for the electron capture and its corresponding β-decay, respectively. The quantity t_1/2 is the half-life of the β-decay for a specific reaction channel and can be extracted from experimental data or estimated analytically, τ_ec is the electron capture timescale. Electron captures and β-decays are subject to selection rules (Povh et al. 2004) that specify allowed changes of nuclear spins for nuclei involved in electron captures and β-decays. We selected the dominant channel in the energy level diagrams presented in National Nuclear Data Center (2022) for each reaction. In Appendix B, we give some details on the calculation of the electron capture rates. In Table 2, we present values of τ_ec for the first electron capture involved in the four shells of the outer crust.

2.4. Continuity equation

During the active phases of accretion, freshly accreted matter is flowing and sinking towards the core with a velocity:

$$\begin{array}{c}\hfill v(z)=\frac{\dot{M}}{4\pi R^2{m}_{B}n(z)},\end{array}$$(28)

with Ṁ the mass accretion rate, R the total radius of the star, m_B the rest mass of a baryon, and z = R − r the proper distance from the surface in the plane-parallel approximation. It should be noted that even if the accretion stops (quiescent phase), nuclear reactions still occur leading to a slow change of the equation of state in the mixed layer. This results in a shrinking of the reaction layer, and therefore to a certain (low) velocity. In our considerations we use the pressure as an independent variable (in place of z or r). For a given piece of matter, the pressure changes due to accretion with the increase proportional to the accreted mass; during quiescent phases, the pressure does not change as the column of matter above is fixed when the reactions take place.

The most important quantity introduced in the previous section is the parent nucleus abundance, X. The layers of electron captures are defined as the thickness in pressure, for which X ∈ [0 − 1]. The product of the parent nucleus abundance, X, with the baryon density corresponds to the number density of parent nuclei, and is governed by the continuity equation. This equation is written in the local frame with τ and z corresponding to the proper time and proper distance in the neutron star crust. In what follows, we use the Newtonian approach

$$\begin{array}{c}\hfill \frac{\partial P}{\partial r}=-\frac{\mathit{GM}}{R^2}\rho ,\end{array}$$(29)

where ρ = ϵ/c² ≈ nm_B, and M is the total mass of the star; in the following the mass and the radius of the star are chosen to be 1.4M_⊙ (solar mass) and 11km, respectively. The continuity equation is therefore given by:

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}ln(n(X(P,t),P)X(P,t))& +\frac{1}{{\tau }_{\mathrm{acc}}(t)}\frac{\partial lnX(P,t)}{\partial \tilde{P}}\hfill \\ \hfill & =-{\mathcal{R}}_{\mathrm{ec}}(X(P,t),P),\hfill \end{array}$$(30)

with $\stackrel{\sim }{P}=P/P_{\mathrm{th}}$, with n = n(P, X) and X = X(P, t), and τ_acc is the accretion timescale,

$$\begin{array}{c}\hfill {\tau }_{\mathrm{acc}}(t)=\frac{4\pi R^4{P}_{\mathrm{th}}}{GM\dot{M}(t)}.\end{array}$$(31)

Typical values of the accretion timescale are presented in Table 2. In the general relativity approach, with a metric ds² = e^2ϕdt² − e^2λdr² − dΩ², the continuity equation has the same form as in the Newtonian one in the local frame with coordinates (te^ϕ, re^λ), and the boundary condition at the surface is e^ϕ(R)+λ(R) = 1. By choosing the pressure, P, as the independent variable instead of z, the relativistic equivalent of Eq. (30) would include Ṁ multiplied by a factor e^ϕ(R) = (1 − 2GM/Rc²)^1/2 in the expression of the accretion timescale (Eq. (31)); in our case (M = 1.4M_⊙, R = 11km), this factor is ∼0.8. The mass accretion rate would be defined as $\dot{M}=\dot{N}/m_B$, with $\dot{N}=dN/dt$, as measured by a distant observer. In the rest of the paper, we neglect relativistic effects and treat the Newtonian approach of the continuity equation.

2.5. Astrophysically motivated accretion cycle

We intend to study a realistic cycle of accretion (sequence of active and quiescent phases) motivated by X-ray observations of sources in accreting low-mass X-ray binaries. Several sources have been observed alternating between active accretion stages lasting from weeks to years and quiescent stages lasting from months to decades. The information on characteristic duration of active accretion can be extracted from the observations of X-ray outbursts (a lasting surge in luminosity) of sources EXO 0748–676 (see Parikh et al. 2020 and references therein), KS 1731–260 (see Merritt 2017 and references therein), XTE J1701–462 (Fridriksson et al. 2010), and IGR J17480–2446 (see Ootes et al. 2019 and references therein). This lasted 24 years, 12.5 years, 1.6 years, and 10 weeks, respectively; all four sources are now in quiescence. The source MXB 1659–29 (see Parikh et al. 2019 and references therein) presented two well monitored outbursts lasting, respectively, 2.5 and 1.7 years, interspersed by a quiescent period of 14 years. This source was first observed in 1976 during an outburst estimated to last between 2 and 2.5 years as well. Unless otherwise stipulated, in this paper we study the following cycle: an active phase period, t_a, of four years and a quiescent period of t_q = 10t_a.

The average accretion rate during active phases were estimated for the sources XTE J1701–462 (Fridriksson et al. 2010) and IGR J17480–2446 (Degenaar & Wijnands 2011) to be 1.7 × 10⁻⁸ M_⊙ per year and 3 × 10⁻⁹ M_⊙ per year, respectively. The observation of X-ray bursts indicate a luminosity variability in the active phase from which the accretion rate as a function of time can be inferred, however, we do not intend to discuss this variability; rather, we aim to use a mean value of the accretion rate in the active phase. The active accretion rate, denoted Ṁ_a, is described by an exponential onset and offset such that:

$$\begin{array}{c}\hfill {\dot{M}}_a(t)={\dot{M}}_{\max }(\frac{1}{1+e^{a(t-t_{\mathrm{a}}+t_{\mathrm{o}})}}+\frac{1}{1+e^{-a(t-t_{\mathrm{o}})}}-1),\end{array}$$(32)

with Ṁ_max the accretion rate outside the onset and offset of active accretion and a a constant controlling the steepness of the Ṁ increase and decrease during onset and offset. For a ≃ 118 (in yr⁻¹), the onset and offset are very steep. We consider the accretion is switched on as soon as Ṁ reaches ${10}^{-10}{\dot{M}}_{\max }$. At t_o, half of the maximum value of the accretion rate ${\dot{M}}_{\max }/2$ is reached. After t_o, because a is so large and the onset so steep, we quickly reach Ṁ_max to a very high precision. Unless otherwise stipulated, in this paper we set ${\dot{M}}_{\max }={10}^{-8}$M_⊙ per year, a reasonable value with regard to observed mean accretion rate during active accretion, and t_o = 0.05t_a.

The number of accreted baryons during an accretion cycle denoted N_b is defined by the time integral of the accretion rate during the cycle

$$\begin{array}{c}\hfill N_b={\int }_0^{t_a+t_q}\frac{{\dot{M}}_a(t)}{m_B}\mathrm{d}t.\end{array}$$(33)

Given that we consider a finite reaction rate, the number of accreted nuclei over one accretion cycle is not a priori the same as the number of nuclei undergoing electron captures. Therefore, we define ${\dot{N}}_r$ the rate of baryons involved in an electron capture per unit time and per unit volume

$$\begin{array}{c}\hfill {\dot{N}}_r(P,t)={\mathcal{R}}_{\mathrm{ec}}(P,X)X(P,t).\end{array}$$(34)

The number of reactive baryons per unit time in the shell gives

$$\begin{array}{c}\hfill N_r(t)=\frac{4\pi R^4{P}_{\mathrm{th}}}{GM{m}_B}{\int }_1^{{\tilde{P}}_0}{\mathcal{R}}_{\mathrm{ec}}(\tilde{P},X)X(\tilde{P},t)\mathrm{d}\tilde{P},\end{array}$$(35)

where ${\stackrel{\sim }{P}}_0=P(X=0)/P_{\mathrm{th}}$. The number of reactive baryons in the shell over the cycle of accretion

$$\begin{array}{c}\hfill N_r={\int }_0^{t_a+t_q}{N}_r(t)\mathrm{d}t.\end{array}$$(36)

2.6. Heat release from electron captures

We consider the first electron capture presented in Eq. (1), while neglecting the neutrino heat loss (cooling). We also consider that processes under consideration proceed at constant pressure and temperature in a piece of matter. Following the derivation of Chapter X in Landau & Lifshitz (1969), we define an elementary reaction e⁻ + (A, Z)→(A, Z − 1), with the coefficients ν_e = ν_A, Z = 1 and ν_{A, Z − 1} = −1. The number of electron captures in the time interval dt is denoted δN_ec; the changes of numbers of particles i are δN_i = −ν_iδN_ec, while the thermodynamic potential Φ changes by δΦ = ∑_iδN_iμ_i = −∑_iν_iμ_iδN_ec. The heat released at constant pressure and temperature is (Chapter X in Landau & Lifshitz (1969)):

$$\begin{array}{c}\hfill \delta Q_P=T\delta S=-T^2\left(\frac{\partial }{\partial T}\frac{\delta \mathrm{\Phi }}{T}\right),\end{array}$$(37)

with T the temperature and S entropy. In our approximation, μ_i values do not depend on T (electrons are fully degenerate and nuclear masses do not depend on T), so that

$$\begin{array}{c}\hfill \delta Q_P=-\sum _i{\nu }_i{\mu }_i\delta N_{\mathrm{ec}}.\end{array}$$(38)

The heating rate is therefore

$$\begin{array}{c}\hfill {\dot{Q}}_P=T\dot{S}=-\sum _i{\nu }_i{\mu }_i{\dot{N}}_{\mathrm{ec}}.\end{array}$$(39)

The heat release per nucleus for the first and second electron capture of the pair, denoted respectively q₁ and q₂, neglecting neutrino heat loss, are

$$\begin{array}{cc}\hfill q_1(P,X)& =\mathcal{M}(A,Z)c^2-\mathcal{M}(A,Z-1)c^2+{\mu }_e(P,X),\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill q_2(P,X)& =\mathcal{M}(A,Z-1)c^2-\mathcal{M}(A,Z-2)c^2+{\mu }_e(P,X).\hfill \end{array}$$(41)

The heat released by the pair of electron captures per nucleus is denoted q and is given by

$$\begin{array}{c}\hfill q(P,X)=2{\mu }_e(P,X)-(W_1+W_2)+E_{\mathrm{exc}},\end{array}$$(42)

with E_exc the excitation energy of the first electron capture daughter nucleus; note that except for ⁵⁶Mn^*, we only consider ground-state to ground-state reactions. Therefore, the heating rate per unit volume is given by

$$\begin{array}{c}\hfill \dot{q}(P,t)=q(P,X){\mathcal{R}}_{\mathrm{ec}}(P,X)\frac{n(P,X)}{A}X(P,t).\end{array}$$(43)

The heat released in the layer of electron captures per unit time is

$$\begin{array}{cc}\hfill \dot{\mathcal{Q}}(t)& =\frac{4\pi R^4{P}_{\mathrm{th}}}{GM{m}_B}{\int }_1^{{\tilde{P}}_0}\frac{\dot{q}(\tilde{P},t)}{n(\tilde{P},X)}\mathrm{d}\tilde{P}.\hfill \end{array}$$(44)

The total heat released in a layer per one accreted nucleon during one accretion cycle is

$$\begin{array}{c}\hfill Q=\frac{1}{N_b}{\int }_0^{t_a+t_q}\dot{\mathcal{Q}}(t)\mathrm{d}t,\end{array}$$(45)

with N_b the number of accreted nucleons during this cycle.

2.7. Neutrino energy loss

In the assessment of the heat release in the outer crust of a neutron star in an accreting binary system, the neutrinos must be taken into consideration: in conditions relevant for adult neutron stars, neutrinos escape the system taking away energy.

Similarly to Eq. (42), we can define the energy allocated to neutrinos in the pair of electron captures per nucleus

$$\begin{array}{cc}\hfill q_{\nu }(P,X)& =m_e{c}^2(\frac{{\mathcal{G}}_{\nu }({\overline{E}}_F,{\overline{W}}_1)}{\mathcal{G}({\overline{E}}_F,{\overline{W}}_1)}+\frac{{\mathcal{G}}_{\nu }({\overline{E}}_F,{\overline{W}}_2)}{\mathcal{G}({\overline{E}}_F,{\overline{W}}_2)}).\hfill \end{array}$$(46)

The energy of neutrinos per electron capture and per parent nucleus, in the units of m_ec², denoted 𝒢_ν, is given by

$$\begin{array}{cc}\hfill {\mathcal{G}}_{\nu }({\overline{E}}_F,\overline{W})& ={\int }_{\overline{W}}^{{\overline{E}}_F}\overline{E}\sqrt{{\overline{E}}^2-1}{(\overline{E}-\overline{W})}^3\mathrm{d}\overline{E}\hfill \\ \hfill & =[{\mathcal{G}}_{\nu b}({\overline{E}}_F,\overline{W})-{\mathcal{G}}_{\nu b}(\overline{W},\overline{W})];\hfill \end{array}$$(47)

the function 𝒢_νb can be obtained analytically as

$$\begin{array}{cc}\hfill {\mathcal{G}}_{\nu b}(x,w)& =\frac{1}{240}(\sqrt{x^2-1}[40x^5-144w{x}^4\hfill \\ \hfill & +10x^3(18w^2-1)+16w{x}^2(3-5w^2)\hfill \\ \hfill & -15x(1+6w^2)+16w(6+5w^2)]\hfill \\ \hfill & -15(1+6w^2)[log(\sqrt{x^2-1}+x)]),\hfill \end{array}$$(48)

see also Bisnovatyi-Kogan (2001).

The energy radiated away by neutrinos per unit time and per unit volume ε_ν is given by

$$\begin{array}{c}\hfill {\epsilon }_{\nu }(P,t)=q_{\nu }(P,X){\mathcal{R}}_{\mathrm{ec}}(P,X)\frac{n}{A}X(P,t).\end{array}$$(49)

The heat deposited in the layer of electron capture per time unit is defined as the integral over the volume of the difference between the heating rate and the cooling rate due to neutrino emission

$$\begin{array}{cc}\hfill \dot{\mathcal{H}}(t)& =\frac{4\pi R^4{P}_{\mathrm{th}}}{GM{m}_B}{\int }_1^{{\tilde{P}}_0}\frac{\dot{q}(\tilde{P},t)-{\epsilon }_{\nu }(\tilde{P},t)}{n(\tilde{P},X)}\mathrm{d}\tilde{P}.\hfill \end{array}$$(50)

The total heat deposited per one accreted baryon during one accretion cycle gives

$$\begin{array}{c}\hfill H=\frac{1}{N_b}{\int }_0^{t_a+t_q}\dot{\mathcal{H}}(t)\mathrm{d}t.\end{array}$$(51)

2.8. Heating from (A, Z − 1)^* decay versus neutrino cooling

The daughter nucleus (A, Z − 1) from the first electron capture in Eq. (1) has odd numbers of neutrons and protons. Therefore, it has a lower energy threshold for the second electron capture, W₂ < W₁. Moreover, it has numerous excited states above the ground state one. This means that there is an energy range of electrons W₂ < E_e < E_F, whose capture leads to the excited states (A, Z − 1)^*. These excited states with energies 0 < E_exc < E_e − W₂ de-excite via electromagnetic transition, heating the matter. The multi-component plasma calculations of the accreting crust evolution show that the de-excitation heating balances the neutrino heat losses calculated in Sect. 2.7 (Gupta et al. 2007). Therefore, neglecting neutrino losses and using the heating rate calculated in Sect. 2.6 is a reasonable approximation. Such an approximation was used in (Haensel & Zdunik 2008; Fantina et al. 2018).

3. Results

3.1. Thickness of the layer of electron capture

3.1.1. Stationary solution

As a first approach to the behavior of the parent nucleus abundance, X, in each shell of the outer crust, we establish the solution of the continuity equation for a constant accretion rate; the solution depends only on the pressure. Equation (30) is reduced to an ordinary differential equation treated with a Runge-Kutta numerical approach, with a constant accretion rate of $\dot{M}={\dot{M}}_{\max }={10}^{-8}{\mathrm{M}}_{\odot }$ per year. The boundary condition X(P = P_th) = 1 is used: at the pressure threshold, all nuclei are parent nuclei. The accreted material is continuously pushed deeper in the crust by freshly accreted material and the solution does not depend on time.

When a piece of matter composed of parent nuclei reaches the pressure threshold of the first electron capture, parent nuclei are allowed to transform into daughter (and immediately grand-daughter) nuclei. In the instantaneous approximation, all parent nuclei react at the threshold pressure of the first electron capture. However, when considering finite reaction rates, parent nuclei electron captures take place according to the finite reaction rate and as the piece of matter is still pushed deeper in the crust. Contrary to the instantaneous approximation, in our approach (later referred to as the mixed layer approach), the sinking piece of matter still contains parent nuclei that shall be changed to grand-daughter nuclei at values of the pressure higher than the threshold pressure. The stationary solution corresponds to a particular distribution of the parent nuclei number as a function of the pressure: it appears during active phases of accretion if a steady state has been reached between a continuous supply of parent nuclei (and compression) provided by accretion and the reactions operating on the parent nuclei.

The stationary solution X_stat(P) is presented for parent proton number Z₀ = 26,  24,  22,  20, in Fig. 1. We define the thickness in pressure of the layer of electron capture, that is to say, the region in the shell where both parent and grand-daughter nuclei are found, using the value of the pressure at which the shell is made of only 1% of parent nuclei. It is denoted ΔP_layer = P(X = 10⁻²)−P_th; this quantity is also presented in Table 3 for the stationary solution. For all four pairs of electron captures, the thickness of the layer is on the order of 10⁻⁷MeV/fm³. The slope of decrease in the parent abundance for the stationary solution is correlated with the ratio τ_acc/τ_ec (see Table 2): the higher the ratio, the steeper the decrease.

Fig. 1. Stationary solution for the parent nucleus abundance, X, as a function of the pressure for the four reactions in the outer crust.

It is important to note that the thickness of the electron capture layer is much smaller than the thickness of the shell composed of the grand-daughter nuclei. The thickness in pressure for all considered shells is ∼10⁻⁵ − 10⁻³MeV/fm³; the reaction layer thickness represents at most one percent of the shell thickness, as said shells are located deeper in the crust and the layer thickness is on the same order of magnitude throughout the outer crust. This implies that no shell of the outer crust with parent number nucleus, Z′, is also composed of its shallower neighbor parent nucleus, Z′−4.

Using the ultra-relativistic approximation of the electron capture reaction rate presented in Eq. (24), we can provide a simple analytical approximation of the increase in pressure δP required to reach the region in which most of the parent nuclei are transformed into grand-daughter nuclei. The stationary form of the continuity equation for ultra-relativistic electrons is written as

$$\begin{array}{c}\hfill \frac{{\tau }_{\mathrm{ec}}}{{\tau }_{\mathrm{acc}}}{P}_{\mathrm{th}}\frac{\partial ln(X)}{\partial (\delta P)}=-\frac{{\overline{W}}_1^5}{3}(\frac{E_F}{W_1}-1{)}^3.\end{array}$$(52)

In this approximation, we include only the leading term in Eq. (24) – it is valid for E_F/W₁ − 1 ≪ 1, equivalent to the condition δP/P ≪ 1, (see Eq. (25)). This condition is fulfilled for all considered reaction layers but the first one, for which it is on the order of 0.1–0.2 (see Tables 1, 3). Under the assumption that the parent nuclei abundance is close to one, the factor α (see Eq. (13)) related to the linear mixing rule function is considered to be constant; thus, the right hand side of Eq. (52) is independent of X. Solving the approximation of the continuity equation in the ultra-relativistic approach (P ∼ E_F⁴ and δP/P ≃ 4(E_F/W − 1)), we obtain the approximate solution to Eq. (52) X(δP) = exp[ − (δP/δP_s)⁴], where δP_s is e-folding value of δP at which X = 1/e defined by:

$$\begin{array}{c}\hfill \delta P_{\mathrm{s}}=4P_{\mathrm{th}}(\frac{3}{{\overline{W}}_1^5}\frac{{\tau }_{\mathrm{ec}}}{{\tau }_{\mathrm{acc}}}{)}^{1/4}.\end{array}$$(53)

This approximation works well for the a small decrease of X: for X = 0.95, it can be quite accurately determined by

$$\begin{array}{c}\hfill \delta P(X=0.95)=0.476\delta P_{\mathrm{s}}=2.51P_{\mathrm{th}}(\frac{1}{{\overline{W}}_1^5}\frac{{\tau }_{\mathrm{ec}}}{{\tau }_{\mathrm{acc}}}{)}^{1/4},\end{array}$$(54)

a quantity presented in Table 3.

Similarly to the estimation of the pressure thickness in Eq. (53), we can estimate the radius thickness of the electron capture layers without solving the continuity equation as:

$$\begin{array}{c}\hfill \delta R=4\frac{R^2{P}_{\mathrm{th}}}{GM{\rho }_{\mathrm{th}}}(\frac{3}{{\overline{W}}_1^5}\frac{{\tau }_{\mathrm{ec}}}{{\tau }_{\mathrm{acc}}}{)}^{1/4},\end{array}$$(55)

with ρ_th = n_thm_B. For a neutron star mass of M = 1.4 M_⊙ and radius R = 11km: the layers for Z₀ = 26, Z₀ = 24, Z₀ = 22, and Z₀ = 20 are about 6m, 0.4m, 0.2m, and 0.2m thick, respectively. This is much less than the thickness of the shells which in all studied cases are larger than 100m (see Eq. (10) in Suleiman et al. 2022). The maximum size of the reaction layer (corresponding to stationary solution) is significantly smaller than the thickness of the shell also for other, astrophysically relevant accretion rates. The presented numbers correspond to $\dot{M}={10}^{-8}{\mathrm{M}}_{\odot }/\mathrm{yr},$ but even a decrease of Ṁ by two orders of magnitude only leads to a factor of ∼0.3 in Eqs. (53)–(55). For a small (5%) decrease in X, we get:

$$\begin{array}{c}\hfill \delta R(X=0.95)=2.5\frac{R^2{P}_{\mathrm{th}}}{GM{\rho }_{\mathrm{th}}}(\frac{1}{{\overline{W}}_1^5}\frac{{\tau }_{\mathrm{ec}}}{{\tau }_{\mathrm{acc}}}{)}^{1/4},\end{array}$$(56)

in the case of which we obtained about 3m for the first reaction layer and ∼0.1 − 0.2 m for deeper reaction layers².

3.1.2. Time-dependent active phase solution for Z₀ = 26

To solve the complete continuity equation in Eq. (30) during active accretion, we used a time-dependent accretion timescale with ${\dot{M}}_{\max }={10}^{-8}{\mathrm{M}}_{\odot }$. The following accretion cycle parameters were used in Eq. (32): t_a = 4years and t_o = 0.2 years. The boundary condition is the same as the stationary case and the initial condition is X(P, t = 0) = X_stat(P). This simplistic approach is reasonable if the equation is solved for multiple accretion cycles and we present results for the fourth accretion cycle for the shell with Z₀ = 26. It turned out that after two cycles, we reached a stable solution, in which X(P, t) = X(P, t + τ_cycle), where τ_cycle = t_a + t_q is the full period of our accreting source. After reaching this stable solution all cycles look identical.

Results for the parent nucleus abundance X(P, t) for Z₀ = 26 are presented in Fig. 2. This figure shows that the pre-accretion profile X(P, t − t₀ = 0) (shown in black) is not step-like: the quiescent time from the previous cycle is not sufficient to reach a very small abundance of parent nuclei in the whole layer. After the accretion is switched on, the parent abundance profile increases (in dark blue) and reaches the stationary profile (in orange). At t − t₀ = 3.8years (dashed brown), the offset has already started and we have ${\dot{M}}_a(t-t_0=3.8)={\dot{M}}_{\max }/2$: we observe that a decrease by half of the accretion rate barely changes the profile of X(P). However, at t − t₀ = 3.9years (in teal), for which the accretion rate can be considered negligible with ${\dot{M}}_a(t-t_0=3.9)/{\dot{M}}_{\max }\simeq {10}^{-5}$, X(P) has significantly evolved from the stationary profile.

Fig. 2. Time-dependent solution, X(P, t), represented as snapshot at various times during the four years of active accretion cycle, starting at t0. (a) ${\dot{M}}_{\max }={10}^{-8}$M⊙ per year. (b) ${\dot{M}}_{\max }={10}^{-10}$M⊙ per year.

The time needed to reach the stationary profile for a given reaction layer can be estimated from the formula in Eq. (31) with P_th replaced by δP defined by Eq. (53):

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{stat}}& ={\tau }_{\mathrm{acc}}\frac{\delta P}{P_{\mathrm{th}}}.\hfill \end{array}$$(57)

The formula is most accurate close to the threshold; therefore, we can defined the time needed to reach the stationary profile at X = 0.95 and use the following formula:

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{stat}}(X=0.95)& =2.5(\frac{1}{{\overline{W}}_1^5}{\tau }_{\mathrm{ec}}{\tau }_{\mathrm{acc}}^3{)}^{1/4}\hfill \\ \hfill & =0.4\mathrm{yr}{R}_{10}^3{\tau }_{\mathrm{ec},5}^{1/4}{W}_{10}^{-5/4}{\left(\frac{P_{\mathrm{th},-5}}{{\dot{M}}_{-8}M/M_{\odot }}\right)}^{3/4}\hfill \end{array}$$(58)

where R₁₀ = R/10km, τ_ec, 5 = τ_ec/10⁵s, W₁₀ = W/10MeV, P_th, −5 = P_th/10⁻⁵Mev/fm³, and ${\dot{M}}_{-8}=\dot{M}/{10}^{-8}{M}_{\odot }/\mathrm{yr}$. The timescale needed to reach the stationary state is on the order of the active phase (accretion) in low mass X-ray binaries. To decide whether this state is reached for a given reaction layer, one needs to study the full solution of the continuity equation.

To assess the role of the value of the maximum accretion rate on our results, we present in Fig. 2 the parent nucleus abundance X(P, t) for ${\dot{M}}_{\max }={10}^{-10}{\mathrm{M}}_{\odot }$ per year instead of ${10}^{-8}{\mathrm{M}}_{\odot }$ per year; the same accretion cycle parameters are used. The stationary profile for ${\dot{M}}_{\max }={10}^{-10}{\mathrm{M}}_{\odot }$ is different from the stationary profile for ${\dot{M}}_{\max }={10}^{-8}{\mathrm{M}}_{\odot }$ presented in Fig. 1. The solution X(P, t) does not reach the stationary profile before the accretion is switched off. The profile of X(P, t) evolves slower than for the value of the accretion rate $\dot{M}={10}^{-8}$M_⊙ per year; the layer of electron captures is also much thinner. The accretion rate is directly related to the number of parent nuclei deposited onto the neutron star surface, and the speed at which nuclei are pushed to higher pressures. Therefore, this quantity dictates how fast the solution X(P, t) increases towards the stationary profile.

The time of active accretion t_a also plays a role on the thickness of the electron capture layer: from Fig. 2 with an accretion rate $\dot{M}={10}^{-8}$M_⊙ per year, we can infer that if the accretion is switched off before or around 0.6years (blue curve in the figure), the stationary profile is not reached during the active phase period; in this case, the stationary state is reached after ∼1 yr.

3.1.3. Quiescence following a stationary profile

We now consider the particular case of quiescence ($\dot{M}(t)=0$) for the continuity equation Eq. (30):

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{dt}}ln(n(X(P,t),P)X(P,t))=-{\mathcal{R}}_{\mathrm{ec}}(X(P,t),P),\end{array}$$(59)

which follows a stationary profile, namely, with the initial condition X(P, t = 0) = X_stat(P). Results for 40years of quiescence are presented in the movie available online.

In the mixed layer approach, X(P, t) in the quiescent phase is not zero in the shell, contrary to the instantaneous approach: electron captures still occur even when the accretion-related pressure increase is switched off. All X(P) profiles evolve very rapidly during the first ∼10years after the beginning of quiescence. However, the profiles evolve increasingly slowly (similarly to an exponential decrease) in the remaining years of the quiescent phase towards a step-like function that is not actually reached.

3.1.4. Accreted versus reactive baryons

The number of accreted baryons N_b (see Eq. (33)) and the number of baryons effectively involved in electron captures (reactive baryons) N_r (see Eq. (36)), are equal in the instantaneous approach, but not in the mixed layer approach. For the mixed layer approach, when considering only the stationary case of active accretion, we also have N_b = N_r: the stationary solution represents an equilibrium between the number of parent nuclei pushed beyond the threshold pressure of the reaction by continuous accretion and the number of parent nuclei transformed into daughter nuclei when they react deeper in the shell.

When considering the time-dependent solution X(P, t) during the active phase, not all parent nuclei accreted necessarily have time to undergo their electron capture. For an active phase of fouryears with ${\dot{M}}_{\max }={10}^{-8}$M_⊙ per year, the number of accreted baryons is N_b = 4.3 × 10⁴⁹. During the same period of active accretion, the time-dependent solution for Z₀ = 26 leads to N_r ≃ 3.9 × 10⁴⁹: this corresponds to ∼90% of the baryons accreted during the active phase having reacted.

During the quiescent phase, in the assumption that it follows from a stationary profile, the number of reactive baryons during quiescence is directly related to the speed of evolution of the X(P) profile: we have N_r(Z₀ = 26)≃6.1 × 10⁴⁸, N_r(Z₀ = 24)≃4.2 × 10⁴⁸, N_r(Z₀ = 22)≃5.5 × 10⁴⁸ and N_r(Z₀ = 20)≃2.1 × 10⁴⁹. Those number represent 14%, 10%, 13%, and 48% of the number of accreted baryons respectively for Z₀ = 26, 24, 22, and 20. For Z₀ = 20, the number of baryons that react during quiescence is not negligible compared to the number of reactive baryons during fouryears of the active phase in the stationary profile.

By adding the number of reactive baryons in the active phase and in quiescence for the full time-dependent solution for Z₀ = 26, we obtain N_r^cycle = N_b although a step-like function is not reached at the end of quiescence for the profile X(P, t − t_q). Let us denote N_nr as the number of nuclei that have not reacted at the end of the quiescent phase. At the end of cycle i, we have $N_{\mathrm{nr}}^i$ parent nuclei that have not reacted in a small layer close to the reaction threshold (black line in Fig. 2a) of size P_nrⁱ − P_th. The value of P_nrⁱ and the number of nuclei that have not reacted during this cycle (proportional to P_nrⁱ − P_th³ can be estimated comparing reaction the time (1/ℛ_ec) for P_nrⁱ with the quiescence timescale, τ_q. For example, given the solution presented in Fig. 2b, the value of τ_qℛ_ec(P) is on the order of 1 for $\overline{P}\sim 1.03$.

3.2. Heat deposition in layers of electron capture

3.2.1. Heat release in active phase: Stationary solution

During active phases of accretion, the amount of heat released in the mixed layer approach and in the instantaneous approach are different. The reason lies in the pressure dependence, as demonstrated via Eq. (42): the chemical potential of electrons ensures that electron captures occurring at larger pressure (deeper in the layer) release more heat. In the instantaneous approach, all the heat is released exactly at the pressure threshold of the first electron capture and the heat release per parent nucleus is given by

$$\begin{array}{c}\hfill q_{\mathrm{th}}=W_1-W_2+E_{\mathrm{exc}}.\end{array}$$(60)

In Fig. 3, we present the heating rate per unit volume $\dot{q}(P)$ for the stationary solution of Z₀ = 26. To emphasize the role of the pressure dependence in the amount of heat released in layers of electron capture, we compare it to the heating rate per unit volume established from a “quasi-instantaneous” approach $\dot{q}(q_{\mathrm{th}})$ in which we assume that the energy release is equal to that at the threshold pressure. The pressure dependence of q implies that the amount of heat released during the active phase is larger in the mixed layer approach than in the instantaneous approach. To illustrate this point, we computed the heat release per accreted baryon (see Eq. (45)) for the stationary solution in the mixed layer approach for four years of active accretion with constant accretion rate. We then compared the outcome to the instantaneous approach, with the results presented in Table 4 for Z₀ = 26,  24,  22, and 20.

Fig. 3. Heating rate per unit volume $\dot{q}$ in the shell with proton number Z0 = 26 as a function of the pressure during active accretion. This quantity is given for the mixed layer approach in blue and compared to the quasi-instantaneous approach in orange.

The stationary profiles presented in Fig. 1 represent the limiting case for the maximum potential heat release: X(P) can never go beyond the stationary profiles at zero temperature. The stationary case corresponds to the maximum increase of the energy release due to the increase of the pressure and chemical potential at the point where most of the reactions take place compared to the instantaneous approach. This increase in δE can be estimated from Eq. (42) as δE = 2δμ_e ≃ 0.5μ_thδP/P_th with μ_th = W₁. This approximation gives the increase in energy release (per nucleon) equal to ∼10,  0.5,  0.2,  0.2keV for Z₀ = 26,  24,  22, and 20, respectively (on the same order as the numbers given in Table 4).

3.2.2. Heat release in quiescence following a stationary profile

We compute the released heat during the forty years of quiescence phase that occur following a stationary profile and present results in Table 4. The heat release per accreted baryon by the deepest reaction layer (Z₀ = 20) is quite large compared to other shells. When considering heat release per reactive baryons however, the values are very similar to what we see for the stationary heat release. For the deepest shell, the heat release per accreted baryon during quiescence is significant compared to the heat release during an active phase in the stationary profile. There are layers of electron captures present in the crust of accreting neutron stars for which the heat release during quiescence phases is quite large; hence, it is not negligible.

Finally, we note that the energy release in the instantaneous and mixed layer approach are on the same order of magnitude because the energy reservoir is directly connected to the mass (energy) difference of a nuclei under consideration. In this sense, the conclusions concerning the shallow heating problem where a much larger energy release is needed are the same as in the instantaneous approach (see e.g., Table 1 in Chamel et al. 2020).

3.2.3. Heat release in the active phase: Time-dependent solution for Z₀ = 26

We computed the released heat per accreted baryon, by solving the full continuity equation during 4 years of active accretion followed by 40 years of quiescence with a time-dependent accretion rate following Eq. (32) and ${\dot{M}}_{\max }={10}^{-8}$M_⊙ in the case of Z₀ = 26. During the active phase, we obtain a heat release of 42.9keV per accreted baryon. We have shown that the pressure dependence of Eq. (42) leads to an increase in heat release ∼22% with respect to the instantaneous approach when considering stationary accretion in Section 3.2.1. However, over the same period of time in active phase, including a realistic time-dependent accretion cycle only leads to an increase of the heat release per accreted baryon by ∼10% when compared to the instantaneous approach. This is due to the decreased number of baryons that effectively react (N_r = 3.9 × 10⁴⁹) compared to the number of accreted baryon (N_b = 4.3 × 10⁴⁹). If we compute the heat release per reactive baryons (and not per accreted baryons) for the time-dependent solution, we obtain 47.3keV per baryon during the active phase. This value is closer to the heat release per accreted baryon in the stationary case.

During quiescence, an additional 4 × 10⁴⁸ baryons are involved in electron captures, leading to a total heat release (active phase and quiescence) of 46.4keV per accreted baryon. This represents an increase of ∼19% compared to the heat release in the instantaneous approach.

3.2.4. Heat deposited

Finally, in Fig. 4, we present the heating rate per unit volume $\dot{q}(P)$ and the neutrino loss rate ε_ν(P) for the stationary solution of Z₀ = 26. The neutrino loss (see Eq. (49)) and total heat deposited (see Eq. (51)) for 4 years of stationary active accretion and for 40 years of quiescence with ${\dot{M}}_{\max }={10}^{-8}$M_⊙, are presented in Table 4. Overall, considerations of neutrino loss in an accreted outer crust with ⁵⁶Fe ashes implies that the heat deposited is approximately one quarter of the heat released. The role of excited states of daughter and grand-daughter nuclei was discussed in Gupta et al. (2007), who reported that the portion of heat release allocated to excitation energy compensated significantly the neutrino heat loss (see discussion in Section 2.8). In the current paper, we only consider excited states for the daughter nucleus of the first electron capture.

Fig. 4. Heating rate per unit volume $\dot{q}$ and neutrino loss εν as a function of the pressure established from the stationary solution of Z0 = 26.

4. Conclusion

We studied the composition and heat release associated to electron captures in a fully accreted outer crust with ⁵⁶Fe ashes when the kinetics of those reactions are taken into account. Making use of the mixed layer approach, we established the relation between the baryon density and the pressure of a lattice of nuclei permeated by a gas of degenerate electrons for a mixture of two nuclear species. Using recent nuclear physics experiments of nuclear mass measurements as well as the theoretical approach of Skyrme force based on the energy density functional described in Goriely et al. (2010), we calculated the threshold density and pressure for each of the four pairs of electron captures in our accreted outer crust. Nuclear data sheets, containing experimental data on nuclei lifetime with respect to β-decays, were used to derive the reaction rates of the four electron captures that dictate the kinetics relevant to our study. We explored the literature to find values for the accretion rate and typical durations for the active and quiescence phases that stem from observations of X-ray transient sources. We derived the continuity equation followed by the parent nuclei abundance and solved it numerically for realistic accretion parameters. The impact of electron capture kinetics on the composition in the shell was evaluated and we computed the thickness of the layers comprising both parent and grand-daughter nuclei for each pair of electron captures in the outer crust. Finally, we studied the impact of kinetics on the heat release, and on the neutrino heat loss in the shells.

We find that the thickness in pressure for the layer containing a mixture of parent and grand-daughter nuclei of electron captures are on the order of 10⁻⁷MeV/fm³ for the four pairs of electron captures in the outer crust. Our calculations were done neglecting thermal effects. Finite temperature would increase the reaction rates for the electrons due to the expansion of the available momentum space. Captures then become possible even below the zero temperature threshold of the reaction (Ushomirsky et al. 2000), shifting the reaction layers to lower densities.

During the active phase of accretion, in the case of a constant accretion rate and studying the stationary solution, we find that including a finite reaction rate leads to an increase in heat release compared to the instantaneous approach that considers an infinite reaction rate for electron captures because electron captures can occur deeper in the shell. For a time-dependent active accretion of four years with a maximum accretion rate of 10⁻⁸M_⊙ per year, we find that the number of accreted baryons N_b is larger than the number of baryons that are involved in electron captures N_r during the active phase of accretion; this is contrary to the stationary solution, for which N_b = N_r. This effect can then significantly decrease the heat release per accreted baryon, as the time-dependent solution X(P, t) does not remain in the stationary profile for the entire phase of active accretion.

We find that heat release occurs also during quiescence phases in the mixed layer approach (contrary to the instantaneous approach) and that said heat for some shells is not negligible, with respect to the heat released during the active phase. The neutrino loss significantly reduces the heat deposited in the matter, however, a more complete consideration of excited states is likely to compensate this loss. Those new results imply that the modeling of accreting neutron star thermal relaxation should be reevaluated by including time-dependent sources that would also cover the quiescent phases of the accretion process.

Data availability

Movie is available at https://www.aanda.org

Acknowledgments

We thank M. Fortin for useful discussions on accreting neutron star sources. We thank the referee for very pertinent remarks which have significantly improved the quality of the paper. The authors acknowledge the financial support of the National Science Center, Poland grants No. 2018/29/B/ST9/02013 and No. 2019/33/B/ST9/00942, as well as the National Science Foundation grant Number PHY 21-16686.

References

Appendix A: Derivation of the β-decay and electron capture rates

From Fermi’s golden rule (Dirac 1927), the probability of transition from an initial state denoted i to a final state denoted f is given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{i\to f}=\frac{2\pi }{ħ}{|{\mathcal{M}}_{i\to f}|}^2\rho (E_f),\end{array}$$(A.1)

with ℳ the matrix element of the reaction and ρ(E_f) the density of final states with energy E_f.

The matrix element for an electron capture ℳ_ec and its corresponding β-decay ℳ_β are established from initial and final state products with the Hamiltonian of weak interaction denoted ${\widehat{H}}_w$ such that

$$\begin{array}{cc}& {\mathcal{M}}_{\mathrm{ec}}={\overrightarrow{p}}_{{\nu }_e},s_{\nu }{\overrightarrow{p}}_{\mathrm{f}},M_{\mathrm{f}}{J}_{\mathrm{f}}{\widehat{H}}_w{\overrightarrow{p}}_{\mathrm{i}},M_{\mathrm{i}}{J}_{\mathrm{i}}{\overrightarrow{p}}_e,s_e,\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& {\mathcal{M}}_{\beta }={\overrightarrow{p}}_{{\overline{\nu }}_e},s_{{\overline{\nu }}_e}{\overrightarrow{p}}_e,s_e{\overrightarrow{p}}_{\mathrm{f}},M_{\mathrm{f}}{J}_{\mathrm{f}}{\widehat{H}}_w{\overrightarrow{p}}_{\mathrm{i}},M_{\mathrm{i}}{J}_{\mathrm{i}},\hfill \end{array}$$(A.3)

with ${\mathit{p}}_{{\nu }_e}$, ${\mathit{p}}_{{\overline{\nu }}_e}$ and p_e the momenta of the neutrino, anti-neutrino and electron, s designating the spin of leptons, and M_fJ_f and M_iJ_i the nuclear spin and its projection on the spin quantization axis of final and initial nuclei.

By averaging over initial states and summing over final states, the reaction rates denoted ℛ are

$$\begin{array}{cc}\hfill {\mathcal{R}}_{\mathrm{ec}}& =\frac{2\pi }{ħ}\langle \sum _{{\overrightarrow{p}}_{{\nu }_e},{\overrightarrow{p}}_{Z-1}}\sum _{J_{Z-1}}{|{\mathcal{M}}_{\mathrm{ec}}|}^2\delta (E_e-W-E_{{\nu }_e}){\rangle }_{{\overrightarrow{p}}_e,s_e,{\overrightarrow{p}}_Z,J_Z},\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill {\mathcal{R}}_{\beta }& =\frac{2\pi }{ħ}\langle \sum _{{\overrightarrow{p}}_{{\overline{\nu }}_e},{\overrightarrow{p}}_e,{\overrightarrow{p}}_Z}\sum _{s_e}\sum _{J_Z}{|{\mathcal{M}}_{\beta }|}^2\delta (W-E_e-E_{{\overline{\nu }}_e}){\rangle }_{{\overrightarrow{p}}_{Z-1},J_{Z-1}},\hfill \end{array}$$(A.5)

where the Dirac delta functions express the conservation of energy, and E_e, E_ν, and $E_{\overline{\nu }}$ are the energy of the electron, of the neutrino and the anti-neutrino; for details on the absence of a sum over the (anti-)neutrino spin, we refer to Goldhaber et al. (1958).

In the heavy nucleus approximation, the transfer of momentum between the initial and final nuclei is neglected so that

$$\begin{array}{cc}\hfill {\mathcal{R}}_{\mathrm{ec}}& =\frac{2\pi }{ħ}\langle \sum _{{\overrightarrow{p}}_{{\nu }_e}}\sum _{J_{Z-1}}{|{\mathcal{M}}_{\mathrm{ec}}|}^2\delta (E_e-W-E_{{\nu }_e}){\rangle }_{{\overrightarrow{p}}_e,s_e,J_Z},\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill {\mathcal{R}}_{\beta }& =\frac{2\pi }{ħ}\langle \sum _{{\overrightarrow{p}}_{{\overline{\nu }}_e}}\sum _{{\overrightarrow{p}}_e}\sum _{s_e}\sum _{J_Z}{|{\mathcal{M}}_{\beta }|}^2\delta (W-E_e-E_{{\overline{\nu }}_e}){\rangle }_{J_{Z-1}}.\hfill \end{array}$$(A.7)

under the assumption that only allowed reactions are studied and that the lepton states are given by plane wave functions, the matrix element of each reaction can be simplified by operating a Taylor expansion of the lepton plane waves and keeping only the zeroth order term (Povh et al. 2004). The matrix element therefore does not depend on the momenta of leptons. Moreover, we make the assumption that the matrix elements relevant for the reactions do not depend neither on the spin of leptons nor on the spin projection M of nuclei. Consequently, we obtain such simplified expressions

$$\begin{array}{cc}\hfill {\mathcal{R}}_{\mathrm{ec}}& =\frac{2\pi }{ħ}{|{\mathcal{M}}_{\mathrm{ec}}|}^2\langle \sum _{{\overrightarrow{p}}_{{\nu }_e}}\delta (E_e-W-E_{{\nu }_e})\sum _{J_{Z-1}}{\rangle }_{{\overrightarrow{p}}_e,s_e,J_Z},\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill {\mathcal{R}}_{\beta }& =\frac{2\pi }{ħ}{|{\mathcal{M}}_{\beta }|}^2\langle \sum _{{\overrightarrow{p}}_{{\overline{\nu }}_e}}\sum _{{\overrightarrow{p}}_e}\delta (W-E_e-E_{{\overline{\nu }}_e})\sum _{s_e}\sum _{J_Z}{\rangle }_{J_{Z-1}}.\hfill \end{array}$$(A.9)

Both the the electrons and (anti-)neutrinos have a continuous spectrum, which translates the sum over their momenta into three dimensional phase-space integrals I_β and I_ec that appear in the reaction rates and are given by

$$\begin{array}{cc}\hfill I_{\mathrm{ec}}& =\frac{2\pi }{ħ}\int \frac{{\mathrm{d}}^3{p}_e}{{(2\pi ħ)}^3}\int \frac{{\mathrm{d}}^3{p}_{\nu }}{{(2\pi ħ)}^3}\delta (E_e-W-E_{\nu }),\hfill \end{array}$$(A.10)

$$\begin{array}{cc}\hfill I_{\beta }& =\frac{2\pi }{ħ}\int \frac{{\mathrm{d}}^3{p}_e}{{(2\pi ħ)}^3}\int \frac{{\mathrm{d}}^3{p}_{\nu }}{{(2\pi ħ)}^3}\delta (W-E_e-E_{\overline{\nu }}).\hfill \end{array}$$(A.11)

These integrals lead to the functions defined in Eq. (23).

Appendix B: Calculating electron capture timescales

In this section of the Appendix, we present details on how the values of the electron capture timescale presented in Table 2 are established for each of the four reactions studied. Unless explicitly stated otherwise, experimental data from National Nuclear Data Center (2022) is used.

B.1. Reaction ⁵⁶Fe → ⁵⁶Mn

In the β-reaction ${}_{25}^{56}$Mn → ${}_{26}^{56}$Fe, the allowed transition is energetically prohibited, so that the inverse of this reaction produces an excited state ⁵⁶Mn^* with a 110keV excitation energy. A doubly forbidden reaction is energetically favorable. However, the timescale of a n-forbidden reaction is much larger than that of an allowed reaction. Very roughly, the quantity log₁₀(ft_1/2) is increased by some 2n + 1 compared to the allowed decay, see Sec. 2.1.4 of Grotz & Klapdor (1990). The timescale for γ-ray deexcitation of ⁵⁶Mn^* is less than 10⁻⁹s, see Sec. 3.4 in Povh et al. (2004). Moreover, we are concerned with the crust of an accreting neutron star which increases continuously the pressure in a piece of matter subject to electron capture. Assuming an accretion rate during the active phase 10⁻⁸M_⊙ per year, the time required for the chemical potential to reach the threshold of the allowed reaction is ∼0.3years, which is small compared to the timescale of the doubly forbidden reaction. Therefore, we should consider an allowed reaction with additional heat (excitation energy) of 110keV from the γ-decay.

Experimental data is available for the quantity ft_1/2 of the β-decay ${}_{25}^{56}$Mn$(3+){\to }_{26}^{56}$Fe(2 + ), and we are interested in establishing the electron capture rate of the reaction ${}_{26}^{56}$Fe$(0+){\to }_{25}^{56}$Mn^*(1 + ). We can now calculate the reaction timescale as

$$\begin{array}{c}\hfill \frac{1}{{\tau }_{\mathrm{ec}}}=\frac{ln(2)}{f{t}_{1/2}}\frac{21}{10}=1.154\times {10}^{-7}{\mathrm{s}}^{-1},\end{array}$$(B.1)

with ft_1/2 = 10^7.101 s.

B.1.1. From ⁵⁶Cr → ⁵⁶V

Experimental data is available for the quantity ft_1/2 of the β-decay ${}_{23}^{56}$V$(1+){\to }_{24}^{56}$Cr(0 + ). We can define the inverse of the timescale of the electron capture as

$$\begin{array}{c}\hfill \frac{1}{{\tau }_{\mathrm{ec}}}=\frac{ln(2)}{f{t}_{1/2}}\frac{9}{2}=7.312\times {10}^{-5}{\mathrm{s}}^{-1},\end{array}$$(B.2)

with ft_1/2 = 10^4.63 s.

B.1.2. From ⁵⁶Ti → ⁵⁶Sc

Experimental data is only available for the quantity T_1/2 of the β-decay, the half-time for the decay of the nucleus through all channels of the reactions. Therefore, we assume that the channel we are studying, that is the β-decay ${}_{21}^{56}$Sc$(1+){\to }_{22}^{56}$Ti(0 + ), is the dominant channel such that T_1/2 ≃ t_1/2, with

$$\begin{array}{c}\hfill T_{1/2}=26\mathrm{ms}.\end{array}$$(B.3)

The function f which appears in the reaction rate of beta decays in Eq. (20) is given by

$$\begin{array}{cc}\hfill f& ={\int }_1^{\overline{W}}\sqrt{{\overline{E}}^2-1}\overline{E}{(\overline{E}-\overline{W})}^2\mathrm{d}\overline{E}.\hfill \end{array}$$(B.4)

However, among the assumptions used on the matrix element of the β-decay, we have neglected the impact of the Coulomb interaction between electrons and proton on the shape of the wave functions of particles under consideration. To correct for that, it is customary to introduce the Fermi function (see, e.g., Sect. 2.1.3 in Grotz & Klapdor (1990)). The non-relativistic formulation of the Fermi function is easier to implement than the relativistic one, however we have considered relativistic particles so far. To use the non-relativistic formulation, we have verified that the quantity log(ft_1/2) calculated with the non-relativistic formulation of the Fermi function is close the the experimental data value for the electron captures involved in the first two shells of the accreted crust; we found a relative difference of at most 8%, and therefore kept the non-relativistic formulation of the Fermi function in the evaluation of ft_1/2. For the β-decay ${}_{21}^{56}$Sc$(1+){\to }_{22}^{56}$Ti(0 + ) reaction rate we obtain

$$\begin{array}{c}\hfill \frac{1}{{\tau }_{\mathrm{ec}}}=\frac{ln(2)}{f{t}_{1/2}}\frac{9}{2}=1.362\times {10}^{-4}{\mathrm{s}}^{-1},\end{array}$$(B.5)

with ft_1/2 = 10^4.36 s.

B.1.3. From ⁵⁶Ca → ⁵⁶K

Experimental data for this reaction is not available but calculations in the quasi-particle Random Phase Approximation are presented in Minato et al. (2021) for the quantity T_1/2 of the β-decay ${}_{19}^{56}$K ${\to }_{20}^{56}$Ca. The energy level diagrams are not available either in National Nuclear Data Center (2022), therefore we assume a (0 + ) → (0 + ) reaction. Similarly to the treatment of the previous reaction, we obtain the rate for the electron capture

$$\begin{array}{c}\hfill \frac{1}{{\tau }_{\mathrm{ec}}}=\frac{ln(2)}{f{t}_{1/2}}\frac{1}{2}=6.307\times {10}^{-6}{\mathrm{s}}^{-1}.\end{array}$$(B.6)