Blaise Tshibangu-Mbuebue

Rober Rojas Laguna

Min Won Lee
email: min.lee@univ-paris13.fr

Jorge Rodríguez-Asomoza

Ignacio Enrique Zaldívar-Huerta

Rojas Laguna
email: rlaguna@ugto.mx

R ; Lee

Numerical Study of a Reconfigurable Multiband Microwave Photonic Filter Using a Tunable Fabry-Perot Filter

Keywords: microwave photonic filter, microwave signals, Fabry-Perot filter, VPIphotonics software

Microwave photonic filters (MPFs) with the capabilities of bandwidth reconfigurability and frequency tunability in the GHz range are of great interest in high-speed communications systems. In this paper, we propose a new reconfigurable multiband microwave photonic filter (MPF) using a tunable Fabry-Perot Filter (FPF). It is demonstrated by numerical simulations that the modification of the intermodal separation (δλ) of a multimode laser diode (MLD) by tunable FPF allows for the reconfigurable multi-passband of the MPF. Our simulation results show that our new filter system is promising to communications systems.

Introduction

At the present time, the delivery of services, such as the Internet, High-Definition Television (HDTV), data, etc., to the final user has begun to migrate to systems based on the use of optical fibers. In this technology, electrical signals are transmitted using photonic devices and optical fiber. This multidisciplinary field is possible thanks to high-speed photonic devices operating at microwave or millimeter-wave frequencies in microwave photonic systems [1]. In this regard, the MPFs are to be applied within microwave photonic technology. Potential applications of an MPF include the distribution and transmission of microwave signals via broadband optical links. Multiband MPFs have become excellent candidates with their inherent advantages from photonic technologies, such as their low loss, high bandwidth, high tunability, and high reconfigurability [START_REF] Capmany | A tutorial on microwave photonic filters[END_REF][START_REF] Yao | Photonics to the Rescue: A Fresh Look at Microwave Photonic Filters[END_REF]. They are also attractive in optical communication systems because they can be used as electrical carriers to transmit information [START_REF] Correa-Mena | Experimental Transmission of Digital Data Coded on Electrical Carriers at 2.1 GHz and 4.2 GHz by Using a Microwave Photonic Filter[END_REF]. Currently, the tunability and reconfigurability of the frequency response of an MPF is a great topic of interest for researchers. For instance, a novel MPF was proposed based on a tunable ultra-sharp roll-off optical filter to demonstrate a flexible MPF with a large reconfigurable bandwidth and a wide frequency tuning range [START_REF] Ding | Wide Bandwidth Reconfigurable and Frequency Tunable Microwave Photonic Filter Based on Tunable Ultra-sharp Roll-off Optical Filter[END_REF]. In [START_REF] Liu | Optical Spectral Slicing Based Reconfigurable and Tunable Microwave Photonic Filter[END_REF], a reconfigurable and tunable multiband MPF was also developed, and flexible amplitudes and phase controls of the filter were achieved via optical spectral slicing by means of a programmable sampling function. A reconfigurable and application-specific tunable MPF with finite impulse response was also reported [START_REF] Mirza | Design and optimization of a microwave photonic filter exploiting differential mode group delay of a multi-mode fiber[END_REF]. This filter exploits the differential mode group delays of linearly polarized modes in a multi-mode fiber to achieve single and multiple high-frequency passband responses. In our previous report, we proposed a tunable filter with a section of birefringent optical fiber placed between two crossed polarizers in order to modify the optical spectrum of an MLD. Thanks to the fiber, the intermodal separation (δλ) is varied. Therefore, it allows for reconfigurability in the frequency of an MPF. However, the drawback of this filter is the birefringent fiber itself. When a particular passband is needed, it is necessary to use a particular length of fiber. In this new work, we propose a reconfigurable MPF using a tunable FPF in a dynamic way. The free spectral range (FSR) of the FPF is modified by an applied voltage, thus adjusting the intermodal separation δλ of the MLD. In this way, reconfigurability in the frequency of the MPF is achieved. In this work, we study the feasibility of the reconfigurable MPF in the frequency range of 0-12 GHz. Due to the COVID-19 pandemic, we study this only by means of simulation using VPIphotonics software (Version 11.1, VPIphotonics GmbH, Berlin, Germany) [8]. The rest of this paper is organized as follows. The operation principles of an FPF and the MPF used in this work are described in Section 2. The simulation results of reconfigurability of the passbands by using a tunable FPF are shown in Section 3. Finally, conclusions are given in Section 4.

Brief Principle of Operation of an FPF and the MPF

In order to familiarize the reader with the basic operation of an FPF and an MPF, this section is devoted to explaining the basic principle of an FPF first, followed by that of the MPF.

Principle of Operation of an FPF

An FPF consists of a cavity formed by an index material refraction n embedded between two parallel flat and highly reflective mirrors, separated by a distance d. This filter possesses one input and one output port, the mirrors of which constitute the resonating cavity creating the multiple-beam interference process. Modification of the spacing between the mirrors can be done manually and/or electro-mechanically. The fastest tuning is achieved by means of a piezo actuator. The modification of spacing generates a tunable filter. The function transfer of the FPF in terms of optical power is given as [START_REF] Kasap | Optoelectronics and Photonics: Principles and Practices[END_REF].

T( f ) = (1 -A -R) 2 (1 -R) 2 + 4R s sin 2 π f FSR ( 1
)
where A represents the signal losses within each mirror, R is the optical power reflectivity, FSR = 1/2τ, with τ = nd/c, which is the time it takes for the signal to propagate from one cavity mirror to the other, and c is the speed of light in a vacuum. Figure 1 shows the transfer function of the FPF that is composed of a sequence of maximum resonances that satisfy the transmission condition [START_REF] Venghaus | Wavelength Filters in Fibre Optics[END_REF]. The distance between two successive peaks is known as the free spectral range (FSR), also called axial mode spacing. In this work, varying the FSR of the simulated FPF allows for the selection of specific optical modes from the multi-modes optical source. As a consequence, the adjustment of δλ parameter in a dynamic way is achieved. When the SM-SF is dispersive to the wavelength used, the frequency response of the system is composed of a series of filtered signals or microwave passbands windows [START_REF] González-Mondragón | Multiple passband microwave photonic filter with adjustable bandwidth[END_REF]. These passband windows are generated by considering the δλ parameter of the MLD, the length (L) of the optical fiber, and its associated chromatic dispersion parameter (D). The corresponding value to the central frequency of the nth filtered microwave band-pass windows can be determined by the equation [START_REF] González-Mondragón | Multiple passband microwave photonic filter with adjustable bandwidth[END_REF].

Principle of Operation of the MPF

f n = n DLδλ ( 2
)
where n = 1, 2, 3 . . . From Equation (2), it is evident that modifying the δλ parameter allows for reconfigurability in the frequency of the filtered microwave band-pass windows. To illustrate this parameter, Figure 3 shows a typical optical spectrum of a commercial MLD (NX7538BF-AA, Renesas Electronics, Santa Clara, CA, USA) characterized by δλ = 0.40 nm. The inset depicts a magnified central zone showing the δλ parameter in detail.

Simulation and Results

VPIphotonic software is a simulation package that possesses a powerful simulation environment for the design, testing, and optimization of virtually any type of optical link. Its library is composed of active and passive integrated photonics devices, different types of optical sources and fibers, microwave electrical devices, optical and virtual electrical analyzers, among other things. It is very important to indicate that the parameters used in these simulations are adjusted to the optical characteristics of the devices available at our laboratory, considering that in the future it will be possible to carry out the experiments in the future. Unless otherwise indicated, the parameters L = 13 km and D = 16.75 ps/nm•km are maintained in the following. Figure 4 illustrates the layout implemented in VPIphotonics that corresponds to the scheme of Figure 2. For δλ = 1.61 nm, four passbands are observed at f 1 = 2.77 GHz, f 2 = 5.64 GHz, f 3 = 8.47 GHz, and f 4 = 11.26 GHz, as seen on the red curve in Figure 7. Its bandwidths at -3 dB are ∆f bp1 = ∆f bp2 = ∆f bp3 = ∆f bp4 = 488.28 MHz. Now, the length of the optical fiber is set to L = 25 km and a set of simulations are conducted again with the same parameters. Figure 8 shows the results for these cases. In the case of δλ = 0.82 nm (FSR = 100 GHz), three passbands are observed at f 1 = 2.92 GHz, f 2 = 5.87 GHz, and f 3 = 8.78 GHz (blue curve). For δλ = 1.61 nm (FSR = 200 GHz), seven passbands are observed at f 1 = 1.44 GHz, f 2 = 2.92 GHz, f 3 = 4.41 GHz, f 4 = 5.87 GHz, f 5 = 7.30 GHz, f 6 = 8.78 GHz, and an incipient f 7 = 10.25 GHz (red curve). In all of the cases, their measured bandwidths at -3 dB are ∆f bp = 253.9 MHz. The passband frequencies are listed together with those calculated by Equation ( 2) in Table 1. For the sake of comparison, the previous results are also listed. It is clearly shown in the table that the numerical and theoretical values are in a good agreement.

f 2 (GHz) f 3 (GHz) f 4 (GHz) f 5 (GHz) f 6 (GHz) f 7 (GHz) δλ = 0.

Conclusions

In this paper, we have successfully demonstrated, at the simulation level, a proposal of reconfigurability in the frequency of a passband MPF. This reconfigurability was based on the use of a tunable FPF. It was demonstrated that the modification of the FSR parameter of the FPF allows for changes to the intermodal separation (δλ) of the MLD used as an optical source. This proposal was supported by a series of simulations by means of the VPIphotonics software. Reconfigurability of the passband was validated in the frequency range of 0-12 GHz. Simulation results are in good agreement with the theoretical values computed by the use of Equation [START_REF] Capmany | A tutorial on microwave photonic filters[END_REF]. Furthermore, it was demonstrated that this proposal of reconfigurability can be easily adapted to other values of L and δλ parameters. It is remarkable that, as the length of the optical link was increased, the number of passbands was increased too. However, its electrical bandwidth was reduced. Nevertheless, this reduction is not critical because it is enough to code, for example, an analog signal in these passband windows for its transmission. As was established at the beginning of this paper, the filtering of passband windows is attractive in optical communication systems because they can be used as electrical carriers to transmit information coded on them. Now, with the obtained results, the versatility of the MPF analyzed in this work has been improved by selecting a particular passband window by appropriate manipulation of the FPF. As future work, once our laboratories open, we could corroborate our proposal with experimental results. For future experimental demonstration, we are going to use a tunable Fabry-Perot filter, which is thermally stable. In order to improve the frequency stability, it will be placed in a thermally insulated box. We will also study the frequency stability of our tunable MPF over time. In practice, the FSR will be controlled by the voltage applied to the piezoelectric electrode of the FPF. In summary, the main advantage of this proposal resides in the fact that the intermodal separation can readily be adjusted in a dynamical way by using a tunable FPF and, as a consequence, the reconfigurability in the frequency of a passband MPF. Furthermore, if the optical signal delivered by the photodetector is fed back by means of electrical components (passband filter and amplifier) to the electro-optic modulator, an optoelectronic oscillator configuration can be generated, and its tunability can also be achieved by means of the use of the Fabry-Perot filter.

Figure 1 .

 1 Figure 1. Periodic transmittance of an FPF.

Figure 2

 2 Figure2shows the basic scheme of the passband MPF used in this work. Essentially, it is formed by a multimode laser diode (MLD), an electro-optic Mach-Zehnder intensity modulator (MZ-IM), an optical link of single-mode standard fiber (SM-SF), and a fast photodetector (PD).

Figure 2 .

 2 Figure 2. Basic scheme of the passband MPF used in this work.

Figure 3 .

 3 Figure 3. Optical spectrum corresponding to the MLD NX7538BF-AA.

Figure 4 .

 4 Figure 4. Layout of the simulated MPF in VPIphotonics software to obtain the filtering of microwave signals.

Figure 5

 5 Figure 5 shows the result of simulation, considering δλ = 0.40 nm. Note the presence of a well-formed passband centered at f = 11.26 GHz exhibiting a bandwidth at -3 dB of ∆f bp = 488.28 MHz.

Figure 5 .

 5 Figure 5. Simulated frequency response of the MPF corresponding to δλ = 0.40 nm value.

Figure 6 .

 6 Figure 6. Layout of the simulated MPF in VPIphotonics software where a tunable FPF has been placed to modify the δλ parameter.

Figure 7

 7 Figure 7 shows the frequency response of the MPF with FSR = 100 GHz (δλ = 0.82 nm) in the blue curve and FSR = 200 GHz (δλ = 1.61 nm) in the red curve. In the case of δλ = 0.82 nm, a new passband f 1 = 5.64 GHz is observed; its value corresponds to a submultiple of the initial frequency of 11.26 GHz. Thus, two passbands at f 1 = 5.64 GHz and f 2 = 11.26 GHz are generated and their bandwidths at -3 dB are measured as ∆f bp1 = ∆f bp2 = 488.28 MHz.

Figure 7 .

 7 Figure 7. Simulated frequency response of the MPF corresponding to δλ = 0.82 nm (blue curve), and δλ = 1.61 nm (red curve).

Figure 8 .

 8 Figure 8. Simulated frequency response of the MPF corresponding to the use of δλ = 0.82 nm (blue curve), and δλ = 1.61 nm (red curve).

Table 1 .

 1 Summary of results.

		L = 13 km, D = 16.75 ps/nm•km
	-	f 1 (GHz)	f 2 (GHz)	f 3 (GHz)	f 4 (GHz)	-	-	-
	δλ = 0.40 nm	11.26	-	-	-	-	-	-
	Computed by Equation (2)	11.48	-	-	-	-	-	-
	δλ = 0.82 nm	5.64	11.26	-	-	-	-	-
	Computed by Equation (2)	5.60	11.20	-	-	-	-	-
	δλ = 1.61 nm	2.77	5.64	8.47	11.26	-	-	-
	Computed by Equation (2)	2.85	5.70	8.55	11.40	-	-	-
		L = 25 km, D = 16.75 ps/nm•km
	-	f 1 (GHz)

Acknowledgments: B.T.-M. wishes to thank to the Mexican Consejo Nacional de Ciencia y Tecnología (CONACyT) for the student scholarship number 846218. M.W.L. and I.E.Z.-H. thank the ECOS-Nord project (No. 299017).

Conflicts of Interest:

The authors declare no conflict of interest.


G Grams

R Somasundaram

J Margueron

S Reddy

Properties of the neutron star crust: Quantifying and correlating uncertainties with improved nuclear physics

A compressible liquid-drop model (CLDM) is used to correlate uncertainties associated with the properties of the neutron star (NS) crust with theoretical estimates of the uncertainties associated with the equation of state (EOS) of homogeneous neutron and nuclear matter. For the latter, we employ recent calculations based on Hamiltonians constructed using Chiral Effective Field theory (χEFT). Fits to experimental nuclear masses are employed to constrain the CLDM further, and we find that they disfavor some of the χEFT Hamiltonians. The CLDM allows us to study the complex interplay between bulk, surface, curvature, and Coulomb contributions, and their impact on the NS crust. It also reveals how the curvature energy alters the correlation between the surface energy and the bulk symmetry energy. Our analysis quantifies how the uncertainties associated with the EOS of homogeneous matter implies significant uncertainties for the composition of the crust, its proton fraction, and the volume fraction occupied by nuclei. We find that the finite-size effects impact the crust composition, but have a negligible effect on the net isospin asymmetry of matter. The isospin asymmetry is largely determined by the bulk properties and the isospin dependence of the surface energy. The most significant uncertainties associated with matter properties in the densest regions of the crust, the precise location of the crust-core transition, are found to be strongly correlated with uncertainties associated with the Hamiltonians. By adopting a unified model to describe the crust and the core of NSs, we tighten the correlation between their global properties such as their mass-radius relationship, moment of inertia, crust thickness, and tidal deformability with uncertainties associated with the nuclear Hamiltonians.

I. INTRODUCTION

The understanding of neutron star (NS) properties from fundamental nuclear physics inputs requires the precise determination of the relation between nuclear physics uncertainties and dense matter predictions. This became possible recently due to advances in theoretical efforts to predict properties of nuclei and dense nuclear matter, and advances in experimental nuclear physics that are now providing more stringent constraints. In addition, recent observations of NS radii by 1 NICER [1,2], and tidal deformabilities by the LIGO-Virgo collaboration [3] have also reached the accuracy to sharply constrain the dense matter EOS. These developments motivate the construction of models that can provide a unified description of the EOS of the crust and the core.

Our current understanding of NS crusts suggests that it is composed of finite nuclei, usually referred to as nuclear clusters since their properties are modified by the dense matter environment and differ from those of isolated nuclei probed in the laboratories, see for instance Ref. [START_REF] Haensel | Neutron stars 1: Equation of state and structure[END_REF]. The outer crust is dominated by the presence of an electron gas filling the whole volume in beta-equilibrium with neutrons and protons bound inside nuclear clusters. The nuclear symmetry energy controls the energy difference between neutrons and protons, and thus the isospin asymmetry inside the nuclear clusters. Electric charge neutrality is ensured by the presence of electrons, and its rapid increase with density favors the appearance of increasingly neutron-rich nuclear clusters [START_REF] Chamel | [END_REF]. From the experimental viewpoint, neutron-rich nuclei can be produced in nuclear facilities and provide strong constraints on the properties of nuclei present in the outer crust, see Refs. [6,7] for recent updates. Nuclear clusters in the inner crust are significantly more neutron-rich and coexist with a neutron fluid and their properties cannot be directly probed by experiments. Their description relies on theoretical models that are sensitive to the properties of bulk nuclear matter and the density and isospin dependence of the nuclear surface tension [8]. The neutron fluid in the inner crust is very likely to be in a superfluid state at low temperature and superfluidity is known to impact NS spin and thermal evolution Recently, several conceptual milestones have been reached in the prediction of neutron star matter from microscopic ab initio approaches which are based on realistic nuclear Hamiltonians constrained by nucleon scattering data. In particular, Hamiltonians derived using chiral effective field theory (χEFT) incorporate the symmetries of QCD and provides a systematic expansion of the operators in powers of the nucleon momenta. χEFT has two distinct features: (1) it consistently include three and higher-body interactions along with the two-body interactions that are well constrained by experiments, (2) it provides a robust method to estimate errors associ-ated with the truncation of the momentum expansion when the nucleon p is small compared to the breakdown scale Λ χEFT of the χEFT. Beginning with the pioneering work of Hebeler and Schwenk [9], several groups have used χEFT to predict the EOS of homogeneous neutronrich matter [10][START_REF] Hebeler | Annual Review of Nuclear and Particle[END_REF][12][13][14]. At sufficiently low densities, neutron matter is well understood because three-body interactions are small, and the two-body neutron-neutron interaction is strongly constrained by the neutron-neutron scattering phase shifts [15,16]. However, with increasing density and correspondingly larger nucleon momenta, higher-dimension operators including three-body forces begin play an increasingly important role. For both these reasons the associated uncertainty which can be estimated quite systematically grows. The densities up to which χEFT remains useful is still a matter of debate, current expectations are that it breaks-down between saturation density (n sat ) and twice-n sat [17][18][19].

Nuclear clusters in NS crust result from the equilibrium between attractive volume interaction, and repulsive surface and Coulomb interactions, at leading order. Despite recent progress in the description of finite nuclei based on chiral nuclear interaction, it is still computationally not feasible to calculate directly the properties of nuclear clusters in NS crust. Several well-motivated approximations could however be employed to predict and understand the properties of this complex system. Among them, the liquid-drop model (LDM) is a macroscopic approach that allows us to combine together the nuclear matter predictions from χEFT with finite-size (FS) terms generated from a leptodermous expansion of the total energy. The compressible liquid-drop model (CLDM) includes variations of the cluster density from one nucleus to another, through the density dependence of the bulk contribution to the total energy. The LDM and CLDM describe the collective degrees of freedom. They are different from the microscopic approach, typified by the Shell model or the energy-density functional approach, which centers around the single-particle degrees of freedom. Since the microscopic approach is the most general one, it contains the macroscopic one as an average.

In the leptodermous expansion of the total energy [20], the different contributions are sorted by decreasing powers of A 1/3 . The order A belongs to the domain of nuclear matter studies, as can be directly related to the metamodel approach for nuclear matter, the orders A 2/3 and A 1/3 characterizes the finite-size contributions, and finally, orders A 0 and below belong to the single-particle contributions, e.g. shell effects or pairing contribution, described by microscopic theories. In our study, we investigate several FS terms which are sorted by decreasing powers of A 1/3 as in the leptodermous expansion. In the following, we illustrate the fact that the leptodermous expansion provides an interesting scheme where the effects of different terms play a lower role as they contribute to a higher order. Our analysis is however stopped at a given order in the leptodermous expansion, hereafter called FS4, and we disregard higher order corrections for the present work. For instance, the neutron skin is not implemented and the in-medium surface modification of cluster energy is also currently disregarded. Some more microscopic effects are also not implemented, such as for instance the actual cluster density profile, non-uniform in Thomas-Fermi and Hartree-Fock approaches as well as quantum shell effects, see for instance the following Refs. [21,22]. We discuss the impact of such approximations in Sec. IV and VI. Being systematical in the implementation of the leptodermous expansion allows us however to compare our findings with previous ones, as well as to evaluate the impact of these different terms on the NS crust properties.

The impact of nuclear physics uncertainties on the NS crust has been analyzed in earlier work. Steiner [23], for instance, has constructed several NS crust using inputs from current experimental information while allowing exploration of the EOS uncertainties, in particular the one induced by the symmetry energy. Other approaches are constructed from currently available EOSs, which may not respect the low-density neutron matter expectations because they are fit to the properties of nuclei near saturation densities. In the present work, we explore both the nuclear experimental uncertainties from the knowledge of measured nuclear masses, as well as the theoretical uncertainties in the nuclear matter equation of state from many-body approach based on chiral NN and 3N interactions. This is first systematic investigation that accounts for these two sources of uncertainties. Our main findings, obtained by incorporating these two sources of information in our CLDM are:

• A suggested upper limit for the energy density of nuclear matter at saturation density: max sat ≈ -2.30 MeV fm -3 . Above this limit, our model can not equilibrate the bulk and FS terms over the nuclear chart.

• The CLDM disfavors some χEFT Hamiltonians, even if the FS terms in the CLDM vary independently to the bulk term while they are correlated from first principle. This is a strong rejection since in our modeling the FS terms are optimized to fit nuclear masses, while in reality they shall be fixed. So even if -by chance -nature chooses to fix these FS terms to be equal to our optimization, these models would still be rejected.

• The correlation between the surface energy (isoscalar and isovector properties) and the bulk symmetry energy depends on the considered FS model. The inclusion of the curvature contribution suggests however a typical value for the symmetry energy (≈ 32 MeV) where the model dependence is minimal.

We find that the present CLDM predicts the following properties for the neutron star crust:

• The crust composition (A cl , Z cl ) is mainly determined by the considered FS model, which -for most part of them -are controlled by the experimental nuclear energies. We find however that the cluster asymmetry I cl is much less impacted than A cl and Z cl .

• The crust composition (A cl , Z cl , I cl ), the proton fraction Y e , and the volume fraction u are essentially determined by the bulk contribution to the energy, fixing the nuclear cluster and the neutron fluid contributions. The bulk term is varied across the chiral Hamiltonians in our study.

• In the densest region of the crust both the Hamiltonians and the surface energy isospin parameter p surf have a dominant role in the determination of the crust-core properties and of the matter composition. These two properties are crucial for the determination of the crustal moment of inertia.

• While not negligible, the influence of the loss function in the fit to finite nuclei as well as of the effective mass is more sub-dominant.

Finally, we analyze the global NS properties and produce uncertainties for a few observables associated to light and canonical mass NS (1.0M and 1.4M ). We obtained the following results:

• Global properties, e.g. mass, radius, moment of inertia, tidal deformability, are critically determined by the chiral Hamiltonian properties in uniform matter.

• Other ingredients discussed here -FS terms, p surf , effective mass, loss function -play a much smaller role than the present uncertainties from astronomical observations, e.g. radius uncertainty from NICER or tidal deformability uncertainty from LIGO-Virgo gravitational wave detectors.

• Assuming the absence of phase transition in massive NS, we found that the present uncertainty in the nucleonic effective mass modifies the spherical NS maximum mass, M TOV , by about 0.15M at maximum.

The paper is organized as follows. Sec. II is reserved to homogeneous matter where the meta-model and the fit to χEFT are described. We detail the CLDM used to describe finite-size effects on the clusterized matter in Sec. III. We compute and quantify uncertainties on the NS crust in Sec. IV, while details on NS macroscopic properties are discussed on Sec. V. Finally, the conclusions of the present work are drawn on Sec. VI.

II. HOMOGENEOUS MATTER

For homogeneous matter, we consider the six Hamiltonians, H1-H5 and H7 (H6 being disregarded for not fitting well the binding energy of 3 He), which have been generated by many-body perturbation theory (MBPT) based on chiral NN and 3N interactions [24] and the two recent χ EF T predictions from Ref. [25]: DHS L59 and DHS L69 . These eight predictions for nuclear matter are used to calibrate a set of eight nuclear meta-models (MM), and the version we consider here is a small extension of the original one [26]. This extension has already been presented in Ref. [27]. The low density correction to the energy is now controlled by the function b(δ) = b sat + b sym δ 2 instead of the parameter b.

Let us briefly summarize the main ingredients of the MM. The energy density is the sum of a kinetic and potential terms, e MM = t * + e pot , where the kinetic term reads,

t * (n, δ) = t sat 2 n n sat 2/3 m m * n (n, δ) (1 + δ) 5/3 + m m * p (n, δ) (1 -δ) 5/3 , (1)
with t sat = 3 2 k 2 F /(10m) and given τ 3 = 1(-1) for neutrons (protons), the effective mass reads

m m * τ (n, δ) = 1 + κ sat n sat + τ 3 δ κ sym n sat n. (2)
In the present work we explore the effect of including or not the effective mass on the kinetic energy. We show results with the effective mass equal to the bare mass m * τ = m τ , and also with the effective mass given by Eq. (2). In order to obtain the parameters κ sat and κ sym for the Hamiltonians we consider here, we first derive the density dependent Landau effective mass for neutrons from its single-particle spectrum, as it was done in Ref. [27]. Then Eq. ( 2), with τ = n, was fit to this quantity in symmetric matter (SM) and neutron matter (NM) for each individual Hamiltonian. As far as the finite residuals of the fits are concerned, the uncertainties on the parameters κ sat and κ sym are negligible for all Hamiltonians. There is however a spread in the results of the fit across different Hamiltonians reflecting an intrinsic uncertainty in the χEFT predictions. Ref. [27] performed a Bayesian quantification of this uncertainty, whereas in this work we probe the χEFT uncertainty by making distinct predictions for each individual Hamiltonian and then monitoring the dispersion among these predictions.

The isospin splitting of the effective mass is defined as,

∆m * sat = m * n (n = n sat , δ = 1) -m * p (n = n sat , δ = 1) = -2κ sym (1 + κ sat ) 2 -κ 2 sym m N , (3)
where we fixed m n = m p = m N . Sometimes, the proton effective mass in NM is not calculated from microscopic approaches, but SM and NM calculations are performed.

In these cases, it is interesting to consider the difference between the neutron effective masses in NM and SM defined as,

Dm * sat = m * n (n = n sat , δ = 1) -m * n (n = n sat , δ = 0). ( 4
)
The effective mass, the splitting of the effective mass ∆m * sat and Dm * sat for the eight Hamiltonians are shown in Table I. We find that Dm * sat , which is obtained directly from the MBPT single-particle spectrum as explained before, is in good agreement with our previous study in Ref. [27]. The splitting of the effective mass can be inferred using the MM parametrization of the effective mass, see Eq. (2). We found a very good correlation between Dm * sat and ∆m * sat , suggesting that ∆m * sat is about 50% larger than Dm * sat . Finally, for DHS L59 and DHS L69 we only consider the bare mass case.

In the MM, the potential term is expressed as,

e pot (n, δ) = N j=0 1 j! v sat,j + v sym2,j δ 2 x j +e pot,low-n (5)
with the low-density correction expressed as

e pot,low-n = v low-n (δ) x N +1 e -b(δ)n/n emp sat , (6)
with v low-n (δ) = v low-n sat + v low-n sym δ 2 . The parameters v low-n sat and v low-n sym are fixed by the condition that e pot (n = 0, δ) = 0 for δ = 0 and δ = 1, and the parameters b sat and b sym are adjusted to reproduce the very low density dependence of the MBPT calculations, see Ref. [27] for more details. With the MM, the binding energy e MM (n, δ) can be obtained for any arbitrary value of the density n and the isospin asymmetry parameter δ. There are two interesting limits which are SM and NM defined as e SM (n) = e MM (n, δ = 0) and e NM (n) = e MM (n, δ = 1).

Table II shows the residual χ 2 values for the fit, where the χ 2 loss function is defined as

χ 2 = 1 2 i e data,i -e MM,i σ i 2 , (7)
where σ i was taken to be a 10% uncertainty on the data, due to imperfect modeling [28]. The fit was performed by minimizing the χ 2 with data in 0.4 < k F n < 1.6 by using the standard Levenberg-Marquardt algorithm implemented in python's scipy package. The impact of the effective mass is very small. At most it improves the reduced χ 2 /N by 10% in SM, no effect in NM. In addition, in the density region out of the data, the impact of the effective mass is also very small. This will also be confirmed in Sec. V where the impact of the effective mass on the mass-radius relation will be presented. In Table II, for Hamiltonians H1-H5 and H7, the number of data points N is 11 [START_REF] Hebeler | Annual Review of Nuclear and Particle[END_REF] for SM (NM) when 0.4 < k F n < 1.0 and 22 (24) for SM (NM) when 0.4 < k F n < 1.6. For the other two Hamiltonians, N is 17 (9) in SM (NM) when 0.4 < k F n < 1.6 and 2 (-) for SM (NM) when 0.4 < k F n < 1.0. A detail comparison of the MBPT calculations (dots) and the MM fit (lines) is shown in Fig. 1, the well known SLy4 Skyrme model prediction is also shown for reference. For the bare mass case (top panels) we also represent the recent MBPT calculations DHS L59 and DHS L69 . The three models DHS L59 , DHS L69 and H7 present stiffer NM energy at high densities compared to H1-H5. There is a very good agreement between the data and the MM down to k F = 0.4 fm -1 (in density n ≈ 0.004 fm -3 ), as shown in panels (b) and (d) for instance. The vertical gray band shows the region where the fit is less accurate. In the following, the MM calibrated on the Hamiltonian Hn will be labelled as Hn MM .

Note in Fig. 1(b) that the MM reproduces very well the NM energies as predicted by χEFT, at variance with SLy4 which overestimates the energy per particle at low density. This is indeed a general feature of Skyrme interactions. We recently analyzed the impact of this systematical differences between χEFT and Skyrme SLy4 in low density NM on the crust EOS within the CLDM [29]. We found that some observables are very sensitive to these differences, e.g. energy density, pressure, sound speed, while other are much less impacted, e.g. cluster configuration (A cl , Z cl , I cl ), which are mostly determined by experimental nuclear masses. Having a good description of NM as predicted by χEFT is however important to predict NS crust properties.

The symmetry energy is defined as the energy difference between NM and SM, and the quadratic contribution to the symmetry energy reads,

e sym (n) = e NM (n) -e SM (n) , (8)
e sym,2 (n) = 1 2 ∂ 2 e(n, δ) ∂δ 2 δ=0 , (9)
The topological properties of the energy per particle around saturation density, with empirical expectation n emp sat ≈ 0.155(5) fm -3 , are encoded into the nuclear empirical parameters (NEP), e.g. E sat , E sym , E sym,2 , which are defined as,

e SM (n) = E sat + 1 2 K sat x 2 + 1 6 Q sat x 3 + 1 24 Z sat x 4 + . . . , (10)
e sym (n) = E sym + L sym x + 1 2 K sym x 2 + 1 6 Q sym x 3 + 1 24 Z sym x 4 + . . . , (11)
e sym,2 (n) = E sym,2 + L sym,2 x + 1 2 K sym,2 x 2 + 1 6 Q sym,2 x 3 + 1 24 Z sym,2 x 4 + . . . , (12)
where the density expansion parameter is defined as x = (n -n sat )/(3n sat ).

The low order NEP of the 8 Hamiltonians are given in Tab. III. The values are rounded and the uncertainties in these parameters are smaller than the rounding. These uncertainties are thus not given. The dispersion of the NEPs across the different Hamiltonians however serves to capture uncertainties intrinsic to the χEFT expansion. The saturation energy and density of H2 MM -H4 MM agree reasonably well with the empirical ones determined from experimental data, see for instance Ref. [26], while the saturation energy of H5 MM and H7 MM is higher than the expected one. The saturation density of H7 MM is also quite lower than the empirical one. For H1 MM the saturation density is higher than the expected one, and the saturation energy is lower. In the following, we will confirm that the fit to the experimental masses is poorer for H5 MM and H7 MM compared to the other Hamiltonians, as we can already anticipate from their empirical properties.

It is well known that the high order empirical parameters Q sat/sym and Z sat/sym are not constrained by χ EF T calculations in uniform matter, since χ EF T is limited to low densities, n < 2n sat [27]. In the present study, we fix the following values for all Hamiltonians:

Q sat = -220 MeV, Z sat = -200 MeV, Q sym = 700
MeV, and Z sym = 500 MeV. This choice allows all the considered Hamiltonians to reach the astrophysical constraint related to the observed maximum mass of NSs, which is 2.0(1) M [30,31]. The recent measurement for PSR J0740+6620, suggesting M TOV ≥ 2.14 (10) M is also compatible with this constraint [32].

The symmetry energy e sym (n) is shown in Fig. 2. Left panel shows in light (dark) blue band the constrains from isobaric analog state IAS (IAS + neutron skin, ∆r np ), from Ref. [33]. In yellow are shown the PREX-II predictions for the symmetry energy, where we vary E sym = 38.1 ± 4.7 MeV and L sym = 106 ± 37 MeV as suggested by the publication [34]. The symbols represent the MBPT calculations [24,25] as in Fig. 1. Right panel shows the symmetry energy, e sym , normalized by the free fermi gas symmetry energy. Note the higher value for the symmetry energy predicted by SLy4 compared to the Hamiltonians at low densities. This is a direct effect of the high NM energy predicted by SLy4, as seen in Fig. 1(b) and(d). The solid lines shows the MM with the bare mass while the dashed lines include the corrections induced by the effective mass, see Eqs. ( 1) and (2). Note that the effect of the effective mass is very small on these curves.

There is a disagreement between χ EF T predictions and the recent constrain from PREX-II [34] for the symmetry energy and its density dependence, as shown in Fig. 2. To reproduce PREX-II predictions not only a large value for L sym is necessary, like DHS L59 and DHS L69 , but also for E sym .

The β-equilibrium with e and µ satisfies the following equations,

µ n = µ p + µ e , µ e = µ µ . ( 13
)
The total pressure, including baryon and lepton contributions, at β-equilibrium is shown in Fig. 3 for the eight Hamiltonians and SLy4. The prediction by the LIGO-Virgo collaboration inferred from GW170817 for the pressure at 2 n sat [35] is also shown. Fig. 3 shows that the eight models are in well agreement with GW170817 constrain. The impact of the effective mass m * (dashed lines versus solid lines with the bare mass) remains small for the pressure. In summary, we have shown that the MM can accurately reproduce the MBPT predictions for each of the considered eight Hamiltonians. In addition, we impose a prescription for the high order empirical parameters that allows the present extension of the MBPT to reach the astrophysical constraint for the TOV mass. We have also checked that all the present Hamiltonians are consistent with the inferred pressure at 2 n sat by the LIGO-Virgo collaboration.

III. THE COMPRESSIBLE LIQUID-DROP MODEL FOR FINITE NUCLEAR SYSTEMS

Finite nuclear systems results from the equilibrium between the bulk attraction, as seen in homogeneous matter, and the surface repulsion, which originates mostly from the surface tension and the Coulomb repulsion. In the present study, we consider various extensions of the compressible liquid-drop model (CLDM), see for instance the seminal BBP model [36], which can be justified from the leptodermous expansion [20].

In the crust of NS, we consider the following composition: the nuclear clusters are composed of neutrons and protons, which are described by A cl and I cl , being the mass number and the isospin asymmetry. The neutron and proton particle numbers in the nuclear clusters are obtained from:

N cl = A cl (1 + I cl )/2 and Z cl = A cl (1-I cl )/2.
In addition, clusters are embedded in electron and neutron gases, described by their uniformly distributed densities n e and n ng . We have implicitly assumed the r-representation [37] for the Wigner-Seitz cell. In this representation, the particles in the cluster volume V cl are in equilibrium with the ones in the outside volume V WS -V cl , where V WS is the Wigner-Seitz volume. There are therefore five variables in total (four particles and one volume), but the actual variables can be any combination of these variables. In the present study, we consider the following ones: A cl , I cl , n cl , n e and n ng , where n cl = A cl /V cl as in Ref. [38].

The total cluster energy in NS crust is expressed as the sum of the independent contributions from the clusters E cl , the electrons E e and the neutron gas E ng ,

E tot (A cl , I cl , n cl , n e , n ng ) = E cl (A cl , I cl , n cl ) + E e (n e ) +E ng (n ng ) . (14)
The cluster contribution in the CLDM is expressed as a bulk energy contribution, determined from homogeneous matter, and a finite-size contribution, including Coulomb, surface, curvature terms at leading orders. The cluster binding energy contributing to Eq. ( 14) is given by,

E cl (A cl , I cl , n cl , n e ) = E bulk (I cl , n cl )+E FS (A cl , I cl , n cl , n e ) , (15)
n cl being the cluster density. Note that in the CLDM the cluster density n cl is fixed to be uniform, at variance with the droplet model [39], Thomas-Fermi approaches [40,[START_REF] Haensel | Physics of Neutron Star Interiors[END_REF] or microscopic models [21,22,[START_REF] Negele | [END_REF]43] where it decreases at the surface and depends on the radial coordinate r as n cl (r). The non-uniformities are however incorporated in a effective way by adjusting the FS terms to the experimental nuclear masses. Also, the neutron and proton radii are identically equal to the cluster radius R cl (there is no neutron skin in the present model).

In the NS crust, the total density n B = (A cl +N g )/V WS is further imposed, contributing to fix one constraint among the five independent variables. This constraint is treated with the Lagrange multiplier technique, as suggested in Ref. [44] for instance.

The energy of an isolated nucleus such as the ones present on Earth is simply defined as,

E nuc (A cl , I cl , n cl ) = E bulk (A cl , I cl , n cl ) +E FS (A cl , I cl , n cl , n e = 0) . ( 16
)
Note that there are only three independent variables in this case. There are no electron and neutron gases surrounding the nuclear cluster. In addition, note that the Wigner-Seitz volume V WS is undefined (it is indeed infinite for isolated nuclei) but the cluster volume V cl is.

In the present CLDM, the global asymmetry of the cluster I cl = (N cl -Z cl )/A cl coincides with the cluster bulk asymmetry δ cl = (n cl,n -n cl,p )/n B since the neutron or proton skin are not considered here. Note however that neutron skin has been considered in Ref. [23] by introducing a fit parameter ζ relating I cl and δ cl , as δ cl = ζI cl . If ζ is unity, there is no skin, while if ζ < 1, then all nuclei with N cl > Z cl will have a neutron skin (R n > R p ). In reality, the parameter ζ is function of A cl and Z cl , see Ref. [45], as well as of the nuclear interaction, as illustrated by the correlation between the neutron skin in 208 Pb and the slope of the symmetry energy L sym [46,47]. It is thus a strong approximation to impose the relation δ cl = ζI cl with ζ constant that we prefer not to consider here. The skin contribution modifies the Coulomb term by using the proton radius instead of the cluster radius, which modified the Coulomb energy by a factor proportional to the difference between the bulk asymmetry δ cl and the global asymmetry I cl . This term increases with the nuclear size and asymmetry and can be important for nuclear clusters present close to the crust-core transition. However, this modification is small compared to the leading order terms considered here and we follow the procedure of recent works [38,48] and neglect the presence of neutron skin in the present study. In a future development, a consistent derivation in the spirit of Ref. [45] for instance will be considered.

The CLDM we consider is comparable to the pioneering BBP model [36] and is well suited to analyze the origin of the uncertainties in the predictions of the NS crust. More microscopic models for the crust have indeed been developed, see for instance Ref. [START_REF] Negele | [END_REF] and recent efforts in Refs. [22]. While being less accurate than microscopic models in reproducing finite nuclei, the present CLDM allows us a better understanding of the various features influencing the properties of the NS crust, which are difficult to analyze in a microscopic model. There are however missing features, such as shell and pairing effects, but these feature are sub-dominant in the leptodermous expansion: they represent a refinement in the description of experimental binding energies which is of the order of a few MeV in total energy, compared to the leading order contributions which are of the order of hundred of MeV. The uncertainties are indeed still large at the leading order, as we will see in the following.

The contributions to the energy E bulk , E FS , E e and E ng will be detailed in the following sub-sections.

A. The cluster bulk contribution

The cluster bulk contribution to the energy per particle is the leading order term in the leptodermous expansion (order A cl , the mass term). It is related to the homogeneous matter calculation, represented here by the MM energy density MM (n n , n p ) = n B e MM (n n , n p ), given by

e bulk (I cl , n cl ) = 1 n B MM (n n , n p ) -n n m n c 2 -n p m p c 2 , (17)

B. The finite-size contribution

The finite-size term in Eq. ( 16) incorporates the nuclear contributions to the cluster energy at all orders in the leptodermous expansion. In the present study, we limit ourself to the leading order terms: the Coulomb term is in

Z 2 cl /A 1/3 cl ≈ A 5/3 cl
(dominant term at large A which prevents super-heavy nuclei to exist), the surface is in A

E FS (A cl , I cl , n cl ) = E Coul (A cl , I cl ) + E surf (A cl , I cl ) +E curv (A cl , I cl ), (18)
where the Coulomb term for a spherical and uniform distribution of protons is given by the direct and exchange contributions,

E Coul = C Coul (E Coul,Dir + E Coul,Ex ) , (19)
with

E Coul,Dir = 3 5 Z 2 cl e 2 R p f Coul (u) (20) = a c 1 -I cl 2 2 f Coul (u)A 5/3 (21)
E Coul,Ex = - 3 4 3 2π 2/3 Z 4/3 cl e 2 R p h Coul (u) (22)
= -5a c 4 3 2π

2/3 1 -I cl 2 4/3 h Coul (u)A , ( 23
)
where e 2 ≈ c/137 and the functions f Coul and h Coul are defined as:

f Coul (u) = 1 -(3/2)u 1/3 + (1/2)u and h Coul (u) = 1 + u 1/3
, with u the volume fraction of the cluster, defined as,

u = V cl V WS = n e n cl,p = 2n e (1 -I cl )n cl , (24)
where n cl,p = n cl (1 -I cl )/2 is defined as the proton density in the cluster. In the function f Coul , the first term corresponds to the proton-proton repulsive interaction, the second term is the "lattice contribution" including the electron-proton and electron-electron interaction, under the hypothesis of a globally neutral Wigner-Seitz cell. Then the third term in f Coul stands for the finitesize correction which becomes important when the cluster volume is comparable with the Wigner-Seitz volume. This term is important near the crust-core transition and pushes the transition to nuclear matter towards higher densities. The first term is the only one remaining in the case of isolated nuclei, corresponding to the limit u = 0. Since there is no proton-electron contributions to the exchange Coulomb energy, the first (second) term in h Coul corresponds to the proton-proton (electron-electron) contribution.

The coefficient C Coul in Eq. ( 19) is a variational parameter which is fine tuned over the nuclear mass table. It describes -in an effective way -the effect the diffusive nuclear surface on the Coulomb energy, which is neglected in the sharp drop off density profile that we consider here. Since the diffusive surface is expected to be a small correction, the fit value is expected to remain close to 1, C Coul ≈ 1.

Note that the direct Coulomb term scales like A most of nuclei in the nuclear chart, the Coulomb interaction remains small compared to the nuclear one. For this reason, the bulk term in A cl and the Coulomb direct contribution in A 5/3 cl are considered at the same order in the leptodermous expansion. It is also interesting to note that the exchange Coulomb term contributes to two order lower compared to the direct term in the A 1/3 cl leptodermous expansion. It is thus expected to effectively contribute at one order below the curvature contribution.

Neglecting the difference between the neutron and proton radii -no skin approximation -, we have

R p = R cl = r cl (n cl )A 1/3 cl with r 3 cl (n cl ) = 3/(4πn cl )
, and the Coulomb factor reads,

a c (n cl ) = 3 5 4π 3 n cl 1/3 e 2 . ( 25
)
Note that the Coulomb factor a c defined from Eq. ( 25) depends on the cluster density. The Coulomb parameter, a c , is however often taken as a constant, see for instance Ref. [38], either as a free parameter to be fitted or as function of the constant n sat . The different assumptions for a c give differences on the description of isolated nuclei. Note that if a c is taken to be constant (often taken to be of the order ∼ 0.7 MeV), neither the Coulomb nor the surface energy contributes to nuclear pressure, the pressure derives from the bulk term only. In NS crust however, the Coulomb term contributes to the pressure thanks to its dependence in the volume fraction u. We show the difference of having a c constant or not in the sequence of the paper, FS1 refers to a c = a c (n sat ), while in FS2 and others we have a c = a c (n cl ).

The surface energy is proportional to the surface tension σ surf (I cl ) and scales as

A 2/3 cl . It reads E surf (A cl , I cl , n cl ) = 4πR 2 cl σ surf (I cl ) (26) = 4πr 2 cl σ surf (I cl )A 2/3 cl , (27)
with σ surf (I cl ) as expressed, as suggested in [49], as

σ surf (I cl ) = σ surf,sat 2 p surf +1 + b surf Y -p surf p + b surf + (1 -Y p ) -p surf , (28)
where Y p = Z cl /A cl = (1 -I cl )/2 is the cluster proton fraction and σ surf,sat is a parameter that determines the surface tension in symmetric nuclei. The parameter p surf entering into the expression of the surface tension (28) plays an important role at large isospin asymmetries. It is usually fixed to be p surf = 3 since the seminal contribution [49], but a small variation around 3 plays an important role at large asymmetries, which occurs around the core-crust transition densities in NSs [48].

For small asymmetries,

σ surf (I cl ) ≈ σ surf,sat -σ surf,sym I 2 cl ( 29
)
with

σ surf,sym = σ surf,sat 2 p surf p surf (p surf + 1) 2 p surf +1 + b surf . ( 30
)
One can thus relate the parameter b surf to the surface symmetry energy σ surf,sym . We have

b surf = 2 p surf p surf (p surf + 1) σ surf,sat σ surf,sym -2 . ( 31
)
In the following, we prefer to use the parameter σ surf,sym instead of b surf , since it directly reflects the isospin dependence of the surface tension for small isospin asymmetries, as shown in Eq. (29). For this reason the domain of variation of σ surf,sym is better constrained than the one for the parameter b surf , which ease the determination of the prior for this parameter, see section III C. The surface energy parametrization (29) do not include a correction due to the presence of the neutron gas, while it has been adjusted on Hartree-Fock microscopic calculations including this effect [50]. The authors of this parametrization already realized this "paradox" and concluded that the contribution of the neutron gas and that induced by the change of the surface density as clusters get more and more neutron rich play an opposite role.

The in-medium effects (neutron gas and nuclear surface modifications) can however be important for the description of the crust-core transition, as suggested from the following authors [21,22,[START_REF] Haensel | Physics of Neutron Star Interiors[END_REF]. Microscopic calculations shall be used to further refine the functional dependence of the the surface energy term, like in [51,52] for instance. Tab. IV suggests standard values for the FS parameters obtained by averaging over the usual parameters, see Ref. [48] for a sample of these parameters associated to various Skyrme interactions.

The curvature energy is controlled by the curvature tension σ curv (I), and follows [53],

E curv (A cl , n cl , I cl ) = 8πr cl σ curv (I cl )A 1/3 cl , (32)
with

σ curv (I cl ) = ασ curv,sat σ surf (I cl ) σ surf,sat β curv - 1 -I cl 2 . ( 33
)
The parameter α is fixed to be α = 5.5, since we allows the variation of the parameter σ curv,sat in the fit to the binding energy over the nuclear chart. The standard values for the curvature parameters σ curv,sat and β curv are also given in Tab. IV.

In the present paper, we explore the role of various approximations in the FS terms on the NS crust properties. We order these approximations by their expected impact and selected four of them, that we call FS1 to FS4, see Tab. V:

• FS1 is the simplest approximation for the FS term that we consider, where only the surface and direct Coulomb contributions are included and for which the Coulomb and surface parameters, which are related to r cl , are taken constant and fixed by setting n cl = n emp sat . • FS2 is an improved version of FS1, where the Coulomb and surface parameters are varied and fixed from the actual value of the cluster density n cl .

• In FS3, we additionally incorporate the effect of the curvature contribution.

• Finally in FS4, we add the exchange Coulomb contribution to the Coulomb energy term.

In the future, we plan to incorporate more contributions and go beyond FS4. The approximations captured into the FS1-4 models represent however a systematical development where the refinements are expected to play a smaller and smaller role. This will be confirmed in Sec. IV. The FS1-4 models allows also to understand the differences between the models proposed for the crust. In their seminal paper in 1971, Bethe, Baym, and Pethick [54] have introduced the first version of the CLDM with Coulomb and surface terms only, similar to our FS2 with the addition of neutron skin effect and without optimization to the nuclear chart. Another well-known model for the crust was proposed by Douchin and Haensel in 2001 [55], still considered a CLDM model with the additional contribution of the curvature term, as in our model named FS3. The surface and curvature terms were however determined from many-body methods [56]. They have also considered different geometries in the pasta phase and incorporated the effect of skins in the CLDM, which go beyond the present approximations. In 2008, Steiner considered a CLDM [23] with surface and Coulomb terms comparable to our FS2 approximation, but introduced in addition an effective way to describe skins that we have previously discussed. Newton et al. suggested in 2012 [53] a CLDM with a full Coulomb term (direct and exchange) as well as surface and curvature contributions, as our FS4 approximation. They additionally studied different geometries in the core-crust region (pasta phases). In 2016, Tews adopted the same CLDM as the one suggested by Steiner [23] and he additionally considered a bulk contribution determined by chiral EFT calculations [57]. The FS terms were also adjusted to reproduce the binding energy over the nuclear chart. These two features make this study comparable to our model for the crust. In 2017, Viñas et al. [40] extended the CLDM by considering the Thomas-Fermi approximation and introduced surface, Coulomb and curvature terms, all FS terms being function of the proton radius. This treatment is similar to our FS4 approximation, but with additional contribution due to skin and pasta phases. In 2019, Carreau et al. use a CLDM with surface and Coulomb terms in [38], which can be compared to our FS1 approximation, and then included in 2020 curvature and shell effects [58], which can be compared to our FS3, but fixing a c constant and including shell effects.

Fantina et al. [21] constructed a unified EOS within Brussels-Montreal Skyrme (BSk) functionals, using experimental nuclear masses at the outer-crust, with the Brussels-Montreal mass model when there is no data, and a extended Thomas-Fermi plus Strutinsky integral method at the inner crust [59]. This model is an interesting alternative to describe neutron star crust, where it is possible to model shell effects. Recently this approach inspired a new model combining the advantage of CLDM with the modelling of shell effects [58].

In conclusion, there is not a strict equivalence between the FS1-4 approximations we suggest and the various models for the crust which have been investigated, but we believe our series of approximation is incremental and well suited to the understanding of the role of various terms, and associated uncertainties, on the properties of the crust of NS. The comparison to other crust models suggests also that our FS approximation series could be extended in order to incorporate neutron skins, shell effects, and different geometries in the core-crust region.

C. Fit to experimental nuclear masses

In this section we discuss the impact of the fit on isolated nuclei. We analyse two definitions for the loss function, ∆ E or ∆ E/A , which are defined as: where X = E or X = E/A: X i exp is the experimental values and X i nuc is the energy given by the CLDM model. In the present fit, we have considered N = 3375 nuclei taken from the 2016 Atomic Mass Evaluation (AME) [60], with nuclei in the range: A = 12 -295 and Z = 6 -118.

∆ X = 1 N N i=1 (X i exp -X i nuc ) 2 1/2 , (34)
It is interesting to analyze the overall impact of the fine tuning of the mass models on the surface and curvature contributions to the nuclear energy. To do so, we represent in Fig. 4(a) the following quantity (e surf + e curv )/(4πr 2 sat A 2/3 cl ) as function of the nuclear asymmetry I cl for a typical value A cl = 100. The homogeneous contribution is also fixed to be given by SLy4 MM . The different colors shows the result of the different FS terms, as indicated in the legend, the solid lines correspond to the standard values, given in Tab. IV, while the dashed lines results from the minimization of the loss function ∆ E/A . The overall trend is similar for all the considered cases. The contribution of the curvature term for FS1 and FS2 is null. The difference between FS1 and FS2 reflects the role played by the parameter r 2 cl in FS2 while it is fixed to be r 2 sat in FS1: since the cluster density n cl decreases as I cl increases, r 2 cl becomes larger than r 2 sat at large I cl , as shown in Fig. 4(a). The consistent treatment of the cluster radius in the FS terms, included in FS2 and further approximation, tends to increase the surface term at large asymmetry I cl compared to FS1, where the values radius parameter is constant.

Employing the standard parameter set (solid lines), FS3 and FS4 predict larger values than FS1 and FS2, since the contribution of the curvature e curv adds up to the surface one e surf . The exchange Coulomb term has no effect on the surface tension, therefore FS3 and FS4 shows identical curves in the figure. The dispersion between the dotted lines is smaller than between the solid lines as a result of the minimization which is performed for each FS model. Finally, in the fit the curvature contribution is absorbed by a reduction of the surface one, such that the sum remains identical, as seen from the comparison of FS3 to FS4 (dashed lines). The fit with the Coulomb exchange in FS4, absent in FS3, tends to slightly increase the surface and curvature terms around I cl ≈ 0. Fig. 4(b) shows the influence of the surface parameter p surf on the same quantity as in Fig. 4(a), where we show only FS2 and FS3 for clarity. Around I cl ≈ 0 the impact of varying p surf is null, while it plays an major role for I cl > 0.3, as it was already remarked in Ref. [38]. As a consequence p surf cannot be determined from the confrontation to the experimental nuclear chart [38]. It however plays an important role in the densest layers of the NS crust at the vicinity of the core, where matter is the most neutron rich. The value of the parameter p surf is thus an important source of uncertainties which cannot be controlled by the nuclear experimental data. This will be illustrated in Sec. IV D.

We now vary both the bulk and the FS terms, where the bulk terms we consider are the ones which reproduce the MBPT predictions based on chiral interactions, see Sec. II. The impact of varying the bulk contribution on the minimization based on either the two loss functions ∆ E and ∆ E/A is shown in Fig. 5. The horizontal axis is chosen to be sat = n sat E sat and the different colors represent the best fit obtained for each FS approximation (FS1 to FS4). The solid (dashed) lines are obtained with the bare (effective) nucleon mass. There is a small improvement using the effective mass instead of the bare mass, but it however remains small compared with the impact of changing the Hamiltonian.

For instance, the Hamiltonians for which sat > -2.30 MeV fm -3 are less good in reproducing the experimental nuclear masses than the others. The reason is that the reproduction of the experimental binding energies requires a delicate balance between the attractive bulk term sat and the repulsive FS term FS . As sat increases, the FS term decreases such that the sum remains constant. The compensation can happen until the FS term becomes almost zero. Above this limit, if sat continues to increase the energy density becomes too large to be able to reproduce the experimental masses, and the quality of the fit gets worst and worst, as illustrated in Fig. 5.

The groups FS3-4 -where the curvature term has been incorporated -reduce the loss function compared to the groups FS1-2. The best models are H1, H2, H3, and H4, among which H2 and H3 are even better. In the following analysis, the predictions based on H1-H4 will be marked with a light grey band, while the ones based on H2-H3 will be identified with a darker grey band, as illustrated in Fig. 5.

We also analyse in Fig. 5 the CLDM based on the Hamiltonians DHS. Note that it is still sat which drives the goodness of the model: the closer it is to -2.6/-2.8 MeV fm -3 the better is the agreement with experiments. There is however a reduction of the goodness of these models as the value of L sym departs from the one of the other models by 10 to 15 MeV, see Table III. Since the DHS Hamiltonians do not show a good reproduction of experimental masses, we neglect them on the NS crust study. In the next section we show the Hamiltonians which better reproduced nuclear masses, H1-H4, in gray band together with H5 and H7 for comparison.

Finally, we explore the correlation between the surface energy parameters σ surf,sat and σ surf,sym and the symmetry energy E sym for the Hamiltonians H1-H7, the FS models FS1-FS4, the prescription of the nucleon mass, and the loss function used for the minimization, as shown in Fig. 6. Having or not the curvature contribution in the CLDM is the main source of difference between these correlations. With the contribution of the curvature term the parameters σ surf,sat and σ surf,sym are almost linearly correlated with E sym . The correlation between σ surf,sym and E sym was already discussed in Ref. [23], but we note here that this correlation is model dependent: without the curvature term it become an anti-correlation, which can even be almost flat if the minimization of the mass table is based on the loss function ∆ E/A .

It is further interesting to note that even if these correlations are model dependent, they are crossing for the value of the symmetry energy around 32 MeV, comparable to the one of the Hamiltonians H2 and H3. Therefore, selecting H2-H3 Hamiltonians reduces strongly the model dependence of the parameters σ surf,sat and σ surf,sym . In the future, it will be interesting to check if this property remains while adding more terms to the FS approximation.

The conclusion of this confrontation of the CLDM to the experimental nuclear masses is that the Hamiltonians H1-H4 which satisfy the condition sat < -2.30 MeV fm -3 reproduce well the experimental masses over the mass table. A better reproduction over the mass table is obtained for the Hamiltonians for which sat = -2.70 (20) MeV fm -3 . We obtained the best results for the confrontation of the CLDM to the experimental masses for the Hamiltonians H2 and H3, which will represent in the following the best models for the NS crust properties.

D. The electron and neutron gas contributions

At variance with isolated nuclei, the energetics of the NS crust incorporates the contributions from the electron and neutron gas, that we briefly present here for completeness.

Similarly to the contribution to the bulk properties for the nuclear clusters, the energy of the neutron gas is given by the MM as,

e ng (n ng ) = 1 n ng MM (n ng = n n , n p = 0) . (35)
Note that by considering the same model providing the core properties and the bulk and gas contribution in the crust, our CLDM provides a unified description of the entire NS [21,55,61].

Since the electron interaction between electrons and between electrons and protons have already been absorbed in the Coulomb term, the remaining electron gas contribution is purely kinetic and reads for a relativistic gas,

e = C e x e (1 + 2x 2 e ) 1 + x 2 e -asinh(x e ) . ( 36
)
where x e = c k Fe /m e c 2 , the electron Fermi momentum k F e = (3π 2 n e ) 1 3 , and the overall constant C e = (m e c 2 ) 4 /(8π 2 ( c) 3 ). The electron chemical potential is expressed as,

µ e = m e c 2 1 + x 2 e , (37)
and the pressure,

P e = -e + n e µ e . (38)

IV. NS CRUST PROPERTIES

In this section, we derive the equilibrium configurations in the crust of NS described by the CLDM (15). We first derive the equilibrium equations and then present and discuss our results for the NS equation of state.

A. Equilibrium equations in the crust

The CLDM cluster energy (15) is minimized under the constraint of the baryon density n B defined as,

n B = A cl + N g V WS = n cl u + n ng (1 -u) , (39)
=

n ng + 2n e 1 -I cl 1 - n ng n cl . (40)
In Eq. ( 40), the density n B is expressed in terms of four of the five independent variables: A cl , I cl , n cl , n e and n ng . We use the Lagrange multipliers technique, as suggested in Ref. [44], to minimize the µcanonical energy density µcan in the Wigner-Seitz cell,

µcan = tot -µ B n B . (41)
The total energy reads,

tot (A cl , I cl , n cl , n e , n g ) = u cl + (1 -u) g +n e (m p c 2 -m n c 2 ) + n B m n c 2 + ρ e , (42)
where u = V cl /V WS is the volume occupied by a nucleus in a Wigner-Seitz cell, ρ e is the electron energy density (with rest mass), g = MM (n ng = n g , n p = 0) is the neutron gas energy density and cl = MM (n n,cl , n p,cl ) + FS (A cl , I cl , n e , n cl ) is the cluster energy density, with the finite-size contributions given by Eq. ( 18) and discussed in the previous section.

In fact µcan coincides with the pressure, µcan = -P tot . So minimizing µcan is equivalent to maximizing the total pressure P tot . Moreover, minimizing the total energy E tot ≡ tot V WS at fixed baryon density n B is equivalent to minimizing the total Gibbs energy G tot at fixed total pressure, as discussed in Ref. [59].

We define the following thermodynamical quantities (q = n, p),

P cl ≡ n 2 cl ∂E cl /A cl ∂n cl A cl ,I cl , (43)
P g ≡ -g + n g µ g , (44)
µ cl,q ≡ µ nuc,q + P g n B , (45)
µ e ≡ ∂E e ∂N e N cl ,Z cl + 2n e (1 -I cl )A cl ∂E Coul ∂n e , (46)
with

µ nuc,n ≡ ∂E nuc ∂N cl Z cl ,Ne , µ nuc,p ≡ ∂E nuc ∂Z cl N cl ,Ne . (47)
From Eq. ( 43), we deduce,

P cl = P bulk + P Coul + P surf + P curv , (48)
with obvious definitions for these partial contributions. We impose the stationary of the canonical potential (41) with respect to the five independent variables, considering µ B as a constant parameter, and obtain the following equilibrium relations [58],

2E Coul = E surf + 2E curv , (49)
P cl = P g , (50)
µ cl,n = µ g , (51)
µ cl,n = µ cl,p + µ e + ∆mc 2 , (52) µ
B = µ g + 2n e n cl A cl (1 -I cl ) -2n e ∂E surf ∂n g A cl ,I cl ,n cl , (53)
where ∆m = m p -m n . These equilibrium relations have a physical understanding: Eq. ( 49) is an extension of the virial theorem [36] including the curvature energy, Eq. ( 50) reflects the mechanical equilibrium, Eq. ( 51) describes the chemical equilibrium between the cluster and the gas in the r-representation, and finally Eq. ( 52) assures the β-equilibrium. The last Eq. ( 53) describes the baryon chemical potential which fixes the baryon density. Equations ( 49)-( 53) are solved by using the robust Newton-Rhapson method. Note that since the surface energy (27) is independent of the gas density n ng , Eq. ( 53) is subsumed to µ B = µ g .

In the following, we analyse the role played by the different terms in the CLDM and discuss their role in the uncertainties on the various observables in the NS crust.

B. Impact of the FS terms and of the loss function

We first discuss the role of the FS terms and of the loss function, see Eq. (34), used for the confrontation to experimental nuclear masses.

The neutron gas P g and the electron P e pressures are shown in Fig. 7 as function of the baryon density n B . We remind that the equilibrium condition in the WS cell imposes that P g = P cl , and the total pressure is P tot = P g + P e . Fig. 7(a) shows the impact of the FS terms by fixing the FS parameters to their standard values (without optimization), see Table IV, and the bulk term is derived from SLy4 MM . The neutron gas pressure is largely impacted by the FS approximation at low density. The differences between the curves reflects the different density at which the outer-inner crust phase transition takes place. The main differences are between the groups FS1-FS2 and FS3-FS4, reflecting the important role of the curvature term. As the density increases, the effect of FS terms on the cluster pressure is weakened. Note that at low densities the electron pressure dominates over the neutron gas pressure. As a consequence, the total pressure will be weakly impacted by the FS terms, but the equilibrium configurations at the bottom of the inner crust could potentially be. Fig. 7(b) is similar to Fig. 7(a), but we show the electron and gas pressure after the optimization of the FS parameters, i.e. including the impact of the experimental nuclear masses. The dashed (dotted) lines stand for the minimization employing the loss function ∆ E (∆ E/A ). We observe that the optimization of the FS parameters on the experimental nuclear masses considerably reduces the dispersion between the various FS terms, as shown in Fig. 7(a). There are still some differences, but they are smaller once the experimental nuclear masses are considered. For instance, the transition density between the outer and the inner crust is better determined. The impact of the loss function, while not negligible, is also small compared to the uncertainties from the FS terms.

We now analyze more globally the properties of the crust -from outer to inner -in terms of composition (A cl and Z cl ), isospin asymmetries (I cl and Y e ) and volume fraction u, see Fig. 8. Here the bulk term is still SLy4 MM . Note that the right panels (b, d, f, h) shows the same quantities as the left panels (a, c, e, g), but zooming in the crust-core transition region. The end of the curves indicates the boundary of the inner crust where the phasetransition to uniform matter occurs. The crust-core transition is defined as the density at which the energy density in the crust given by Eq. ( 42) matches with e MM the energy density in uniform matter given from the MM at the same density.

As expected, there is a hierarchy in the contributions of the FS terms -from FS1 to FS4 -since the largest difference is observed between FS1 and FS2, then between FS2 to FS3, and finally from FS3 to FS4. The different observable are however not identically sensitive to the FS terms. For instance A cl , Z cl , and I cl are impacted, while Y e , u and n cc (the crust-core transition density) are almost insensitive to the FS terms. For reference, the magic number Z cl =28 and the quasi-magic number Z cl =40 are shown in Fig. 8(c). The number of protons Z cl in the inner crust is more stable employing FS4, which includes curvature and exchange Coulomb terms, as well as a proper treatment of the cluster density n cl in the equilibrium equations, compared to the lower order FS terms. We will see in the following that this stability is also due to the value taken of p surf . The impact of the loss function is also minimal for FS4, compared to the other FS terms. The composition of the neutron star crust is important for the determination of transport properties [62,[START_REF] Schmitt | Reaction Rates and Transport in Neutron Stars[END_REF], which has a direct impact on shear frequencies of neutron stars [57].

Note that even if the effect of the FS terms on the crust-core transition density n cc is small, it can still be discussed from the right panels in Fig. 8. There is a reduction of n cc from FS1 to FS2, consequence of the consistent treatment of the cluster density in the Coulomb and surface terms, which decreases at high density in FS2 while it is fixed to be n sat at all densities in FS1. The cluster energy thus gets higher in FS2, compared to FS1, as the density increases the crust-core transition occurs at a lower density. From FS2 to FS3, the positive curvature term contributing to the energy is compensated in the fit to experimental nuclear masses by a reduction of the surface term for densities given by the experimental data (close to n emp sat ). As a consequence FS2 and FS3 are quite similar. However, the surface and curvature terms have a different density dependence in FS2 and FS3. Then as the global density increases in the crust, implying a decrease of the cluster density, the FS term in FS3 becomes lower than the FS term in FS2, and as a consequence n cc slightly increases. Finally, the exchange Coulomb term is attractive, opposite to the di-rect Coulomb term, but being a small contribution, it very slightly pushes up n cc .

C. Impact of the bulk terms

In this section, we analyse the impact of the bulk terms by varying the Hamiltonians and by analyzing the effect of the effective nucleon mass. We contrast the predictions obtained from the Hamiltonians H1-H4, which represent the best models reproducing the experimental nuclear masses, with the Hamiltonians H5 and H7 for which the reproduction of the experimental data are poorer.

The crust predictions by the chiral Hamiltonians H1-H5 and H7 within the FS4 model are displayed in Fig. 9 similarly to Fig. 8. We notice that the Hamiltonians H5 and H7 -already excluded from the finite nuclei analysis -significantly depart from the predictions by the Hamiltonians H1-H4, which define the light gray band in Fig. 9. They predict lower values for A cl , Z cl , I cl and slightly larger values for Y e and u. Note also that while the nucleon effective mass prescription plays a role, it is much smaller than the uncertainty in the Hamiltonians. The biggest uncertainties are for the composition (A cl , Z cl , I cl ), as seen on top panels. Note that the uncertainty in the value of Z cl originating from the Hamiltonian H1-H4 is small in the outer crust (about 2-4), larger at the bottom of the inner crust (about 8-10), and then becomes very large close to the core-crust transition (about 20 or more). The uncertainties from our best models H2-H3 is however smaller. For instance, in the inner crust, these models predict quite stable values for Z cl up to about 0.01 fm -3 : 34 < Z cl < 40. The vertical gray band in Fig. 9(h) shows the width for n cc corresponding to the two best Hamiltonians H2 and H3: n cc = 8.4 -9.3 10 -2 fm -3 . In Fig. 10 we illustrate the dependence of the EOS on the surface parameter p surf for a reduced set of models. We thus consider our two best models H2 and H3 with FS4 and we fix the loss function to be ∆ E . We vary p surf from 2.5 to 3.5 as suggested in Refs. [38] and [48].

Let us remind that the precise value for p surf cannot be determined from finite nuclei, since its influence on the surface energy becomes important for isospin asymmetries which are way beyond the experimental ones. It was first claimed in Ref. [48] and is similarly illustrated in Fig. 4(b). Fig. 10 shows that p surf has however a remarkable impact in the high density region of the inner crust, close to the core-crust transition. It largely controls for instance the slope of Z cl close to the core-crust transition: a large value of p surf (here 3.5) predicts a decrease of Z cl below 30 at high density while a low value (here 2.5) predicts a increase of Z cl above 40. The stability of Z cl at high density is found only for a very specific value of p surf .

The impact of p surf on the core-crust density n cc is large and comparable with the impact of the Hamiltonian as shown in Fig. 9.

In conclusion, the impact of p surf on the composition of the crust (A cl , Z cl , I cl ) in the high density region is large and quite uncontrolled. More microscopic calculations are required in order to better determine the value of p surf . MM m * models.

V. GLOBAL NS PROPERTIES

The structure of non-rotating neutron stars is provided by the solution of the spherical hydrostatic equations in general relativity, also named the Tolmann-Oppenheimer-Volkof equations [START_REF] Tolman | [END_REF]65],

dm(r) dr = 4πr 2 ρ(r), (54)
dP (r) dr = -ρ(r)c 2 1 + P (r) ρ(r)c 2 dΦ(r) dr , dΦ(r) dr = Gm(r) c 2 r 2 1 + 4πP (r)r 3 m(r)c 2 1 - 2Gm(r) rc 2 -1
, where G is the gravitational constant, c the speed of light, P (r) the total pressure, m(r) the enclosed mass, ρ(r) = ρ B (r) is the total mass-energy density and Φ(r) the gravitational field. P and ρ have contributions from both the baryons (P B , ρ B ) and the leptons (P L , ρ L ). The four variables (m, ρ, P , Φ) are obtained from the solution of the three TOV equations ( 54) and the EOS. In the present calculation, the crust and core EOS are unified, i.e., the same nuclear interaction describes crust and core, as seen in the previous sections. The tidal deformability Λ GW induced by an external quadrupole field is expressed in terms of the Love number k 2 as Λ GW = 2k 2 /(3C 5 ), where the compactness is C = GM/(Rc 2 ), and k 2 is calculated from the pulsation equation at the surface of NS [67,68],

k 2 = 8C 5 5 (1 -2C) 2 (2 -y R + 2C(y R -1)) × 2C(6 -3y R + 3C(5y R -8)) +4C 3 13 -11y R + C(3y R -2) + 2C 2 (1 + y R ) +3(1 -2C) 2 (2 -y R + 2C(y R -1)) ln(1 -2C) -1 , (55)
where y R is the value of the y function at radius R, y R = y(r = R), and y(r) is the solution of the following differential equation,

r dy dr + y 2 + yF (r) + r 2 Q(r) = 0 , (56)
with the boundary condition y(0) = 2 and the functions F (r) and Q(r) defined as,

F (r) = 1 -4πr 2 G[ρ(r) -P (r)]/c 4 1 -2M (r)G/(rc 2 ) , (57)
r 2 Q(r) = 4πr 2 G c 4 5ρ(r) + 9P (r) + ∂ρ(r) ∂P (r) [ρ(r) + P (r)] × 1 -2M (r)G/(rc 2 ) -1 -6 1 -2M (r)G/(rc 2 ) -1 - 4G 2 r 2 c 8 M (r)c 2 + 4πr 3 P (r) 2 1 -2M (r)G/(rc 2 ) -2 . ( 58
)
The NS moment of inertia is obtained from the slow rotation approximation [69,70] as

I = 8π 3 R 0 drr 4 ρ(r) 1 + P ρ(c)c 2 ω Ω e λ-Φ , ( 59
)
where ω is the local spin frequency, which represents the general relativistic correction to the asymptotic angular momentum Ω and λ is defined as exp(-2λ) = 1 -Gm/(rc 2 ). The crust moment of inertia I crust is deduced by considering on the contribution of the crust to the general expression (59).

As usual, for a given EOS the family of solutions is parameterized by the central density or pressure or enthalpy. The EOS are characterized by their evolution in the mass-radius diagram, where both masses and radii of compact stars could in principle be measured, see also Ref. [71].

We show the mass and radius relations for families of NS in Fig. 11. Horizontal error bars show the constrains from NICER analyzis of the pulsar J0030+0451 [1] and J0740+6620 [2], and the constrain from Capano et al. [66] combining multi-messengers signals with nuclear physics. Gray band shows the maximum mass constrain of 2 M ± 0.1 NS [30,31]. Wheat band marks the lowest NS mass observed of 1.17M [72]. Fig. 11(a) shows the influence of the Hamiltonians H1-H4 and the loss function, used in the fit to experimental nuclear masses, on the mass-radius diagram. All curves are calculated with FS4 model and include the effective mass m * . The square-symbols on the curves shows the density above which causality is violated for a given EOS. We note that the loss function has a very small influence for the massradius relation, especially in the domain of observed NS (we observe a very small effect for H3 and H4 for masses below 1.1M ). The Hamiltonian H1 predicts the smallest radius (about 11.5km for 1.4M ), while H3 the highest one (about 12.3km for 1.4M ). Note that this prediction is also influenced by the extrapolation of the chiral Hamiltonian at high density, which we do not discuss in this paper.

In Fig. 11(b) we show the impact of the nucleon mass prescription (bare versus effective mass), considering FS4 finite-size term and ∆ E loss function. The impact of the nucleon mass prescription is smaller than the one of the Hamiltonian, but it is the largest one among all other terms discussed here. Note that the use of the bare nucleon mass makes all EOSs a bit stiffer compared to the use of the effective mass, but the overall impact on the radius is not larger than about 100 m (less than 1% of the radius). The maximum mass is also influenced by the effective mass prescription. We remark that it is increased by 0.15M between the effective mass and the bare mass.

In Fig. 11(c) we show the mass-radius predictions for the two best Hamiltonians, H2 and H3, for which the different FS approximation are shown, fixing the nucleon effective mass and the ∆ E loss function. We note that the finite-size terms has essentially no effect on the global mass-radius relation.

Finally, in Fig. 11(d) we show the effect of varying the surface parameter p surf for H2 and H3, fixing FS4, m * and the ∆ E loss function. Although this parameter has shown to be important for the composition and crust-core transition density, for the macroscopic quantities shown Fig. 11(d) we see no effect of p surf .

While the global mass-radius relation is rather insensible to the details of the energetic modeling of the crust, some more specific properties could be. We show for instance in Fig. 12 the correlation of the crust thickness, CT , with the symmetry energy evaluated at half saturation density for three NS masses. We note that the higher E sym (n = n sat /2) the smaller the crust thickness. By changing the effective mass and the loss function, we also note small effects: the effective mass reduces CT by about 5% compared to the bare mass, while ∆ E reduces by about 1% CT compared to ∆ E/A .

We have performed a more systematical study of the uncertainties related to the crust properties which is shown in Table VI. We present in this table the values obtained by averaging over the different parameters, as well as the standard deviations reflecting their present uncertainties.

We compare low mass NS (1.0M ) with canonical mass NS (1.4M ) for the following quantities: radius, crust thickness, tidal deformability, and crustal moment of inertia. In more details, we have varied the bulk term from H1 to H4, with and without effective mass within the two loss functions, keeping however FS4 fixed. We note a ≈ 3% uncertainty on the radii. Note that this error could be larger if one assumes a non-unified EOS, as shown in Ref. [61]. For the crust thickness we found an error of ≈ 6%. The size of the crust impacts directly the crust moment of inertia, therefore we see an error of the same order for I crust /I. Here we bring attention again for the importance of using a unified EOS, since this error could be as high as 10% [61] in case of a not consistent treatment for the NS crust and core. Regarding the tidal deformability of the 1.4 M NS, the models presented in this work predict Λ 1.4 = 337(±56), therefore all models are inside the constrain from the GW170817 event [3].

VI. CONCLUSIONS

In the present paper we produce a set of compressible liquid-drop models (CLDM) at different orders in the finite-size terms. These CLDM are qualified as unified models for the NS crust and core since the bulk contribution in these two regions are derived from the same model. Here we use the meta-model fitted to a set of Hamiltonians constructed from chiral EFT interactions. Based on these ingredients, we build a set of unified EOS which allows us to explore a number of ap-proximations leading to uncertainties in the crust EOS. We have then investigated the impact of various nuclear approximations and model uncertainties on the prediction for NS crust properties and global quantities. While the present crust model could be extended in the future, we reached some conclusions on the respective impact of various sources of uncertainties which are already interesting for the understanding of NS modeling and of their predictive power.

The confrontation of the CLDM to the experimental nuclear masses presented in this work allowed us to exclude four Hamiltonians: H5, H7, DHS L59 and DHS L69 . We also predicted for the first time the existence of an upper limit for the energy density at saturation max sat , induced by the finite-size terms. Based on the CLDM and chiral EFT Hamiltonians, we obtained the following value for this upper limit: max sat ≈ -2.30 MeV fm -3 . The good models considered here (H2 and H3) lead to an even more refined value for the energy density at saturation which is best sat = -2.70 (20) MeV fm -3 . We have also discussed the effect of including the curvature term in the CLDM: this term leads to a substantial reduction of the residuals between the CLDM and the experimental data, and it also impact the correlation of the surface parameters σ surf,sat and σ surf,sym versus the symmetry energy E sym . In the latter case, we found that the model dependence of this correlation is minimal for our best Hamiltonians H2 and H3. The origin of this effect is still to be understood.

The analysis of NS crust shows that the finite-size terms impact directly the crust composition (A cl , Z cl , and I cl ) and very marginally Y e and the volume fraction u. The impact of the finite-size terms on the core-crust transition density n cc is also small, while non negligible. n cc is however very largely impacted by the Hamiltonian and the surface energy parameter p surf . We finally make predictions on the EOS properties based on the best models we have (H1-H4 for the wider ones and H2-H3 for the smaller ones). Our results somehow establishes new boundaries based on our present knowledge on nuclear physics (χEFT and experimental nuclear masses) and NS crust modeling (represented here by the different finitesize terms employed in the CLDM). These boundaries are important for the modeling of NS phenomenology, such as for instance NS cooling or the understanding of QPO in x-ray flares, oscillation modes in the NS crust as well as spin glitches observed in young NS. In the future, the present approach could be enriched by incorporating the contribution of pairing, neutron skins, different geometries, and of shell effects, and better determination of the parameter p surf shall also be explored.

We remind that several approximations were considered in the present analysis, and it would be interesting to add further terms in the leptodermous expansion to check, for instance, if this expansion continues to provide a pertinent ordering of the missing terms. Among important corrections for future extensions of the present work, let us mention the neutron skin, the in-medium effects on the surface energy, as well as the superfluid contribution from the neutron gas. Thomas-Fermi and Hartree-Fock microscopic calculations also show that the density distribution in the nuclear clusters is not uniform and that shell effects play a role, see for instance the following Refs. [21,22] among others. Finally proton shell effects and neutron pairing could play opposite roles as illustrated in Ref. [43] based on Hartree-Fock-Bogoliubov approach: In the bottom layers of the crust where neutrons are nonexistent or very few, the shell effects induced by protons dominate and impose magicity for the proton number Z, while in the deep layers of the crust where neutrons dominate, superfluidity is maximal for mid-shell nuclei reducing the stabilizing influence of the proton shells.

Finally, the uncertainties on NS macroscopic properties are dominated by the Hamiltonians themselves. A realistic estimation of this error requires however a uni-fied description of nuclear matter in the crust and the core of NS, as implemented in this work.

FIG. 1 .

 1 FIG. 1. Left: Energy per particle as function of the baryon density, in NM and SM. Right: Energy per particle normalized to the free Fermi gas energy (EFFG = (3/5)EF with EF = 2 k 2 F /(2m) and kF the Fermi momentum) in SM and NM function of the neutron Fermi momentum kF n for the six Hamiltonians H1-H7 (except H6). Top panels include DHS Hamiltonians for comparison. Data from the original model are plotted in dots (squares) for H1-H5 and H7 (DHSL59 and DHSL69), with the same color as the MM version. Bottom panels include the nucleon effective mass.

  FIG. 2. Left: Symmetry energy w.r.t the baryon density for H1-H5 and H7, DHSL59, DHSL69 and SLy4. Yellow band shows constrain from neutron skin experiments by PREX-II. Blue bands show constrains from isobaric analog state IAS (IAS + neutron skin, ∆rnp). Right: symmetry energy normalized to the free Fermi gas symmetry energy (Esym,FFG = (3/5)(EF,NM -EF,SM)) w.r.t the Fermi momentum kF . Continuous (dashed) lines consider bare (effective) nucleon mass.

FIG. 3 .

 3 FIG. 3. Pressure as function of the baryon density in βequilibrium for the six Hamiltonians, DHSL59, DHSL69 and SLy4. Continuous (dashed) lines consider bare (effective) nucleon mass. Error bar shows the constrain for the pressure at 2nsat inferred by the LIGO/Virgo collaboration with the GW170817 observation at 90% credible level.

  where the neutron and proton masses m n and m p are fixed to their bare mass, m n c 2 = 939.565346 MeV and m p c 2 = 938.272013 MeV, and n n and n p are the uniform neutron and proton densities in the cluster.

2 / 3 cl

 23 , and the curvature term is in A 1/3 cl . They are expressed as

FIG. 4 .

 4 FIG. 4. Sum of surface and curvature energies normalized by 4πr 2 0,sat A 2/3 cl function of the cluster asymmetry I cl in isolated nuclei. The left panel shows a comparison between the four finite-size models, in the case they are fixed by the standard values given in TableIV(solid lines) and in the case they are fit to the nuclear masses (dotted lines). The right panel shows the influence of the surface parameter p surf at large I cl , after the fit to the nuclear masses.

FIG. 6 .

 6 FIG.6. Top: The isoscalar surface tension parameter σ surf,sat is represented against Esym for the six Hamiltonians H1-H5 and H7, and for the FS models FS1 to FS4. The two Hamiltonians DHS are included with dots (squares) for the minimization to the total energy (energy per particle). Bottom: Same for σsym. Silver band shows the values for the two Hamiltonians which best reproduce nuclear masses, H2 and H3.

FIG. 7 .

 7 FIG. 7. Neutron gas and electron contribution to the pressure on the NS crust. Left panel shows the difference originated by the four finite-size models (with standard parameters, see TableIV). Right panel shows the influence of the different optimization procedures (see Section III C).

FIG. 8 .

 8 FIG. 8. Crust composition (top), asymmetry (center) and volume fraction occupied by the cluster (bottom) for the different ingredients in the CLDM within SLy4MM interaction. FS1 (red), FS2 (magenta), FS3 (green) and FS4 (black) as explained in Tab. V. Dotted lines on panel (c) represents the magic numbers Z cl = 28 and Z cl = 40. Continuous lines represent the parameters fitted to the total energy, while dashed lines represent the fitting to the energy per particle. Left: outer and inner crust. Right: zoom at crust-core transition.

FIG. 9 .

 9 FIG. 9. Crust composition (top), asymmetry (center) and volume fraction occupied by the cluster (bottom) for the different ingredients in the CLDM. In the left panels we show the 6 H from outer to inner crust, on the right we show the 4 H allowed by the finite nuclei analysis with a zoom close the crust-core transition. Dotted lines on panel (c) represent the magic numbers Z cl = 28, Z cl = 40 and Z = 50. Dashed lines include the nucleon effective mass. All figures use FS4 ∆ E model.

FIG. 10 .

 10 FIG. 10. Same as Fig. (9). Comparison of three different values of the surface parameters p surf , within the H2 FS4,∆ E MM m * and H3 FS4,∆ E MM m * models.

FIG. 11 .

 11 FIG. 11. Mass and radius relation. Observational constrains from NICER to the pulsar J0030+0451 and J0740+6620, and the constrain from Capano et al. [66] combining multi-messengers signals with nuclear physics are shown in error bars. Wheat band shows the lowest mass NS observed with a mass of 1.17M . Marks on curves represent the density where causality is broken for a given EOS. Top panels show the results for H1-H4 and (a) Influence of the optimization procedure (b) effect of the inclusion of the nucleon effective mass. Bottom panels show the results for the two Hamiltonians the better reproduce nuclear masses, H2 and H3, and (c) the impact of finite-size terms, (d) the impact of the surface parameter p surf . See text for more details.

FIG. 12 .

 12 FIG. 12. Correlation of crust thickness for a 1.0 M , 1.2 M and 1.4 M NS and the symmetry energy evaluated at nsat/2. Dashed lines include the nucleon effective mass.

TABLE IV .

 IV Standard FS parameters for the CLDM considered in this work. Note the associated value b surf = 29.9.

	σ surf,sat	σ surf,sym	p surf	σcurv,sat	βcurv
	MeV fm -2	MeV fm -2		MeV fm -1
	1.1	2.3	3.0	0.1	0.7

5/3 cl

and therefore dominates the CLDM energy at large A cl . Since the Coulomb term is repulsive, this induces a limitation in the maximum A cl for finite nuclei. However, for

TABLE V .

 V Definition of nuclear macroscopic models used in this work. In the left column we have the different terms on the CLDM implemented in this work. On the first line we have the label of the model. The table show which term and how many parameters are necessary to fix each model.

	Model	Variables FS1 FS2 FS3 FS4
	Bulk from MM	(I cl , n cl )	×	×	×	×
	FS Surface	(nsat)	×	-	-	-
	FS Coulomb (Dir.)	(nsat)	×	-	-	-
	FS Surface	(n cl )	-	×	×	×
	FS Coulomb (Dir.)	(n cl )	-	×	×	×
	FS Curvature	(n cl )	-	-	×	×
	FS Coulomb (Ex.)	(n cl )	-	-	-	×
	Number of param.		3	3	5	5

ACKNOWLEDGMENTS G.G., J.M. and R.S. are supported by CNRS grant PICS-08294 VIPER (Nuclear Physics for Violent Phenomena in the Universe), the CNRS IEA-303083 BEOS project, the CNRS/IN2P3 NewMAC project, and benefit from PHAROS COST Action MP16214. S.R. is supported by Grant No. DE-FG02-00ER41132 from the Department of Energy , and the Grant No. PHY-1430152 (JINA Center for the Evolution of the Elements), and PHY-1630782 (Network for Neutrinos, Nuclear Astrophysics, and Symmetries (N3AS)) from the National Science Foundation. This work is supported by the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) of the Université de Lyon for its financial support within the program Investissements d'Avenir (ANR-11-IDEX-0007) of the French government operated by the National Research Agency (ANR).