Betelgeuse Fainter in the Submillimeter Too: An Analysis of JCMT and APEX Monitoring during the Recent Optical Minimum

We derive the SCUBA-2 point-source photometry for each observation by integrating under the point-spread function (PSF) scaled to the peak flux at sub-pixel precision. The PSF is generated by convolving the standard SCUBA-2 PSF consisting of two Gaussian components (Dempsey et al. 2013) and the 7'' smoothing Gaussian. This approach is more effective in order to measure the compact-flux component, avoiding contamination by any negative bowling effects, and by the extended circumstellar-envelope component, which methods such as aperture photometry may cause. The pre-optical dimming (2012–2013) JCMT/SCUBA-2 850 μm fluxes are consistent with those derived by O'Gorman et al. (2017) using the Atacama Large Millimeter/submillimeter Array (ALMA) main array at the same wavelength.

The JCMT/SCUBA-2 light curve, along with the archival APEX/LABOCA data, is shown in Figure 1. For comparison, the AAVSO⁹ optical (V-band) light curve for the same date range is also shown here.

Zoom In Zoom Out Reset image size

To determine whether there is any evidence of variation in the sub-mm fluxes of Betelgeuse, we perform Bayesian inference by forward-modeling the fluxes with three different models and evaluate which one best reproduces the fluxes simultaneously for SCUBA-2 850, SCUBA-2 450, and LABOCA. Because the stellar signal is strong, we assume that the uncertainties of the fluxes are dominated by the systematic, uncorrelated scatter in the determination of the absolute calibration (Dempsey et al. 2013). The models we consider are as follows.

• C1: constant flux. This model has three free parameters: c, the SCUBA-2 850 μm flux at 850 μm; α₄₅₀, the ratio between the 450 and 850 fluxes; and α₈₇₀, the ratio between the SCUBA-2 850 and LABOCA 870 μm fluxes. These α parameters allow us to straightforwardly handle the conversion from the calibration of one instrument to another, which would otherwise inject additional uncertainty on the value of c.
• C2: the behavior is broken into two distinct epochs each with constant flux. Unlike model C1, only the α parameters are assumed to be fixed at all times. There are then three other parameters, t_break, the date at which the behavior changes, and c and c', the (otherwise constant) SCUBA-2 850 μm flux before and after t_break, respectively.
• L: the flux varies linearly with time, with all bands following the slope. This is similar to model C1, except that an additional parameter m is added so that the fluxes are proportional to mt + c. This effectively assumes that the fractional rate of change at each band is the same.

In the absence of further constraints, the priors on all parameters are assumed to be flat, drawn from uniform random distributions that cover the required ranges that are listed in Table 2.

Table 2.  Results of Fitting Using Nested Sampling

Model	lnZ	α870	α450	c	c'	m	tbreak
				(Jy)	(Jy)		MJD
C1	−36.11 ± 0.06	1.23 ± 0.05	2.73 ± 0.13	0.525 ± 0.015	⋯	⋯	⋯
C2	−29.24 ± 0.07	1.15 ± 0.06	2.81 ± 0.13	0.567 ± 0.018	0.466 ± 0.021	⋯	${58320}_{-900}^{+380}$
L	−44.84 ± 0.11	1.15 ± 0.05	2.79 ± 0.13	2.22 ± 0.41	⋯	$(-2.96\pm 0.72)\times {10}^{-5}$	⋯
Priors	⋯	0–2	0–10	C1, C2: 0–1	0–1	−3–3	56000–59000
				L: 0–100

Download table as:  ASCIITypeset image

We now have a model-selection problem including multiple comparisons. To identify the most appropriate model of the three, we must compute Bayes' factors (K) for each comparison. The Bayes factor provides a way to select between families of models by considering the total evidence (or integrated likelihood) over the whole prior volume and comparing the ratio of evidences for different families of models. This requires that we have a reliable way of estimating the evidence, Z for each model (for further details on Bayes' factors, see e.g., Goodman 1999a, 1999b; Morey et al. 2016). To do this, we use the Python nested-sampling package dynesty, which supports dynamic nested sampling, allowing optimization for either evidence estimation or posterior exploration (Skilling 2004, 2006; Higson et al. 2019; Speagle 2020). For each model, we conduct two runs, one optimized to estimate the evidence for the model selection, and the other optimized to evaluate the posterior and hence provide parameter estimates. All runs used 1500 initial live points and dlogz_init = 0.01. The estimated evidence and parameter values for each model are given in Table 2, along with their respective 1σ uncertainties.

Using the evidence computed using dynesty, we now compute Bayes' factors for each pair of models. Unlike frequentist approaches, Bayes' factors automatically penalize models with more free parameters but do not privilege any model by considering it the null hypothesis. By assuming that all models are a priori equally likely, the Bayes factor collapses to the ratio of the evidences—the impact of this assumption is difficult to quantify, but this is a necessary starting point. The frequentist equivalent to the process would be to separately compare model C2 and L to C1, considering model C1 as the null hypothesis. However, this would require an additional penalty as we are performing multiple comparisons, while the use of Bayes factors avoids this issue.

The ratios of the evidence support model C2 over the other two models by a significant margin. For the comparison between C2 and C1, the Bayes' factor K = 710, while the K comparison between C2 and L is ∼6 × 10⁶. This represents decisive evidence in favor of model C2 (using the scale of, e.g., Kass & Raftery 1995), suggesting not only that the sub-mm flux is variable but, crucially, that it is lower during the recent minimum than it was before. Given that the sub-mm emission is dominated by the star, this implies that the luminosity has decreased by nearly 20% (a 3.6σ change) and argues against models in which the recent optical dimming is caused by the formation of a new cloud of dust along the line of sight to the star, as increases in the dust mass cannot decrease the sub-mm flux as the dust emission at this wavelength is optically thin.¹⁰

To understand the role of dust versus the photosphere in the dimming across the optical–sub-mm wavelength range, we perform dust radiative-transfer modeling. We do not aim to quantitatively fit the observations, merely to qualitatively illustrate the trends that could be expected under different assumptions. Using the Monte Carlo radiative-transfer code hyperion (Robitaille 2011), we compute models of a spherically symmetric wind with constant outflow speed of 14 km s⁻¹ (comparable to Loup et al. 1993). We use dust that consists of compact spherical grains composed of 30% alumina using optical constants from Begemann et al. (1997) and 70% Mg–Fe silicate for optical constants from Dorschner et al. (1995), with grain sizes from 0.1 to 1 μm following an MRN size distribution (Mathis et al. 1977). We use one million photons per iteration, which is sufficient as the shell is optically thin and spherically symmetric, with grid cells distributed logarithmically in radius. Using a dust mass-loss rate of 3 × 10⁻⁹ M_⊙ yr⁻¹ (e.g., Verhoelst et al. 2009), we achieve an acceptable fit to the optical–mid-infrared photometry (see Table A1 in the Appendix). We then compute a second model, in which an additional shell of dust has been added between 2 and 4 R_*, close to the observed dust-condensation radius (Haubois et al. 2019), which adds an additional 1 mag of optical extinction.

Comparing these two models (see Figure 2), we see two distinctive features. First, while the addition of new dust substantially changes the optical and infrared emission, it makes very little difference at wavelengths longer than 100 μm. Second, Betelgeuse's emission at sub-mm wavelengths is entirely dominated by the stellar photosphere, with small contributions from circumstellar dust and the radio pseudo-photosphere (i.e., optically thick excess radio emission whose size changes with wavelength, not included in our models; e.g., Richards et al. 2013; O'Gorman et al. 2017). The radio pseudo-photosphere, which dominates the spectrum at cm wavelengths and beyond (O'Gorman et al. 2015), clearly makes only a minor contribution at 850 μm and a negligible one at 450 μm. Therefore, while we are unable to probe changes in dust mass at the level that would produce the optical dimming, it is clear that only changes in the photosphere can reproduce the sub-mm dimming—this would otherwise require that the dust emission or the radio pseudo-photosphere had completely disappeared. Using the same logic, we can exclude changes in the sub-mm line emission as a cause of the variations. Similar to CW Leo (Dharmawardena et al. 2019b), the CO contribution (the brightest line in these spectral regions) to the measured continuum is at the 5%–10% level; short of all the molecular line emission disappearing, the continuum could not change by this amount.

Zoom In Zoom Out Reset image size

Under the naïve assumption that the change in sub-mm flux corresponds to a change in temperature at constant radius, and that the original temperature was T_eff = 3650 K (the 2004 temperature from Levesque & Massey 2020), this would suggest a present value of 3450 K. Alternatively, this could be reconciled with the sub-mm data if starspots with T_eff = 3250 K cover 50% of the visible surface, or if spots with T_eff = 3350 K cover 70% of the visible surface. Indeed, as can be seen in Figure 3, all such scenarios would also be accompanied by a ∼0.9 mag dimming in the V-band, and are qualitatively similar to the 2020 spectrum of Levesque & Massey (2020), including the convergence of the spectra at the blue end with the warmer/2004 spectrum. However, as a dynamic, non-equilibrium system, changes in temperature and the resulting molecular bands in pulsating stars may take some time to settle after a significant perturbation (e.g., McDonald & van Loon 2007; Lebzelter et al. 2010, 2014). Hence, 1D models have difficulty capturing the observed behavior and detailed 3D modeling is required. This might explain why the expected changes in TiO bands were not visible in the spectra of Levesque & Massey (2020), as it is typical to compare to static, 1D models. If the change were instead due to a change in radius at constant temperature, the change in radius would be small (∼10%), as ${\rm{\Delta }}R\propto \sqrt{{\rm{\Delta }}L}$.

Zoom In Zoom Out Reset image size