On the Granulation and Irregular Variation of Red Supergiants

The mechanisms and characteristics of the irregular variations of red supergiants (RSGs) are studied based on the RSG samples in the Small Magellanic Cloud (SMC), Large Magellanic Cloud (LMC), and M31. With the time-series data from the All-Sky Automated Survey for SuperNovae and Intermediate Palomar Transient Factory survey, we use the continuous time autoregressive moving average model to estimate the variability features of the light curves and their power spectral density (PSD). The characteristic evolution timescale and amplitude of granulations are further derived from fitting the posterior PSD with the COR function, which is a Harvey-like granulation model. The consistency of theoretical predictions and results is checked to verify the correctness of the assumption that granulations on RSGs contribute to irregular variation. The relations between granulation and stellar parameters are obtained and compared with the results of red giant branch (RGB) stars and Betelgeuse. It is found that the relations are in agreement with predictions from basic physical process of granulation and fall close to the extrapolated relations of RGB stars. The granulations in most of the RSGs evolve at a timescale of several days to a year with a characteristic amplitude of 10–1000 mmag. The results imply that the irregular variations of RSGs can be attributed to the evolution of granulations. When comparing the results from the SMC, LMC, and M31, the timescale and amplitude of granulation seem to increase with metallicity. The analytical relations of the granulation parameters with stellar parameters are derived for the RSG sample of each galaxy.

In general, RSGs show some degree of variability in the visual and infrared bands (Levesque et al. 2007;Yang et al. 2018). According to the characteristics of the light curves, RSGs are divided into three types: semiregular variables, variables with a long secondary period (LSP), and irregular variables (Yang & Jiang 2011, 2012Ren et al. 2019). The irregular variation, present in all three types of RSGs, becomes very popular in comparison with the bluer variables such as Cepheids or the less luminous red variables such as Miras. Meanwhile, the study of irregular variation is hindered by the difficulty of characterizing the features because no simple period can be retrieved from the irregular variation. The widely accepted mechanism for irregular variation of RSGs is large convection cells first proposed by Schwarzschild (1975), since RSGs have a convective envelope due to the high opacity in the low-temperature extended atmosphere. The study suggested that the surface of RSGs should be dominated by a few large convective cells (Schwarzschild 1975), which was later confirmed by the 3D simulation (Freytag et al. 2002), and that the convective outer layers of RSGs cause both the spectrum and the photometry to be variable (Freytag et al. 2002;Josselin & Plez 2007). Observationally, the 1/f (colored) noise component found in the power spectra means that the irregular variation may be driven by stochastic mechanisms (Kiss et al. 2006;Yang & Jiang 2008) related to the large convection cells. Kiss et al. (2006) found a strong 1/f noise component in almost all of the stars in their sample, which suggests that convection plays an important role in all types of RSGs. With increasing power of convection, the light variation of an RSG becomes more and more irregular.
Betelgeuse (α Ori), as the nearest (∼130-200 pc; see van Leeuwen 2007 andHarper et al. 2008) and brightest RSG, has been studied in detail in many aspects, among which is the granulations as the observable appearance of convective cells. Chiavassa et al. (2009Chiavassa et al. ( , 2010Chiavassa et al. ( , 2011 compared their 3D simulations of RSG convection with the interferometric observations of Betelgeuse. They concluded that the granulation pattern is detected on the surface of Betelgeuse in both the optical and the H band based on the excellent fits to the observed visibility points and closure phase. They further determined the size of the convection cell to be from about 5to 30 mas (i.e., about 200-1200 R e at a distance of 200 pc), which implies a small number of granulations if considering the diameter of Betelgeuse to be about 1800 R e . Other works (Monnier 2003;Haubois et al. 2009;Ohnaka et al. 2009Ohnaka et al. , 2011Kiss et al. 2010) also found evidence for the existence of large convective cells in Betelgeuse.
There are only a few RSGs like Betelgeuse that can be resolved by interferometric observation thanks to a very nearby location. In this sense, Betelgeuse is an exception; however, the granulations of normal RSGs cannot be studied in the Betelgeuse way. Fortunately, light variation provides an alternative observable to investigate the granulations and their relation to stellar parameters. Mathur et al. (2011) systematically analyzed the Kepler light curves of ∼1000 red giants and successfully determined the granulation timescale τ gran and power P gran . They also derived the relation of the granulation characteristics to stellar parameters. Similarly, Derekas et al. (2017) analyzed the Kepler light curve, but of the Cepheid variable V1154 Cyg, and retrieved the effective timescale of granulation to be about 3.5 days, in agreement with the extrapolation of red giants. For comparison, they estimated the timescale of the granulation τ gran of Betelgeuse to be 280±160 days from the power density spectrum of the visual light curve assembled by the AAVSO organization, consistent with the timescale of convective motions in Betelgeuse derived from line bisector variation (Gray 2008). These studies justified the possibility of estimating the features of RSG granulations from their light curves and their relation to stellar parameters, which is the intention of this work.
The difference for RSGs with previous studies of red giants and Cepheids lies in the lack of the Kepler observation, i.e., long-term continuous and accurate photometry. However, various time-domain surveys have collected large data sets, in particular toward the nearby galaxies where RSGs are detectable for their high luminosity. In addition, extragalactic RSGs share the same distance, which brings about the convenience to estimate stellar parameters with high relative precision. On the contrary, RSGs in our Galaxy suffer serious interstellar extinction, which messes with the measurement of stellar distance and parameters. Following our previous study of the RSG light variation, the nearby galaxies, namely the Large Magellanic Cloud (LMC), Small Magellanic Cloud (SMC), M31, and M33, are chosen as the targets. On the light variation data sets, we choose the All-Sky Automated Survey for SuperNovae (ASAS-SN; Shappee et al. 2014;Kochanek et al. 2017), which covers the SMC and LMC, and the Intermediate Palomar Transient Factory (iPTF) survey Rau et al. 2009), which observed M31 and M33. We first describe the selection of the RSG sample and data for light variation in Section 2, then the method to derive the granulation parameters in Section 3 followed by the analysis of the relation between granulation and stellar parameters in Section 4.

The Preliminary Sample
The most distinguishable features of RSGs are their high luminosity and low effective temperature. Based mainly on nearinfrared brightness and colors, RSGs are selected in the nearby galaxies. Historically, there have been many collections of RSG samples in the Magellanic Clouds (e.g., Feast et al. 1980;Catchpole & Feast 1981;Wood et al. 1983;Pierce et al. 2000;Massey 2002;Massey & Olsen 2003;Levesque et al. 2006;Massey et al. 2007;Kastner et al. 2008;van Loon et al. 2008van Loon et al. , 2010. Recently, Yang et al. (2018) and Yang et al. (2019Yang et al. ( , 2020 selected 773 and 1405 RSG candidates in the LMC and SMC, respectively, by combining most of the available infrared observations and carefully checking the brightness and colors; this is the largest RSG sample in the Magellanic Clouds to date. For M31 and M33, we adopted the RSG samples in our previous work (Ren et al. 2019), which contain 420 and 717 RSGs in M31 and M33, respectively. Ren et al. (2019) found that the minimum luminosity of RSGs in M31 and M33 is about 1 mag above the theoretical limit of RSGs corresponding to 9 M e . This is understandable because the observations are limited in particular for the relatively distant M31 and M33. Besides, the selection criteria are not uniform. Judging from the number of RSGs in the samples, i.e., 773, 1405, 420, and 717 RSGs in the LMC, SMC, M31, and M33, respectively, we can tell that the samples are incomplete. Nevertheless, they are the largest and seemingly most reliable RSG samples in these galaxies. In addition, the samples, though incomplete, provide many cases to study the granulation and irregular variation in RSGs. Therefore, we take them as the initial samples.

The Time-series Data
As mentioned earlier, the time-series data are taken from the ASAS-SN for RSGs in the LMC and SMC and the iPTF survey for RSGs in M31 and M33. The ASAS-SN consists of 24 telescopes all around the globe that observed variable objects in the SMC and LMC at a cadence of 4-5 days for approximately 1600 days (Shappee et al. 2014;Kochanek et al. 2017). The ASAS-SN has magnitude limits of g∼18 and V∼17.3 and photometric precision of 0.08 mag at V∼16 (Jayasinghe et al. 2018). The cadence, duration, and sensitivity all fit the purpose of studying the irregular variation of RSGs in the Magellanic Clouds. The data are retrieved with a 3″ radius positional match. For M31 and M33, the iPTF survey monitored for about 2000 days in the g and R bands, reaching 20.5 mag at a 5σ level Rau et al. 2009). With a distance modulus of about 24.5 mag, this sensitivity should be compatible with the brightness of RSGs if there is not much dimming by interstellar extinction. Since Ren et al. (2019) already performed photometry of the RSGs in M31 and M33 from the iPTF images, their results are adopted. However, it is noted in the later analysis that the data for RSGs in M33 lead to a very uncertain determination of the granulation timescale, which may be caused by insufficient sampling, since the light curves of RSGs in M33 have a similar time span as M31 but only onethird of the observational points. The RSG sample of M33 is thus discarded in further study and will have to wait for better data collection. In practice, the g-band photometry is much less sampled and usually of lower quality in comparison with V or R band due to the red color of RSGs and the facility sensitivity skewed to longer wavelengths; therefore, the time-series data to be used are from the V-band observation for the SMC and LMC and the R-band observation for M31.

The Photometric Data
The stellar parameters, mainly the luminosity and effective temperature, are derived from multiband photometry. We take the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) measurement in the J, H, and K S bands. According to the study of Ren et al. (2019), 18 RSGs in M31 with less than 250 photometric data points are classified as undefined type (U-type) and discarded in the following analysis because the observational data are not enough for the study of variation. Finally, we collect 128, 385, and 359 light curves of RSGs in the SMC, LMC, and M31, respectively. Since the granulation signal could be detected in all categories of RSGs, all of the RSGs with sufficient time-series data are analyzed further, no matter whether they have been classified into semiregular variables, irregular variables, or variables with LSP (Yang & Jiang 2011, 2012Ren et al. 2019).

The Granulation Model
The major parameters that characterize the granulations are the timescale τ gran and the amplitude σ gran . The common way to obtain the granulation parameters uses the background of the power spectrum. Harvey (1985) constructed a simple but reasonable model to describe the granulations that is widely used, e.g., in red giants. The idea is that the autocovariance of granulation evolution over time can be approximated by an exponential decay function with a characteristic timescale τ gran and variance σ 2 . This assumption leads to a Lorentz profile power spectrum, in which P(ν) is the total power at frequency ν. It should be noted that the exponent in the denominator ("2" here) is changeable (Harvey 1985). Some modified Harvey functions are used to fit the background of the power spectral density  Mosser et al. 2011) to fit the PSD of red giants and retrieve the granulation parameters from the Kepler observations. Following them, we also use the COR function to fit the background of the PSD and obtain the granulation parameters in RSGs. The COR function takes the following form: where ξ is a normalization factor that depends on the value of α. Moreover, from Karoff et al. (2013), ξ is related to α as

Continuous Time Autoregressive Moving Average Models
The ASAS-SN and iPTF surveys are ground-based, which makes the sampling irregular. Moreover, there are usually gaps in the time-series data. Both the irregular sampling and gaps in data would deform the PSD when estimated directly from the light curves of irregular variation (Kelly et al. 2009(Kelly et al. , 2014. Besides, the fitting would be affected by the pulsation and noise of measurement when we use the COR model to fit the background of the PSD. In order to solve this problem, Kelly et al. (2014) developed a flexible and scalable method, the Continuous-time AutoRegressive Moving Average (CARMA) model, to estimate the variability features of the light curves and their PSDs. The CARMA(p q , ) process has been used to model X-ray variability of active galactic nuclei (Kelly et al. 2009(Kelly et al. , 2011) and X-ray binaries (Nowak 2000); it is the first time to fit the granulation background using CARMA models. We make use of this model to compensate for the irregularly sampled data. A CARMA(p q , ) model can be defined as the solution to the stochastic differential equation in which y(t) is the time-series data, a a -,..., p 0 1 are autoregressive coefficients, β 1 ,..., β q are moving average coefficients, and ò(t) is Gaussian white noise with zero mean and variance 1. For a stationary stochastic process, p<q, and it is necessary that the roots r r , ..., p 1 of the autoregressive polynomial According to the Wiener-Khinchin theorem, the PSD of a wide-sense stationary stochastic process is the Fourier transform of its autocovariance function. Then there is a PSD corresponding to a CARMA(p q , ) process: The details of the mathematics and algorithm can be found in Kelly et al. (2014). In this algorithm, the values of the parameters (p q , ) must be chosen for time-series analysis, which is usually done by using the Akaike information criterion (AIC; Akaike 1973). The AIC is an estimator of the relative information lost of models for a given set of data, and AIC provides a method for model selection. Hurvich & Tsai (1989) made a correction to the AIC for small samples (hereafter AICc). In our work, the CARMA(p, q) model that minimizes the AICc yields p=3, q=0. Afterward, the CARMA(3, 0) parameters and the posterior PSD are calculated as the product of the likelihood function with a prior probability distribution using a Markov Chain Monte Carlo (MCMC) method.
To specify the prior distribution, we need to build the relations between the coefficients of the CARMA(3, 0) model and the physical parameters of convection (i.e., timescale τ and characteristic amplitude σ). This is done by considering a CARMA(1, 0) model that is also known as the Ornstein- As can be seen, P(ν) is a Lorentzian function, and ( ) R t 1,0 decays with a timescale τ=1/α 0 , which is related to the physical parameters τ and σ. Since the PSD of a CARMA(3, 0) process can be expressed as a sum of three Lorentzian functions, the convection characteristics can be deduced from the model. We set a prior for each O-U process and assume a uniform prior on α 0 =1/τ. The lower and upper bounds chosen for the uniform prior distribution of 1/τ are set to 1/0.01 and 1/2000, which implies that the timescale of convection can range from 0.01 to 2000 days. For a characteristic amplitude σ, we set a uniform distribution on autocovariance at lag t=0 (see Equation (7)); then, R 3,0 (0) is a uniform distribution from 0 to 100 s mag 2 , where s mag is the standard deviation of the time-series data.
Once the posterior distributions of coefficients for the CARMA(3, 0) model are obtained, we can determine the autocovariance function (Equation (7)) and posterior distribution of the PSD (Equation (8)). After obtaining the coefficients for the CARMA(3, 0) model and its posterior PSD, the residuals of the light curve need further analysis to evaluate the CARMA(3, 0) model. The time-series data are fitted by using a Kalman filter that estimates the mean and variance of the measured time series at time t j when the covariance of timeseries data and measurements at time t i (i<j) are given. Then we analyze the fitting residuals, including the distribution and autocorrelation function (ACF) of residuals. If the CARMA(3, 0) model is correct, the fitting residuals should have a normal distribution, and we can perform a Kolmogorov-Smirnov (K-S) test on the residuals. In addition, the residuals χ i and the squared residuals ( ( )) c c -E i i 2 should form the white-noise sequence. If the distribution of residuals deviates from the normal distribution or a large number of ACFs of residuals are outside the region, assuming a white-noise process, it means that the CARMA(3, 0) model does not capture all correlation structures in the time-series data.  The K-S test is performed on the residuals, yielding a p-value of 0.49, which is greater than the significance level (say 5%) and cannot reject the hypothesis that the residuals come from the normal distribution (orange line). The ACFs of the residuals and their square are shown in the bottom panels compared with the 95% confidence region assuming a white-noise process (gray shaded region).

Determination of Granulation Parameters
The COR function is used to fit the median value of the posterior PSD and obtain the granulation parameters τ gran , σ, and the exponent α in Equation (3), which are shown in Figure 2 for the example star. It can be seen in Figure 2 that the posterior PSD (blue line) traces well the background of the observational PSD and is barely affected by the pulsation (red plus sign) and measurement noise (black dashed-dotted line). In order to compare the characteristic timescales resulting from different α, Mathur et al. (2011) defined an effective timescale τ eff as the e-folding time of the ACF. 1 For this purpose, the derived τ gran is converted to τ eff based on the ACF. The uncertainty of the granulation parameter is calculated by the MCMC method from taking the 68% confidence interval of the posterior PSD into account.
The granulation parameters of RSGs in the SMC, LMC, and M31 are listed in Tables 1-3, and the distributions of τ eff , σ V , and σ R are shown in Figure 3. The characteristic timescale of granulations, τ eff , ranges mainly from several to a few hundred days. This is significantly longer than in red giants, where τ eff is generally not longer than a few days (Mathur et al. 2011). Similarly, the characteristic amplitude of granulation, σ gran , ranges from several to several tens of mmag, also larger than in red giants (Mathur et al. 2011). Such differences between RSGs and red giants are understandable and will be discussed later.
The granulation parameters exhibit systematic differences between the three galaxies. The median value of τ eff increases from 31.60 days in the SMC to 44.62 days in the LMC to 62.34 days in M31. An apparent reason may be the metallicity effect, since the three galaxies have very different metallicities, then τ eff increases with metallicity. The effect of metallicity on τ eff was previously observed in red giants by Corsaro et al. (2017). They found that the timescale of granulations lengthens with metallicity for the red giants in three open clusters spanning a metallicity range from [Fe/H] ; −0.09 to 0.32. Here the metallicity ranges from about −0.7 for the SMC to −0.46 for the LMC to +0.3 for M31, and τ eff changes by more than a factor of 2. Collet et al. (2007) found that increasing metallicity would enlarge T g eff and t µ T g gran eff (see Section 4.1 for details), where T eff is the effective temperature and g is the surface gravity, hence a long timescale. Basically, metallicity increases opacity, which in turn increases the mixing length l, and the granulation timescale would then become longer, since t µ l c s gran as far as the sound velocity c s remains the same. The median σ gran changes from 47.00 mmag in the SMC to 59.22 mmag in the LMC, both in the V band, and 64.03 mmag in M31 in the R band. The amplitude in the R band should be converted to that in the V band for comparison because the V-band amplitude is generally larger. There is no systematic study of the ratio of amplitude in the V and R bands (hereafter C V R ) for RSGs. Instead, we replace with the mean ratio of red       . 2 Then the median σ gran of RSGs in M31 should be 96.05 mmag in the V band. For Betelgeuse, the V-band observation is available, and the amplitude σ V is calculated to be 114 mmag (see Section 3.4 for details), which is consistent with this converted value. Consequently, the characteristic amplitude of granulation increases with metallicity, which also agrees with the discovery in red giants in three open clusters with different metallicities by Corsaro et al. (2017). Similar to the argument for the granulation timescale, metallicity increases opacity and mixing length to affect the convective dynamics (Tanner et al. 2013a). The increasing mixing length makes the granules occupy a larger surface area (Collet et al. 2007;Tanner et al. 2013b) and hence a higher characteristic amplitude of granulation (Ludwig 2006;Kjeldsen & Bedding 2011). Figure 4 shows the change of granulation effective timescale and characteristic amplitude with metallicity after excluding the outliers in the granulation parameters (see Section 4.3.1 and 4.4), where the median and 1σ values of τ and σ V are displayed for each galaxy at the given metallicity. It can be seen from Figure 4 that the granulation effective timescale and characteristic amplitude increase systematically with metallicity. Meanwhile, the dispersion is significant, which could be caused by the dispersion of metallicity of individual RSGs. Moreover, the granulation parameters depend not only on metallicity but also on stellar parameters. The influence of metallicity will be clearer after the effects of other stellar parameters such as mass and age are peeled off.

Betelgeuse
Betelgeuse is appended to the sample of RSGs in the three nearby galaxies as one of the nearest Galactic RSGs for comparison, as well as a check of the method. The timescale of granulation on Betelgeuse is derived to be several weeks to 1 yr in the optical or several years in the H band by comparing the 3D radiative hydrodynamics (RHD) simulation with interferometric observations (Chiavassa et al. 2010) and 280±160 days by fitting the ACF of the light curve (Derekas et al. 2017). By analyzing the images through spectropolarimetry, López Ariste et al. (2018) found that individual convective structures can be tracked over 1 yr, but the shape and position of the convection cells change from several weeks to several months. All of these results indicate that the evolutionary timescale of the granulation of Betelgeuse ranges from several days to 1 yr. Based on the V-band light curve from the AAVSO database, we calculated the timescale and amplitude of granulation on Betelgeuse ( Figure 5). It should be noted that the data after JD-2,458,665 in Figure 5 are removed because the recent persistent fainting of Betelgeuse (Guinan et al. 2019) is mysterious and very likely unrelated to granulation. As can be seen from Figure 6, our method by fitting the posterior PSD of the light curve yielded τ eff =138 −3 +8 days, which is consistent with previous results in the optical bands. Besides, the amplitude in 19 mmag, which agrees with the expectations as discussed above.

GP Regression
Since the granulation signals can be modeled by the autocovariance and corresponding PSD (see Section 3.2), the Gaussian process (GP) can also be introduced to investigate the granulation. A key fact of GPs is that they can be completely defined by their mean and covariance functions. So the GP regression has mathematical ideas very similar to the CARMA model and has been used to model granulation and oscillations in red giant stars ( when α=4. In order to compare the CARMA model with the CELERITE method, the PSD of the light curves is calculated by using CELERITE, and the corresponding granulation parameters are derived. The results are also listed in Tables 1-3 as a reference. Figure 7 compares the results from the CARMA and CELERITE models, which shows that σ gran agrees with each other, but τ eff presents some dispersion, especially for the RSGs in M31. This may be explained by the dependence of the accuracy of τ eff on the regular sampling, while the accuracy of σ gran does not have such dependence. In general, the two methods show consistent results, since both originated from the same mathematical idea. In this work, all further analysis is based on the results from the CARMA model.

Basic Assumptions
The relations between granulation and stellar parameters are expected from both the models and observations of convection and the resultant granulation. As verified in Huber et al. (2009) andBedding et al. (2011), the timescale of granulation is proportional to the pressure scale height and inversely proportional to the sound speed. Then, the timescale of granulation can be expressed as t µ H p (Kjeldsen & Bedding 1995). Besides, another granulation parameter σ can be expressed as s µ c n s (Kjeldsen & Bedding 2011), σ being the rms intensity fluctuation or granulation characteristic amplitude (hereafter σ gran ) and n being the number of granulations. Since µ n ( ) R H p 2 , s µ µ µ T gR T R M L MT gran eff 1.5 eff 1.5 0.5 eff 0.5 . In this work, we investigate the relations between the granulation parameters and stellar radius, effective temperature, luminosity, surface gravity, and mass.

Effective Temperature
The effective temperature is calculated from the nearinfrared intrinsic color index. A preparatory step is to transform the 2MASS color index ( )  (Russell & Dopita 1990), and 9.00 (Zaritsky et al. 1994), respectively. The symbols indicate the median value of granulation parameters in three galaxies, and the error bars show the 1σ range of the granulation parameter distributions.
standard Johnson system (Bessell & Brett 1988) following Carpenter (2001): The difference between the K and K 2MASS bands is very small and ignored. Further, the interstellar extinction is corrected to obtain the intrinsic color index ( ) -J K 0 by adopting The effective temperature T eff is then calculated by its relation with ( ) -J K 0 (Neugent et al. 2020) based on RSGs with known physical properties (Levesque et al. 2005(Levesque et al. , 2006;

Luminosity and Radius
For RSGs in nearby galaxies, Davies et al. (2013) derived a function to calculate the luminosity as follows: where m λ is the apparent magnitude in band λ and μ is the distance modulus. We choose the K band because RSGs are bright in this band and it is little affected by interstellar extinction. The distance moduli of the SMC, LMC, and M31 are taken to be 18.91 (Hilditch et al. 2005 It should be noted that the errors in L and T eff are transferred into the uncertainty of R; consequently, the error of R could be much larger than the L and T eff derived directly from observation.

Mass and Surface Gravity
The stellar mass, M, is obtained from the mass-luminosity relation ( ) Stothers & Leung 1971). The surface gravity can then be calculated by = g GM R 2 . Again, we should be careful about these secondarily derived parameters. The result depends on the accuracy of the firsthand parameter L and secondhand parameter R, as well as the mass-luminosity relation. Nevertheless, these relations have solid fundamental physics and are reasonable.

Elimination of Outliers
The CARMA(3, 0) model may identify a fake high frequency in light variation due to wrongly fitting the gaps in the light curve. Figure 8 presents an example, No. 458 in the LMC, for which the CARMA(3.0) model captures a high-frequency variation around day 500 in the gap duration. In such a case, a very short timescale of granulation would be inferred incorrectly. The granulation parameters derived from this type of data appear as outliers to the granulation-stellar parameter relations. The selection of a narrower boundary for a uniform prior would help to avoid the occurrence of extreme values, but on the other hand, this selection of a prior will force some incorrect results to fall into a reasonable range.
To reduce the effect of the outliers, the Kendall-Theil method resistant to the effects of outliers (Theil 1950;Sen 1968) is used to calculate the coefficients for robust linear regression between granulation and stellar parameters twice.   . Relation of the granulation effective timescale τ eff with stellar radius, effective temperature, surface gravity, luminosity, and mass (top to bottom) for RSGs in the LMC (red dots), compared with the RGB stars (gray dots) and Betelgeuse with error bars. The blue dashed-dotted lines are the first robust linear fit with the 95% confidence region (blue shaded region), and the black open circles are outliers beyond the 95% confidence region from the first fitting. The black solid lines are the second robust linear fit with the 95% confidence region (yellow shaded region).
The slope of the line is the median of all possible pairwise slopes between points; meanwhile, the intercept is calculated so that the line will run through the median of input data. The points beyond the 95% confidence region in the relations from the first fitting are regarded as outliers and dropped in the final linear fitting, which is marked as Y in the columns "Out-lier_taueff" and "Outlier_sigma" in Tables 1-3. Actually, since the Kendall-Theil method is resistant to the effects of outliers, there is little difference between the first and second fitting parameters. Figure 10. Same as Figure 9, but for the granulation amplitude σ gran .

Relations between τ eff and Stellar Parameters
The effective timescale τ eff of granulation is found to be tightly correlated with stellar parameters such as luminosity, effective temperature, radius, surface gravity, and mass. The relations for RSGs in the LMC are taken as an example in Figure 9, where the red dots represent the RSGs to derive the linear relation and the blue shaded region marks the 95% confidence region, while the outliers are represented by black open circles and dropped in fitting. The timescale increases with radius and luminosity and decreases with effective temperature and surface gravity, which agrees with the analysis based on the granulation assumption in Section 4.1. The timescale also increases with mass, which can be understood from t µ L MT gran eff 3.5 given that L∝M 4 . For the RSGs in the SMC and M31, the relations with the stellar parameters are analogous to the LMC with somewhat different coefficients, which are not present here.
Betelgeuse is compared by taking the effective timescale τ eff and stellar parameters from Chiavassa et al. (2009) andDerekas et al. (2017), as well as the timescale from this work. As can be seen from Figure 9, Betelgeuse lies at the high-mass, lowtemperature end, in comparison with the RSGs in the LMC. On one hand, Betelgeuse is in agreement with the relations of τ eff versus stellar parameters derived from the LMC RSGs. On the other hand, its τ eff seems to be larger than that at a given stellar parameter in the LMC. In Section 3.3, it is found that τ eff increases with metallicity. Betelgeuse is an RSG in our Galaxy with [Fe/H] ∼0.1 (Lambert et al. 1984), much more metalrich than RSGs in the LMC, and a larger τ eff is expected in agreement with what is revealed here.
Since the granulation scaling relations for RSGs have never been explored before, here we only take the red giant branch (RGB) stars for comparison. Granulation is common in red giants and considered to be the mechanism for their smallamplitude light variation. The granulation effective timescale τ eff is calculated by de Assis Peralta et al. (2018), and the stellar parameters are calculated by Yu et al. (2018) for a large sample of RGB stars (Yu et al. 2018), which are compared with the RSGs in the LMC in Figure 9. It can be seen that the relations of τ eff with stellar radius R and surface gravity log g coincide with the extrapolation of the relations for RGB stars. On luminosity, the RSGs lie slightly above but close to the extrapolation of the relation for RGB stars. Meanwhile, a very good correlation is present between τ eff and mass for RSGs with a mass range from about 10 to 20 M e , which is not visible for RGB stars within a relatively narrow range of mass. A remarkable difference occurs in the relation with effective temperature in that the relation of τ eff with T eff of RSG stars lies significantly above that of RGB stars. As mentioned above, the relation between τ eff , surface gravity, and effective temperature can be expressed as t µ T g . The difference of the surface gravity would lead to a vertical shiftg log in the t eff versus log T eff diagram.
Because the g log value is∼0.0 for RSGs and ∼2.0 for RGB stars, the τ eff of an RSG would be∼2.0 dex larger than that of an RGB star at the same T eff , which is consistent with the difference observed in Figure 9. In short, the relations of the granulation characteristic timescale τ eff are consistent with the RSG Betelgeuse and RGB stars and in agreement with the granulation theory.

Relations between σ gran and Stellar Parameters
The characteristic amplitude of granulation is expected to relate to the stellar parameters as analyzed in Section 4.1. In the same way as dealing with the parameter τ eff , we use the RSGs in the LMC as an example. The results are displayed in Figure 10. The amplitude increases with stellar radius, luminosity, and mass but decreases with effective temperature and surface gravity, which agree with the analysis based on the granulation assumption in Section 4.1. Indeed, the relations with individual stellar parameter are not very straight. For example, the relation with effective temperature can be derived from s µ L MT . Whether σ gran increases or decreases with T eff would depend on the relation of M/R with T eff . Assuming that , it can be inferred from the inverse correlation between σ gran and T eff that the exponent β<0.67.
The relations of σ gran with stellar parameters are compared with RGB stars in Figure 10. The deviation of RGB stars from the derived relation for RSGs looks much more apparent than in the case of τ eff presented in Figure 9. That RGB stars lie below the σ V versus T eff relation can be easily understood from the smaller surface gravity of RSGs, since σ gran is inversely proportional to g as µT gR eff 1.5 . The other deviations cannot be directly explained in this way, depending on the mutual relations between stellar parameters, which is beyond the scope of this paper, but we present more discussion in Section 4.4.2.

Scaling Relations between Granulation Parameters L and T eff
The brief analysis in Section 4.1 reveals that the granulation parameters τ eff and σ gran are related to multiple stellar parameters, such as L, T eff , R, log g, M, etc. Because these stellar parameters are not independent, the relations with all of these parameters would be unnecessary. Although there could be various options to combine some parameters, we decide to choose L and T eff as the principle ones because only L and T eff are derived directly from the observations, while the other parameters depend on L and T eff and would have more uncertainties (see Section 4.2). Without knowing the analytical form of the relations, a polytropic equation is adopted, i.e., t = b g L T 10 a eff eff , and σ gran takes the same form with different power-law indexes.   Figure 11. Scaling relations of the granulation effective timescale τ eff (left) and amplitude σ gran (right) with stellar luminosity and effective temperature for RSGs in the SMC (top), LMC (middle), and M31 (bottom). For RSGs in M31, the σ V of RGB stars (gray dots) and Betelgeuse (gray square) is scaled down to σ R by dividing by 1.5.
The indexes a, β, and γ of the scaling relations are determined by the MCMC method using PyMC3 (Salvatier et al. 2016). After excluding the outliers in granulation parameters, 94, 342, and 226 RSGs are contained to determine these unknowns for τ eff in the SMC, LMC, and M31, respectively, and 73, 271, and 188 RSGs for σ gran . For RSGs in the SMC, we notice that seven stars (Nos. 3839, 6999, 11765, 19164, 20019, 38653, and 43337) deviate from the scaling relation significantly. We rechecked their light curves, which were found to be inconsistent with RSGs. A further examination found that the seven objects are relatively blue and faint, with four of them being classified as Cepheid variables (Gaia Collaboration et al. 2018;Holl et al. 2018). Thus, they are discarded in the final fitting. This case analysis implies that there may be a potential to identify an RSG from the consistency with the scaling relation or not.
The fitted indexes are listed in Table 4 for τ eff and Table 5 for σ gran , together with the correlation coefficient. Figure 11 illustrates the scaling relations of the granulation parameters with stellar luminosity and effective temperature for RSGs in the SMC, LMC, and M31. The granulation amplitudes of RGB stars were derived from the Kepler or CoRoT observations, whose effective wavelength of filters is analogous to the V band; thus, their σ gran is divided by a factor of 1.5 when comparing with the σ R of RSGs in M31. It can be seen that the correlations are apparent in all cases, for both τ eff and σ gran , and for all three galaxies. Among them, the relation for τ eff in M31 is relatively loose, which can be understood by the presence of bigger gaps and less regular sampling in the time-series data from iPTF, which leads to the more uncertain determination of τ eff . Besides, the coefficients β σ and γ σ in the scaling relations for the SMC and LMC (σ V ) are larger than for M31 (σ R ). In other words, σ V is more sensitive to stellar parameters than σ R because the intensity depends on the opacity at a given wavelength. The metallicity effect is visible in comparison with the same RGB star sample. The relative position of RGB stars moves from slightly above the fitting line to beneath; i.e., τ eff and σ gran systematically rise from the SMC through the LMC to M31. This effect is explained in Section 3.3.
From the fundamental analysis in Section 4.1, In comparison with this fundamental analysis, the fitted relations with L coincide with the expectation that β ranges from about 0.4 to 0.9 for both τ eff and σ gran . On the other hand, the fitted relations with T eff show some discrepancy. The fitted value of γ for τ eff ranges from about −2 to −4, more or less in agreement with the expected −3.5, though slightly different. Meanwhile, the value of γ for σ gran ranges from about −2 to −5, significantly different from −0.5, the expected value, which means σ gran is much more sensitive to T eff than expected. In Section 4.4.2, we give a possible explanation for this large discrepancy.

Offsets between RSGs and RGB Stars
The relations of granulation parameters with individual stellar parameter exist offsets between RSGs and RGB stars, which can be seen in Figures 9 and 10. In particular, the relation with effective temperature exhibits the most significant offset. According to the analysis of Section 4.3.2, the offset of τ eff with T eff could be attributed to the difference in surface gravity. In order to clarify the effect of g, we plot τ eff versus T g eff in Figure 12(a) for the RSGs in the LMC and the RGB stars by following the fundamental relation of granulation t µ T g eff eff . The offsets between RSGs and RGB stars disappear in the τ eff versus T g eff diagram, which confirms previous judgement.
The case of the characteristic amplitude σ gran is more complicated, as discussed in Section 4.3.3. Following the basic assumption of Kjeldsen & Bedding (2011), s µ c n s , then s µ T gR gran eff 1.5 . Figure 12(b) presents the relation of σ gran with µT gR eff 1.5 for the RSGs in the LMC and the RGB stars, but the offset still exists between RSGs and RGB stars, which may indicate that the s µ T gR , which does bring about the agreement between RSGs and RGB stars. Moreover, substituting the relation between L and M yields s µ L T gran 0.75 eff 2.5 , a very different relation with T eff in comparison with previous works. This roughly coincides with the scaling relation found in Section 4.4.1, where σ gran is found to be proportional to ( ) ( ) ---T eff 2 5 .

Summary and Conclusion
The irregular variation and its relation with the granulation of RSGs are analyzed for the RSG samples in the SMC, LMC, and M31 with the time-series data from ASAS-SN and iPTF. The CARMA(3, 0) model is used to capture the correlation structures of the 128, 385, and 359 light curves of RSGs in the SMC, LMC, and M31, respectively, which forms the basis to use the COR method to fit the posterior PSD and obtain the characteristic evolution timescale τ eff and amplitude σ gran of granulations. It is found that the granulations in most of the RSGs evolve at a timescale of several days to 1 yr with a characteristic amplitude of 10-1000 mmag. Both τ eff and σ gran increase from the SMC to the LMC to M31, which may indicate the influence of metallicity on the granulation characteristics.
The relations between granulation and stellar parameters are analyzed. Both the effective timescale τ eff and amplitude σ gran of granulations are tightly correlated with stellar effective temperature, luminosity, radius, surface gravity, and mass. These relations agree with the expectations from the granulation model. In addition, they are in accord with the extrapolation of the relations derived from the Kepler light curves of RGB stars. The numerals of the Galactic RSG Betelgeuse, whose granules were detected by interferometric observations, support these relations as well. A scaling relation is derived for the granulation parameters with stellar luminosity and effective temperature for the RSG sample in the SMC, LMC, and M31.
The results of the ranges of granulation parameters and their relation with stellar parameters agree with the assumption that irregular variation of RSGs can be attributed to the evolution of granulations on their surfaces. The large granulation characteristic amplitude implies that the surfaces of RSGs are dominated by a few but huge granules.