The macroscopic model for a neutron star NS as a liquid drop at the equilibrium is used to extend the TolmanOppenheimerVolkoff TOV equations taking into account the gradient terms responsible for the system surface. The parameters of the Schwarzschild metric in the spherical case are found with these surface corrections in the leading zero order of the leptodermic approximation aRltlt1 where a is the NS effectivesurface ES thickness and R is the effective NS radius. The energy density mathcalE is considered in a general form including the functions of the particle number density and of its gradient terms. The macroscopic gravitational potential Phirho is taken into account in the simplest form as expansion in powers of rhooverlinerho  where overlinerho is the saturation density up to second order in terms of its contributions to th separation particle energy and incompressibility. Density distributions rho across the NS ES in the normal direction to the ES which are derived in the simple analytical form at the same leading approximation was used for the derivation of the modified TOV MTOV equations by accounting for their NS surface corrections. The MTOV equations are analytically solved at first order and the results are compared with the standard TOV approach of the zero order.
['General Relativity and Quantum Cosmology', 'Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Solar and Stellar Astrophysics', 'Nuclear Theory', 'F.2']
Leptodermic corrections to the TOV equations and nuclear astrophysics within the effective surface approximation
{'LEPTODERMIC CORRECTIONS TO THE TOV EQUATIONS AND NUCLEAR ASTROPHYSICS WITHIN THE EFFECTIVE SURFACE APPROXIMATION': "A.G. Magner, 1, 2, ∗ S.P. Maydanyuk, 1 A. Bonasera, 2  \nH. Zheng, 3 T. Depastas, 2 A.I. Levon, 1 and U.V. Grygoriev 1  \n1 Institute for Nuclear Research, 03028 Kyiv, Ukraine  \nCyclotron Institute, Texas A&M University, College Station, Texas 77843, USA  \n2 3 School of Physics and Information Technology,  \nShaanxi Normal University,Xi'an 710119, China  \n(Dated: November 27, 2024)  \nThe macroscopic model for a neutron star (NS) as a liquid drop at the equilibrium is used to extend the Tolman-Oppenheimer-Volkoff (TOV) equations taking into account the gradient terms responsible for the system surface. The parameters of the Schwarzschild metric in the spherical case are found with these surface corrections in the leading (zero) order of the leptodermic approximation a/R << 1, where a is the NS effective-surface (ES) thickness, and R is the effective NS radius. The energy density E is considered in a general form including the functions of the particle number density and of its gradient terms. The macroscopic gravitational potential Φ( ρ ) is taken into account in the simplest form as expansion in powers of ρ -ρ , where ρ is the saturation density, up to second order, in terms of its contributions to the separation particle energy and incompressibility. Density distributions ρ across the NS ES in the normal direction to the ES, which are derived in the simple analytical form at the same leading approximation, was used for the derivation of the modified TOV (MTOV) equations by accounting for their NS surface corrections. The MTOV equations are analytically solved at first order and the results are compared with the standard TOV approach of the zero order.", 'I. INTRODUCTION': "The macroscopic effective-surface approach (ESA) [1-8] taking into account the gradient  \nterms of the energy density was applied [9] for the neutron stars. Following the works of Tolman, Oppenheimer, and Volkoff (TOV) [10, 11] we consider the neutron star (NS) as a liquid drop at equilibrium under the gravitational,  \nand other realistic forces; see also the Tolman and Rowlinson books [12, 13]. As shown in Ref. [12], the TOV equations were derived [11] 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 [14, 25]. Many applications in terms of the Schwarzschild metric solutions in the spherical case can be found in Refs. [14, 16-19, 21-24, 34]. For the basic relations of the nuclear and neutron liquid-drop models to NS properties we should also mention a pioneer work [15] 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. [11, 12, 14, 20, 22, 23, 27, 2936, 69] for the description of the finite NS dense matter systems in terms of the leptodermic approximation, a/R /lessmuch 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 interparticle forces in the non-relativistic and relativistic cases can be found in the recent review Ref. [23].  \nUsing the step-function for the particle number density ρ , which is a constant inside and  \nzero outside of the system, in line of the Schwarzschild derivations [26], Tolman obtained [12] 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 ρ , e.g., as shown in Fig. 1.  \nIn this work we derive the modified TOV (MTOV) equations by taking into account the gradient surface corrections within the leptodermic approximation Refs. [9, 20, 22, 27-36]. We consider first a dense one-component system as in Ref. [9], and take the energy density functional E in a rather general form, linear in ( ∇ ρ ) 2 terms, with the coefficients which are smooth functions of the particle density ρ . In such a system, the particle density ρ 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  \nTABLE I. Equations of Ref. [9] (the first column), marked by the asterisk in the text, are presented for the quantity with symbols in the second column, and explained in the third column. The vdW-Skyrme is the van der Waals and Skyrme inter-particle interaction.  \naccount through the macroscopic gravitational potential Φ determined by the second order expansion over deflection of the density ρ from its saturation value ρ (see, e.g., Ref. [9]).  \nWithin the leptodermic approximation, a/R /lessmuch 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 [1-8]; see also Ref. [9] for applications  \nto the neutron stars. The ESA exploits the property of saturation of the nuclear particle density ρ 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 ρ , corresponding approximately to the infinite dense matter [37]. As a result, relatively narrow  \nFIG. 1. Particle number density, ρ , in units of the saturation value, ρ , as function of the radial variable r (in km) for a neutron star in a simple compressed liquid-drop model at the stable equilibrium. Solid line is related to the asymmetric solution, Eq. (32) ∗ [the asterisk means the equation taken from Table I of Ref. [9], line '1'], which is the leading approximation over a small leptodermic parameter, a/R /lessmuch 1. Dashed line presents the same but for the Wilets symmetric solution [1, 9] (line '2'). Parameters, for example: the effective radius R = 10 km, diffuseness of the NS crust a = 1 km, and temperature T = 0; see, e.g., Refs. [9, 27, 69]. The dots ES1 and ES2 show the ES at the effective radius R for the solid and dashed lines, respectively (from Ref. [9]).  \n<!-- image -->  \nedge 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).  \nThe coordinate system related to the effective surface is defined in such a way that one of  \nthe spatial coordinates ( ξ ) is the distance from a given point to the ES, and others ( η ) are tangent to the ES. This nonlinear coordinate system ξ, η 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 ρ 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 [38, 39]. Neglecting the other smallerorder 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 ξ ). 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.  \nWhen 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 [4, 40, 41]. 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 [43]. A relatively large change of the density ρ on a small distance a with respect to the curvature radius R takes place for a finite liquidmatter drop (nuclei, molecular water drops, neutron stars). We emphasize that the macroscopic ESA theory [9] allows to consider the deformed and superdeformed systems for which the same coordinate system ξ, η can be used locally near the ES. Inside of such dense systems, the density ρ is changed a little around constant saturation density ρ . Therefore, including also the second order expansion of the gravitational potential near ρ , one obtains essential effects of the surface capillary pressure.  \nThe accuracy of the ESA was checked in Ref. [5] for the case of static and dynamic nuclear physics by comparing the results with the existing nuclear theories like Hartree-Fock (HF) [44] and ETF [38, 39], based on the Skyrme forces [17, 39, 44-51], but for the simplest case without spin-orbit and asymmetry terms of the energy density functional. The direct varia-  \ntional principle for finding numerically the parameters of the tested particle number density functions in simple forms of the Woods-Saxonlike in Ref. [52] or their powers (Ref. [53]) were applied by using the realistic Skyrme energy functional [17, 50]; see also recently published variational ETF approach [54, 55]. The extension of the ES approach to the nuclear isotopic symmetry and spin-orbit interaction has been done in Refs. [6-8]. The Swiatecki derivative terms of the symmetry energy for heavy nuclei [56-63] were taken into account within the ESA in Ref. [8]. The discussions of the progress in nuclear physics and astrophysics within the relativistic local density approach, can be found in reviews Refs. [23, 64]; see also Refs. [65-67].  \nIn the present work we apply the effective surface approach for the neutron stars of Ref. [9], based on Refs. [5-8] 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.", 'II. TOLMAN-OPPENGEIMERVOLKOFF APPROACH': "Starting from the Einstein-Gilbert equations of the GRT for the spherically symmetric system, one has the Schwarzschild gravitational metric in a general form:  \nd s 2 = -e λ d r 2 -r 2 d θ 2 -r 2 sin 2 θ d φ 2 + c 2 e ν d t 2 , (1)  \nwhere c is the speed of light. Functions λ = λ ( r ) and ν = ν ( r ) are given by the differential equations (see Ref. [12] at zero cosmological constant Λ). First, outside of the system on big distances from masses, or for relatively small gravitational forces [12, 14], 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, θ, φ , and t , due to the GRT invariance, for the line element in this space, Eq. (1), one obtains the simple differential equations for λ and ν of the Schwarzschild metric for the region inside of the liquid-drop system (see, Refs. [12, 14]),  \n8 πG c 4 P = e -λ ( ν ' r + 1 r 2 ) -1 r 2 , 8 πG c 4 E = e -λ ( λ ' r -1 r 2 ) + 1 r 2 , dP dr = -1 2 ( E + P ) ν ' , (2)  \nwhere G is the gravitational constant. Primes above λ and ν mean the corresponding derivatives over the radial coordinate r , and E is the energy density. With the EoS, E = E ( ρ ) [or, P = P ( ρ )], one arrives at the complete system of equations for λ , ν , E (or density ρ ), and P .  \nUsing the substitution e -λ = f , one can solve the second equation in Eq. (2), as a linear firstorder differential equation over f ( r ), in terms of the quadratures. Then, the derivative ν ' 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 λ ( r ) and ν ( r ). Substituting the derivative ν ' into the third equation of Eq. (2), one finally obtains the famous TOV equations [10, 11, 28],  \nd P d r = -G ( E + P )( mc 2 +4 πr 3 P ) rc 2 ( rc 2 -2 Gm ) , d m d r = 4 πr 2 E c 2 . (3)  \nHere, m ( r ) is the gravitational mass interior to the radius r , and E is a given energy density. Boundary conditions for these equations are given by  \nm ( r = 0) = 0 , (d P/ d r ) r =0 = 0 . (4)  \nThe 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. [9].  \nAs shown in Ref. [12], taking approximately into account the step-function density, ρ = ρ for r < R inside of the system, and ρ = 0 for r > R outside of it, one can solve TOV equations (3) explicitly analytically (see Refs. [9, 12]). (For instance, ρ is the constant density ρ 00 using units G = c = 1, inside of the system in the notations of Tolman). Our leading order solutions for the density ρ are obtained in Ref. [9] by  \nusing 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 /lessmuch 1 (see Fig. 1). Taking also into account the first boundary condition in Eq. (4), the second equation in Eq. (3) can be integrated analytically,  \nm ( r ) = 4 π E 0 3 c 2 r 3 , (5)  \nwhere E 0 is a constant at r < R ,  \nE 0 = E ( ρ = ρ ) Θ ( ρ -ρ ) = A ( ρ )Θ( R -r ) , (6)  \nand ρ is the saturation approximation for the density ρ . We introduced here the Heaviside step function, Θ( x ), Θ( x ) = 1 for x ≥ 0 and 0 for x < 0. The function A ( ρ ) in Eq. (6) is given by Eq. (9) ∗ of Ref. [9] for A ( ρ ) at the saturation density ρ = ρ ; see Table I. The density ρ and energy density E 0 can be taken approximately as those found in Ref. [9]; see Eqs. (3) ∗ , (9) ∗ , (10) ∗ , and (11) ∗ , including the macroscopic gravitational field. Formally, one can consider the expression (9) ∗ for E 0 ( ρ ) inside of the system [or Eq. (16) ∗ ] as the macroscopic EoS additional to the TOV equations (3). 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  \nM = m ( R ) = 4 π E 0 3 c 2 R 3 , (7)  \nwhere R is the effective NS radius ( r 1 in the notations of Tolman's book, Ref. [12]).  \nFor the Schwarzschild metric in the interior region, one has  \nd s 2 = -d r 2 1 -r 2 /R 2 S -r 2 d θ 2 -r 2 sin 2 θ d φ 2 + c 2 [ A S -B S √ 1 -r 2 /R 2 S ] 2 d t 2 , (8)  \nwhere R S , A S , and B S are the specific integration constants, and r < R S . This metric, Eq. (8), is nonsingular at r = 0, in contrast to the original Schwarzschild metric [12, 14], 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. [12]). 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 S in Eq. (8) is expressed in terms of the energy density E inside of the system, E 0 , far away from the ES as  \nR S = ( 8 πG E 0 3 c 4 ) -1 / 2 . (9)  \nOther constants, A S and B S , as well as R S , are the following parameters of the transformed Schwarzschild metric, Eq. (8) (see the derivations in Ref. [12]),  \nA S = 3 2 √ 1 -R 2 R 2 S , B S = 1 2 . (10)  \nSubstituting Eq. (5) into the first TOV Eq. (3), one can see that the variables P and r are separated. Therefore, we analytically find the solution for r ≤ R ≤ R S :  \nP = E 0 3 3 ζ √ 1 -( r/R S ) 2 -1 1 -ζ √ 1 -( r/R S ) 2 , (11)  \nFIG. 2. Contour plots for the pressure P ( r ), Eq. (11), in units of the central value, P ( r = 0) as function of the radial coordinate r and the dimensionless parameter A S of the Schwarzschild metric, Eq. (8). The numbers in squares show the values of this pressure. White color presents regions where we have indetermination, infinity by infinity, with the finite limit 1 at r → 0. Red color shows negative values of the ratio P ( r ) /P ( r = 0) [a positive pressure P ( r )]. The value '0.00' displays the zero value in the horizontal and vertical coordinate lines on right of plots. The effective NS radius R = 10 km is the same as in Fig. 1 (from Ref. [9]).  \n<!-- image -->  \nwhere  \nζ = ∣ ∣ P 0 + E 0 / 3 P 0 + E 0 ∣ ∣ , (12)  \n∣ ∣ and P 0 ≡ P ( ρ ) = P ( r = 0). Here, R S is the radius parameter, Eq. (9), in the Schwarzschild metric, Eq. (8). Notice that the solution, Eq. (11), coincides, up to a constant, with that presented in Tolman's book at ζ = B S /A S = 1 / (2 A S ) (see Refs. [9, 12] and Eq. (9) for R S ). The cosmological constant also is assumed to be zero, and the units for the dimensionless pres-  \n8  \nFIG. 3. The pressure P ( r ), Eq. (11), in units of the central value, P ( r = 0), as function of the radial coordinate r for the Schwarzschild metric, Eq. (8), at two values of the parameter A S , A S = 1 . 15 (solid line) and A S ≈ 0 . 96 (dashed line). The latter corresponds to the two values for the NS masses, M = 1 . 4 M /circledot and M = 2 . 0 M /circledot ( R S ≈ 13 . 0 km and 15 . 6 km, respectively; see Refs. [70-73] and text). The effective NS radius R = 10 km is the same as in Figs. 1 and 2.  \n<!-- image -->  \nsure are used, E 0 ∝ R -2 S , Eq. (9). Thus, one obtains the relation of these parameters of the Schwarzschild metric, Eq. (8), to the initial condition P = P 0 at r = 0 through Eq. (12) for ζ . 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 S and B S . Finally, our result for the pressure, Eq. (11), coincides identically with that of Ref. [12].  \nFigure 2 shows the contour plots for the pressure P ( r ) [in units of P ( r = 0)], as function of r (km) and parameter A S , Eq. (10). It is more convenient to use the finite dimensionless values  \nof A S , 0 ≤ A S ≤ 1 . 5, instead of dimensional R S for the fixed effective NS radius R = 10 km, 0 ≤ r ≤ R . For a very rough estimate neglecting the relatively small gravitational mass defect, E 0 ∼ c 2 ρ (0) , where the central mass density ρ (0) ≈ M/ (4 πR 3 / 3) ∼ M /circledot / (4 πR 3 / 3), ( M ≈ M V ) and M /circledot = 1 . 989 ∗ 10 30 kg is the Sun mass, one finds significantly larger Schwarzschild radius parameter R S than the effective radius R , R S ≈ 15 . 6 km for M = 1 . 4 M /circledot , as example. By definition, R S is always larger than (or equal to) the effective NS radius R , R ≤ R S , while the gravitational radius r g should be smaller than R , r g = 2 GM/c 2 ≤ R [14]. These conditions are fulfilled in our calculations for the values R S ≈ 15 . 6 km, and r g ≈ 4 . 1 km in the case of the effective radius R = 10 km. With these evaluations, one has the constant A S ≈ 1 . 15. The maximum value A S ≈ 1 . 5 corresponds to zero pressure P ( r ).  \nFigure 3 shows the two cuts of the contour plots (Fig. 2) at the value, A S ≈ 1 . 15 ( R S ≈ 15 . 6 km, solid black line), and the value, A S = 0 . 96 ( R S = 13 . 0 km, dashed red line). The physical mass region, M/M /circledot ≈ 1 . 4 -2 . 0, in line with the experimental data for NS masses and radii (see, Refs. [70-73]), 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 S , for instance, taking into account the gravitational defect of mass [14] in the evaluations of the central energy density E 0 , one may have different values of A S presented in Fig. 2. Thus, as seen  \nfrom the solid and dashed lines of Fig. 3 and assumed in the derivations of the TOV equations [11, 12], the pressure P turns, indeed, to zero at the boundary of the NS system, r = R .", 'III. MODIFIED TOV APPROACH': 'We now derive the leptodermic corrections to the TOV equations due to the gradient terms of the energy density E , in line of Ref. [9]. The details of these derivations are shown in Appendix A.', 'A. General points': "The key quantity of our derivations is the energy density E ( ρ ) [5, 8, 9],  \nE ( ρ ) = A ( ρ ) + C ( ∇ ρ ) 2 . (13)  \nwhich determines the total system energy,  \nE = ∫ d V E [ ρ ( r )] . (14)  \nThe volume element, d V , for a spherical system is determined by the Schwarzschild metric, d V = exp( λ/ 2)d r , where d r = r 2 sin θ d θ d ϕ . In Eq. (13), A ( ρ ) is the volume non-gradient part, Eq. (9) ∗ ,  \nA = -b V ρ + ε ( ρ ) + G 4 ρ Φ( ρ ) , (15)  \nwhere b V is the non-gravitational contribution into the energy of separation of particle from the matter. In Eq. (15),  \nε ( ρ ) = K 18 ρ 2 ρ ( ρ -ρ ) 2 , (16)  \nwhere K is the non-gravitational part of the incompressibility modulus, e.g., due to the Skyrme nuclear interaction (for nuclear matter K ∼ 200 MeV), Φ( ρ ) is the statistically averaged (macroscopic) gravitational potential; see Eq. (8) ∗ and discussions around it in Ref. [9]. For simplicity, we will consider the gradient term of the interparticle interaction with a main constant C of Eq. (12) ∗ for a smooth coefficient B ( ρ ) in front of the gradient-density square term.  \nWe also use the condition for a minimum of the energy per particle, W , at a stable equilibrium,  \n( d W d ρ ) ρ = ρ = 0 , W = E ρ . (17)  \nTherefore, when we can neglect external forces, there is no linear terms in expansion of A , Eq. (15), over powers of the difference ρ -ρ near the saturation value ρ . Near the saturation value, ρ → ρ , one has W = A /ρ . Expanding the potential Φ( ρ ) near the saturation density ρ up to second order, and taking into account Eq. (17), one indeed finds that there is no linear term, namely, ∂ Φ /∂ρ = 0 at ρ = ρ , i.e., Φ = Φ 0 + (1 / 2)Φ 2 ( ρ -ρ ) 2 with some constants Φ 0 and Φ 2 . Combining different terms in Eq. (15), one writes  \nA = --b ( G ) V ρ + ε G ( ρ ) , (18)  \nwhere b ( G ) V = b V -m Φ 0 , m is the particle mass, Φ 0 = Φ( ρ ). The particle separation energy, b ( G ) V , and energy density component, ε G ( ρ ) (see Eq. (16) for ε ( ρ )) are modified both by the grav-  \nl field Φ( ρ ),  \nε G ( ρ ) = ε ( ρ ) + m 2 Φ 2 ( ρ -ρ ) 2 = K G 18 ρ 2 ρ ( ρ -ρ ) 2 . (19)  \nHere, K G is the total incompressibility modulus modified by the gravitational field Φ( ρ ) ( K G > 0) as  \nK G = K +9 mρ 2 Φ 2 . (20)  \nand Φ 2 is the second derivative of Φ( ρ ), Φ 2 = ∂ 2 Φ /∂ρ 2 , taken at the saturation density value ρ = ρ .  \nFor non-rotated one-component NS system we have to fix the particle number N ,  \nN = ∫ d V ρ ( r ) . (21)  \nTherefore, introducing the chemical potential µ as the Lagrange multiplier, one obtains from Eqs. (14) and (13) the equation for the equilibrium,  \nδ E δρ ≡ ∂ A ∂ρ -2 C ∆ ρ = µ . (22)  \nEquation (21) determines the chemical potential µ 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. [6-8]). They determine the two (neutron and proton) chemical potentials.  \nThe solution of the Lagrange equation (22) at leading order of the leptodermic expansion over the parameter a/R , ρ ( r ) = ρy (( r -R ) /a ), as function of the radial coordinate r , is shown in Fig. 1 in the main case (i), B = C = const , 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  \nvalue ρ ; see Refs. [9, 27, 69]. We compare the two dimensionless solutions for the particle density ρ ( r ) /ρ . One of them (solid line) is related to the vdW-Skyrme interaction constant C '1' for a dense liquid drop; see Eq. (32) ∗ for the universal dimensionless particle number density of a dense liquid drop, y = ρ/ρ , as function of x = ( r -R ) /a . Another limit case corresponds to the symmetric Wilets solution (ii) [1, 9] for zero constant C and no spin-orbit term but with taking into account the only kinetic gradient term in a gas system; see Eq. (12) ∗ for B . 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. As seen from Fig. 1 and expected, the range of the exponential decrease of ρ ( r ) becomes much smaller for the vdW-Skyrme liquid-drop case (i) than for the Wilets gas case (ii). But the surface effects on the local density ρ ( r ) in both cases are still very important to get a continuous analytical behavior through the ES. The smallness factor, proportional to a/R , with respect to the volume contribution, appears because of the integration over the normal-to-ES (radial in this example) variable x = ( r -R ) /a in this ES region. Thus, for the calculations of the surface corrections to the NS total energy E , Eq. (14), and particle number N , Eq. (21), it would seem possible to assume that they are not (very) important for so small relative diffuseness thickness a/R = 0 . 1. However, as shown in Ref. [9], the influence of the NS surface becomes really important for the local quantities as the particle  \ndensity ρ ( r ), Fig. 1, the EoS, and therefore, for the pressure P ( r ) calculating through Eqs. (13) and (17),  \nP = ρ 2 δ W δρ . (23)  \nMoreover, the gravitational mass defect due to the radial curvature of the nonlinear Schwarzschild metric in calculations of the NS mass leads to its dramatic nonmonotonic behavior of M = M ( R ); see Ref. [9].  \nFigure 4 shows the pressure P as function of the dimensionless density y [panel (a)] and radial variable r [panel (b)] through the density, y (( r -R ) /a ), at leading order over the parameter a/R ; see Eq. (60) ∗ for the full pressure ( P ), Eq. (61) ∗ for the volume ( P V ), and Eq. (62) ∗ for the surface ( P S ) components of the pressure (Table I and Ref. [9]). As seen from Fig. 4 (a), all pressures, P , P V and P S , take zero value at y = 0. The surface pressure, P S , takes zero value also asymptotically in the limit y → 1, in contrast to the volume part, P V , as well as the total pressure P . The full dot shows the ES at y = y 0 = 1 / 3 for the (i) case. As expected, the volume pressure P V is monotonically increasing function of the density y . In Fig. 4 (b), the pressure P ( ρ ( r )) (solid line) has a typical leptodermic behavior with a relatively short thickness of the order of a , where the pressure decreases sharply from the almost constant asymptote inside of the system at r ∼ > R -a to zero. Dashed red and dotted lines show the surface ( P S ) and volume ( P V ) components, respecively. Sum of these components is the total pressure P with  \n<!-- image -->  \nFIG. 4. Pressure P , Eq. (60) ∗ , in units of ρK G / 9 through the NS diffuse surface as function of the density variable y , (a), for the universal Equation of State, and the radial coordinate r , (b), through the particle number density, ρ = ρy (( r -R ) /a ). The density y ( x ) is given, e.g., by Eq. (32) ∗ (solid lines). Dashed and dotted lines show the surface (S) and volume (V) components, P S [Eq. (62) ∗ ] and P V [Eq. (61) ∗ ], respectively. The crust thickness a = 1 . 0 km and the effective radius parameter R = 10 km are the same as in Figs. 1 and 2. The full dots and arrows present the ES for y = y 0 , (a), and r = R , (b).  \n<!-- image -->  \nzero value at the ES because of the equilibrium; see the full point (ES) and arrow which present the effective surface and radius R , respectively.", 'B. Surface corrections to the TOV equations': "For the energy density E , Eq. (13), one can write  \nE = E 0 + E 1 , (24)  \nwhere E 0 is given by Eq. (6); see Eqs. (3) ∗ and (12) ∗ (Table I). The second component can be considered as a small leptodermic correction,  \nE 1 = C ( ∇ ρ ) 2 . (25)  \nIndeed, E 1 enters the solutions of Eq. (2) in terms of the integral over the radial coordinate r . Our derivations at leading order over leptodermic parameter a/R are similar to those for  \nthe equation for the density ρ and its solutions in Ref. [9].  \nStarting from the same equations (2), as presented in the previous section (see Ref. [12]), we split all their solutions into two components, the leading component and leptodermic corrections,  \nP = P 0 + P 1 , λ = λ 0 + λ 1 , ν = ν 0 + ν 1 , (26)  \nwhere P 0 is the zero order solution shown in Eq. (11). It is related to the constant saturation value for the density ρ inside of the system and zero outside, and corresponding energy density constant E 0 ; see Eq. (6). These quantities correspond to the zero order Schwarzschild metric parameter, λ 0 ( r ), which is the function of r (see Eq. (1) and Ref. [12]),  \ne -λ 0 = 1 -r 2 /R 2 S , r < R S . (27)  \nThe component ν 0 ( r ) for the Schwarzschild  \nmetric (1) can be found similarly as done in Ref. [12]. The leptodermic corrections, P 1 and ν 1 , Eq. (26), can be expressed in terms of the λ 1 .  \nAs shown in Appendix A, using approximately the property of the energy density component E 1 ∝ ( ∂ρ/∂r ) 2 , which is concentrated near the ES, for the correction λ 1 , at first order in the leptodermic expansion, one finally obtains from the second equation of Eq. (2),  \nλ 1 ( r ) = -8 πR S 1 -R 2 /R 2 S D 1 ( r ) , (28)  \nwhere  \nD 1 ≡ ∫ d r E 1 ( r ) ≈ a R RK G 9 ρ ∫ d ξ ( ∂ρ ∂ξ ) 2 . (29)  \nSubstituting the asymptote of the solutions, ρ ∼ exp( -ξ/a ), for larger distances ξ/a /greatermuch 1 (see Ref. [9]), one has an exponential disappearance of the integral D 1 , Eq. (29),  \nD 1 ∝ exp ( -2 ξ/a ) . (30)  \nFor the interior region, one finds also the exponentially decreasing asymptotes of D 1 , Eq. (29),  \nD 1 ≈ -a R ρRK G 9 Y ( y ) . (31)  \nThe function Y is given by  \nY ( y ) = ∫ d y √ y (1 -y ) = 2 3 y 3 / 2 -2 5 y 5 / 2 + C y , (32)  \nwhere y = ρ/ρ (see Ref. [9]), C y is an integration constant. To remove singularity at the origin, we shall assign the value -4 / 15 for the constant C y . Notice that in the derivations for λ 0 , Eq.(27)  \n[Eq. (96.7) in Tolman's book [12]] similar corrections were neglected by the same arguments. The solution λ 1 ( r ), Eq. (28), is zero outside of the system because of zero E 0 , Eq. (6), and P 0 , Eq. (11). Therefore, λ 1 ( r ) is concentrating near the NS surface, as it should be true (see also Appendix A).  \nSubstituting the derivative d ν/ d r from the first equation of Eq. (2) to its last equation, for a given energy density E , one obtains the modified TOV equation for the pressure P ,  \nd P d r = -E + P 2 r [( 8 πPr 2 +1 ) e λ -1 ] , (33)  \nwhere λ = λ 0 + λ 1 ; see Eqs. (27) and (28). This is valid approximately up to high, e.g. second order in the leptodermic expansion. Let us use now Eqs. (24) for E and (26) for P . Using also Eq. (A1) for P 0 at zero order of the leptodermic parameter a/R , one can derive the linearized equation for P 1 beyond the analytical solutions (27) and (11) of the standard TOV equations for λ 0 and P 0 , respectively.  \nFor the leptodermic correction P 1 to the pressure P , one has the same standard form as Eq. (A2) for λ 1 ,  \nd P 1 d r + βP 1 = β 1 . (34)  \nHere, β ( r ) and β 1 ( r ) are new functions of the radial variable r ,  \nβ = 8 πr 2 ( E 0 + P 0 ) + 8 πP 0 r 2 + ˜ r 2 2 r (1 -˜ r 2 ) , (35)  \nβ 1 = -8 πP 0 r 2 +˜ r 2 2 r (1 -˜ r 2 ) E 1 + ( E 0 + P 0 ) (8 πP 0 r 2 +1) 1 -˜ r 2 λ 1 , (36)  \nwhere ˜ r = r/R S , as in Appendix A. In these equations, P 0 ( r ) is the analytical solution of the TOV equations at zero order of the leptodermic parameter for the pressure, Eq. (11); see Refs. [9, 12]. Using the same method of the arbitrary constant variations, one obtains the solution in terms of the quadratures (see Appendix A),  \nP 1 = -C ( P ) 1 ( r )exp ( -∫ d rβ ( r ) ) , (37)  \nwhere  \nC ( P ) 1 ( r ) = ∫ d rβ 1 ( r )exp [∫ d rβ ( r ) ] . (38)", 'C. Leading solutions for the pressure': "Substituting the function β ( r ), Eq. (35), into the integral for the pressure correction P 1 , Eqs. (37) and (38), one finds  \nexp [∫ d rβ ( r ) ] = ( 1 -ζ √ 1 -˜ r 2 ) 2 √ ˜ r (1 -˜ r 2 ) , (39)  \nwhere ζ = 1 / (3 √ 1 -˜ R 2 ). Using then Eq. (36) for β 1 ( r ) in Eq. (38) for C ( P ) 1 , one obtains the sum of the two terms,  \nC ( P ) 1 ( r ) = -∫ d r [ 8 πP 0 r 2 +˜ r 2 2 r (1 -˜ r 2 ) E 1 + ( E 0 + P 0 ) ( 8 πP 0 r 2 +1 ) 1 -˜ r 2 λ 1 ] exp [∫ d rβ ( r ) ] . (40)  \nThe integrand of one of them is proportional to E 1 , Eq. (25), while the second one is proportional to λ 1 , Eq. (28). The first order corrections, E 1 and λ 1 , are determined by Eqs. (25) and (28), respectively [see also the approximate Eq. (A6)]. Both these functions exponentially decrease for r outside of a relatively small range  \na of the NS crust at small leptodermic parameter, a/R . Therefore, the first integral in Eq. (40) can be taken approximately at leading order, a/R , in the leptodermic expansion as for the derivation of λ 1 . Notice that the second integral in Eq. (40), with the integrand ∝ λ 1 ( r ), over radial coordinate r can be neglected as the second order of ( a/R ) 2 in the leptodermic expansion. Therefore, it does not contribute at first order over a/R . Indeed, this integral is reduced approximately to ∫ d r Y ( y ),  \nNotice that by the same reason, the λ 1 surface correction to the mass defect owing to the curved space can be neglected in calculations of the NS mass M (Ref. [9]) at the first order over a/R . Transferring the last integral over r to the variable ρ , d r = d ρ (d r/ d ρ ), and using Eqs. (24) ∗ , and (19) for the inverse derivative of ρ over r (d r = d ξ ), one finds  \n∫ d rλ 1 ( r ) ≈ -a R K G ρ 9 4 πR 2 √ 1 -R 2 /R 2 S × ∫ d r Y ( y ) . (41)  \n∫ d r Y ( y ) ≈ -a R R ∫ d y √ y (1 -y ) Y ( y ) . (42)  \nThe integral on the right hand side is finite because its integrand has no singularity in the limit y → 1 by the L'hospital theorem, Y / (1 -y ) → Y ' ( y ) / ( -1) → (1 -y ) → 0 . There is also no singularity at small y , ∫ d y/ √ y = 2 √ y .  \nTaking all smooth quantities in front of E 1 off the first integral in Eq. (40) over the radial variable r at r = R , one approximately finds  \nC ( P ) 1 ( r ) ≈ -˜ R 2 2 R ( 1 -˜ R 2 ) D 1 , (43)  \nwhere D 1 is given by Eq. (29) with the same asymptotic properties (30) and (31) outside of the narrow NS crust, and ˜ R = R/R S . With Eq. (39) for the integral exponent taken at r = R and the coefficient C ( P ) 1 , Eq. (43), one arrives at the final analytical approximation for the surface pressure correction P 1 ,  \nP 1 = -3 R 3 / 2 4 R 5 / 2 S √ 1 -R 2 /R 2 S D 1 ; (44)  \nsee Eq. (37). Here, R is the effective radius, R S is the Schwartzschild radius, Eq. (9), related to the initial value of the energy density, E 0 , and D 1 is given above by Eq. (29) with the internal (31) and external (30) asymptotes. Therefore, these leptodermic corrections for C ( P ) 1 , and hence, for P 1 , (as ∝ D 1 ), are exponentially decreasing functions of radial coordinate r , r →∞ outside [see Eq. (30)] and to zero inside of the system [Eqs. (31) and (32)]. In particular, inside of the system, Eq. (31) for D 1 , one arrives at  \nP 1 = a R ρK G 12 R 5 / 2 R 5 / 2 S √ 1 -R 2 /R 2 S Y ( y ) . (45)  \nThus, as expected, the leading leptodermic correction P 1 to the zero order pressure P 0 is peaked at the NS effective-surface range due to the function D 1 [Eqs. (29)-(30)]. Therefore, one finds the equation of state in the form:  \nP ≡ P ( ρ ) = P 0 ( y ) + P 1 ( y ) , (46)  \nwhere y = ρ/ρ , P 0 and P 1 are given approximately, at the first order leptodermic approach, by Eqs. (11) and (44), respectively. Substituting the solutions y = y ( r ) [ x = x ( y ), Eq. (31) ∗ with explicitly given Eq. (32) ∗ ] into Eq. (46) we  \nFIG. 5. Contour plots for the dimensionless pressure p 0 ( r, A S ), Eq. (47), as function of the radial coordinate r and the dimensionless parameter A S of the Schwarzschild metric, Eq. (8). Other notations and parameters are the same as in Fig. 2.  \n<!-- image -->  \nobtain r coordinate dependence of the pressure P ( ρ ( r )) at first order in a/R .  \nIt is also convenient to introduce the following dimensionless pressures as functions of r and A S variables, which are explicitly related to the constants c and G ,  \np 0 ≡ 8 πP 0 GR 2 c 4 = R 2 R 2 S 3 √ 1 -r 2 /R 2 S / 2 -A S A S -√ 1 -r 2 /R 2 S / 2 . (47)  \nThe explicit relation of the ratio R/R S to A S is given by R/R S = √ 1 -4 A 2 S / 9; see Eqs. (10) and (9). Then, for the correction P 1 , one similarly writes  \np 1 ≡ 8 πP 1 GR 2 c 4 = 3 κR 9 / 2 R 9 / 2 S √ 1 -R 2 /R 2 S Y ( y ) a R , (48)  \nwhere  \nκ ≡ ρK G 12 E 0 = -K +9 mρ 2 Φ 2 12( b V -m Φ 0 ) , (49)  \nFIG. 6. The same as in Fig. 5 but for the total pressure p ( r, A S ), Eq. (50), as function of the radial coordinate r and the dimensionless parameter A S of the Schwarzschild metric, Eq. (8). The numbers in squares show again the values of this pressure. White color presents regions where we have zero values. The effective NS radius R = 10 km is the same as in Figs. 1-4. The leptodermic parameter has a typical value a/R = 0 . 1. For example, the parameter κ ≡ ρK G / 12 E 0 = 10 for large gravitational forces ( κ /greatermuch 1).  \n<!-- image -->  \nK is a non-gravitative (e.g., nuclear) incompressibility; see also Eqs. (6) for E 0 and (18) for a smooth part of the energy density A ( ρ ) in terms of the gravitational potential derivatives of Φ. For the total dimensionless pressure p , at first order in the leptodermic parameter a/R , one finally obtains  \np ≈ p 0 + p 1 , (50)  \nwhere p 0 and p 1 are given by Eqs. (47) and (48).  \n16  \nFIG. 7. The dimensionless pressure p = P (8 πGR 2 /c 4 ), Eq. (50), as function of the radial coordinate r for the Schwarzschild metric in the form of Eq. (8) at the same value of the parameter A S , A S ≈ 1 . 15, and different values of the incompressibility parameter κ = 0 [black solid for only zero order pressure p 0 , Eq. (47)], 1 [red dashed line, Eq. (50)], 5 (green dotted line), and 10 (blue dashdotted line). The NS mass is the same as for the solid line in Fig. 3, M = 1 . 4 M /circledot . Other parameters are the same as in Fig. 6.  \n<!-- image -->  \np  \nFIG. 8. The same as in Fig. 7 but for the NS mass presented as the dashed line in Fig. 3, M = 2 . 0 M /circledot ( A S = 0 . 96). Other parameters are the same as in Figs. 6 and 7.  \n<!-- image -->", 'IV. DISCUSSIONS OF THE RESULTS': '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 S , relatively. The first order contribution p 1 of the pressure p is proportional to the dimensional incompressibility parameter κ , Eq. (49). As seen from the comparison of Figs. 6 ( κ = 10) and 5 ( κ = 0), the effect of the gravitational field on the pressure is significant near the ES ( r ∼ > R ) for this change because of the ES correction p 1 , Eq. (48) at a large dimensionless incompressibility κ , Eq. (49).  \nFigures 7 and 8 show the cuts of the contour plot of Fig. 6 at A S ≈ 1 . 15 corresponding approximately to the typical NS mass M = 1 . 4 M /circledot and A S ≈ 0 . 96 for M = 2 . 0 M /circledot , 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. [7072]; see also the comparison with the theoretical results in Ref. [73]. We specify different values of the incompressibility parameter κ from a weak (almost nuclear κ = 1) to a strong value ( κ = 10) of the gravitational potential. As mentioned above, the value κ = 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 κ = 10 in Fig. 8. For nuclear physics, one has  \nκ ≡ K/ (12 b V ) = 1 . 25 for K = 240 MeV, and b V = 16 MeV (or K = 225 MeV and b V = 15 MeV from Ref. [21]).  \nFigure 9 shows the EoS (a,c), p ( y ), for small [ κ = 1, (a)] and large [ κ = 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 . 4 M /circledot , red solid) and relatively large ( M = 2 . 0 M /circledot , 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 . 4 M /circledot ) and green long dashed ( M = 2 . 0 M /circledot ) lines, respectively. The surface effects become much more significant with increasing gravitational parameter κ . 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. [70-73]). 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  \nFIG. 9. The dimensionless pressures, p , Eq. (50), (thick solid black line for the surface component p 1 and thin red one for the total volume and surface quantity p ) as functions of the dimensionless particle number density, y , (density ρ in units of ρ ) for the solution (32) ∗ of leading order over the leptodermic parameter a/R , (a,c), and over the radial coordinate r , (b,d), for different dimensionless incompressibilites κ = 1, (a,b), and 10, (c,d), [see Eq. (49)], for the NS mass 1 . 4 M /circledot . Dashed (rare green p 1 and frequent blue p ) lines show these pressures for the NS mass 2 . 0 M /circledot . The arrows show the ES value y 0 = 1 / 3 in (a,c) and the effective radius R in (b,d).  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nr (km)  \nthat 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. [9].', 'V. CONCLUSIONS': '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 neu-  \ntron star was considered as a dense finite liquid (including amorphous solid) system at the equilibrium. The mean gravitational potential Φ( ρ ) was taken into account in the simplest form as expansion over powers of differences ρ -ρ , where ρ 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 /lessmuch 1) parti-  \number density ρ . The density ρ was derived as function of the normal-to-ES coordinate ξ of the orthogonal local nonlinear-coordinate system ξ, η through the effective surface. With its help, one finds the equation of state (EoS) for the pressure P = P ( ρ ). 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 [12] and numerical calculations [9]. 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.  \nAs 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 ξ, η for deformed and superdeformed shapes of the NS surface, one can  \napply our method to the pair of neutron star pulsars at their large angular momenta.', 'ACKNOWLEDGMENTS': 'The authors greatly acknowledge C.A. Chin, V.Z. Goldberg, S.N. Fedotkin, J. Holt, C.M. Ko, E.I. Koshchiy, J.B. Natowitz, A.I. Sanzhur, and G.V. Rogachev for many creative and useful discussions. A. Bonasera, T. Depastas, and A.G. Magner are partially supported by the US Department of Energy under Grant no. DE-FG0393ER-40773.', 'Appendix A: Derivations for λ 1': "Substituting Eqs. (24) and (26) into the second equation in Eq. (2), one has at zero order leptodermic parameter,  \n8 π E 0 = e -λ 0 ( 1 r d λ 0 d r -1 r 2 ) + 1 r 2 . (A1)  \nThen, using this zero order equation, from Eq. (2), one obtains in the linear order in λ 1 ( G = c = 1),  \n8 π E 1 = 1 r e -λ 0 d λ 1 d r -λ 1 ( 8 π E 0 -1 r 2 ) . (A2)  \nThis equation has the standard form with respect to the surface correction λ 1 as the first order linear differential equation,  \nd λ 1 d r + αλ 1 = α 1 , (A3)  \nwhere  \nα ( r ) = re λ 0 ( -8 π E 0 + 1 r 2 ) = 1 r (1 -˜ r 2 ) -8 π E 0 r 1 -˜ r 2 , (A4)  \nα 1 = 8 πr E 1 1 -˜ r 2 ≈ 8 πr C 1 -˜ r 2 ( d ρ d r ) 2 ≈ 8 πrε G 1 -˜ r 2 . (A5)  \nHere,  \nE 1 ≈ C ( d ρ d r ) 2 ≈ ε G ; (A6)  \nsee Eq. (27) for λ 0 and Eq. (26) for E 1 at the leading order of the leptodermic expansion for r < R S (˜ r < 1), as well as Eq. (24) ∗ for the last approximation. According to Eq. (19) for ε G , one has ε G ∝ y (1 -y ) 2 , i.e., ε G has a sharp maximum of the width a near the boundary ES (see Fig. 1).  \nUsing the method of an arbitrary constant variation, one obtains from Eq. (A2) the solution for λ 1 in the form:  \nλ 1 = C ( λ ) 1 ( r )exp ( -∫ d rα ( r ) ) , (A7)  \nwhere  \nC ( λ ) 1 ( r ) = ∫ d rα 1 ( r )exp (∫ d rα ( r ) ) . (A8)  \nFor the integrals in exponent of Eq. (A7) and (A8), one obtains  \n∫ d rα ( r )=ln ( ˜ r √ 1 -˜ r 2 ) +8 π E 0 R 2 S ln √ 1 -˜ r 2 , (A9)  \ni.e.,  \nexp (∫ d rα ( r ) ) = ˜ r ( 1 -˜ r 2 ) -1 / 2+4 π E 0 R 2 S . (A10)  \nTherefore, for C ( λ ) 1 ( r ), Eq. (A8), one has  \nC ( λ ) 1 ( r )=8 π ∫ r d r E 1 ( r ) ( 1 -˜ r 2 ) 4 π E 0 R 2 S -3 / 2 . (A11)  \nTaking into account that α 1 ∝ E 1 ∝ (d ρ/ d r ) 2 , Eq. (A5), at leading order of the leptodermic expansion, one can take the smooth coefficient in the integrand of Eq. (A11) off the integral at the ES, r = R . The reason is that the derivative square, ( ∂ρ/∂r ) 2 , is a bell-like function of r with a relatively small width of the order of the crust thickness a . We assume here that r is not close to the Schwarzschild radius R S , where the integrand in Eq. (A11) is singular. In line of Ref. [9], the remaining integration within the same precision can be taken analytically. Finally, for C ( λ ) 1 ( r ) one finds ( r < R S )  \nC ( λ ) 1 ( r ) = 8 πR ( 1 -˜ R 2 ) 4 π E 0 R 2 S -3 / 2 D 1 , (A12)  \nwhere D 1 is given by Eq. (29); see also Eq. (30) for the asymptote outside and Eq. (31) [in terms of Y ( y ), Eq. (32)] inside of the system. The effective radius R is presented in units of R S and c = G = 1, for shortness. From Eqs. (A12) and (A10), one arrives at the final expression (28) for λ 1 .  \nUsing Eqs. (24) for E and (26) for P and λ , and Eq. (A1) for P 0 at the zero order of leptodermic parameter a/R , one can derive the linearized equation (34) for P 1 beyond the standard TOV equations (A1) and (3) for P 0 . Equation (34) for the leptodermic pressure correction P 1 has the same standard form as Eq. (A2) for λ 1 . Using also the same method of the arbitrary constant variations, one obtains the solution in terms of the quadratures, Eqs. (37) for P 1 and (38) for C ( P ) 1 . Substituting the function β ( r ), Eq. (35), into the integrals in Eqs. (37)  \nand (38), one finds Eqs. (39) for ∫ d rβ ( r ) and (40) for C ( P ) 1 . As shown in section III, within the first order over a/R , one has to neglect λ 1 which gives the second order leptodermic correc-  \n- [1] L. Wilets, Phys.Rev. 101 , 1805 (1956); Rev. Mod. Phys. 30 , 542 (1958).\n- [2] V.M. Strutinsky, and A.S. Tyapin, Exp. Theor. Phys. (USSR) 18 , 664 (1964).\n- [3] A.S. Tyapin, Sov. Journ. Nucl. Phys. 11 , 401 (1970), 13 , 32 (1971), 14 , 50 (1972).\n- [4] V.M. Strutinsky, A.G. Magner, and M. Brack, Z. Phys. A 319 , 205 (1984).\n- [5] V.M. Strutinsky, A.G. Magner, and V. Yu. Denisov, Z. Phys., 322 , 149 (1985).\n- [6] A.G. Magner, A.I. Sanzhur, and A.M. Gzhebinsky, Int. J. Mod. Phys. E 92 , 064311 (2009).\n- [7] J.P. Blocki, A.G. Magner, P. Ring, and A.A. Vlasenko, Phys. Rev. C 87 , 044304 (2013).\n- [8] J.P. Blocki, A.G. Magner, and P. Ring, Phys. Rev. C 92 , 064311 (2015).\n- [9] A.G. Magner, A.I. Sanzhur, S.P. Maydanyuk, H. Zheng, A. Bonasera, T. Depastas, A.I. Levon, U.V. Grygoriev, arXiv:2403.01445 (nucl-th), accepted by the Int. Journ. Mod. Phys. E, 2024.\n- [10] R.C. Tolman, Phys. Rev. C 55 , 364 (1939).\n- [11] J.R. Oppenheimer and G.M. Volkoff, Phys. Rev. C 55 , 374 (1939).\n- [12] R.C. Tolman, Relativity, Thermodynamics, and Cosmology , (Dover Publications, New York, 1987; Oxford, the University Press, 1934,  \ntion to the pressure correction P 1 . Finally, with Eq. (43) for C ( P ) 1 , for the leptodermic correction P 1 , one arrives at Eq. (44).  \n1946, 1949, 1987).  \n- [13] J.S. Rowlinson and B. Widom, Molecular Theory of Capillarity (Clarendon Press, Oxford, 1982).\n- [14] L.D. Landau and E.M. Lifshitz, Theoretical Physics, v.2, The Classical Theory of Fields (Elsevier, Amsterdam,1951-1971; FIZMATGIZ, Moscow, 2003).\n- [15] G. Baym, H. Bethe, and C.J. Pethick, Nucl. Phys. A 175 (1971) 225.\n- [16] R.B. Wiringa, V. Fiks, A. Fabrocini, Phys. Rev. C 38 , 1010 (1988).\n- [17] E. Chabanat, P. Bonche, P.Haensel, J. Meyer, and R. Shaeffer, Nucl. Phys. A 627 , 710 (1997).\n- [18] I. Sagert, M. Hempel, C. Greiner, and J. Schaffner-Bielich, Eur, J. Phys. 27 , 577 (2006).\n- [19] Bao-An Li, Lie-Wen Chen, and Che Ming Ko, Phys. Rep., 464 , 113 (2008).\n- [20] A.Y. Potekhin, A. F. Fantina, N. Chamel, J.M. Pearson, and S. Goriely, Astronomy & Astrophysics, 560 , A48 (2013).\n- [21] G. Giuliani, H. Zheng, A. Bonasera, Prog. Part. Nucl. Phys. 76 , 116 (2014).\n- [22] Y. Lim and J.W. Holt, Eur.Phys. J. A 55 , 209 (2019).\n- [23] Boyang Sun, Saketh Bhattiprolu, and James  \nM. Lattimer, archiv:2311.00843v1, 2023.  \n- [24] R. Belverde, F. Cipolletta, C. Cherubini, S.M. de Carvalho, S. Filippi, R. Negreiros, J. P. Pereira, J.A. Rueda, and R. Ruffini, Physics and Astrophysics of Neutron Stars, The second ICRANet C'esar Lattes Meeting, AIP Conf. Proc. 1693,030001-1-030001-19; doi;10.1063/1.4937184 (AIP Publishing LLC 978-0-7354-5/ $ 30.00), p.030001-1 (2015).\n- [25] L.D. Landau and E.M. Lifshitz, Theoretical Physics, v.6, Fluid Mechanics (Elsevier, Amsterdam, 1959-1987).\n- [26] K. Schwarzschild, Berl. Ber., p. 424 (1916).\n- [27] P. Haensel, A.Y. Potehin, D.G. Yakovlev, Astrophysics and space science library, vol. 326, Neutron Stars 1. Equation of State and Structure (Springer, New York, 2007).\n- [28] J.M. Lattimer, J. Phys., Conf. Ser., 2536 , 012009 (2023).\n- [29] J.M. Lattimer and M. Prakash, The astrophysical Journal, 550 , 426 (2001).\n- [30] C.J. Horowitz and J. Piekarewicz, Phys. Rev. Lett., 86 , 5647 (2001).\n- [31] C.J. Horowitz and J. Piekarewicz, Phys. Rev. C, 64 , 062802(R) (2001).\n- [32] P.Haensel, Neutron Star Crusts, N. Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18, PL-00-716 Warszawa, Poland, 2000.\n- [33] Z. Arzoumanian, S. Bogdanov, J. Cordes, K. Gendreau, D. Lai, J. Lattimer, B. Link, A. Lommen, C. Miller, P. Ray, R. Rut-  \n- ledge, T. Strohmayer, C. Wilson-Hodge, and K. Wood, arXiv:0902.3264 , 2009; https://doi.org/10.48550/arXiv.0902.3264.\n- [34] A.F. Fantina, N. Chamel, J.M. Pearson, and S. Goriely, Astronomy & Astrophysics, 559 , A128 (2013).\n- [35] A.F. Fantina, and F. Gulminelli, J. Phys., Conf. Ser. 2586 , 012112 (2023); doi:10.1088/1742-6596/2586/1/012112 .\n- [36] H. Dinh, A.F. Fantina, and F. Gulminelli, Eur. Phys. J. A 59 , 292(2023); https://doi.org/10.1140/epja/s10050-02301199-x .\n- [37] H. Bethe, Annu. Rev. Nucl. Sci., 21 , 93 (1971).\n- [38] M. Brack and R. K. Bhaduri, Semiclassical Physics (Addison-Wesley, Reading MA) 1997; 2nd edition (Westview Press, Boulder) 2003.\n- [39] M. Brack, G. Guet, H.-B. H˚akansson, Phys. Rep. 123 , 275 (1985).\n- [40] A. Bohr and B. Mottelson, Nuclear structure (W. A. Benjamin, New York, 1975), Vol. II.\n- [41] A. G. Magner, V. M. Strutinsky, Z. Phys. A 324 , 633 (1985).\n- [42] J.S. Rowlinson, Journ. Stat. Phys. 20 , 197 (1979).\n- [43] A.N. Tikhonov, Systems of differential equations containing a small parameter in front of the derivatives , Mat. Sb. , 31 : 3 (1952) pp. 575-586 (In Russian).\n- [44] D. Vautherin, D. Brink, Phys. Rev. C 5 626 (1972).\n- [45] T.H.R. Skyrme, Philos. Mag., 1 8th ser., 1043  \n(1956).  \n- [46] R.C. Barrett, D.F. Jacson, Nuclear sizes and structure Oxford, Clarendon Press 1977.\n- [47] P. Ring, P. Schuck, The nuclear manybody problem, Berlin, Heidelberg, New York, Springer-Verlag 1980.\n- [48] J.P. Blaizot, Phys. Rep. 64 , 172 (1980).\n- [49] H. Krivine, J. Treiner, and O. Bohigas, Nucl. Phys. A, 336 , 155, (1980).\n- [50] E. Chabanat, P. Bonche, P.Haensel, J. Meyer, and R. Shaeffer, Nucl. Phys. A 635 , 231 (1998).\n- [51] B. Gramaticos, A. Voros, Ann.Phys., 123 , 359 (1979); 129 , 153 (1980.)\n- [52] M. Brack and R. K. Bhaduri, Semiclassical Physics. Frontiers in Physics No. 96, 2nd ed. (Westview Press, Boulder CO, USA, 2003).\n- [53] V.M. Kolomietz and A.I. Sanzhur, Eur. Phys. J. A 38, 345-354 (2008).\n- [54] V.M. Kolomietz, A.I. Sanzhur, and S. Shlomo, Phys. Rev. C 97, 064302 (2018).\n- [55] V.M. Kolomietz and S. Shlomo, Mean Field Theory, World Scientific, Singapore (2020).\n- [56] W.D. Myers and W.J. Swiatecki, Ann. Phys. (NY) 55 , 395 (1969); 84 , 186 (1974).\n- [57] W.D. Myers, W.J. Swiatecki, and C.S. Wang, Nucl. Phys. A 436 , 185 (1985).\n- [58] P. Danielewicz and J. Lee, Int. J. Mod. Phys. E 18 , 892 (2009).\n- [59] M. Centelles, X. Roca-Maza, X. Vi˜nas, and M. Warda, Phys. Rev. Lett., 102 , 122502 (2009).\n- [60] M. Warda, X. Vi˜nas, X. Roca-Maza, and  \n- M. Centelles, Phys. Rev. C 80 , 024316 (2009); 81 , 054309 (2010); 82 , 054314 (2010); arXiv:0906.0932 [nucl-th] (2009).\n- [61] X. Roca-Maza, M. Centelles, X. Vi˜nas, and M. Warda, Phys. Rev. Lett. 106 , 252501 (2011).\n- [62] X. Vi˜nas, M. Centelles, X. Roca-Maza, and M. Warda, Eur. Phys. J. A 50 , 27 (2014).\n- [63] J. Piekarewicz and M. Centelles, Phys. Rev. C 79 , 054311 (2009).\n- [64] T. Niksic, D. Vretenar, and P. Ring, Prog. Part. Nucl. Phys., 66 , 519 (2011).\n- [65] B.D. Reed, F.J. Fattoyev, C.J. Horowitz, and J. Piekarewicz, Phys. Rev. Lett. 126 , 172503 (2021).\n- [66] N. Chamel and P. Haensel, Liv. Rev. Relativ. 11 , 10 (2008).\n- [67] N. N. Shchechilin, N. Chamel, J. M. Pearson, Phys. Rev. C 108 , 025805 (2023).\n- [68] V.M. Strutinsky, N.Ya. Ljashtchenko, N.A. Popov, Nucl. Phys. 46 171 (1963).\n- [69] S.L. Shapiro, S.A. Teukolsky, Black holes, white dwarfs, and neutron stars: The physics of compact objects (Wiley-VCH Verlag GmbH& Co.KGaA, Weinheim, 2004).\n- [70] G. Raaijmaket et al., The Astrophysical Journ. Lett. 918 : L29 (2021), 1 (13 pp).\n- [71] C.D. Capano et al., Nature Astronomy 4 , (2020) 625.\n- [72] T.E. Riley et al., The Astrophysical Journ. Lett. 918 : L27 (2021), 1 (30 pp).\n- [73] O. Lorenco, M. Bhuyan, C.H. Lenzi, M. Dutra, C. Gonzalez-Boquera, M. Centelles, X. Vi˜nas,  \nPhys. Lett. B 803 , (2020) 135306.  \n- [74] H. Zheng and A. Bonasera, Phys. Rev. C 83 , 057602 (2011).  \n- [75] L.D. Landau and E.M. Lifshitz, Theoretical Physics, v.5, Statistical Physics (Elsevier, Amsterdam, 1951-1980)."}