NEUTRON STARS AS DENSE LIQUID DROP AT EQUILIBRIUM WITHIN THE EFFECTIVE SURFACE APPROXIMATION

R.C. Tolman suggested [1] first to study the simplest model for a neutron star (NS) considering it as a dense liquid-matter drop at its equilibrium under the gravitational, nuclear and other realistic forces; see also his book [2], chapt. 7, sect. 96. Within this model, the Einstein-Gilbert equations of the General Relativistic Theory (GRT) for the spherical symmetry case has been reduced to the three independent equations for four unknown quantities, namely, two parameters, $\lambda$ and $\nu$, of the Schwarzschild metric (see also Ref. 3, the chapters 11 and 12)

and the pressure $P$, and energy density $\mathcal{E}$. To get the complete system of equations Tolman suggested also to derive independently the equation of state (EoS), $\mathcal{E}=\mathcal{E}(\rho)$, or $P=P(\rho)$. The EoS can be found from the condition of a static equilibrium for a liquid-matter drop under gravitational, and nuclear (and Coulomb) forces, which all should be taken into account by the simplified GRT equations for the gravitational field. Following Tolman’s ideas, Oppenheimer and Volkoff[4] derived the simple (TOV) equations by using essentially the macroscopic properties of the system [Refs. 3 (chapt. 12) and 5 (chapt. 1)] of the system as a perfect liquid drop at equilibrium. As shown in Ref. 2, one can derive the analytical solution for the pressure $P$ as function of the radial coordinate. So far, the TOV equations [4] are considered with the equation of state, $\mathcal{E}=\mathcal{E}(\rho)$, which was obtained independently of the macroscopic assumptions used in the TOV derivations. For instance, the EoS with a polytropic expression for the pressure[6, 7] $P=P(\rho)$ as function of the density $\rho$ is assumed to be similar to that for a particle gas system with fitted parameters. The gradient terms are usually neglected in the energy density, and then, no equilibrium for a finite gas system with no external fields. However, the pressure $P$ was calculated as function of the radial $r$ coordinate by using the TOV equations for another dense liquid-drop system. Such a system, including the gravitation, can be described mainly at a stable equilibrium by short-range forces, in contrast to a gas matter with long-range inter-particle interaction. Then, the NS mass $M$ was obtained numerically as function of the NS radius $R$, $M=M(R)$, with important restrictions due to those of the NS radius. In any case, we should emphasize the importance of the gradient terms in the energy density $\mathcal{E}(\rho)$ for a macroscopic condition of the system equilibrium, in spite of a relative small NS crust thickness. As a hint to these conclusions, see many applications of the TOV equations, e.g., in Refs. 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18. Concerning the relation of the nuclear and neutron liquid-drop models (LDMs) to NS properties, we should also mention the work[19] by Baym, Bethe, and Pethick. They discuss the liquid matter drop with the leptodermic property having a sharp decrease of the particle number density in a relatively small edge region considered as its surface. Clear specific definitions and complete updated results for the energy density with density gradient (surface) terms and for equations of the infinite matter state with many inter-particle forces in the non-relativistic and relativistic cases can be found in the recent review Ref. 18.

Taking Tolman’s ideas, one can try first to extend the EoS to those for a dense macroscopic[3, 5] system of particles. In this leptodermic system, the particle density $\rho$ is function of the radial coordinate with exponentially decreasing behavior from an almost constant saturation value inside of the dense system to that through the effective surface (ES) of NS in a relatively small crust range $a$; see Fig. 1. The ES is defined as the points of spatial coordinates with a maximum density gradient. To obtain the analytical solutions for the particle density and EoS, we use the leptodermic approximation[3, 2, 4, 18, 21, 22, 23, 24, 25, 26, 27, 19, 20, 17, 6], $a/R\ll 1$, where $R$ is the (curvature) radius of the ES.

Within this effective surface approximation (ESA), $a/R\ll 1$, simple and accurate solutions of many nuclear and liquid-drop problems involving the particle number density distributions were obtained for nuclei[28, 29, 30, 31, 32, 33, 34, 35]. The ESA exploits the property of saturation of the nuclear particle density $\rho$ inside of the system, which is a characteristic feature of dense systems as molecular systems¹¹1 For the one-dimensional and more complicate dense molecular (e.g., liquid-drop) systems, van der Waals (vdW) suggested the phenomenological capillary theory[36] which predicted the results for the particle number density $\rho$ and surface tension coefficients $\sigma$. These results are similar to those obtained later in Refs. 28, 29; see also comments below., liquid drops, nuclei, and presumably, NSs. The realistic energy-density distribution is minimal at a certain saturation density of particles (nucleons, neutrons, or nuclei) corresponding approximately to the infinite matter[37]. As a result, relatively narrow edge region exists in finite nuclei or NS (crust) in which the density drops sharply from its almost central value to zero. We assume here that the part inside of the system far from the ES can be changed a little (saturation property of the dense system as hydrostatic liquid drop, nucleus or NS in the final evolution state).

The coordinate system related to the effective surface is defined in such a way that one of the spatial coordinates ($\xi$) is the distance from the given point to the ES; see Fig. 1, and A. This non-linear coordinate system is conveniently used in the region of nuclear and NS edges. They allow for an easy extraction of relatively large terms in the density distribution equations for the variation equilibrium condition. This condition means that the variation of the total energy $E$ over the density $\rho$ is zero under the constraints which fix some integrals of motion beyond the energy $E$ by the Lagrange method. The Lagrange multipliers are determined by these constraints within the local energy-density theory, in particular, the extended Thomas-Fermi (ETF) approach from nuclear physics[39, 38]. Neglecting the other smaller-order perpendicular- and all parallel-to-ES contributions, sum of such terms leads to a simple one-dimensional equation (in special local coordinates with the coordinate normal-to-surface $\xi$); see Fig. 1.

Such an equation mainly determines approximately the density distribution across the diffused surface layer of the relatively small order of the ratio of the diffuseness parameter $a$ to the (mean curvature) ES radius $R$, $a/R\ll 1$, for sufficiently heavy systems. Notice that within this manuscript, as in Refs. 28, 29, 30, 31, 32, 33, 34, 35, the “diffuseness parameter”, “ the crust range”, and “the thickness of the system edge” have the same meaning for neutron stars. A small parameter, $a/R$, of the expansion within the ESA can be used for analytical solving the variational problem for the minimum of the system energy with constraints for a fixed particle number, and other integrals of motion, such as angular momentum, quadrupole deformation, etc. When this edge distribution of the density is known, the leading static and dynamic density distributions which correspond to diffused surface conditions can be easily constructed. To do so, one has to determine the dynamics of the effective surface which is coupled to the volume dynamics of the density by certain LDM boundary conditions[40, 31, 41]. This ESA approach is based on the catastrophe theory for solving differential equations with a small parameter of the order of $a/R$ as the coefficient in front of the high order derivatives in normal to the ES direction[42]. A relatively large change of the density $\rho$ on a small distance $a$ with respect to the curvature radius $R$ takes place for the liquid-matter drop (nuclei, water drops, neutron stars). Inside of such dense systems, the density $\rho$ is changed slightly around a constant saturation density $\overline{\rho}$. Therefore, one obtains essential effects of the surface capillary pressure. Another important idea of Tolman is that we should consider as simpler as possible the neutron star in terms of the liquid drop at a static equilibrium under the stability condition, i.e., having a minimum of the total energy under constraints mentioned above.

The accuracy of the ESA was checked in Ref. 32 for the case of nuclear physics by comparing the results with the existing nuclear theories like Hartree-Fock[43] (HF) and extended Thomas-Fermi[39, 38] approaches, based on the Skyrme forces [43, 44, 45, 46, 47, 39, 50, 48, 9, 49], but for the simplest case without spin-orbit and asymmetry terms of the energy density functional. The direct variational principle for finding numerically the parameters of the tested particle density functions in simple forms of the Woods-Saxon-like in Ref. 38 or their powers (Ref. 51) were applied by using the realistic Skyrme energy functional[9, 49]. The main focus in Ref. 51 aimed to the surface corrections to the nuclear symmetry energy of spherical nuclei, see also recently published variational ETF approach[52, 53]. The extension of the ES approach to the nuclear isotopic symmetry and spin-orbit interaction has been done in Refs. 33, 34, 35. The Swiatecki derivative terms of the symmetry energy for heavy nuclei[54, 55, 56, 57, 58, 59, 60, 61] were taken into account within the ESA in Ref. 35. The discussions of the progress in nuclear physics and astrophysics within the relativistic local density approach, can be found in reviews Refs. 18, 62; see also Refs. 13, 63, 64. Notice that the simplest nuclear Thomas-Fermi approach (without gradient terms in the energy density) was applied for neutron stars in Ref. 65.

In the present work, we extend the ES approximation of Refs. 32, 33, 34, 35, to the spherical neutron stars by the method working also for deformed systems. In Sect. 2, we present the basic local energy-density formalism taking into account particle density gradients which are responsible for surface terms in finite macroscopic dense systems. The particle number density solutions inside of the nucleus are obtained analytically within the ESA expansion in a small leptodermic parameter $a/R$ at leading order in Sect. 3. The leading density distributions for the density $\rho$ at the lowest order in a relatively small diffuse surface-layer size, $a/R\ll 1$, are obtained analytically in Sect. 4. The NS mass with the surface correction is derived in Sect. 5 taking into account the Schwarzschild metric. The surface energy in terms of the tension coefficients of the vdW macroscopic capillary theory[42, 40, 31, 41] is obtained analytically through vdW-Skyrme forces parameters in Sect. 6. Equation of state in the form $P=P(\rho)$, for the system surface corrections, is derived through the ES at the same leading order in Sect. 7. The TOV approach is presented in Sect. 8 for a step-like particle number density. The results of our calculations and comparison of the macroscopic and some semi-microscopic calculations of the EoS are discussed in Sect. 9. The main results are summarized in Sect. 10. Some details of the mathematical textbook relations will be shown in A.

The total energy $E$ for static problems can be written as

where ${\cal E}(\rho)$ is the energy density[32, 35],

The integration is carried out over the volume of a system, $\hbox{d}\mathcal{V}=e^{\lambda/2}\hbox{d}{\bf r}$, where $\lambda$ is the coefficient of the Schwarzchild metric (1) for a spherical system. As shown in Ref. 3, the multiplier $e^{\lambda/2}$ takes into account the gravitational defect of the NS mass. In Eq. (3), $\mathcal{A}(\rho)$ and $\mathcal{B}(\rho)$ are smooth functions of the density $\rho$ which are coefficients in expansion of the energy density over gradients of $\rho$. A non-gradient part $\mathcal{A}(\rho)$ of the energy density $\mathcal{E}$ can be written as

where $b_{V}$ is the non-gravitational energy component for the separation of particle from the matter, $m$ is the particle mass, and $\Phi$ is the macroscopic gravitational potential determined in more details below. The second and third terms take into account the non-gravitational (like nuclear) and gravitational contributions into the incompressibility. In order to define properly other quantities in Eq. (4), $\varepsilon(\rho)$ and $\Phi(\rho)$, we now use the condition for a minimum of the energy per particle, $\mathcal{W}$, at a stable equilibrium,

Near the saturation value, $\rho\rightarrow\overline{\rho}$, one has $\mathcal{W}=\mathcal{A}/\rho$. Therefore, there is no linear terms in expansion of $\mathcal{A}$, Eq. (4), over powers of the difference $\rho-\overline{\rho}$ near the saturation value for an isolated system including the gravitation. For instance, for $\varepsilon(\rho)$ in Eq. (4), one can write, generally speaking, for the considered dense system (similarly as in Refs. 32, 33),

where $K$ is the non-gravitational part of the incompressibility modulus, e.g., due to the Skyrme nuclear interaction (for nuclear matter $K\sim 200$ MeV). In Eq. (4), as mentioned above, $\Phi(\rho)$ is the main statistically (macroscopically) averaged part of the gravitational potential. In what follows, we may restrict ourselves by only quadratic terms in expansion of the smooth function $\mathcal{A}(\rho)$ due to the saturation property of the NS inside of the dense finite drop of the matter.

As well known[2, 3], the Einstein-Gilbert GRT equations in the non-relativistic limit and/or for a week gravitational field can be transformed to the continuity Poisson equation for the gravitational potential $\Phi$. As shown in Ref. 3, in this limit, $g_{00}=e^{\nu}$ of a more general Schwarzschild metric, Eq. (1), can be presented as $e^{\nu}\approx 1+\nu=1+2\Phi/c^{2}$, where $\Phi$ is the solution of the Poisson equation[3, 7] $\nu=2\Phi/c^{2}$. We will use a more general definition for the gravitational potential, $\Phi=(c^{2}/2)\ln g_{00}=(c^{2}/2)\nu$. Instead of the Newtonian limit expansion, for the statistically averaged potential $\Phi$, we will use another macroscopic approximation within the leptodermic approach. Therefore, after such a macroscopic averaging, $\Phi$ can be considered as function of the radial coordinate $r$ through the particle number density $\rho=\rho(r)$ in the form of expansion of $\Phi(\rho)$ over powers of $\rho-\overline{\rho}$ near the saturation density $\overline{\rho}$ for leptodermic systems. This is similar to the macroscopic Coulomb potential in nuclear physics, considered up to the residual inter-particle interaction[29, 30, 33] but accounting now for second order terms in $\rho-\overline{\rho}$. Up to second order, one can take into account the macroscopically averaged gravitational contribution to the incompressibility modulus $K$,

where $\Phi_{0}=\Phi(\overline{\rho})$, $\Phi_{n}=\partial^{n}\Phi(\overline{\rho})/\partial\rho^{n}$ (n=1,2) are the derivatives of the gravitational potential $\Phi(\rho)$ at the saturation density $\overline{\rho}$. According to Eq. (4) for the non-gradient part $\mathcal{A}$ of the energy density $\mathcal{E}$, and the saturation condition (5), one finds for the zero constant $\Phi_{1}$ in the expansion (7) of the gravitational potential $\Phi(\rho)$. With this condition and the expansion (7), one obtains

Therefore, for $\mathcal{A}(\rho)$, Eq. (4), we arrive at

It was convenient to introduce the two new quantities, $b^{(G)}_{V}$ is the total separation energy per particle, and $\varepsilon_{G}(\rho)$ is the total second-order energy density component, both modified by the gravitational field,

Here, $K_{G}$ is the total incompressibility modulus modified by the gravitational field ($K_{G}>0$),

Notice that our second order approximation for the gravitational potential $\Phi$ expansion, Eq. (8), agrees with the energy density presentation up to second order gradients, Eq. (3). The latter is, in turn, rather general for all known Skyrme forces in nuclear physics taking into account the volume and surface terms. Thus, we account for the contribution of a strong gravitational field in terms of separation energy $b^{(G)}_{V}$, Eq. (9), and the total incompressibility $K_{G}$, Eq. (11), which will be agreed with the Schwarzschild metric space curvature within the ES approximation. As we will see, this second-order potential term is consistently related with the density gradient squared term of the energy density $\mathcal{E}(\rho)$, Eq. (3), by the equilibrium equation at leading order of the leptodermic parameter $a/R$.

The coefficient $\mathcal{B}(\rho)$ in front of the gradient squared term in Eq. (3) is given by

These terms are associated with the nuclear Skyrme interaction. First term is related to the interaction term which is a main reason of the diffuse surface thickness for a non-rotating NS. In this sense it has a more general meaning including the main effective interaction in a dense molecular system, studied by van der Waals[36] (vdW); see also the footnote on page 4. Therefore, we will call this component as the vdW-Skyrme interaction. The second term is coming from the spin-orbit interaction. This interaction might be important, for instance, for any dense rotating liquid-drop systems with a sharp surface edge, i.e., with large density-gradient terms in the surface region for a finite leptodermic systems; see also the same general arguments for using the ETF approach in Refs. 38, 39. The last term is a gradient correction to the kinetic energy ($\hbar^{2}$ correction of the nuclear kinetic energy in the ETF approach[39, 38]). It is introduced here for a comparison of the dense and gas systems. Thus, the forms (9) for $\mathcal{A}$ and (12) for $\mathcal{B}$ are rather general among the simplest analytical solutions for dense finite systems.

We should add the constraint for the variational procedure to get equations for the static equilibrium. For non-rotated one-component NS system we have to fix the particle number $N$,

The integration is carried out over the volume occupied by the gravitational mass and determined by the GRT metric. Therefore, introducing the chemical potential $\mu$ as the Lagrange multiplier, one obtains from Eqs. (2) and (3) the equation for the equilibrium,

where $\mathcal{A}(\rho)$ and $\mathcal{B}(\rho)$ are given by Eqs. (9) and (12), respectively. Equation (13) determines the chemical potential $\mu$ in terms of the particle number $N$. For nuclear liquid drop, one has two equations related to the two constraints for the fixed neutron and proton numbers (Refs. 33, 34, 35). They determine two (neutron and proton) chemical potentials. The Coulomb interaction can be taken into account for a proton part of the nucleus through the Coulomb potential; see Refs. 30, 33. The terms like $\propto\Delta\rho$ with constant of the proportionality are omitted because they do not affect the equilibrium density because of their disappearance from the variational equations for the total energy owing to the well-known Ostragradsky-Gauss theorem. High order derivative terms in the energy per particle $\mathcal{E}/\rho$ [see Eq. (3) for $\mathcal{E}$], are neglected for simplicity in analogy of the GRT, Ref. 2. Equation (3) for $\mathcal{E}$ overlaps most of the Skyrme forces[9, 49]. As function of a local particle density $\rho$, the $\mathcal{E}(\rho)$ corresponds to a saturation condition, Eq. (5). The most remarkable in this form is that the energy density $\mathcal{E}(\rho)$, Eq. (3), is a sum of the three terms related to the incompressible constant, $-b^{(G)}_{V}$, Eq. (9), and the compressible energy, $\propto(\rho-\overline{\rho})^{2}$, both including the gravitational components, Eq. (8), and the surface gradient terms. In nuclear physics, one has, instead of the gravitational term but similarly, the Coulomb potential (even without quadratic terms). Similarly, the symmetry energy, and surface ($\nabla^{2}$) terms can be taken into account in nuclear physics and astrophysics. Within the nuclear ETF approach, the terms being proportional to $\Gamma$ of the gradient part in (3) comes from the $\hbar^{2}$ correction to the TF kinetic energy density[38], $\Gamma=\hbar^{2}/18m$, where $m$ is the nucleon mass. In the gradient squared part, Eq. (12), the $\mathcal{C}$ constant term is typical for the potential component of the Skyrme energy density functional and $\mathcal{D}$ is the constant of the spin-orbit term[39] in the ETF model in nuclear physics. For a nucleus, the spin-orbit coefficient $\mathcal{D}$ is given by $\mathcal{D}=-(9m/16\hbar^{2})W_{0}^{2}$, where $W_{0}=100-130$ MeV fm⁵ is the nuclear spin-orbit constant; see Refs. 39, 9, 49. The Coulomb part of the nuclear energy density (3) is considered similarly as suggested in Refs. 29, 30. Meaning of all terms in the energy density (3) will be specified more below.

For the spherical and deformed system of $N$ particles, we may find the equilibrium particle density $\rho$ from the variational problem for the energy functional (2) with respect to the variations of the density $\delta\rho$, which obey the constraints in the form:

The last constraint fixes a certain deformation parameter. The function $q({\bf r})$, for instance, can be a multipole moment, or it can be chosen in such a way that Q determines the distance between the centers of masses of the two halves of the stretched nucleus for fission problems[66]. Thus, for the energy functional (2), (3) and the constraints (15), one finds the same Lagrange variational equation (14) but $\mu$ would be equal to the sum of the chemical potential and $q\mu_{Q}$, $\mu\Rightarrow\mu+q\mu_{Q}$. Lagrange multipliers $\mu$, and $\mu_{Q}$ are determined by the two constraints (15). Formally, we may consider the same equation (14) for spherical and deformed systems taking into account $q$ in $\mu$.

In the system volume, the terms of Eq. (14) containing derivatives of $\rho$ are small. These derivatives in a normal direction become large near the nuclear edge. For a rather wide class of deformed shapes of the dense finite system, as well as for the axially-symmetric, in particular, spherical one, we may assume that the thickness $a$ of a system edge is small as compared with its mean curvature radius $R$, considering $a/R$ as a small parameter. In this respect, we define the effective surface as the points of maximum of the gradient of a particle density $\nabla\rho$, as shown in Fig. 1 by dots in the right panel. Locally near a point of the ES, one may introduce the local coordinate system $\xi,\eta$, where $\xi$ is normal to the ES, and $\eta$ presents two other orthogonal coordinates; see for instance, as shown in Fig. 1 and A. In particular, for the spherical coordinates, $\eta$ can be two spherical angle coordinates. See A, also for simple geometric relations for the non-linear $\xi,\eta$ coordinates. It naturally appears the mean curvature $H$ in terms of a small leptodermic parameter of the ES approximation, $aH$; see Eq. (79). For heavy enough nuclei near the spherical shape with the effective radius $R=r_{0}A^{1/3}$, one has $aH=a/R\sim A^{-1/3}\ll 1$ because $a/r_{0}\sim 1$ for realistic nuclear parameters, $r_{0}=(4\pi\overline{\rho}/3)^{-1/3}\approx 1.14$ fm at the density of infinite nuclear matter, $\overline{\rho}=0.16$ fm^-3, and the typical diffuseness parameter $a\approx 0.8$ fm, see below more precise definition for the diffuseness parameter through the decrement of decreasing of the exponential asymptotes of particle density. For NSs, it is well known that the typical surface diffuseness (the NS crust) $a\approx 1$ km, and the mean effective radius $R\approx 10$ km; see, e.g., Ref. 20. Therefore, the ratio $a/R$ can be also considered as a small parameter, and the leptodermic approximation, $a/R\ll 1$, can be used too.

The largest terms in Eq. (14) within the region of a sharp density descent are the second-order derivative of particle density $\rho$ in the $\xi$ direction, normal to the ES, $\hbox{d}^{2}\rho/\hbox{d}\xi^{2}$, and its derivative square, $(\hbox{d}\rho/\hbox{d}\xi)^{2}$, both of the order of $(\overline{\rho}/a)^{2}\propto(R/a)^{2}$; see the expression for the Laplacian and gradient in $\xi,\eta$ coordinates in A. Our approach is based on expansion in the same small parameter $aH$ (for a nucleus, $aH=a/R\sim 1/k_{F}R\sim A^{-1/3}$, where $k_{F}$ is the Fermi momentum in units of $\hbar$). This leptodermic approach is used in the liquid-drop and the extended Thomas-Fermi approach. Following Refs. 29, 31, 32, 33, 34, 35, we shall call this statistical²²2The ETF in nuclear physics is the statistical semiclassical approach because for convergence of the expansion in $\hbar$ of the partition function we have first to average statistically it removing strong (e.g., shell) oscillations, and thus, get a smooth behavior of the partition function[38]. semiclassical ETF approach as the ES approximation. We have to evaluate also the derivatives of the particle density and gravitational terms in the energy density $\mathcal{E}(\rho)$, Eq. (3), and Lagrange equation (14). For simplicity, the gravitational terms [see Eqs.(8) and (14)] are assumed to be of the same order as other non-gradient terms. Notice that according to these estimations, the transformation of radial derivatives due to the Schwarzschild radial curvature can be taken into account through the leptodermic parameter $a/R$ because they are important near the ES. The gravitational corrections are in contrast to those due to the Coulomb potential for the nuclear liquid drop, where the Coulomb corrections to a constant are assumed to be negligible as $a/R$. The latter corrections[30] are of the relative order $\sim(e^{2}/\hbar c)(\hbar c/b_{V}r_{0})(Z^{2}/A^{4/3})/2\hbox{\kern 1.00006pt\lower 2.58334pt\hbox{$\sim$} \kern-11.19997pt\raise 2.58334pt\hbox{$<$} \kern 1.00006pt}0.2$.

As function $\mathcal{E}(\rho)$ obeys the saturation property (5), in the system volume at $\xi\hbox{\kern 1.00006pt\lower 2.58334pt\hbox{$\sim$} \kern-11.19997pt\raise 2.58334pt\hbox{$<$} \kern 1.00006pt}\xi_{vol}$ ($\xi<0$ and $|\xi|\gg a$), where $\rho-\overline{\rho}$ is small asymptotically far from the ES, one can expand the non-gradient function $\mathcal{A}$ of the energy density $\mathcal{E}(\rho)$, Eq. (3), up to second order including the compression term of the second order in powers of $\rho-\overline{\rho}$, $\mathcal{A}=\mathcal{A}_{V}$; see Eqs. (9), (10), and (6),

The zero and second order terms of expansion of the gravitational potential, Eq. (8), are included in the zero and second order terms through constants $b^{(G)}_{V}$ [Eq. (9)] and $K_{G}$ [Eq. (11)], respectively.

Introducing the dimensionless quantities for convenience to exclude the transformation of units,

one can present Eq. (16) in the following way:

where

Discarded terms are of the order of $((\rho-\overline{\rho})/\overline{\rho})^{3}\sim(a/R)^{3}$. Neglecting gradient terms in Eq. (14) in the system volume and using the approximation (8) for the gravitational potential $\Phi$, for simplicity, up to quadratic terms over $\rho-\overline{\rho}$, one can reduce equation (14) as

where $\mathcal{M}$ is an integration constant. Solving this equation with respect to the particle number density $\rho$, one obtains

Therefore, $\mathcal{M}$ is the surface (capillary) correction to the leading component of the chemical potential $\mu$ [introduced above as the Lagrange multiplier, determined through the constraint (13)]. As seen from Eq. (21), one finds the relatively small surface correction, $9\mathcal{M}/K_{G}\sim a/R$, to the saturation constant value $\overline{\rho}$; see Refs.  29, 30, 31, 33, 34, 35. Therefore, $\mathcal{M}$ is also the capillary pressure correction; see also Ref. 42.

In this section, we will analytically solve equation (14) for the distribution of the particle density $\rho$ through the ES at leading order in small parameter $a/R$. Equation (14) is a typical differential equation in the catastrophe theory. In such differential equations, a small coefficient (of the order of $a/R$ in our case) appears in front of the highest order derivative. These terms become important because of the product of this small coefficient (of the order of $a/R$) and large (of the order of $R/a$ and high) can be of zero order in this leptodermic parameter. At the leading order, we need to keep only leading terms in small parameter $a/R$ in the Lagrange equation (14). These are the second order derivatives $\rho^{\prime\prime}$, and first-order derivative squares $(\rho^{\prime})^{2}$, over their $\xi$ variable, according to Eqs. (80) and (77), $\rho^{\prime}=\partial\rho/\partial\xi$. We may develop some iteration procedure to find the solutions with improved precision in $a/R$ in terms of the solutions of the leading order. As mentioned above, the transformation of radial derivatives due to a smooth Schwarzschild radial curvature for the spherical case can be taken into account near the ES through the re-definition of the leptodermic parameter $a/R$.

Multiplying Eq. (14) by the derivative $\partial\rho/\partial\xi$, from Eqs. (14), (9), and (12) at leading order, one obtains

where

and $\epsilon(y)$ is a given function of $y$, which is determined by a short- and long-range inter-particle interaction. In particular, one can use a simple quadratic approximation (19) far from the critical point where $K_{G}=0$. In Eq. (22), the term with the coefficient $\mathcal{C}$ in parentheses in front of the derivative squared of the density $\rho$ is the main term, e.g., as related to the vdW capillary theory. In particular, it is the main term, associated with the Skyrme interaction, with respect to the spin-orbit term $\mathcal{D}\rho$ and, moreover, the gradient correction $\Gamma/(4\rho)$ to the kinetic energy of the system in nuclear physics. Let us first assume the same for the NS liquid drop. Therefore, $a$ should be determined through the constant $\mathcal{C}$. Otherwise, we have to include average of whole coefficient in front of the derivative square term in Eq. (22).

With employing the boundary conditions $\rho\to 0$ and $\rho^{\prime}\to 0$ for $\xi\to\infty$, one can easily integrate equation (22). Finally, we come to a simple ordinary 1st order differential equation for $\rho$, depending only on the normal-to-ES coordinate $\xi$, at leading order in a small parameter $a/R$, that is

Introducing several dimensionless quantities,

one can relate the crust thickness $a$ to the vdW interaction constant $\mathcal{C}$ by

Using these definitions, one can present Eq. (24) in the following dimensionless form:

where

From the asymptotes of the explicit analytical expressions for the particle density $y(x)$ at large $x$, one may see a typical behavior, $y(x)\propto e^{-x}$, where $x=\xi/a$, that is a reason to call $a$, Eq. (26), as the crust diffuseness parameter. Thus, we may check that $a$ has exactly the same meaning as that introduced above in terms of the small parameter $a/R$. As noted above, in the derivation of Eq. (27) within leading order of the ES expansion, we reduced the second order differential equation (14) to the first order one by integrating with using the standard boundary conditions $y\to 0$ and $\partial y/\partial\xi\to 0$ and $y\epsilon(y)\to 0$ for the definition of the $\epsilon(y)$ asymptotes. Finally, all parameters of the energy density in Eq. (3) are reduced within this leading order ES approximation to the only two dimensionless constants; see Eq. (28). For nuclear physics,

where $\varepsilon_{F}=\hbar^{2}k^{2}_{F}/2m\approx 37$ MeV is the Fermi energy for the Fermi momentum, $k_{F}=(3\pi^{2}\;\overline{\rho}/2)^{1/3}$ (in units of $\hbar$). We are coming to the realistic nuclear data for the infinite-matter particle density, $\overline{\rho}=0.16$ fm^-3.

Single boundary condition which we need for the unique solution of the first order differential equation can be easily found within the leading ES approximation from the definition of the ES as the points of maximal values of the derivative, $\partial\rho/\partial\xi$ at $\xi=0$, namely, $y^{\prime\prime}(0)=0$ at $x=0$. Differentiating equation (27) over $x$, one obtains the algebraic equation for the position of the ES, $y=y_{0}$ at $x=0$,

defined by function $\epsilon(y)$; see, e.g., Eq. (19).

For any function $\epsilon(y)$, Eq. (27) can be easily integrated. The integration constant is determined by the boundary condition (30). The leading-order solution of Eq. (27) can be found explicitly in the inverse form:

For the expression (19) for $\epsilon(y)$, this integral can be calculated analytically in terms of the elementary functions for any parameters $\beta$ and $\gamma$ within the condition $a/R\ll 1$. Note that the solution (31) to equation (27) is considered at leading order in small parameter $a/R$. This solution satisfies asymptotically the condition of its matching with the volume result, Eq. (21), of the same order at the point $x$ corresponding $\xi\approx\xi_{vol}$ defined in Section 3, $y\to 1$ exponentially for $x\to-\infty$.

Let us consider the examples of simple solutions $x=x(y)$, Eq. (31), of Eq. (27) with the boundary condition, Eq. (30), for the quadratic expression $\epsilon(y)=y(1-y)^{2}$.

(i) For $\gamma=0$ and $\beta=0$, neglecting the gradient correction to the kinetic energy density ($\Gamma=0$ in Eq. (12)), and spin-orbit ($\mathcal{D}=0$ in Eq.(12)) terms, one finally obtains[32, 33]

This solution is related to the gradient squared term due to the main nuclear Skyrme, or molecular van der Waals interaction term [$\mathcal{C}\neq 0$ in Eq.(12)]. It is very asymmetric with respect to the ES value, $y_{0}=y(0)=1/3$; see Fig. 1(b).

(ii) Keeping only $\Gamma\neq 0$ term in Eq. (24), but neglected gradient terms of the vdW-Skyrme ($\mathcal{C}=0$) and spin-orbit ($\mathcal{D}=0$) interaction, one obtains³³3Dimensionless equation in this case is $\hbox{d}y/\hbox{d}x=-(y\epsilon(y))^{1/2}$ for $a=\Gamma[9/(8K_{G})]^{1/2}$. Equation for the ES, $x=0$, $\epsilon(y_{0})+y_{0}\hbox{d}\epsilon(y_{0})/\hbox{d}y=0$, has the solution $y_{0}=1/2$. Integrating analytically the equation for $y(x)$ in this case, one obtains the expression (33). the Wilets solution[28], derived early for the semi-infinite system,

This solution is symmetric with respect to the ES because of $y(0)=1/2$. This Fermi-gas form of the solutions is used for the variational problems within the ETF model, in contrast to the dense liquid-drop solution (32), related to the vdW-Skyrme interaction; see Fig. 1(b).

Figure 1(b) shows the particle density, $\rho(r)=\overline{\rho}y((r-R)/a)$, as function of the radial coordinate $r$ for typical parameters, the thickness of the NS diffused layer (crust) $a=1$ km, and the effective radius $R=10$ km, in units of the saturation value $\overline{\rho}$. We compare the two dimensionless solutions for the main asymmetric particle density $\rho(r)/\overline{\rho}$. One of them (solid line) is related to the corresponding vdW-Skyrme interaction “1”; see Eq. (32) for the universal dimensionless particle density of a dense liquid-drop, $y=\rho/\overline{\rho}$, as function of $x=(r-R)/a$. Another limit case corresponds to the symmetric Wilets solution based on Eq. (33) for zero constants $\mathcal{C}$ and $\mathcal{D}$ in Eq. (12) for $\mathcal{B}$ (only the gradient correction to the kinetic energy is taken into account, i.e. for a gas system). This solution “2” is often used for the variational version of this problem. These two curves are very different in the surface layer, especially outside of the system far away from the NS, and almost the same inside of the NS. Notice that so far the effective surface method can be applied for the deformed NS under the condition of smallness of crust thickness $a$ over the effective NS radius $R$ determined by the mean surface curvature $H$, $aH=a/R\ll 1$; see A. However, in the following sections, our macroscopic method will be applied for the simplest solved analytically case of the spherical symmetry.

Let us derive the NS mass M within our macroscopic ESA approach,

where $m$ is the particle (nucleon, nucleus, or molecule) mass. The integration is carried out over the spatial volume,

This element of the spatial volume is associated with the Schwarzschild interior spherical coordinate system (see Refs. 2, 3),

where $R_{\rm SM}$, $A_{\rm SM}$, and $B_{\rm SM}$ are some constants; see Ref. 2, and more details below. Here and below we will use the subscript SM to relate the quantities to the Schwarzschild metric, Eq. (36). The Schwarzschild radius $R_{\rm SM}$ is given by[2, 3]

where $G$ is the gravitational constant, $\mathcal{E}_{0}=\mathcal{A}(\overline{\rho})=-b^{(G)}_{V}\overline{\rho}$, Eq. (9). Then, according to Eqs. (35) and (36), for the radial part of the Jacobian factor $J$ in this transformation, one has

The Jacobian $J$ takes approximately into account the gravitational defect mass[3].

Adding and subtracting $\overline{\rho}$ in the integrand of Eq. (34), one can identically rewrite the NS mass $M$ in the following way:

Here, the first volume component of the NS mass, $M_{V}$, is given by

where

for $0<z<1$. The second term, $M_{S}$, in Eq. (39) is the surface component of the NS mass,

where $S$ is the surface area of the spherical system, $S=4\pi R^{2}$. The integral over the radial coordinate $r$ is taken effectively over a small diffuse crust region of the order of $a$, where $\rho-\overline{\rho}$ essentially differs from zero. A smooth Jacobian $J(r)$, Eq. (38), multiplied by $r^{2}$, for $r$ far away from the Schwarzschild radius $R_{\rm SM}$, can be taken approximately off the integral in Eq. (42) over the radial coordinate $r$ at the ES $r=R$, $J(r)\approx J(R)$. Using the differential equation (24) for the density $\rho$ at leading order in $a/R$ in the local coordinate system $\xi,\eta$ (A), one can transform the integration variable $\xi$ to $\rho$, $\hbox{d}r=\hbox{d}\xi=(\hbox{d}\xi/\hbox{d}\rho)\hbox{d}\rho$. Thus, the surface correction, $M_{S}$, Eq. (42), in terms of the dimensionless density $y=\rho/\overline{\rho}$, is resulted in

where $J(R)$ is given by Eq. (38) for $J(r)$ at $r=R$, and $\epsilon(y)$ can be parameterized by using the vdW-Skyrme interaction model but with NS parameters. For the simplified (quadratic) expression of $\epsilon(y)$, Eq. (19), the integral in Eq. (43) can be taken analytically in terms of the elementary functions. In particular, for the vdW-Skyrme case (i), $\beta=\gamma=0$, Eqs. (32) and (38) at $r=R$, one explicitly finds from Eq. (43)

Finally, we arrive from Eqs. (39), (40), and (44) at the simple result in the case (i):

where $M_{V}$, Eq. (40), is the volume part of the total NS mass $M$ [see Eq. (39)], i.e., the mass $M$ at $a=0$, and $f(z)$ is given by Eq. (41).

For completeness, one also can present the NS mass with the surface correction in the Wilets case (ii), $\mathcal{C}=\mathcal{D}=0$, Eq. (33),

where $\gamma$ is the dimensionless constant of the gradient correction to the kinetic energy density; see Eq. (28).

Figure 2 shows the mass distribution $M(R)$, Eq. (45), as function of the effective radius $R$ by using only two physical parameters, the relative crust thickness $a/R$, and the central NS particle-number density $\overline{\rho}$ in units of the nuclear matter-saturation density $\rho_{0}=0.16$ fm^-3, $\overline{\rho}/\rho_{0}=1-4$, as examples. The surface effect is measured by the relative difference between the full NS mass, $M=M_{V}+M_{S}$ [see Eq. (45) with Eqs. (40) and (41), solid lines] and its corresponding volume component $M_{V}$ ($a=0$, dashed curves). As seen from Fig. 2, the surface component is rather notable even at a small leptodermic parameter $a/R=0.08$. The significant influence of the gravitational forces through the Schwarzschild metric (36) is shown by comparing the solid [Eq. (41) for $f$] and dotted ($f=1$) green lines, both taken at $a/R=0.08$. It depends also essentially on the value of $a/R$ through the equilibrium condition; cf. dashed-dotted ($a=0.06$) and solid ($a=0.08$) green lines. The NS mass $M(R)$ is the non-monotonic function of $R$ for a constant central NS density, $\overline{\rho}$, because of the surface component, $M_{S}$, Eq. (44). This is in contrast to the monotonic behavior of the volume mass, $M_{V}\propto R^{3}$, in the Cartesian case of the small Newtonian gravitational limit at $a=0$; see also the dashed ($a=0$) at finite $f(R/R_{\rm SM})$ [Eq. (41)] curves. For any given value of $\overline{\rho}$, one finds a rather pronounced maximum in dependence of the full mass $M(R)$, Eq. (45). As well known[2, 3], the physical values of the NS radius for the Schwarzschild metric, Eq. (36) inside of the system, have to obey the condition, $r_{g}<R<R_{\rm SM}$. This means on the left side of the sharp addiction decline, i.e. more pronounceable near the maxima in Fig. 2. These maxima are increased from about 1.2 to 2.5 $M_{\odot}$ for decreasing $\overline{\rho}/\rho_{0}$ from 4 to 1 at $a/R=0.08$ as example, respectively. The NS mass for each of these curves at a given central value, $\rho_{c}=\overline{\rho}$, disappears sharply in the limit $R\rightarrow R_{\rm SM}$, and does not exist at $R>R_{\rm SM}$. As mentioned above, we should emphasize that our derivations for the surface component of the NS mass, $M_{\rm S}$, fail near the point $R=R_{\rm SM}$ because of the obvious singularity of the Jacobian $J(r)$, Eq. (38), at $r=R=R_{\rm SM}$. As seen from Fig. 2, our results are in reasonable agreement with the experimental data[68, 67, 69]; see also Ref. 70 for the theoretical results about the neutron stars. For smaller NS masses, $M\approx(1.2-1.5)M_{\odot}$ for the NS pulsars GW170817 (Ref. 68) and J0030+0451 (Ref. 67) with the slightly different radii, $R=10.3-11.8$ km (orange spot) and $R=11.6-13.1$ km (red spot), one finds respectively good agreement with our results for the central density $\overline{\rho}=3\rho_{0}$, with $a/R=0.08$. Larger mass of the NS pulsar J0740+6620 (Ref. 69), $M=(2.0-2.1)M_{\odot}$ and radius $R=11.3-13.6$ km (green spot), corresponds to the same central density, $\rho_{c}=\overline{\rho}$, but with smaller leptodermic parameter $a/R=0.06$.

Figure 3 shows contour plots for the NS mass, $M$ (in Solar units) as functions of its radius, $R$, and central density, $\rho_{c}$, ($\overline{\rho}$ in units of $\rho_{0}=0.16$fm^-3) for the leptodermic parameter values, $a/R=0.08$ (left) and $0.06$ (right). As seen from the left figure, one for instance finds the line for the NS mass $1.5M_{\odot}$ for radius $R=12.6-14.0$ km and central density, $\overline{\rho}/\rho_{0}=1.8-2.7$. In the right figure we can see the line for the NS mass $2.0M_{\odot}$ for radius $R=10.9-14.0$ km and central density, $\overline{\rho}/\rho_{0}=1.8-4.3$. They are close to the experimental data shown in Fig. 2. The NS mass structure is changed significantly with radius $R$ and relative central density $\rho_{c}$, also with the leptodermic parameter. The calculations become unstable with approaching the line $R=R_{\rm SM}$, as shown by white areas. The physical region is below this line because of the restriction by the condition $r_{g}<R<R_{\rm SM}$.

Let us derive now the NS total energy $E$ which is similar to the main two (volume and surface) terms of the Weizsäker mass formula, basic in nuclear physics, but with the analytical expression for its tension coefficient and accounting for the gravitational forces. According to the basic definitions, the total energy $E$, Eq. (2), can be presented as

where $\hbox{d}\mathcal{V}$ is given by Eqs. (35) and (38), and $E_{V}$ is the volume part of the total energy. Following the mass derivations in the previous section, one has

where $b^{(G)}_{V}$ is given by Eq. (9), and $f(z)$ is shown in Eq. (41). In Eq. (47), $E_{S}$ is the surface part,

where $S=4\pi R^{2}$ is the surface area value, and $\sigma$ is the tension coefficient. The leading expression for $\sigma$ in our leptodermic expansion over small parameter $a/R$ is given by⁴⁴4 Similar expressions as Eqs. (50) and (51) at $\mathcal{D}=\Gamma=0$ were obtained early by van der Waals; see Ref. 36.

where $J(R)$ is the Schwarzschild metric Jacobian $J(r)$ at $r=R$; see Eq. (38). In these derivations, we present locally the integration over $\hbox{d}\mathcal{V}$ as that over the system surface and the normal coordinate $\hbox{d}\xi$, $\hbox{d}\mathcal{V}=J\sqrt{g}\hbox{d}\eta\hbox{d}\xi$, where $\sqrt{g}$ is the Newtonian volume element related to the metrics given in Eq. (77) [see Eq. (76) for $\sqrt{g}$] and Eq. (38) for the Jacobian $J(r)$. For the spherical system, one has $\sqrt{g}=1$. Then, we expand the integration limits of $\xi$ to $\pm\infty$ because $(\partial\rho/\partial\xi)^{2}$ has a sharp maximum at the ES ($\xi=0$). Factor 2 appears because of Eq. (24) which leads to the approximate (at leading order in small parameter $a/R$) equivalence of the second order [$\varepsilon_{G}(\rho)$] correction to the bulk energy density $-b^{(G)}_{V}\overline{\rho}$, and the surface (gradient) parts of the total energy $E$. As we have square of density derivative over $\xi$ in Eq. (50), one can use Eq. (24) for the density $\rho$ at the leading order. Equation (50) can be transformed to the expression with the dimensionless integral over $y$, that is more convenient for calculations,

where $J(R)$ is given by Eq. (38) at $r=R$. In these derivations we transform the variable $\xi$ to $\rho$, $\hbox{d}\xi=(\hbox{d}\xi/\hbox{d}\rho)\hbox{d}\rho$, as above, and use Eqs. (23), (24), and (26)-(28). For the approximation (19) to $\epsilon(y)$, the integral in Eq. (51) can be taken in terms of the elementary functions. The simplest analytical expression for the tension coefficient $\sigma$ is obtained for $\beta=\gamma=0$ [case (i)]. Using Eq. (32) for $y=y(x)$, one finally arrives at

Similarly, the expression for $\sigma$ in the Wilets case (ii) can be derived too, $\sigma=a\overline{\rho}K_{G}J(R)\sqrt{\gamma}/36$. Thus, one obtains the explicitly analytical expressions for the tension coefficient $\sigma$ as functions of constants of the energy density $\mathcal{E}$, Eq. (3), and Schwarzschild metric Jacobian at the ES, $J(R)$. For nuclear physics, the tension coefficient $\sigma$ is related to the vdW-Skyrme interaction constants, in good agreement with the van der Waals capillary theory[36].

Equation (2) is the two main (volume and surface) terms of the usual leptodermic expression for the total energy of the dense liquid drop. In nuclear physics, this expression was obtained with including additional terms of the symmetry energy [with volume (and Swiatecki linear term) and surface (isovector) corrections] and curvature and Coulomb components (see, for instance, Refs. 34, 35 for the ESA).

Notice that due to the integration over the normal-to-ES radial variable $\xi=r-R$, the smallness factor, proportional to $a/R$, appears for the surface with respect to the volume contributions to the NS mass, Eq. (44), and total energy $E$, Eq. (47). Thus, these surface components, $M_{S}$ and $E_{S}$ are of the order of a small parameter $a/R$. However, these surface contributions, $M_{S}$ and $E_{S}$, dramatically influences on the total NS characteristics because the vdW capillary surface pressure (including the gravitational forces) equilibrates the volume pressure acting from the interior of the liquid-drop to be a leading reason of its equilibrium stability. We will study this stability condition in more details in the next section.

The equation of state (EoS), $P=P(\rho)$, can be determined[71] by using the energy density $\mathcal{E}(\rho)$, Eq. (3),

The energy per particle $\mathcal{W}$ is a function of the particle number density $\rho$; see Eqs. (5) for $\mathcal{W}$, and (3) for the energy density $\mathcal{E}$. According to Eq. (14) for the equilibrium of the dense liquid-drop system, one obtains the pressure $P$ in the form:

where $\mu$ is the chemical potential, and $\mathcal{E}(\rho)$ is the energy density, Eq. (3). The leading term of the pressure (54) can be presented as