A Statistical Approach to Neutron Star Crust-Core Transition Density and Pressure

1. Introduction

The nuclear symmetry energy belongs to the key factors that define the problem of the exact neutron star’s internal structure and, to some extent, determine its solution. How and in what range symmetry energy controls the emergence of different phases of nuclear matter is one of the main topics of current theoretical research in nuclear physics and astrophysics. The uncertainties in the internal structure of a neutron star, which is expected to exhibit nuclear matter at different physical states, are mainly due to the limited knowledge of the equation of state (EoS) of such a matter being in extreme physical conditions of density, temperature, and isospin asymmetry. Without experimental data extracted at such extreme conditions, it is necessary to use models that meet the results of ground-based experiments and reproduce nuclear matter’s saturation properties. Such models yield considerable uncertainty when extrapolated and applied to densities relevant to neutron stars. There are many dubious points in the modeling of neutron stars. One of the most critical concerns is the precise description of the crust-core crossing boundary and, thus, the extent of the crust. The neutron star matter EoS allows for obtaining a neutron star’s hydrostatic model, and its general stratification distinguishes three layers: the outermost is the atmosphere, then the crust, which splits into the inner and outer parts. The inner crust extends outward to the well-determined neutron drip density ${\rho }_{drip}=4\times {10}^{11}$ g/cm³. The very inner part of a neutron star is a liquid core comprising interacting neutrons in $\beta$-equilibrium with the admixture of protons and electrons. Theoretical considerations point to the complex structure of a neutron star’s inner crust. It consists of atomic nuclei with significant neutron excess immersed in a gas of free neutrons and relativistic degenerate electrons. Depending on the density, atomic nuclei have different shapes, being spherical in most of the inner crust. Calculations suggest that non-spherical configurations of nuclei in the crust’s deepest layers become energetically favorable, forming the so-called pasta phase [1,2,3,4,5]. This complex structure of the inner crust transforms into its equally complicated EoS. The missing precise physical model that allows for constructing an accurate EoS adequate to describe the asymmetric nuclear matter in the full range of densities characteristic for a neutron star and which correctly reproduces its properties forces the use of approximate methods. However, these methods allow for determining only approximately a neutron star model and, among other properties, the crust-core boundary’s location. Often, the physical model is described by various statistical characteristics. One of the most general is the joint probability distribution of all variables needed to describe the physical phenomenon under study. However, finding such a distribution by proposing a theoretical model is generally impossible. Therefore, the initial stage in constructing a physical model is selecting a regression model(s) between variables suspected of being essential in describing a physical phenomenon. Finding this regression model(s), in turn, helps capture the correct form of the physical theoretical model. Searching for different regression functions between different variables may be the initial stage of its search in case of ignorance of the fundamental formulation of the physical model. The obvious help is appropriate statistical analysis. One of the most elegant statistical methods is the maximum likelihood method (MLM) [6] and the resulting Akaike information criterion (AIC) [7,8,9,10]. In general, AIC helps search for the true statistical model from which the data visible in the observation are generated. Between the models accepted for analysis, the statistical model (e.g., regression model) closest to the unknown accurate statistical model gives the highest probability of producing the observed data. Section 4.3 is devoted to the AIC criterion, which selects a regression model between the crust-core transition pressure $P_t$ and the characteristics of the system’s energy. The statistics that measure the goodness of fit of the dependent variable to the data for a specific group of independent variables in a linear regression model is the coefficient of determination $R^2$. In this paper the adjusted $R^2$, $R_{adj}^2$, is also used [10], Section 4.4. $R_{adj}^2$ helps to eliminate the overestimation of the model obtained by applying $R^2$, i.e., $R_{adj}^2$ may have a maximum and may start to decline as the number of effects in the regression model increases. The maximum of $R_{adj}^2$ indicates the place where the expansion of the regression model should be stopped so as not to overfit the model in the sample when compared to the unknown model in the population (theoretical model). In this paper, the $R_{adj}^2$ is considered an auxiliary criterion in searching for the optimal regression model. Another method used in this paper for selecting the appropriate regression model is the backward elimination method [10], (Section 4.5), which allows for choosing a regression model with factors that have a significant statistical impact on the goodness of fit of the dependent variable to the data. It is good if all these methods point to the same group of factors and produce the same regression model, although there is generally no guarantee that this will happen. In this paper, the AIC criterion for selecting the regression model is preferred, as it gives the highest probability of the appearance of the particular data.

The regression model between the transition pressure $P_t$ and the characteristics of the energy of the system estimated in this paper using the AIC method is a particular characteristic of the sought true physical model that is expected to describe nuclear and astrophysical observations correctly. One of the quantities characterizing a physical system is the crust-core transition density $n_t$. The proposed approach allows for determining the most probable value of the transition density ${\tilde{n}}_t$ related to the selected regression model for the analyzed sample of the RMF models (Section 5.1 and Section 5.2).

2. The Inner Edge of a Neutron Star Inner Crust

The location of the crust-core boundary in a neutron star can be specified if accurate models describing the matter of the crust and core are known. Generally, a hydrostatic equilibrium equation supplemented with a proper form of the EoS can provide valuable clues about the neutron star’s internal structure. However, in a neutron star’s inner crust, one can deal with a form of nuclear matter whose a priori predictions are not obvious. Model calculations indicate the possibility of a very complex, nonhomogenous phase called nuclear pasta, which further complicates the form of the equation of state of this matter. Due to its highly complex structure, determining the EoS of matter in this layer of a neutron star is problematic and burdened with very high uncertainty. Thus, it has become necessary to develop alternative methods that lead to the estimation of the transition density at which homogeneous matter becomes unstable against small density fluctuations, indicating the beginning of the formation of the nucleus clusters. In the following, the location of the inner boundary of the neutron star’s inner crust is determined based on thermodynamic methods [13,14,15,16], which require that the system meets the stability condition given by the pair of inequalities:

$$-{\left(\frac{\partial P}{\partial v}\right)}_{\mu }>0,-{\left(\frac{\partial {\mu }_{asym}}{\partial q_c}\right)}_v>0,$$

(1)

if not, it loses stability against small density fluctuations. In the above inequalities, v and $q_c$ are volume and charge per baryon number, P is the total pressure of the system, and ${\mu }_{asym}={\mu }_n-{\mu }_p$ is the difference of neutron and proton’s chemical potentials. The energy of nuclear matter considered in terms of binding energy (EoS) is given by the relation

$$E(n_b,\delta )=\frac{\epsilon (n_b,\delta )}{n_b}-M,$$

(2)

where the energy density $\epsilon (n_b,\delta )$ of the system is a function that depends on baryon density $n_b=n_n+n_p$ and the isospin asymmetry parameter $\delta$, M is the nucleon mass. It is expected that the function $E(n_b,\delta )$ can be represented by its Taylor series, which under expansion to the fourth-order around $\delta =0$ takes the following form

$$\begin{array}{ccc}\hfill E(n_b,\delta )& =& \sum _{n=0}^{\infty }{E}_{2n}\left(n_b\right){\delta }^{2n}\hfill \\ & =& E_0\left(n_b\right)+E_2\left(n_b\right){\delta }^2+E_4\left(n_b\right){\delta }^4+\dots .\hfill \end{array}$$

(3)

Coefficients of the series (3) are functions of baryon density and denote the binding energy of the symmetric matter $E_0\left(n_b\right)$, the symmetry energy $E_2\left(n_b\right)\equiv E_{sym,2}\left(n_b\right)$ and the fourth-order symmetry energy $E_4\left(n_b\right)\equiv E_{sym,4}\left(n_b\right)$. The simplest case considers only the second-order term in (3), and it is known as the parabolic approximation. Using the dependence $\delta =1-2Y_p$, where $Y_p=n_p/n_b$ is the relative proton concentration, the following relation for the isospin-dependent part of the binding energy can be obtained

$$E_{N,asym}(n_b,Y_p)=E_{sym,2}\left(n_b\right){(1-2Y_p)}^2+E_{sym,4}\left(n_b\right){(1-2Y_p)}^4.$$

The energy per baryon of relativistic electrons has the form

$$E_e\left(n_b\right)=\frac{3}{4}\hslash c{\left(3{\pi }^2{n}_b\right)}^{1/3}{Y}_e^{1/3}.$$

The charge neutrality condition demands that $Y_e=Y_p$. Thus, the total energy per baryon of the matter in the core is given by

$$E_{Tot}=E_0\left(n_b\right)+E_{N,asym}(n_b,Y_p)+E_e(n_b,Y_p).$$

Minimization of $E_{Tot}(n_b,Y_p)$ with respect to $Y_p$ gives the $\beta$ equilibrium condition

$$\begin{array}{ccc}\hfill {\mu }_e={\mu }_n-{\mu }_p& =& -\frac{\partial E_{Tot}(n_b,Y_p)}{\partial Y_p}=4(1-2Y_p)E_{sym,2}\left(n_b\right)+\hfill \\ & +& 8{(1-2Y_p)}^3{E}_{sym,4}\left(n_b\right).\hfill \end{array}$$

(4)

For the chemical potential of relativistic electrons ${\mu }_e=\hslash c{\left(3{\pi }^2{n}_b\right)}^{1/3}{Y}_e^{1/3}$, the above-given condition allows one to determine the equilibrium proton fraction $Y_p^{eq}$

$$\begin{array}{ccc}\hfill \hslash c{\left(3{\pi }^2{n}_b\right)}^{1/3}{Y}_p{\left(n_b\right)}^{1/3}& =& 4(1-2Y_p\left(n_b\right))E_{sym,2}\left(n_b\right)\hfill \\ & +& 8{(1-2Y_p\left(n_p\right))}^3{E}_{sym,4}\left(n_b\right).\hfill \end{array}$$

(5)

The condition $-{\left(\frac{\partial {\mu }_{asym}}{\partial q_c}\right)}_v>0$ is usually satisfied, whereas the inequality $-{\left(\frac{\partial P}{\partial v}\right)}_{\mu }>0$ can be expressed by requiring the expression $V_{ther}$ to be positive

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

(6)

where $E(n_b,Y_p)$ is the binding energy of nuclear matter. Solving Eq.(5) and (6) allows for determining the value of the transition density $n_t$ and the corresponding proton concentration value $Y_p^{eq}\left(n_t\right)=Y_t$. Using the thermodynamic relation

$$P=n_b^2\frac{\partial E(n_b,Y_p)}{\partial n_b}$$

(7)

to calculate the pressure of the npe system of particles results in a total pressure that is the sum of contributions from nucleons ($P_N$) and electrons ($P_e$), $P_{Tot}=P_N+P_e$. The calculations made for the transition density $n_t$ and the corresponding $Y_t$ value can lead to the equation for the pressure at the crust-core boundary.

$$\begin{array}{ccc}\hfill P_t\left(n_t\right)& =& n_t^2{\left.\frac{d{E}_0\left(n_b\right)}{d{n}_b}\right|}_{n_t}+n_t^2{(1-2Y_t)}^2\left({\left.\frac{d{E}_{sym,2}\left(n_b\right)}{d{n}_b}\right|}_{n_t}+{(1-2Y_t)}^2{\left.\frac{d{E}_{sym,4}\left(n_b\right)}{d{n}_b}\right|}_{n_t}\right)+\hfill \\ & +& n_t{Y}_t(1-2Y_t)\left(E_{sym,2}\left(n_t\right)+2E_{sym,4}\left(n_t\right){(1-2Y_t)}^2\right).\hfill \end{array}$$

(8)

In general, it is expected that higher-order terms in the expansion (3) have to be included to obtain a more accurate description of the binding energy of systems with a significant value of the isospin asymmetry. In this case, an improvement in the accuracy of the obtained solution is expected. In further analysis, each function $E_0\left(n_b\right)$, $E_{sym,2}\left(n_b\right)$ and $E_{sym,4}\left(n_b\right)$ is represented by a Taylor series expansion around $n_0$. This procedure can be presented in the general form as

$$E_j\left(n_b\right)=\sum _{i=0}^{\infty }{C}_j^i{\left(\frac{n_b-n_0}{3n_0}\right)}^i.$$

(9)

The index j distinguishes between symmetric $\delta =0$ and asymmetric $\delta \ne 0$ nuclear matter. The case of symmetric nuclear matter is denoted by $j=0$, and $E_0\left(n_b\right)$ means the binding energy of symmetric nuclear matter. The case $j=2$ corresponds to the second-order symmetry energy $E_{sym,2}\left(n_b\right)$ and $j=4$ the fourth-order symmetry energy $E_{sym,4}\left(n_b\right)$. The expansion coefficients

$$C_i^j={\left(3n_0\right)}^i\frac{1}{i!}\frac{d^i{E}_j\left(n_b\right)}{d{n}_b^i}{|}_{n_0}$$

(10)

represent the following characteristics of nuclear matter: $C_0^0\equiv E_0\left(n_0\right)$ - the binding energy per nucleon of symmetric nuclear matter at saturation density $n_0$, the nuclear matter incompressibility $C_2^0\equiv K_0$, $C_0^2\equiv E_{sym,2}\left(n_0\right)$ - the symmetry energy at the saturation density, $C_1^2\equiv L_{sym,2}$ - the second order symmetry energy slope, $C_2^2\equiv K_{sym,2}$ - curvature of the second-order symmetry energy, $C_1^4\equiv L_{sym,4}$ - the fourth-order symmetry energy slope, $C_2^4\equiv K_{sym,4}$ - curvature of the fourth-order symmetry energy.

Applying the Taylor series expansions of functions $E_0\left(n_b\right)$, $E_{sym,2}\left(n_b\right)$ and $E_{sym,4}\left(n_b\right)$, it is possible to obtain an approximate value of the pressure at the crust-core boundary

$$\begin{array}{ccc}\hfill P_{app}\left(n_t\right)& \approx & \frac{n_t^2(n_t-n_0)}{9n_0^2}\left(K_0+K_{sym,2}{\delta }_t^2+K_{sym,4}{\delta }_t^4\right)+\hfill \\ & +& L_{sym,2}\left(\frac{n_t(n_t-n_0)Y_p\left(n_t\right){\delta }_t}{3n_0}+\frac{n_t^2{\delta }_t^2}{3n_0}\right)\hfill \\ & +& n_t{Y}_p\left(n_t\right){\delta }_t\left(E_{sym,2}+2E_{sym,4}{\delta }_t^2\right)+\hfill \\ & +& L_{sym,4}\left(\frac{2n_t(n_t-n_0)Y_p\left(n_t\right){\delta }_t^3}{3n_0}+\frac{n_t^2{\delta }_t^4}{3n_0}\right).\hfill \end{array}$$

(11)

Another approximate form of the expression defining the pressure can be obtained, assuming that $\delta$ equals 1, which leads to $Y_p=0$ and corresponds to the case of pure neutron matter.

$$\begin{array}{ccc}\hfill P_{app}\left(n_t\right)& \approx & {\left(\frac{n_t}{3n_0}\right)}^2(n_t-n_0)\left(K_0+K_{sym,2}+K_{sym,4}\right)+\hfill \\ & +& \frac{n_t}{3n_0}{n}_t\left(L_{sym,2}+L_{sym,4}\right).\hfill \end{array}$$

(12)

3. Determination of the EoS

The determination of the EoS is based on the Lagrangian density function being the sum of free baryon and meson fields part ${\mathcal{L}}_0$ and the part ${\mathcal{L}}_{int}$ describing the interaction. Individual parts are given in the following forms:

$$\begin{array}{ccc}\hfill {\mathcal{L}}_0& =& \overline{\psi }(i{\gamma }^{\mu }{\partial }_{\mu }-M)\psi +\frac{1}{2}({\partial }^{\mu }\sigma {\partial }_{\mu }\sigma -m_{\sigma }^2{\sigma }^2)-\frac{1}{4}{F}^{\mu \nu }{F}_{\mu \nu }\hfill \\ & +& \frac{1}{2}{m}_{\omega }^2{\omega }_{\mu }{\omega }^{\mu }-\frac{1}{4}{\mathbf{B}}^{\mu \nu }{\mathbf{B}}_{\mu \nu }+\frac{1}{2}{m}_{\rho }^2{\overrightarrow{\rho }}_{\mu }·{\overrightarrow{\rho }}^{\mu },\hfill \end{array}$$

(13)

where $\sigma$, ${\omega }_{\mu }$ and ${\overrightarrow{\rho }}_{\mu }$ represent the scalar-isoscalar $\sigma$, vector-isoscalar $\omega$, and vector-isovector $\rho$ meson fields, respectively and $\psi$ is the isodoublet nucleon field, $F_{\mu \nu }$ and ${\mathbf{B}}_{\mu \nu }$ are field tensors defined as $F_{\mu \nu }={\partial }_{\mu }{\omega }_{\nu }-{\partial }_{\nu }{\omega }_{\mu }$ and ${\mathbf{B}}_{\mu \nu }={\partial }_{\mu }{\overrightarrow{\rho }}_{\nu }-{\partial }_{\nu }{\overrightarrow{\rho }}_{\mu }$,

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{int}& =& \overline{\psi }\left(g_{\sigma }\sigma -(g_{\omega }{\omega }_{\mu }+\frac{1}{2}{g}_{\rho }\overrightarrow{\tau }·{\overrightarrow{\rho }}_{\mu }){\gamma }^{\mu }\right)\psi -\frac{A}{3}{\sigma }^3-\frac{B}{4}{\sigma }^4+\frac{C}{4}{\left(g_{\omega }^2{\omega }_{\mu }{\omega }^{\mu }\right)}^2\hfill \\ & +& g_{\sigma }{g}_{\omega }^2\sigma \left({\omega }_{\mu }{\omega }^{\mu }\right)({\alpha }_1+\frac{1}{2}{\alpha }_1^{\prime }{g}_{\sigma }\sigma )+g_{\sigma }\sigma g_{\rho }^2\left({\overrightarrow{\rho }}_{\mu }{\overrightarrow{\rho }}^{\mu }\right)({\alpha }_2+\frac{1}{2}{\alpha }_2^{\prime }{g}_{\sigma }\sigma )+\hfill \\ & +& \frac{1}{2}{\alpha }_3^{\prime }{\left(g_{\omega }{g}_{\rho }\right)}^2\left({\omega }_{\mu }{\omega }^{\mu }\right)\left({\overrightarrow{\rho }}_{\mu }{\overrightarrow{\rho }}^{\mu }\right).\hfill \end{array}$$

(14)

The Lagrangian density function ${\mathcal{L}}_{int}$ contains the Yukawa couplings between the nucleons and the meson and collects various nonlinear meson interaction terms. The individual coupling constants determine the strength of the meson interactions. The equations of motion derived based on the given above Lagrangian density function $\mathcal{L}={\mathcal{L}}_0+{\mathcal{L}}_{int}$ were solved in the mean-field approximation. In this approach, meson fields are separated into classical components and quantum fluctuations; the quantum fluctuation terms vanish, and only classical components remain. The mean field limit, in the case of a static and spherically symmetric system, leads to the following relations:

$$\begin{array}{c}\hfill \sigma \to \langle \sigma \rangle \equiv s\\ \hfill {\omega }^{\mu }\to \langle \omega \rangle \equiv \langle {\omega }_0\rangle {\delta }^{\mu 0}\equiv {\omega }_0\\ \hfill {\overrightarrow{\rho }}^{\mu }\to \langle {\rho }_3\rangle \equiv \langle {\rho }_{0,3}\rangle {\delta }^{\mu 0}\equiv r_{0,3}.\end{array}$$

(15)

The mesons are coupled to the nucleon sources, which are also replaced by their expectation values in the mean-field ground state. The solution of the equations of motion allows one to calculate the energy density of the system

$$\begin{array}{ccc}\hfill \epsilon & =& \frac{1}{2}{m}_{\sigma }^2{s}^2+\frac{A}{3}{s}^3+\frac{B}{4}{s}^4-\frac{1}{2}{m}_{\omega }^2{\omega }_0^2-\frac{C}{4}{\left(g_{\omega }^2{\omega }_0^2\right)}^2+g_{\omega }{\omega }_0{n}_b-\frac{1}{2}{m}_{\rho }^2{r}_{0,3}^2+g_{\rho }{r}_{0,3}{n}_{3b}\hfill \\ & -& g_{\sigma }s{\left(g_{\omega }{\omega }_0\right)}^2\left({\alpha }_1+\frac{1}{2}{\alpha }_1^{\prime }{g}_{\sigma }s\right)-g_{\sigma }s{\left(g_{\rho }{r}_{0,3}\right)}^2\left({\alpha }_2+\frac{1}{2}{\alpha }_2^{\prime }{g}_{\sigma }s\right)-\frac{1}{2}{\alpha }_3^{\prime }{\left(g_{\omega }{\omega }_0\right)}^2{\left(g_{\rho }{r}_{0,3}\right)}^2\hfill \\ & +& \sum _{j=n,p}\frac{g}{2{\pi }^2}{\int }_0^{k_{Fj}}{k}^2\sqrt{k^2+M_{\mathrm{eff},j}^2}dk,\hfill \end{array}$$

(16)

where $M_{\mathrm{eff}}=M-g_{\sigma }s$ denotes the effective nucleon mass and $n_{3b}=\langle \overline{\psi }{\gamma }^0{\tau }_3\psi \rangle =n_p-n_n$, and g represents the number of degrees of freedom. The nonlinear meson interaction terms necessary for constructing a correct nuclear matter EoS alter both the isoscalar and isovector sectors [17,18]. The calculations were carried out in the framework of relativistic mean field (RMF) theory. In this approach, the nuclear many-body problem is considered a relativistic system of baryons and mesons. In the original Walecka model, only scalar-isocalar $\sigma$ (attractive) and vector-isocalar $\omega$ (repulsive) mesons [19,20] were involved to account for the saturation properties of symmetric nuclear matter. This model was then extended with the vector-isovector meson $\rho$ and subjected to further modifications, leading to more sophisticated models containing various nonlinear self and mixed meson interaction terms [21]. Specifying this model in such an extended form allows one to successfully reproduce some ground-state properties of finite nuclei and nuclear matter. The implemented modifications increase the usefulness of the models in satisfactory descriptions of the properties of asymmetric nuclear matter [22,23]. The properties of nuclear matter determined based on RMF models rely on selected groups of parameters that are the research subject presented in papers [23,24]. The acceptance of a given parameterization depends on the degree of compliance of the determined properties of symmetric and asymmetric nuclear matter with the constraints resulting from the analysis of experimental data. The choice of experimental constraints in the case of symmetrical matter ($\delta =0$) considers the nuclear matter incompressibility at saturation density $K_0$ in the range $190-270$ MeV [25,26,27] the skewness coefficient $-Q$ in the range $200-1200$ MeV [28], the pressure $P\left(n_b\right)$ in density ranges $(2n_0,5n_0)$ and $(1.5n_0,2.5n_0)$ [29,30]. Considering the asymmetric nuclear matter [31], experimental constraints apply to the coefficients characterizing the density dependence of the symmetry energy. One can specify the following limitation ranges: symmetry energy coefficient $E_{\mathrm{sym}}\left(n_0\right)$ - (25 MeV - 35 MeV) and (30 MeV - 35 MeV) [32], symmetry energy slope $L_0$ calculated at $n_0$ - (25 MeV - 115 MeV) [33,34], volume part of isospin incompressibility $K_{\tau ,v}^0$ at $n_0$ - ( -700 MeV - -400 MeV) [23,35,36] and the ratio of the symmetry energy in $n_0/2$ to its value in $n_0$ - (0.57 - 0 .86) [37].

The RMF models applied in the analysis performed in this paper can be characterized and distinguished by different types of nonlinear couplings between mesons. It becomes possible to divide all models into three groups. Group I includes the BSR [38] and FSUGZ03, FSUGZ06 [39] models with the following types of mixed meson couplings $\sigma -{\omega }^2$, ${\sigma }^2-{\omega }^2$, $\sigma -{\rho }^2$, ${\sigma }^2-{\rho }^2$, ${\omega }^2-{\rho }^2$. Group II of the BKA [40], G2 [41] and G$2^{★}$[42] models includes $\sigma -{\omega }^2$, ${\sigma }^2-{\omega }^2$, $\sigma -{\rho }^2$ non-linear terms. Group III FSUGold [18], FSUGold4 [43], IU FSU, XS [44] and TM1 [45] is characterized by ${\omega }^2-{\rho }^2$. Values of parameters for individual models and saturation properties of symmetric and asymmetric nuclear matter are collected in the papers [46,47]. The energy density of the system given by Eq.(16) encodes the correct form of the symmetry energy.

4. Regression Analysis

Various concepts that belong to the category of measuring the goodness of fit of quality of statistical modeling have been developed, including $R^2$, adjusted $R^2$, which represents some attempt to adjust for the number of parameters in the model, AIC, and statistical backward elimination. The approaches are not always equivalent, and using different methods allows for a better understanding of which factors in the regression models are the most important.

4.1. The Regression Function

Regression analysis is a valuable method for estimating the relationships between a dependent random variable Y (response) and independent variables $X_1,X_2,...X_k$ (factors). Let the regression model has the form:

$$\begin{array}{ccc}\hfill Y& =& \mathrm{E}\left[Y\right|X_1,X_2,...,X_k]+E\hfill \\ & =& {\alpha }_0+{\alpha }_1{X}_1+{\alpha }_2{X}_2+...+{\alpha }_k{X}_k+E,\hfill \end{array}$$

(17)

where E denotes the random error, and $\mathrm{E}\left[Y\right|X_1,X_2,...,X_k]$ is the conditional expectation value of Y, ${\alpha }_i$, $i=0,1,2,...,k$ are the structural parameters. The parameter ${\alpha }_0$ is called the intercept. Considering a sample of N models chosen randomly from a population of models, the regression model (17) can be estimated by:

$$\begin{array}{c}\hfill Y={\widehat{\alpha }}_0+{\widehat{\alpha }}_1{X}_1+{\widehat{\alpha }}_2{X}_2+...+{\widehat{\alpha }}_k{X}_k+\widehat{E}.\end{array}$$

(18)

Here ${\widehat{\alpha }}_0$, ${\widehat{\alpha }}_1$, ${\widehat{\alpha }}_2,...,{\widehat{\alpha }}_k$ are the estimators of the structural parameters ${\alpha }_0,{\alpha }_1,{\alpha }_2,...,{\alpha }_k$ of a particular regression model and $\widehat{E}$ is the estimator of the error term E in Eq.(17). The error term $\widehat{E}$ variance is denoted by $MSE$ (mean squared error) throughout this paper.

Given a linear regression model with k factors, the null hypothesis

$$\begin{array}{c}\hfill H_0:{\alpha }_1={\alpha }_2=...={\alpha }_k=0,\end{array}$$

(19)

is a question about the irrelevance of the correlation between the dependant variable Y and the group of independent variables $X_i$, $i=1,2,...,k.$

4.2. The Consistency Assumption for Considered Models

This paper assumes that every theoretical point in the sample of $N=23$ models is estimated consistently, that is, without any bias, at least asymptotically. Therefore, every theoretical point on the scatter diagram coincides with the estimate obtained for n of hypothetical experiments testing this model. It follows that in the limits $n\to \infty$ and for all the population of models, the finite sample error $\widehat{E}$ tends to E. Therefore, the requirement to use the method is the assumption that it is possible to determine the values of the estimators’ model parameters from the experiment. Each model introduced into the analysis satisfies as many experimental constraints as possible. This group is an optimal sample of models in this paper.

4.3. The Akaike Information Criterion Analysis

The Akaike information criterion (AIC) [9] is very useful in mining the most probable appearance of the observed sample with the simultaneous model extensions limiting. Let the data $y=(y_1,y_2,...,y_N)$ be generated by the true but unknown regression model g for the random variable Y (to simplify the notation, only the values $y_i$ of the response variable Y are written). Consider a regression model $f\equiv f(Y,A_k)$ with a vector parameter $A_k$, as a candidate for describing the investigated interdependence between the dependent variable Y and the group of factors. $A_k$ is a free parameter of the regression model f as all its components ${\alpha }_1$, ${\alpha }_2$, ... , ${\alpha }_k$ can be put to zero in the null hypothesis (19). To select a better regression model f for the response variable Y and explanatory variables $X_1,X_2,...X_k$ with a parameter $A_k$, the following form of AIC is used

$$\begin{array}{c}\hfill AIC(f,A_k)=-2\ln L\left({\widehat{A}}_k\right)+2(k+1).\end{array}$$

(20)

Here, $L\left(A_k\right)\equiv L\left(y\right|A_k)$ denotes the likelihood function corresponding to the model f for a N-dimensional sample, ${\widehat{A}}_k$ is a maximum likelihood method (MLM) estimator of the parameter $A_k$, and $k+1$ is the number of the estimated structural parameters in the regression model, i.e., the vector of slope coefficients $A_k$ = $({\alpha }_1,{\alpha }_2,...,{\alpha }_k)$ plus the intercept ${\alpha }_0$. The maximization of the log-likelihood function $\ln L\left(y\right|A_k)$ is in mean equivalent to maximizing the expectation value ${\mathrm{E}}_g\left[\ln f(Y,A_k)\right]$ calculated for the true model g [9]. As the unknown parameter $A_k$ is replaced by its MLM estimator ${\widehat{A}}_k$, thus, instead of ${\mathrm{E}}_g\left[\ln f(Y,A_k)\right]$ the expectation value $Q_k\equiv {\mathrm{E}}_{g,h_{A_k}}\left[\ln f(Y,{\widehat{A}}_k)\right]$ is maximized, where $h_{A_k}$ is the distribution $h_{A_k}\left({\widehat{A}}_k\right)$ of the estimator ${\widehat{A}}_k$. The maximization of $Q_k$ is equivalent to the minimization of $-2N{Q}_k$, where N is the dimension of the sample. Because $AIC(f,A_k)$ is approximately an unbiased estimator of $-2N{Q}_k$,9], the model that minimizes $AIC(f,A_k)$ is the candidate for the searched model. This can be confirmed by considering the Kullback-Leibler (K-L) distance between the models f and g [9]:

$$\begin{array}{ccc}\hfill D(g,f)& =& {\mathrm{E}}_g\left[\ln g\left(Y\right)\right]-{\mathrm{E}}_g\left[\ln f(Y,A_k)\right].\hfill \end{array}$$

(21)

As ${\mathrm{E}}_g\left[\ln g\left(Y\right)\right]$ is constant, the minimization of $AIC(f,A_k)$ implies the selection of the model that minimizes the K-L distance of the, chosen for the statistical analysis, model f from the unknown true model g. The model selected by AIC gives the highest probability of occurrence of the observed sample. Details concerning the AIC model selection procedure can be found in [9].

4.4. The Characteristics of the Regressions. Coefficient of Determination $R^2$ and adjusted $R_{adj}^2$

In statistics, the sums of squares, used in a regression analysisol, measure the variability in data. It reveals the dispersion of data points concerning the mean and how much the response variable differs from the predicted values. For a given dependent variable Y, it is convenient to define the total sum of squares SSY= ${\sum }_{i=1}^N{(Y_i-\overline{Y})}^2$, which is the sum of squares of deviations of the observed $Y_i$ from their mean $\overline{Y}$. $SSY$ is often partitioned to the sum of squares due to regression, $SSR={\sum }_{i=1}^N{({\widehat{Y}}_i-\overline{Y})}^2$ and due to error, $SSE={\sum }_{i=1}^N{(Y_i-{\widehat{Y}}_i)}^2$, where $SSE$ is called the residual (error) sum of squares. This leads to the ANOVA equation for the linear regression, $SSY=SSR+SSE$[10]. The characteristics of the regressions in use are the mean square due to regression $MSR=SSR/d{f}_{SSR}$, the mean squared error $MSE=SSE/d{f}_{SSE}$ and the coefficient of determination:

$$\begin{array}{c}\hfill R^2=\frac{SSR}{SSY}=1-\frac{SSE}{SSY}\in \langle 0,1\rangle ,\end{array}$$

(22)

which measures the ratio of the variability of the dependent variable explained by the regression to the overall variability of this variable. Here, $d{f}_{SSR}=k$ and $d{f}_{SSE}=N-k-1$ are the number of degrees of freedom for $SSR$ and $SSE$, respectively. The descriptive limits of the correlation strength in this paper are assumed to be: $0.1<R^2<0.25$ for weak correlation, and, if $R^2\ge 0.64$ ($\left|R\right|\ge 0.8$) for strong correlation. In one-dimensional linear regression $Y=a+bX+\widehat{E}$, the sign of the Pearson linear correlation coefficient $r_{YX}$ between Y and X,10] equals the sign of b, $r_{YX}=\mathrm{sgn}\left(b\right)\left|R\right|$. An additional statistics to the coefficient of determination $R^2$, which considers the number of parameters in the model is the adjusted coefficient of determination $R_{adj}^2$, defined as [11]

$$\begin{array}{c}\hfill R_{adj}^2=1-\frac{N-i}{SSY}MSE=1-\frac{N-i}{N-K}(1-R^2),\end{array}$$

(23)

where the mean square error $MSE=\frac{SSE}{N-K}$ and N is the number of observations used to match the model, K is the number of parameters in the model, including the intercept, $K=k+1$, and i equals 1 when the model has an offset (intercept) and 0 otherwise. $R_{adj}^2$ starts decreasing when the model has too many parameters. The moment in which $R_{adj}^2$ starts to drop is a signal that the model no longer needs to be developed.

4.5. The Backward Selection Method

4.5.1. $F_{partial}$ Statistics

One might be tempted to compare models with different numbers of parameters. To this aim, a convenient notation is $SS{E}_k\equiv SSE(X_1,X_2,$$...,X_k)$. Let the model includes $k^{\prime }$ factors $X_1$, $X_2$, ..., $X_k$, $X_{k^{\prime }}$, $k^{\prime }=k+1.$ To determine the significance of introducing an additional factor $X_{k^{\prime }}$, the partial $F_{\mathrm{p}}$-statistics ($\mathrm{p}$ for "partial") for the models ($X_1,X_2,...,X_k,X_{k^{\prime }}$) and ($X_1,X_2,...,X_k$) can introduced [10]:

$$\begin{array}{ccc}\hfill F_{\mathrm{p}}& =& \frac{(SS{R}_{k^{\prime }}-SS{R}_k)/(d{f}_{SS{R}_{k^{\prime }}}-d{f}_{SS{R}_k})}{MS{E}_{k^{\prime }}},\hfill \end{array}$$

(24)

where, according to the introduced notation, $MS{E}_{k^{\prime }}=SS{E}_{k^{\prime }}/d{f}_{SS{E}_{k^{\prime }}}.$ The statistics $F_{\mathrm{p}}$ is a random variable on the sample space, having the F-distribution $F_{k^{\prime }-k,N-k^{\prime }-1}$ with $k^{\prime }-k=d{f}_{SS{R}_{k^{\prime }}}-d{f}_{SS{R}_k}$ and $N-k^{\prime }-1=d{f}_{SS{E}_{k^{\prime }}}$ . For a particular sample, $F_{\mathrm{p}}$ takes the observed value $F_p^{obs}$. In this case, the empirical significance level (the p-value) can be calculated:

$$\begin{array}{c}\hfill p=\mathrm{Prob}(F_{\mathrm{p}}\ge F_{\mathrm{p}}^{obs}).\end{array}$$

(25)

If in the observed sample, for a chosen significance level $\alpha$, $p>\alpha$, then there is no reason to reject the null hypothesis $H_0:{\alpha }_{k^{\prime }}=0$, which now reads: "the lower model fits the observed data as well as the higher model," i.e., the sample gives no incentives to extend the model. Here, the adjectives lower and higher reflect the number of factors. The lower model is rejected if $p\le \alpha$. With the hierarchical development of the regression model in the paper, the value of the partial statistics $F_{\mathrm{p}}$ in the observed sample and the corresponding p-values, (Eq.(25)) for the factor added last can be determined at all stages of model construction. It can be shown that the $F_{\mathrm{p}}$ test (24) for the significance of a one-variable extension of a model with k variables coincides with Student’s t-test for the null hypothesis for the structural parameter ${\alpha }_{k^{\prime }}=0$ (the p-values of the tests are the same).

4.5.2. Backward Elimination Method

The backward elimination regression method is a statistical procedure of the model selection used to reduce the number of less significant variables [10]. The selection should start with a possible most complete model and simplify it until it turns out that all the remaining variables have a substantial impact on the accuracy of fitting the model regression to empirical data. As a tool, the partial $F_{\mathrm{p}}$ value (24) (or the empirical significance level p-value (25)) for each variable in the model is calculated. When comparing the highest value of the empirical significance level p with the value of the previously chosen significance level $\alpha$ (for the variable to remain in the model, e.g., $\alpha =0.01,0.05$), it is possible to decide whether to remove or keep the considered variable. The procedure can be repeated after deciding to neglect a given variable until obtaining a model where all the values of the estimators of the model’s structural parameters are statistically significant.

5. Discussion

5.1. The Results of the Selection of the Regression Models

The analysis performed in this paper uses a sample of the most reliable RMF models that describe nuclear matter whose high credibility follows from the fact that they meet the largest number of experimental constraints. Based on these models, nuclear matter EoSs given in terms of binding energy $E(n_b,\delta )$ (2) were constructed. In the first step of the analysis, the function $E(n_b,\delta )$ is approximated by its Taylor series expansion around $\delta =0$. This leads to the separation of the symmetric $E_0(n_b,0)\equiv E_0\left(n_b\right)$ and asymmetric $E_{asym}(n_b,\delta )=E_{sym,2}\left(n_b\right){\delta }^2+E_{sym,4}\left(n_b\right){\delta }^4+\dots$ parts of the EoS and allows one to consider the asymmetric part of the EoS at different levels of approximation. The coefficients of the expansion depend on baryon density. The analysis was carried out for the symmetry energy given by the parabolic approximation and for the case when the description of asymmetric matter additionally considers the fourth-order symmetry energy term. The transition pressure at the neutron star crust-core boundary following (8) decisively depends on the functions $E_0\left(n_b\right)$, $E_{sym,2}\left(n_b\right)$, and $E_{sym,4}\left(n_b\right)$. The approximate expression for the transition pressure given in terms of the defined expansion coefficients has the form given by Eq. (11). All variables that enter this formula form the set of explanatory variables. In the parabolic approximation, it contains the following terms:

$$\begin{array}{c}\hfill \left(K_0,E_{2,Y{\delta }_2}\equiv E_{sym,2}\left(n_0\right)Y_{t,2}{\delta }_{t,2},L_{2,{\delta }_2^2}\equiv L_{sym,2}{\delta }_{t,2}^2,\right.\\ \hfill \left.L_{2,Y{\delta }_2}\equiv L_{sym,2}{Y}_{t,2}{\delta }_{t,2},K_{2,{\delta }_2^2}\equiv K_{sym,2}{\delta }_{t,2}^2\right)\end{array}$$

(26)

and in the fourth-order approximation:

$$\begin{array}{c}\hfill \left(K_0,E_{2,Y{\delta }_{24}}\equiv E_{sym,2}\left(n_0\right)Y_{t,24}{\delta }_{t,24},L_{2,{\delta }_{24}^2}\equiv L_{sym,2}{\delta }_{t,24}^2,\right.\\ \hfill L_{2,Y{\delta }_{24}}\equiv L_{sym,2}{Y}_{t,24}{\delta }_{t,24},K_{2,{\delta }_{24}^2}\equiv K_{sym,2}{\delta }_{t,24}^2,\\ \hfill E_{4,Y{\delta }_{24}}\equiv E_{sym,4}\left(n_0\right)Y_{t,24}{\delta }_{t,24}^3,L_{4,{\delta }_{24}^4}\equiv L_{sym,2}{\delta }_{t,24}^2,\\ \hfill \left.L_{4,Y{\delta }_{24}}\equiv L_{sym,4}{Y}_{t,24}{\delta }_{t,24}^3,K_{2,{\delta }_{24}^2}\equiv K_{sym,4}{\delta }_{t,24}^4\right).\end{array}$$

(27)

These variables serve as input parameters in the regression analysis. The nuclear matter at the crust-core transition boundary is highly isospin asymmetric, and thus, an additional approximation consisting in taking ${\delta }_t=1$, which corresponds to pure neutron matter, was also adopted. The description of nuclear matter was based on a selected group of RMF models. Although this is an optimal sample of models that meets many experimental constraints, none is the final true physical model, i.e., having all the necessary components in the correct form. Since the true physical model is unknown, its search can start at a selected basic stage. This means to provide for this physical model statistical evidence based on a regression analysis, which will reproduce the given sample of RMF models with the highest probability. The procedure of evaluating regression models, called model selection, has been applied. The selected model should be the one that provides an adequate representation of the data. However, it must be emphasized that the selected model is not desirable to be represented by the maximal number of explanatory variables. The selection analysis identifies the explanatory variables for the selected regression model. Different selection procedures, such as the AIC method and the $R_{adj}^2$, and the backward elimination method yielded the chosen regression model (Section 4.3-Section 4.5). The analysis covers several cases. The first concerns approximations used to describe the symmetry energy, namely the parabolic approximation $E_{sym,2}\left(n_b\right)$ (Table 1), and the one that also considers the contribution from the fourth-order term $E_{sym,2}\left(n_b\right)+E_{sym,4}\left(n_b\right)$ (Table 2). In each table, the collected results of regression models for a different number of explanatory variables are given. The results for the pure neutron matter (the isospin asymmetry $\delta =1$) obtained for the parabolic approximation are given in Table 3. In Table 4, the results for the fourth-order case are gathered.

Selection analysis indicates that when the parabolic approximation describes the symmetry energy for both considered values of the $\delta$ parameter, the minimum AIC value applies to the maximal model, meaning that the selected regression model covers the entire set of explanatory variables (26) (Table 1 and Table 3). In the case $\delta \ne 1$ a global AIC minimum appears (Table 1) for five explanatory variables denoted by $char=cha{r}_2\equiv (K_0,E_{2,Y{\delta }_2},L_{2,{\delta }_2^2},L_{2,Y{\delta }_2},K_{2,{\delta }_2^2})$. For pure neutron matter ($\delta =1$), there are three explanatory variables $char=cha{r}_{2;\delta =1}\equiv (K_0,L_{sym,2},K_{sym,2})$ (Table 3). This situation changes when the symmetry energy function is the sum of $E_{sym,2}\left(n_b\right)$ and $E_{sym,4}\left(n_b\right)$. In this case, for $\delta \ne 1$, a global AIC minimum appears (Table 2) for a model with six explanatory variables $char=cha{r}_{24}\equiv (K_0,E_{2,Y{\delta }_{24}},L_{2,{\delta }_{24}^2},L_{2,Y{\delta }_{24}},K_{2,{\delta }_{24}^2},$$E_{4,Y{\delta }_{24}^3})$ selected from the set (27). For pure neutron matter ($\delta =1$), there are three AIC selected explanatory variables $char=cha{r}_{24;\delta =1}\equiv (L_{sym,2},K_{sym,2},K_{sym,4})$ (Table 4). The results obtained by the AIC method coincide with the results for $R_{adj}^2$ in three out of four cases, and the selected model is characterized by the maximal value of $R_{adj}^2$ (see Table 1,Table 2,Table 3). An exception is for ${\delta }_{24}=1$ (Table 4), for which there is a minor compatibility violation at the third significant figure. However, for the AIC selected model, there still is a local maximum of $R_{adj}^2$.

When multiplied by the $Y_t$ factor, the roles of explanatory variables from sets (26) and (27) are practically negligible due to the small $Y_t$ value. Therefore, for the regression analysis with only one factor, the explanatory variable multiplied by $Y_t$ is not considered. Otherwise, an artificial effect of a statistical nature may happen, suggesting a good fit to the data for a model with an insignificant variable.

Results of the employed AIC and $R_{adj}^2$ model selection techniques are presented in Figure 1 and Figure 2. Both figures depict values of AIC and $R_{adj}^2$ against the number $\mathrm{n}$ of explanatory variables.

5.2. The Most Probable Value of the Transition Density

The exact numerical values of the transition pressure $P_t\left(n_t\right)$ is calculated according to Eq.(8) (see Table 6) and $P_{app}\left(n_t\right)$ is its approximated form given by Eq.(11). The following equality follows:

$$P_t\left(n_t\right)=P_{app}\left(n_t\right)+R,$$

(28)

where R is the remainder of the Taylor series expansion. This equation is valid for every RMF model from the considered sample. Treating the regression model as an alternative way to represent the data makes it possible to approximate the transition pressure in the sample by the sum of the function $P_{fit}$ plus the error (residual) term $\widehat{E}$ (see Eq.(17)):

$$P_t\left(n_t\right)=P_{fit}(char;\widehat{{\alpha }_j}{|}_{j=0}^k)+\widehat{E},$$

(29)

where $P_{fit}(char;\widehat{{\alpha }_j}{|}_{j=0}^k)$ is the regression function with a general form $P_{fit}(char,{\widehat{\alpha }}_i)=\widehat{{\alpha }_0}+f(char;\widehat{{\alpha }_j}{|}_{j=1}^k)$, $\widehat{{\alpha }_0}$ denotes the intercept term. The estimate of the variance of E is $MSE$, which is the variance of $\widehat{E}$ and $\sqrt{MSE}$ is its standard deviation (see Section 4.4) . The regression function in the parabolic case has the form (see (26)):

$$\begin{array}{ccc}\hfill P_{fit}& =& {\widehat{\alpha }}_0+{\widehat{\alpha }}_1{K}_0+{\widehat{\alpha }}_2{E}_{2,Y{\delta }_2}+{\widehat{\alpha }}_3{L}_{2,{\delta }_2^2}+{\widehat{\alpha }}_4{L}_{2,Y{\delta }_2}+{\widehat{\alpha }}_5{K}_{2,{\delta }_2^2}\hfill \end{array}$$

(30)

and when the fourth-order term is included in the description of the symmetry energy, $P_{fit}$ is given by (see (27)):

$$\begin{array}{ccc}\hfill P_{fit}& =& {\widehat{\alpha }}_0+{\widehat{\alpha }}_1{K}_0+{\widehat{\alpha }}_2{E}_{2,Y{\delta }_{24}}+{\widehat{\alpha }}_3{L}_{2,{\delta }_{24}^2}+{\widehat{\alpha }}_4{L}_{2,Y{\delta }_{24}}+{\widehat{\alpha }}_5{K}_{2,{\delta }_{24}^2}\hfill \\ & +& {\widehat{\alpha }}_6{E}_{4,Y{\delta }_{24}^3}+{\widehat{\alpha }}_7{L}_{4,{\delta }_{24}^4}+{\widehat{\alpha }}_8{L}_{4,Y{\delta }_{24}^3}+\widehat{{\alpha }_9}{K}_{4,{\delta }_{24}^4}.\hfill \end{array}$$

(31)

The basis for determining the most probable value of the transition density is the assumption of the validity of the following relation, which is the consequence of the two possible representations of the data, namely (28) and (29),

$$P_{app}\left(n_t\right)=P_{fit}(char;\widehat{{\alpha }_j}{|}_{j=0}^k),$$

(32)

where the function on the RHS is given by relation (30) in the parabolic approximation or (31) in the fourth-order case. This requires the appearance of the constant $\widehat{{\alpha }_0}$, which results from a different form of the coefficients in $P_{app}\left(n_t\right)$ and the constant coefficients $\widehat{{\alpha }_i}$, $i=1,2,...,k$, in $P_{fit}(char;\widehat{{\alpha }_j}{|}_{j=0}^k)$. In addition, the residual standard deviation $\sqrt{MSE}$ is a mean estimate of the remainder R in the range of the considered transition density $n_t$. The parameters of the regression model including $\widehat{{\alpha }_0}$ depend only implicitly on the transition density $n_t$ and in the limit $n_t\to 0$, ${\widehat{\alpha }}_0\to 0.$

The regression model selected by the AIC, which gives the most probable appearance of the sample, allows one based on Eq.(32) to calculate the most probable value of the transition density $n_t={\tilde{n}}_t$. The transition density ${\tilde{n}}_t$ resulting from the performed regression analysis is considered as the best approximation of the transition density implied by the true regression model, and as a consequence, it should characterize the true physical model. The regression function for the AIC and $R_{adj}^2$ selected regression model in the fourth-order approximation has the form

$$\begin{array}{ccc}\hfill P_{fit}& =& {\widehat{\alpha }}_0+{\widehat{\alpha }}_1{K}_0+{\widehat{\alpha }}_2{E}_{2,Y{\delta }_{24}}+{\widehat{\alpha }}_3{L}_{2,{\delta }_{24}^2}+{\widehat{\alpha }}_4{L}_{2,Y{\delta }_{24}}+{\widehat{\alpha }}_5{K}_{2,{\delta }_{24}^2}+{\widehat{\alpha }}_6{E}_{4,Y{\delta }_{24}^3}\hfill \\ & =& -4.9549+0.002479K_0+8.009438E_{2,Y{\delta }_{24}}+0.08421L_{2,{\delta }_{24}^2}\hfill \\ & -& 3.4876L_{2,Y{\delta }_{24}}-0.004782K_{2,{\delta }_{24}^2}-49.4E_{4,Y{\delta }_{24}^3}\hfill \end{array}$$

(33)

The above regression model is also confirmed by the backward analysis procedure applied to the set of factors (27) as all values of the estimators of the structural parameters ${\alpha }_j$, $j=0,1,2,...,k=6$ of the regression model are statistically significant at level $\alpha =0.05$. It is assumed that the significance levels of introducing a variable into the model and keeping it in the model are the same. The other characteristics of the selected regression model with regression function given by Eq.(33) are given in Table 5.

In the parabolic case, the selected by the AIC and $R_{adj}^2$ regression model has the following regression function:

$$\begin{array}{ccc}\hfill P_{fit}& =& {\widehat{\alpha }}_0+{\widehat{\alpha }}_1{K}_0+{\widehat{\alpha }}_2{E}_{2,Y{\delta }_2}+{\widehat{\alpha }}_3{L}_{2,{\delta }_2^2}+{\widehat{\alpha }}_4{L}_{2,Y{\delta }_2}+{\widehat{\alpha }}_5{K}_{2,{\delta }_2^2}\hfill \\ & =& -2.9636+0.002974K_0+2.4789E_{2,Y{\delta }_2}+0.04254L_{2,{\delta }_2^2}\hfill \\ & -& 1.1365L_{2,Y{\delta }_2}-0.003259K_{2,{\delta }_2^2}.\hfill \end{array}$$

(34)

This regression model is confirmed by the backward analysis procedure applied to the set of factors (26) at the significance level $\alpha =0.065$. The values of the estimators of the structural parameters other than ${\alpha }_1$ for the factor $K_0$ and ${\alpha }_4$ for the factor $L_{2,Y{\delta }_2}$ remain in the model at the significance level lower than $\alpha =0.05$. The other characteristics of the selected regression model with regression function given in Eq.(34) are given in Table 5.

Figure 3 shows the values of means $\overline{{\tilde{n}}_t}$ of the most probable density values ${\tilde{n}}_t$ (see Table 1, Table 2, Table 3 and Table 4), obtained for $\delta ={\delta }_t$ and $\delta =1$ for the two considered cases of symmetry energy approximations as a function of the number of explanatory variables $\mathrm{n}$ that characterize a given regression model selected by the AIC and $R_{adj}^2$ methods.

The crucial relation for the further construction of a true physical model is ${\tilde{n}}_t\left(L_{sym,2}\right)$. This relation for the model selected by the AIC and $R_{adj}^2$ in the fourth-order approximation with the regression function (33) is shown in Figure 4.

6. Conclusions

A regression model selected according to the AIC and $R_{adj}^2$ model selection procedures can be used to identify and support a correct, from the physical and experimental point of view, model. The results obtained include, among others, the importance of approximations used to describe the symmetry energy. Suppose the isospin asymmetry parameter has the value resulting from the adopted model $\delta ={\delta }_t$. In this case, the selected regression models corresponding to the extreme values of AIC and $R_{adj}^2$, both for the parabolic approximation and considering the fourth-order term, include as explanatory variables $K_0$ and $L_{sym,2}$, $K_{sym,2}$, $E_{sym,2}$ multiplied by functions of ${\delta }_t$ (see (26) and (27)). The preferred regression model that incorporates contributions from the fourth-order symmetry energy term weakly depends on the fourth-order symmetry energy characteristics (depending only on $E_{sym,4}{Y}_t{\delta }_{24}^3$ ). A different result is obtained, assuming the regression analysis is done for pure neutron matter. Then, in the parabolic approximation, the set of independent variables is the maximum set and includes the factors $K_0$, $K_{sym,2}$ and $L_{sym,2}$. However, in the case when the symmetry energy is given by the formula $E_{sym,2}\left(n_b\right)+E_{sym,4}\left(n_b\right)$, the set of explanatory variables does not include the characteristics of symmetric nuclear matter namely its incompressibility $K_0$. In this case, the regression model is described by the following independent variables: $L_{sym,2},L_{sym,4},K_{sym,2}$ and $K_{sym,4}$. The selected regression models for the parabolic approximation, regardless of the considered values of $\delta$, always contain the maximal set of independent variables. An additional conclusion concerns the value of the most probable transition density, which for pure neutron matter when taking into account the contribution of the fourth-order symmetry energy, is of significantly lower value $\overline{{\tilde{n}}_t}=0.05876\pm 0.00572$ fm⁻³. This means that, in this case, the crust-core boundary is moved to much lower densities. After obtaining the transition density in the fourth-order approximation for the AIC and $R_{adj}^2$ selected model (with the regression function (33)), it is possible to examine the relationship between ${\tilde{n}}_t$ and $L_{sym,2}$. The obtained results confirm the existence of anti-correlation between these quantities (see Figure 4). Moreover, the mean square error $MSE$ for the exponential fit is much lower than in the linear case. Thus, the exponential fit is better. This could suggest the existence of a valuable shift of the transition boundary towards densities approaching the saturation density in the case of nuclear matter models characterized by a low $L_{sym,2}$ value.