An approach to constrain models of accreting neutron stars with the use of an equation of state

We investigate X-ray bursts during the thermal evolution of an accreting neutron star which corresponds to the X-ray burster GS\ 1826-24. Physical quantities of the neutron star are included using an equation of state below and above the nuclear matter density. We adopt an equation of state and construct an approximate network that saves the computational time and calculates nuclear energy generation rates accompanying the abundance evolutions. The mass and radius of the neutron star are got by solving the stellar evolution equations from the center to the surface which involve necessary information such as the nuclear energy generation in accreting layers, heating from the crust, and neutrino emissions inside the stellar core. We reproduce the light curve and recurrence time of the X-ray burst from GS 1826-24 within the standard deviation of 1$\sigma$ for the assumed accretion rate, metallicity, and equation of state. It is concluded that the observed recurrence time is consistent with the theoretical model having metallicity of the initial CNO elements $Z_{\rm CNO}$ = 0.01. We suggest that the nuclear reaction rates responsible for the $rp$-process should be examined in detail, because the rates may change the shape of the light curve and our conclusion.


Introduction
The type I X-ray burst was discovered from the Low Mass X-ray Binary in 1975 [1] and 112 bursters have been observed until now [2]. These observations provide firm ground to the arguments that bursts are identified to be a phenomenon associated with shell flashes on the accreting neutron stars in close binary systems (e.g. [3]). The shell flashes are now believed to be initiated with thermonuclear runaway as shown in Fig. 1 [4]. The X-ray burster GS 1826-24 has been observed in 1988 by X-ray astronomical satellite Ginga [5]. Since an accretion rate for an X-ray burst to start generally varies even from the same burster, the shape of the light curve changes significantly [6]. However, the accretion rate (dM/dt) of GS 1826-24 has been known to be almost constant during a year and therefore the shape of the light curve is nearly invariant [7]. From the theory of the X-ray burst, it is convenient to make models of light curves under the constant dM/dt [8,9]. As a consequence, it is believed to be appropriate to check the validity of models [10,11]. Therefore, we can use the X-ray burster GS 1826-24 to compare numerical results of light curves calculated from 3α Hot CNO cycle Breakout reaction 12 19 Ne accelerated by α(2α, γ) 12 C reactions [4]. Due to the increasing temperature, the rapid proton capture process (rp-process) ensures by overcoming the Coulomb barrier.
the assumed dM/dt with the observational data. In the present paper, we adopt the data of RXTE reported in Ref. [12] to construct models. In the meanwhile, it was investigated to reproduce observational results by analyzing numerical results concerning X-ray bursts of GS 1826-24 [10,11,13,14]. In Ref. [10], a large nuclear reaction network which includes 1300 nuclides was used to follow the detailed nucleosynthesis, which is based on the numerical method by Woosley et al. [15]. None the less, the numerical method does not include an equation of state at the high density above nuclear matter which relates to light curves. This is because their calculations have been limited to around the accretion layers above the crust, where the neutron star mass and radius are selected as parameters to rebuild light curves. Furthermore, artificial luminosity due to shallow and/or crustal heating is given at the bottom of the calculated region which may locate above the outer crust of the neutron star [16]. Generally speaking, it is difficult to constrain physical processes occurred inside the neutron stars as far as only accretion layers responsible for bursts are taken account of [17]. In this connection, the ignition properties have been investigated from the point of "mixed bursts" [18], where impact on the burst ignition related to stable and/or unstable nuclear burnings has been examined with the use of a linear stability analysis; The discrepancies between the results and other multi-zone calculations may be ascribed to the strength of the Hot CNO cycle [19]. Since many reaction rates responsible for the rp-process are still very uncertain, the suggested problem remains unresolved. On the other hand, effects of the inner core below the neutron star crust are not clear due to the uncertainty of nuclear physics at the high density, where an equation of state (EoS) is still very uncertain. For example, cooling of neutron stars [20,21] or supernova explosions [22] have been studied taking account of some EoSs having significant uncertainty.

2/11
In the present study, to perform the calculation of X-ray bursts, the general relativistic equations for the stellar structure and evolution are solved exactly from the center to the outermost layer, where we can study thermal interactions between the accreted layers and the core of the neutron star [17]. Besides, we develop an approximate network that includes 88 nuclides; it reproduces the nuclear energy generation and rapid proton capture process (rpprocess) during the burst. By including necessary physical inputs such as an equation of state, nuclear burning, opacities, and neutrino losses, we propose some models with the mass and radius of the neutron star responsible for the X-ray burster GS 1826-24. In section 2, basic equations of neutron star evolution are given with physical inputs. Our approximate network to follow the nuclear processes is explained in section 3. The results of the calculations are compared with the observations using the statistical analysis in section 4. We discuss our results toward further understanding of X-ray bursts in section 5.

Basic equations and physical inputs for accreting neutron stars
We have performed numerical calculations of the thermal evolution of neutron stars in hydrostatic equilibrium by using the spherical symmetric stellar evolutionary code [17], which includes full general relativistic effects formulated by Throne [23]. The basic simultaneous differential equations are written as follows: Here, M tr and M r in a radius r are gravitational and rest masses; ρ and ρ 0 denote the total mass-energy and rest mass densities; P and T are the pressure and local temperature; ε n and ε g are the energy generation rates by nuclear burning and gravitational energy release, respectively. Besides, ε ν represents the energy loss rate by neutrino emission; ∇ rad and ∇ ad are the radiative and adiabatic gradients; φ is the gravitational potential in unit mass. We adopt EoS by Lattimer and Swesty [24] with the incompressibility K of 180 or 220 MeV in the inner layers (ρ ≥ 10 12.8 g cm −3 ) and connect it to EoS in Refs. [17,28] for the outer layers (ρ < 10 12.8 g cm −3 ). It is noted that the EoS we have adopted is still applied to study the neutron star properties above the nuclear density [21,22], because it may be consistent with both experimental and astrophysical constraints compared to other sophisticated EoS [20]. Opacities appropriate for the neutron star evolution are the same as in Ref. [27]. The neutrino emission process is taken into account as a slow cooling process: electron-positron pair, photo, and plasmon processes [29,30]; bremsstrahlung process; modified Urca (MURCA) process. Dominant processes are MURCA [31] and bremsstrahlung [29]. The corresponding energy loss rates are summarized in Ref. [30]: Although these neutrino emission rates have been studied, there remains the uncertainty of a factor of ten because of insufficient understanding of the symmetry energy and nucleon effective mass in the dense matter [32]. Moreover, we do not include a strong process such as pion condensation. While pion condensation accompanies strong neutrino emissions, the effects of the super-fluidity may reduce the neutrino emissions below the critical temperature T cr in proportion to exp(−aT cr /T ) with a constant a. Since the present purpose is to show our approach to constrain the properties coming from the neutron star core using X-ray burst observations for high accretion rates dM/dt > 10 −9 M ⊙ yr −1 , we neglect the strong neutrino emission process [17]. We construct initial models of a neutron star accreted at the constant rates dM/dt = (2 − 4) × 10 −9 M ⊙ yr −1 with mass fractions, X( 1 H) = 0.73, X( 4 He) = 0.25, X( 14 O) = 0.007, and X( 15 O) = 0.013, where X(O 14 ) and X(O 15 ) are assumed to be equilibrium values of Hot CNO cycle (see Fig. 1). The initial model is in a steady-state, in which the nonhomologous part of the gravitational energy release can be neglected [17], where nuclear burning is temporally switched off. As a result, we can get initial models forṀ −9 = 2, 2.5, 3, and 4, whereṀ −9 is the accretion rate in units of 10 −9 M ⊙ yr −1 . Corresponding to the mass-radius relation against K, the gravitational mass and radius are obtained as follows: (M/M ⊙ , R(km)) = (1.57, 11.9) for K = 180 MeV, (1.58, 12.6) and (2.0, 11.3) for K = 220 MeV. The mass of 2 M ⊙ may be rather high for the low mass X-ray binary; we adopt this model for reference as shown in Table 2.

Approximate network
We have constructed an approximate network that nearly reproduces the nuclear energy generation rates and nuclear abundances related to X-ray bursts [4,33]. The small network was developed first by Wallace and Woosley [4] with 10 nuclei up to 56 Ni and then by Hanawa  et al. [33] with the inclusion of 16 nuclei up to 68 Se. The construction of these approximate networks is based on the calculations of the nucleosynthesis by large networks. We have made a revised version of the approximate network (hereafter APRX3) which includes 88 nuclei as shown in Table 1. To construct APRX3, we used a large network [34], where the updated nuclear reaction rates have been adopted from Ref. [35] except for some data: Two reactions 64 Ge(p, γ) 65 As and 65 As(p, γ) 66 Se and their reverse rates are taken from Ref. [36]. The β + rates are from Refs. [37][38][39][40]. APRX3 is built to reproduce the nuclear reaction paths examined by Schatz [41]. In particular, it is pointed out that the rp-process proceeds to the formation of nuclei of Z < 50 and the nucleosynthesis is closed down at the Sn-Sb-Te cycle: 103 Sn(β) 103 In(p, γ) 104 Sn(β) 104 In(p, γ) 105 Sn(p, γ) 106 Sb(p, γ) 107 Te(γ, α) 103 Sn.
Our large nuclear reaction network contains 897 nuclei up to 130 Sm (hereafter FNRN) and is adequate to study X-ray bursts [34]. To simulate bursts we assume one zone model of the ignition pressure to be log P = 22.8 for the neutron star with M = 1.4M ⊙ and R = 10 km [33,34]. The initial mass fractions are the same as described in the previous section. The networks reproduce changes in hydrogen abundance and nuclear energy generation rate as shown in Figs. 2 and 3, respectively, where APRX2 is another approximate network that contains 61nuclei [42]. It can be seen that X 1 (H) calculated by APRX2 deviates significantly from that of FNRN at t > 850 s. APRX3 can reproduce both X( 1 H) and ε n of FNRN with enough accuracy. Therefore, we adopt APRX3 in the present study. In the above calculations, the computational time of FNRN is 6.3 times longer than that of APRX3. We note that our parallel computer machine has the power of fourteen cores per cpu.

Comparison with observations
The observed luminosity L obs b can be expressed in terms of the observed flux F obs b during a burst and the distance d to GS 1826-24 as follows: where ξ b is a factor of anisotropy for the burst; if the radiation is isotropic, ξ b = 1. It is noted that the calculated peak luminosity L pk divided by (1 + z) 2 corresponds to L obs b . On the other hand, d 2 ξ b is not obtained from observations. Therefore, L obs b cannot be compared to the calculated value of L thr b . In our case, we define χ 2 at a given trial value of dξ Here, t pk is the time from the beginning to the peak of the luminosity, n obs is the total number of observations consisting of the light curve (n obs = 10) [12], L obs b,i (dξ 1/2 b , t pk ) is the luminosity under the assumption of two assumed parameters of dξ 1/2 b and t pk , and L thr b is the theoretical luminosity calculated from our models. Two values of the standard deviations, σ i,obs and σ i,thr are those of observation and calculation, respectively. We note that σ i,thr is taken to be the deviation from the average luminosity over the total number of the successive luminosities n burst . The minimum quantity of χ 2 is determined by changing both dξ 1/2 b and t pk .
In Table 2, the results of calculations are summarized for 27 models, where we can see the effects of EoS, the accretion rate, and the metallicity on the light curve. Models LS180 and LS220 adopt EoS [24] with the respective incompressibility K of 180 and 220 MeV; the accretion rateṀ −9 =2, 2.5, 3, and 4. The metallicity of CNO elements is selected to be Z CNO = 0.005, 0.01, and 0.02. Eight light curves are shown for models L2n20Z1 and L2n20Z05 in Fig. 4. We note that the recurrence time ∆t depends onṀ , M , and Z CNO . In Fig. 5, it is seen that the light curves depend significantly on the metallicity, which constrains the model; Both the height of the peak luminosity and the shape of the light curve are different in obedience to Z CNO . We note that it is difficult to fit both the height of the light curve and the recurrence time in case of M = 2 M ⊙ .
For comparison, we show the results in the table with the use of K = 180 MeV whose value is incapable of reproducing the experimental one of K [25,26], where the calculated luminosity is lower than that of the observed one. Figure 6 designates the recurrence time between bursts against the accretion rate for the calculated models. It depends significantly on the metallicity. As shown in Table 2 we obtained the recurrence time ∆t = 3.56 ± 0.07 hr for L2n20Z1 which is in the range of the observed value of ∆t = 3.536 ± 0.04 hr [12]. In addition to the recurrence time, it is clear that the overall shape of the light curve for the lower-mass model (L2n20z1) is much better than the higher mass model such as L2n20Z1 as shown in Fig. 7.
Since other models cannot fit the peak value of the light curve observation, among models investigated in the present paper, the best fit model which reproduces the observation of the light curve in 2007 [12] becomes L2n20Z1, where we can see the degree of agreement between calculations with 1σ band and observations. We remark that while the observational average of the light curve was obtained from 10 bursts (n burst = 10), we select 21 successive bursts (n burst = 21) to get average luminosity. We may eliminate numerical ambiguity by including more bursts compared to actual observations. It should be noted that the model of

Discussion
We have tried to reproduce the light curve observation of X-ray bursts from GS 1826-24. To properly calculate the X-ray bursts, it is needed to include physical quantities such as EoS, nuclear energy generation rates, opacities, and neutrino loss rates inherent to accreting neutron stars. Zamfir et al. [16] constrained neutron star mass and radius to be R < 9.0 − 13.2 km and M < 1.2 − 1.7 M ⊙ , respectively. However, they used a theoretical light curve model in Ref. [10] who assumed the radius R = 11.2 km of M = 1.4 M ⊙ and the gravitational redshift is z = 0.26. In this paper, we focus on specific EoS which determines the mass and radius. As a result, by using the statistical method of χ 2 , we reproduce the observational light curve with the standard deviation of 1σ level. If we adopt another EoS, we may get different results. To see the dependence on EoS, we have shown the case of K = 180 MeV in Table 2 , where the peak luminosity is rather low and the recurrence time becomes shorter. The incompressibility K is the most uncertain quantity in EoS and the softer one  seems to be inconsistent with observations. It should be continued to study the effects of EoS properties on X-ray bursts. Concerning the burster GS 1826-24, the accretion rate is inferred to be around 10 −9 M ⊙ yr −1 which is a rather high accretion rate to satisfy the condition of regular bursts. In this case, X-ray bursts will depend on the accretion layers and the connection between the accretion layers and the stellar core is weak [17]. By comparing the theoretical and observational light curves, some physical processes such as a strong cooling process should be examined based on theories and experiments of nuclear physics [27]. From Fig. 5, the shape of the theoretical light curve depends on Z CNO . If we choseṀ −9 = 2.0, Z CNO = 0.01 and M = 1.58 M ⊙ , the tail of the light curve agrees rather well compared to the observation as seen in Fig. 7. This is due to the energy generation of the rp-process. On the other hand, many reaction rates concerning the nucleosynthesis will be responsible 8/11 Table 2 Physical quantities for 27 models. We adopt the observational data in 2007 [12]. Numerical value inside the bracket indicates the standard deviation 1σ.

MODEL
EoS for the shape. The uncertainty of the reaction rates has been guessed in the range of factor 100, because of the poor understanding of the reaction rates in the proton-rich nuclei [13,15]. As a consequence, it is difficult to select key reactions important for the construction of the light curve [14]. Finally, we consider the ratio of the burst anisotropy ξ b to the persistent one ξ p between bursts. The anisotropy related to the persistent luminosity is obtained fromṀ and the gravitational redshift z g : 4πd 2 ξ p F x =Ṁ c 2 z g /(1 + z g ) with the speed of light c [10]. Our best fit model L2n20Z1 gives z g = 0.26,Ṁ = 2 × 10 −9 M ⊙ yr −1 , and the persistent flux F x = 5.4 × 10 −9 erg cm −2 s −1 . With the use of the peak flux and luminosity, d 2 ξ b is calculated from F peak = L peak /4πd 2 ξ b . As the result, the ratio ξ b /ξ p = 1.01. Furthermore, the inclination angle relative to the observer becomes θ ≈ 58.6 • on the assumption of a thin flat disk [43]. 9/11 It has been reported that while the angle becomes θ ≈ 80 • using the MESA code [13], the KEPLER best fit requires θ ≈ 65 • [10]. Although it is complicated to constrain the disk model, it may be worthwhile to study the dependence of the anisotropy on EoS.
We remark future studies concerning type I X-ray bursts which should be carried out based on the present work.
1. Detailed constraints for the equation of state beyond the nuclear saturation density will become possible by comparing light curves between theories and observations, which include the uncertainty of disk models.
2. Whole reaction rates related to the rp-process should be restudied from the point of nuclear experiments and theories.
3. Bursts with strong cooling coupled to parameters associated with the superfluidity model would be examined for lower accretion rates dM/dt ≤ 10 −9 M ⊙ yr −1 and/or long intervals between bursts.
4. Long term simulations of X-ray bursts beyond one billion years could be desirable by using the approximate network since the beginning of bursts after the formation of low mass X-ray binary cannot be identified.