Accreting neutron stars from the nuclear energy-density functional theory. II. Equation of state and global properties

The accretion of matter onto the surface of a neutron star in a low-mass X-ray binary triggers X-ray bursts, whose ashes are buried and further processed thus altering the composition and the properties of the stellar crust. In this second paper of a series, the impact of accretion on the equation of state and on the global properties of neutron stars is studied in the framework of the nuclear energy-density functional theory. Considering ashes made of $^{56}$Fe, we calculated the equations of state using the same Brussels-Montreal nuclear energy-density functionals BSk19, BSk20, and BSk21, as those already employed for determining the crustal heating in our previous study for the same ashes. All regions of accreting neutron stars were treated in a unified and thermodynamically consistent way. With these equations of state, we determined the mass, radius, moment of inertia, and tidal deformability of accreted neutron stars and compared with catalyzed neutron stars for which unified equations of state based on the same functionals are available. The equation of state of accreted neutron stars is found to be significantly stiffer than that of catalyzed matter, with an adiabatic index $\Gamma \approx 4/3$ throughout the crust. For this reason, accreting neutron stars have larger radii. However, their crustal moment of inertia and their tidal deformability are hardly changed provided density discontinuities at the interface between adjacent crustal layers are properly taken into account. The enhancement of the stiffness of the equation of state of accreting neutron stars is mainly a consequence of nuclear shell effects, thus confirming the importance of a quantum treatment as stressed in our first study. With our previous calculations of crustal heating using the same functionals, we have thus obtained consistent microscopic inputs for simulations of accreting neutron stars.


Introduction
The accretion of matter onto a neutron star (NS) in a lowmass X-ray binary induces thermonuclear explosions with total energies 10 39 − 10 40 erg observed as X-ray bursts lasting a few tens of seconds and with a typical recurrence time of hours to days (see e.g. Strohmayer & Bildsten (2006) for a review and Galloway et al. (2020) for a recent compilation of observations). The unstable carbon burning in deeper layers is thought to power less frequent superbursts with energies ∼ 10 42 erg. On a longer timescale, the ashes from Xray bursts may be further reprocessed by electron captures, neutron emissions, and neutron captures, and possibly pycnonuclear fusion reactions as they slowly sink inside the NS (see, e.g., Meisel et al. (2018) for a recent review). Depending on the duration of accretion episodes, the original crust of the NS can be partially or fully replaced (Tauris et al. (2012); see also the discussion in Fantina et al. (2018)) and this may radically change its properties, as first shown by Haensel & Zdunik (1990a).
Simulations of accreting NSs usually rely on the equation of state (EoS) of Haensel & Zdunik (1990b) for the accreting crust together with the heat sources obtained by Haensel & Zdunik (1990a) (or the more recent calculations of Haensel & Zdunik (2003, 2008). These results, which are based on the liquid-drop model, are then combined with some EoS for the deepest layers of the inner crust and for the core. However, significant errors on the global structure of the NS and on its tidal deformability may arise from the use of thermodynamically inconsistent EoSs (Fortin et al. 2016;Ferreira & Providência 2020;Suleiman et al. 2021).
In our previous study, we calculated the heat released in accreting NS crusts and the location of the heat sources using accurately calibrated Brussels-Montreal nuclear energydensity functionals (EDFs) taking into account nuclear shell effects perturbatively . We found that the composition of the crust of an accreted NS can differ substantially from that of an isolated NS made of cold catalyzed matter. In this second paper, we present the corre-sponding EoSs. Both the outer and inner regions of accreted NS crusts are treated consistently within the same microscopic framework using the same EDFs as in Fantina et al. (2018). Together with the EoSs of the liquid core previously calculated in Goriely et al. (2010), unified EoSs describing all regions of accreted NSs with the same EDFs are constructed.
After briefly recalling the main features of the microscopic model of accreted NS crusts and of the adopted EDFs in Sect. 2, the calculated EoSs are presented in Sect. 3. The global properties of accreted NSs are discussed in Sect. 4 and conclusions are given in Sect. 5.

Main assumptions
We recall here the main assumptions on how we model the crust of an accreting NS, by which we mean the region of the NS below the envelope and whose density is higher than 10 6 g cm −3 . Because of the relatively low temperatures generally prevailing in accreting NS crusts, namely T < 10 9 K (see e.g. Chamel & Haensel (2008)), the thermal contributions to the thermodynamic potentials can be neglected. For such temperatures, except for the shallowest regions, the innermost layers are expected to be in the solid state Carreau et al. 2020). The composition and structure of each layer were determined in the one-component plasma approximation: at each pressure P , matter is fully described by the proton number Z and the total nucleon number A cell in a single Wigner-Seitz (WS) cell. In the outer crust, A cell reduces to the mass number A of the corresponding nuclei, while in the inner crust A cell also accounts for free neutrons. Unlike cold catalyzed matter, accreting NS crusts are not in full thermodynamic equilibrium and reactions can thus occur. The compression of the ashes from X-ray bursts may trigger a series of electron captures, by which the nucleus (A, Z) transforms into a nucleus (A, Z − 1) with the emission of an electron neutrino. The one-component plasma approach induces discontinuities in the neutron density and the neutron chemical potential, while in fact the onset of neutron emission occurs at lower density and pressure than those predicted by the mean-nucleus approximation (see  for details). However, we did not consider neutron diffusion (Shchechilin & Chugunov 2019;Chugunov & Shchechilin 2020;Gusakov & Chugunov 2020), for consistency with our previous calculations of the heating . With increasing density and decreasing proton number, neighbouring nuclei may overcome their Coulomb barrier and fuse. However, large uncertainties reside in the determination of the rates of these pycnonuclear reactions (see e.g. Yakovlev et al. 2006). We thus simply assumed that these processes occur whenever Z reaches the minimum value Z min = 8 (see also the discussion in Sect. 3.2 in Fantina et al. (2018)). As in Fantina et al. (2018), we supposed that the crust is fully accreted.

Brussels-Montreal nuclear energy-density functionals
Our models of accreted NS crusts are based on the nuclear EDF theory Stone & Reinhard 2007), which has proved successful in describing finite nuclei, such as those encountered in the outer crust of a NS , but also provides a consistent treatment of the inner crust, where neutron-proton clusters coexist with a neutron liquid. In this quantum mechanical approach, bound and unbound neutrons are consistently calculated. They do not need to be considered separately as in the liquid-drop treatment of Haensel & Zdunik (1990a,b, 2003, 2008. Moreover, the EDF approach is suitable to describe homogeneous matter, thus allowing for a unified description of all regions of a NS. Here, we employ the same Brussels-Montreal EDFs labelled BSk19, BSk20, and BSk21 ) as in Fantina et al. (2018).
The Brussels-Montreal EDFs BSk19, BSk20, and BSk21 are based on generalised Skyrme effective nucleon-nucleon interactions (Chamel et al. 2009;Goriely et al. 2010), supplemented with a microscopic contact pairing interaction (Chamel 2010). These EDFs were fitted to the 2149 measured masses of nuclei with neutron and proton numbers, N, Z ≥ 8, given in the 2003 Atomic Mass Evaluation (AME) (Audi et al. 2003), with a root-mean square (rms) deviation as low as 0.58 MeV for the three functionals, while ensuring at the same time an optimal fit to charge radii. Incidentally, these models also yield an equally good fit to the experimental masses of nuclei with N, Z ≥ 8 given in the latest AME Wang et al. 2021). The masses of bound nuclei were obtained by adding to the Hartree-Fock-Bogoliubov (HFB) energy a phenomenological Wigner term and correction term for the collective energy.
The BSk19, BSk20, and BSk21 EDFs were simultaneously constrained to reproduce various properties of homogeneous nuclear matter as obtained from many-body calculations using realistic two-and three-nucleon interactions. Specifically, these EDFs were fitted to three different neutron-matter EoSs, reflecting the current uncertainties of the high-density behaviour of dense matter. BSk19 was adjusted to the 'soft' EoS of neutron matter of Friedman & Pandharipande (1981) obtained from the realistic Urbana v14 nucleon-nucleon force with the three-body force TNI, BSk20 was fitted to the EoS of Akmal et al. (1998) labelled "A18 + δv + UIX", whereas BSk21 was constrained to reproduce the 'stiff' EoS labelled "V18" from Li & Schulze (2008). All three neutron-matter EoSs are consistent at the low densities prevailing in nuclei and in NS crusts and are compatible with the more recent calculations based on auxiliary field diffusion Monte Carlo method and chiral effective field theory (see Fantina et al. (2014) and Fig. 13 in Fantina et al. (2018)).
Furthermore: (i) the symmetry energy coefficient was set to J = 30 MeV for the three functionals (see e.g. Goriely et al. 2010;Tsang et al. 2012;Baldo & Burgio 2016;Tews et al. 2017;Burgio & Fantina 2018, for a discussion of experimental and theoretical estimates of J), (ii) the incompressibility K v of symmetric nuclear matter at saturation was required to fall in the experimental range 240±10 MeV (Colò et al. 2004;Stone et al. 2014), (iii) the ratio of the isoscalar effective mass to bare nucleon mass in symmetric nuclear matter at saturation was set to the realistic value of 0.8 (Goriely et al. 2003), (iv) the isovector effective mass was found to be smaller than the isoscalar effective mass, in agreement with both experiments and many-body calculations (see Goriely et al. 2010, for a discussion), (v) a qualitatively realistic distribution of the potential energy among the four spin-isospin channels in nuclear matter was obtained, and (vi) spurious spin and spin-isospin instabilities in nuclear matter that generally plagued earlier Skyrme functionals were eliminated for all densities prevailing in NSs (Chamel et al. 2009;. Finally, all these functionals are consistent with constraints on the EoS of symmetric nuclear matter inferred from heavy-ion collisions (Danielewicz et al. 2002;Lynch et al. 2009).
The BSk19-21 EDFs have been already applied to calculate the structure, the composition, and the EoS of nonaccreted NSs under the cold catalyzed matter hypothesis Pearson et al. 2012) and analytical representations of these EoSs have also been obtained by Potekhin et al. (2013). Even though the EoS based on the BSk19 EDF fails to explain the existence of the most massive observed NSs ) (see also Fantina et al. (2013) for a discussion of various astrophysical constraints), such a soft EoS seems to be favoured by the analyses of K + production (Fuchs et al. 2001;Sturm et al. 2001;Hartnack et al. 2006) and π − /π + production ratio (Xiao et al. 2009) in heavy-ion collisions. This EoS also appears to be supported by the LIGO-Virgo data from the gravitational-wave signal GW170817 emitted by the merger of intermediate mass NSs (Perot et al. 2019). We have thus employed these three EDFs to calculate the EoS of all regions of accreted NSs, as described in the following, and we have studied their impact on NS properties.

Outer crust
We assumed that the outer crust is made of fully ionised atoms embedded in a charge-neutralising degenerate electron background. Since in hydrostatic equilibrium the pressure P varies continuously throughout the star, the suitable thermodynamic potential for determining the equilibrium structure of the crust is the Gibbs free energy per nucleon g(A, Z, P ). The composition and EoS in each layer of pressure P were thus found by minimising g, given by E = ρc 2 being the internal energy density and n the average baryon number density, keeping the mass number A fixed (recalling that in the outer crust the number of nucleons in the cell is equal to the nucleus mass number, i.e. A cell = A).
For the specific expressions of E and P , we followed Pearson et al. (2011). To recall the main points, the internal energy per nucleon E/n in the WS cell is given by where M ′ (A, Z) is the nuclear mass, c is the speed of light, and E L is the lattice energy density, with e the elementary charge, and we adopted the WS estimate C = −1.4508 for the lattice constant (Salpeter 1961). The last term in Eq.
(2), E e , is the energy density of electrons, which accounts for the energy density of a free uniform electron gas of number density n e (see Weiss et al. (2004) for complete expressions), including the exchange and polarisation contributions, for which we used Eqs.
(1), (6), and (8) in . We therefore neglected the electron correlation correction, as well as the finite-size and zero-point contributions, which were shown to be negligible ). The electron number density is related to the baryon number density by n e = (Z/A)n to ensure electric charge neutrality.
Similarly, for the pressure we have with P e the pressure of the electron (including exchange and charge-polarisation corrections) and P L = E L /3 the lattice pressure, the nuclei exerting no pressure in cold matter. The only microscopic inputs of this model are the nuclear masses M ′ (A, Z), that were obtained from the atomic masses M (A, Z) after subtracting out the binding energy of atomic electrons (see Eq. (A4) in Lunney et al. (2003)). For consistency with our previous study , we made use of the experimental atomic masses from the 2016 AME (Wang et al. 2017), whenever available 1 , complemented with the Brussels-Montreal HFB atomic masses from the corresponding functional, also available from the BRUSLIB database (Xu et al. 2013).
With increasing pressure, Z remains unchanged until P reaches the threshold pressure P β for the onset of electron captures and such that g(A, Z, P β ) = g(A, Z − 1, P β ). Accurate analytical formulae for P β and the corresponding density ρ β were obtained by . The daughter nucleus (A, Z −1) is generally unstable and further electron captures occur until the proton number reaches some value Z 0 < Z such that g(A, Z 0 , P β ) < g(A, Z 0 − 1, P β ), which corresponds to a local minimum of the Gibbs free energy per nucleon at pressure P β .
The onset of neutron drip was determined as in , thus ensuring the thermodynamic consistency across the boundary between the outer and inner crusts: the proton number of the dripping nucleus is given by the highest value of Z (below that of the initial ashes) for which the ∆N neutron separation energy is negative, i.e.
m n being the neutron mass. The properties of the accreted crust at the neutron-drip transition are summarised in Table 1, where we list, for the considered models, the proton number Z of the dripping nucleus, the number ∆N of emitted neutrons, the mass-energy density ρ drip , and the corresponding pressure P drip .
For the inner crust, we calculated the EoS using the 4th order Extended Thomas-Fermi (ETF) approach, with proton shell corrections added perturbatively via the Strutinsky integral (SI) theorem (see Onsi et al. (2008) and Pearson et al. (2012) for details). Neutron shell corrections, which are expected to be much smaller than proton shell corrections (Pearson & Chamel 2022), were neglected. This so-called ETFSI approach is a computationally high-speed approximation of the fully self-consistent HFB method. We employed the same EDFs as those underlying the HFB nuclear mass models used for the outer crust, thus ensuring a consistent description of both the outer and the inner regions of the crust. We further assumed that nuclear clusters remain spherical and we calculated the Coulomb energy in the WS approximation. Indeed, the presence of so-called nuclear 'pasta' has recently been found to be marginal within the ETFSI approach (Pearson & Chamel 2022). In order to further reduce the computation time, the nucleon density distributions in the WS cell were parametrised as n q (r) = n Bq + n Λq f q (r), where q = n, p for neutrons and protons respectively, r is the radial distance from the centre of the spherical cell, n Bq is the background density, and the second term accounts for the presence of inhomogeneities, with n Λq a constant term and f q (r) a damped Fermi function ensuring that density derivatives vanish at the surface of the cell, thus providing a smooth matching of the nucleon distributions between adjacent cells (Onsi et al. 2008).
Before the onset of pycnonuclear reactions, A cell remains unchanged with increasing pressure whereas Z can decrease if electron captures occur. The EoS can in principle be determined by minimising the baryon chemical potential, which is equal to the Gibbs free energy per nucleon g, with respect to all the parameters of the WS cell keeping A cell fixed. We note that, because of the presence of the free neutron gas, A cell > A with the number A of nucleons in clusters defined by with R cell the WS cell radius. However, this procedure would be computationally very costly. Instead, we proceeded as follows: (1) We first minimised the internal energy E/n per nucleon at constant average baryon number density n, for given values of the proton number Z and nucleon number A cell (see Pearson et al. 2012, for details); (2) For this set of parameters (n, Z, A cell ), we calculated the corresponding pressure P and the Gibbs free energy per nucleon g; (3) We repeated the calculations for different values of n and Z, thus allowing us to construct g as a function of Z and P for the given (fixed) value of A cell .

EoS for accreted neutron-star crusts
We calculated the composition and the EoSs of fully accreted NS crusts for the EDFs BSk19, BSk20, and BSk21, considering X-ray burst ashes made of 56 Fe as in our previous study ).

Pressure-density relation
Results are plotted in Fig. 1. The EoS tables containing the number density n, the mass-energy density ρ, and the pressure P for the three EDFs are available at the CDS. As we can see from the left panel of Fig. 1, the three EoSs exhibit a similar behaviour, with steps associated to changes in the composition induced by electron captures and/or pycnonuclear reactions ). On the other hand, remarkable differences are observed between the EoS of accreted and catalyzed crusts, the former being notably stiffer than the latter, particularly in the density region between the neutron drip and ≈ 10 13 g cm −3 . This is clearly seen from the right panel of Fig. 1, where the EoS based on the BSk19 EDF is shown for both catalyzed (blue dashed line) and accreted (red solid line) crusts, as an illustrative example (the other EDFs yielding qualitatively similar results). One can also notice that, at higher densities, ρ > ∼ 10 13 g cm −3 (corresponding to P > ∼ 5×10 31 dyn cm −2 ), the two EoSs merge. The global properties of accreted NSs can thus be fairly accurately calculated by combining the EoS of accreted crust with that of catalyzed matter for densities beyond ρ > ∼ 10 13 g cm −3 (deep inner crust and core, see Sect. 4). These calculations assumed that the crust is fully accreted and stationary. This assumption is valid if the star has accreted for a long enough time, approximately given by (Suleiman et al. 2022) where P 31 is the pressure in 10 31 dyn cm −2 at the crustcore transition,Ṁ −10 is the accretion rate in 10 −10 M ⊙ /yr (M ⊙ being the mass of our Sun), and R 6 = R/(10 km). Setting P 31 = 60 (Pearson et al. 2012), M = 1.4 M ⊙ , and R = 12 km, we find τ acc ∼ 430 Myr/Ṁ −10 . On longer timescales, the EoS remains unchanged with further accretion and is therefore well defined. On the other hand, on a shorter accretion timescale, the NS crust is only partially replaced by the accreted material and a different treatment is required (Suleiman et al. 2022).
To assess the importance of shell effects on the EoS, we also display the results obtained using the ETF approach (black solid line in the right panel of Fig. 1). As previously discussed in Fantina et al. (2018), the compressed matter element undergoes in this case successive electron captures until the proton number is low enough to allow for pycnonuclear reactions to occur. As a consequence, variations of Z are more continuous and jumps in density are much smaller with respect to those obtained with the full ETFSI treatment. Even though the composition was found to be radically different, we showed that the proton fractions Z/A cell are comparable (see Fig. 4 in Fantina et al. (2018)). Since the pressure is mainly determined by electrons and ρ ≈ n e m u A cell /Z, with m u the unified atomic mass unit (recalling that n e = nZ/A cell due to electric charge neutrality), the EoS of accreted NS crusts is thus expected to be very similar to that of catalyzed crusts when shell corrections are neglected, as can be seen in Fig. 1. This shows that nuclear shell effects are important not only for determining the composition but also for the EoS of accreted crusts.

Adiabatic index
To better see the differences between the EoS of accreted and catalyzed crusts, we calculated the adiabatic index, Results are plotted in Fig. 2 for models based on the BSk19 EDF. We focus on pressures around P ≈ 10 30 − 10 31 dyn cm −2 (or, equivalently, ρ ≈ 10 12 − 10 13 g cm −3 ) corresponding to the largest deviations in Fig. 1.
Assuming that the pressure in the shallow layers of the inner crust remains mainly determined by electrons, we have P ≈ P (n e ) = P (nZ/A cell ). The adiabatic index can thus be expressed as where in the last equality we have considered the pressure of an ultrarelativistic electron Fermi gas being the Dirac-Planck constant. In catalyzed crusts, Z/A cell decreases continuously with increasing n so that d(Z/A cell )/d log(n) < 0, see e.g. Fig. 4 in Fantina et al. (2018). After the onset of neutron drip, Z/A cell falls rapidly leading to a sharp drop of Γ, as shown by the blue dashed line in Fig. 2. In accreted crusts, both Z and A cell remain unchanged until electrons are captured by clusters or pycnonuclear fusion occurs. Therefore, the second term in Eq. (9) vanishes beyond the neutron-drip point and Γ remains equal to 4/3 as in the densest layers of the outer crust. For the same reason, the adiabatic indices obtained within the ETF (black solid line) and ETFSI (red solid line) approaches are comparable. With further increase of the density, the neutron contribution to the pressure can no longer be ignored and Γ > ∼ 4/3. This shows that the polytropic approximation according to which Γ is replaced by a constant is much more accurate for accreted crusts than for catalyzed crusts.

Global structure of neutron stars
The structure of a non-rotating NS is determined by the well-known Tolman-Oppenheimer-Volkoff (TOV) equations (Tolman 1939;Oppenheimer & Volkoff 1939) dP (r) dr = − Gρ(r)M(r) r 2 1 + P (r) c 2 ρ(r) and The numerical solution of the TOV Eqs. (11) and (12) is found by integrating from the centre of the star at r = 0 up to the surface at (circumferential) radius r = R such that P (R) = 0. The gravitational mass of the star is then given by M ≡ M(R). In practice, the TOV equations were integrated down to the density ρ = 10 6 g cm −3 .
Numerical results are shown in Fig. 3 for both accreted and catalyzed NSs for the EoSs based on the EDFs BSk19, Article number, page 5 of 10  BSk20, and BSk21. We checked that the approximation proposed by Zdunik et al. (2017) for determining the massradius relation from the core EoS only remains equally accurate for the EoSs calculated here. As can be seen in Fig. 3, accreted NSs (solid lines) have larger radii than their nonaccreted relatives made of catalyzed matter (dashed lines). This stems from the fact that the EoS of accreted NSs is stiffer, as discussed in Sect. 3. The differences in the radii is approximately given by (Zdunik et al. 2017) where R is the radius of the catalyzed NS and Q tot is the total energy per nucleon released in accreted crust. For comparison, we show in Fig. 3 the recent mass and radius estimates from NICER observations of PSR J0030+0451 (Miller et al. 2019;Riley et al. 2019) and PSR J0740+6620 (Miller et al. 2021;Riley et al. 2021). Contours for the source PSR J0740+6620 include XMM-Newton priors on the mass. The contours taken from Miller et al. (2019) include both the two and three oval-spot analyses of the signal. PSR J0030+0451 and PSR J0740+6620 are two millisecond pulsars and are therefore expected to have been recycled by accretion from a companion in their past evolution. The EoS of accreted NSs should thus be more realistic than that of catalyzed matter. However, the deviations δR remain much smaller than the current observational uncertainties on the radii. Still, these observations clearly favour the EoSs based on the BSk20 and BSk21 EDFs. Incidentally, the latter is the one that was previously found to be in best agreement with nuclear data ) and various other astrophysi- Fig. 3. Mass-radius relation for accreted and catalyzed NSs for the EoSs based on the BSk19, BSk20, and BSk21 EDFs. Results obtained from Eq. (13) are indistinguishable from the exact ones. The coloured areas represent the recent mass-radius measurements (at 1σ) from NICER of PSR J0030+0451 and PSR J0740+6620. The analyses of Riley et al. (2019Riley et al. ( , 2021 are shown in green and those of Miller et al. (2019Miller et al. ( , 2021 in blue. See text for details. cal observations . Even though the EoS based on the BSk19 EDF fails to reproduce the mass and radius of the very massive pulsar PSR J0740+6620, it remains compatible with those of the less massive pulsar PSR J0030+0451. This comparison is consistent with the analysis of GW170817.

Moment of inertia
The moments of inertia I of accreted and catalyzed NSs were calculated in the slow rigid rotation approximation (Hartle 1967), as implemented in Haensel & Prószyński (1982)  whereω(r) is the local spin frequency as measured in a local inertial frame, j(r) represents the contribution from the sphere of radius r to the stellar angular momentum J = j(R), Ω is the (uniform) angular frequency of the star expressible as and we have introduced the metric function Φ(r) defined by Results for the BSk21 EDF are shown in the top left panel of Fig. 4 for different NS masses. The differences ∆I in the moment of inertia between accreted and catalyzed NSs are indistinguishable and are of the order of ∆I ∼ 10 −7 I.
Of particular interest for the analysis of pulsar frequency glitches (Andersson et al. 2012;Chamel 2013) is the fractional moment of inertia of the crust I crust /I, with r cc denoting the radial coordinate at the crust-core boundary determined by the pressure P cc = 4.294×10 32 dyn cm −2 (Pearson et al. 2012). Results are shown in the top right panel of Fig. 4. Whether the crust is accreted or not hardly changes the crustal moment of inertia. The deviations ∆I crust plotted in the bottom panel of Fig. 4 amount to 0.1% of I crust at most.

Tidal deformability
To assess the importance of accretion on the tidal deformability of a NS, we calculated the tidal deformability function: with the tidal Love number function k 2 (r) given by (Hinderer 2008;Hinderer et al. 2010) where we have introduced the dimensionless compactness parameter and the function y(r) is obtained by integrating the following differential equation with the boundary condition y(0) = 2 (Postnikov et al. 2010): Q(r) = 4πGr 2 /c 4 1 − 2GM(r)/(rc 2 ) 5E(r) + 9P (r) c s = c dP/dE denoting the sound speed. Results for the functions y(r) and λ(r) for a NS with a gravitational mass M = 1.4 M ⊙ are plotted in Figs. 5 and 6 respectively, focussing on the outer region where the EoS of accreted NSs differs from that of catalyzed NSs. Calculations were made for the EoS based on the BSk21 EDF. Accreted NSs are found to have a slightly smaller tidal deformability than catalyzed NSs. However, the deviations almost completely disappear when density discontinuities at the interfaces between adjacent crustal layers are properly taken into account following the prescription of Damour & Nagar (2009), as can be seen by comparing the solid lines and the dashed lines. The (observable) tidal deformability of the star given by Λ = λ(R) is displayed in the left panel of Fig. 7  They are of the same order 10 −5 − 10 −6 as the precision of the calculations with the prescription of Damour & Nagar (2009). The Love numbers and tidal deformability coefficients previously calculated for catalyzed matter with the EoSs based on the BSk19, BSk20, and BSk21 EDFs as well as the comparison with the LIGO/Virgo measurement of GW170817 (Perot et al. 2019) therefore remain valid for accreted NSs. Ignoring density discontinuities leads to errors of the order of a few % (below 2% for NSs with masses M > ∼ M ⊙ ); although these errors remain rather small, they are larger than the relative contribution of the crust to Λ (Perot et al. 2020). A proper treatment of discontinuities is therefore necessary to correctly assess the role of the crust.

Conclusions
In this paper, we presented the EoSs for accreted NSs, calculated employing the Brussels-Montreal EDFs BSk19, BSk20, and BSk21, and we evaluated their impact on global NS properties. Both the outer and inner crusts and the core were described with the same functional, thus ensuring a unified and thermodynamically consistent treatment. The EoS of accreted NS crusts is found to be remarkably different from that of catalyzed crusts, the former being notably stiffer than the latter, particularly in the density region between the neutron drip and ≈ 10 13 g cm −3 . This reflects in the adiabatic index, whose value ≈ 4/3 remains rather constant throughout the accreted crust, implying that the polytropic approximation is much more accurate for accreted crusts than for catalyzed ones. The different stiffness between accreted and catalyzed EoSs also impacts the global NS structure. Indeed, accreted NSs with low and intermediate masses have larger radii than the non-accreted counterparts. However, the deviations remain smaller than current observational uncertainties on the radii. Likewise, the crustal moment of inertia and the tidal deformability of accreted NSs are almost indistinguishable from those of non-accreted NSs, provided density discontinuities at the interface between adjacent crustal layers are properly taken into account.
Although we considered here three specific models, to a good approximation the EoSs calculated for catalyzed matter remain applicable to compute the global properties of