Relativistic mean-field study of the neutron star inner crust using the asymmetric finite difference method

Neutron stars represent the densest observable matter in the Universe and serve as unique natural laboratories for testing nuclear physics under extreme conditions [1, 2]. A typical cold neutron star is structured into three distinct regions beneath its atmosphere: the inhomogeneous outer crust, the inner crust, and a dense, homogeneous liquid core. The outer crust consists of a body-centered cubic (BCC) lattice of nuclei embedded in a degenerate electron gas, which ensures global charge neutrality. With increasing depth, the baryon number density rises accordingly, and the nuclei become progressively more neutron rich [3−5]. When the density reaches the neutron drip threshold, nuclei are no longer able to bind all neutrons [6, 7]. This leads to the neutron drip phenomenon, which marks the transition to the inner crust. In this region, nuclear matter comprises extremely neutron-rich clusters immersed in a dilute gas of free neutrons and electrons. As the density further increases toward the crust–core transition, the interplay between nuclear surface tension and Coulomb repulsion favors the formation of complex non-spherical structures, commonly referred to as nuclear pasta phases [8].

The matter in the inner crust of a neutron star can be regarded as a nearly perfect crystal and should, in principle, be described by the band theory of solids. Such a treatment naturally incorporates periodic boundary conditions and the long-range correlations of the neutron gas. Band-structure effects in the inner crust have been discussed in many works based on nonrelativistic density functional theory, in which pairing correlations of nucleons were also taken into account [9−12]. For practical calculations, however, the Wigner–Seitz (WS) cell approximation is commonly adopted to simplify the underlying periodic structure [13]. The nucleon–nucleon interactions for neutron-rich nuclei in the inner crust are usually treated within nuclear density functional theory, such as the Skyrme–Hartree–Fock model [14], the Gogny Hartree–Fock model [15], and the relativistic mean-field (RMF) theory [16−18]. All these approaches have successfully described the properties of nuclear matter and finite nuclei.

In general, the nuclei appearing in the inner crust contain a very large number of neutrons, so that many occupy continuum states. Consequently, the choice of boundary conditions for these unbound states is essential for solving the corresponding Schrödinger or Dirac equations. In the nonrelativistic framework, the first derivative of odd-parity wave functions vanishes at the cell radius, whereas even-parity wave functions vanish at the cell boundary [13]. However, for the Dirac equation, the boundary conditions at the cell radius are not uniquely defined and involve a certain degree of arbitrariness. As demonstrated by Cao et al. in Ref. [19], a direct extension of nonrelativistic boundary conditions to the relativistic regime violates the orthogonality of single-particle states. The semiclassical Thomas–Fermi (TF) approximation has also been widely applied to calculate the equation of state (EOS) of the inner crust within RMF models [20−22], in which shell effects are completely neglected. This approach further enables computationally efficient, large-scale investigations of neutron-star matter.

Traditional techniques, such as the shooting method [23] and basis-expansion approaches [24−27], are well established as robust tools for standard nuclear-structure problems. However, their application to the neutron-star inner crust faces significant challenges. As noted in Ref. [28], conventional methods can be highly sensitive to the choice of box size or require an extremely large basis space to accurately describe the weakly bound and continuum states characteristic of the inner crust.

Alternatively, discretizing the Hamiltonian on a coordinate grid—the finite-difference method (FDM) [29−31]—provides an efficient and accurate approach. Because this method does not rely on basis expansions, it is particularly well suited for describing states with extended spatial distributions. Despite these advantages, the direct application of FDM to the Dirac equation suffers from the fermion-doubling problem, which is also well known in lattice quantum chromodynamics (QCD) [32, 33]. In particular, the standard central-difference discretization generates unphysical, rapidly oscillating spurious states that contaminate the physical spectrum. Although remedies such as the Wilson mass term are commonly employed in lattice QCD to suppress these modes, they inevitably introduce artificial modifications to the Hamiltonian.

In this work, the asymmetric finite-difference (AFD) method is adopted to overcome these numerical difficulties and to solve the radial Dirac equations for neutron-rich nuclei in WS cells. As proposed in Ref. [28], this method employs distinct finite-difference formulas for the large and small components of the Dirac spinor, depending on the sign of the quantum number κ. This prescription effectively eliminates spurious states for massive fermions without introducing artificial terms such as the Wilson mass, while preserving the hermiticity of the Hamiltonian matrix. As a result, the AFD method enables an accurate and fully quantum-mechanical description of extremely neutron-rich nuclei in the inner crust of neutron stars.

The paper is organized as follows. In Sec. II, we briefly outline the RMF framework and introduce the AFD method for solving the Dirac equations in the WS cell. Section III presents the results and discussion, focusing on the binding energies per nucleon, density distributions, and the neutron root-mean-square radii and neutron chemical potentials of inner-crust nuclei, using TM1-family parameter sets. A summary is given in Sec. IV.

In this work, the properties of nuclear systems are investigated within the RMF formalism. In this framework, the interaction among nucleons is described by the exchange of mesons, including the isoscalar-scalar meson (σ), the isoscalar-vector meson ($ \omega^\mu $), and the isovector-vector meson ($ \rho^{a\mu} $). In addition, the electromagnetic interaction between protons is accounted for through the photon field ($ A^\mu $). A cross-coupling term between the $ \omega^\mu $ and $ \rho^{a\mu} $ mesons is included to adjust the density dependence of the nuclear symmetry energy and its slope [34, 35]. The Lagrangian density of the RMF model is given by:

Here, $ W^{\mu\nu} $, $ R^{a\mu\nu} $, and $ F^{\mu\nu} $ denote the antisymmetric field tensors associated with the $ \omega^\mu $, $ \rho^{a\mu} $, and $ A^\mu $ fields, respectively. Within the mean-field approximation, the meson fields are treated as classical, and the corresponding field operators are replaced by their expectation values. Consequently, the only nonvanishing classical fields are $ \sigma = \langle\sigma\rangle $, $ \omega = \langle\omega^0\rangle $, and $ \rho = \langle\rho^{30}\rangle $. From the Lagrangian density given in Eq. (1), the equations of motion for the meson fields follow.

The source terms appearing on the right-hand sides of the meson field equations are constructed from the nucleon wave functions, with $ \rho_i^s $ denoting the scalar density and $ \rho_i $ the vector density for the species $ i=n,\; p $. To evaluate these densities, we first solve the Dirac equation for the nucleons. For a spherically symmetric nucleus, the single-particle nucleon wave function $ \phi_\alpha $ can be expressed as a two-component Dirac spinor

In this expression, $ G_{n\kappa}(r) $ and $ F_{n\kappa}(r) $ denote the radial wave functions of the large and small components, respectively. $ \Phi_{\kappa m}(\hat{{\bf{r}}}) $ denotes the spinor spherical harmonics. The single-particle state α is characterized by the set of quantum numbers $ {n, \kappa, m} $, where n is the radial quantum number, κ is the relativistic angular-momentum quantum number defined as $ \kappa = (-1)^{j+l+1/2} (j+1/2) $, and m is the projection of the total angular momentum j on the z axis. The wave function is normalized to unity, and its radial dependence is governed by the following set of coupled equations

where $ M^* = M + g_\sigma \sigma $ denotes the effective nucleon mass. Once the nucleon wave functions have been determined, the source terms that enter the meson field equations can then be evaluated. The scalar and vector densities are then constructed from the wave-function components as follows:

Here, the summation runs over all occupied states, and $ w_\alpha $ denotes the corresponding occupation probability. In the inner crust of neutron stars, dripped neutrons are expected to form spin-singlet Cooper pairs through the attractive interaction in the $ ^1S_0 $ channel. These pairing correlations can significantly affect a range of neutron-star phenomena, including the crust's microscopic properties, its thermal behavior, and its response to external perturbations. However, in the present work we do not include pairing correlations. Consequently, the occupation probability $ w_{\alpha} $ is determined using the uniform filling approximation.

The meson-field equations and the Dirac equations for nucleons are mutually coupled: the mean-field potentials determine the nucleon wave functions, which in turn generate the source densities for the meson fields. Accordingly, this coupled system is solved self-consistently by iteration, which is continued until both the potentials and the densities converge to within a prescribed numerical accuracy.

To solve the Dirac equation using the AFD method, we perform the calculation within a finite box of radius $ R_c $. This domain is discretized into N equidistant grid points, $ r_i = i \times h $ (for $ i=1, ..., N $), with a step size $ h = R_c / N $, so that the radial wave functions $ G(r) $ and $ F(r) $ can be represented as N-dimensional vectors for a fixed κ. In the AFD method, the first-order derivative $ \mathrm{d}/\mathrm{d}r $ is approximated by asymmetric finite-difference formulas. In this work, we employ the five-point approximation. For a general function $ f(r) $, the forward difference formula at the grid point $ r_i $ is given by

The corresponding backward difference formula is given by:

This discretization transforms the derivative operator into a sparse matrix. For example, the backward-difference formula is represented by a lower-banded matrix.

The forward-difference formula likewise yields an upper-banded matrix structure. Applying the discretization scheme described above to the radial Dirac equations (7) and (8) converts the original differential eigenvalue problem into a $ 2N \times 2N $ matrix eigenvalue problem. This problem can be written in the following matrix form:

In this formulation, the diagonal blocks A and C contain the potential and mass terms, while the off-diagonal blocks $ {\bf{B_1}} $ and $ {\bf{B_2}} $ include the derivative and centrifugal terms proportional to $ \kappa/r $. The boundary condition imposes that the wave functions vanish at $ r = 0 $ and outside the box, $ r \gt R_c $ [28]. To eliminate spurious states while preserving the hermiticity of the Hamiltonian matrix, the choice of the finite-difference scheme is crucial. The finite-difference formulas are selected according to the sign of κ. For $ \kappa \gt 0 $, the large component G is discretized using a backward-difference form, while the small component F is treated with a forward-difference form; the assignment is reversed for $ \kappa \lt 0 $. This distinction necessitates treating states with different signs of κ separately in numerical implementations. The present approach follows the methodology of Ref. [28]. A similar consideration applies to the boundary conditions, where Cao et al. [19] demonstrated that strict preservation of the orthogonality of single-particle states can be achieved using a scheme that depends on the sign of κ, specifically by imposing a vanishing boundary condition on the appropriate component.

In this work, the properties of extremely neutron-rich nuclear systems are investigated within the WS cell approximation. It should be noted that the present setup does not aim to provide a fully realistic description of the neutron-star inner crust, since the electron gas and β-equilibrium conditions are not explicitly taken into account. To systematically explore the influence of the density dependence of the symmetry energy on inner-crust structures, a family of parameter sets based on the TM1 interaction is adopted [34], with symmetry-energy slope L in the range $ 40 $–$ 110 $ MeV. The symmetry-energy slope L governs the density dependence of the symmetry energy and is strongly correlated with the neutron-skin thickness of finite nuclei and the radii of neutron stars. However, current experimental and observational constraints still involve large uncertainties. The PREX-II experiment favors a relatively large value, exceeding $ 100 $ MeV, whereas astrophysical observations of neutron stars and gravitational waves suggest a value around $ 40 $–$ 80 $ MeV. Therefore, in the present work, we choose L in the range $ 40 $–$ 110 $ MeV [36].

In this scheme, the couplings of the isoscalar-scalar and isoscalar-vector mesons are kept fixed at their original TM1 values, thereby preserving the saturation properties of symmetric nuclear matter. The isovector meson couplings, $ g_\rho $ and $ \Lambda_{\mathrm{v}} $, are adjusted to generate different values of the symmetry-energy slope L at nuclear saturation density, while the symmetry energy is constrained to be the same at the subsaturation reference density of $ n_{\text{fix}} = 0.11 \text{ fm}^{-3} $. Specifically, calculations are performed using parameter sets with $ L = 40 $, $ 60 $, and $ 80 $ MeV, together with the original TM1 parameter set corresponding to $ L = 110.8 $ MeV.

In addition, to examine the impact of the nucleon effective mass, the recently developed TM1m* parameter set is also employed [37], since an EOS with a larger effective nucleon mass makes supernova explosions more likely to occur. This interaction features a larger effective mass ($ M^*/M \approx 0.79 $) than the standard TM1 value ($ M^*/M \approx 0.63 $) and has been adjusted to better reproduce the ground-state properties of finite nuclei. The TM1m* interaction shares the same symmetry-energy slope ($ L = 40 $ MeV) as the corresponding model in the TM1-based series discussed above. The detailed coupling constants of all parameter sets used in this work are listed in Table 1, and the remaining parameters are identical to those of the original TM1 set [38].

In this work, we adopt the nine WS cell configurations with different mass numbers A for Zr and Sn listed in Ref. [19]. These configurations correspond to baryon densities ranging from $ \rho = 0.0016 \rho_0 $ to $ 0.12 \rho_0 $, where $ \rho_0 = 0.17 \text{ fm}^{-3} $. In the original many-body study of the neutron-star inner crust [13], it was shown that, at a given baryon density, local minima of the nuclear binding energy arise from shell closures at proton numbers 40 and 50 when pairing correlations are neglected. The ground-state solutions obtained in that work, for densities from $ 4.6\times10^{11} $ g/cm$ ^3 $ to $ 3.4\times10^{13} $ g/cm$ ^3 $, correspond to nuclei ranging from $ ^{180} $Zr to $ ^{1800} {\rm{Sn}}$, as listed in Table 2. Subsequent studies investigated these nuclei to elucidate the structure of the neutron-star inner crust; see, e.g., Refs. [19, 39, 40]. Therefore, we also adopt these isotopes here for comparison. The detailed parameters for each layer, including the proton number Z, the neutron number N, and the cell radius $ R_c $, are summarized in Table 2.

Based on the WS cell configurations detailed in Table 2, we calculate the binding energy per nucleon $ E/A $ for the neutron-rich nuclei in the inner crust. These calculations are performed to assess the energetics of the system under extreme neutron-rich conditions. For comparison, we also perform calculations within the TF approximation, following the self-consistent relativistic TF approach; details are given in Ref. [37]. In this approximation, shell effects are neglected, and nucleons are treated as locally uniform Fermi gases. The local scalar and vector densities are determined from the local Fermi momenta, while the neutron and proton chemical potentials are spatially constant in the WS cell. For a fixed cell radius $ R_c $ and given proton and neutron numbers, we start from trial mean fields (or, equivalently, trial density distributions), determine the corresponding chemical potentials, calculate the local densities, and solve the meson field equations, Eqs. (2-5), self-consistently. The iteration is repeated until the densities, chemical potentials, and mean fields converge. Thus, the TF results presented in this work correspond to the same prescribed WS cells as those adopted in the AFD calculations and serve as semiclassical reference results for comparison. The values of $ E/A $ for the nine inner crust layers, employing the TM1 family parameter sets with different symmetry energy slopes L and effective masses, are summarized in Table 3.

When examining the TM1 family of parameter sets, we observe a systematic trend in which the binding energy per nucleon decreases with increasing symmetry-energy slope L. This suggests that higher values of L correspond to a lower total energy for the specified inner-crust configurations, since the symmetry-energy slope plays a negligible role below the saturation density $ \rho_0 $ compared with its role above $ \rho_0 $. A smaller L yields a larger symmetry energy in the density range relevant to finite nuclei. This means that the interaction in neutron-rich nuclei becomes more repulsive, thereby reducing the binding energy.

It should be emphasized that the present discussion of the L dependence of $ E/A $ is conducted for the fixed WS cell configurations adopted in Table 2. In realistic inner-crust calculations, however, one should impose the β-equilibrium condition and determine the energetically optimal configuration at each baryon density. In that case, the optimal proton fraction, as well as the corresponding neutron and proton numbers, would generally depend on L. Therefore, the influence of L in the inner crust is reflected not only in the energy variation of a given configuration, but also in possible changes of the equilibrium composition. When comparing two interactions with the same slope, $ L=40 $ MeV, the TM1m* set consistently yields lower $ E/A $ values than the corresponding TM1-based set. This discrepancy is particularly significant in layers with higher nucleon densities, indicating that the larger effective mass in the TM1m* set plays a vital role in reducing the system's energy.

It should be noted that the lower $ E/A $ obtained with the TM1m* interaction should not be attributed solely to the kinetic-energy reduction associated with a larger effective mass. Since TM1m* also differs from the TM1-based $ L = 40 $ set in its self-consistent mean-field parameters, the reduction of the total energy should be understood as a combined effect of the larger effective mass and the corresponding rearrangement of the scalar and vector mean fields, which modifies both the kinetic- and potential-energy contributions.

Regarding the numerical methods employed, results obtained from the TF approximation are generally lower than those derived from the AFD method. For the TM1-based parameter sets (excluding TM1m*), the results from TF and AFD are quite similar in the lower-density region (Layers 1–5), where $ E/A $ is less than zero. However, the differences become more pronounced for the TM1m* set, as shown in Ref. [37], as well as in higher-density layers. This variance is typically attributed to the inclusion of quantum shell effects and gradient terms in the AFD method, which are often averaged out or neglected in the semiclassical approximation.

In Fig. 1, we show the proton and neutron density distributions of $ ^{950} {\rm{Sn}}$, an extremely neutron-rich nucleus, calculated with the original TM1 interaction ($ L=110.8 $ MeV) using both the AFD method and the TF approximation. The solid lines depict the AFD results, while the dashed lines represent those from the TF approximation. The nucleon densities obtained with the AFD method clearly exhibit pronounced fluctuations, particularly in the central region of the nucleus ($ r \lt 7 $ fm). These oscillations arise from quantum shell effects, reflecting the discrete nature of the single-particle wave functions [37]. In contrast, the TF approximation yields smooth and nearly flat density distributions. At larger distances ($ r \gt 10 $ fm), a finite, approximately constant neutron density is observed, indicating the presence of a surrounding neutron gas in the inner crust environment. Meanwhile, the proton density rapidly drops to zero around $ r=9 $ fm owing to the small Z.

Figure 1.  Proton and neutron density distributions for $ ^{950} {\rm{Sn}}$ calculated with the TM1 interaction ($ L=110.8 $ MeV). The solid lines show the AFD results, while the dashed lines indicate the TF approximation.

In Fig. 2, the nucleon density distributions for the nine layers of the inner crust are presented, calculated using different parameter sets. Initially focusing on the neutron distributions within the TM1-based family, which retains identical isoscalar properties, we observe a systematic dependence of the distributions on the symmetry-energy slope L. In the central region of the WS cell, the neutron density exhibits an inverse correlation with L: the parameter set with $ L=40 $ MeV shows the highest central density, whereas the original TM1 set ($ L=110.8 $ MeV) yields the lowest. However, at the nuclear surface this trend reverses, with larger L generally corresponding to higher neutron densities. As discussed in Ref. [34], this redistribution occurs because a larger L yields a larger symmetry energy in the dense central region, which favors reducing the central neutron density and enhancing the distribution at the surface. In the outer tail region, however, the parameter sets with smaller L eventually sustain a higher neutron density than those with larger L. While this effect is subtle and almost imperceptible in the lower-density layers [Figs. 2(a)–2(e)], it becomes distinctly evident in the boundary drop-off region of Fig. 2(i). Regarding protons, their distributions are relatively less sensitive to variations in L. Minor differences among the curves are still apparent, since protons and neutrons are coupled in the self-consistent mean field through the isovector meson. For smaller L, the proton densities are suppressed as the neutron number increases.

Figure 2.  (color online) Proton and neutron density distributions in WS cells for nine inner-crust layers, with panels (a)-(i) corresponding to Layers 1-9. Results are obtained using the AFD method with different parameter sets, namely $ L=40, 60, 80,110.8 $ MeV from the TM1 family and the TM1m* interaction.

In examining the impact of the effective nucleon mass, we compare the TM1m* set with the TM1-based $ L=40 $ MeV set, since both have the same symmetry-energy slope. The first three panels, Figs. 2(a)-(c) for Layers 1–3, reveal a crossover behavior in the central region: the TM1-based $ L=40 $ MeV set shows a higher central peak, but its density decreases more rapidly with increasing radius. As a result, the TM1m* density eventually exceeds that of the $ L=40 $ MeV set in the outer part of the central region. In Figs. 2(d)-(e) for Layers 4 and 5, the density profiles of both parameter sets are quite similar in the central area.

However, a distinct trend emerges in Figs. 2(f)-(i) for Layers 6–9 for Sn isotopes, which correspond to extremely neutron-rich conditions. In these layers, the TM1m* profile consistently shows a lower central density than the TM1-based $ L=40 $ set throughout the central region, and its curve generally lies between those of $ L=40 $ and $ L=60 $. This suggests that, for a fixed symmetry-energy slope, increasing the effective nucleon mass tends to diminish the nucleon density in the central region, producing an effect comparable to that of increasing L while keeping the effective mass fixed. A larger effective mass also yields a weaker spin-orbit coupling.

Furthermore, the neutron density distributions for $ ^{1800} {\rm{Sn}}$ exhibit sharp decreases near the boundary, owing to the boundary condition mentioned in Sec. II B. This decrease is independent of the choice of boundary conditions, as a similar feature is also observed in Fig. 1 of Ref. [19].

The calculated neutron root-mean-square radii, $ R_n $, for the nine layers of the inner crust are presented in Table 4. Within the TM1 family, $ R_n $ decreases with increasing L. This trend is most pronounced in the lower-density layers (e.g., Layers 1–6), underscoring the influence of the symmetry energy on the spatial distribution of neutrons. This indicates that, despite surface-region expansion at larger L, the higher neutron density in the center region for the low-L sets (as discussed in Fig. 2) dominates the total radius.

Moreover, the contrast in $ R_n $ between large and small L is more pronounced with the AFD method than with the TF approximation, because in AFD the neutron gas density decays to zero near the boundary. Notably, the magnitude of this discrepancy depends strongly on the layer's average density. The difference between the TF and AFD results is largest in the low-density layers (e.g., Layer 1) and gradually diminishes with increasing layer number. By Layer 9, where the density is highest, the results from the two methods are much closer, suggesting that the quantum shell effects captured by the AFD method are more pronounced in dilute environments.

To assess the impact of the effective nucleon mass, we compare the TM1m* parameter set with the TM1-based $ L=40 $ configuration. The TM1m* interaction yields smaller neutron radii than the TM1-based $ L=40 $ set, suggesting that a larger effective mass leads to a more concentrated neutron distribution within the WS cell.

Table 5 presents the neutron chemical potentials $ \mu_n $ for the nine inner-crust layers, calculated with different parameter sets, which are crucial for the β-equilibrium conditions. It is noteworthy that $ \mu_n $ decreases significantly with increasing symmetry-energy slope L. For instance, in Layer 9, the chemical potentials calculated with the $ L=40 $ MeV set are more than twice those obtained with the original TM1 set, in both the TF approximation and the AFD method. This implies that more high-angular-momentum orbitals are occupied at smaller L. Regarding the effective nucleon mass, the TM1m* interaction yields systematically lower chemical potentials than the $ L=40 $ MeV set, despite both having the same symmetry-energy slope. This reduction is attributed to the larger effective mass in TM1m*, which suppresses the kinetic-energy contribution to the Fermi energy.

When comparing the two numerical methods, the neutron chemical potentials obtained with the AFD method are generally larger than those from the TF approximation, indicating that quantum shell effects typically enhance the chemical potential in these inner-crust environments. However, a slight reversal is observed for TM1m* in the dilute Layers 1 and 2, where the TF values marginally exceed the AFD ones. This anomalous behaviour aligns with Ref. [37], which reported that parameter sets with large effective masses (such as TM1m*) can exhibit an opposite trend in the deviations between semiclassical and quantum calculations, compared with other models like the original TM1 and the TM1-based $ L=40 $. As suggested in Ref. [37], this difference is largely attributable to specific nuclear-matter properties, particularly the effective mass. Interestingly, this distinctive feature appears to be confined to the low-density regime, as the TM1m* results return to the general trend in the denser layers.

In this work, we systematically explore the ground-state properties of extremely neutron-rich nuclear systems in the inner crust of neutron stars, employing the RMF formalism within the Wigner–Seitz approximation. By adopting a series of TM1-family parameter sets distinguished by varying symmetry-energy slopes L and effective nucleon masses, we analyze the effects of the symmetry-energy slope L and the effective nucleon mass $ M^* $ on inner-crust properties. To address the spurious-state problem in the numerical procedure, we implement the AFD method to solve the radial Dirac equations and compare the results with those derived from the TF approximation.

The investigation demonstrates that, in contrast to the smooth density profiles yielded by the TF approximation, the AFD method captures significant density fluctuations in the central region and leads to noticeable differences in the neutron chemical potentials. These shell effects may strongly influence the microscopic structure, equilibrium composition, superfluid-related properties, and transport properties of representative Wigner–Seitz cells in the inner crust, particularly in the low-density region. This indicates that quantum shell effects are non-negligible for the microscopic properties of the neutron-star inner crust.

Building on these results, the study highlights the crucial influence of nuclear equation-of-state parameters on the inner-crust structure. We find that the symmetry-energy slope L and the effective nucleon mass play similarly critical roles in regulating the system’s properties. Specifically, an increase in either the slope L or the effective nucleon mass leads to stronger neutron binding, resulting in a more compact neutron distribution within the cell and a significant increase in the binding energy. These findings elucidate that under extremely neutron-rich conditions, the density dependence of the symmetry energy and the effective nucleon mass constitute the core physical mechanisms determining the microscopic structure and energetic characteristics of the inner crust.

The β-equilibrium conditions will be considered in future work to construct the inner-crust EOS. Furthermore, periodic boundary conditions in the Dirac equation will be treated using the AFD method to investigate band-structure effects in the inner crust. We also note that neutrons in the inner crust are expected to exhibit pairing correlations, particularly in the weakly bound and continuum regions, which may affect the detailed shell structure and related microscopic properties. A self-consistent inclusion of pairing would require an extension to an RMF+BCS or relativistic Hartree–Bogoliubov treatment in coordinate space. Such an extension is compatible with the present framework and will be considered in future work.