Forms of the Symmetry Energy Relevant to Neutron Stars

Abstract

The symmetry energy is an invaluable tool for studying dense nuclear matter. Unfortunately, its definition is somewhat implicit, and therefore, phenomenological methods are necessary to describe experimental facts. This paper discusses the differences arising from the use of Taylor series expansion and Padé approximation to determine theoretically the symmetry energy and the possible consequences for neutron stars. For this purpose, a form of the nuclear matter equation of state that explicitly depends on the symmetry energy is used. The obtained results point out that the applied approximations lead to modifications of the equilibrium proton fractions and equation of state, especially in their high-density limit. However, this effect is small near the saturation density $n_0$.

1. Introduction

Symmetry considerations are often indispensable in obtaining solutions and explanations of many important problems in science and technology. Moreover, even approximate or broken symmetries acquire a special significance in physical sciences. This type of issue that includes analyses of the properties of dense and very dense states of isospin asymmetric nuclear matter is commonly considered as a serious problem of modern physics due to the importance of modeling quantum nuclear physics systems and compact astrophysical objects. The complexity of such many-body high-density systems practically prevents obtaining realistic analytical solutions, which necessitates that one use various effective methods of analysis. A huge amount of data from ground-based experiments should be compared with astrophysical data, obtained from observations of much denser objects than those available in laboratories on Earth, to get as much information as possible about nuclear matter in a wide range of densities. Therefore, a separate and highly nontrivial problem is the extrapolation of laboratory results to densities relevant to neutron star cores. The starting point for the analysis of a nuclear system is the equation of state (EoS) of isospin asymmetric nuclear matter. The EoS, defined as the binding energy per baryon, is a function of three variables: temperature, density, and isospin asymmetry, which reduce to density and isospin asymmetry when a model of a cold neutron star is considered. The understanding of the asymmetric nuclear matter EoS is currently insufficient to solve various problems in nuclear physics and astrophysics. Previous studies based on experimental data focused on symmetric nuclear matter, and the obtained results allowed fairly good determination of its properties near the saturation density [1]. The most valuable information about systems with different neutron and proton densities is encoded in the symmetry energy, which determines both the properties of asymmetrical nuclei and the matter of neutron stars. The standard approach to their analysis requires separating of the EoS into the symmetric and isospin-dependent parts. The commonly applied method uses the approximation of the function, which expresses this EoS, by a partial sum of terms of its Taylor series at the point of zero asymmetry. This is also a way to define the symmetry energy. Many theoretical approaches to the issue of the symmetry energy dependence on density have been developed so far. Their diversity comes from the lack of restrictive constraints on the form of the EoS resulting from ground-based experiments, as well as from astrophysical observations. However, even using very sophisticated models, the correct reconstruction of neutron star parameters and the explanation of experimental data obtained for finite nuclei are still debatable. Different types of isospin sensitive observables, which include among others: atomic masses, neutron skins of heavy nuclei, heavy-ion multifragmentation reactions, Pygmy dipole resonances, measurements of the electric dipole polarizability of ²⁰⁸Pb [2], frequency of isovector giant dipole resonance in ²⁰⁸Pb, isospin diffusion in heavy ion reactions, excitation energies of isobaric analog states, isoscaling of fragments from intermediate energy heavy-ion collisions used in ground-based experiments to obtain data that allow constraining the symmetry energy below and in the vicinity of the saturation point $n_0$. Analysis of these data led to relatively consistent results for the form of the symmetry energy for sub-saturation densities [3,4]. Heavy-ion collisions (HICs) represent a class of experiments that allow one to study dense nuclear matter under conditions changing in a controlled manner. Although the experimental creation and then the analysis of nuclear matter in a wide temperature, density, and neutron-proton asymmetry range become feasible, the results strongly depend on the achievable composition of projectile and target nuclei and the energy of beams. Recent investigations of the asymmetric nuclear matter mostly focused on the symmetry energy. Nevertheless, the symmetry energy is still unconstrained for densities higher then the saturation density. This applies to both the experimental and theoretical aspects of this issue. The direct measurement of the symmetry energy is not possible, and therefore, many theoretical assumptions are necessary for the interpretation of the experimental data. The basic problem that influences the interpretation of the existing data is the fact that HICs are highly non-equilibrium processes, and their analysis necessarily involves a transport model. Another problem is the identification of observables sensitive to the symmetry energy in the presence of many uncertainties including those introduced by transport models. Analyses based on the theory of transport suggest using isospin ratio observables for extracting information about the symmetry energy [5,6]. The measurements of the observables in question are conditional on the energies achievable during experiments. The choice of observables suitable for obtaining information on the properties of asymmetric nuclear matter in the high-density range is connected with the phenomena of collective flows and meson production. The upper limit of the density that can be achieved in central heavy-ion collisions at beam energies of several hundred MeV per nucleon in the range $2÷3n_0$ [7]. A description of collisions during subsequent phases of the system evolution is carried out based on the transport model, which typically uses the power-law function that represents the symmetry energy dependence on density. The potential part of the symmetry energy has the form:

$$E_{2,sym}^{pot}=C{\left(\frac{n_b}{n_0}\right)}^{\gamma },$$

(1)

with the exponent $\gamma$ and the experiment-dependent parameter C, which determines the value of the potential energy at saturation density $n_0$. In the ultrarelativistic quantum molecular dynamics (UrQMD) transport model [8,9], it is taken at the level of $20÷22$ MeV. The kinetic part of the symmetry energy is given in the following form:

$$E_{2,sym}^{kin}\left(n_b\right)=\frac{1}{3}{ϵ}_F{\left(\frac{n_b}{n_0}\right)}^{2/3},$$

(2)

where ${ϵ}_F$ is the Fermi energy (${ϵ}_F$∼28 MeV). Calculations done within the UrQMDmodel that used the presented power-law parametrization of the symmetry energy allow for improved interpretation of the preliminary data from the FOPI-LANDexperiment [10,11]. The obtained results allow limiting the possible functional dependence of the symmetry energy on density to the value of the coefficient $\gamma =0.9\pm 0.4$. The new ASY-EOSexperiment that was carried out at the energy of 400 MeV/nucleon for ${}^{197}$Au+${}^{197}$Au, ${}^{96}$Zr+${}^{96}$Zr, and ${}^{96}$Ru+${}^{96}$Ru [7] improved the precision of measurements of the experimental flow parameters, and the statistics allowed restricting the value of the exponent $\gamma$ to $\gamma =0.72\pm 0.19$ [7]. The power-law form of the symmetry energy is a concise, but convenient representation of the EoS of asymmetric nuclear matter. Exemplification of this fact was presented in paper [12]. The parametrization of the symmetry energy in the form $E_{2,sym}\left(n_b\right)=C{\left(\frac{n_b}{n_0}\right)}^{\gamma }$ together with the definition of the symmetry energy slope L, after considering that symmetric nuclear matter has a well determined value of energy at the saturation density and imposing the requirement of vanishing pressure of symmetric nuclear matter at saturation, allow for the determination of the exponent $\gamma$. It is given by a combination of parameters that were used for the general description of the function representing the form of the EoS of pure neutron matter (PNM) [12].

The analysis of the symmetry energy in a much more general context indicates the need to study it in a broader range of asymmetry and density. Such physical conditions are realized in the interiors of neutron stars. Thus, astrophysical observations, as well as theoretical modeling of neutron star matter are an additional source of information about the symmetry energy and allow verification and clarification of existing experimental and theoretical constraints. Phenomena and astrophysical processes whose correct understanding requires knowledge of the symmetry energy are mainly the formation and evolution of neutron stars. The formation of a neutron star strongly depends on the details of the EoS of dense, hot nuclear matter. The essential problem is to understand under what thermodynamic conditions do the stages of neutron star evolution occur. Different effects of the symmetry energy, which is the main part of the corresponding EoS, are crucial, but their impact on the individual stages of a neutron star evolution varies and depends on both the density and the value of isospin asymmetry. These changing physical conditions allow one to trace the process of a neutron star formation during different phases starting with the core collapse, the post-bounce behavior of a supernova core, the post-bounce shock evolution, the early phases of proto-neutron star thermal evolution through the neutrino emission, which is related to neutrino opacity, or the late phase of proto-neutron star cooling. A separate issue is the analysis of the observational properties of neutron stars carried out in order to identify those that depend on the form of the symmetry energy [13]. An attempt to answer the question about the dependence of symmetry energy on density based on information obtained from observations of neutron stars is subject to considerable uncertainty. This uncertainty increases significantly due to the possibility of the emergence of additional particles in the very inner part of a neutron star core. Equilibrium conditions, which settle the chemical composition of neutron star matter, suggest that hyperons are expected to emerge and remain stable at sufficiently high densities, leading to the formation of the internal hyperon core. From the point of view of symmetry energy, it is important to analyze such properties of the neutron star as its radius and chemical composition, which result from the dependence of symmetry energy on density. Thus, analysis of the concentration of hyperons and their onset points is particularly promising in an attempt to answer the question about the form of the symmetry energy. The ability to activate various hyperonic neutrino processes can alter the efficiency of neutron star cooling. Significant theoretical efforts have been made to build models of strangeness-rich neutron star matter that meet the limitations of both ground-based experiments and astrophysical observations. The development of models that allow the occurrence of direct URCA reactions in neutron stars with hyperons goes well with the limitation of the observation of SAXJ1808.4-3658 [14,15,16,17]. The reported results indicated the direct URCA process involving hyperons as a probable dominant process. A possible very accurate scenario of neutron star cooling was considered in the paper [18]; contributions to the neutrino emissivity that come from the processes involving hyperons were taken into account. Observational determination of neutron star masses and radii would give the opportunity to obtain constraints on the EoS of their matter. However, only neutron star mass is determined directly from observations and gives real bounds on the form of the EoSs by eliminating those that give a maximum mass less than the largest mass determined observationally. Thus, the existence of a class of massive neutron stars with precisely determined masses is an extremely effective factor constraining the EoS. Previously obtained results include MSPJ1614-2230 with the updated value of mass $1.928\pm 0.017M_{\odot }$ [19,20,21] and pulsar J0348 + 0432 with the mass $2.01\pm 0.04M_{\odot }$. These two mass measurements verify models that take into account additional components of matter. Recent astronomical observations have shifted the maximum mass limit to ${2.14}_{-0.09}^{+0.10}{M}_{\odot }$ [22], which further complicates the theoretical description of dense nuclear matter with additional degrees of freedom (e.g., hyperons or $\Delta$ isobars). Their appearance in the core of a neutron star modifies the mass of the star, controls its evolution, softens the EoS, and at the same moment, reduces the value of the maximum mass possible to obtain in a given model. This reveals an inconsistency between models and the astrophysical observations setting a discriminating limit to the theoretical description of nuclear forces. Analysis of this problem indicates a special relationship between the properties of neutron star matter with $\Delta$ isobars and hyperons and the form of the symmetry energy. Their threshold densities have been shown to correlate with the symmetry energy slope parameter L. Theoretical models indicate that for experimentally acceptable values of L, $\Delta$ isobars appear in the density range of $2÷3n_0$ [23]. This fact has a fundamental impact not only on the chemical composition of stars, but also causes a significant softening of EoSs. The problem of the formation of $\Delta$ isobars in the matter of neutron stars due to its similarity to the problem associated with the appearance of hyperons requires multi-faceted studies that take into account the various types of existing experimental limitations. However, it is the radius of a neutron star that is considered as the most promising probe of the symmetry energy. The problem is that it can be estimated only indirectly. The assessment of the value of radius is based among others on observations of low-mass X-ray binaries [24]. The available data come from thermonuclear X-ray bursts and spectra. However, recently, the NICER mission succeeded in measuring simultaneously both the mass and radius of PSR J0030+0451 with a radius $R={12.7}_{-1.2}^{+1.1}$ km and ${13.0}_{-1.1}^{+1.2}$ km at the $68\%$ confidence level [25,26,27]. The importance of the symmetry energy can also be perceived in the context of its capabilities to probe the internal structure of a neutron star. Theoretical models of neutron stars can provide hints of the internal structure of neutron stars, yet their verification is out of the reach of present experiments. A related very important issue is the stratification of matter inside a neutron star. The focus is on crust-core transition density and pressure that directly determine the properties of neutron star crust including its thickness, fractional mass, and moment of inertia. Another phenomenon that enables us to study the effect of the symmetry energy on the structure of a neutron star is connected with the determination of the r- mode instability window [28]. A very important source of information on EoS of dense asymmetric nuclear matter is a precise measurement of neutron star tidal polarizability being understood as the quadrupole deformation of a neutron star caused by a quadrupolar gravitational field from the companion star in a binary system. It is defined as:

$$\mathsf{\Lambda }=\frac{2}{3}{k}_2{\left(\frac{c^2}{G}\frac{R}{M}\right)}^5,$$

(3)

where $k_2$ is the tidal Love number, c the velocity of light, G the gravitational constant, and R and M the global neutron star parameters. The factor $\mathsf{\Lambda }$ is linked to global neutron star parameters and through this can impose constraints on the EoS. The initial analysis of data from GW170817 predicted $\tilde{\mathsf{\Lambda }}\le 800$, where $\tilde{\mathsf{\Lambda }}$ denotes the combined dimensionless tidal deformability of the binary neutron star system. The subsequently modified results constrained this value to $70\le \tilde{\mathsf{\Lambda }}\le 700$. Results obtained within the chiral effective field theory can be quoted here [29,30]. The obtained EoS based on the chiral interactions provides constraints on the symmetry slope L (43 < L < 52 MeV). It also offers the opportunity of constraining the mass-radius relation in the range of $\tilde{\mathsf{\Lambda }}$ consistent with the results obtained from the GW170817 event.

The analysis carried out in this paper is based on a field-theoretical Walecka-type model that gives the opportunity to construct classes of solutions representing different forms of the symmetry energy. The linking between the symmetry energy and neutron star properties can be studied changing the symmetry energy slope parameter L. Correct adjustment of parameters in the isovector sector of the model enables examining the dependence of the symmetry energy on density for a wide density range and at the same time with an acceptable modeling of nuclei. The influence of the symmetry energy on neutron star properties strongly depends on the form of symmetry energy itself due to the necessity of approximation of the function that represents the EoS of asymmetric nuclear matter. Here, a key issue is the approximation method being used. There are mathematical premises allowing one to arrive at the conclusion that Padé approximants often perform much better than the Taylor expansion. The use of Padé approximation in the hope of obtaining more faithful formulas modifies the basic characteristics of neutron stars. The paper is organized as follows. In the next section, a general introduction to the symmetry energy with stress on its dependence on the approximation method is given. Then, isospin sensitive properties of neutron stars are discussed. The results presented in Section 6 are based on the formalism that was introduced in Section 5.

2. Symmetry Energy: The Taylor Series

In order to derive typical properties of nuclear matter, preliminary analysis consisting of adopting the appropriate representation of the function describing the EoS $E(n_b,{\delta }_a)$ of the system must be performed. The EoS $E(n_b,{\delta }_a)$ depends on the baryon density $n_b$ and the isospin asymmetry parameter ${\delta }_a=(n_n-n_p)/(n_n+n_p)$, where $n_n$ and $n_p$ are neutron and proton number densities. An important result is obtained from the Taylor expansion of the function $E(n_b,{\delta }_a)$ at the point ${\delta }_a=0$. This leads to the typical decoupling of the symmetric matter from the isospin dependent part:

$$E(n_b,{\delta }_a)=E(n_b,{\delta }_a=0)+\frac{1}{2}\frac{{\partial }^{2}E(n_b,0)}{\partial {\delta }_a^2}{\delta }_a^2+\frac{1}{24}\frac{{\partial }^{4}E(n_b,0)}{\partial {\delta }_a^4}{\delta }_a^4+\dots .$$

(4)

The symmetric part of the function $E(n_b,{\delta }_a)$ is given by $E(n_b,{\delta }_a=0)$, whereas the asymmetric one is denoted by $E_{2,sym}\left(n_b\right)$ in the case when it is limited to the second-order term of the expansion Equation (4). Subsequently, the Taylor series expansion of the symmetric and asymmetric parts at equilibrium density $n_0$ results in the following formulas:

$$E(n_b,{\delta }_a=0)=E\left(n_0\right)+\frac{1}{2}\frac{{\partial }^{2}E(n_0,0)}{\partial n_b^2}{(n_b-n_0)}^2+\dots ,$$

(5)

$$E_{2,sym}^{full}\left(n_b\right)=E_{2,sym}\left(n_0\right)+\frac{\partial E_{2,sym}\left(n_0\right)}{\partial n_b}(n_b-n_0)+\frac{1}{2}\frac{{\partial }^2{E}_{2,sym}\left(n_0\right)}{\partial n_b^2}{(n_b-n_0)}^2+\dots .$$

(6)

In order to find different types of correlations [31] between relevant physical parameters that can be probed in experiments (or observations), several factors that characterize the EoS can be defined:

• the binding energy at saturation density $n_0$: $E(n_0,0)$
• the symmetry energy $E_{2,sym}\left(n_b\right)$ with the symmetry energy coefficient $E_{2,sym}\left(n_0\right)$
• the slope of the symmetry energy $E_{2,sym}$: $L=3n_0\frac{d{E}_{2,sym}}{d{n}_b}\left(n_0\right)$
• the incompressibility of symmetric nuclear matter $E(n_b,0)$: $K_0=9n_0^2\frac{d^{2}E(n_b,0)}{d{n}_b^2}\left(n_0\right)$
• the curvature of the symmetry energy: $K_{sym}=9n_0^2\frac{d^2{E}_{2,sym}\left(n_b\right)}{d{n}_b^2}\left(n_0\right)$.

One may wonder why the term “symmetry energy” is used. The basic assumption is the $SU\left(2\right)$ isospin symmetry invariance of strong interactions. If the squared isospin $I^2$ and its third component $I_3$ of $SU\left(2\right)$ commute with the interaction Hamiltonian, then the $I_3=\frac{1}{2}(n_p-n_n)$ and I parameterize the eigenstates. The total energy $E(n_p,n_n)$ is a combination of Clebsch–Gordan coefficients and the reduced matrix elements of the Hamiltonian, reflecting the isospin symmetry. This makes it possible, for a given energy density functional that reflects the nucleon binding energy, to calculate functions $E_{2,sym}\left(n_b\right)$ that reproduce $E_{2,sym}\left(n_0\right)$ at, say, densities close to the saturation density $n_0$. Of course, the definitions of experimentally relevant variables or parameters depend on the approach being used. For example, the inclusion of the fourth-order term $E_{4,sym}\left(n_0\right)$ instead of $E_{2,sym}\left(n_0\right)$ alone is possible [32].

3. Symmetry Energy: The Padé Approximation

The computational problems in the analysis of the symmetry energy are involved in its implicit-like definition. The commonly accepted one refers to the Taylor series of the energy density $E(n_b,{\delta }_a)$ of nuclear matter expanded at the parameter ${\delta }_a=0$. According to Taylor’s theorem, if we are in its disc of convergence, the series Equation (4) gives a useful representation of the function $E(n_b,{\delta }_a)$ in the vicinity of the point $(n_0,0)$ with a quantitative estimate on the error of approximation. However, this formalism can be questioned in the case of high density and isospin asymmetry [33,34]. Such conditions prevail in the interiors of neutron stars. The use of a finite number of Taylor terms does not guarantee stable nor systematically improvable approximation to the exact form of the energy density, i.e., it could happen that the finite partial sum explodes with increasing order in the Taylor series extension [35]. At the same time, Padé approximants represents more faithfully the initial function than its truncated Taylor series. Such behavior was noticed in the perturbative calculations of the ground-state energies of closed-shell nuclei models [35]. Since the late 19th Century, the Padé approximation [36] has become an independent research field of analysis and an important chapter in the theory of rational approximation of analytic functions. Padé approximants can be perceived as formal Gaussian quadrature methods [37]. It is possible to use Taylor series’ coefficients to build Padé approximants and thus study the global properties of the original analytic function (e.g., analytical continuation, character, and distribution of singularities). Padé approximants can be used in a larger region of the function domain (where the Taylor series might not be converged) [35,38]. For a given formal power series $c\left(z\right),$ $c\left(z\right)\equiv {\sum }_{i=0}^{\infty }{c}_i{z}^i$ representing the function $f\left(z\right)$ and for any pair of non-negative integers $(l,m)$, the ${[l,m]}^{\mathrm{th}}$ Padé approximation of $f\left(z\right)$ is defined as the rational function $P_m^l\equiv \frac{A_l}{B_m}=\frac{{\sum }_{i=0}^l{a}_i{z}^i}{{\sum }_{i=0}^m{b}_i{z}^i}$, provided all terms in the formal power series $\left({\sum }_{i=0}^m{b}_i{z}^i\right)\left({\sum }_{i=0}^{\infty }{c}_i{z}^i\right)-{\sum }_{i=0}^l{a}_i{z}^i$ vanish up to the term $z^{l+m}$. Under some technical assumptions, the Padé approximation is unique [36], and $b_0$ can be chosen as equal to one. It is essential to control the truncation error, $R_{l,m}\left(z\right)=f\left(z\right)-P_m^l$, which arises when the function $f\left(z\right)$ is replaced by its approximation. The correct behavior $R_{m_k,m_k}\left(z\right)\to 0$ is guaranteed by the Padé conjecture [39,40]. It assumes the convergence of the infinite diagonal $(l_k=m_k)$ subsequence $\left\{P_{m_k}^{m_k}\right\}$ of the full sequence ${\left\{P_m^m\right\}}_{m\in \mathrm{N}}$ of Padé approximants for increasing sequences $\left(l_k\right)$ and $\left(m_k\right)$. Then, if $l_k+m_k\to \infty$ as $k\to \infty$, $k=1,2,\dots$, and for $(l_k=m_k)$, the limit $R_{m_k,m_k}\left(z\right)\to 0$ follows from the Padé conjecture in the formulation given by Baker, Gammel, and Wills [39]. The Padé conjecture laid the foundation for various scientific applications and numerical experiments [41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56]. The validity of this conjecture is also assumed in this paper. Finally, various generalizations of the classical Padé approximation (rational approximation to a function with close contact at one point) are also considered, e.g., among others: the Padé–Hermite approximation (simultaneous rational approximation to a vector of functions), the multiple-point Padé approximation (general interpolation processes by means of rational functions with free poles) [57,58].

In order to compare the effectiveness of Taylor and Padé methods, several algorithms have been developed [59]. An analysis of the performance of both approximations for the isospin-dependent part of the EoS is one of the purposes of this paper. As the Taylor series Equation (4), which was obtained by expanding the function $E(n_b,{\delta }_a)$ with respect to the asymmetry parameter $z\equiv {\delta }_a$, has vanishing odd-numbered coefficients, it is seen that $P_2^2$ is the lowest nontrivial Padé approximant. Vanishing of the coefficients $c_1$ and $c_3$, ($c_1=c_3=0$), results in the following formula:

$$\begin{array}{ccc}{P}_2^2(n_b,{\delta }_a)\hfill & =\hfill & \frac{c_0{c}_2+(c_2^2-c_0{c}_4){\delta }_a^2}{c_2-c_4{\delta }_a^2}\hfill \\ \hfill & =\hfill & c_0+\frac{c_2^2}{c_2-c_4{\delta }_a^2}{\delta }_a^2,\hfill \end{array}$$

(7)

where $c_i$, $i=0,2,4$ are the coefficients of the Taylor series, and both are functions of baryon number density $n_b$. Therefore, the following relations hold: $c_0\left(n_b\right)\equiv E(n_b,0)$, $c_2\left(n_b\right)\equiv E_{2,sym}\left(n_b\right)$, and $c_4\left(n_b\right)\equiv E_{4,sym}\left(n_b\right)$. This leads to the result:

$$E_{asym}^P(n_b,{\delta }_a)=\frac{E_{2,sym}\left(n_b\right)}{1-\frac{E_{4,sym}\left(n_b\right)}{E_{2,sym}\left(n_b\right)}{\delta }_a^2}.$$

(8)

Thus, the asymmetry energy $E_{asym}^P(n_b,{\delta }_a)$ obtained from the second summand in Equation (7) is a combination of two factors, the symmetry energy $E_{2,sym}\left(n_b\right)$ and the asymmetry factor:

$${\Delta }_a^P(n_b,{\delta }_a)\equiv {\left(1-\frac{E_{4,sym}\left(n_b\right)}{E_{2,sym}\left(n_b\right)}{\delta }_a^2\right)}^{-1}.$$

(9)

4. Isospin-Sensitive Properties of Neutron Stars

Neutron star matter has one of the largest isospin asymmetry values; therefore, the symmetry energy, which enters the EoS of asymmetric nuclear matter, is considered as the decisive factor controlling the process of evolution, properties, and structure of a neutron star. The analysis of the significance of the symmetry energy for neutron stars should also take into account the approximation method used for its calculation. Calculations performed for early phases of neutron star evolution, which include the gravitational collapse of the core of massive stars and, after that, the proto-neutron star cooling, confirm the significance of the symmetry energy [60]. The isospin asymmetry changes in a wide range during the formation of a neutron star, so not at every stage of evolution is the influence of symmetry energy the same. Numerical simulations of the gravitational collapse of the massive star core indicate a small influence of the symmetry energy in the early phase of evolution involving core bounce and a subsequent period of time lasting about 200 ms. This fact indicates a low isospin asymmetry value. During successive phases of a proto-neutron star evolution, consisting of cooling through the emission of neutrinos, the increasing role of the symmetry energy is evident. Theoretical models of the symmetry energy that differ in the value of the symmetry energy slope L lead to changes in both the cooling time scale and the brightness of neutrinos [60].

A model of the internal structure of a neutron star must accurately reconstruct the generally accepted stratification of its matter, and the importance of the symmetry energy with regard to this problem is to determine the transition density between neutron star crust and core. The appropriate location of the transition density allows determining the physical properties of the crust such as its thickness, fractional mass, and moment of inertia [13]. The approach that allows one to determine the transition density relies on the thermodynamical stability condition of the n-p-e matter against the growth of small density fluctuations and can be formulated by the requirement of the positiveness of the expression [13,61]:

$$V_{therm}(n_b,Y_p)=\left(2n_b\frac{\partial E(n_b,Y_p)}{\partial n_b}\right)+n_b^2\frac{{\partial }^{2}E(n_b,Y_p)}{\partial n_b^2}-\left(n_b\frac{{\partial }^{2}E(n_b,Y_p)}{\partial n_b\partial Y_p}\right)\left(\frac{{\partial }^{2}E(n_b,Y_p)}{\partial Y_p^2}\right),$$

(10)

where $E(n_b,Y_p)$ is the binding energy of asymmetric nuclear matter and ${\delta }_a=1-2Y_p$. The solution of the equation $V_{therm}=0$ determines the value of the core-crust transition density. The symmetry energy directly enters Equation (10); thus, an important part of this study is to determine the extent to which the approximation method influences the crust-core transition density and pressure. Different results are obtained when parabolic or Padé approximations are applied. Reformulation of the EoS of dense asymmetric nuclear matter being in $\beta$-equilibrium, so that it explicitly depends on the symmetry energy, requires determination of the equilibrium proton fraction $Y_p^{eq}$, which can be calculated on the basis of the following equation [62]:

$$\hslash c{\left(3\pi n_b^2{Y}_p\right)}^{1/3}=4E_{2,sym}\left(n_b\right)(1-2Y_p).$$

(11)

The solution of this equation shows that the equilibrium concentration of protons $Y_p^{eq}$ also explicitly depends on the form of the symmetry energy $E_{sym}\left(n_b\right)$.

5. The Formalism

5.1. Nuclear Matter

The core of a neutron star contains matter that constitutes an extremely complex many-body system. One of the approaches used to describe this type of matter theoretically is based on relativistic mean-field theory (RMF). This is an effective theory describing the interaction between nucleons through the exchange of different kinds of mesons. This model includes the isospin doublet nucleon field ${\psi }_N$ and meson fields: scalar $\sigma$ and vector $\omega$ and $\rho$. The dynamics of the system is governed by the following Lagrangian density function:

$$\begin{array}{ccc}\mathcal{L}\hfill & =\hfill & {\overline{\psi }}_N\left(i{\gamma }^{\mu }{D}_{\mu }-m_N\right){\psi }_N+\frac{1}{2}{\partial }_{\mu }\sigma {\partial }^{\mu }\sigma -\frac{1}{2}{m}_{\sigma }^2{\sigma }^2+\frac{1}{2}{m}_{\omega }^2\left({\omega }_{\mu }{\omega }^{\mu }\right)+\hfill \\ \hfill & +\hfill & \frac{1}{2}{m}_{\rho }^2\left({\rho }_{\mu }^a{\rho }^{a\mu }\right)-\frac{1}{4}{\mathsf{\Omega }}_{\mu \nu }{\mathsf{\Omega }}^{\mu \nu }-\frac{1}{4}{R}_{\mu \nu }^a{R}^{a\mu \nu }+U_{eff}(\sigma ,\omega ,\rho ),\hfill \end{array}$$

(12)

where ${\mathsf{\Omega }}_{\mu \nu }={\partial }_{\mu }{\omega }_{\nu }-{\partial }_{\nu }{\omega }_{\mu }$ and $R_{\mu \nu }^a={\partial }_{\mu }{\rho }_{\nu }^a-{\partial }_{\nu }{\rho }_{\mu }^a$ are the field tensors and $D_{\mu }={\partial }_{\mu }+i{g}_{\omega }{\omega }_{\mu }+i{g}_{\rho }{I}_3{\tau }^a{\rho }_{\mu }^a$ is the covariant derivative. The effective potential $U_{eff}(\sigma ,\omega ,\rho )$ describes nonlinear self-couplings and mixed-couplings between scalar and vector mesons:

$$\begin{array}{ccc}{U}_{eff}(\sigma ,\omega ,\rho )\hfill & =\hfill & \frac{1}{3}{g}_2{\sigma }^3+\frac{1}{4}{g}_3{\sigma }^4-\frac{1}{4}{c}_3{\left({\omega }_{\mu }{\omega }^{\mu }\right)}^2-\hfill \\ \hfill & -\hfill & {\mathsf{\Lambda }}_V{\left(g_{\omega }{g}_{\rho }\right)}^2\left({\omega }_{\mu }{\omega }^{\mu }\right)\left({\rho }_{\mu }^a{\rho }^{a\mu }\right).\hfill \end{array}$$

(13)

The form of the potential is dictated by the necessity for a realistic description of the properties of finite nuclei and nuclear matter. For the $\sigma$ scalar meson field, its nonlinear self-coupling terms modify the EoS of symmetric nuclear matter at a low density at the vicinity of saturation density $n_0$, while the $\omega$ meson self-interaction term softens the EoS in the high-density limit. The mixed isoscalar-isovector coupling ${\mathsf{\Lambda }}_V{\left(g_{\omega }{g}_{\rho }\right)}^2\left({\omega }_{\mu }{\omega }^{\mu }\right)\left({\rho }_{\mu }^a{\rho }^{a\mu }\right)$ [63,64,65] supplements the isovector part of the model, remodeling the symmetry energy dependence on density. This coupling is determined by the combination of parameters ${\mathsf{\Lambda }}_V,g_{\rho }$ and is fixed to meet the value of the symmetry energy $E_{2,sym}\left(n_b\right)=25.68$ MeV at the baryon density $n_b$, which corresponds to $k_F=1.15f{m}^{-1}$. The TM1parameter set together with ${\mathsf{\Lambda }}_V$ and the proper value of the parameter $g_{\rho }$ are presented in

Table 1. The TM1 parameter set described in [66].

$\sigma$	$m_{\sigma }=511.198$ MeV	$g_{\sigma }=10.029$	$g_2=$ 7.2327 fm${}^{-1}$	$g_3=0.6183$	-
$\omega$	$m_{\omega }=783$ MeV	$g_{\omega }=12.614$	$c_3=71.308$	-	-
$\rho$	$m_{\rho }=770$ MeV	${\mathsf{\Lambda }}_V$	0.0	0.0165	0.03
		$g_{\rho }\left({\mathsf{\Lambda }}_V\right)$	9.2644	10.037	11.10

[66].

The Lagrangian density function $\mathcal{L}$ leads to field equations that are solved in the mean-field approximation. In this approach, the meson field operators are replaced by their expectation values: $\sigma \to <\sigma >\equiv s_0$, ${\omega }_{\mu }\to <{\omega }_{\mu }{\delta }_{\mu 0}>\equiv w_0$, ${\rho }^{\mu a}\to <{\rho }^{\mu a}{\delta }_{0\mu }{\delta }^{3a}>\equiv r_0$. This leads to the equations of motion, which have to be solved self-consistently:

$$m_{\sigma }^2{s}_0+g_2{s}_0^2+3s_0^3=g_{\sigma }{n}_s,$$

(14)

$$m_{\omega }^2{w}_0+c_3{w}_0^3+{\mathsf{\Lambda }}_V{\left(g_{\omega }{g}_{\rho }\right)}^2{r}_0^2{w}_0=g_{\omega }{n}_b,$$

(15)

$$m_{\rho }^2{r}_0+{\mathsf{\Lambda }}_V{\left(g_{\omega }{g}_{\rho }\right)}^2{w}_0^2{r}_0=g_{\rho }{\delta }_a{n}_b,$$

(16)

$$\left(i{\gamma }_{\mu }{\partial }^{\mu }-m_{eff}\left(s_0\right)-g_{\omega }{\gamma }^0{w}_0-g_{\rho }{I}_{3N}{\gamma }^0{\tau }^3{r}_0\right){\psi }_N=0.$$

(17)

The scalar density $n_s$ is given by:

$$n_s=\sum _{N=n,p}\frac{g_s}{{\left(2\pi \right)}^3}{\int }_0^{k_{F,N}}\frac{m_{eff}\left(s_0\right)d^{3}k}{\sqrt{(k^2+m_{eff}^2\left(s_0\right))}},$$

(18)

where $g_s$ is the spin-degeneracy factor, $k_{F,N}$ the nucleon Fermi momentum, $k_{F,N}={(6{\pi }^2/g_s)}^{1/3}{n}_b^{1/3}$, $N=n,p$, and $m_{eff}=m_N-g_{\sigma }{s}_0$ [67]. The obtained solutions allow the determination of the energy $\epsilon$ and pressure P of the system:

$$\begin{array}{ccc}\epsilon \hfill & =\hfill & \frac{1}{2}{m}_{\omega }^2{w}_0^2+\frac{1}{2}{m}_{\rho }^2{r}_0^2+\frac{1}{2}{m}_{\sigma }^2{s}_0^2+\frac{3}{4}{c}_3{w}_0^4+3{\mathsf{\Lambda }}_V{\left(g_{\omega }{g}_{\rho }\right)}^2{w}_0^2{r}_0^2+\hfill \\ \hfill & +\hfill & \frac{1}{3}{g}_2{\sigma }^3+\frac{1}{4}{g}_3{\sigma }^4+\sum _{N=n,p}\frac{2}{{\left(2\pi \right)}^3}{\int }_0^{k_{F,N}}{d}^{3}k\sqrt{k^2+m_{eff}^2\left(s_0\right)},\hfill \end{array}$$

(19)

$$\begin{array}{ccc}P\hfill & =\hfill & \frac{1}{2}{m}_{\omega }^2{w}_0^2+\frac{1}{2}{m}_{\rho }^2{r}_0^2-\frac{1}{2}{m}_{\sigma }^2{s}_0^2+\frac{1}{4}{c}_3{w}_0^4+{\mathsf{\Lambda }}_V{\left(g_{\omega }{g}_{\rho }\right)}^2{w}_0^2{r}_0^2-\hfill \\ \hfill & -\hfill & \frac{1}{3}{g}_2{\sigma }^3-\frac{1}{4}{g}_3{\sigma }^4+\sum _{N=n,p}\frac{2}{3{\left(2\pi \right)}^3}{\int }_0^{k_{F,N}}{d}^{3}k\frac{k^2}{\sqrt{k^2+m_{eff}^2\left(s_0\right)}}.\hfill \end{array}$$

(20)

Given the solutions of the field equations, it is possible to calculate the symmetry energy for different approximations. In the case of parabolic approximation (PA), the following result is obtained:

$$E_{2,sym}\left(n_b\right)=\frac{k_F^2}{6\sqrt{k_F^2+m_{eff}^2}}++\frac{g_{\rho }^2{n}_b}{2(m_{\rho }^2+{\mathsf{\Lambda }}_V{\left(g_{\rho }{g}_{\omega }\right)}^2{w}_0^2)}.$$

(21)

The correction resulting from the Padé approximation requires knowledge of the fourth-order term of the Taylor series. The analytical form of this term obtained in [32] was used in this calculation. The quartic symmetry energy $E_{4,sym}\left(n_b\right)$ term can be decomposed into the sum of the kinetic and potential parts: $E_{4,sym}\left(n_b\right)=E_{4,sym}^{kin}\left(n_b\right)+E_{4,sym}^{pot}\left(n_b\right).$ As a result, the additional contributions to the symmetry energy take the following form:

$$E_{4,sym}^{kin}\left(n_b\right)=\frac{k_F^2}{648}\frac{4m_{eff}^4+11m_{eff}^2{k}_F^2+10k_F^4}{{\left(\sqrt{k_F^2+m_{eff}^2}\right)}^5}$$

(22)

$$\begin{array}{ccc}{E}_{4,sym}^{pot}\left(n_b\right)\hfill & =\hfill & \frac{g_{\rho }^8{n}_b^3}{2Q_{\rho }^4}+\frac{{\mathsf{\Lambda }}_V{g}_{\omega }^4{w}_0^2}{Q_{\omega }}-\hfill \\ \hfill & -\hfill & \frac{g_{\sigma }^2{n}_b{m}_{eff}{k}_F^2}{24Q_{\sigma }{\left(\sqrt{k_F^2+m_{eff}^2}\right)}^3}\times \frac{m_{eff}{k}_F^2}{3{\left(\sqrt{k_F^2+m_{eff}^2}\right)}^3},\hfill \end{array}$$

(23)

where $Q_{\sigma },Q_{\omega }$ and $Q_{\rho }$ are given by relations [32]:

$$Q_{\sigma }=m_{\sigma }^2+g_{\sigma }^2\left(\frac{3n_s}{m_{eff}}-\frac{3n_b}{\sqrt{k_F^2+m_{eff}^2}}\right)+2g_2{s}_0+3g_4{g}_3{s}_0^3,$$

(24)

$$Q_{\omega }=m_{\omega }^2+3c_3{\left(g_{\omega }{w}_0\right)}^2$$

(25)

$$Q_{\rho }=m_{\rho }^2+{\mathsf{\Lambda }}_V{\left(g_{\omega }{g}_{\rho }\right)}^2{w}_0^2.$$

(26)

Correction to the symmetry energy resulting from the use of the Padé approximation contains the term $E_{4,sym}/E_{2,sym}$. This term is the most important factor that remodels the symmetry energy. The theoretical knowledge of the quartic term $E_{4,sym}$ is limited, and unfortunately, it is not experimentally constrained; consequently, there are no constraints on the ratio $E_{4,sym}/E_{2,sym}$. Therefore, any method of estimation of the value of this ratio is much sought after to make a connection with experimental data on the high-density limit of the symmetry energy coming from flow experiments [7]. A rough linear approximation is possible within the considered model. The obtained linear approximation for the $E_{4,sym}/E_{2,sym}$ ratio together with the exact results calculated within this approach are depicted in Figure 1. The least-squares function approximation for an arbitrary finite number of points is used for which the squared correlation coefficient is approximately equal to 0.998.

5.2. Matter of a Neutron Star

The chemical composition of a neutron star is implied by the assumptions of charge neutrality and $\beta$-decay equilibrium. To meet these conditions, the presence of electrons and muons is required in the model of neutron star matter. Consequently, the Lagrangian density must be extended by adding the free lepton Lagrangian ${\mathcal{L}}_l$:

$${\mathcal{L}}_l=\sum _{l=e,\mu }{\overline{\psi }}_l(i{\gamma }^{\mu }{\partial }_{\mu }-m_l){\psi }_l.$$

(27)

The presence of leptons results in additional contributions: $E_l$ to the energy density of the system and $P_l$ to pressure of neutron star matter, where:

$$E_l=\frac{2}{{\left(2\pi \right)}^3}{\int }_0^{k_F}{d}^{3}k\sqrt{k^2+m_l^2}$$

(28)

and:

$$P_l=\frac{2}{3{\left(2\pi \right)}^3}{\int }_0^{k_F}{d}^{3}k\frac{k^2}{\sqrt{k^2+m_l^2}}.$$

(29)

The relations between chemical potentials of the constituents of neutron star matter ${\mu }_{asym}={\mu }_n-{\mu }_p={\mu }_e-{\mu }_{{\nu }_e}$ are obtained by considering the process $p+e^{-}\leftrightarrow n+{\nu }_e$. If neutrinos are not trapped in neutron star matter (${\mu }_{{\nu }_e}=0$), then ${\mu }_{asym}$ reduces to ${\mu }_{asym}={\mu }_n-{\mu }_p={\mu }_e$.

6. Results

6.1. The Effect of the Symmetry Energy for Neutron Star Matter

The possible modifications of the symmetry energy can be studied using a comparative analysis based on the predefined methods of approximation with factoring in each of their specific features. In the paper [68], nuclear symmetry energy was described in two ways. The first one, known as the parabolic approximation (PA), was identified with the second-order term in the series expansion of the binding energy about a point ${\delta }_a=0$, whereas the second was calculated based on the Padé approximation. The two ways of presenting the symmetry energy lead to different results caused by distinct convergence properties of approximations. The knowledge about the asymmetric nuclear matter is crucial for modeling of neutron star interiors. Accordingly, the analysis of nuclear matter with varying degrees of asymmetry allows one to study how this type of change induces changes in the properties of neutron stars. The factor that modifies the symmetry energy in the case of the Padé approximation is a function of baryon density and isospin asymmetry. The results given in [68] have been obtained for any value of isospin asymmetry, whereas in the case of neutron star matter, the isospin asymmetry follows from the imposed conditions of $\beta$-equilibrium and charge neutrality. The symmetry energy is implicit in the EoS of asymmetric nuclear matter, and variation of its form manifests itself as a modification of the value of isospin asymmetry. In the given model, it is perceived as a change caused by variations in the parameter ${\mathsf{\Lambda }}_V$. This step is in accordance with the adjustment of the value of the symmetry energy slope L, with particular regard to the value in line with experimental constraints. Considering data from terrestrial experiments and limitations from astrophysical observations, one can estimate values of the symmetry energy factor $E_{2,sym}\left(n_0\right)=31.6\pm 2.66$ MeV and the slope $L\left(n_0\right)=58.9\pm 16$ MeV [13]. The obtained results are presented in figures that are given below. In all these figures, the values of ${\mathsf{\Lambda }}_V$ correspond to the values of the slope parameter L. The following cases are considered: (${\mathsf{\Lambda }}_V=0,L=110.79$ MeV), (${\mathsf{\Lambda }}_V=0.0165,L=73.34$ MeV), (${\mathsf{\Lambda }}_V=0.03,L=55.76$ MeV). The calculated EoSs for different values of ${\mathsf{\Lambda }}_V$ are presented in Figure 2. The results revealed only moderate sensitivity of the EoS to the form of the symmetry energy. The stiffest EoS was obtained for ${\mathsf{\Lambda }}_V=0$.

The subsequent figures show the changes in the chemical composition of neutron star matter caused by the change in the value of ${\mathsf{\Lambda }}_V$. In Figure 3, the concentrations of protons $Y_p=n_p/n_b$ are given. Proton concentration is one of the most important isospin-sensitive characteristics of neutron star matter, and its value is of special interest as it affects the thermal evolution of neutron stars. It has been shown that if the proton concentration exceeds the critical value $Y_p^{crit}$ given by the relation:

$$Y_p^{crit}=\frac{1}{1+{(1+Y_{e,l}^{1/3})}^3},$$

(30)

where $Y_{e,l}$ equals $n_e/(n_e+n_{\nu })$, then the direct URCA process [69] can occur, activating the fast cooling of neutron stars.

The threshold proton concentration $Y_p^{crit}$ is determined at a given density $n_{crit}$, which fixes the critical neutron star mass in which the proton concentration reaches the level sufficient for the onset of the direct URCA process. The proton fraction highly exceeds the URCA threshold for the model with ${\mathsf{\Lambda }}_V=0$. However, this happens for a class of models with experimentally unacceptable high values of the slope parameter L. $Y_p$ lies in the vicinity of the critical level $Y_p^{crit}$ for other values of ${\mathsf{\Lambda }}_V$. Concentrations of leptons (electrons and muons) that result from the conditions of charge neutrality and $\beta$-equilibrium are presented in Figure 4. Models with reduced values of the slope parameter lead to low concentrations of leptons.

Given the EoS, the neutron star mass M and radius R can be calculated basing on the condition of hydrostatic equilibrium (the Tolman–Oppenheimer–Volkoff (TOV)) equations:

$$\begin{array}{ccc}\frac{dP\left(r\right)}{dr}\hfill & =\hfill & \frac{-G\left(\epsilon \left(r\right)+P\left(r\right)\right)\left(m\left(r\right)+4\pi r^{3}P\left(r\right)\right)}{r^2\left(1-\frac{2Gm\left(r\right)}{r}\right)}\hfill \\ \frac{dm\left(r\right)}{dr}\hfill & =\hfill & 4\pi r^2\epsilon \left(r\right)\hfill \\ \frac{d{n}_b\left(r\right)}{dr}\hfill & =\hfill & 4\pi r^2{\left(1-\frac{2Gm\left(r\right)}{r}\right)}^{-1/2}.\hfill \end{array}$$

(31)

Solutions to the TOV equations in the form of R-M diagrams are presented in Figure 5. Two specific configurations were selected for further analysis: the first with the central density of the order of $2n_0$ and the second one being the maximum mass configuration. The first choice was dictated by an attempt to link the obtained experimental constraints on the form of the symmetry energy with neutron star models. The reported results obtained in HIC experiments [7] suggested that the upper-density limit achievable in flow experiments for the symmetry energy density dependence was within the range $2÷3n_0$. Thus, it is instructive to study some isospin sensitive properties of a neutron star whose central density equals $2n_0$. The selected configurations are marked by dots in Figure 5. The mass difference $\Delta M$ between models with varying values of ${\mathsf{\Lambda }}_V$ were not substantial, and this supported the conclusion that the form of the symmetry energy did not influence the value of mass notably, but it significantly changed the value of the stellar radius (Figure 5).

The solution of TOV equations opens up opportunities for theoretical modeling of the internal structure of a neutron star. Information on the chemical composition of neutron star matter helps to transfer this knowledge to the neutron star interior and enables us to study the change in isospin asymmetry caused by different forms of symmetry energy. The results given in Figure 6 present ${\delta }_a$ dependence on dimensionless baryon density u and ${\delta }_a$ radial dependence for the selected stellar model. The general conclusion was that neutron star matter models with a nonzero value of ${\mathsf{\Lambda }}_V$ were characterized by higher values of isospin asymmetry. The comparison of the radial dependence of ${\delta }_a$ calculated for selected configurations also supported this conclusion. The shown configurations differed in the values of isospin asymmetry.

In order to determine the total relative excess of neutrons in a star, a new parameter ${\Delta }_a$, defined analogously to isospin asymmetry parameter ${\delta }_a$, is introduced:

$${\Delta }_a=\frac{N_n-N_p}{N_n+N_p},$$

(32)

where $N_n$ and $N_p$ are given by $N_i={\int }_0^{R}4\pi \frac{n_i\left(r\right)}{\sqrt{1-\frac{2m\left(r\right)}{r}}}{r}^{2}dr$, $m\left(r\right)$ is the mass enclosed inside the volume of radius r, and $n_i,(i=n,p)$ are the number densities of neutrons and protons. The relations between the mass M and radius R of neutron stars and the value of ${\Delta }_a$ are plotted in Figure 7.

In this figure, dots mark the selected and previously discussed stellar configurations. These configurations differed significantly in the values of ${\Delta }_a$. In general, stars with a central density equal to $2n_0$ have a much higher neutron excess than stars with maximum masses. In the case of a model characterized by a soft form of the symmetry energy (${\mathsf{\Lambda }}_V=0.03$), the excess of neutrons in the star with a central density of $2n_0$ was so large that the star could be virtually considered to be made of neutrons alone, and the content of such a star could be treated as pure neutron matter. This result is confirmed in the right panel of Figure 7, which shows the ${\Delta }_a$ dependence on a stellar radius for the maximum mass configuration and for the configuration with central density equal $2n_0$. The values of ${\Delta }_a$ for the discussed mass configurations and for the fixed ${\mathsf{\Lambda }}_V$ are gathered in

Table 2. Parameters of the selected neutron star configurations.

	${\mathsf{\Lambda }}_{\mathit{V}}$	${\mathbf{\Delta }}_{\mathit{a}}$	${\mathit{n}}_{\mathit{c}}/{\mathit{n}}_0$	M[M${}_{\odot }$]	R (km)	${\mathit{\delta }}_{\mathit{a}}\left(0\right)$
M${}_{max}$	0.03	0.77	6.05	2.11	11.6	0.69
$2\times n_0$	0.03	0.85	2.08	1.13	12.42	0.79
M${}_{max}$	0.0	0.6	5.71	2.17	12.09	0.41
$2\times n_0$	0.0	0.79	2.05	1.33	13.05	0.63

.

6.2. The Effect of the Approximation Method for Neutron Star Matter

The striking feature of the symmetry energy is its sensitivity to the adopted method of approximation. In order to demonstrate this effect, it is necessary to reformulate the EoS to have an explicit dependence on the symmetry energy [62]. Then, one can compare the results obtained for the Taylor and the Padé approximation methods used for the determination of the binding energy of nuclear matter in $\beta$-equilibrium. There is a general approach to the analysis of the EoS of isospin asymmetric nuclear matter, which allows one to present its binding energy $E_N(n_b,Y_p)$ in the parabolic approximation (PA) as $E_N(n_b,Y_p)=E(n_b,Y_p=1/2)+E_{sym}\left(n_b\right){(1-2Y_p)}^2$. This is a generally accepted method that leads to the decomposition of the EoS into the symmetric and isospin-dependent parts. Minimization of the function $E_T(n_b,Y_p)=E_N(n_b,Y_p)+E_e(n_b,Y_e)$, where $E_e(n_b,Y_e)$ is the energy of free electrons, with respect to $Y_p$, provides the condition for $\beta$-equilibrium and leads to the result:

$${\mu }_n-{\mu }_p=4(1-2Y_p)E_{sym}\left(n_b\right).$$

(33)

The pressure of cold $\beta$-equilibrated neutron star matter can be expressed as [70]:

$$P^{eq}\left(n_b\right)=P(n_b,1/2)+n_b(1-2Y_p^{eq})E_{2,sym}\left(n_b\right)\left((1-2Y_p^{eq})\frac{dln{E}_{2,sym}\left(n_b\right)}{dln\left(n_b\right)}+Y_p^{eq})\right),$$

(34)

where $P(n_b,1/2)$ denotes the symmetric matter pressure.

The equilibrium proton fraction was obtained for the TM1 parameterization for ${\mathsf{\Lambda }}_V=0.0165$ and for ${\mathsf{\Lambda }}_V=0.03$. The results showed differences that arose when different methods of approximations were applied (parabolic approximation and Padeapproximants).

The equilibrium proton fractions obtained on the basis of the PA, as well as the Padé approximants are presented in Figure 8. In addition, the dependence on the form of the symmetry energy was introduced through the parameter ${\mathsf{\Lambda }}_V$. Both approaches led to almost identical results near the saturation density, but the obtained solutions differed substantially for higher densities, relevant to neutron star interiors. Equilibrium proton concentrations $Y_p^{eq}$ obtained based on Equation (33) for different ${\mathsf{\Lambda }}_V$ values maintained the tendency revealed in the calculations presented in the previous subsection. Increasing ${\mathsf{\Lambda }}_V$ led to a lower proton concentration. The results that were obtained allows one to calculate the equilibrium pressure and determine the EoS. The EoSs acquired for different values of ${\mathsf{\Lambda }}_V$ are given in Figure 9. The internal structure of a neutron star depends on the location of the transition between its core and crust. The corresponding transition density and pressure depend on the form of the symmetry energy and are also influenced by the approximation method. The obtained results are presented in

Table 3. Crust-core transition density $n_{tr}$ and corresponding values of equilibrium proton fraction $Y_p^{eq}\left(n_{tr}\right)$ and pressure $P^{eq}\left(n_{tr}\right)\left[MeVf{m}^{-3}\right]$ obtained for the parabolic approximation (PA) and Padé approximants.

PA	${\mathsf{\Lambda }}_{\mathit{V}}=0$	${\mathsf{\Lambda }}_{\mathit{V}}=0.0165$	${\mathsf{\Lambda }}_{\mathit{V}}=0.03$	Padé	${\mathsf{\Lambda }}_{\mathit{V}}=0$	${\mathsf{\Lambda }}_{\mathit{V}}=0.165$	${\mathsf{\Lambda }}_{\mathit{V}}=0.03$
n${}_{tr}$ (fm${}^{-3})$	0.0876	0.0932	0.0947	n${}_{tr}$ (fm${}^{-3}\left]\right)$	0.0877	0.0931	0.0946
Y${}_p^{eqT}\left(n_{tr}\right)$	0.037	0.041	0.0430	Y${}_p^{eqP}\left(n_{tr}\right)$	0.041	0.047	0.05
P${}^{eqT}\left(n_{tr}\right)$	1.256	0.378	0.410	P${}^{eqP}\left(n_{tr}\right)$	1.285	0.400	0.423

.

Determining the core-crust transition density $n_{tr}$ is related to the issue of low-density symmetry energy; therefore, the approximation effect, although noticeable, is not too great. However, its role increases with the increase in symmetry energy stiffness. Constraints that have been imposed on the symmetry energy affect the chemical composition of neutron star matter. Through the relation given by Equation (33), they determine the corresponding range of change in proton concentration in $\beta$ equilibrium. The equilibrium concentration of protons is decisive in the cooling process of neutron stars. As the stiffest form of the symmetry energy allows for the fast cooling process, thus the Padé approximation that gives stiffening of the symmetry energy is more favorable than PA.

Theoretical calculations performed on the basis of the effective Skyrme energy density functional indicate the existence of significant correlations between the symmetry energy parameter $E_{2,sym}\left(n_0\right)$ and its slope L [71]. A similar conclusion can be drawn for the complete and partial correlation coefficients for the same variables [31]. This kind of analysis is performed in order to discover possible interdependence among variables. If it is statistically significant, then one can ask about the sensitivity between variables. The aim is to estimate the impact of the type of approximation on the sensitivity of the relative proton concentration $Y_p$ and $E_{2,sym}\left(n_b\right)$ to the symmetry energy slope L. Sensitivity analysis determines how the change in values of an independent variable affect a specific dependent variable for a fixed level of other effects. It can provide information on how the uncertainty in the output of a mathematical model (possibly numerical) can be allocated to different sources of uncertainty in its inputs. Therefore, it may suggest a better or more convenient parameterization of experimental data being modeled. This type of analysis is carried out to detect possible correlations between variables. Therefore, if it is statistically significant and one decides to use a specific form of regression between variables, one can focus on determining the sensitivity between variables. In this spirit, an analysis of the complete and partial correlation coefficients for the variables $E_{2,sym}$ and L was considered in [31]. In [31,71], the linear regression was assumed. In this case, the sensitivity can be calculated. In this paper, the impact of the method of approximation (the Padé approximation and the truncated Taylor series) on the sensitivity of proton concentrations $Y_p$ and $E_{2,sym}$ to the value of the symmetry energy slope L was numerically estimated. The analysis was performed for the baryon density $n_b$ close to $2n_0$. It was found that the sensitivity of proton concentrations of $Y_p$ to the change of L was meaningful, but only weakly depended on the approximation.

When L increased from 96.5 to 223, then for the Padé approximation, $Y_p$ increased from approximately 0.14 to 0.46, while for the truncated Taylor series, it increased from 0.11 to 0.41, cf. Figure 10. The sensitivity of $E_{2,sym}$ to the change of L was moderate, and in the Padé case, it was negligibly higher then in the Taylor case (cf. Figure 10). Numerically, if L increased from 96.5 to 223, then for Padé approximation $E_{2,sym}$ increased from approximately 48.76 to 74.82, while for the truncated Taylor series, it increased from 48.74 to 74.81. For high values of $n_b$, the difference in the sensitivities to the symmetry energy slope L between both approaches increased and even changed sign. This suggested that the range of applicability of the used approximations was exceeded.

7. Conclusions

The goal of this paper was to analyze the way the symmetry energy influenced the isospin dependent part of the EoS of dense nuclear matter. The differences resulting from the approximation method were particularly informative. To demonstrate this, it was necessary to determine the EoS of nuclear matter in $\beta$-equilibrium in a way that explicitly took into account the form of the symmetry energy. Such calculations were performed, and selected isospin-dependent characteristics of neutron star matter were calculated using two methods of approximation: the Taylor expansion and the Padé approximants. The obtained results pointed out that different approximations led to modifications of the equilibrium proton fractions and EoSs, especially in their high-density limit. This is important in the case of attempts to extrapolate the results of experiments for densities relevant for neutron star cores. Thus, the interpretation of experimental data for completeness should also include the aspects of approximation methods, which could lead to inherently different results. The Padé approximation resulted in a form of the symmetry energy that depends on isospin asymmetry ${\delta }_a$ and on the ratio $E_{4,sym}\left(n_b\right)/E_{2,sym}\left(n_b\right)$. However, there were virtually no experimental constraints on $E_{4,sym}\left(n_b\right)$. The results obtained in this paper were based on the model for which the factor $E_{4,sym}/E_{2,sym}$ was estimated at around $0.07$ [68]. Even in this case, differences in equilibrium concentrations of protons, as well as differences in equilibrium EoSs were quite big. Various theoretical models gave estimates of $E_{4,sym}$ that were significantly higher than those obtained for the considered model [72]. The study showed that an increase in the value of $E_{4,sym}$ led to the stiffening of the symmetry energy, and in this case, the obtained form of the symmetry energy better met the experimental constraints [68].

Neutron stars, being objects with high isospin asymmetry, are described by the EoS, of which a significant part is the symmetry energy. Thus, it is instructive to know how the neutron star interior is modified by changing the symmetry energy stiffness. In particular, by changing the concentration of protons, it affected the chemical composition of the star, leading to objects with a higher asymmetry value for the softer form of the symmetry energy.

Summing up: The experimentally determined limitations on the symmetry energy dependence on density and the related attempt to transfer the obtained results to the ground of the EoS of the asymmetric nuclear matter should take into account both the form of the symmetry energy and the method of approximation used. This would lead to a more precise determination of the real isovector part of the EoS.