LEPTODERMIC CORRECTIONS TO THE TOV EQUATIONS AND NUCLEAR ASTROPHYSICS WITHIN THE EFFECTIVE SURFACE APPROXIMATION

The macroscopic effective-surface approach (ESA) wilets ; strtyap ; tyapin ; strmagbr ; strmagden ; MS09 ; BM13 ; BM15 taking into account the gradient terms of the energy density was applied MS24 for the neutron stars. Following the works of Tolman, Oppenheimer, and Volkoff (TOV) RT39 ; OV39 we consider the neutron star (NS) as a liquid drop at equilibrium under the gravitational, and other realistic forces; see also the Tolman and Rowlinson books RT87 ; RW82 . As shown in Ref. RT87 , the TOV equations were derived OV39 from the equations of the General Relativistic Theory (GRT) assuming the spherical symmetry and the macroscopic approximation to the energy-momentum tensor. The expression for the energy-momentum tensor was essentially used in the TOV derivations for a macroscopic system as the perfect liquid drop LLv2 ; LLv6 . Many applications in terms of the Schwarzschild metric solutions in the spherical case can be found in Refs. LLv2 ; WF88 ; CB97 ; SH06 ; Ko08 ; FCPG13 ; AB14 ; BC15 ; LH19 ; SBL23 . For the basic relations of the nuclear and neutron liquid-drop models to NS properties we should also mention a pioneer work BBP71 done by Baym, Bethe, and Pethick. They discussed 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. See also Refs. LLv2 ; RT87 ; OV39 ; SBL23 ; LP01 ; HP01PRL ; HP01PRC ; PH00Poland ; ABCG09 ; PFCPS13 ; FCPG13 ; FG23 ; DFG23 ; HPY07 ; LH19 ; ST04 for the description of the finite NS dense matter systems in terms of the leptodermic approximation, $a/R\ll 1$, where $a$ is the crust thickness, and $R$ is the mean curvature radius of the ES. 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. SBL23 .

Using the step-function for the particle number density $\rho$, which is a constant inside and zero outside of the system, in line of the Schwarzschild derivations KS16 , Tolman obtained RT87 the explicit analytical solutions for the pressure. They were found in terms of the Schwarzschild metric coefficients for the spherical system symmetry, and the macroscopic energy-momentum tensor for the perfect liquid drop. Analytical expressions for the pressure in terms of the simple one-dimensional integrals can be obtained also for the leptodermic particle-number density $\rho$, e.g., as shown in Fig. 1.

In this work we derive the modified TOV (MTOV) equations by taking into account the gradient surface corrections within the leptodermic approximation Refs. LJ23 ; LP01 ; HP01PRL ; HP01PRC ; PH00Poland ; ABCG09 ; PFCPS13 ; FCPG13 ; FG23 ; DFG23 ; HPY07 ; LH19 ; MS24 . We consider first a dense one-component system as in Ref. MS24 , and take the energy density functional $\mathcal{E}$ in a rather general form, linear in $(\nabla\rho)^{2}$ terms, with the coefficients which are smooth functions of the particle density $\rho$. In such a system, the particle density $\rho$ in the leading order approximation is function of the radial coordinate with exponentially decreasing behavior from almost a constant saturation value inside of the dense system to that of zero through the NS effective surface in a relatively small crust range $a$; see, e.g., Fig. 1. The ES is defined as the inflection points of spatial coordinates, i.e., with density gradients maximum. To obtain the analytical solutions for the particle number density and Equation of State (EoS), we use the leptodermic approximation. In this approximation, the gravitational forces can be taken into account through the macroscopic gravitational potential $\Phi$ determined by the second order expansion over deflection of the density $\rho$ from its saturation value $\overline{\rho}$ (see, e.g., Ref. MS24 ).

Within the leptodermic approximation, $a/R\ll 1$, simple and accurate solutions of many nuclear and liquid-drop problems involving the particle number density distributions were obtained by using the ES approach for nuclei wilets ; strtyap ; tyapin ; strmagbr ; strmagden ; MS09 ; BM13 ; BM15 ; see also Ref. MS24 for applications to the neutron stars. 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 such as molecular systems, liquid drops (amorphous solids), nuclei, and presumably, NSs. The realistic energy-density distribution per particle is minimal at a certain saturation density $\overline{\rho}$, corresponding approximately to the infinite dense matter bete . 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 a given point to the ES, and others ($\eta$) are tangent to the ES. This nonlinear coordinate system $\xi,\eta$ is conveniently used in a region near the nuclear and NS edge. They allow for an easy extraction of the relatively large terms in density distribution equations for the condition of a variational equilibrium. 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, dense molecular, and metallic-cluster physics brguehak ; brbhad . 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$). Such an equation mainly determines approximately the density distribution across the diffused surface layer. The latter is of the relatively small order over the ratio of the diffuseness parameter $a$ to the (mean curvature) ES radius $R$, $a/R<<1$, for sufficiently heavy systems. A small parameter, $a/R$, of the expansion within the ES approach can be used for analytical solving the variational problem for the minimum of the energy density functional per volume 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 the 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 liquid-drop model boundary conditions bormot ; strmagbr ; magstr . Our ESA 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 Ti52 . A relatively large change of the density $\rho$ on a small distance $a$ with respect to the curvature radius $R$ takes place for a finite liquid-matter drop (nuclei, molecular water drops, neutron stars). We emphasize that the macroscopic ESA theory MS24 allows to consider the deformed and superdeformed systems for which the same coordinate system $\xi,\eta$ can be used locally near the ES. Inside of such dense systems, the density $\rho$ is changed a little around constant saturation density $\overline{\rho}$. Therefore, including also the second order expansion of the gravitational potential near $\overline{\rho}$, one obtains essential effects of the surface capillary pressure.

The accuracy of the ESA was checked in Ref. strmagden for the case of static and dynamic nuclear physics by comparing the results with the existing nuclear theories like Hartree-Fock (HF) vauthbrink and ETF brguehak ; brbhad , based on the Skyrme forces vauthbrink ; skyrme ; barjac ; ringshuk ; blaizot ; brguehak ; gramvoros ; krivin ; CB97 ; CB98 , 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 number density functions in simple forms of the Woods-Saxon-like in Ref. BB03 or their powers (Ref. kolsan ) were applied by using the realistic Skyrme energy functional CB97 ; CB98 ; see also recently published variational ETF approach KS18 ; KS20 . The extension of the ES approach to the nuclear isotopic symmetry and spin-orbit interaction has been done in Refs. MS09 ; BM13 ; BM15 . The Swiatecki derivative terms of the symmetry energy for heavy nuclei myswann69 ; myswiat85 ; danielewicz2 ; vinas1 ; vinas2 ; vinas4 ; vinas5 ; Pi09 were taken into account within the ESA in Ref. BM15 . The discussions of the progress in nuclear physics and astrophysics within the relativistic local density approach, can be found in reviews Refs. SBL23 ; NVR11 ; see also Refs. CH08 ; Pi21 ; Pe23 .

In the present work we apply the effective surface approach for the neutron stars of Ref. MS24 , based on Refs. strmagden ; MS09 ; BM13 ; BM15 of nuclear physics, to derive the MTOV equations taking into account the surface leptodermic corrections. In Sect. II we present the TOV derivations at the zero order approach. In Sect. III we derive the MTOV equations with accounting for the first order surface corrections to the TOV equations. The derivations of the analytical solutions to the MTOV equations are obtained too. Sect. IV is devoted to the discussions of the results. These results are summarized in the conclusion section V. Some details of the calculations are given in Appendix A.

Starting from the Einstein-Gilbert equations of the GRT for the spherically symmetric system, one has the Schwarzschild gravitational metric in a general form:

where $c$ is the speed of light. Functions $\lambda=\lambda(r)$ and $\nu=\nu(r)$ are given by the differential equations (see Ref. RT87 at zero cosmological constant $\Lambda$). First, outside of the system on big distances from masses, or for relatively small gravitational forces RT87 ; LLv2 , one can derive these equations using the macroscopic properties of a dense spherical system for the energy-momentum tensor. Transforming then the radial coordinate $r$ of the spherical four-coordinate system $r,\theta,\phi$, and $t$, due to the GRT invariance, for the line element in this space, Eq. (1), one obtains the simple differential equations for $\lambda$ and $\nu$ of the Schwarzschild metric for the region inside of the liquid-drop system (see, Refs. RT87 ; LLv2 ),

where $G$ is the gravitational constant. Primes above $\lambda$ and $\nu$ mean the corresponding derivatives over the radial coordinate $r$, and $\mathcal{E}$ is the energy density. With the EoS, $\mathcal{E}=\mathcal{E}(\rho)$ [or, $P=P(\rho)$], one arrives at the complete system of equations for $\lambda$, $\nu$, $\mathcal{E}$ (or density $\rho$), and $P$.

Using the substitution $e^{-\lambda}=f$, one can solve the second equation in Eq. (2), as a linear first-order differential equation over $f(r)$, in terms of the quadratures. Then, the derivative $\nu^{\prime}$ can be found from the first equation of Eq. (2). Integrating it over $r$, one obtains the Schwarzschild metric (1) in terms of these solutions for $\lambda(r)$ and $\nu(r)$. Substituting the derivative $\nu^{\prime}$ into the third equation of Eq. (2), one finally obtains the famous TOV equations RT39 ; OV39 ; LJ23 ,

Here, $m(r)$ is the gravitational mass interior to the radius $r$, and $\mathcal{E}$ is a given energy density. Boundary conditions for these equations are given by

The integrations are terminated when $P=0$, which defines the surface $r=R$. A specified value of the central pressure $P_{0}=P(r=0)$ determines the NS volume mass $M=m(r=R)$. Notice that the surface corrections to the NS volume mass were obtained in Ref. MS24 .

As shown in Ref. RT87 , taking approximately into account the step-function density, $\rho=\overline{\rho}$ for $r<R$ inside of the system, and $\rho=0$ for $r>R$ outside of it, one can solve TOV equations (II) explicitly analytically (see Refs. RT87 ; MS24 ). (For instance, $\overline{\rho}$ is the constant density $\rho_{00}$ using units $G=c=1$, inside of the system in the notations of Tolman). Our leading order solutions for the density $\rho$ are obtained in Ref. MS24 by using the ESA within the leading leptodermic approximation; see Eqs. (31)^∗ for more general inter-particle interaction and (32)^∗ for the specific vdW-Skyrme forces. They have a sharp transition from the saturation to zero value, which agrees with above mentioned approximation at zero order of the leptodermic expansion in small parameter $a/R\ll 1$ (see Fig. 1). Taking also into account the first boundary condition in Eq. (4), the second equation in Eq. (II) can be integrated analytically,

where $\mathcal{E}_{0}$ is a constant at $r<R$,

and $\overline{\rho}$ is the saturation approximation for the density $\rho$. We introduced here the Heaviside step function, $\Theta(x)$, $\Theta(x)=1$ for $x\geq 0$ and $0$ for $x<0$. The function $\mathcal{A}(\rho)$ in Eq. (6) is given by Eq. (9)^∗ of Ref. MS24 for $\mathcal{A}(\rho)$ at the saturation density $\rho=\overline{\rho}$; see Table 1. The density $\overline{\rho}$ and energy density $\mathcal{E}_{0}$ can be taken approximately as those found in Ref. MS24 ; see Eqs. (3)^∗, (9)^∗, (10)^∗, and (11)^∗, including the macroscopic gravitational field. Formally, one can consider the expression (9)^∗ for $\mathcal{E}_{0}(\overline{\rho})$ inside of the system [or Eq. (16)^∗] as the macroscopic EoS additional to the TOV equations (II). Therefore, the total mass $M$ neglecting the surface effects and the influence of the Schwarzschild metric curvature, is given within this zero order ESA approximation by

where $R$ is the effective NS radius ($r_{1}$ in the notations of Tolman’s book, Ref. RT87 ).

For the Schwarzschild metric in the interior region, one has

where $R_{\rm S}$, $A_{\rm S}$, and $B_{\rm S}$ are the specific integration constants, and $r<R_{\rm S}$. This metric, Eq. (II), is nonsingular at $r=0$, in contrast to the original Schwarzschild metric RT87 ; LLv2 , which is valid outside of the gravitating system. Therefore, it can be applied for the part of the space inside of this system (see Ref. RT87 ). A smooth match of the external and internal Schwarzschild metrics through the effective radius, $r=R$, is assumed to be carried out. The radius $R_{\rm S}$ in Eq. (II) is expressed in terms of the energy density $\mathcal{E}$ inside of the system, $\mathcal{E}_{0}$, far away from the ES as

Other constants, $A_{\rm S}$ and $B_{\rm S}$, as well as $R_{\rm S}$, are the following parameters of the transformed Schwarzschild metric, Eq.  (II) (see the derivations in Ref. RT87 ),

Substituting Eq. (5) into the first TOV Eq. (II), one can see that the variables $P$ and $r$ are separated. Therefore, we analytically find the solution for $r\leq R\leq R_{S}$:

where

and $P_{0}\equiv P(\overline{\rho})=P(r=0)$. Here, $R_{\rm S}$ is the radius parameter, Eq. (9), in the Schwarzschild metric, Eq. (II).

Notice that the solution, Eq. (11), coincides, up to a constant, with that presented in Tolman’s book at $\zeta=B_{\rm S}/A_{\rm S}=1/(2A_{\rm S})$ (see Refs. RT87 ; MS24 and Eq. (9) for $R_{\rm S}$). The cosmological constant also is assumed to be zero, and the units for the dimensionless pressure are used, $\mathcal{E}_{0}\propto R^{-2}_{\rm S}$, Eq. (9). Thus, one obtains the relation of these parameters of the Schwarzschild metric, Eq. (II), to the initial condition $P=P_{0}$ at $r=0$ through Eq. (12) for $\zeta$. Notice that according to Eq. (11) for the pressure $P$, one has $P=0$ at the boundary $r=R$ because of Eq. (10) for the constants $A_{\rm S}$ and $B_{\rm S}$. Finally, our result for the pressure, Eq. (11), coincides identically with that of Ref. RT87 .

Figure 2 shows the contour plots for the pressure $P(r)$ [in units of $P(r=0)$], as function of $r$ (km) and parameter $A_{\rm S}$, Eq. (10). It is more convenient to use the finite dimensionless values of $A_{\rm S}$, $0\leq A_{\rm S}\leq 1.5$, instead of dimensional $R_{\rm S}$ for the fixed effective NS radius $R=10$ km, $0\leq r\leq R$. For a very rough estimate neglecting the relatively small gravitational mass defect, $\mathcal{E}_{0}\sim c^{2}\rho^{(0)}$, where the central mass density $\rho^{(0)}\approx M/(4\pi R^{3}/3)\sim M_{\odot}/(4\pi R^{3}/3)$, ($M\approx M_{V}$) and $M_{\odot}=1.989*10^{30}$ kg is the Sun mass, one finds significantly larger Schwarzschild radius parameter $R_{\rm S}$ than the effective radius $R$, $R_{\rm S}\approx 15.6$ km for $M=1.4M_{\odot}$, as example. By definition, $R_{\rm S}$ is always larger than (or equal to) the effective NS radius $R$, $R\leq R_{\rm S}$, while the gravitational radius $r_{\rm g}$ should be smaller than $R$, $r_{\rm g}=2GM/c^{2}\leq R$ LLv2 . These conditions are fulfilled in our calculations for the values $R_{\rm S}\approx 15.6$ km, and $r_{\rm g}\approx 4.1$ km in the case of the effective radius $R=10$ km. With these evaluations, one has the constant $A_{\rm S}\approx 1.15$. The maximum value $A_{\rm S}\approx 1.5$ corresponds to zero pressure $P(r)$.

Figure 3 shows the two cuts of the contour plots (Fig. 2) at the value, $A_{\rm S}\approx 1.15$ ($R_{\rm S}\approx 15.6$ km, solid black line), and the value, $A_{\rm S}=0.96$ ($R_{\rm S}=13.0$ km, dashed red line). The physical mass region, $M/M_{\odot}\approx 1.4-2.0$, in line with the experimental data for NS masses and radii (see, Refs. GR21 ; CC20 ; TR21 ; OL20 ), is related to a small range in Fig. 3. These curves have the same starting value 1 and final value 0 for the ratio of pressures $P(r)/P(0)$. For more realistic values of $R_{\rm S}$, for instance, taking into account the gravitational defect of mass LLv2 in the evaluations of the central energy density $\mathcal{E}_{0}$, one may have different values of $A_{\rm S}$ presented in Fig. 2. Thus, as seen from the solid and dashed lines of Fig. 3 and assumed in the derivations of the TOV equations OV39 ; RT87 , the pressure $P$ turns, indeed, to zero at the boundary of the NS system, $r=R$.

We now derive the leptodermic corrections to the TOV equations due to the gradient terms of the energy density $\mathcal{E}$, in line of Ref. MS24 . The details of these derivations are shown in Appendix A.

Figures 5 and 6 show the dimensionless pressure $p_{0}$, Eq. (47), at zero order and the total pressure $p$, Eq. (50), at the first order over a leptodermic parameter $a/R=0.1$ as function of the radial coordinate $r$ and Schwarzschild metric constant $A_{\rm S}$, relatively. The first order contribution $p_{1}$ of the pressure $p$ is proportional to the dimensional incompressibility parameter $\kappa$, Eq. (49). As seen from the comparison of Figs. 6 ($\kappa=10$) and 5 ($\kappa=0$), the effect of the gravitational field on the pressure is significant near the ES ($r\hbox{\kern 1.00006pt\lower 2.58334pt\hbox{$\sim$} \kern-11.19997pt\raise 2.58334pt\hbox{$>$} \kern 1.00006pt}R$) for this change because of the ES correction $p_{1}$, Eq. (48) at a large dimensionless incompressibility $\kappa$, Eq. (49).

Figures 7 and 8 show the cuts of the contour plot of Fig. 6 at $A_{\rm S}\approx 1.15$ corresponding approximately to the typical NS mass $M=1.4M_{\odot}$ and $A_{\rm S}\approx 0.96$ for $M=2.0M_{\odot}$, respectively, as function of the radial coordinate $r$ (in km). These NS masses $M$ and approximately the radii by the order of the magnitude were found in the recent experiments presented in Refs. GR21 ; CC20 ; TR21 ; see also the comparison with the theoretical results in Ref. OL20 . We specify different values of the incompressibility parameter $\kappa$ from a weak (almost nuclear $\kappa=1$) to a strong value ($\kappa=10$) of the gravitational potential. As mentioned above, the value $\kappa=0$ is related to the pressure $p_{0}$ over zero order leptodermic parameter $a/R$; see Fig. 5. For example, one can consider stronger gravitational field up to the value $\kappa=10$ in Fig. 8. For nuclear physics, one has $\kappa\equiv K/(12b_{\rm V})=1.25$ for $K=240$ MeV, and $b_{\rm V}=16$ MeV (or $K=225$ MeV and $b_{V}=15$ MeV from Ref. AB14 ).

Figure 9 shows the EoS (a,c), $p(y)$, for small [$\kappa=1$, (a)] and large [$\kappa=10$, (c)] gravitational forces, and the same $p(y(r))$ (b,d) but with $y=y(r)$ for the vdW-Skyrme interaction, Eq. (32)^∗. Again, we compare two characteristic cases, small ($M=1.4M_{\odot}$, red solid) and relatively large ($M=2.0M_{\odot}$, blue frequent dashed) NS masses. The EoS pressures $p=p(y)$, Fig. 9(a,c), are increasing functions of the density $y$ while the pressures $p=p(y(r))$, (b,d), are decreasing as functions of the radial variable $r$ because of the leptodermic properties of the particle number density $y(r)$ for the vdW-Skyrme forces, Eq. (32)^∗. The surface ($S$) components, $p_{1}$, Eq. (48), are shown by black solid ($M=1.4M_{\odot}$) and green long dashed ($M=2.0M_{\odot}$) lines, respectively. The surface effects become much more significant with increasing gravitational parameter $\kappa$. These results for the pressure are basically in agreement with those presented in Fig. 4 for the macroscopic EoS. However, the surface correction $p_{1}$ leads to a more wide coordinate dependence of the total pressure $p$ near the ES (Fig. 9 (b,d)). As mentioned above, the NS masses $M$ and effective radii $R$ correspond largely to very accurate experimental data; see Refs. GR21 ; CC20 ; TR21 ; OL20 ). Notice also that the mass dependence is also essential for all results presented in Fig. 9 for the volume and surface contributions especially for a large gravitational field. We should emphasize that our results for the volume ($p_{V}$) contribution to the total pressure $p$ are in good agreement with the numerical calculations based on the TOV equations discussed in Ref. MS24 .

The effective surface approximation based on the leptodermic expansion over a small parameter - a ratio of the crust thickness $a$ to the effective radius $R$ of the system - is extended for the description of neutron star properties. The neutron star was considered as a dense finite liquid (including amorphous solid) system at the equilibrium. The mean gravitational potential $\Phi(\rho)$ was taken into account in the simplest form as expansion over powers of differences $\rho-\overline{\rho}$, where $\overline{\rho}$ is the saturation density, up to second order in terms of the particle separation energy and incompressibility. Taking into account the gradient terms of the energy density in a rather general form we analytically obtained the leading (over a small parameter $a/R\ll 1$) particle number density $\rho$. The density $\rho$ was derived as function of the normal-to-ES coordinate $\xi$ of the orthogonal local nonlinear-coordinate system $\xi,\eta$ through the effective surface. With its help, one finds the equation of state (EoS) for the pressure $P=P(\rho)$. The TOV equations were re-derived within the same ESA approximation taking into account the leptodermic surface corrections. Our results for the leading order volume pressure components are in good agreement with those of Tolman RT87 and numerical calculations MS24 . We obtained the first order leptodermic surface corrections owing to the Schwartzschild metric, and corresponding modified TOV equations, and have solved them analytically within the same order. These first order corrections, being of the order of $a/R$, are also proportional to the dimensionless incompressibility depending essentially on the gravitational interaction constant, along with the nuclear interaction part. The surface pressure corrections are significantly increased with increasing NS mass, and even more pronouncely, with the gravitational contribution to the total incompressibility.

As perspectives, one can generalize our analytical approach to the the many-component and rotating systems. It is especially interesting to extend this approach to take into account the symmetry energy of the isotopically asymmetric systems, and collective NS rotations. As this approach was formulated in the local nonlinear system of coordinates $\xi,\eta$ for deformed and superdeformed shapes of the NS surface, one can apply our method to the pair of neutron star pulsars at their large angular momenta.