The Impacts of Neutron-Star Structure and Base Heating on Type I X-Ray Bursts and Code Comparison

Type I X-ray bursts are rapidly brightening phenomena triggered by thermonuclear burning on accreting layer of a neutron star (NS). The light curves represent the physical properties of NSs and the nuclear reactions on the proton-rich nuclei. The numerical treatments of the accreting NS and physics of the NS interior are not established, which shows uncertainty in modelling for observed X-ray light curves. In this study, we investigate theoretical X-ray-burst models, compared with burst light curves with GS~1826-24 observations. We focus on the impacts of the NS mass, the NS radius, and base-heating on the NS surface using the MESA code. We find a monotonic correlation between the NS mass and the parameters of the light curve. The higher the mass, the longer the recurrence time and the greater the peak luminosity. While the larger the radius, the longer the recurrence time, the peak luminosity remains nearly constant. In the case of increasing base heating, both the recurrence time and peak luminosity decrease. We also examine the above results using with a different numerical code, HERES, based on general relativity and consider the central NS. We find that the burst rate, burst energy and burst strength are almost same in two X-ray burst codes by adjusting the base-heat parameter in MESA (the relative errors $\lesssim5\%$), while the duration time and the rise time are significantly different between (the relative error is possibly $\sim50\%$). The peak luminosity and the e-folding time are ragged between two codes for different accretion rates.


INTRODUCTION
Type I X-ray bursts are periodic eruptions caused by unstable thermonuclear burning on the surface of neutron stars (NSs) in low-mass X-ray binary (LMXB) systems (Joss 1977;Parikh et al. 2013). The NS accreted matter from the companion star, overflowed through the Roche lobe, and formed an envelope on the surface of the neutron star. Under the action of gravity, the accreted matter was continuously compressed and heated, thereby increasing the temperature and density, when the energy generation rate is greater than the cooling rate, thermonuclear unstable combustion will occur, resulting in type I X-ray bursts (Woosley & Taam 1976;Lewin et al. 1993;Bildsten 2000;Galloway & Keek 2021). The accreted matter mainly provides energy for type I X-ray bursts through 3α reaction, CNO cycle and rp-process, etc (Wallace & Woosley 1981;Taam 1985;Bildsten 1998;Galloway et al. 2008). Burning produces a heavy accumulation of ash, and as new material continues to pile on top of it, the accreted material undergoes gradual compression until it reaches a condition for ignition, producing another burst sequence.
Since the first discovery of the X-ray burst in 1975 (Belian et al. 1976;Grindlay et al. 1976), more than 7000 events from 118 1 bursting sources have been observed so far (Galloway et al. 2020). By comparing with observations, theoretical models can be calibrated and the physical properties of NSs can be constrained. (Cromartie et al. 2020). GS1826-24 is one of the preferred sources because of a nearly uniform accretion rate and regular burst behaviour, which is called "clock burst" or "textbook" burst (Ubertini et al. 1999;Bildsten 2000). X-ray burst models require input parameters regarding the accreted fuel composition (X,Y,Z), mass accretion rate (Ṁ ), base heating (Q b ), mass (M ) and radius (R) of NS, as well as nuclear reaction rates. Heger et al. (2007) studied the effect of metallicity (Z) and mass accretion rate on the theoretical light curves, by the comparison with the light curve of GS1826-24, they estimated the initial metallicity and the accretion rate of GS1826-24. Meisel (2018 investigated the sensitivity of models to varied accretion rate, base heating, metallicity and the nuclear reaction rate 15 O(α, γ) 19 Ne, by modelobservation comparisons, they constrained the shallow heating in GS 1826-24 should be below 0.5 MeV/u. The influence of nuclear reaction rate uncertainties on NS properties also has been studied from X-ray burst model-observation comparisons (Meisel et al. 2019).
The above X-ray burst simulations are based on KE-PLER (Heger et al. 2007) or MESA (Meisel 2018;Meisel et al. 2019), which consider the NS envelope using inner boundary conditions with fixed NS mass and radius (1.4M and 11.2 km). The effects of the NS mass and radius on thermonuclear flashes are investigated by Joss & Li (1980) and Ayasli & Joss (1982) using the stellar evolution code ASTRA. They adopted M = 1.4M and R = 6.57 km as a standard case, and vary mass to M = 0.705M as the low mass case and radius to R = 13.14 km as the large radius case. The results show that the recurrence time, accumulated mass, burst energy, burst strength and peak luminosity have obvious change. However, the results are not consistent with the recent NS mass-radius constraint (Steiner et al. 2010;Abbott et al. 2018).
Recently, Dohi et al. (2020Dohi et al. ( , 2021Dohi et al. ( , 2022) studied X-ray bursts using a general relativistic stellar-evolution code 1 https://personal.sron.nl/∼jeanz/bursterlist.html with several NS equation of states (EOSs). They focused on the microphysics inside NSs (e.g., the mass and radius with different EOSs and the NS cooling process). By comparing with the burst parameters of GS 1826-24, they constrained the EOS and the NS mass and raidus. Meanwhile, Johnston et al. (2020) apply Markov chain Monte Carlo methods to 3840 Kepler X-ray burst models and obtain system parameter estimates for GS 1826-24. They estimate a metallicity of Z CNO = 0.010 +0.005 −0.004 , hydrogen fraction of X 0 = 0.74 +0.02 −0.03 , mass M > 1.7M , radius R = 11.3 +1.3 −1.3 , etc. So far, the NS mass and radius are unknown for burst sources, but mass and radius change the burst properties. It is worth for us to extract the information on the macroscopic properties of NS from the observation of X-ray bursts.
As X-ray burst simulations with a general relativistic stellar-evolution code solve the stellar evolution equations from the center to the surface with the EOS, neutrino emission, crust heating as well as the nuclear energy generation in accreting layers are important for the comparison to X-ray burst observations (Dohi et al. 2020(Dohi et al. , 2022. The MESA and KEPLER code only consider the accreting layers above NS solid crust, where base heating parameter Q b is adopted at the inner boundary to mimic the energy transfer from the NS interior. However, the value of base heating is not well constrained by observation. Keek & Heger (2016) assumed Q b = 0.1 MeV u −1 , the deep crust heating theory suggests that the generated heat may be larger, up to Q b = 2 MeV u −1 (Haensel & Zdunik 1990, 2003, although most of the heating in the deep crust is conducted into the core and carried off by neutrinos, a considerable amount of local heating will occur, which may increase Q b . A yet-unknown shallow heating may also increase Q b (Brown & Cumming 2009;Deibel et al. 2015;Lu et al. 2022). On the other hand, Q b may be reduced by the competing effect of neutrino cooling (Cumming et al. 2006), the Urca neutrino cooling process in the outer crust may also complicate the estimation of Q b (Schatz et al. 2014). Thus it is significant for us to study the effect of Q b on X-ray bursts, by model-observation comparison, we may give a constraint on its value.
In addition, in Newtonian codes such as MESA or KEPLER, to accurately model bursts, it is important to account for the General Relativity (GR) effects when comparing models with observations. The MESA code adopts the Post-Newtonian correction to include the effects of GR (Paxton et al. 2015;Meisel 2018). The KE-PLER code uses Newtonian gravity and ignores the GR effects, thus GR corrections are adopted in X-ray burst simulations (Keek & Heger 2011;Johnston et al. 2018Johnston et al. , 2020 Figure 1. A schematic picture of the NS structure. The computational domain of several X-ray burst codes for the thermal evolution of accreting NS is shown. The label with the asterisk (*) considers the effects of convection and nuclear reaction networks, and with (s) treats the envelope as in the steady state.
simulate a sequence of X-ray bursts, and compare the results from a general relativistic stellar-evolution code, HERES , focusing on the burst observables. The structure of the paper is as follows. In Section 2, we describe the Post-Newtonian hydrodynamic MESA model and the GR hydrostatic HERES model. In Section 3, we present the results of computations wherein the effects of masses, radii and base heating upon the X-ray burst properties are taken into account, and we compare the results from the MESA code and HERES code  in Section 4. Finally, we summarize our results and briefly discuss their implications.

MODEL
There are several X-ray burst models which have different features due to the nature of burst code. In Fig. 1, we show the schematic of NS structure and several corresponding burst codes. In several stellar evolutionary models, the most used code is the MESA code (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018, which solves the (post-)Newtonian hydrodynamics within the accreted regions. The formulation of MESA is quite similar to some codes of KEPLER (Woosley et al. 2004) and SHIVA (José & Hernanz 1998;José et al. 2010). These codes can be accessible even with the thermal evolution of relativistic compact objects by using the "GR correction" (e.g., Keek & Heger (2011)) 2 , but since the boundary condition on the crust surface is inevitably introduced as "Q b " value, it is hard for them to probe the NS physics. The approximate treatment of strong gravity in NSs as above may not be valid except for the surface, and therefore consistent treatment based on general relativistic formulation is indispensable for more exact calculation of burst light curves.
The sophisticated public code which takes into account the above is the dSTAR code (Brown 2015) originally developed by Brown (2000). dSTAR simultaneously solves the TOV and energy transport equations without convection and reaction networks for X-ray bursts. It covers the regions except for the NS core and can probe the crust physics such as crustal heating, shallow heating, Urca cooling, and so on (Deibel et al. 2015(Deibel et al. , 2016Meisel & Deibel 2017). For the envelope, it treats as in the steady state, which can construct a relation between surface temperature and the crust temperature at the shallowest point (Brown et al. 2002). Still, dSTAR leaves the boundary condition on the core surface, which must be changed by the NS physics such as the EOS and ν cooling effects.
To include more possible NS physics, we have recently developed the code of HERES 3 (Dohi et al. 2020). HERES is essentially the same as dSTAR in that both follow quasi-thermal evolution of accreting NSs, but the covered regions for calculation are extended to the center of the NSs. Unlike the other codes, no artificial boundary condition such as Q b is required. Another two codes of NSCool (Colpi et al. 2001;Page & Reddy 2013) 4 and PC18 (Potekhin & Chabrier 2018) 5 (see also Potekhin & Chabrier (2021)) are similar to HERES in regard with the formulation but without convection and reaction network for X-ray bursts. Therefore, HERES is currently the unique code that can probe the NS physics from X-ray burst light curves. In this work, we adopt two distinct codes of MESA and HERES. Next, we briefly explain the properties of each code.

Post-Newtonian Hydrodynamic model with large reaction network (MESA)
We use an open-source stellar evolution code (MESA, version 9793) (Paxton et al. 2015) to perform calculations on Type I X-ray bursts. The MESA equation of state (EOS) is based on the 2005 OPAL EOS tables (Rogers & Nayfonov 2002), besides, the SCVH tables (Saumon et al. 1995), HELM (Timmes & Swesty 2000) and PC (Potekhin & Chabrier 2010) EOSs are employed for various conditions (Paxton et al. 2011). It is worth to mention that a new Skye EOS for fully ionized is designed by Jermyn et al. (2021), which has been tested in action in the MESA stellar evolution code by computing white dwarf cooling curves. OPAL opacity tables are used with the proto-solar abundances from Asplund et al. (2009). Following the approach described in Meisel (2018), we model a series of NS envelopes by considering inner boundary conditions for different NS masses and radii. The most pertinent details are repeated here. The luminosity at the base of the envelope is set to L base = M Q b , where Q b is the base heat, a parameter adopted by many models to simulate the heat flow from the neutron star's crust into the envelope (Brown & Cumming 2009;Galloway & Keek 2021;Keek & Heger 2017). The change of mass, radius and base luminosity by using the commands "relax M center ", "relax R center " and "relax L center ", respectively (Paxton et al. 2011). GR effects were accounted for using a post-Newtonian modification to the local gravity (Paxton et al. 2011(Paxton et al. , 2015, where the MESA setting "use GR factors = .true." was chosen. The envelope thickness is approximately 0.01 km, and the initial metal abundance uses the solar metal abundance Z = 0.01, 0.02 (Grevesse & Sauval 1998). We use the rp.net, which contains 304 isotopes (see Lund Fisker et al. (2007)), and nuclear reaction rates use the reaction rates from the REACLIB V2.2 library (Cyburt et al. 2010). Adaptive time and spatial resolution were employed according to the MESA controls varcontrol target=1d-3 and mesh delta coeff=1.0 (Paxton et al. 2013). In order to achieve convergent solutions, some models need slightly different settings. In table 1 from the Appendix, we provide our burst models, which describe the input parameters and some of the outputs in more detail.

General-Relativistic hydrostatic Evolutional model
with an approximate reaction network (HERES) As explained above, MESA has two issues of the treatment of NS gravity and artificial boundary condition introduced as Q b . As the definitions of Q b are different in previous works, e.g., Q b is defined by Keek & Heger (2016) to mean the amount of heat generated by crustal heating at the base of the envelope, and the typical value for Q b of 0.1 MeV/u was adopted. Q b is defined by Meisel (2018) to mean not only the crust heating but also the shallow heating. Q b = 0.1, 0.5, 1.0 MeV/u were adopted to mimic the shallow heating of unknown origin. Hereafter, we define the net base heat as Q e , which represents the energy exchange between the interior NS and the accreting layer, its value could be changed by the unknown shallow heating or the ν cooling processes inside NSs, related to the EOS and mass (see Table 2 in Dohi et al. (2021)). In such formulation, it is principally impossible to treat the heat flux coming from the interior of NSs, which drastically changes the overall temperature through electron (and radiative) thermal conductivity. Thus, we should validate the MESA burst models in particular for the physical effects inside NSs.
As the most realistic burst model which covers whole NS regions, we utilize some of them presented by Dohi et al. (2021), which follows quasi-hydrostatic evolution by using HERES. We take an approximate reaction network with 88 isotopes for mixed hydrogen and helium burning (APRX3 in Dohi et al. (2020)) and the same data of reaction rates as MESA 6 . In the energy transport equation, we implement the Schwarzschild convection. Note that convection is required for causing the mixed hydrogen/helium burning though it is somewhat artificial due to one-dimensional formulation. The initial models for our X-ray burst calculation are set to be the steady-state models (Liu et al. 2021) with gravitational compressional heating (see Matsuo et al. (2018) for details).
Let us explain model parameters in HERES. The accretion rate and compositions of accreted matter are the same as Sec 2.1. We utilize the nuclear EOS of Togashi, which has been based on the variational approach with the use of the bare nuclear potentials for two-body interaction and phenomenological three-body interaction (Togashi et al. 2017). For the heating source, standard crustal heating rates of Haensel & Zdunik (1990) are implemented. For the cooling source, we consider the slow ν cooling processes mainly composed of the modified Urca process and bremsstrahlung. The occurrence of fast ν cooling processes such as the nucleon direct Urca process, i.e., neutrino emissions induced by (inverse) β decay, could affect burst light curves (Dohi et al. 2022), but for any mass, it is prohibited with the Togashi EOS due to the quite low symmetry energy (the slope parameter L is 30 MeV) (Dohi et al. 2019).

THE IMPACT OF NEUTRON STAR MASS, RADIUS AND BASE HEATING ON TYPE I X-RAY BURST
We build a series of scenarios (models 1-12 in Table 1) with variation in mass (models 1-4), radius (models 5-8) and base heating (models 9-12), then type I X-ray bursts on the surface layer of accreting NSs are simulated by using MESA with the above inputs.
The light curves of X-ray bursts are usually characterised by several parameters, e.g., the recurrence time ∆t, which represents the time from one burst to the next. The burst duration τ , which is defined to be the time after the peak at half value of L peak . The rise time t rise is defined from transience to peak point. The e-folding time τ e is defined after peak point. The peak luminosity L peak is taken from the light curve maximum. The burst energy E b is obtained by integrating over the light curve The burst strength α is defined by the ratio of the accretion energy to the burst energy where z g is the gravitational redshift.
In order to compare with observations, we stack a sequence of bursts from each model and obtain the average light curve, burst parameters and 1σ error for them. Since the wait time for the next burst is usually shortened as the ash from the previous burst is mixed with the new fuel, i.e., compositional inertia (Taam 1980;Woosley et al. 2004), we remove the data of the first four bursts and start processing from the fifth burst. The convergence of MESA light curves is almost archived at ∼5 bursts, which is fewer than ∼10 bursts in KEPLER without nuclear preheating (see Figure A1 in Johnston et al. (2020)) 7 . 3.1. Variations in NS mass, radius and base heating, and X-ray burst parameters X-ray burst with various values of NS mass, radius and base heating are calculated. In the left panel of Figure 2, we show the luminosity of the burst sequence with different NS mass models. We calculate the averaged light curves by aligning bursts in each sequence 7 HERES light curves are converged around 10-30 burst times, which are more than those in MESA and KEPLER. This is because the HERES adopts the aniso-thermal structure as the initial model (Matsuo et al. 2018), which spends the convergence time due to the existence of thermal flux.
by their peak luminosities, the results are shown in the right panel. Similarly, the luminosity of the burst sequence with different NS radius models are shown in the left panel of Figure 3 and the averaged light curves are shown in the right panel. We find that with the increase of mass, L peak increases, ∆t increases, decay time decreases. However, as radius is increased, ∆t also increases, L peak remains constant, decay time increases.
The results for different base heating cases are shown in Figure 4, as Q e is increased, L peak decreases, ∆t decreases and decay time decreases. In the following, we calculated burst parameters such as ∆t, L peak , α, E burst , t rise , τ and τ e , one can find the values in detail in Table  1. Meanwhile, the ignition pressure P ign for each model is obtained to understand the variation of parameters. The parameters change with the variation in mass, radius and base heating are shown in Figure 5. For the models 1-4 in Table 1 (left panel in Figure 5), as M increases, ∆t, α, L peak , E b and P ign increase. For a fixed NS radius, M increases, the surface gravitational acceleration (g s ) will become larger, resulting in an increase in ignition pressure(P ign ). One can also find the ignition pressure from the bottom panel in Figure 5, which is increased with mass increases. According to one-zone model, the column density σ is expressed in two ways (Bildsten (1998); Dohi et al. (2022)): where g s = GM R 2 (1 − 2GM Rc 2 ) −1/2 , one can see the surface gravity acceleration g s on the mass-radius plane in de-   tail from Figure 3 of Dohi et al. (2021). As a result, for the fixed accretion rate and NS radius, with increasing mass, the recurrence time is proportional to P ign /g s . The increase of ignition pressure overtakes the increase of surface gravity acceleration, which leads to the increase of ∆t. The peak luminosity can be scaled as the Eddington limit (Lewin et al. 1993), which is proportional to M , but independent of R.
where κ is the electron scattering opacity. Therefore, the peak luminosity is increased with M increases. Assumed all accreted matter is processed in flashes, the burst strength is the ratio of the average luminosity emitted in the persistent X-ray emission (L p ) to that emitted in X-ray bursts (L b ) (Lewin et al. 1993): where ε G = GM/R is the gravitational energy release per gram, ε N is the nuclear energy. According to equation (5), for a fixed NS radius, α increases as mass increases. Our results from MESA simulation are almost consistent with the above simple one-zone model assumption.
In the middle panel of Figure 5, it shows that P ign and α of the bursts are inversely proportional to the radius, and ∆t and E b are proportional to the radius. While the peak luminosity L peak remains constant. As the radius increases, the gravitational acceleration on the surface of the neutron star becomes smaller, the ignition pressure decrease ∆t is longer due to the increased NS surface area. Burst energy is larger due to the longer e-folding time, to the reduced gravitational redshift from the NS surface. Burst strength α is reduced due to the lower surface gravitational potential, which also can be easily understood from equation 5. According to equation 4, as the peak luminosity does not dependent on radius, the peak luminosity almost constant with radius increases.
The results for the parameter variation with base heating are shown in the right panel of Figure 5. With the in- crease of Q e , the peak luminosity of the burst decreases continuously, and the interval between bursts becomes smaller. This is because the first hot CNO cycle, i.e., 12 C(p, γ) 13 N(p, γ) 14 O(β + ) 14 N(p, γ) 15 O(β + ) 15 N(p, α) 12 C, lasts longer with smaller Q e ; The timescale of hot CNO cycle is almost determined by abundances of 14 O and 15 O, which could trigger new (α, p) reaction paths, indirectly leading to proton-rich nucleosynthesis 8 . If Q e is smaller (in the range of 0 < Q e <0.5 MeV/u), i.e., interior NS is colder, excessive protons turn into helium, which burns to 12 C by the 3α reaction at faster rates because it takes more time to accumulate the seeds, 14 O and 15 O, leading to higher ∆t. Then, e-folding time tends to be shorter because protons being critical fuel of the rp process are more exhausted. As a result, the larger energy is produced during the hot CNO cycle due to its longer duration if Q e is smaller, leading to a higher peak luminosity. We note that Q e dependence of the hot CNO cycle timescale is, in a sense, similar to the 15 O(α, γ) 19 Ne rate dependence of that, which has been studied by Fisker et al. (2006Fisker et al. ( , 2007.

Model-observation comparisons
The light curves with variations in NS mass, radius and base heating are compared with observation in Figure 6, where the observed light curve of GS 1826-24 in 2007 is adopted. We include the burst anisotropy 8 At low temperature of T 4 × 10 8 K, β decays are dominant, but at high temperature, the second hot CNO cycle, 14 O(α, p) 17 F(p, γ) 18 Ne(β + ) 18 F(p, α) 15 O, occurs instead of 14 O(β + ). The resultant breakout reactions to αp and rp processes are therefore 15 O(α, γ) 19 Ne and 18 Ne(α, p) 21 Na. ξ b in the distance, dξ 1/2 b is calculated from F peak = L peak /4πd 2 ξ b . With use of the χ 2 method in Dohi et al. (2020), we can get the best fit dξ 1/2 b for each modelobservation comparison. From the left panel of Figure  6, we can see that the peak luminosity increases as mass increases, the peak luminosity is too high to fit the observation for M ≥ 1.7M . In the middle panel of Figure 6, the peak luminosity almost constant as radius increases, the light curve can be well fitted with radius in the range ∼11.2-13 km. In the right panel of figure 6, the peak luminosity decreases as base heating increases, the light curve can be well fitted with the variation of base heating in the range Q e = 0.1 − 0.4 MeV/u. Besides, we also compare the recurrence time with observation in the uppermost panel in figure 5. The burst models of 1-10 are consistent with observed values. The recurrence time is too short to interpret observations for burst models 11-12 with Q e = 0.3, 0.4 MeV/u. However, as the source distance is uncertain, which is crucial to determine the shape of the light curve. The input parameters such as metallicity, accretion rate also affect the burst light curve. It is better for us to use the MCMC method (e.g. Johnston et al. (2020)) to determine the system parameters. In our calculation, models 9 and 10 are consistent with the observation of GS 1826-24 in 2007 (whether the light curve or the recurrence time).

CODE COMPARISON
In order to validate the models for X-ray burst calculation such as MESA which solves the Newtonian hydrodynamics with the accreted layers, we adopt the realistic code HERES which solves the whole NS as a comparison. In table 2 from the appendix, we show our calculation models with HERES code. By using adopted mass, radius and accretion rate under X/Y = 2.9 and Z CNO = 0.02, we obtain several burst parameters, such as burst strength α, burst duration τ , recurrence time ∆t, total burst energy E burst , peak luminosity L peak , rise time t rise . The 1σ errors are also presented for each output parameters. The base heating inferred from the 1.4M NS with Togashi EoS is Q e =0.35 MeV/u ). The light curves calculated with the HERES code are shown in Figure 7. It shows that the recurrence time decreases as accretion rate increases, the peak luminosity almost constant. Meanwhile, we adopt the same mass, radius, metallicity, X/Y , base heating, accretion rate for MESA X-ray burst calculations. The input parameters and some of the output parameters are shown in table 1 from model 13 to model 16 in appendix. Figure 8 shows the comparison of the mean light curves between MESA and HERES calculation. The difference between two light curves is very small with accretion rateṀ = 2.5 × 10 −9 M yr −1 andṀ = 3.0 × 10 −9 M yr −1 . There are big differences for peak luminosity and luminosity at the tail parts of the light curve between two codes under accretion ratė M = 2.0×10 −9 M yr −1 andṀ = 4.0×10 −9 M yr −1 . The main difference is due to the higher hydrostatic force, i.e., higher compressional heating in HERES models (Matsuo et al. 2018), which leads to higher peak luminosity than MESA. Note that the contribution of compressional heating to total luminosity (∼ 10 38 erg/s) is around 10%. We show in Figure 9 that the differences of compressional heating luminosity L g between MESA and HERES codes with 4 different accretion rates. At low accretion rateṀ = 2.0 × 10 −9 M yr −1 , the peak luminosity of L g obtained from HERES code is much higher than that in MESA code. While for the rest three accretion rates, compressional heating luminosities are almost the same between two codes. The difference due to reaction networks, i.e., nuclear burning energy rates and compositions, appears in the tail parts, where the luminosity is higher in MESA models regardless ofṀ . Next, we also calculate models 17-20 with different base heating based on model 14 in the right upper panel in Figure 8, it shows that the lower the base heating, the higher the peak luminosity and the luminosity at the tail parts, which leads to a big deviation from the case with Q e = 0.35 MeV/u. Finally, we compared the predicted burst parameters (α, E b , ∆t, L peak , t rise , τ e ) of the two codes for a range of accretion rates in Figure 10. In both codes, the burst strength α, total burst energy E b and recurrence time ∆t are highly consistent. The differences of the peak luminosity and the tail parts of the light curve between two codes are obvious (e.g., the maximum relative errors of t rise and L peak are about ∼ 50%). Thus, for the first time, the consistency of the two codes are identified by our comparison. The differences for the peak luminosity at low accretion rate caused by the high compressional heating luminosity as shown in Figure 9, while the high luminosity at the tail parts regardless of accretion rate possibly be caused by the nuclear reaction energy and compositions. The comparison of the nuclear reaction network adopted in MESA (rp.net) and HERES (APRX3) are shown in table 3. It is worth noting that our input values and the values of the burst parameters obtained from MESA code which adopts a post-Newtonian modification for GR effects are unified to the local frame, in order to compare with the observations which was detected by a distant observer, the red shift of the parameters should be considered. In Appendix A, we show the detailed formulas to transfer the local frame quantities to the frame of a distant observer.

CONCLUSIONS
In this work, we present a set of simulations of X-ray burst with variation in NS mass, radius and base heating by using the open source code MESA. The light curves and burst parameters are obtained for each model. We find that the recurrence time, burst strength, peak luminosity and total burst energy are increased as mass increases. As radius increases, the recurrence time and total burst energy increase, the burst strength decreases, while the peak luminosity remains constant. The recurrence time, burst strength, peak luminosity and total burst energy decrease as base heating increases. The above phenomenon can be well explained with use of the simple one-zone model. One can see section 3.1 for the detailed explanation.
The codes, such as KEPLER and MESA, solve the Newtonian hydrodynamics only within the accreted regions to simulate X-ray burst. As a result, it is hard to probe the NS physics. HERES solves the TOV and energy-transport equations, hence it can include all possible physics. To assess the validity of the boundary condition on the crust and the GR correction for the Newtonian hydrodynamics calculations, we made a comparison between multi-zone burst models from MESA and HERES code for the first time. The results show that the average light curves are highly consistent under accretion rateṀ = 2.5 × 10 −9 M yr −1 anḋ M = 3.0 × 10 −9 M yr −1 . While under accretion ratė M = 2.0×10 −9 M yr −1 andṀ = 4.0×10 −9 M yr −1 , the peak luminosity and cooling tail are obviously different between two codes. However, the burst strength, total burst energy and recurrence time are consistent between two codes regardless of accretion rate. It is worth noting that the light curves are inconsistent when we choose other values of Q e .
We demonstrate that the NS mass, radius and base heating have a non-negligible effect on the X-ray burst simulation. The validity of the boundary condition and GR correction for MESA code are verified by the code comparison between MESA and HERES. The variation trend of the output parameters with different NS mass, radius, base heating and accretion rate can help us to understand the properties of NS via X-ray burst observations.
The difference of X-ray burst codes appears in not only light curves but also rp-process nucleosynthesis. In fact, Parikh et al. (2008) showed the difference in final products among three burst models with post-process calculation. A similar comparison with the use of MESA and HERES may also give information on some model parameters and will be present in near future.

APPENDIX
A. CORRECTING THE QUANTITIES FROM LOCAL FRAME TO THE FRAME OF A DISTANT OBSERVER To compare the burst parameters with observations, it is crucial to correct the GR quantities from local reference frame of the NS surface to the frame of a distant observer, which we mark with the superscript "∞".
The timescale will be redshifted by The redshifted luminosity can be written as Because the burst energy E b is obtained by integrating over the time (see equation (1)), from equations (A1)(A2), the redshifted burst energy is given by Similarly, the redshifted mass accretion rate is given byṀ The local burst strength from equation (2) can be redshifted by One can transfer the time scales (e.g. recurrence time ∆t, rise time t rise , duration time τ , e-folding time τ e ), the luminosities ( e.g. peak luminosity L pk ), the burst energy E b , mass accretion rateṀ and the burst strength α from the local reference frame to an observer frame by the above formulas. Note-(a) The base heat calculated from the luminoty value on the crust surface is Q b = 0.35 MeV/u ).