On period distribution of RR Lyr type variables in the globular cluster M3

TL;DR

This study models the period distribution of RR Lyrae stars in globular cluster M3 using stellar evolution and pulsation simulations, achieving partial agreement with observed period distributions and mode fractions.

Contribution

It introduces a combined stellar evolution and hydrodynamics approach to predict RR Lyrae period distributions in M3, including analytic relations for instability strip edges.

Findings

Good match with observed mode fractions (~75% fundamental mode)

Predicted mean period of fundamental mode pulsators is 0.79 days

Models reproduce the period distribution shape but overestimate mean periods

Abstract

Evolutionary calculations of population II stars with chemical composition of the globular cluster M3 were carried out under various assumptions about the initial stellar mass ($0.809M_\odot\le M_\mathrm{ZAMS} \le 0.83M_\odot$) and the mass loss rate parameter in the Reimers formula ($0.45\le\eta_\mathrm{R}\le 0.55$). In general, 30 evolutionary tracks of the horizontal branch stars were computed. Selected models of evolutionary sequences were used as initial conditions for solution of the equations of hydrodynamics that describe radial stellar oscillations. Hydrodynamic models of RR Lyr type stars were computed for the core helium burning stage as well as for the preceding pre--ZAHB stage. Analytic relations for the effective temperature of the instability strip edges as a function of stellar luminosity are obtained. Theoretical histograms of the period distribution of RR Lyr type variables were produced for each evolutionary sequence using Monte--Carlo simulations based on the consistent stellar evolution and nonlinear stellar pulsation calculations. A satisfactory agreement with observations (i.e. the greater number of RRab variables) was found for the evolutionary sequence $M_\mathrm{ZAMS} = 0.811M_\odot$, $\eta_\mathrm{R}=0.55$ with the number fraction of fundamental mode pulsators $\approx 75\%$. At the same time the mean period of fundamental mode pulsators ($\langle\Pi\rangle_0=0.79$ day) is substantially greater compared to the observational estimate of $\langle\Pi\rangle_\mathrm{ab}$.

Tables

Table 1. Table 1: Characteristic lifetimes of RR Lyr stars

$M_{\mathrm{ZAHB}}/M_{\odot }$	$t_{\mathrm{RR}}/t_{\mathrm{HB}}$	$t_{\mathrm{RR},\mathrm{pre}-\mathrm{ZAHB}}/t_{\mathrm{RR}}$
0.58	0.003	0.42
0.59	0.004	0.10
0.60	0.01	0.12
0.61	0.07	0.04

Equations

$$\mu ={\alpha }_{\mu }\rho \Lambda E_{\mathrm{trb}}^{1/2},$$

$$\log T_{\mathrm{eff},\mathrm{b}}=4.0561-0.1040\log (L/L_{\odot }),$$

$$\log T_{\mathrm{eff},\mathrm{r}}=3.9847-0.1096\log (L/L_{\odot }).$$

Full text

On period distribution of RR Lyr type variables in the globular cluster M3

© 2019 г.  Yu. A. Fadeyev*E–mail: fadeyev@inasan.ru**

Institute of Astronomy, Russian Academy of Sciences, Pyatnitskaya ul. 48, Moscow, 119017 Russia

Received March 14, 2019; after revision April 1, 2019;

accepted April 1, 2019

Abstract — Evolutionary calculations of population II stars with chemical composition of the globular cluster M3 were carried out under various assumptions about the initial stellar mass ($0.809M_{\odot}\leq M_{\mathrm{ZAMS}}\leq 0.83M_{\odot}$) and the mass loss rate parameter in the Reimers formula ($0.45\leq\eta_{\mathrm{R}}\leq 0.55$). In general, 30 evolutionary tracks of the horizontal branch stars were computed. Selected models of evolutionary sequences were used as initial conditions for solution of the equations of hydrodynamics that describe radial stellar oscillations. Hydrodynamic models of RR Lyr type stars were computed for the core helium burning stage as well as for the preceding pre–ZAHB stage. Analytic relations for the effective temperature of the instability strip edges as a function of stellar luminosity are obtained. Theoretical histograms of the period distribution of RR Lyr type variables were produced for each evolutionary sequence using Monte–Carlo simulations based on the consistent stellar evolution and nonlinear stellar pulsation calculations. A satisfactory agreement with observations (i.e. the greater number of RRab variables) was found for the evolutionary sequence $M_{\mathrm{ZAMS}}=0.811M_{\odot}$, $\eta_{\mathrm{R}}=0.55$ with the number fraction of fundamental mode pulsators $\approx 75\%$. At the same time the mean period of fundamental mode pulsators ($\langle\Pi\rangle_{0}=0.79$ day) is substantially greater compared to the observational estimate of $\langle\Pi\rangle_{\mathrm{ab}}$.

Keywords: stars: variable and peculiar

introduction

RR Lyr type pulsating variables observed in globular clusters are low–mass ($M\approx 0.6M_{\odot}$) population II stars on the evolutionary stage of steady–state thermonuclear core helium burning (Iben, 1974; Caputo, 1998). In the Hertzsprung–Russel diagram (HRD) the RR Lyr type variables locate on the horizontal branch within the pulsational instability strip and have effective temperatures $6000\leavevmode\nobreak\ \mathrm{K}\lesssim T_{\mathrm{eff}}\lesssim 7600\leavevmode\nobreak\ \mathrm{K}$ (Stellingwerf, 1984; Bono, Stellingwerf, 1994). Most of RR Lyr stars are the fundamental mode or the first overtone pulsators (i.e. variables of RRab and RRc type, respectively).

It has been shown by Oosterhoff (1939) that galactic globular clusters may be divided into two groups depending on the mean period and the number fractions of RRab and RRc variables. The criteria of the globular clusters of the first group are $\langle\Pi\rangle_{\mathrm{RRab}}<0.6$ day and $\langle\Pi\rangle_{\mathrm{RRc}}<0.35$ day with the relative number of RRc variables $f_{\mathrm{RRc}}<30\%$. On the contrary, in clusters of the second group $\langle\Pi\rangle_{\mathrm{RRab}}>0.6$ day, $\langle\Pi\rangle_{\mathrm{RRc}}>0.35$ day and $f_{\mathrm{RRc}}>30\%$ (van den Bergh, 1957). The Oosterhoff dichotomy seems to be due to different physical conditions corresponding to pulsational mode switching (van Albada, Baker, 1973) but the cause of the dichotomy remains ubclear.

The globular cluster M3 (NGC 5272) is one of the most thoroughly investigated clusters due to its relatively small distance and insignificant interstellar reddening. Recent distance estimates range from 4.44 kpc (Watkins, van der Marel, 2017) to 10.05 kpc (Marconi et al., 2003) and the interstellar reddening is $E=(B-V)=0.013$ (VandenBerg et al., 2016). In comparison with other globular clusters M3 is conspicuous due to the large number of RR Lyr pulsating variables. By now about 274 RR Lyr variables were discovered in this cluster (Bakos et al., 2000). The observational period distribution of RR Lyr type variables (Corwin, Carney, 2001) undoubtedly indicate that M3 belongs to the first group of Oosterhoff classification.

Theoretical analysis of horizontal branch stars is based on methods of population synthesis using results of stellar evolution calculations and observations of RR Lyr type variables. In the framework of this approach the relations between the pulsation period $\Pi$ and fundamental parameters of the star (the mass $M$, the muninosity $L$, the effective temperature $T_{\mathrm{eff}}$, the metal abundance $Z$) are obtained from the formulae approximating the extensive grids of hydrodynamic models of RR Lyr stars (Caputo et al., 1998; Marconi et al., 2003; 2015). However attempts to reproduce the period distribution of RR Lyr stars in the globular cluster M3 have failed (Catelan, 2004).

The goal of the present study is to consider the theoretical period distribution of RR Lyr variables in the globular cluster M3 using the consistent calculations of stellar evolution and nonlinear stellar pulsations. In this work we employ various assumptions on the initial stellar mass and the mass loss rate. The study implies calculation of the evolutionary sequences of horizontal branch stars and the use of selected evolutionary models as initial conditions for solution of the equations of radiation hydrodynamics and time–dependent convection describing radial stellar oscillations. Consistent evolutionary and nonlinear pulsation calculations allow us to determine the stellar age when the star crosses the instability strip edge as well as the age of the pulsational mode swithing. Finally, for each evolutionary sequence we produce the theoretical period distribution which is compared with observations.

evolutionary tracks of horizontal branch stars

Following the methodology of our preceding work (Fadeyev, 2018) the evolutionary sequences of horizontal branch stars were computed from the zero–age main sequence (ZAMS) to core helium exhaustion. The nuclear network includes 29 nuclides from hydrogen ${}^{1}\mathrm{H}$ to aluminium ${}^{27}\mathrm{Al}$ which are coupled by 51 reactions. The reaction rates were calculated using the JINA Reaclib data (Cyburt et al., 2010). The initial fractional mass abundance of helium was assumed to be $Y_{0}=0.25$ (Salaris et al., 2004) with initial metallicity $Z_{0}=0.001$ (Catelan, 2004).

Convective mixing was treated in the framework of the standard theory (Böhm–Vitense, 1958) with mixing length to pressure scale height ratio $\alpha_{\mathrm{MLT}}=\Lambda/H_{\mathrm{P}}=2.0$. In general, we computed 30 evolutionary sequences with initial masses $0.809M_{\odot}\leq M_{\mathrm{ZAMS}}\leq 0.83M_{\odot}$ for the mass loss rate parameter $0.45\leq\eta_{\mathrm{R}}\leq 0.55$ in the Reimers formula (1975). All evolutionary calculations were carried out with the program MESA version 10398 (Paxton et al., 2018).

The main difficulty encountered in evolutionary calculations of the horizontal branch stars on the stage of core helium burning is due to the jump of the helium abundance on the outer boundary of the convective core. Finite–difference representation of the evolutionary model is responsible for discrete ingestion of the helium–rich material into the convective core. As a result, the energy generation rate in the triple–$\alpha$ reactions may substantially increase so that the evolution time on the horizontal branch significantly increases due to spurious loops in the HRD (Constantino et al., 2015; 2016). In the present study to avoid appearence of spurious breathing pulses we employed the method that allows us to constrain the ingestion rate of unburned helium–rich material on the outer boundary of the convective core (Spruit, 2015; Constantino et al., 2017).

Selected models of evolutionary sequences located in the vicinity of the instability strip in the HRD were used as initial conditions for solution of the equations of radiation hydrodynamics describing radial stellar oscillations. Basic equations and the choice of parameters of the time–dependent convection theory (Kuhfuß, 1986) are discussed in our preceding papers (Fadeyev, 2013; 2015). In the present study the solution of the transport equations for turbulent convection was carried out with $\alpha_{\mu}=0.3$ in the expression for the turbulent viscosity

[TABLE]

where $\rho$ is the gas density, $\Lambda$ is the mean free path of the turbulent element (i.e. the mixing length), $E_{\mathrm{trb}}$ is the mean kinetic energy of turbulence. The parameter $\alpha_{\mu}$ determines efficiency of ineraction between turbulent elements and the gas flow. For pulsating stars this parameter ranges within $0.1<\alpha_{\mu}<0.5$ and the pulsation period is almost independent of $\alpha_{\mu}$ (Wuchterl, Feuchtinger, 1998; Olivier, Wood, 2005; Smolec, Moskalik, 2008).

For each crossing of the instability strip we computed from 10 to 15 hydrodynamic models. The edge of the instability strip was determined using two adjacent models one of which is unstable against radial pulsations and another shows decaying oscillations. The growth (or decay) rate of oscillations is $\eta=\Pi^{-1}d\ln E_{\mathrm{K,max}}/dt$, where $\Pi$ is the period of radial stellar pulsations, $E_{\mathrm{K,max}}$ is the maximum kinetic energy of pulsation motions. The kinetic energy reaches its maximum value twice per pulsation period. The evolutionary time $t_{\mathrm{ev}}$ (i.e. the star age) correspoinding to the instability edge ($\eta=0$) was evaluated by linear interpolation of $\eta(t_{\mathrm{ev}})$.

Effective temperatures of the blue $T_{\mathrm{eff,b}}$ and the red $T_{\mathrm{eff,r}}$ edges of the instability strip can be expressed as a function of the bolometric stellar luminosity $L$ by following relations

[TABLE]

The constant coefficients in (2) and (3) were obtained by the least–squares method within $1.638<\log L/L_{\odot}<1.883$ for 58 and 84 crossings of the blue and red edges of the instability strip, respectively, by 30 evolutionary tracks. The larger number of red edge crossings is due to decaying oscillations of energy generation in the triple–$\alpha$ reactions after the helium flash on the tip of the red giant branch. Formulae (2) and (3) give the effective temperature on the instability edge with mean r.m.s. deviations of $\sigma(\log T_{\mathrm{eff,b}})\approx 10^{-3}$ and $\sigma(\log T_{\mathrm{eff,r}})\approx 1.2\times 10^{-3}$. The scatter of $\eta=0$ points around the regression line is mainly due to the linear intrerpolation errors. The average width of the instability strip is $\Delta\log T_{\mathrm{eff}}=0.081$.

Location of the evolutionary track in the HRD relative to the instability strip depends on the initial stellar mass $M_{\mathrm{ZAMS}}$ and the mass lost during the preceding stage of the red giant (i.e. the parameter $\eta_{\mathrm{R}}$). Dependence of the evolutionary track on the initial mass is illustrated in Fig. 1 for evolutionary sequences $M_{\mathrm{ZAMS}}=0.81M_{\odot}$ and $M_{\mathrm{ZAMS}}=0.83M_{\odot}$ computed for the mass loss rate parameter $\eta_{\mathrm{R}}=0.5$. The role of mass loss is illustrated in Fig. 2 for evolutionary sequences $M_{\mathrm{ZAMS}}=0.82M_{\odot}$ computed for $\eta_{\mathrm{R}}=0.45$, 0.5 and 0.55.

Solid lines in Figs. 1 and 2 correspond to the evolutionary stage of core helium burning. For all evolutionary sequences considered in the present study the time of core helium burning with an accuracy of a few percent is $t_{\mathrm{HB}}\approx 10^{8}$ yr. The preceding evolutionary stage when the star leaves the tip of the red giant branch and approaches the horizontal branch is significantly shorter. For example, the time interval between the maximum energy generation during the helium flash and commencement of core helium burning is $\approx 1.4\times 10^{6}$ yr. In Figs. 1 and 2 this stage of evolution is shown by dotted lines and in close vicinity of the instability strip its duration is $t_{\textrm{pre-ZAHB}}\lesssim 10^{6}$ yr. The point between the solid and dotted lines corresponds to the horizontal branch of the zero age (ZAHB).

periods of radial pulsations

Computations of each hydrodynamic model were carried out on the time interval encompassing hundreds of pulsation cycles. After completion of hydrodynamic computations the pulsation periods of the fundamental mode $\Pi_{0}$ and the first overtone $\Pi_{1}$ were evaluated using the discrete Fourier transform of the pulsating stellar envelope kinetic energy. Both periods are time–independent because after attainment of the limiting amplitude nonlinear effects remain negligible. Therefore, to enhance the accuracy of period determination we calculated the power spectrum of the kinetic energy for the whole time interval of the solution of the equations of hydrodynamics. While approaching the limiting cycle the amplitude of one mode in the spectrum reduces whereas another mode becomes prevalent.

During the instability strip crossing the star undergoes $\sim 10^{9}$ oscillations on the stage of core helium burning and $\sim 10^{7}$ oscillations on the pre–ZAHB stage. Therefore, we may assume with high accuracy that the mode switching occurs instantly. In the present study the stellar age $t_{\mathrm{ev}}$ corresponding to the mode switching was determined as a mean evolutionary time of two adjacent hydrodynamic models pulsating in different modes.

The pulsation period as a continuous function of evolutionary time $\Pi(t_{\mathrm{ev}})$ was determined using the cubic interpolation splines. Results of approximation are shown in Figs. 3 and 4 for evolutionary sequences displayed in Figs. 1 and 2. For the sake of graphical clarity the evolutionary time $t_{\mathrm{ev}}$ is set to zero at ZAHB, the time scale for $t_{\mathrm{ev}}<0$ being nearly two orders of magnitude shorter in comparison with that for $t_{\mathrm{ev}}>0$. Each plot in Figs. 3 and 4 represents the temporal dependence of the pulsation period within the instability strip and a discontinuous jump corresponds to the mode switching.

As clearly seen in Figs. 3 and 4 the lifetime of the star within the instability strip $t_{\mathrm{RR}}$ significantly reduces with decreasing mass of the horizontal branch star (i.e. with decreasing $M_{\mathrm{ZAMS}}$ or increasing $\eta_{\mathrm{R}}$). Approximate estimates of $t_{\mathrm{RR}}$ in units of the core helium burning time $t_{\mathrm{HB}}$ are given in the second column of the Table 1 for several values of the zero age horizontal branch stellar mass $M_{\mathrm{ZAHB}}$. The instability strip lifitime on the pre–ZAHB stage $t_{\mathrm{RR,pre-ZAHB}}$ is almost independent of stellar mass, so that decrease of $M_{\mathrm{ZAHB}}$ is accompanied by increasing number fraction of pre–ZAHB RR Lyr variables (see the third column in Table 1).

period distribution

Theoretical histograms of RR Lyr period distribution were produced for all computed evolutionary sequences using $10^{7}$ Monte–Carlo simulations. To this end the age of the synthetic star $t_{\mathrm{ev}}$ was randomly selected from a uniform distribution on the interval $t_{\textrm{pre-ZAHB}}\leq t_{\mathrm{ev}}\leq t_{\mathrm{HB}}$, where $t_{\textrm{pre-ZAHB}}=-10^{6}$ yr.

First of all, we have to note that histograms of evolutionary sequences computed with $\eta_{\mathrm{R}}\leq 0.5$ should be excluded from our consideration due to contradiction with observations because they exhibit significant excess of the first overtone pulsators. Moreover, as clearly seen in Fig. 2, the maximum effective temperature of the horizontal branch stars seems to be insufficiently high in comparison with its observational estimates (Catelan et al., 2001; Cacciari et al., 2005).

Fig. 5 shows the normalized histograms of the period distribution produced from evolutionary and stellar pulsation computations of three evolutionary sequences with initial masses $M_{\mathrm{ZAMS}}=0.81M_{\odot}$, $0.818M_{\odot}$, $0.828$ and the mass loss rate parameter $\eta_{\mathrm{R}}=0.55$. The masses of RR Lyr stars are $M=0.578M_{\odot}$, $0.589M_{\odot}$ and $0.603M_{\odot}$, respectively. As clearly seen, to avoid the contradiction with observations and to obtain excess of fundamental mode pulsators we have to assume that the initial mass is $M_{\mathrm{ZAMS}}=0.81M_{\odot}$.

For a more detailed comparison between results of computations and observations Fig. 6 shows normalized histograms of evolutionary sequences with initial masses $0.809M_{\odot}\leq M_{\mathrm{ZAMS}}\leq 0.812M_{\odot}$ where each plot is accompanied by the number fraction of fundamental mode pulsators $f_{0}$. As clearly seen, the satisfactory agreement with observations is obtained for evolutionary sequences $M_{\mathrm{ZAMS}}=0.810M_{\odot}$ and $M_{\mathrm{ZAMS}}=0.811M_{\odot}$, that is for RR Lyr stars with masses $0.578M_{\odot}\leq M\leq 0.580M_{\odot}$. The number fractions of fundamental mode pulsators (71% and 75%) are enough close to 78% of RR Lyr variables in the globular cluster M3 (Catelan, 2004). At the same time one should note that the mean period of fundamental mode pulsators $\langle\Pi\rangle_{0}=0.79$ day is significantly larger in comparison with the mean period of RRab variables $\langle\Pi\rangle_{\mathrm{ab}}=0.56$ day.

The existence of the narrow interval of intial masses $M_{\mathrm{ZAMS}}$ with satisfactory agreement between theory and observations is due to the fact that the mode switching is sensitively dependent on the stellar mass. This dependence is illustrated in Fig. 7 where the periods of the fundamental mode and the first overtone at the mode switching are plotted as a function of $M_{\mathrm{ZAMS}}$. The plots correspond to the final core helium burning stage when radial pulsation mode changes from the first overtone to the fundamental mode. As is seen, in the initial mass interval $0.809\leq M_{\mathrm{ZAMS}}\leq 0.812M_{\odot}$ the period of each mode at the mode switching varies by $\approx 15\%$. A strong sensitivity of the pulsation period at the mode switching to the stellar mass is the main cause of different shapes of histograms displayed in Fig. 6.

conclusions

In the present study we carried out consistent calculations of stellar evolution and nonlinear stellar pulsations for 30 evolutionary sequences of horizontal branch stars. Initial masses of evolutionary models correspond to the age of RR Lyr stars from $1.195\times 10^{10}$ yr for $M_{\mathrm{ZAMS}}=0.83M_{\odot}$ to $1.309\times 10^{10}$ yr for $M_{\mathrm{ZAMS}}=0.809M_{\odot}$. Almost all theoretical histograms exhibit excess of short–period RRc type variables and contradict to observations of the globular cluster M3. A satisfactory agreement with observations was obtained for two evolutionary sequences $M_{\mathrm{ZAMS}}=0.810M_{\odot}$ and $0.811M_{\odot}$ computed with the mass loss rate parameter $\eta_{\mathrm{R}}=0.55$. The number fraction of fundamental mode pulsators (71% and 75%) was found to be insignificantly smaller than the observed number fraction of RRab variables (78%).

The evolutionary sequences of horizontal branch stars computed in the present study will be used to improve the theoretical estimate of the mean period of RR Lyr variables pulsating in the fundamental mode. Bearing in mind the fact that the mode switching is strongly sensitive to the stellar mass we should consider in the future work the evolutionary sequences computed in vicinity of the values $M_{\mathrm{ZAMS}}=0.811M_{\odot}$, $\eta_{\mathrm{R}}=0.55$ and $Y_{0}=0.25$.

Solution of the problem is complicated by requirements of more thorough evolutionary calculations. The instability strip crossing on the stage of core helium burning occurs when the fractional mass abundance of helium in the convective core is $Y<0.1$, so that effects of unsteady mass ingestion on the outer boundary of the convective core become stronger and lead to perceptible changes in evolutionary calculations. As was shown above reduction of such effects is necessary for correct determination of the mode switching and thereby for better agreement between theory and observations.

The author is indebted to L.R. Yungelson for critical comments and useful discussions.