A Statistical Approach to Neutron Stars’ Crust–Core Transition Density and Pressure

In this paper, a regression model between neutron star crust–core pressure and the symmetry energy characteristics was estimated using the Akaike information criterion and the adjusted coefficient of determination Radj2. The most probable value of the transition density, which should characterize the crust–core environment of the sought physical neutron star model, was determined based on the obtained regression function. An anti-correlation was found between this transition density and the main characteristic of the symmetry energy, i.e., its slope L.


Introduction
Nuclear symmetry energy is a key factor that defines the problem of the neutron star's exact internal structure and, to some extent, determines 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's matter EoS allows a neutron star's hydrostatic model to be obtained, and its general stratification distinguishes three layers: the outermost is the atmosphere and then the crust, which splits into the inner and outer parts.The inner crust extends outward to the well-determined neutron drip density ρ drip = 4 × 10 11 g/cm 3 .The very inner part of a neutron star is a liquid core comprising interacting neutrons in β-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 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 of a neutron star and that correctly reproduces its properties forces the use of approximate methods.
However, these methods allow for only approximately determining a neutron star model and, among other properties, the crust-core boundary's location.Often, the physical model is described with 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 the 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, the AIC helps search for the actual statistical model from which the data visible in the observation are generated.Between the models accepted for analysis, the statistical model (e.g., the regression model) closest to the unknown accurate statistical model gives the highest probability of producing the observed data.Section 4.2 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 2 adj , is also used [10]; see Appendix A.1.R 2 adj helps to eliminate the overestimation of the model obtained by applying R 2 ; i.e., R 2 adj may have a maximum.It may decline as the number of regression model's effects increases.The maximum of R 2 adj indicates 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).This paper considers the R 2 adj to be 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] (Appendix A.2), 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 system's energy characteristics estimated in this paper using the AIC method is a particular characteristic of the sought actual physical model 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 ñt related to the selected regression model for the analyzed sample of the RMF models (Sections 5.1 and 5.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 to estimate the transition density at which homogeneous matter becomes unstable against small density fluctuations, indicating the beginning of the formation of the nucleus clusters.Below, the location of the inner boundary of the neutron star's inner crust is determined based on thermodynamic methods [11][12][13][14], which require that the system meets the stability condition given by the pair of inequalities: Otherwise, 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 system's total pressure, and µ asym = µ n − µ p is the difference in neutrons' and protons' chemical potentials.The energy of nuclear matter considered in terms of binding energy (EoS) is given by the relation where the energy density ε(n b , δ) of the system is a function that depends on baryon density n b = n n + n p and the isospin asymmetry parameter δ; M is the nucleon mass.It is expected that the function E(n b , δ) can be represented by its Taylor series, which, under expansion to the fourth order around δ = 0, takes the following form Coefficients of the series (3) are functions of baryon density and denote the binding energy of the symmetric matter E 0 (n b ), the symmetry energy E 2 (n b ) ≡ E sym,2 (n b ), and the fourthorder symmetry energy E 4 (n b ) ≡ E sym,4 (n b ).The simplest case considers only the secondorder term in (3), and it is known as the parabolic approximation.Using the dependence δ = 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: The energy per baryon of relativistic electrons has the form 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 The minimization of E Tot (n b , Y p ) with respect to Y p gives the β equilibrium condition For the chemical potential of relativistic electrons µ e = hc(3π 2 n b ) 1/3 Y 1/3 e , the condition given above allows one to determine the equilibrium proton fraction Y eq p hc(3π The condition − ∂µ asym ∂q c v > 0 is usually satisfied, whereas the inequality − ∂P ∂v µ > 0 can be expressed by requiring the expression V ther to be positive where E(n b , Y p ) is the binding energy of nuclear matter.Solving Equations ( 5) and ( 6) allows for determining the value of the transition density n t and the corresponding proton concentration value Y eq p (n t ) = Y t .Using the thermodynamic relation to calculate the pressure of the n-p-e 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.
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 (n b ), E sym,2 (n b ), and E sym,4 (n b ) is represented by a Taylor series expansion around n 0 .This procedure can be presented in the general form as The index j distinguishes between symmetric δ = 0 and asymmetric δ = 0 nuclear matter.The case of symmetric nuclear matter is denoted by j = 0, and E 0 (n b ) means the binding energy of symmetric nuclear matter.The case j = 2 corresponds to the second-order symmetry energy E sym,2 (n b ) and j = 4 the fourth-order symmetry energy E sym,4 (n b ).The expansion coefficients represent the following characteristics of nuclear matter: C 0 0 ≡ E 0 (n 0 ) is the binding energy per nucleon of symmetric nuclear matter at a saturation density n 0 , the nuclear matter incompressibility C 0 2 ≡ K 0 , C 2 0 ≡ E sym,2 (n 0 ) is the symmetry energy at the saturation density, C 2 1 ≡ L sym,2 is the second-order symmetry energy slope, C 2 2 ≡ K sym,2 is the curvature of the second-order symmetry energy, C 4  1 ≡ L sym,4 is the fourth-order symmetry energy slope, and C 4 2 ≡ K sym,4 is the curvature of the fourth-order symmetry energy.By applying the Taylor series expansions of the functions E 0 (n b ), E sym,2 (n b ), and E sym,4 (n b ), it is possible to obtain the approximate value of the pressure at the crust-core boundary In the case when the density dependence of the symmetry energy is given by the parabolic approximation, Equation (11) reduces to the following form Another approximate form of the expression defining the pressure can be obtained, assuming that δ equals 1, which leads to Y p = 0 and corresponds to the case of pure neutron matter The above equation reduces to a very simple form for the parabolic approximation of the symmetry energy: Only when the transition density reaches values equal to the saturation density n 0 does the dependence of pressure P t on parameters characterizing the incompressibility of nuclear matter disappear, and a straightforward relation P app ≈ 1/3n 0 L sym,2 is obtained.

Determination of the EoS
The determination of the EoS is based on the Lagrangian density function that is the sum of free baryon and meson fields part L 0 and the part L int describing the interaction.The individual parts are given in the following forms: where σ, ω µ , and ρ µ represent the scalar-isoscalar σ, vector-isoscalar ω, and vector-isovector ρ meson fields, respectively ,and ψ is the isodoublet nucleon field, F µν and B µν are field tensors defined as The Lagrangian density function 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 above Lagrangian density function L = L 0 + L int were solved in the mean-field approximation.This approach separates meson fields into classical components and quantum fluctuations; the quantum fluctuation terms vanish, and only classical parts remain.The mean field limit, in the case of a static and a spherically symmetric system, leads to the following relations: The mesons are coupled to the nucleon sources, which are also replaced by their expectation values in the mean-field ground state.The solution to the equations of motion allows one to calculate the energy density of the system where M eff = M − g σ s denotes the effective nucleon mass and n 3b = ψγ 0 τ 3 ψ = 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 [15,16].The calculations were carried out in the framework of relativistic mean field (RMF) theory.This approach considers the nuclear many-body problem a relativistic system of baryons and mesons.In the original Walecka model, only scalar-isoscalar σ (attractive) and vector-isoscalar ω (repulsive) mesons [17,18] were involved in accounting for the saturation properties of symmetric nuclear matter.This model was then extended with the vector-isovector meson ρ and subjected to further modifications, leading to more sophisticated models containing various nonlinear self and mixed meson interaction terms [19].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 [20,21].The properties of nuclear matter determined based on RMF models rely on selected groups of parameters that are the research subject presented in papers [21,22].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 (δ = 0) considers the nuclear matter's incompressibility at saturation density K 0 in the range of 190-270 MeV [23][24][25], the skewness coefficient is Q in the range 200-1200 MeV [26], and the pressure P(n b ) is in density ranges of (2n 0 , 5n 0 ) and (1.5n 0 , 2.5n 0 ) [27,28].Considering the asymmetric nuclear matter [29], 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 sym (n 0 ) −(25-35 MeV) and (30-35 MeV) [30], symmetry energy slope L 0 calculated at n 0 −(25-115 MeV) [31,32], volume part of isospin incompressibility K 0 τ,v at n 0 −(−700-−400 MeV) [21,33,34], and the ratio of the symmetry energy in n 0 /2 to its value in n 0 −(0.57-0.86)[35].
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 [36] and FSUGZ03, FSUGZ06 [37] models with the following types of mixed meson couplings: Group II of the BKA [38], G2 [39] and G2 [40] models includes the σ − ω 2 , σ 2 − ω 2 , σ − ρ 2 non-linear terms.Group III FSUGold [16], FSUGold4 [41], IU FSU, XS [42] and TM1 [43] is characterized by ω 2 − ρ 2 .The values of parameters for individual models and saturation properties of symmetric and asymmetric nuclear matter are collected in the papers [44,45].The energy density of the system given by Equation (18) encodes the correct form of the symmetry energy.

Regression Analysis
Various concepts that belong to the category of measuring the goodness of fit of the 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 necessary information on the AIC method used in the paper to select the regression model is presented below.The basic information on other methods is given in the Appendix A. 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.

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 → ∞ and for all the population of models, the finite sample error Ê 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.

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 limiting model extensions.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 ≡ 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 α 1 , α 2 , . . ., α k can be made zero in the null hypothesis (A3) (Appendix A).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 Here, L(A k ) ≡ L(y|A k ) denotes the likelihood function corresponding to the model f for a N-dimensional sample, 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 = (α 1 , α 2 , . . ., α k ) plus the intercept α 0 .The mean of the maximization of the log-likelihood function lnL(y|A k ) is equivalent to maximizing the expectation value E g [ln f (Y, A k )] calculated for the true model g [9].As the unknown parameter A k is replaced by its MLM estimator A k , thus, instead of 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]: As E g [ln g(Y)] is constant, the minimization of AIC( f , A k ) implies the selection of the model that minimizes the K-L distance chosen for the statistical analysis is model f from the unknown true model g.Details concerning the AIC model-selection procedure can be found in [9].

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 , δ) (2) were constructed.In the first step of the analysis, the function E(n b , δ) is approximated by its Taylor series expansion around δ = 0.This leads to the separation of the symmetric E 0 (n b , 0) ≡ E 0 (n b ) and asymmetric E asym (n b , δ) = E sym,2 (n b )δ 2 + E sym,4 (n b )δ 4 + . . .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 (n b ), E sym,2 (n b ), and E sym,4 (n b ).The approximate expression for the transition pressure given in terms of the defined expansion coefficients has the form given by Equation (11).All variables that enter this formula form the set of explanatory variables.In the parabolic approximation, it contains the following terms: and in the fourth-order approximation: These variables serve as input parameters in the regression analysis.The nuclear matter at the crust-core transition boundary is highly isospin-asymmetric.Thus, an additional approximation consisting of taking δ 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., one with all the necessary components in the correct form.Since the true physical model is unknown, the search for it can start at a selected basic stage.This means providing statistical evidence for this physical model 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 it is not desirable that the selected model is 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 2 adj and the backward elimination method, yielded the chosen regression model (Section 4.2 andAppendices A.1 and A.2).The analysis covers several cases.The first concerns approximations used to describe the symmetry energy, namely the parabolic approximation E sym,2 (n b ) (Table 1) and the one that also considers the contribution from the fourth-order term E sym,2 (n b ) + E sym,4 (n b ) (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 δ = 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 δ parameter, the minimum AIC value applies to the maximal model, meaning that the selected regression model covers the entire set of explanatory variables (21) (Tables 1 and 3).In the case that δ = 1, a global AIC minimum appears (Table 1) for five explanatory variables denoted by char ).For pure neutron matter (δ = 1), there are three explanatory variables char = char 2;δ=1 ≡ (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 (n b ) and E sym,4 (n b ).In this case, for δ = 1, a global AIC minimum appears ( ) selected from the set (22).For pure neutron matter (δ = 1), there are three AIC selected explanatory variables char = char 24;δ=1 ≡ (L sym,2 , K sym,2 , K sym,4 ) (Table 4).The results obtained using the AIC method coincide with the results for R 2 adj in three out of four cases, and the selected model is characterized by the maximal value of R 2 adj (see Tables 1-3).An exception is for δ 24 = 1 (Table 4), for which there is a minor compatibility violation in the third significant figure.However, for the AIC selected model, there still is a local maximum of R 2 adj .When multiplied by the Y t factor, the roles of explanatory variables from sets ( 21) and ( 22) are practically negligible due to the small Y t value.Therefore, the explanatory variable multiplied by Y t is not considered for the regression analysis with only one factor.Otherwise, an artificial effect of a statistical nature may occur, suggesting a good fit to the data for a model with an insignificant variable.
Results of the employed AIC and R 2 adj model-selection techniques are presented in Figures 1 and 2. Both figures depict values of AIC and R 2 adj against the number n of explanatory variables.

The Most Probable Value of the Transition Density
The exact numerical values of the transition pressure P t (n t ) are calculated according to Equation (8) (see Table 5) and P app (n t ) is its approximated form given in Equation (11).The following equality is as follows: 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 with the sum of the function P f it plus the error (residual) term E (see Equation (A2)): where P f it (char; α j | k j=0 ) is the regression function with a general form P f it (char, ), α 0 denotes the intercept term.The estimate of the variance of E is MSE, which is the variance of E, and √ MSE is its standard deviation (see Appendix A.1) .The regression function in the parabolic case has the form (see ( 21)) and when the fourth-order term is included in the description of the symmetry energy, P f it is given by (see ( 22)) The basis for determining the most probable value of the transition density ñt is the assumption of the validity of the relation, which is the consequence of the two possible representations of the transition pressure P t (n t ), given by Equations ( 23) and ( 24), where the function on the RHS is given in relation (25) in the parabolic approximation or (26) in the fourth-order case.This requires the appearance of the constant α 0 , which results from a different form of the coefficients in P app (n t ) and the constant coefficients α i , i = 1, 2, . . ., k, in P f it (char; α j | k j=0 ).In addition, the residual standard deviation √ 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 α 0 , depend only implicitly on the transition density n t and in the limit n t → 0, α 0 → 0.
The regression function for the AIC and R 2 adj selected regression model in the fourthorder approximation has the form The above regression model is also confirmed by the backward analysis procedure applied to the set of factors (22) as all values of the estimators of the structural parameters α j , j = 0, 1, 2, . . ., k = 6 of the regression model are statistically significant at the level α = 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 the regression function given in Equation ( 28) are given in Table 6.The most probable transition density n t = ñt was determined by solving Equation (27).The RHS of this equation is the appropriate regression model with a specific number of factors.At the same time, the LHS in the case when the fourth-order contribution to the symmetry energy is included, following Equation (11), is expressed by elements characterizing the dependence of nuclear matter on density.Since the coefficients of the factors describing nuclear matter depend on the transition density, solving equation Equation( 27) makes determining the transition density value possible.For example, the form of Equation ( 27) for the AIC-selected regression model is presented as Solving the above equation for n t allows one to calculate its value for the selected regression model.This procedure was carried out for each of the 23 RMF models, and then the average value ñt was calculated from the obtained 23 ñt values.Similar calculations were performed for the parabolic approximation of the symmetry energy E sym,2 (n b ) and the case of δ = 1.The regression model selected by the AIC gives the most probable appearance of the sample [9].Equation ( 27) is a relationship imposed on the model specified by the AIC that guaranteed that the observed sample appeared with maximum probability for a fixed number of factors.Thus, the value of n t = ñt determined from Equation ( 27) is the most probable value for the determined number of factors selected by the AIC regression model.Because the AIC's globally selected regression model is the best estimate of the true regression model, the transition density ñt resulting from the performed regression analysis is considered the best approximation of the transition density implied by the true regression model.As a consequence, ñt should characterize the true physical model.
In the parabolic case, the regression model selected via the AIC and R 2 adj has the following regression function: This regression model is confirmed through the backward analysis procedure applied to the set of factors (21) at the significance level α = 0.065.The values of the estimators of the structural parameters other than α 1 for the factor K 0 and α 4 for the factor L 2, Yδ 2 , remain in the model at a significance level lower than α = 0.05.The other characteristics of the selected regression model with the regression function given in Equation ( 30) are given in Table 6.
To calculate the uncertainty of the estimation of a particular value of ñt , two components, namely the error of the estimation of the conditional expectation value P f it (char; α j | k j=0 ) (which appeared to be decisive) and the error propagation from P f it (char; α j | k j=0 ) to ñt , were calculated.The obtained uncertainty of the estimation of a particular ñt coming from these two sources is, on average, approximately ±0.024 when the fourth-order symmetry energy term is included in the analysis (k = 6) and ±0.03 for the parabolic case (k = 5).
Table 5.The characteristics of the regression models selected by the AIC and R 2 adj with the regression functions (30) in the parabolic case and (28) in the fourth-order case.SSR, SSE, and SSY are the sum of squares due to regression, the error sum of squares, and the total sum of squares of the response Y ≡ P t , respectively, and SSY = SSR + SSE.MSE (which is the variance of the error term Ê) is the mean squared error (Appendix A.1), [10].σα 0 to σα 5 are the standard errors of α0 to α5 in the parabolic approximation case, and σα 0 to σα 6 are the standard errors of α0 to α6 in the fourth-order approximation case.Figure 3 shows the values of the means ñt of the most probable density values ñt (see Tables 1-4), obtained for δ = δ t and δ = 1 for the two considered cases of symmetry energy approximations as a function of the number of explanatory variables n that characterize a given regression model selected using the AIC and R 2 adj methods.The crucial relation for the further construction of a true physical model is ñt (L sym,2 ).This relation for the model selected by the AIC and R 2 adj in the fourth-order approximation with the regression function ( 28) is shown in Figure 4. Green squares illustrate the n t values calculated for individual RMF models.The n t values are shown in Table 7.The AIC method, searching for phenomena related to the location of the neutron star's crust-core transition described by the RMF models shifts the crust-core boundary to higher densities.In the paper [14], the neutron star's core-crust transition densities obtained within the dynamical and thermodynamical methods using the full EoS and its PA with the MDI and Skyrme interactions have been analyzed.It should be emphasized that the most probable transition density ñt values determined in this paper, based on the proposed probabilistic method for the symmetry energy supplemented by the fourth-order contribution, given as a function of the symmetry energy slope L sym,2 , follow a very similar course as in the case of models with modified Gogny (MDI) interaction and 51 Skyrme interactions.The results reported in other papers show that the transition density n t decreases as L increases, as it is anticorrelated with L. This has been verified with many different models [46][47][48][49].Table 6.Numerical values of the transition pressure P t calculated for individual models for δ = δ t , i.e., for the value resulting from Equations ( 5) and ( 6) and for neutron matter (δ = 1).The table contains pressure values for the parabolic approximation and the case when a fourth-order contribution is included in the description of the symmetry energy.The subscript 2 refers to the quantities calculated based on the parabolic approximation, and the subscript 24 indicates the sum of the second-and fourth-order contributions.The pressure P t is given in Mev/fm 3 .

Model
. The linear and exponential fits to the sample of (L sym,2 ; ñt ) points for the AIC selected model for the fourth-order case with the regression function given by Equation ( 28).The most probable values of ñt (black circles) are calculated from Equations ( 27) and ( 28) for the values of x ≡ L sym,2 for the sample of N = 23 RMF models.As for the mean square errors, MSE(exponential) = 2.683 × 10 −5 < MSE(linear) = 4.5767 × 10 −5 ( f m −6 ); thus, the exponential fit is better [50].Green squares illustrate the n t values calculated for individual RMF models.
Table 7. Numerical values of the transition density n t calculated for individual RMF models for δ = δ t , i.e., for the value of δ resulting from Equations ( 5) and ( 6).The case of pure neutron matter (δ = 1) is also included.The table contains transition density values for the parabolic approximation and the case when the fourth-order contribution in the description of the symmetry energy is included.The subscript 2 refers to the quantities calculated based on the parabolic approximation, and the subscript 24 indicates the sum of the second-and fourth-order contributions.The transition density n t is given in fm −3 .

Conclusions
A regression model selected according to the AIC and R 2 adj model-selection procedures can be used to identify and support a correct model from the physical and experimental points of view.The results obtained include, among others things, the importance of approximations used to describe the symmetry energy.Suppose the isospin asymmetry parameter has the value resulting from the adopted model δ = δ t .In this case, the selected regression models corresponding to the extreme values of AIC and R 2 adj , both for the parabolic approximation and considering the fourth-order term, include as explanatory variables K 0 and L sym,2 , K sym,2 , and E sym,2 multiplied by functions of δ t (see (21) and ( 22)).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 δ 3  24 ).A different result is obtained assuming the regression analysis is performed 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 (n b ) + E sym,4 (n b ), 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 .Regardless of the considered values of δ, the selected regression models for the parabolic approximation 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 ñt = 0.05876 ± 0.00572 fm −3 .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 model with AIC and R 2 adj selected (with the regression function (28)), it is possible to examine the relationship between ñ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.

δ 2 Parabolic with δ 2 Figure 1 .
Figure 1.The minimal values of the AIC for a given number of explanatory variables in the regression model.The regression model with a globally minimal AIC gives the highest probability of the appearance of the sample of RMF points.The lines connecting the symbols are a guide for the eyes only.

δ 2 Parabolic with δ 2 Figure 2 .
Figure 2. The adjusted coefficient of determination R 2 adj vs the numbers n of the explanatory variables used in a given regression model.The lines connecting the symbols are a guide for the eyes only.

3 )Parabolic with δ 2 Parabolic with δ 2 Figure 3 .
Figure 3.The most probable transition density ñt vs. the numbers n of the explanatory variables used in a given regression model selected by AIC.The maximal value of ñt refers to the global minimal AIC.The lines connecting the symbols are a guide for the eyes only.

Table 1 .
(27) characteristics of the regression models in the parabolic approximation case with δ 2 .R 2 is the coefficient of determination, R 2 adj is the adjusted coefficient of determination (Appendix A.1), AIC is the Akaike information criterion given in Equation(19), and the ñt values in the table are the means of the most probable density values assuming Equation(27).The most likely value is for the model with a globally minimal AIC given in boldface characters.Yδ 2 L 2, Yδ 2 L 2, δ 2 Yδ 2 L 2, Yδ 2 L 2, δ ( K 0 E 2,

Table 2 .
(27) characteristics of the regression models when the fourth-order contribution is included with δ 24 .R 2 is the coefficient of determination, R 2 adj is the adjusted coefficient of determination (Appendix A.1), AIC is the Akaike information criterion given in Equation(19), and ñt values in the table are the means of the most probable density values assuming Equation(27).The most likely value is for the model with a globally minimal AIC given in boldface characters.Yδ 24 L 2, Yδ 24 K 4, δ

Table 3 .
(21)case when the parabolic approximation gives the symmetry energy.The regression models are determined for δ 2 = 1.ñt in the table are the means of the most probable density values assuming Equation(27).The most likely value of ñt is the one obtained for the model with a globally minimal AIC value.The variables in this table are from set(21)in the case of δ 2 = 1.The most likely value is for the model with a globally minimal AIC given in boldface characters.

Table 4 .
(22)case when the symmetry energy is represented by the functions E sym,2 (n b ) + E sym,4 (n b ) and for δ 24 = 1.ñt denotes the means of the most probable density values assuming Equation(27).The most likely value is for the model with a globally minimal AIC value.The variables in this table are the ones from set(22)in the case of δ 24 = 1.The most likely value is for the model with a globally minimal AIC given in boldface characters.