{'fraction_non_alphanumeric': 0.039744432960631404, 'fraction_numerical': 0.026825143286667293, 'mean_word_length': 4.417918045304148, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Direct diagnostics of forming massive stars: stellar pulsation and periodic variability of maser sources

The 6.7 GHz methanol maser emission, a tracer of forming massive stars, sometimes shows enigmatic periodic flux variations over several 10-100 days. In this Letter, we propose that this periodic variations could be explained by the pulsation of massive protostars growing under rapid mass accretion with rates of Mdot>10^-3 Msun/yr. Our stellar evolution calculations predict that the massive protostars have very large radius exceeding 100 Rsun at maximum, and we here study the pulsational stability of such the bloated protostars by way of the linear stability analysis. We show that the protostar becomes pulsationally unstable with various periods of several 10-100 days, depending on different accretion rates. With the fact that the stellar luminosity when the star is pulsationally unstable also depends on the accretion rate, we derive the period-luminosity relation log (L/Lsun) = 4.62 + 0.98log(P/100 day), which is testable with future observations. Our models further show that the radius and mass of the pulsating massive protostar should also depend on the period. It would be possible to infer such protostellar properties and the accretion rate with the observed period. Measuring the maser periods enables a direct diagnosis of the structure of accreting massive protostars, which are deeply embedded in dense gas and inaccessible with other observations.


INTRODUCTION
Massive stars have significant impacts on the interstellar medium through various feedback processes such as supernova explosions, stellar winds, and ultraviolet radiation. These feedback processes would be also important for shaping the stellar initial mass function because stars are mostly formed in clusters including massive stars (e.g., Lada & Lada 2003). However, our understanding of the massive star formation is still limited with observational difficulties, one of which arises from the fact that forming massive stars are deeply embedded in obscuring dense gas.
Observing maser emission is one of the possible methods to see the vicinity (< 10 3 AU) of forming massive stars because of the high brightness. It is known that the 6.7 GHz methanol maser emission can be thought to be often associated with circumstellar disks around forming massive stars (e.g., Norris et al. 1993; Bartkiewicz et al. 2009;Sanna et al. 2010). Some of the methanol masers show periodic flux variations over 10 days (e.g., Goedhart et al. 2004Goedhart et al. , 2009). Since the 6.7 GHz methanol masers are radiatively pumped by infrared emission of warm dusts (∼ 150 K; Cragg et al. 2005), the periodicity could reflect the luminosity variation of nearby forming massive stars or accretion disks.
Several authors have proposed different explanations for this periodic flux variations, e.g., colliding-wind binary (van der Walt 2011), and periodic accretion onto binary systems (Araya et al. 2010). However, these binarybased models do not explain why the periodic variations shorter than 10 days have not been found yet, despite the fact that a number of OB star eclipsing binaries have the periods of < 10 days (Harries et al. 2003;Hilditch et al. 2005) In this Letter, we propose an alternative picture that the periodic variability of the maser emission can be due to the pulsation of protostars growing via rapid mass accretion withṀ * 10 −3 M ⊙ yr −1 , which is expected in the massive star formation (e.g., Osorio et al. 1999;McKee & Tan 2003;Zhang et al. 2005;Beltrán et al. 2006;Krumholz et al. 2009). Our previous work shows that, by numerically modeling the stellar evolution, the massive protostar should have the large radius exceeding 100 R ⊙ with such the high accretion rate (Hosokawa & Omukai 2009;Hosokawa et al. 2010, hereafter HO09 and HYO10 respectively). We here examine the pulsational stability of the bloated massive protostars with the linear stability analysis, and show that the observed periodicity of ∼ 10 − 100 days can be well explained with the stellar pulsation. Our models predict that, with the observed period of the pulsation, we could infer the radius, mass, and luminosity of the protostars as well as the accretion rate onto the protostar. Measuring the maser periods would thus make a direct diagnosis of the central smallscale ( a few AU) structure of forming massive stars, which is beyond observations in the optical and infrared bands.

MASSIVE STARS
We study the pulsational stability of massive protostars growing at constant accretion rates ofṀ * ≥ -Evolution of the protostellar radius with various accretion ratesṀ * = 10 −4 , 10 −3 , 3 × 10 −3 , and 6 × 10 −3 M ⊙ yr −1 (taken from HO09). Symbols indicate the stellar models we conduct the linear stability analysis. The circles and stars represent the pulsationally stable and unstable models. The shaded layer denotes the instability strip, where the protostar becomes unstable.
10 −4 M ⊙ yr −1 . We here adopt the protostellar models with spherical accretion taken from HO09. In Section. 3, we will discuss that the effects of the accretion geometry (disk accretion) does not change the qualitative results from the spherical accretion case. Figure 1 shows the evolution of the stellar radius as the stellar mass increases with different accretion rates. The outline of the evolution is briefly summarized as follows (see HO09 for the details). Initially, the stellar radius gradually increases with increasing the stellar mass. After this stage, e.g., M * 10 M ⊙ forṀ * = 10 −3 M ⊙ yr −1 , the protostar contracts by losing its energy via radiation (Kelvin-Helmholtz or KH contraction). The stellar central temperature increases during this contraction stage and finally reaches 10 7 K. The hydrogen burning is ignited, and the protostar reaches the zero-age main sequence (ZAMS, except with 6 × 10 −3 M ⊙ yr −1 ). After this point, e.g., M * ≃ 40 M ⊙ forṀ * = 10 −3 M ⊙ yr −1 , the stellar radius increases again as the stellar mass increases. We here note that the maximum stellar radius during this evolution is larger with the higher accretion rate. This is because the accreting gas has the higher specific entropy with the more rapid mass accretion, which leads to the higher average entropy in the stellar interior (also see HO09). As a result, the maximum stellar radius exceeds 100 R ⊙ with the high accretion rateṡ M * 10 −3 M ⊙ yr −1 .
We apply the linear stability analysis to the above protostellar models (see e.g., Cox 1980;Unno et al. 1989;Inayoshi et al. 2013 for the details). We here consider radial (spherical) perturbations with the time dependence of e iσt , where σ = σ R + iσ I is the eigen frequency, σ R is the frequency of the pulsation, and |σ I | is the growth or damping rate of the perturbation. The protostar is pulsationally stable (unstable) with the positive (negative) σ I . If unstable, the perturbation grows until reaching the non-linear regime, where the pulsation energy is dissipated by shock waves near the stellar surface. The dissipated energy is partly converted into the kinetic energy of periodic outflows (e.g., Appenzeller 1970; Yoon & Cantiello 2010). -The period-luminosity (P-L) relation of forming massive protostars. Thin black (blue) curves show the evolutionary tracks with the spherical (disk) accretion at the rates oḟ M * = 10 −4 , 10 −3 , 3×10 −3 , and 6×10 −3 M ⊙ yr −1 . The star symbols on the tracks denote the pulsationally unstable models. The eigen frequencies are plotted for the stable stellar models. The shaded layer shows the instability strip as in Figure 1. The solid line represents the P-L relation given by equation (1), which fits the unstable models. The filled and open symbols represent the observed sources, whose distances are measured with the trigonometric parallax and kinematics. The triangles indicate that the sources are associated with ultra/hypercompact H II regions (Walsh et al. 1998, Reid et al. 2009a, Urquhart et al. 2007). Figure 1 represent the stellar models for which we conduct the linear stability analysis. Our calculations show that the protostar becomes pulsationally unstable only when the stellar radius maximally expands at a given accretion rate. This instability is caused by the κ mechanism in the He + ionization layer, where the radiative energy flux is blocked and converted into the pulsation energy (e.g., Cox 1980;Unno et al. 1989). In the KH contraction stage, the stellar surface temperature increases and the He + ionization layer disappears. The protostar is consequently stabilized and the pulsation ceases. Although the protostar would be also unstable with the lower accretion rates ( 10 −3 M ⊙ yr −1 ), the growth time is much longer than the duration for which the star is in the instability strip (σ −1 I ≫ M * /Ṁ * ∼ 10 4 yr). The perturbation does not grow enough to cause the stellar pulsation in this case. The instability strip thus does not extend for M * 10 M ⊙ , where the star becomes unstable with the lower accretion rate (see Figure  1). Figure 2 shows the evolution of the pulsation period and stellar luminosity in the examined cases. In each case, the stellar luminosity monotonically increases as the stellar mass increases. The pulsation period increases with the stellar mass in the early expansion phase, and decreases in the later KH contraction phase. The protostar becomes pulsationally unstable around the turning point of the period, which corresponds to that of the stellar radius as seen in Figure 1. Figure 2 also present the observed sources with periodic flux variations in the 6.7 GHz methanol masers, whose parameters are summarized in Table 1. The luminosities of these sources are estimated with the far-infrared data of Infrared Astronomical Satellite (IRAS) following Casoli et al. (1986) and Guzmán et al. (2012). We see that the pulsation periods of the unstable models are several 10 − 100 days, the same order of the observed periods of the maser sources. Although the pulsation period is shorter than 10 days with the low accretion rateṡ M * 10 −3 M ⊙ yr −1 , the instability strip does not cover such cases as explained above. This explains well why the maser sources which have such the short periodic variability have not been observed.

PERIOD-LUMINOSITY RELATION
Our calculations predict that both the period and luminosity of the pulsationally unstable protostars increase with the accretion rate. We thus predict a positive correlation between the maser periodicity and the stellar luminosity. We fit the unstable models by a single power law and obtain the period-luminosity (P-L) relation (thick line in Fig. 2) log(L/L ⊙ ) = 4.62 + 0.98 log(P/100 days), which is similar to that of the Cepheids (L ∝ P 4/3 ), which is well-used as a cosmic distance ladder (e.g., Tammann et al. 2003;Sandage et al. 2004). The above P-L relation and the instability strip can explain the distribution of the observed parameters in the periodic methanol masers within errors of an order of magnitude. We show that, in Appendix A, a similar P-L relation is analytically derived.
The instability strip shown in Figures 1 and 2 is not broad. The duration for which a protostar spends in the strip, i.e., the prospective lifetime of each periodic maser source is only ∼ 10 3 years. This actually matches the rarity expected with observations. The number of 6.7 GHz methanol masers observed ever is ∼ 900 (e.g., Pestalozzi et al. 2005;Caswell et al. 2011;Green et al. 2012), and 56 masers of that have been monitored for more than a year at intervals shorter than a week or month (Goedhart et al. 2007(Goedhart et al. , 2009Araya et al. 2010;Szymczak et al. 2011). About 20 % of such the sources show the characteristic periodic variabilities (Table 1). Figure 2 might suggest that our models are more applicable to the sources with the shorter periods, e.g., P 0.5 year, which is only ≃ 5 % of the monitored sources. Since the appearance time of the methanol masers during massive star formation is thought to be ≃ 3 × 10 4 yr (van der Walt 2005), the lifetime of the periodic maser sources should be as short as expected with our calculations.
An uncertainty of the P-L relation arises from possible variations of our stellar evolution tracks. While HO09 study the protostellar evolution with the spherical accretion, for instance, HYO10 present that the evolution slightly changes with different accretion geometries, e.g., accretion via a geometrically thin disk. 6 We also perform the stability analysis of case MD3-D-b0.1 in HYO10, whose model has the largest radius among the cases forṀ * = 10 −3 M ⊙ yr −1 . Figure 2 shows that the protostar also becomes pulsationally unstable in this case, and that the periods are 10 times longer than those  with the spherical accretion at the same rate. This effect might explain the observed sources which have the longer periods than predicted by the P-L relation (1).

CONCLUSION AND DISCUSSION
In this Letter, we have shown that the pulsation of massive protostars could explain the periodic variability of the 6.7 GHz maser sources. Our linear stability analysis clearly shows that a rapidly accreting (Ṁ * 10 −3 M ⊙ yr −1 ) massive protostar becomes pulsationally unstable when the star is bloated before reaching the ZAMS. Typical periods of the pulsation are several 10 − 100 days, which explain the observed periodicity well. The period depends on the adopted accretion rate, getting longer with the higher rate. On the other hand, the protostar with the lowerṀ * becomes unstable but does not produce the periodicity, which could also explain why the periodic variability shorter than 10 days has not been observed. Our stability analysis predicts the period-luminosity (P-L) relation for the pulsationally unstable protostars, log(L/L ⊙ ) = 4.62 + 0.98 log(P/100 day).
Moreover, our stellar evolution models predict that other stellar properties such as the mass, radius, and accretion rate as well as the luminosity should depend on the pulsation period. We show this in Figure 3. Single power-law functions which fit the results are, If the P-L relation (1) is confirmed by further observations, we can infer the above quantities with equations (2) -(4) only from the maser period. As seen in Figure  2, the periods of the observed maser sources are longer than 10 days. Our stellar-pulsation models predict that, in particular to explain the periods longer than 100 days, very rapid mass accretion withṀ * 3 × 10 −3 M ⊙ yr −1 is required. In this way, we can make a direct diagnosis of the small-scale structure of accreting massive protostars and their vicinities ( a few AU), which are difficult to see in the optical and infrared bands, by measuring the maser periods.
The stellar pulsation examined above should cause the periodic variation of the stellar luminosity L * . Although it is still uncertain where the maser emission is excited around protostars, some observations suggest that the maser sources are located on circumstellar disks at r ∼ 10 3 AU (e.g., Bartkiewicz et al. 2009). The temperature of the circumstellar disk irradiated by the stellar radiation obeys T d (r) ≃ 120 K (L * /3 × 10 4 L ⊙ ) 1/4 (r/10 3 AU) −1/2 , where r is the distance from the protostar. The variation of the stellar luminosity δL * changes the dust temperature following δT d /T d ≃ δL * /4L * . Since the response of the dust temperature to the irradiating flux variation is so rapid (∼ a few sec), δT d should synchronize with δL * . As evaluating δL * is beyond our linear analysis, we here consider the flux variation of Mira variables, which are also excited by the κ-mechanism. Interestingly, SiO maser sources associated with some Mira variables are likely pumped by the stellar radiation, and show periodic variations which synchronize with the stellar pulsation (Pardo et al. 2004). With the typical flux variations of ≃ 3 magnitudes of Mira-type stars (Samus et al. 2009), the amplitude of the variation is estimated as δT d /T d ≃ 0.3, which results in δT d ≃ 40 K around T d ≃ 120 K. Cragg et al. (2005) show that the variations of the dust temperature δT d 20 K could produce the observed amplitudes of the maser variability (∼ 10 − 100 Jy). Therefore, the stellar pulsation could produce the flux variations of the observed 6.7 GHz methanol masers.
In the above, we have supposed that the stellar radiation reaches the irradiated disk surface without significant time delays, i.e., the light path from the star to the disk surface is optically thin. This should be valid with the disk accretion, where the density above the disk quickly decreases as the gas falls onto the disk. If outflow is launched from the disk, however, the disk wind would enhance the density above the disk and increase the optical depth. According to a disk wind model by Zhang, Tan & McKee (2013), the optical depth when the protostar is pulsationally unstable is estimated as ≃ 5 (Ṁ * /10 −3 M ⊙ yr −1 ) 0.27 (κ/3 cm 2 g −1 ), whereκ is the mean opacity of dusts. The diffusion time with this optical depth is ≃ 30 days, which is shorter than the typical periods of the maser sources. This effect would smear out only variabilities shorter than a few 10 days. Nonetheless, the launching mechanism of the disk wind and the detailed density structure above the disk are still highly uncertain. Some periodic maser sources which have P 100 days also show rapid fluctuations of light curves over a few 10 days (G009.621, G022.357, and G338.935). The above estimate of the diffusion time do not explain such the rapid changes. Although the variations of these sources might be explained with different models rather than ours (also see below), it would be also possible that the disk wind is weaker and the diffusion time is shorter than evaluated above. Simultaneous monitoring of the maser sources in the infrared band will verify the relation between the periodic variations of the stellar luminosity and the maser variability.
Finally, we compare our pulsation model to alternative models for explaining the periodic variability of the 6.7 GHz methanol maser sources. For instance, van der Walt (2011) proposes that radiation from shocked gas in the colliding binary wind explains both the periodicity and shapes of the light curves of some maser sources. This model supposes a particular condition, the binary system surrounded by an ultracompact H II region (∼ 10 3 AU), and explains the periodic maser sources for which ultracompact H II regions have been observed (triangles in Fig 2). On the other hand, a number of the sources presented in Figure 2 do not show associating H II regions at least at 5 and 8 GHz. Our models could explain these sources, as the protostar becomes pulsationally unstable when the stellar radius is very large and effective temperature is only T eff ≃ 5000 K, with which the UV luminosity is too low to create a detectable H II region. This suggests that measuring the UV luminosities of massive protostars would be a key for discriminating the models. Further observations at the higher frequencies (> a few tens GHz) would detect hypercomact H II regions, which are too small and optically thick to be detected at the lower frequencies. Such the observations with the higher spatial resolution (e.g., with ALMA) are sure to advance our understanding of the enigmatic periodic variability of the maser sources.
We would like to thank Kazuyuki Omukai, Munetake Momose, Yoshinori Yonekura, Kenta Fujisawa, Mareki Honma, Yichen Zhang, Takashi Nakamura, Yudai Suwa, and Daisuke Nakauchi for their fruitful discussions. This work is supported in part by the Grants-in-Aid by the Ministry of Education, Culture, and Science of Japan (23·838 KI; 24·6525 KM) and by the project of Yamaguchi University entitled "The East-Asian VLBI Network and the Circulation of Matter in the Universe".

ANALYTIC DERIVATION OF THE P-L RELATION
In this Appendix, we analytically derive a P-L relation similar to equation (1) with considering the protostellar evolution. Stahler et al. (1986, hereafter SPS86) show that the evolution of an accreting protostar is well understood with the balance between following two timescales: the accretion time scale t acc = M * /Ṁ * , and the Kelvin-Helmholtz timescale t KH = GM 2 * /L * R * . As seen in Figure 1, the protostar gradually expands in the early stage with t acc < t KH , and turns to contract in the later stage with t acc > t KH . The turnaround of the radius occurs when t acc ≃ t KH (also see HO09 and HYO10). According to our linear stability analysis, the protostar becomes unstable at the epoch of the maximum radius, i.e., t acc ≃ t KH (see star symbols in Figure 1).
SPS86 also show that the stellar expansion in the early evolutionary stage is well described as R * ∝ M * 0.27Ṁ 0.41 * . With Kramers' law κ ∝ ρT −7/2 , which approximates the dependences of opacity in the stellar interior in this early stage, the stellar luminosity generally obeys L * ∝ M * 11/2 R * −1/2 (e.g., Cox & Giuli 1968). Using these relations, we can express the luminosity of the unstable protostar (i.e., t acc ≃ t KH ) as a function of the pulsation period P (∝ R 3/2 * M 1/2 * ) and obtain the period-luminosity (P-L) relation L * ∝ P 1.2 .