On Granulation and Irregular Variation of Red Supergiants

The mechanism and characteristics of the irregular variations of red supergiants (RSGs) are studied based on the RSG samples in Small Magellanic Cloud (SMC), Large Magellanic Cloud (LMC) and M31. With the time-series data from All-Sky Automated Survey for SuperNovae (ASAS-SN) 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. The characteristic evolution timescale and amplitude of granulations are further derived from fitting the posterior power spectral density 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 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 the 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 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 visual and infrared bands (Levesque et al. 2007;Yang et al. 2018). According to the characteristics of light curves, RSGs are divided into three types: semi-regular variables, variables with long secondary period (LSP), and irregular variables (Yang & Jiang 2011, 2012Ren et al. 2019). The irregular variation, present in all the 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 to characterize 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 firstly 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 three dimensional simulation (Freytag et al. 2002), and that the convective outer layers of RSGs cause both the spectra 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) find a strong 1/f noise component in almost all the stars in their sample, which suggests that convection play an important role in all the types of RSGs. With increasing power of convection, the light variation of a RSG becomes more and more irregular.
Betelgeuse (α Ori), as the nearest (∼ 130-200 pc, see van Leeuwen (2007) and Harper et al. (2008)) and brightest red supergiant, has been studied in details 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 ∼ 30 mas (i.e. about 200-1200 R at a distance of 200pc), which implies small number of granulations if considering the diameter of Betelgeuse to be about 1800 R . Other works (Monnier 2003;Haubois et al. 2009;Ohnaka et al. 2009Ohnaka et al. , 2011Kiss et al. 2010) also found the evidences for the existence of large convective cells in Betelgeuse.
There are only a few RSGs like Betelgeuse which can be resolved by interferometric observation thanks to 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 a 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 time scale of convective motions in Betelgeuse derived from line bisector variation (Gray 2008). These studies justified the possibility to estimate 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 Cepheid lies in the lack of the Kepler observation, i.e. long-term continuous and accurate photometry. However, various time-domain surveys have collected large datasets, in particular, toward the nearby galaxies where RSGs are detectable for their high luminosity. In addition, extra-Galactic 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 measurement of stellar distance and parameters. Following our previous study of the RSG light variation, the nearby galaxies, namely LMC, SMC, M31 and M33, are chosen as the targets. On the light variation datasets, we choose the All-Sky Automated Survey for SuperNovae (ASAS-SN, Shappee et al. 2014;Kochanek et al. 2017) that covers the Small Magellanic Cloud (SMC) and the Large Magellanic Cloud (LMC), and the Intermediate Palomar Transient Factory (iPTF) survey Rau et al. 2009) that observed M31 and M33. We first describe 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 near-infrared brightness and colors, RSGs are selected in the nearby galaxies. Historically there have been many collections of RSG samples in the MCs (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 LMC and SMC respectively by combining most of the available infrared observations and checking carefully the brightness and colors, which is the largest RSG sample in the MCs up 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) find that the minimum luminosity of RSGs in M31 and M33 is about one magnitude above the theoretical limit of RSGs corresponding to 9M . This is understandable because the observations are limited in particular to 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 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 survey for RSGs in LMC and SMC, and from the iPTF survey for RSGs in M31 and M33. The ASAS-SN consists of 24 telescopes all around the globe, which observed variable objects in SMC and LMC at a cadence of 4-5 days for 1600 days approximately (Shappee et al. 2014;Kochanek et al. 2017). 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 g-and R-band, reaching 20.5 mag at 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 leads to 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 one third 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 -band or R-band due to the red color of RSG and the facility sensitivity skewed to longer wavelength, therefore the time-series data to be used are from the V -band observation for 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 multi-band photometry. We take the Two Micron All Sky Survey (2MASS) (Skrutskie et al. 2006) measurement in the J, H, 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 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 SMC, LMC and M31, respectively. Since the granulation signal could be detected in all categories of RSGs, all the RSGs with sufficient time-series data are analyzed further, no matter whether they have been classified into semi-regular variables, irregular variables or variables with LSP (Yang & Jiang 2011, 2012Ren et al. 2019).

The granulation model
The major parameters to characterize the granulations are the timescale τ gran and the amplitude σ gran . The common way to obtain the granulation parameters uses the background of power spectrum. Harvey (1985) constructed a simple but reasonable model to describe the granulations which 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 (PSD), such as OCT , SYD (Huber et al. 2009), CAN (Kallinger et al. 2010), and A2Z (Mathur et al. 2010). Mathur et al. (2011) and Kallinger et al. (2014) used the the Harvey-like COR function (Mosser & Appourchaux 2009;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 positive parameter to be fitted that characterizes the slope of the decay. For various α, a normalisation factor is required for Eq. 2, thus we change Eq. 2 to where ξ is a normalisation factor that depends on the value of α. Moreover, from Karoff et al. (2013), ξ is related to α as ξ = 2α sin(π/α).

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 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, i.e. the Continuous-time AutoRegressive Moving Average (CARMA for short) model to estimate the variability features of the light curves and their PSD. CARMA(p, q) process has been used to model X-ray variability of AGNs (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, α 0 , . . . , α p−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 1 , ..., r p of the autoregressive polynomial have negative real parts. The autocovariance function of CARMA(p,q) model at lag t is: 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 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 1992). 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 AIC for small samples (hereafter AICc). In our work, the CARMA(p, q) model that minimizes the AICc yields p = 3, q = 0. Afterwards, the CARMA(3,0) parameters and the posterior PSD are calculated as the product of the likelihood function with a prior probability distribution using Markov Chain Monte Carlo (MCMC).
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 which is also known as the Ornstein-Uhlenbeck (O-U) process. The autocovariance function R 1,0 (t) and the PSD P (ν) of an O-U process are: As can be seen, P (ν) is a Lorentzian function, and R 1,0 (t) decays with a timescale τ = 1/α 0 , which are 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 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-day to 2000-day.
For characteristic amplitude σ, we set a uniform distribution on autocovariance at lag t = 0 (see Eq. 7), then R 3,0 (0) is a uniform distribution from 0 to 100 σ 2 mag , where σ 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 (Eq. 7) and posterior distribution of the PSD (Eq. 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 is fitted by using Kalman filter that estimates the mean and variance of the measured time series at time t j when the covariance of time-series data and measurements at time t i (i < j) are given. Then we analyze the fitting residuals, including the distribution of residuals and the autocorrelation function (ACF) of residuals. If the CARMA(3,0) model is correct, then the fitting residuals should have a normal distribution and we can perform Kolmogorov-Smirnov (KS) test on the residuals. In addition, the residuals χ i and the squared residuals (χ i − E(χ i )) 2 should form the white-noise sequence. If the residuals distribution 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 do not capture all correlation structures in the time series data.

Determination of Granulation Parameters
The COR function is used to fit the median value of posterior PSD and obtain granulation parameters: τ gran , σ, and the exponent α in Eq. 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 observational PSD and is barely affected by the pulsation (red cross) and measurement noise (black dash-dot line). In order to compare characteristic timescale resulted 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 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 SMC, LMC and M31 are listed in Table 1, 2 and 3 and the distributions of τ eff , σ V , σ R are shown in Figure 3. The characteristic timescale of granulations, τ eff , ranges mainly from several days to a few hundred. 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 mmag to several tens mmag, also larger than in red giants (Mathur et al. 2011). Such difference between red supergiants and red giants are understandable and will be discussed later.
The granulation parameters exhibit systematic difference between the three galaxies. The median value of τ eff increases from 31.60 days in SMC through 44.62 days in LMC to 62.34 days in M31. An apparent reason may be the metallicity effect since the three galaxies have very different metallicity, then τ eff increases with metallicity. The effect of metallicity on τ eff is previously observed in red giants by Corsaro et al. (2017). They found that the time scale 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 SMC, through -0.46 for LMC to +0.3 for M31, τ eff changes by more than a factor of two. Collet et al. (2007) found that increasing metallicity would enlarge √ T eff /g, and τ gran ∝ √ T eff /g (see Section 4.1 for details), where T eff is effective temperature and g is 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 τ gran ∝ l/c s as far as the sound velocity c s remains the same.
The median σ gran changes from 47.00 mmag in SMC to 59.22 mmag in 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 V band and R band (hereafter C V /R ) for RSGs. Instead, we replace with the mean ratio of red giants from a relatively large sample by Percy et al. (2009), i.e. C V /R = 1.5 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 metallicity 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 occupying 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 (c.f. Section 4.3.1 and 4.4), where the median and 1-sigma 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 depends not only on metallicity, but also on stellar parameters. The influence of metallicity would be more clear 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 RSG for comparison as well as a check of the method. The timescale of granulation on Betelgeuse is derived to be several weeks to one year in 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) find that individual convective structures can be tracked over one year but the shape and position of the convection cells change from several weeks to several months. All of these results indicate that the evolutionary time scale of the granulation of Betelgeuse ranges from several days to one year. Based on the V band light curve from the AAVSO database, we calculated the timescale and the amplitude of granulation on Betelgeuse ( Figure 5). It should be noted that the data after JD-2458665 in Figure 5 is removed because the recent persistent fainting of Betelgeuse (Guinan et al. 2019) is mysterious and very likely unrelated to granulation. Our method by fitting the posterior PSD of light curve yielded τ eff = 138 +8 −3 days, which is consistent with previous results in the optical bands. Besides, the amplitude in the V band σ V = 114 +19 −1 mmag that agrees with expectation as discussed above.

Gaussian Process Regression
Since the granulation signals can be modeled by the autocovariance and corresponding PSD (see 3.2), the Gaussian Process (GP) can be also 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 (Pereira et al. 2019). In the resultant PSD becomes Since Eq. 12 corresponds to the PSD of the kernel in Eq. 11 and is the same as Eq. 3 with α = 4, this kernel can be used to capture the granulation signal in the time domain. After applying a normalization factor K = 2 √ 2π (Pereira et al. 2019) to Eq. 12 and comparing with Eq. 3, we can derive τ = 1/ω 0 and σ = S 0 ω 0 / √ 2 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 Table 1, 2 and 3 as a reference. Figure 8 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 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 the further analysis are 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) and Bedding 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 (Kjeldsen & Bedding 1995). Besides, another granulation parameter σ can be expressed as σ ∝ c s / √ n (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 , σ gran ∝ T 1.5 eff /gR ∝ T 1.5 eff R/M ∝ L 0.5 /M T 0.5 eff . 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 near-infrared intrinsic color index. A preparatory step is to transform the 2MASS color index (J − K) 2MASS into J − K in the 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 foreground interstellar extinction is corrected to obtain the intrinsic color index (J − K) 0 by adopting A K /A V = 0.12 (Cardelli et al. 1989;Wang & Jiang 2014) and E(J − K) = A V /5.79 (Schlegel et al. 1998), in which A V = 0.297 mag for RSGs in SMC and LMC (Zaritsky et al. 2002(Zaritsky et al. , 2004 and A V = 1 mag for RSGs in M31 (Massey & Evans 2016). 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. , 2006Massey et al. 2009):

Luminosity and radius
For RSGs in nearby galaxies, Davies et al. (2013) derived a function to calculate the luminosity as following: 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 this band is little affected by interstellar extinction. The distance modulus of SMC, LMC and M31 are taken to be 18.91 (Hilditch et al. 2005), 18.41 (Macri et al. 2006) and 24.40 (Perina et al. 2009) respectively. For RSGs in SMC and LMC, Davies et al. (2013) obtained a = 0.90 ± 0.11 and b = −0.40 ± 0.01 in Eq. 15 for the extinction-corrected K band magnitude. For RSGs in M31 which were not studied by Davies et al. (2013), we use bolometric correction (Neugent et al. 2020) to calculate absolute magnitude and luminosity: Thus, With luminosity L and T eff calculated, the radius of RSGs can be derived from R = (4πσ) −0.5 L 0.5 T −2 eff . 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 L and T eff derived directly from observation.

Mass and Surface Gravity
The stellar mass, M , is obtained from the mass-luminosity relations L/L = (M/M ) γ where γ ≈ 4 (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 first-hand parameter L, the second-hand 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 faked high-frequency in light variation due to wrongly fitting the gaps in the light curve. Figure 7 presents an example, No.458 in LMC for which the CARMA(3.0) model captures a high-frequency variation around day 500 in the gap duration. In such 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 parameters relations. The selection of more narrow boundary for uniform prior would help to avoid the occurrence of extreme values, but on the other hand, this selection of 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. 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 'Outlier taueff ' and 'Outlier sigma ' in Table 1, 2 and 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.

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 LMC is 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, while decreases with effective temperature and surface gravity, which agree with the analysis based on the granulation assumption in Section 4.1. The timescale also increases with mass, which can be understood from τ gran ∝ L/M T 3.5 eff given that L ∝ M 4 . For the RSGs in SMC and M31, the relations with 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 Derekas et al. (2017) and Chiavassa et al. (2009) as well as the timescale from this work. As can be seen from Figure 9, Betelgeuse lies at the high-mass low-temperature end in comparison with the RSGs in LMC. On one hand, Betelgeuse is in agreement with the relations of τ eff vs. stellar parameters derived from the LMC RSGs. On the other hand, its τ eff seems to be larger than that at given stellar parameter in LMC. In Section 3.3, it is found that τ eff increases with metallicity. Betelgeuse is a RSG in our Galaxy with [Fe/H] ∼ 0.1 (Lambert et al. 1984) much more metal-rich than RSGs in LMC, 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 RGBs for comparison. Granulation is common in red giants and considered to be the mechanism for their small-amplitude light variation. The granulation effective timescale τ eff is calculated by de Assis Peralta et al. (2018) and 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 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 lies slightly above but close to the extrapolation of the relation for RGBs. Meanwhile, a very good correlation is present between τ eff and mass for RSGs with mass range from about 10 M to 20 M , which is not visible for RGBs 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 RGBs. As mentioned above, the relation between τ eff , surface gravity and effective temperature can be expressed as τ eff ∝ √ T eff /g, thus log τ eff ∝ 0.5 × log T eff − log g. The difference of the surface gravity would lead to a vertical shift -log g in the τ eff vs. log T eff diagram. Because the log g value is ∼0.0 for RSGs while ∼ 2.0 for RGBs, τ eff of a RSG would be ∼2.0 dex larger than of a RGB star at the same T eff , which is consistent with the difference observed in Figure 9. In a word, the relations of the granulation characteristic timescale τ eff are consistent with the red supergiant 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 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 σ gran ∝ L 0.5 /M T 0.5 eff . With L ∝ R 2 T 4.0 eff , this expression can be written as σ gran ∝ RT 1.5 eff /M . Whether σ gran increases or decreases with T eff would depend on the relation of M/R with T eff . Assuming that T eff ∝ (M/R) β , 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 RGBs from the derived relation for RSGs looks much more apparent than in the case of τ eff present in Figure 9. That RGB stars lie below the σ V vs. T eff relation can be easily understood from the smaller surface gravity of RSGs since σ gran inversely proportion to g as ∝ T 1.5 eff /gR. 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 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. τ eff = 10 a L β T γ eff , and σ gran takes the same form with different power law indexes.
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 SMC, LMC and M31 respectively, and 73, 271 and 188 RSGs for σ gran . For RSGs in SMC, we notice that 7 stars (No. 3839, 6999, 11765, 19164, 20019, 38653, 43337) deviate from the scaling relation significantly. We re-checked their light curves which were found to be inconsistent with RSGs. A further examination found that the 7 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 final fitting. This case analysis implies that there may be a potential to identify a 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 SMC, LMC and M31 respectively. The granulation amplitudes of RGBs were derived from the Kepler or CoRoT observation 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 σ R of RSGs in M31. It can be seen that the correlations are apparent in all the cases, for both τ eff and σ gran , and for all the three galaxies. Among them, the relation for τ eff in M31 is relatively loose, which can be understood by 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 β σ , γ σ in the scaling relations for 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 given wavelength. The metallicity effect is visible in comparison with the same RGB stars sample. The relative position of RGBs moves from slightly above the fitting line to beneath, i.e. τ eff and σ gran systematically rise from SMC through LMC to M31. This effect is explained in Section 3.3.
From the fundamental analysis in Section 4.1, τ eff ∝ L/(M T 3.5 eff ) and σ ∝ L/(M T 0.5 eff ), then the relations can be expressed with only L and T eff if M is approximated by L 1/4 , i.e. τ eff ∝ L 0.75 T −3.5 eff and σ ∝ L 0.75 T −0.5 eff . In comparison with this fundamental analysis, the fitted relations with L coincide with the expectation in 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 single stellar parameter exist offsets between RSGs and RGB stars, which can be seen in Figure 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 eff /g in Figure 12 (a) for the RSGs in LMC and the RGB stars by following the fundamental relation of granulation τ eff ∝ √ T eff /g . The offsets between RSGs and RGBs disappear in the τ eff versus √ T eff /g 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), σ ∝ c s / √ n, then σ gran ∝ T 1.5 eff /gR. Figure 12 (b) presents the relation of σ gran with ∝ T 1.5 eff /gR for the RSGs in LMC and the RGB stars, but the offset still exists between RSGs and RGBs, which may indicate the σ gran ∝ T 1.5 eff /gR relation is incorrect. From the very fundamental concept of granulation, the characteristic amplitude is an indicator of the granulation cell energy, which should be proportional to the mixing length of convection and the sound speed. Taking that mixing length is on the order of pressure scale height, then σ gran ∝ c s H p ∝ T 1.5 eff /g ∝ T 1.5 eff R 2 /M ∝ L/T 2.5 eff M . Figure 12 (c) plot the relation between σ gran and T 1.5 eff /g, which does bring about the agreement between RSGs and RGBs. Moreover, substituting the relation between L and M yields σ gran ∝ L 0.75 /T 2.5 eff , a very different relation with T eff in comparison with previous works. This roughly coincides with the scaling relation found in previous section 4.4.1, where σ gran is found to be proportional to T (−2)−(−5) eff .

SUMMARY AND CONCLUSION
The irregular variation and its relation with granulation of RSGs are analyzed for the RSG samples in 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, 359 light curves of RSGs in 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 a year with the characteristic amplitude of 10-1000 mmag. Both τ eff and σ gran increases from SMC through 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 accord with the extrapolation of the relations derived from the Kepler light curves of red giant branch stars. The numerals of the Galactic red supergiant 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 SMC, LMC and M31 respectively.
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 surface. The large granulation characteristic amplitude implies that the surface of RSGs are dominated by a few but huge granules.

ACKNOWLEDGMENTS
This paper is published to commemorate the 60th anniversary of the Department of Astronomy, Beijing Normal University. We are grateful to Drs. Ming Yang and Yong Shi, and Ms. Yuxi Wang for very helpful discussions, and the anonymous referee for very good suggestions. This work is supported by the National Natural Science Foundation of China through the project NSFC 11533002. This work has made use of data from the surveys by ASAS-SN and iPTF as well as the AAVSO database.  The K-S test is performed on the residuals, yielding the 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 lower panels compared with the 95% confidence region assuming a white noise process (gray region).   τ eff (days) Figure 9. Relation of the granulation effective timescale τ eff with stellar radius, effective temperature, surface gravity, luminosity and mass (from top panel to bottom panel) for RSGs in LMC (red dots), compared with the RGB stars (gray dots) and Betelgeuse with errorbar. The blue dashdot lines are the first robust linear fit with 95% confidence region (blue region), and the black circles are outliers beyond the 95% confidence region from the first fitting. The black solid lines are the second robust linear fit with 95% confidence region (yellow region).