The Great Dimming of Betelgeuse seen by the Himawari-8 meteorological satellite

Betelgeuse, one of the most studied red supergiant stars, dimmed in the optical by ~1.2 mag between late 2019 and early 2020, reaching an historical minimum called"the Great Dimming."Thanks to enormous observational effort to date, two hypotheses remain that can explain the Dimming: a decrease in the effective temperature and an enhancement of the extinction caused by newly produced circumstellar dust. However, the lack of multi-wavelength monitoring observations, especially in the mid infrared where emission from circumstellar dust can be detected, has prevented us from closely examining these hypotheses. Here we present 4.5-year, 16-band photometry of Betelgeuse between 2017-2021 in the 0.45-13.5 micron wavelength range making use of images taken by the Himawari-8 geostationary meteorological satellite. By examining the optical and near-infrared light curves, we show that both a decreased effective temperature and increased dust extinction may have contributed by almost the same amount to the Great Dimming. Moreover, using the mid-infrared light curves, we find that the enhanced circumstellar extinction actually contributed to the Dimming. Thus, the Dimming event of Betelgeuse provides us an opportunity to examine the mechanism responsible for the mass loss of red supergiants, which affects the fate of massive stars as supernovae.

Examining the images taken by the Himawari-8 satellite, we found that several bright stars, including Betelgeuse, sometimes appear in the images (Extended Data Fig. 2). We therefore measured the amount of light coming from Betelgeuse, typically once per 1.72 days, between January 2017 and June 2021, and we made a 4.5-year catalog of light curves in 16 bands covering 0.45-13.5 µm (Extended Data Fig. 3a). Comparing the light curves, i.e., the time series of spectral energy distributions (SEDs), with a model SED for Betelgeuse, we determined the time variation of the main parameters of Betelgeuse, namely the radius R, effective temperature T eff , and extinction A(V ) in the optical. In addition to these three parameters, which can be determined from optical or near-infrared observation 8,18 , we also succeeded in determining the dust optical depth τ10 at 10 µm using our mid-infrared light curves. With this τ10, we will examine the dust-production history around Betelgeuse. Fig. 2 shows the main results of this paper, namely, the time variations of the stellar parameters, R, T eff , A(V ), and τ10 together with the AB magnitudes for Band 3 at 640 nm, which is located at around the central part of the Johnson's R-band. From the time variations of R, T eff , and A(V ), we found that the Great Dimming of Betelgeuse, the drop of ∼1.2 mag in the V magnitude 4 , is likely to be caused by a cooling of ∼140 K in T eff which results in a dimming of the V magnitude by ∼0.6 mag, and by an increase of ∼0.6 mag in A(V ) which, by definition, darken the V magnitude by ∼0.6 mag. This result supports a previously suggested scenario in which the Dimming was caused by a combination of a decrease in T eff and an increase in A(V ) 1,8 , although several shortcomings in our model SED should be kept in mind; some of the shortcomings are examined in the Methods section and Extended Data Fig. 4.

Time variations of the photospheric parameters
However, the cause of the enhanced extinction A(V ) during the Great Dimming is still under debate. Several works have proposed that occultation by a newly formed circumstellar dust cloud might be a promising cause 1,8 . In contrast to this dust-formation scenario, several other works have argued that the dust formation is unnecessary to explain the Dimming. For example, three-dimensional (3D) effects 7 or molecular opacity 19,20 , rather than dust extinction, can also explain the Dimming. In order to examine whether or not the enhanced extinction A(V ) during the Dimming is due to the circumstellar dust, we focus on the time variation of τ10. Since τ10 was determined using the mid-infrared light curves of Bands 12-14 at around 10 µm, where O-rich dust emissions dominate the flux 21,22 , it reflects the amount of circumstellar dust that emits the infrared light. Therefore, using our infrared τ10 measurements, in contrast to optical diagnostics, we can directly trace the amount of circumstellar dust. Our measurements in the bottom panel of Fig. 2 shows an enhancement of τ10 during the Great Dimming, though there is a slight time delay between the peaks of A(V ) and τ10. This enhancement of τ10 indicates that a clump of gas did in fact produce dust that shaded the photosphere of Betelgeuse and contributed to the Great Dimming.
Two potential dust reservoirs have been found so far for red supergiants. The first is a diffuse and spread-out dust layer located a few to several tens of a stellar radius from the photosphere 21,22 . In addition, recent very-high-spatial-resolution instruments have found circumstellar dust clouds around red supergiants just a few tenths of stellar radii above the photosphere 23,24 . The distribution and amount of dust in the spread-out layer remain mostly unchanged with time 22,25 , and thus it is unlikely that dust condensation in this layer is related to the large time variation of A(V ) and τ10 during the Dimming that we found. In contrast, the latter, lukewarm region 26 coexists with the warm chromosphere, the temperature of which 9,27 is sufficiently high for dust subli- b-e, The results of the SED fitting. b, Radius, R. c, Effective temperature, T eff . d, Dust extinction in the V band, A(V ). e, Dust optical depth at 10 µm, τ 10 . The time of the Great Dimming is indicated by the red vertical line. mation 28 . Therefore, a dust clump that condenses near the photosphere may soon be sublimated by chromospheric heating, or, if the radiative acceleration is insufficient, it may fall back to the photosphere to be sublimated 29 . Considering these points, we conclude that the enhance-ments of A(V ) and τ10 during the Dimming may have occurred very close to the photosphere. This hypothesis is also supported by the fact that the polarized brightness of this circumstellar envelope changed significantly around the Dimming 30,31 .

Transition before the Dimming seen with H 2 O
In order to further examine the dust-condensation process, we here focus on the time variation of the amount of molecular gas around Betelgeuse, in which the dust condensation occurred. We here focus on H2O molecules. In addition to the dust emission seen in the wavelengths around 10 µm, the amount of molecules, particularly H2O, can also be traced using Himawari-8 photometric data. Since H2O molecules in the Earth's atmosphere absorb starlight 32 in the wavelength ranges where the H2O molecular features coming from stars exist 33 , the observation of extraterrestrial H2O molecules usually requires space telescopes 15 . One of the missions of the Himawari-8 satellite is to observe the distribution of water vapor in the Earth's upper atmosphere 10 , and thus Himawari-8 has the wavelength bands where H2O features exist. These bands, Bands 8-10 at around 6-8 µm, are therefore useful for time-series observations of H2O molecules around Betelgeuse. We focus here on the flux variation of Band 8, where the signal-to-noise ratio of the H2O feature is maximized, thanks to the strength of this feature 33,34 and the sensitivity of the Himawari-8 satellite (Extended Data Fig. 1). We note that the 6-8 µm wavelength range is a complicated region in which other molecular bands are present, e.g., SiO molecules. Thus, we should keep in mind the possibility of contamination by these molecules, and further spectroscopic observations are required to fully disentangle the problem. Nevertheless, Band 8 is the least contaminated by other molecules 15,20,34 . Fig. 3 shows the time variation of the flux excess in Band 8, defined as the ratio of the observed flux to predicted photospheric flux minus 1, which is likely to trace the amount of H2O molecular gas around Betelgeuse, maybe in the MOLsphere 15,34 . When gas clumps outside the photosphere are present that emit or absorb H2O molecular features, the excess is positive or negative, respectively. In the first half of the observing period, i.e., before September 2018, to within the error bars the flux excess was almost consistent with zero, which indicates that there was no extra detectable H2O gas around the photosphere. In contrast, most of the data points after September 2018 show negative flux excesses. This indicates the emergence of a clump and/or a layer of cool gas containing H2O molecules, in which the dust that dimmed Betelgeuse may have formed. Considering the clumpy nature of the dust cloud that dimmed Betelgeuse 1 and the existence of the warm gas clump found from UV observations of the Mg II h and k lines around October 2019 9 , the H2O gas that we found is also likely to be a clump.
Finally, and more interestingly, we tentatively identified a rapid transition in the H2O feature from absorption to emission in early April 2019 (the green vertical line in Fig. 3). After visual inspection of the unbinned light curve for Band 8, we found that the time-scale of this transition is likely to be about a week or less. This time-scale is much shorter than either the theoretical hydrodynamic time-scale of a few weeks to a few months 35 or to previously observed time-scales of changes in the H2O features seen in evolved red stars 36,37 . Because the gas clump cannot move outside of the line-of-sight to the photosphere on such a short time-scale, this transition indicates the occurrence of an episodic bursty event. One may speculate that the cause of this unusual event is related to another rare event, the Great Dimming; if so, it would be easy to imagine that such an unusual event could occur. One conceivable scenario involves a shock wave passing through the clump. It has been suggested that such a shock emerged from the bottom of the photosphere in January 2019 and propagated from there to the outer envelope of Betelgeuse 19 . Together with the fact that A(V ) began to increase at the time of this transition (the fourth panel of  might be related to the triggering process of the Great Dimming. Future theoretical investigations and time series of spatially resolved infrared observations would enable us to attack this further arising mystery.

METHODS
In this paper, we first aim to establish and test a brand-new concept, using meteorological satellites as astronomical space telescopes, to obtain multi-band, long-term, and continuous light curves uninterrupted by the Sun. Second, we applied this technique to Betelgeuse to investigate its mysterious Great Dimming.
Images of celestial objects, particularly the Moon, obtained by meteorological satellites have long been used to calibrate and ensure the quality of Earth observations 38,39,40 . In contrast, to our knowledge, this work is the first attempt to use these images to investigate stellar astrophysics. Here, we thoroughly describe the data properties, photometric methods, and their validation.
Basic property of Himawari-8's scanning observation Geostationary meteorological satellites, such as Himawari-8, observe the entire disk of the Earth by sequential scans with several rows, each of which is termed a swath 10 41 . Each swath is scanned with a 1D-like detector array from west to east. Due to this scanning procedure, the properties of images obtained by Himawari-8 have several differences compared to those obtained by typical astronomical observations using 2D devices. Nevertheless, we can use the reduced and distributed data, termed the L1B product, as usual astronomical images, provided that several caveats are taken into consideration.
First, pixel navigation and data sampling are well controlled and regular, and their errors are usually smaller than 0.5 pix for Himawari-8 42 . Thus, we can precisely transform the position in the image onto its location in the sky.
Second, the fill factor of the pixels in the north-south direction is almost 100%, and the east-west sampling interval is smaller than the pixel size 41 . As such, we do not need to be concerned with undersampling due to gaps between pixels. Moreover, after the resampling process during the reduction described below, each image is transformed to the fixed square grid, which is very easy to use.
Third, because each row of un-resampled images is observed by a single detector, inaccurate calibration of the detector occasionally introduces a stripe pattern into the obtained images. The strength of such a stripe has been shown to be small for observations of Earth scenes 43 , but it has not yet been closely examined for faint sources like stars. Therefore, it should be kept in mind that few flux measurements of stars are possibly influenced by this relatively large systematic error.
Finally, other possibly minor effects include quantization noise, nonlinearity of the detector response, and frequency of the calibration. We will further examine some causes of systematic and/or statistical errors and validate the photometric results in Supplementary Discussion. In this work, we limit our analyses to the images observed with the Advanced Himawari Imager (AHI) onboard the Himawari-8 satellite, which is the first "third-generation meteorological satellite" characterized by access to as many as 16 bands 10 . There also are three other third-generation satellites equipped with 16-band imagers: GOES-16 and GOES-17 44 , which are operated by the National Aeronautics and Space Administration (NASA) and the National Oceanic and Atmospheric Administration (NOAA), and GEO-KOMPSAT-2A 45 , which is operated by the Korea Meteorological Administration (KMA). Although they have been operational since November 2016, March 2018, and July 2019, respectively, flux values outside the Earth's disk are found to be masked and unavailable in their published standard L1B products, which makes these satellites unusable for our analyses. Combined analyses with those satellites would improve the precision of our results, if their images of the space around the Earth's disk were to be made publically available.
Reduction of the AHI images is performed by the Japan Meteorological Agency (JMA) 41,46 . This involves geometric re-gridding into a fixed coordinate system using a Lanczos-like resampling; a linear and global flux calibration using deep space as the background, with a solar diffuser for the optical and a blackbody for the infrared as standard flux sources; and nonlinear flux calibrations based on pre-launch ground tests. No corrections are made for atmospheric effects, nor is masking for cosmic rays performed.

Photometry of the Himawari-8's image
About 6-years of 16-band images have been obtained by the Himawari-8 satellite at a frequency of once every ten minutes. We restricted the scope of our analyses to data observed from January 2017 to June 2021 because systematic biases in the flux calibration and in the positioning were identified in observations from July 2015 to November 2016 (see Supplementary Discussion). From this large number of images, we tried to find and measure the light from Betelgeuse, which is moving with the speed of the Earth's rotation, 15 deg hour −1 , and which is found to be located in the field of view about once every two days. First, we calculated the times when Betelgeuse is in the direction of the Earth using the ASTROPY package 47 , and we collected all the images in which Betelgeuse could be detected. Then we searched for the location of the Betelgeuse signal in each of these images in Bands 1-10, where the signal-to-noise ratio (SNR) is high enough to detect the star. For each of those high-SNR bands, we searched for the location of the star by choosing the pixel with the largest counts among those pixels where Betelgeuse is expected to be located. In contrast, we inferred the locations in the other bands (Bands 11-16) by shifting from the position of Betelgeuse in Band 7, where the SNR is the highest, by around 30 arcsec depending upon the differences in the times of observation for the two bands.
The pixel scale of the published images in each band is larger than the actual radius of the point-spread function (PSF) estimated by the radius of the telescope. Thus, the light from a point source occupies only 1×1 pix 2 of the un-reduced images, and it spreads out to ∼2 × 2 pix 2 during the reduction process (see the left panels of Fig. 1). Due to this small PSF radius, PSF photometry does not work well, so instead we used the aperture-photometry technique with a 7 × 7 pix 2 square aperture. The background flux in the aperture was estimated and subtracted by fitting the flux within a square annulus around the star using a 2D polynomial function. We set the inner width of the annulus to be 7 pix and the outer width to be 43, 21, and 21 pix for the bands with the spatial resolutions of 0.5, 1.0, and 2.0 km, respectively. We selected those widths in order to minimize the impact of stripe-shaped fluctuations of the background. We selected the order of the polynomial to be in the range 0-4 for each measurement so that the Akaike Information Criterion calculated from the fitting residuals is minimized. Then, we estimated the stellar flux by summing the background-subtracted spectral irradiances within the 7 × 7 pix 2 aperture.
Since the background shapes are not guaranteed to be continuously connected between adjacent swaths, we excluded from our analyses cases in which the 7 × 7 pix 2 aperture spans multiple swaths. Furthermore, we trimmed the square annulus used for background estimation by requiring that the annulus does not span multiple swaths. Moreover, in order to avoid contamination from terrestrial lights, we introduced distance thresholds from the edge of the Earth's disk (see Extended Data Fig. 1).

Validation of the photometry
Here, we performed a wide variety of validation steps, not only with Betelgeuse but also with other bright stars located in the field of view of the Himawari-8 satellite (Extended Data Fig. 2). From the validation described in Supplementary Discussion and Supplementary Figs. 1-6, we found no evidence of any serious systematic bias in our photometry, and we thus conclude that the measurements are ready for further analyses. Nevertheless, it would be important, in a coming decade, to examine whether or not there are overlooked systematic effects in our method.

SED fitting to the Himawari-8 light curves
The right panel of Fig. 1 shows the spectral energy distribution (SED) averaged in time together with a model SED having the typical atmospheric parameters of Betelgeuse. The photometric bands with wavelengths 8 µm (Bands 1-10) are almost unaffected by the circumstellar dust emission, and thus they can be used to determine the atmospheric parameters of Betelgeuse. In contrast, the other bands receive radiation from silicate-dust emission at around 10 µm, and they can thus be used to measure the amount of dust formed around Betelgeuse. We therefore first fitted the SED for Bands 1-10 with a model photosphere for Betelgeuse to determine the photospheric parameters-R, T eff , and A(V )-and we then determined τ 10 using Bands 12-14, which have the highest SNR for the dust emission. Five examples of the fitting results are shown in Extended Data Fig. 3 Our final time-series photometric catalog of Betelgeuse is an incomplete one; the catalog lacks fluxes in several bands on several observations. This is especially serious for longer-wavelength bands, where the sky-background emission is stronger, and thus a larger number of photometric measurements was removed from the catalog when the star was located near the edge of the Earth's disk. Moreover, the signal-to-noise ratios in the mid-infrared bands are typically small, ∼1. In order to deal with these difficulties, we binned each light curve over 5-12 observations requiring that all the binned data points for Bands 1-10 include five or more points. Note that the number of days in each bin is automatically determined by the observation frequency, which depends on the observation epoch. For example, there is a clump of observations in May-June 2018 with a frequency of 0.85 days, which leads to the smaller number of days in the bins. Finally, we calculated the resulting median and its standard error for each bin. Model SED for the photosphere We fitted each point of this "complete" SED curve for Bands 1-10 with a photospheric model SED using the LMFIT package 48 . To obtain the model spectrum for the photosphere, we interpolated a grid of BT-NextGen model spectra 49 with solar metallicity and the surface gravity log g = 0.0, assuming that this one-dimensional model can be used to model the spectra of Betelgeuse at any phase. We used the Cardelli's reddening law 50 with the selective-to-total extinction ratio R V = 3.1 51 , assuming that this law can be extrapolated to longer wavelengths. We also used the parallax of 6.55 mas measured with the HIPPARCOS satellite 52 . Although the accuracy of the parallax of Betelgeuse is controversial 53,54 , the choice of parallax value only changes the fitted radius in a relative way and does not affect the interpretation of the results. For the same reason, we do not take into account the error in the parallax value. The model spectra are then convolved with the response function of the AHI (September 2013 version retrieved from https://www.data.jma.go.jp/ mscweb/en/himawari89/space_segment/spsg_ahi.html) to calculate the model SEDs for Betelgeuse, and we fitted the observed SED with the model.

Dependence on the model selection
For the grid of model spectra, we first tried the MARCS grid 55 with 15M 54,56,57,58,59 , which is usually used to model the spectra of red supergiants 51,60 . We found that the MARCS grid gives a smaller Akaike Information Criterion than the BT-NextGen grid, probably because the wavelengths and the strengths of the molecular bands in the MARCS grid are more accurate. Moreover, considering the fact that several previous works 51 succeeded in reproducing observed optical spectra of red supergiants dominated by molecular TiO bands with MARCS models, MARCS models appear to be an appropriate grid for the synthesis of TiO bands in red supergiants. However, we found that the results of the fits obtained using the MARCS grid were physically unreasonable for most cases that we considered, e.g., A(V ) < 0 mag, although the reason, possibly related to 3D effects and/or circumstellar gas, is unclear. Thus, instead of the MARCS grid, we used the BT-NextGen grid to compute the model SED. Testing both the MARCS and Next-Gen grids, we found that the resulting model parameters we obtained using the BT-NextGen grid vary with time in the same way as those obtained using the MARCS grid; i.e., the differences between the time variations obtained using both grids are roughly constant (see the left panels of Extended Data Fig. 4). Therefore, the results of this paper using the Next-Gen grid are not affected quantitatively by the choice of model grid, as long as the time variations are interpreted in a relative way.
We have also tested whether or not other model assumptions affect the results. For example, a change in log g from 0.0 to −0.5 makes T eff systematically ∼80 K warmer and A(V ) systematically ∼0.3 mag larger. Similarly, changing R V from 3.1 to 4.4 51 decreases T eff by 10 K and A(V ) by 0.05 mag. From these tests, we concluded that although the absolute values of R, T eff , and A(V ) depend on the model assumptions, their relative variations appear to have been determined to the accuracy ∼10 K in T eff and ∼0.05 dex in A(V ), both of which are much smaller than the fitting errors, ∼40 K and ∼0.2 dex, respectively.

Three-dimensional (3D) effects in the photosphere
The temperatures that previous works 8, 19 61,62 and we have both derived assumed a one-dimensional (1D) geometry for the model atmosphere. Thus these temperatures are affected by three-dimensional (3D) effects in the optical spectrum of the real star 7 . In fact, the SED of a model with a given T eff , especially the strength of the TiO molecular absorption bands that appear in the optical spectra of red supergiants 8 , has a variety of shapes due to 3D effects 35 . Therefore, when we fit the observed SED of the real 3D photosphere of Betelgeuse with a 1D model SED, the fitted values of T eff and A(V ) depend not only on T eff and A(V ) themselves but also on the 3D effects 35 . In particular, a time variation of the patchiness of the surface 63,64 could affect the time variations of T eff and A(V ) 7 .
As a simplified test of the 3D effect on our fitting result, inspired by the optical image of Betelgeuse obtained during the Dimming 1,25 , we considered a model photosphere composed of two regions with different T eff , which is similar to the model used by a previous work 7 . In this model, we assume that the cooler part of the photosphere has a varying temperature and a varying area and that the warmer part has a fixed temperature of 3650 K. In contrast to the previous model 7 , we are able to consider two additional free parameters, A(V ) and R; without these two parameters we could not obtain a satisfactory fit to our multi-wavelength data. Altogether, we considered four free parameters: R, A(V ), and the temperature and area fraction of the cooler part; e.g., when the entire area of the photosphere is occupied by the cooler part, the fraction is 1.0. With this four-parameter model, due to the strong degeneracy between A(V ) and the area, sometimes the fitting failed and gave unphysical results. To remove these cases, we excluded fitted results with areas greater than 1.0 or lower than 0.0 as well as results for which the temperatures reached the boundaries of the acceptable range (we assumed 2800 and 4000 K for the lower and upper boundaries, respectively).
The right panels of Extended Data Fig. 4 compare the fitted results obtained with the three-parameter model adopted in the main text and the fourparameter model described here. We found that with the four-parameter model the contribution of the increased A(V ) to the Dimming decreases but that it is still larger than zero (∼0.4 mag). This result supports the robustness of our results in a qualitative sense. Since the four-parameter model examined here is too simple to take 3D effects fully into account 7 , future theoretical modeling with self-consistent hydrodynamic models 65 will be necessary.

Comparison to previous radial-velocity (RV) measurements
To check our SED fitting results, we compared the radius determined herein to previous radial-velocity (RV) measurements. Because the radius is given by the integration of the RV, we can validate the radius curve using literature RV curves.
First, we found time lags of up to ∼80 days between the times of maximum brightness and minimum radius (Fig. 2). Although such a large lag was not expected in a simple 1D hydrodynamic model of Betelgeuse 54 , a lag between light and RV curves has been observationally found for Betelgeuse 9,19 66 . Moreover, we found clockwise trajectories between the R/T eff and Band-3 flux for most observation periods (Extended Data Fig. 3h-o), which has been predicted theoretically 54,67 . The exception, the relation between R and Band-3 flux between 27 October 2017 and 24 December 2018 shown in panel i, could be attributed to effects other than the fundamentalmode pulsation, such as the long-secondary period 68 ; however, the actual cause remains unclear.
Second, we estimated the RV curve by calculating the differential of our radius curve, and compared it with literature RV curves 19 66 . This comparison reveals an overall good agreement between our data and the results cited above, supporting the validity of our radius measurement, and hence other parameters. However, minor differences of up to ∼5 km s −1 may exist, particularly in September 2018-November 2018 and September 2019-December 2019. Plausible causes of these differences are: (1) our imperfect modeling of Himawari-8's SEDs with a 1D model atmosphere, as discussed in the previous subsection, and (2) differences in the atmospheric layers that were observed in the present and previous studies, and hence the real, physical difference in RVs 19 .

Model SED for the dust emission
After fitting the SED curve for Bands 1-10, we calculated the expected photospheric flux curves for Bands 11-16 and subtracted them from the observed flux curves. These residual fluxes are likely to come from the circumstellar dust emission, and they are useful for estimating the variation of the amount of circumstellar dust.
For the model spectra of the circumstellar dust emission, we used the radiative-transfer code DUSTY V2 69 . For the settings of the DUSTY calculation, we basically followed the procedure described in a previous work 70 , except for several modifications. First, we only fitted τ 10 and fixed the other parameters. We used the "warm silicates" dust composition 71 , a constant grain size of 0.3 µm, and a constant dust temperature of 1300 K at the inner boundary, which roughly corresponds to the sublimation temperature of silicate dust 28 . Although inner-boundary temperatures cooler than 1300 K are sometimes found by observations of red supergiants 70 , we adopted this value since a different inner-boundary temperature only results in a systematic bias in τ 10 as long as the temperature is independent of time. We also tried the assumption of a constant radius at the inner boundary, and we found that the results do not change in a qualitative sense. For the stellar radiation that illuminates the dust, we used the time variation of T eff determined in the previous section to calculate the natal photospheric SED, after applying spline smoothing to reduce the noise in each measurement. Then, we calculated the expected silicate-dust emission in Bands 12-14 on the τ 10 grid and determined the τ 10 value that best matches the observed fluxresidual curve for each of Bands 12-14. Finally, we calculated the mean of the three τ 10 values determined from Bands 12-14 to obtain the final estimate of τ 10 . We note that τ 10 that we determined is in fact not the actual optical depth at 10 µm; it is rather a model parameter that we use to match the observed and model fluxes of Bands 12-14, and we use it as an indicator of the amount of circumstellar silicate around Betelgeuse.

Code Availability
We used the LMFIT package 48 to fit the observed SED in the optical to near infrared with publicly available NextGen 49 and MARCS 55 model grids. We also used the DUSTY V2 code 69 to calculate the dust emission. We have not made the codes for the photometry publicly available because they are not prepared for open use. Instead, we have made the photometric catalog available, as described in the Data Availability statement.

Data Availability
Reduced images taken by the Himawari-8 satellite are publicly available at the Data Integration and Analysis System via https:// diasjp.net/en/. The final photometric product of this paper is available at https://d-taniguchi-astro.github.io/homepage/ Data_Himawari_en.html, which is archived as Supplementary Data in Supplementary Information. Fig. 1    In all the panels, the black circles indicate the results adopted in this paper (the same as Fig. 2), which use the grid of Next-Gen model spectra 49 , assuming three free parameters in the fitting. a-c, The effect of the adopted modelspectrum grid: Next-Gen 49 (blue circles) vs. MARCS 55 (upward-pointing orange triangles). d-g, The effect of the number of parameters: three parameters (blue circles) vs. four parameters (downward-pointing pink triangles).

SUPPLEMENTARY INFORMATION
In this Supplementary Information pdf, we describe our validation steps of the photometry.

Signal-to-noise ratio (SNR) of the imager
According to the results of the pre-launch verification carried out by JMA 75 , the noise level of AHI is around 0.7, 0.2, and 0.02 nW m −2 µm −1 per pixel for visible, near-infrared, and mid-infrared bands, respectively. Assuming independence of noise among pixels, the SNR of our flux measurements for Betelgeuse due to noise would be approximately 2-20 and 0.1-1 for Bands 1-10 and 11-16, respectively, before applying the temporal binning described later. Those estimated SNRs are roughly consistent with the sample-to-sample variation of our photometric results (Extended Data  Fig. 3a).

Calibration accuracy
Long-term in-orbit monitoring of the AHI (https://www. data.jma.go.jp/mscweb/en/oper/eventlog/Update_of_ Calibration_Information_2020.pdf) indicates that the sensor gain in the visible and near-infrared bands (Bands 1-6) is affected by seasonal fluctuations of up to 0.5% and by a quasi-linear fading trend of up to 0.5% year −1 . The estimated stellar parameters are affected by the uncertainties due to these calibration biases. However, the impact of the biases is likely to be negligible, since these are sufficiently smaller than the random ∼5% errors in the observed binned fluxes of Betelgeuse in the visible bands (Extended Data Fig. 3a). Moreover, the weakening trend results only in a weak, long-term monotonic trend in the stellar parameters of Betelgeuse, which does not affect the interpretation of the results.

Saturation
Since the AHI has been designed to observe the Earth, it is necessary to confirm that stellar observations are not affected by the saturation of the detector. We confirmed that the maximum valid fluxes of the AHI in all bands (Extended Data Fig. 1) are sufficiently larger than the observed fluxes of Betelgeuse shown in the right panel of Fig. 1. Therefore, the results of this study are not affected the saturation.

Quantization noise
In the opposite direction to saturation, the photometry of AHI data might be affected by quantization during the reduction. Due to the limited size of data that can be downlinked in real time, satellite data is compressed to 11 to 14 bits 76 . Thus, the quantization noise caused by the relatively large least significant bit (LSB) of each band (Extended Data Fig. 1) would have an impact on the photometric error. However, we confirmed that the quantization noise is at least two times smaller than the photometric error due to background fluctuation in all bands. Moreover, with the exception of Band 16, which we did not use for the analysis of Betelgeuse, the typical flux of Betelgeuse is at least five times larger than the quantization noise. This consideration demonstrate that the quantization noise has a non-negligible, but non-dominant contribution to the photometric error.

Point-spread function (PSF)
The PSFs of the reduced Himawari-8 images were estimated by stacking the flux of Betelgeuse around the brightest pixel. The results in Supplementary Fig. 1 show that the PSF occupies ∼2 × 2 pix 2 , regardless of the wavelength or the pixel scale of the observation band; i.e., 0.5 km for Band 3, 1.0 km for Bands 1, 2, and 4, and 2.0 km for the others. This finding supports the conclusion that the PSF mainly results from the resampling process in the reduction. Because of the method used for resampling 76 , we can safely assume that the signal from a point source like a star does not spread further than ±2 pix away from the brightest pixel.

Positioning accuracy
Since the location of Betelgeuse in the images of the mid-infrared bands was inferred from that of Band 7, large errors in the inferred position, if present, would cause flux loss. To evaluate the accuracy of the inter-band positioning inference, we estimated the position of Betelgeuse in the visible and near-infrared bands from Band 7 in the same manner as for the mid-infrared bands, and we compared the inferred and actual positions ( Supplementary Fig. 2). Focusing on the period after January 2017, both in the east-west (x) and north-south (y) directions, we found that the positioning-inference errors are typically less than 5 arcsec, which corresponds to 0.5 pix in the mid-infrared bands. Thus, our aperture size, 7 × 7 pix 2 , is large enough to cover all of the size of the circumstellar envelope of Betelgeuse (< 0.5 pix 22 ), the PSF (±2 pix at most), and the positioning error (±0.5 pix at most) after 2017. Therefore, it is unlikely that the results in the mid-infrared bands missed flux from Betelgeuse due to PSF and/or positioning errors after January 2017. In 2015 and 2016, however, the positioning errors were considerably larger in some cases, especially in Band 5 (the middle panel of Supplementary Fig. 2). Thus, for safety, we excluded all measurements in 2015 and 2016 from our analyses.

Light curves of bright stars other than Betelgeuse
In order to validate the photometric results, we applied the photometric method we used for Betelgeuse to four bright stars within the field of view of the AHI (Extended Data Fig. 2). Extended data Fig. 3a and Supplementary Figs. 3, 4 show the resulting 16-band light curves for Betelgeuse, two non-variable stars (Rigel and Procyon), and two red giant stars (δ Oph and α Cet), using the observational data from the AHI in its period of operation from July 2015 to June 2021. The light curves of the two non-variable stars show that the fluxes in Bands 1-4 before the beginning of 2016 were systematically enhanced relative to the period after this. Similarly, the light curves of the two red giant stars also show enhancement in the early epoch. Though the cause of this enhancement, which could be related to an update in the reduction process by JMA in 9 May 2016 77 , is unclear, in order to avoid this enhancement, we decided to use the data only after 1 January 2017 for the subsequent analyses.
In order to test quantitatively whether there is time variability in the light curves of the two non-variable stars after January 2017, we tested the null hypothesis that there is no autocorrelation in the light curves with time delays up to 30 days using the Ljung-Box test 78 implemented in the STATSMOD-ELS package 79 . We found that, among the 16 bands for each of the 2 stars, only one band of one star has a p value smaller than 0.05 (Extended Data Fig. 1); i.e., there is no evidence of any detectable autocorrelation in the light curves of the non-variable stars. Considering that these two stars are sufficiently bright in the optical and near-infrared bands (Bands 1-5), this result indicates that our photometry in these bands is not likely to be affected by time-dependent systematic biases. Further, it indicates that photometry in the mid-infrared bands is always likely to have accurate zero points, and hence partially supports the reliability of the mid-infrared photometry. Nevertheless, it would be better to test the mid-infrared photometry with infraredbright, non-variable stars; unfortunately, there are no such stars in the field of view of Himawari-8.

Atmospheric contamination
The radiance of a star measured while it is located close to the edge of the Earth's disk might be affected by emission and Rayleigh scattering from the terrestrial atmosphere, and thus these effects need to be examined. First, in order to evaluate the amount of spurious flux resulting from such contaminating emissions, we applied the photometric procedure used for Betelgeuse to empty points around the Earth, rather than to Betelgeuse. The measured spurious fluxes of the empty points correspond to the amount of bias in the photometry due to unsatisfactory background subtraction. As shown in the left panels of Supplementary Fig. 5, the amount of background bias actually depends on the distance from the edge of the Earth's disk, and the dependence of this atmospheric emission differs among the bands. We confirmed that the measurements contaminated by this bias had been removed from our dataset after setting the thresholds for the Earth-Betelgeuse distances for the observed images in the individual bands (Extended Data Fig. 1).
Next, we evaluated the effects of the terrestrial sky on starlight using the photometric data for the two non-variable stars, Rigel and Procyon, measured with the same distance thresholds. We found that the measured fluxes of the two stars have no significant dependence on the Earth-star distances (the right panels of Supplementary Fig. 5). This fact indicates that our photometric measurements also are not likely to be affected by Rayleigh scattering or by other effects relating to the terrestrial atmosphere.

Contamination by Solar-System objects
Solar-System objects like the Moon and Venus sometimes appear in images from meteorological satellites 10 . If these objects are located too close to the stars whose fluxes are to be measured (Fig. 2), they could affect photometry. Using the ASTROPY package 80 , we calculated the direction of bright Solar-System objects as seen from Himawari-8 when the stars are located within the field of view. Herein, we only considered bright objects that could be detected with Himawari-8: the Moon, Mercury, Venus, Mars, Jupiter, and Saturn. As a result, we confirmed that none of these objects approach too close to the stars, and thus their resulting contamination does not affect the photometry. We note that no other bright stars are present around each of the target stars.
Another source of contamination is scattered light in the optical and near-infrared wavelengths from the Sun. This effect becomes strongest at around midnight when the satellite observes the direction close to the Sun, particularly around the equinoxes 81 . However, due to a close relationship between the observing month and hour for each star, midnights correspond to specific month(s), e.g., January and February for Betelgeuse. Similarly, we confirmed that none of our target stars are observed at around midnight during equinox periods. Moreover, we confirmed that scattered light was significantly reduced by the background-subtraction procedure. Although the background flux is as much as, e.g., 1.0 and 0.6 nW m −2 µm −1 per pixel for Bands 1 and 6, respectively, before the subtraction, they are reduced to 0.02 and 0.003 nW m −2 µm −1 , respectively.

Comparison with literature photometric data of Betelgeuse and Rigel
From the validation steps so far, our light curves appear to be independent of both the date of observation and the distance from the edge of the Earth's disk; i.e., they seem to be accurate in a relative sense. Finally, in order to evaluate the absolute accuracy of the photometry, we compared the Himawari-8's photometric data for Betelgeuse and Rigel with literature measurements.
First, the top panel of Supplementary Fig. 6 compares the optical light curves obtained with Himawari-8's Bands 2 and 3 (centered at 0.51 and 0.64 µm, respectively) to the V-band light curve (0.55 µm) cataloged in several sources: the American Association of Variable Star Observers (AAVSO) database, a work doing daylight observation 82 , and a work using images obtained by the STEREO Mission 83 . We found that our light curves are well parallel to that of the V band. Moreover, as expected from the central wavelengths, the V-band light curve is located between those of Bands 2 and 3 and is relatively close to Band 2.
Second, the bottom panels of Supplementary Fig. 6 compare Himawari-8 and literature SEDs of Betelgeuse and Rigel: the infrared spectra of Betelgeuse and Rigel taken with the Short-Wavelength Spectrometer (SWS) on the Infrared Space Observatory (ISO) satellite 84 (TDT numbers 69201980 and 83301505, respectively); the mean spectral irradiances of Betelgeuse between January 2017 and June 2021 in the six visible and near-infrared wavebands (B, V, R, I, J, and H) taken from the AAVSO database; and the spectral irradiances of Rigel in the same six bands 85 . The figures show overall good agreement between the Himawari-8 SEDs and those from the literature, although there is a discrepancy between the two in the optical and near-infrared SED, especially for Betelgeuse in the near infrared. The reason for this discrepancy is unclear, but differences in the dates of observation and in the bandwidths may contribute to the discrepancy, at least in part.
Considering these results, the absolute accuracy of the photometry is estimated conservatively as ∼20%, although the absolute accuracy of the photometry does not affect the final results of this paper; i.e., the time variation of the brightness and stellar parameters of Betelgeuse.    85 (black diamonds in the right panel), and infrared spectra observed with the ISO satellite 84 (gray lines). We ignored the photometric data for Betelgeuse during the Great Dimming period, August 2019-July 2020 indicated by vertical red dashed lines in a, in the averaging process. We calculated the spectral irradiances for these photometric data from the literature by converting from Vega to AB magnitudes using the zero points given in Table 1 in ref. 86 . We retrieved the AAVSO's photometric data on 16 February 2022.