# Long Term Evolution of Surface Features on the Red Supergiant AZ Cyg

Ryan P. Norris, Fabien R. Baron, John D. Monnier, Claudia Paladini, Matthew D. Anderson, Arturo O. Martinez, Gail H. Schaefer, Xiao Che, Andrea Chiavassa, Michael S. Connelley, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Ryan P. Norris, Fabien R. Baron, John D. Monnier, Claudia Paladini, Matthew D. Anderson, et al. "Long Term Evolution of Surface Features on the Red Supergiant AZ Cyg." The Astrophysical Journal Letters, In press, 10.3847/1538-4357/ac0c7e, hal-03367843

HAL Id: hal-03367843

https://hal.science/hal-03367843

Submitted on 6 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_3.png" /></figure>

CrossMark

# Long Term Evolution of Surface Features on the Red Supergiant AZ Cyg

Ryan P. Norris<fnref n="1" />, Fabien R. Baron<fnref n="2" />, John D. Monnier<fnref n="3" />, Claudia Paladini<fnref n="4" />, Matthew D. Anderson<fnref n="2" />,<fnref n="5" />, Arturo O. Martinez<fnref n="2" />, Gail H. Schaefer<fnref n="5" />, Xiao Che<fnref n="6" />, Andrea Chiavassa<fnref n="7" />, Michael S. Connelley<fnref n="8" />, Christopher D. Farrington<fnref n="5" />, Douglas R. Gies<fnref n="2" />, László L. Kiss<fnref n="9" />,¹⁰,<fnref n="11" />, John B. Lester<fnref n="12" />,<fnref n="13" />, Miguel Montargès<fnref n="14" />, Hilding R. Neilson<fnref n="12" />, Olli Majoinen<fnref n="5" />, Ettore Pedretti<fnref n="15" />, Stephen T. Ridgway¹⁶, Rachael M. Roettenbacher<fnref n="17" />, Nicholas J. Scott<fnref n="18" />, Judit Sturmann<fnref n="5" />, Laszlo Sturmann<fnref n="5" />, Nathalie Thureau<fnref n="19" />, Norman Vargas<fnref n="5" />, and Theo A. ten Brummelaar<fnref n="5" />

Received 2021 January 7; revised 2021 June 15; accepted 2021 June 16; published 2021 October 1

## Abstract

We present H-band interferometric observations of the red supergiant (RSG) AZ Cyg that were made with the Michigan Infra-Red Combiner (MIRC) at the six-telescope Center for High Angular Resolution Astronomy (CHARA) Array. The observations span 5 yr (2011–2016), which offers insight into the short and long-term evolution of surface features on RSGs. Using a spectrum of AZ Cyg obtained with SpeX on the NASA InfraRed Telescope Facility (IRTF) and synthetic spectra calculated from spherical MARCS, spherical PHOENIX, and SAtlas model atmospheres, we derive $T_{\text{eff}}$ is between 3972 K and 4000 K and log g between −0.50 and 0.00, depending on the stellar model used. Using fits to the squared visibility and GAIA parallaxes, we measure its average radius $R = 911^{+57}_{-50} R_\odot$. Reconstructions of the stellar surface using our model-independent imaging codes SQUEEZE and OITOOLS.jl show a complex surface with small bright features that appear to vary on a timescale of less than one year and larger features that persist for more than one year. The 1D power spectra of these images suggest a characteristic size of 0.52–0.69 $R_\star$ for the larger, long lived features. This is close to the values of 0.51–0.53 $R_\star$ that are derived from 3D RHD models of stellar surfaces. We conclude that interferometric imaging of this star is in line with predictions of 3D RHD models but that short-term imaging is needed to more stringently test predictions of convection in RSGs.

*Unified Astronomy Thesaurus concepts: Optical interferometry (1168); Red supergiant stars (1375); Stellar surfaces (1632); Stellar properties (1624); Evolved stars (481)*

## 1. Introduction

As stars of mass $9 M_\odot–25 M_\odot$ transition to core He burning, they evolve toward the red supergiant (RSG) stage (Meynet et al. 2015). Reaching radii up to $1500 R_\odot$, stars in this evolutionary stage are cool (3400–4100 K) and have luminosities of 20,000–300,000 $L_\odot$ (Levesque et al. 2005). RSGs display high mass-loss rates ($10^{-7}–10^{-4} M_\odot \text{ yr}^{-1}$), leading to complex circumstellar environments and a variety of evolutionary outcomes (van Loon et al. 2005; De Beck et al. 2010; Mauron & Josselin 2011). Stars may remain in this stage until

they end in supernovae or they may evolve blueward, possibly to return once more to the RSG stage (Meynet et al. 2015). Understanding the evolution of RSGs has implications for a variety of astrophysical questions, including the chemical evolution of the universe and the diversity of supernovae observed.

RSGs display semi-regular variations that are comprised of a short period of hundreds of days and a longer period lasting thousands of days (Kiss et al. 2006). Although pulsations have been linked to some of these variations, it is likely that the motions of large surface features also play a role. In the early 1970s, Stothers (1972) and Schwarzschild (1975) suggested that the observed periodicity was reminiscent of the variability displayed by convection in Sun-like stars but on a larger scale and longer timescale. Using standard mixing length theory, Stothers (1972) suggested that RSGs possess a small number of

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Department of Physics, New Mexico Institute of Mining and Technology, 801 Leroy Place, Socorro, NM 87801, USA; ryan.norris@nmt.edu</footnote>

<footnote n="2">²Center for High Angular Resolution Astronomy, Department of Physics and Astronomy, Georgia State University, P.O. Box 5060, Atlanta, GA 30302-5060, USA</footnote>

<footnote n="3">³Department of Astronomy, University of Michigan, 1085 S. University Avenue, Ann Arbor, MI 48109-1090, USA</footnote>

<footnote n="4">⁴European Southern Observatory, Santiago, Chile</footnote>

<footnote n="5">⁵The CHARA Array of Georgia State University, Mount Wilson Observatory, Mount Wilson, CA 91023, USA</footnote>

<footnote n="6">⁶Independent Researcher</footnote>

<footnote n="7">⁷Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Lagrange, CS F-34229, Nice, France</footnote>

<footnote n="8">⁸University of Hawaii, Institute for Astronomy, 640 N. Aohoku Place, Hilo, HI 96720, USA</footnote>

<footnote n="9">⁹Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, 1121 Budapest, Konkoly Thege Miklós út 15–17, Hungary
¹⁰Institute of Astronomy, School of Physics, University of Sydney, NSW 2006, Australia</footnote>

<footnote n="11">¹¹ELTE Eötvös Loránd University, Institute of Physics, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary</footnote>

<footnote n="12">¹²David A. Dunlap Department of Astronomy & Astrophysics, University of Toronto, 50 Saint George Street, Toronto, ON, M5S 3H4, Canada</footnote>

<footnote n="13">¹³Department of Chemical & Physical Sciences, University of Toronto Mississauga, Mississauga, ON L5L 1C6, Canada</footnote>

<footnote n="14">¹⁴Institute of Astronomy, KU Leuven, Celestijnenlaan 200D Box 2401, B-3001 Leuven, Belgium</footnote>

<footnote n="15">¹⁵ LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, 92195 Meudon, France
¹⁶NSF’s NOIRLab, P.O. Box 26732, Tucson, AZ 85726-6732, USA</footnote>

<footnote n="17">¹⁷Yale Center for Astronomy and Astrophysics, Yale University, 46 Hillhouse Avenue, New Haven, CT 06511, USA</footnote>

<footnote n="18">¹⁸Space Science & Astrobiology Division, NASA Ames Research Center, Moffett Field, CA 94035, USA</footnote>

<footnote n="19">¹⁹Farr Institute Scotland, Nine Edinburgh Bioquarter, Little France, UK</footnote>
extremely large convection cells (~0.5–1R<sub>⋆</sub>) with lifetimes of 1000 s of days. Meanwhile, Schwarzschild (1975) extrapolated from the behavior of supergranulation in the Sun, predicting a surface occupied by dozens of large granules with lifetimes of 100 s of days. Indeed, both types of features are observed in 3D radiative hydrodynamics (RHD) models of RSGs (Chiavassa et al. 2011a, 2011b). Kravchenko et al. (2019) applied tomography to high resolution spectra of the RSG μ Cep finding that line-of-sight velocity variations exhibited a phase-shift behind temperature and photometric variations, a feature also found in 3D RHD simulations. The authors suggested that this was evidence that the shorter period of hundreds of days was related to convection within the star. The presence of a granulation-related signal in the irregular variability of Betelgeuse was noted by Ren & Jiang (2020) who determined a granulation time of $138_{-2}^{+8}$ days by fitting the posterior of the power density spectrum (PSD) of American Association of Variable Star Observers (AAVSO) observations of the star to a modification of an equation Harvey (1985) first used to model solar granulation. However, convection is likely not the cause of the dominant short-term periodicity in RSGs. For example, Joyce et al. (2020) found that radial fundamental and first overtone pulsations, originating in the κ-mechanism, are the cause of Betelgeuse's short-term (<400 day) variability.

In 2019, the dramatic dimming of Betelgeuse (Guinan & Wasatonic 2020) brought further attention to the variability of RSGs and its myriad causes. Dharmawardena et al. (2020) suggested the dimming was caused by the presence of a large cool spot in the photosphere of the star based on decreased flux at submillimeter wavelengths. Meanwhile, other authors have suggested that spectrophotometry (Levesque & Massey 2020), polarimetry (Cotton et al. 2020), and interferometry (Montargès et al. 2021) provide evidence that the decrease in visual magnitude was the result of obscuring of the surface, likely due to large-grain dust from mass loss. This claim was bolstered by Dupree et al. (2020), who reported the presence of a large bright spot on Betelgeuse prior to its dimming event. Because the location of the spot was similar to the location of the dimmed portion of the star, they concluded that pulsation and convection elevated material outward and that this material later cooled, dimming the star. Indeed, the contribution of convection to mass-loss has been a long standing question. Using tomography, Josselin & Plez (2007) found velocity gradients indicative of vertical motion on timescales that match the lifetime of convection cells predicted by Schwarzschild (1975). Turbulent pressure from such motions would result in lifting that could play a role in mass loss. However, Arroyo-Torres et al. (2015) studied the extended atmosphere of three RSGs using spectro-interferometry and found that hydrostatic model atmospheres, 3D RHD simulations of convection, and 1D pulsation models each produced atmospheres that were more compact than observed. While pointing out the limitations of current 3D convection and 1D pulsation models, they noted that these models do not take into account short-timescale velocity gradients, which may result from granulation-like motion in convection cells.

Because these stars are quite large and because their convection related features are predicted to be large relative to the stellar radius, studying convection in RSGs is amiable to imaging studies. Using aperture masking at the William Herschel Telescope, Buscher et al. (1990) and Tuthill et al. (1997) found evidence of hotspots via detection of asymmetries

in visible wavelength observations of three large RSGs. Because the timescale of the variations matched predictions by Schwarzschild (1975), Tuthill et al. (1997) suggested that the hotspots were related to convection. Later observations of Betelgeuse by Young et al. (2000) showed that detection of these hotspots is wavelength dependent and that at some wavelengths, RSGs appear to be featureless. Observations by Haubois et al. (2009) of Betelgeuse in the H-band, where $H^−$ continuum opacities are at a minimum (1.6 µm), suggested that large features do exist in the photosphere, as predicted by 3D RHD models (Chiavassa et al. 2010; Montargès et al. 2014). Using the Center for High Angular Resolution Astronomy (CHARA) Array, Baron et al. (2014) imaged the surface of T Per and RS Per in the H-band, finding more evidence that large surface features are ubiquitous. Looking at short-timescale evolution, Montargès et al. (2018) monitored the evolving surface features of CE Tau between 2016 November and December, finding evidence of short-term changes. As for longer timescale evolution, López Ariste et al. (2018) collected spectro-polarimetry of Betelgeuse between 2013 and 2018, and found evidence of warm regions of scale 0.6 R<sub>⋆</sub> in the photosphere that persisted for up to 4 yr that the authors tied to giant convection cells. Climent et al. (2020) presented reconstructions the surface of V602 Carinae obtained with VLTI-PIONIER. These data, collected in 2016 and 2019, are limited by temporal and angular resolution in their ability to study the size of individual granules. However, they suggest that the surface does not change significantly over a period of at least 70 days and that 3D RHD models are able to reproduce surface features as observed in images of RSGs. Moreover, a comparison of continuum images to 3D RHD models suggests that the large features present in these images are indeed due to convection.

Our work in many ways complements this study of V602 Carinae. To better grasp the timescale and spatial scale of surface variations on RSGs, we undertook a 5 year study of the RSG AZ Cyg using the Michigan Infrared Combiner (MIRC) on the CHARA Array. In Section 2, we discuss the observations and data reduction of the interferometric and spectroscopic data we used in this study. In Section 3, we present the fundamental stellar parameters of AZ Cyg. In Section 4, we discuss our image reconstruction techniques and present images of AZ Cyg from the 5 yr of observations. In Section 5, we use 1D power spectra of the images to compare the reconstructions to predictions of surface feature size by 2D and 3D models. We conclude the paper in Section 6 with a summary of the work and discussion about the need for observations at shorter timescales.

## 2. Observations

### 2.1. AZ Cyg

AZ Cyg is an M2-4.5 Iab star, with a parallax distance of $d = 2090_{-115}^{+130}$ pc based on Gaia data (Bailer-Jones et al. 2021). Using Fourier analysis of AAVSO data, Kiss et al. (2006) determined periods of 495 ± 40 days and 3350 ± 1100 days for the star. Using self-correlation analysis, Percy & Sato (2009) determined a long secondary period of 2000 days. Chatys et al. (2019) found three periods: of 0.2 mag amplitude over 340 days, 0.28 mag over 495 days, and 0.46 over 3350 days. A study of the star by Kiss et al. (2010) with the Michigan InfraRed Combiner (MIRC) using four telescopes on the
<table><caption>**Table 1**
Observing Log for AZ Cyg</caption><thead><tr><th>Date (UT)</th><th>N<sub>telescopes</sub></th><th>N<sub>block</sub></th><th>Calibrators</th></tr></thead><tbody><tr><td>2011 Jul 21</td><td>6</td><td>1</td><td>7 And</td></tr><tr><td>2012 Jun 6</td><td>5</td><td>1</td><td>7 And</td></tr><tr><td>2012 Jun 11</td><td>6</td><td>1</td><td>σ Cyg</td></tr><tr><td>2014 Jul 16</td><td>5</td><td>1</td><td>σ Cyg, 7 And</td></tr><tr><td>2014 Jul 17</td><td>5</td><td>1</td><td>σ Cyg, 7 And</td></tr><tr><td>2015 Aug 22</td><td>5</td><td>2</td><td>σ Cyg</td></tr><tr><td>2015 Aug 23</td><td>6</td><td>2</td><td>σ Cyg</td></tr><tr><td>2016 Aug 29</td><td>6</td><td>1</td><td>σ Cyg</td></tr><tr><td>2016 Aug 30</td><td>6</td><td>1</td><td>σ Cyg</td></tr><tr><td>2016 Sep 3</td><td>6</td><td>1</td><td>σ Cyg</td></tr><tr><td>2016 Sep 6</td><td>5</td><td>1</td><td>σ Cyg</td></tr><tr><td>2016 Sep 8</td><td>5</td><td>1</td><td>σ Cyg, 7 And</td></tr><tr><td>2016 Sep 9</td><td>6</td><td>1</td><td>σ Cyg, 7 And</td></tr><tr><td>2016 Sep 10</td><td>6</td><td>4</td><td>σ Cyg, 7 And</td></tr><tr><td>2016 Sep 11</td><td>6</td><td>2</td><td>σ Cyg, 7 And</td></tr></tbody></table>

<table><caption>**Table 2**
Calibrators Used in This Study</caption><thead><tr><th>Star Identifier</th><th>UD (H -Band) (mas)</th><th>Source</th></tr></thead><tbody><tr><td>7 And</td><td>0.65 ± 0.03</td><td>Mourard et al. (2015)</td></tr><tr><td>σ Cyg</td><td>0.54 ± 0.02</td><td>Zhao et al. (2008)</td></tr></tbody></table>

CHARA Array suggested the presence of complex structures on the surface of the star. These results inspired this more involved study using the updated MIRC after it was upgraded to combine the light from six telescopes, which enabled imaging of the surface structures.

## 2.2. CHARA/MIRC Observations

Over the course of 5 yr, we observed AZ Cyg using MIRC at the Georgia State University's CHARA Array. The CHARA Array is located on Mount Wilson, CA and consists of six 1 m telescopes with 15 baselines ranging from 34 to 331 m and angular resolutions reaching ~0.5 mas in the H-band (ten Brummelaar et al. 2005).

The MIRC instrument, which has since been upgraded to MIRC-X, was a six beam combiner operating in the H-band (1.5–1.8 µm). The instrument was capable of providing 15 visibilities, 10 independent closure phases, and nine independent closure amplitudes per spectral channel. In the low spectral resolution mode used for our observations ($R=30$), eight 30 nm wide spectral channels are acquired (Monnier et al. 2004, 2012). We recorded 21 data blocks for AZ Cyg, as listed in Table 1, with each data block corresponding to a continuous observation of a target for 10 minutes. Note that we were unable to use data from 2012 in our analysis because of issues with calibration frames taken during the observation run.

To reduce data, we used the latest version of the official MIRC reduction pipeline (as of 2017 June; Monnier et al. 2007). In the pipeline, we applied a 17 millisecond coherence time and a cross-talk correction for visibilities less than 0.1. Corrections for variations in atmospheric coherence time and optical changes in the beam path were obtained using calibrator stars, as noted in Table 1 and described in Table 2. We also corrected for systematic errors following the method of Monnier et al. (2012). We used a 6.6% multiplicative error correction and $2 \times 10^{-4}$ additive error correction. For triple

amplitudes, we used a 10% multiplicative error correction and $1 \times 10^{-5}$ additive error correction. Following Zhao et al. (2011), we applied a 1° error floor for closure phases.

Because we do not expect significant brightness variations or rotations over periods of less than one month, we merged calibrated data observed within three weeks of each other into single data files for each epoch prior to analysis and image reconstruction. The purpose of this was to build up (u, v) coverage and increase the number of data to points to assist in image reconstruction. In Figure 1, we display the (u, v) coverage for the 2011 epoch. Figure A1 in the Appendix contains plots of (u, v) coverage for all epochs discussed in this paper. It was not always possible to record data with all six telescopes because of observing conditions or ongoing work on different telescopes.

## 2.3. SpeX IRTF Spectroscopy

In addition to the interferometric data, we obtained a spectrum of AZ Cyg on 2016 September 6 using the short cross dispersed (SXD) mode on the SpeX spectrograph on the NASA Infrared Telescope Facility (IRTF) 3 meter telescope atop Maunakea (Rayner et al. 2003). SpeX's SXD mode covers 0.7–2.5 µm at $R \sim 2000$ and is matched to a $0.3 \times 15''$ slit. To correct for additive systematic effects such as dark current and sky background, we employed an “AB” method in which we observed the target at two different positions along the slit. To build up signal-to-noise (S/N), we collected 5 AB pairs. The total integration time on the target was 13.902 s.

For purposes of flux calibration and correction of telluric features, we obtained 4 AB pairs of the A0 V standard HD 201320, which we selected using a list provided by the SpeX instrument team. The integration time on each exposure was 29.65 s. We used Spextool, which is an Interactive Data Language (IDL) package (Cushing et al. 2004), to reduce data. During this process, the standard star spectrum was used to generate a model telluric spectrum by adjusting a convolved model spectrum of Vega to match the standard star, as described in Vacca et al. (2003), and dividing the standard

<figure><img src="image_5.png" /><figcaption>**Figure 1.** (u, v) coverage of AZ Cyg 2011 observations.</figcaption></figure>
<table><caption>**Table 3**
Best Fit UD and Power-law LDD for AZ Cyg at 1.50–1.75 µm</caption><thead><tr><th>Parameter</th><th>2011</th><th>2014</th><th>2015</th><th>2016</th><th>All years</th></tr></thead><tbody><tr><td>θ<sub>UD</sub> (mas)</td><td>3.82 ± 0.01</td><td>3.81 ± 0.01</td><td>3.90 ± 0.01</td><td>3.99 ± 0.07</td><td>3.85 ± 0.01</td></tr><tr><td>χ<sup>2</sup><sub>UD</sub></td><td>16.08</td><td>41.05</td><td>23.23</td><td>22.49</td><td>24.68</td></tr><tr><td>θ<sub>LD</sub> (mas)</td><td>3.93 ± 0.01</td><td>4.09 ± 0.01</td><td>4.11 ± 0.01</td><td>4.09 ± 0.01</td><td>4.05 ± 0.01</td></tr><tr><td>α</td><td>0.42 ± 0.01</td><td>0.59 ± 0.02</td><td>0.62 ± 0.01</td><td>0.72 ± 0.02</td><td>0.57 ± 0.01</td></tr><tr><td>χ<sup>2</sup><sub>LD</sub></td><td>5.24</td><td>8.07</td><td>3.40</td><td>3.44</td><td>5.51</td></tr></tbody></table>

Note. Note that fits were made using only $V^2$ and not the other interferometric observables. $\chi^2$ are reduced $\chi^2$.

spectrum by this model so that only telluric lines remained. To get a spectrum that has been corrected for telluric lines, we divided the AZ Cyg spectrum by this derived telluric spectrum. We shifted the telluric spectrum roughly 0.1–0.2 pixels to minimize systematic errors around sharp telluric lines. Telluric correction using this method within the SpeX pipeline calibrates flux to 10% accuracy (Rayner et al. 2009). Calibration frames also included a flat field and Argon arc lamp exposure. In the final step of data reduction, we merged the order separated spectra into a single spectrum. In SpeX spectra, the offset between neighboring orders is between 1% and 3% (Rayner et al. 2009).

Following reduction, we followed the procedure of Villaume et al. (2017) to absolutely flux calibrate our data. In short, we determined a calibration factor,

$$ \langle C \rangle = \frac{\sum_x w_x \text{Max}[f(C_x)]}{\sum_x w_x}. \quad (1) $$

that we used to scale the entire spectrum, while preserving the relative flux calibration between each order. In Equation (1), $f(C_x) = 10^{(X_{AB}-N_{obs})/2.5}$ where $X_{AB}$ is the AB magnitude determined from the observed spectrum, $N_{obs}$ is a random normal distribution determined using 2MASS magnitudes and errors from Cutri et al. (2003), and $w_x$ are weights comprised of $\frac{1}{1\sigma}$ deviations from $f(C_x)$.

## 3. Stellar Parameters

### 3.1. Angular Diameter

We determined uniform disk (UD) and power law limb- darkened disk (LDD) diameters using $\chi^2$ minimization of models to the squared visibilities. The UD complex visibilities are given by the equation $V_{UD}(x) = \frac{2\ln(x)}{2}$ where $x = \pi\nu\theta_{UD}$, $J_1$ is a Bessel function of the first kind, $\nu$ is the spatial frequency at which the visibility function is sampled, and $\theta_{UD}$ is the angular diameter being fit. For the limb- darkened disk diameter, we used the Hestroffer law (Hestroffer 1997):

$$ I_{\lambda}(\mu)/I(1) = \mu^{\alpha}; \quad \alpha \in \mathbb{R}^{+} \quad (2) $$

where $I$ is intensity, $\mu = \sqrt{1 - (2r/\theta_{LD})^2}$ with $r$ the angular distance from the center of the star, $\theta_{LD}$ the angular diameter of the photosphere, and $\alpha$ is the limb darkening parameter. The complex visibility function corresponding to this law, as a function of $x = \pi\nu\theta_{LD}$, is:

$$ V(x) = \sum_{k \ge 0} \frac{(-1)^k \Gamma(\alpha/2 + 2)}{\Gamma(\alpha/2 + k + 2)\Gamma(k + 1)} \left(\frac{x}{2}\right)^{2k} \quad (3) $$

where $\Gamma$ is the Gamma function (Lacour et al. 2008).

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Squared visibilities for each epoch, denoted by different colors.</figcaption></figure>

We present the results of our fits in Table 3. Note that we did not observe any significant dependence on wavelength in the angular diameters when fitting using individual channels. The angular diameter and limb darkening parameter of the star seem to vary temporally. Because we did not limit our fits to the first lobe, the reduced $\chi^2$ values for fits to the uniform disk are high. This is to be expected for an object that deviates from a uniform disk, such as a star with spots on it and a darkened limb. Although lower, reduced $\chi^2$ values for fits to the limb- darkened disk are likewise impacted by the presence of such features. In Figure 2, we plot the visibility curve of each epoch, denoted by different colors. Looking at the first two lobes, which provide the diameter information, one can see that the differences in diameter arise from the second lobe and lower visibilities of the first lobe. This suggests that the variation in measured diameters comes from difficulties in fitting to the limb of the star, which would be sensitive to the presence of bright features.

### 3.2. Fundamental Stellar Parameters

To determine the fundamental stellar parameters of AZ Cyg, we made use of atlases of synthetic spectra based on spherical MARCS models (Gustafsson et al. 2008), SAtlas models (Lester & Neilson 2008), and synthetic spectra from the Lançon et al. (2007) spectral library, which is based on spherical PHOENIX models. Within the Lançon et al. (2007) grid, we selected the spectra with solar metallicity but $^{4}$He, $^{12}$C, $^{14}$N, and $^{16}$O abundances modified to match those of a $M = 20 M_{\odot}$ RSG near the end of its life. The MARCS grid consisted of 280 models, covering $T_{eff} = 3300–4500$ K, $\log g = -0.5 \pm 1.0$, and $[Fe/H] = 0.0 \pm 1.0$. All of the models were calculated for $M = 15 M_{\odot}$. To generate the model spectra from MARCS models we used TURBOSPECTRUM v19.1, a 1D LTE spectral synthesis code that incorporates numerous atomic and molecular features (Alvarez & Plez 1998; Plez 2012). To cover the complete range of our SpeX observations, we calculated spectra with $v_{turb} = 2$ km s$^{-1}$ and 0.1 Å steps from 650 to 5500 nm.
<table><caption>**Table 4**
Stellar Parameters of AZ Cyg Corresponding to the Best Fitting Model Atmospheres</caption><thead><tr><th>Model</th><th><i>T</i><sub>eff</sub><br>(K)</th><th>log <i>g</i></th><th>Radius<br>(<i>R</i><sub>☉</sub>)</th><th>Luminosity<br>(<i>L</i><sub>☉</sub>)</th><th>Mass<br>(<i>M</i><sub>☉</sub>)</th><th>[Fe/H]<br>(dex)</th><th><i>E</i>(<i>B</i> − <i>V</i>)<br>(mag)</th></tr></thead><tbody><tr><td>MARCS</td><td>4000</td><td>−0.50</td><td>1040</td><td>249443</td><td>15</td><td>0.0</td><td>0.56</td></tr><tr><td>PHOENIX</td><td>4000</td><td>0.00</td><td>641</td><td>94759</td><td>15</td><td>0.0</td><td>0.55</td></tr><tr><td>SATLAS</td><td>3972</td><td>−0.07</td><td>700</td><td>109828</td><td>15</td><td>0.0</td><td>0.54</td></tr></tbody></table>

To generate the grid of SA atlas model atmospheres we first calculated a grid of plane-parallel model atmospheres, starting with the initial grid of model atmospheres from Fiorella Castelli's website<fnref n="20" /> that covered the parameter space $T_{\text{eff}} = [3500, 3750, 4000, 4250, 4500]$ K, $log g = [0, 0.5, 1.0]$, $[\text{Fe/H}] = [-5.5:0.5]$ in 0.5 steps], and $v_{\text{turb}} = 2$ km s⁻¹ (Castelli & Kurucz 2004). We used the $\kappa_{\text{Rosseland}}$ coefficients, the “new” opacity distribution functions recomputed with an updated H₂O line list, and Kurucz's molecular and atomic line lists from Castelli's website. Using this grid, we calculated an expanded grid covering $T_{\text{eff}} = 3000-4500$ in 100 K steps and $log g = -1.0 \pm 1.0$ in 0.25 steps at $[\text{Fe/H}] = 0.0$ and $v_{\text{turb}} = 2$ km s⁻¹. Using this grid of plane-parallel models, we then calculated spherical models with radius $R = 100-900$ R☉, mass $M = 8, 10, 15, 20, 30$ M☉, and luminosity determined by using the corresponding radius and temperatures from 3000 to 4300 K in 100 K steps. Prior to calculation we checked whether the L, R, and M fell into RSG parameter space, which we defined as $10,000 L_{\odot} < L < 500,000 L_{\odot}$ and $-1.0 < \log g < 1.0$. This resulted in a grid of 542 model atmospheres. We then used the SYNTHE suite<fnref n="21" /> to generate model spectra at resolution $R = 500,000$ from 680 to 2700 nm using the Kurucz line lists as inputs (Kurucz 1993; Sbordone 2005).

We obtained stellar parameters by minimizing the $\chi^2$ between the observed spectra and our model spectra (after convolving to the resolution of SpeX) using the Amoeba function within IDL, which uses the Nelder-Mead downhill simplex method for multi-dimensional minimization. We included a radial velocity correction of $-3.88$ km s⁻¹ based on the measurement of the radial velocity from Gaia Data release 2 and corrected for reddening using a modified version of the Cardelli reddening law (Cardelli et al. 1989) that was proposed by Massey et al. (2005), setting $R_v = 4.2$ rather than the normal 3.1 to account for circumstellar dust. We allowed the color excess $E(B-V)$ to vary as a parameter in our fits. We limited our fit to 850-2290 µm to avoid strong TiO and CN lines that emerge further in the atmosphere of the star and would drive the temperature lower (Davies et al. 2013). We report the results of our fitting in Table 4 and present an example of the fits in Figure 3, with the rest given in Figure A2 in the Appendix. Bear in mind that the grids and input parameters for each model atmosphere atlas are somewhat different. Thus, in Table 4 luminosities and radii for the MARCS and PHOENIX rows are derived from the effective temperature, $log g$, and mass of the best fitting model. For the SA atlas row, effective temperature and $log g$ are derived from the luminosity, radius, and mass of the best fitting model. Note also that the luminosities for the SA atlas models were originally derived from a grid of radii and effective temperatures; as a

result, the luminosity grid does not correspond to rounded numbers.

We find that $T_{\text{eff}}$ ranges from 3972 to 4000 K and $log g$ from $-0.5$ to 0. The color excess we find, ranges from $E(B-V) = 0.54$ to 0.56. Gaia (Gaia Collaboration 2018; Gaia Collaboration et al. 2018) reports $T_{\text{eff}} = 3294$ K for AZ Cyg, which is close to the value of $T_{\text{eff}} = 3300$ K given by Lancon & Rocca-Volmerange (1992). The temperature scale of Levesque et al. (2005), based on fits of TiO rich optical spectra to spherical MARCS models, gives $T_{\text{eff}} = 3535-3660$ K for spectral types M2-4.5. Meanwhile, van Belle et al. (2009) gives $T_{\text{eff}} = 3513 \pm 168$ K to $3755 \pm 194$ K for spectral types M2-M4.7, based on angular diameters derived from the Palomar Testbest Interferometer data and fits of SEDs to templates from Pickles (1998). The cooler temperatures reported previously likely result from the inclusion of regions of cool molecules further in the atmosphere of the star. Davies et al. (2013) found that measurements of RSG temperatures using spectral regions with strong molecular lines resulted in cooler temperatures. Using the limb-darkened diameters and Gaia distances from Bailer-Jones et al. (2021) ($d = 2090_{-115}^{+130}$ pc), we find an average radius of $R = 911_{-50}^{+57}$ R☉.

# 4. Image Reconstruction

## 4.1. Reconstruction Codes

Interferometric imaging constitutes an ill-posed inverse problem, due to both sparse $(u, v)$ coverage and the incompleteness of the phase information carried by closure phases. The problem is typically solved by regularized maximum likelihood (Thiébaut & Young 2017). The reconstructed image $\mathbf{x}^{\text{opt}} \in \mathbb{R}^{N \times N}$ is defined as the most probable $N \times N$-pixel image given the data, under a set of prior assumptions about the image (such as a given level of smoothness or sparsity). The data are assumed to be normally distributed, following the OIFITS standard (Pauls et al. 2005). The priors are then imposed both by the optimization engine for image positivity, and by a regularization function $R(\mathbf{x})$, so that $\mathbf{x}^{\text{opt}}$ is given by:

$$ \mathbf{x}^{\text{opt}} = \underset{\mathbf{x} \in \mathbb{R}^{N \times N}}{\operatorname{argmax}} (\chi^2(\mathbf{x}) + \mu R(\mathbf{x})). \quad (4) $$

where $\mu$ is a hyperparameter that sets the relative strength of the regularization with respect to the $\chi^2$. In this paper, we use two very different imaging codes to produce the reconstructions. Both exemplify two different optimization techniques to solve Equation (4).

OITOOLS.jl<fnref n="22" /> uses a quasi-Newton method with positivity constraint (Thiebaut 2002), a fast minimization technique which makes use of the analytic expressions of the gradients $\frac{\partial \chi^2}{\partial x}$ and $\frac{\partial R(\mathbf{x})}{\partial x}$. However, OITOOLS.jl suffers from a major

<footnote n="20">²⁰ http://wwwuser.oats.inaf.it/castelli/</footnote>

<footnote n="21">²¹ Available in a Linux port at https://wwwuser.oats.inaf.it/castelli/sources.html.</footnote>

<footnote n="22">²² https://github.com/fabienbaron/OITOOLS.jl</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 3. Upper: best fitting PHOENIX model spectrum (red) compared to observed spectrum. Lower: residuals (in sigmas) of the observed and model spectrum.</figcaption></figure>

drawback because of its restricted choice of regularization function, which needs to be both differentiable and convex.

In contrast, SQUEEZE<fnref n="23" /> (Baron et al. 2010, 2012) employs slower, gradient-less, stochastic methods. In this paper, we used its default simulated annealing mode. The image **x** is modeled as a superposition of a large number of flux elements (typically 1000–10,000 elements). These elements are free to randomly move within the image, increasing or decreasing the $\chi^2 + \mu R$ functional as they do so and as the posterior distribution is explored. The evolution of the image as the elements move constitutes a Markov Chain. The chain stabilizes after a period of burn-in where the elements are settling: the posterior distribution is still being explored, but only close to the optimal solution. Averaging the chain over the last hundreds of iterations then provides a mean image that, while it does not strictly minimize Equation (4), still characterizes the posterior around the solution. A standard deviation is also produced, which quantifies the spread around the mean and thus reflects the posterior sharpness.

## 4.2. Regularization

The best non-committal regularization in optical interferometry is generally admitted to be total variation (Renard et al. 2011), which acts on image **x** by imposing greater sparsity on its spatial gradient $\nabla x$. The total variation functional $R_{TV}$ is the $\ell_1$ norm of the spatial gradient, and can be approximated by:

$$ R_{TV}(\mathbf{x}) = \|\nabla \mathbf{x}\|_1 \\ \approx \sum_{i,j} \sqrt{(x_{i,j} - x_{i-1,j})^2 + (x_{i,j} - x_{i,j-1})^2 + \epsilon} \quad (5) $$

where $\epsilon$ is chosen to be a small number close to numerical precision. This approximation is differentiable and can be readily used in OITOOLS.jl.

The effectiveness of the total variation can be traced to two main effects. First, it is very effective in regularizing the outside of the star (i.e., pixels with zero flux). As such, it acts as a soft-support constraint, as opposed to a hard-support constraint enforced (e.g., by a mask). Second, total variation prevents high frequency noise. Total variation, by design, also tends to favor zones of uniform flux separated by sharp boundaries. Since 3D simulations of stellar convection do show convection cells separated by black outline, total variation may be a good choice as long as $\mu$ remains reasonable.

In addition to total variation, we employ a Laplacian regularizer and $\ell_0$ regularizer in this study. The Laplacian regularizer is the $\ell_1$ norm of the Laplacian: $R_{LA}(\mathbf{x}) = \|\nabla^2 \mathbf{x}\|_1^2$. Because it favors contours, it is particularly useful in imaging stellar surfaces with non-uniform flux. The $\ell_0$ regularizer is a pseudo-norm, which is not a norm because it is non-convex:

$$ \ell_0(\mathbf{x}) = \sum_{i}^{N} 1_{\mathbb{R}}(x_i), \quad (6) $$

where $1_{\mathbb{R}}_{\ge 0}$ is the indicator function of $\mathbb{R}_{\ge 0}$, the set of real non-negative numbers. The $\ell_0$ pseudo-norm counts the number of non-zero elements in the vector it is applied to and hereby heavily penalizes stray flux outside of the stellar disk.

## 4.3. Results

Prior to reconstructing images of AZ Cyg, we ran an experiment to test the impact of different combinations of regularizers and hyperparameters on reconstructions of RSGs. To start, we selected a snapshot of a 3D RHD simulation from Chiavassa et al. (2011a), scaled the snapshot to the angular

<footnote n="23">²³ https://github.com/fabienbaron/squeeze</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Reconstruction of a simulated observation of the 3D RHD simulation of Betelgeuse from Chiavassa et al. (2009), using the (u, v) coverage of the AZ Cyg 2011 observation. The left-hand panel presents an image of the source simulation, the middle panel presents the source image scaled down and convolved to the resolution of CHARA, and the right-hand panel presents the reconstructed image using the simulated observation. In each image, the beam size in the lower right corner is the reconstruction resolution. The reconstructed image shows that it is possible to reconstruct the surface features predicted by 3D RHD models at the resolution of CHARA using typical observational coverage.</figcaption></figure>

diameter of AZ Cyg, and copied the (u, v) coverage from each epoch to make four OIFITS files of simulated observations of a RSG. Using each simulated observation, we ran a large stack of reconstructions with SQUEEZE using different regularizers and hyperparameters. In each reconstruction, we used a mask of a 4–4.25 mas disk and an initial image corresponding to the best fitting limb-darkened disk for a given epoch. Each reconstruction ran on five independent chains, and we aligned and averaged the results of all of the chains to produce a single final mean image. We used the $l_1$ norm to compare each mean image to the source snapshot, convolved to the resolution of the reconstruction and rebinned to match a 4.0 mas star. Gomes et al. (2017) found the $l_1$ norm to be the optimal metric for comparing reconstructions to source images.

$$l_1(\mathbf{x}, \mathbf{y}) = \|\mathbf{x} - \mathbf{y}\|_1 = \sum_i |x_i - y_i|. \quad (7)$$

For the reconstruction of images using real data from each epoch, we selected the regularizers and hyperparameters used in the reconstruction with the lowest $l_1$ norm of the simulated data corresponding to that epoch. In Figure 4, we depict the source image and the best fitting reconstruction for the simulated observation based on the 2011 observations. In Table 5 we present the reconstruction parameters we used for each epoch, based on the results of this experiment.

To produce our images, we ran SQUEEZE using the parameters from Table 5 on 10 chains for 500 iterations using 4500 elements on a 64 × 64 grid at a scale of 0.125 mas/pixel, which is roughly 1/8th the resolution of our data. In each reconstruction we used an initial image of a 4.0–4.25 mas disk, which we smoothed with a Gaussian to enforce artificially sharp boundaries at the edges, and the best fitting limb-darkened disk for a given epoch was used as a prior image.

The use of multiple chains offers a chance to investigate the impact of variations in the reconstruction process on the final image. For each epoch, we determined mean and standard deviation images from the result of all chains. We found that the standard deviation image was not very informative, likely because of the use of a mask and initial image, and that the result from each chain was only marginally different.

<table><caption>**Table 5**

Regularization Parameters Used for Reconstructing the Images in Figures 5, 7, and A4</caption><thead><tr><th>Year</th><th>Regularizers</th><th>Hyperparameters</th><th>χ<sup>2</sup></th></tr></thead><tbody><tr><td>2011</td><td>Total Variation (TV), l<sub>0</sub></td><td>1000,0.3</td><td>2.06</td></tr><tr><td>2014</td><td>Laplacian(LA), TV</td><td>2500,500</td><td>2.36</td></tr><tr><td>2015</td><td>LA, TV</td><td>2500,1500</td><td>4.09</td></tr><tr><td>2016</td><td>LA, TV</td><td>300,2000</td><td>4.80</td></tr></tbody></table>

**Note.** Note that the $\chi^2$ describes the reduced $\chi^2$ of fits of the reconstructions to all the interferometric observables.

We also tested our images for artifacts by simulating an observation of a uniform disk using the (u, v) coverage of our observations and reconstructing an image using exactly the same parameters. We depict the result of that experiment in Figure 5. Although the reconstructed images do have some noisy features, we note that in every case the patterns observed in the reconstructions of the true data are not present and that the pixel intensities are much lower than those found in the images from the real data.

In the second row of Figure 5, we present SQUEEZE reconstructed images of AZ Cyg from 2011 to 2016. In Figure 6, we present the squared visibilities and closure phases of the 2011 data compared to synthetic squared visibilities and closure phases generated from the mean reconstructions. A figure presenting a comparison of data from all of the epochs is available in Figure A3 in the Appendix. The reduced $\chi^2$ of the reconstructions, based on comparison of the squared visibilities, triple amplitude, and closure phases of the data, are listed in Table 5. The reduced $\chi^2$ in Table 5 and the residuals in the bottom panels of Figure A3 indicate that there is a close match between the reconstructed images, and the data and low frequencies. However, there are notable deviations at higher frequencies and in datasets with longer temporal coverage.

Although the spectral resolution of the observations is rather low, we also reconstructed images for each wavelength channel to test for wavelength dependence. These images were produced using the same reconstruction parameters that we used in reconstructing images using data from all of the
<figure><img src="image_9.png" /><figcaption>AZ Cyg (OITOOLS.jl) 2011</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>AZ Cyg (OITOOLS.jl) 2014</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>AZ Cyg (OITOOLS.jl) 2015</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>AZ Cyg (OITOOLS.jl) 2016</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>AZ Cyg (SQUEEZE) 2011</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>AZ Cyg (SQUEEZE) 2014</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>AZ Cyg (SQUEEZE) 2015</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>AZ Cyg (SQUEEZE) 2016</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>AZ Cyg 2011</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>AZ Cyg 2014</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>AZ Cyg 2015</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>AZ Cyg 2016</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>UD 2011</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>UD 2014</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>UD 2015</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>UD 2016</figcaption></figure>

**Figure 5.** First row: images of AZ Cyg from 2011 to 2016 produced with OITOOLS.jl. Second row: mean images of AZ Cyg from 2011–2016 produced with SQUEEZE. Third row: standard deviation of the reconstructed images in row two. Fourth row: corresponding reconstruction of a simulated uniform disk (UD). The intensities of the images in each column are scaled to the maximum pixel value of the corresponding mean image in row one. The beam size in the right-hand corner corresponds to resolution given by the maximum projected baseline of that observation.

wavelength channels. We present these images in Figure A4 in the Appendix.

To test the results of our reconstructions, we also used OITOOLS.jl to reconstruct the data using total variation and a centering prior. We used the same masks and initial images as used with SQUEEZE. We display the reconstructions made using OITOOLS.jl in the top row of Figure 5. We note that prominent features found in the images produced using SQUEEZE are present in the reconstructions made with OITOOLS.jl, which

suggests that these features are not the result of a bias in the reconstruction process. In the following analysis section, we will only use the images produced using SQUEEZE.

As further tests of our images, we split up the observations from the 2016 epoch into two parts and we then reconstructed images with these split data. We find that the prominent features in the image produced with the complete dataset remain in the split datasets. We depict the result of this experiment in Figure 7.
<figure><img src="image_25.png" /><figcaption>AZ Cyg 2011 Squared Visibilities</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>AZ Cyg 2011 Closure Phases</figcaption></figure>

**Figure 6.** Upper: comparison of AZ Cyg observations 2011 (black) to squared visibilities and closure phases calculated from the mean SQUEEZE image for the same epoch (red). Lower: residuals (in sigmas) between the quantities in each upper graph.

<figure><img src="image_27.png" /></figure>

2016 Aug 28-Sep 06

<figure><img src="image_28.png" /></figure>

2016 Sep 8-11

<figure><img src="image_29.png" /></figure>

2016 Aug 28-Sep 11

**Figure 7.** SQUEEZE reconstructions of AZ Cyg in 2016 using data from August 28 to September 6, September 8–11, and entirety of the AZ Cyg 2016 observations. The beam size in the right-hand corner corresponds to resolution given by the maximum projected baseline of that observation.

## 5. Discussion

In Figure 5, we can see that the photosphere of AZ Cyg is non-uniform, and consists of features of varying size and intensity. The accompanying reconstructions of uniform disks based on simulated observations using the same (u, v) coverage do not show similar features as the reconstructions of AZ Cyg. Thus, we posit that the features in the reconstructions of the star are not artifacts of observational coverage.

In Figure A4, it is clear that the images have a wavelength dependence. The spectrum of AZ Cyg is dominated by CN lines from roughly 1.450–1.560 µm, and there are several CO bands between 1.560 and 1.640 µm. Thus, it is likely that reconstructions in these channels display the appearance of the star at higher levels in the atmosphere where these lines form. Meanwhile, the channels between 1.670 and 1.760 µm cover a continuum region of the spectrum. Thus, these reconstructions in these channels present the surface of the star.

In each epoch, there are a number of relatively small (<30% $R_*$), bright features of various intensities, as well as larger medium intensity regions and smaller dark regions. The size of the small, bright features appears to be similar in every year. However, given that these features are less than the nominal resolution of the CHARA Array, it could be that this is the result of a limitation in our data. Meanwhile, the medium intensity regions seem to be consistent in their size, occupying >0.5 $R_*$, as was found in Chiavassa et al. (2010). The smaller, dark features are of comparable size to the bright regions.

The 3D RHD models, as well as spectroscopic and photometric evidence, suggest that the surface of RSGs is occupied by large, long lived convection cells with spatial extent ≥0.5 $R_*$. We expect there to be smaller granules with a lifetime on the order of months within these cells. In the reconstruction of the 3D RHD model presented in Figure 4, we see that the large convection cells in the source image appear as bright to medium intensity regions in the reconstruction, with
<table><caption>**Table 6**
Surface Contrast AZ Cyg</caption><thead><tr><th>Year</th><th>1.47–1.52 µm</th><th>1.52–1.56 µm</th><th>1.56–1.60 µm</th><th>1.60–1.63 µm</th><th>1.63–1.67 µm</th></tr></thead><tbody><tr><td>2011</td><td>20.69 ± 1.20</td><td>16.41 ± 1.07</td><td>15.07 ± 1.85</td><td>12.84 ± 0.08</td><td>12.81 ± 0.08</td></tr><tr><td>2014</td><td>6.27 ± 0.92</td><td>6.27 ± 0.92</td><td>6.27 ± 1.30</td><td>6.86 ± 0.52</td><td>4.83 ± 0.35</td></tr><tr><td>2015</td><td>10.50 ± 2.86</td><td>11.95 ± 1.09</td><td>11.11 ± 0.79</td><td>11.36 ± 0.70</td><td>12.35 ± 0.68</td></tr><tr><td>2016</td><td>18.27 ± 2.25</td><td>18.13 ± 2.74</td><td>14.23 ± 1.06</td><td>13.57 ± 2.79</td><td>14.55 ± 2.00</td></tr></tbody></table><table><thead><tr><th>Year</th><th>1.67–1.70 µm</th><th>1.70–1.73 µm</th><th>1.73–1.76 µm</th><th>1.47–1.76 µm</th><th>OITOOLS (1.47–1.76 µm)</th></tr></thead><tbody><tr><td>2011</td><td>12.37 ± 1.99</td><td>13.33 ± 3.22</td><td>11.93 ± 0.96</td><td>13.89 ± 0.18</td><td>8.84</td></tr><tr><td>2014</td><td>6.42 ± 1.77</td><td>6.69 ± 2.10</td><td>5.17 ± 1.10</td><td>9.59 ± 1.65</td><td>4.39</td></tr><tr><td>2015</td><td>10.13 ± 0.53</td><td>10.56 ± 1.04</td><td>9.81 ± 1.54</td><td>12.16 ± 0.79</td><td>9.40</td></tr><tr><td>2016</td><td>14.26 ± 1.22</td><td>11.72 ± 1.90</td><td>10.81 ± 2.10</td><td>14.46 ± 1.30</td><td>6.85</td></tr></tbody></table>

the bright regions being granules of rising hot gas. The darker regions in the reconstruction correspond to the thin lanes of infalling gas between the large convection cells. Although the source image has four or five large cells, in the reconstruction these regions blur together and it is difficult to distinguish more than three separate regions.

Carrying this same analysis to the reconstructions of AZ Cyg, we see a strikingly similar surface pattern. In 2011, the surface of AZ Cyg exhibited several distinct regions that may be convection cells bounded by dark lanes. One of these is in the northwest quadrant of the star, another is in the northeastern quadrant and one is in the southeastern quadrant. In addition, there are several bright regions that may be granules. The size of these granule-like features is at the limit of the interferometer, so these images provide an upper limit of $0.25 R_*$ on the granules.

The image from 2011 is separated too far in time from the 2014 image to make any temporal comparison. However, comparing the 2014–2016 images, allows us to make a rough estimate of the lifetime of these features. In 2014, there is a region in the northwest that may be a convection cell, a feature in the south that may be a separate cell, and a feature in the east that may be another. In 2015, the bright granules are no longer in the same location, as expected, but a medium intensity region remains in the northwest and another remains in the south. The similarities between 2015 and 2016 are even more striking. The pattern in the southwest in the 2015 image remains in the 2016 image. This suggests that the lifetime of the larger features is at least one year and possibly two. The rotation periods of most RSGs are on the order of 20 yr. Thus, it is not likely that yearly changes on the surface are the result of long lived features moving out of view. To better understand the lifetime of surface features will require data obtained at a shorter time interval, perhaps monthly.

To quantitatively analyze our images, we first determined the root mean square (rms) relative intensity contrast (which we will call contrast throughout the rest this paper) of the surface granulation using the definition of Tremblay et al. (2013), which has found use by Wittkowski et al. (2017), Paladini et al. (2018) and Montargès et al. (2018) in studies of evolved stars.

$$ \frac{\delta I_{\text{rms}}}{\langle I \rangle} = \left\langle \frac{\sqrt{\langle I(x, y, t)^2 \rangle_{x,y} - \langle I(x, y, t) \rangle_{x,y}^2}}{\langle I(x, y, t) \rangle_{x,y}} \right\rangle_t \quad (8) $$

As noted in Tremblay et al. (2013), this quantity is a measure of the deviation of a star from plane-parallel approximation and those 3D models with the highest contrast have the least

<figure><img src="image_30.png" /><figcaption>**Figure 8.** Surface contrasts of AZ Cyg for different wavelength channels in each epoch of observation as determined using the mean reconstructed images in Figure 5.</figcaption></figure>

efficient convection (Trampedach et al. 2013). The contrast is correlated with the Mach number, which is the ratio of the flow and sound speeds, such that stars with high contrast ratios will have larger Mach numbers. This is explained by mixing length theory, which suggests that larger convective velocities correspond to larger temperature fluctuations between convective features and their surroundings. The intensity contrast ratio is correlated with density, with lower densities corresponding to higher Mach numbers. Because lower temperature stars have lower densities, there is also a correlation between temperature and intensity contrast, with lower temperatures corresponding to larger intensity contrast ratios. This relationship is expressed quantitatively (for giant stars) by Trampedach et al. (2013):

$$ I_{\text{rms}}/\% = (20.81 \pm 1.81) \times \log T_{\text{eff}} - (1.11 \pm 0.14) \times \log g - 55.46 \pm 0.29. \quad (9) $$

We present the average intensity contrasts in Table 6 and Figure 8. To get these values, we took the average of the contrast measured on the mean image of each chain. This permits our measured values to incorporate the spread resulting from the reconstruction process. We also performed this measurement on the OITOOLS.jl reconstruction, which does not have an error value because the reconstruction process results in only one image.

We find that our contrast ranges from roughly 5.0%–14.5% in the continuum wavelength channels. Contrasts in observations from 2014 are notably lower than those in other years. It is unclear if this is the result of an actual difference in the star, or the lower quality of that year’s data. Using Equation (9) with
<figure><img src="image_31.png" /></figure>

2011

<figure><img src="image_32.png" /></figure>

2014

<figure><img src="image_33.png" /></figure>

2015

<figure><img src="image_34.png" /></figure>

2016

<caption>Figure 9. Power spectra of the mean reconstructed images in Figure 5, with the limb-darkened disk edge filtered out.</caption>

the parameters from Table 4 gives 19.35% ± 6.52% through 20.05 ± 6.53. Paladini et al. (2018) found contrasts of 11.9% ± 0.4% to 13.1% ± 0.2% for an AGB of $T_{\text{eff}} = 3200$ K and $\log g = -0.4$. Montargès et al. (2018) found contrasts of 5%–6% for the RSG CE Tau with $T_{\text{eff}} = 3820 \pm 135$ K and $\log g = 0.05_{-0.17}^{+0.11}$. Wittkowski et al. (2017) found a contrast of 10% ± 4% for the RSG V776 Cen with $T_{\text{eff}} = 4290 \pm 760$ K. Thus, our contrasts fall within the range of those found around similar stars. Based on the models described by Tremblay et al. (2013), which we note do not cover RSG parameter space, these intensity contrasts correspond roughly to Mach numbers between 0.25 and 0.75 at unity Rosseland optical depth ($\langle\tau_R\rangle_{x,y} = 1$). In models of red giant stars, Tremblay et al. (2013) found relative intensity contrasts of roughly 18% and 22%. Although these models have similar temperature to AZ Cyg, they are of higher gravity ($\log g = 1.0$ and $1.5$). A lower gravity would suggest a lower density and thus larger Mach number and larger intensity contrast, this trend is also described by Equation (9). We find a lower rather than higher contrast, which suggests some difference between red giants and red supergiants. If we focus on the Mach number as a guide, then one explanation is that the convective flow velocity behaves differently. A lower contrast suggests a lower Mach number and thus a lower convective flow velocity, which could result if the depth of the cells was larger in RSGs than in red giants. A lower convective velocity would also suggest a larger advective timescale for the same pressure scale height and thus longer cell lifetimes.

## 5.1. Spatial Scales

Tremblay et al. (2013) and Trampedach et al. (2013) also offer a way of measuring the characteristic size of features using the power spectra of the intensity maps they generated. This method was adapted by Paladini et al. (2018) to measure the characteristic size of granulation on the AGB $\pi^1$ Gruis using images generated from Very Large Telescope Interferometer (VLTI) data. We followed a similar method to derive the size of features on AZ Cyg. We first took the 1D spatial power spectra of each mean image. We then filtered the disk out of the images by setting the pixel values of an image from the disk boundary outwards equal to the median flux of the pixels within the star. This had the effect of smoothing out the sharp boundary at the edge of the disk and it enabled detection of prominent power carrying features within the disk. To make identification of the peak easier, we padded this filtered image by 200 pixels. We identified the characteristic size of granulation as the remaining peak in a power spectrum. We present the power spectrum of each epoch in Figure 9. We also present a comparison of power spectra of images produced with SQUEEZE and OITOOLS.jl in Figure 10. We note the characteristic size of granulation in Table 7.

Using mixing length theory with turbulence, Antia et al. (1984) predicted a small number of large granules on the surface of cool, evolved stars, such as red giants and red supergiants. This was in line with the qualitative arguments by Schwarzschild (1975), which extrapolated from solar data to suggest that the atmospheric depth of the maximum of the
<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /><figcaption>**Figure 10.** Power spectra of the mean reconstructed images from OITOOLS in Figure 7 compared to the reconstructions from SQUEEZE in Figure 5, both with the limb-darkened disk edge filtered out.</figcaption></figure>

<table><caption>**Table 7**
Characteristic Granulation Size for AZ Cyg ($R_*$)</caption><thead><tr><th>Year</th><th>1.47–1.52 µm</th><th>1.52–1.56 µm</th><th>1.56–1.60 µm</th><th>1.60–1.63 µm</th><th>1.63–1.67 µm</th><th>1.67–1.70 µm</th></tr></thead><tbody><tr><td>2011</td><td>0.74 ± 0.06</td><td>0.52 ± 0.04</td><td>0.67 ± 0.16</td><td>0.51 ± 0.04</td><td>0.48 ± 0.03</td><td>0.56 ± 0.08</td></tr><tr><td>2014</td><td>0.52 ± 0.09</td><td>0.46 ± 0.11</td><td>0.79 ± 0.30</td><td>0.71 ± 0.14</td><td>0.77 ± 0.06</td><td>1.13 ± 0.03</td></tr><tr><td>2015</td><td>0.66 ± 0.11</td><td>0.78 ± 0.12</td><td>0.83 ± 0.12</td><td>0.80 ± 0.24</td><td>1.00 ± 0.34</td><td>1.14 ± 0.16</td></tr><tr><td>2016</td><td>0.67 ± 0.17</td><td>0.57 ± 0.05</td><td>0.54 ± 0.04</td><td>0.53 ± 0.14</td><td>0.57 ± 0.06</td><td>0.63 ± 0.07</td></tr></tbody></table><table><thead><tr><th>Year</th><th>1.70–1.73 µm</th><th>1.73–1.76 µm</th><th>1.47–1.76 µm</th><th>OITOOLS (1.47–1.76 µm)</th></tr></thead><tbody><tr><td>2011</td><td>0.70 ± 0.18</td><td>0.59 ± 0.04</td><td>0.69 ± 0.01</td><td>1.06 ± 0.01</td></tr><tr><td>2014</td><td>0.94 ± 0.28</td><td>0.90 ± 0.25</td><td>0.67 ± 0.19</td><td>0.8 ± 0.01</td></tr><tr><td>2015</td><td>1.11 ± 0.06</td><td>1.16 ± 0.17</td><td>0.52 ± 0.03</td><td>0.51 ± 0.02</td></tr><tr><td>2016</td><td>0.90 ± 0.23</td><td>0.78 ± 0.26</td><td>0.54 ± 0.02</td><td>0.51 ± 0.02</td></tr></tbody></table>

superadiabatic temperature gradient was related to the vertical extent of a granule and thus also related to the horizontal size of such a feature. Furthermore, 2D and 3D models (Freytag et al. 1997; Trampedach et al. 2013; Tremblay et al. 2013) have found that granule size, $x_{gran}$ is proportional to the scale pressure height, $H_p$, a result of the fact that the mixing length, the length at which a convective element dissipates into its environment, is proportional to $H_p$ in current theories of convection. Interestingly, Freytag et al. (1997) found that the superadiabatic temperature gradient was less useful as a descriptor of $x_{gran}$ because it deviates from a linear relation at lower gravities. Instead, they proposed the relation

$H_p = RT_{eff}/g$ with $R$ the universal gas constant. This can be rewritten (Paladini et al. 2018) as

$$ \frac{x_{\text{gran}}}{R_{*}} \approx 0.0025 \frac{R_{*} T_{\text{eff},*} M_{\odot}}{R_{\odot} T_{\text{eff},\odot} M_{*}} \quad (10) $$

Tremblay et al. (2013) found that the relation between $x_{gran}$ and $H_p$ was not linear but varied based on Mach number. They noted that $x_{gran}/H_p$ might be related to the ratio between horizontal and vertical velocity, with larger velocity ratios corresponding to larger ratios of $x_{gran}/H_p$. To describe the relationship between stellar parameter and $x_{gran}$, they determined a parameterization using least squares fits to the granule
<table><caption>**Table 8**
Calculated Surface Properties of AZ Cyg</caption><thead><tr><th>Year</th><th>x<sub>gran.</sub>, obsv<br>(R<sub>*</sub>)</th><th>x<sub>gran.</sub>, Freytag<br>(R<sub>*</sub>)</th><th>x<sub>gran.</sub>, Tremblay<br>(R<sub>*</sub>)</th><th>x<sub>gran.</sub>, Trampedach<br>(R<sub>*</sub>)</th><th>x<sub>gran.</sub>, Chiavassa<br>(R<sub>*</sub>)</th></tr></thead><tbody><tr><td>2011</td><td>0.69 ± 0.01</td><td>0.10<sup>+0.01</sup><sub>-0.01</sub></td><td>0.04–0.19</td><td>0.09–0.38</td><td>0.51<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>2014</td><td>0.67 ± 0.19</td><td>0.11<sup>+0.01</sup><sub>-0.01</sub></td><td>0.04–0.19</td><td>0.09–0.37</td><td>0.53<sup>+0.04</sup><sub>-0.03</sub></td></tr><tr><td>2015</td><td>0.52 ± 0.03</td><td>0.11<sup>+0.01</sup><sub>-0.01</sub></td><td>0.04–0.19</td><td>0.08–0.37</td><td>0.53<sup>+0.04</sup><sub>-0.03</sub></td></tr><tr><td>2016</td><td>0.54 ± 0.02</td><td>0.11<sup>+0.01</sup><sub>-0.01</sub></td><td>0.03–0.18</td><td>0.08–0.35</td><td>0.53<sup>+0.04</sup><sub>-0.03</sub></td></tr></tbody></table>

Note. When ranges are reported, they represent the lower and upper values determined using the parameters in Table 4.

size described by a power spectrum:

$$ \frac{\text{Char.size}}{[\text{km}]} = 13.5g^{-1}[T_{\text{eff}} - 300 \log g]^{1.75} 10^{0.05[\text{Fe/H}]}, \quad (11) $$

with *g* in cgs and *T*<sub>eff</sub> in K. Also using least square fits to characteristic sizes determined from power spectra of 3D models, Trampedach et al. (2013) found that the granulation size *A*<sub>gran</sub> was described by the relation

$$ \log \frac{A_{\text{gran}}}{[\text{Mm}]} \approx (1.310 \pm 0.0038) \times \log T_{\text{eff}} - (1.0970 \pm 0.003) \times \log g + 0.0306 \pm 0.0359. \quad (12) $$

Finally, using 3D RHD models of red supergiants, Chiavassa et al. (2009) found the equation of Freytag et al. (1997) needed to account for turbulent pressure. Including this factor, they used

$$ H_p = \frac{kT_{\text{eff}}}{g\mu m_{\text{H}}} \left( 1 + \beta\gamma \left( \frac{v_{\text{turb}}}{c_s} \right)^2 \right) \quad (13) $$

with *β* a parameter close to one, *γ* the adiabatic exponent, *c*<sub>s</sub> the sound speed, *μ* the average particle mass, *m*<sub>H</sub> the atomic mass of hydrogen, and *v*<sub>turb</sub> the turbulent velocity. Using this equation and models of RSGs, Chiavassa et al. (2009) found that convective feature size could be obtained by multiplying Equation (10) by five. The feature size predicted by this equation matched that of the large convection features in their 3D RHD models. Chiavassa et al. (2011a) conducted further analysis of the relationship between pressure scale height and large surface features, this time as suggested by the standard deviation of photocenter displacement. They found that both 3D RHD models and interferometric observations of RSGs suggest that the relationship between convection cell size and pressure scale height is different in RSGs than it is in giant stars.

In Table 8, we display the parameters of the surface of AZ Cyg derived from Equations (10)–(13), which we determined using values from Table 4. We report ranges in calculated granule size for Equations (11) and (12) because these equations depend on log *g*, for which we found a wide-range of possible values, depending on the stellar atmosphere code used in the fit (see Table 4). We did not use averages of the parameters in Table 4 because they resulted from very different stellar atmosphere codes. Based on these results, it appears that the features measured in the images of AZ Cyg presented in this paper are the same as those calculated by Equation (13). These features are likely to be the longer lived convection cells. If we assume that these cells last 3350 days,

which is the length of the longest period observed by Kiss et al. (2006) and Chatys et al. (2019), then the same cell could be present during the entirety of the observing cycle. The similarity in measured granule size between 2015 and 2016 suggests that features of that scale do last longer than a year. The difference between those years and the measured size in 2011 may indicate the emergence of a new cell, or it could indicate that there was a change in the size of a cell during its lifetime. Furthermore, long-term observations of the star, which we are pursuing, will be necessary to image the full life-time of a convection cell. Meanwhile, the features derived using Equations (10)–(12) could correspond to the smaller bright features visible in the images in Figure 5, which may represent smaller scale granulation with a shorter lifetime. These features may last on the order months, which indicates a need for short-term observations and also long-term study.

## 6. Conclusion

The observations presented here have allowed us to derive the fundamental stellar parameters of the RSG AZ Cyg. In addition, we have reconstructed H-band images of the star, finding that the surface varies from year to year with large features of scale ~0.5 *R*<sub>*</sub>, as predicted by 3D RHD simulations. We also find small bright spots on the surface, which are perhaps akin to the granulation predicted by such models. However, these features are of a scale that is below the resolution of the CHARA Array, and thus analysis of their size is beyond the scope of this work. Although it seems as if the larger features last for longer than one year, the surface pattern of the small features varies substantially. Observations on a shorter timescale, on the order of months, will be necessary to get a better understanding of the lifetime of these features. Overall, it seems that the surface of this RSG, and likely others, is dominated by a pattern of large, long lived convection cells and smaller features consisting of hot granules of rising gas that persist for months, but last less than one year.

The authors would like to thank the anonymous referee for the helpful comments and edits. R.N. would like to thank Dr. Carlos Lopez Carrillo for the suggestion to use mean images from each chain when measuring contrast and feature size to obtain an understanding of the spread in measurements.

F.B. and R.N. would like to acknowledge funding of this research through NSF Awards #1616483 and #1814777. R.N. would also like to acknowledge start-up funding from the Office for Research at New Mexico Tech. This work is based upon observations obtained with the Georgia State University Center for High Angular Resolution Astronomy Array at Mount Wilson Observatory. The CHARA Array is supported by the National Science Foundation under grant No. AST-
1636624 and AST-1715788. Institutional support has been
provided from the GSU College of Arts and Sciences and the
GSU Office of the Vice President for Research and Economic
Development.

This work is based upon observations obtained with the
Infrared Telescope Facility, which is operated by the University
of Hawaii under contract NNH14CK55B with the National
Aeronautics and Space Administration.

This research has made use of the Jean-Marie Mariotti
Center OIFits Explorer<fnref n="24" /> and Aspro<fnref n="25" /> services.

This research has also made use of the SIMBAD database
and the VizieR catalog access tool, CDS, Strasbourg, France
(DOI:10.26093/cds/vizier). The original description of the
VizieR service was published in Ochsenbein et al. (2000).

This publication makes use of data products from the Two
Micron All Sky Survey, which is a joint project of the
University of Massachusetts and the Infrared Processing and
Analysis Center/California Institute of Technology, funded by
the National Aeronautics and Space Administration and the
National Science Foundation.

This work has made use of data from the European Space
Agency (ESA) mission Gaia (https://www.cosmos.esa.int/
gaia), processed by the Gaia Data Processing and Analysis
Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/
dpac/consortium). Funding for the DPAC has been provided
by national institutions, in particular the institutions participat-
ing in the Gaia Multilateral Agreement.

**Appendix**
**Additional Material**

This appendix contains additional figures as referenced in the
text. Figure A1 presents (u, v) coverage for all of the epochs
described in the paper. Figure A2 presents comparisons of the
observed and model spectra with residuals. Figure A3 presents
comparisons between the squared visibilities and closure
phases of the observed data to the synthetic visibilities and
closure phases generated from mean images reconstructed
using SQUEEZE and residuals. Figure A4 presents mean
images reconstructed for each wavelength channel of MIRC
using SQUEEZE.

<footnote n="24"><sup>24</sup> Available at http://www.jmmc.fr/oifitsexplorer.</footnote>

<footnote n="25">²⁵ Available at http://www.jmmc.fr/aspro.</footnote>
<figure><img src="image_39.png" /></figure>

2011

<figure><img src="image_40.png" /></figure>

2014

<figure><img src="image_41.png" /></figure>

2015

<figure><img src="image_42.png" /></figure>

2016

<caption>**Figure A1.** (u, v) coverage of AZ Cyg for four epochs over a five year time period. Different baseline combinations are identified by colors.</caption>
<figure><img src="image_43.png" /><figcaption>MARCS fit</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>PHOENIX fit</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>SATLAS fit</figcaption></figure>

**Figure A2.** Best fitting model spectrum (red) compared to observed spectrum.
<figure><img src="image_46.png" /><figcaption>AZ Cyg 2011 Squared Visibilities</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>AZ Cyg 2014 Squared Visibilities</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>AZ Cyg 2011 Closure Phases</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>AZ Cyg 2014 Closure Phases</figcaption></figure>

**Figure A3.** Upper: comparison of AZ Cyg observations (black) to squared visibilities and closure phases calculated from the mean SQUEEZE images (red). Lower: residuals (in sigmas) between the quantities in each upper graph.
<figure><img src="image_50.png" /><figcaption>AZ Cyg 2015 Squared Visibilities</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>AZ Cyg 2016 Squared Visibilities</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>AZ Cyg 2015 Closure Phases</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>AZ Cyg 2016 Closure Phases</figcaption></figure>

**Figure A3.** (Continued.)
<figure><img src="image_54.png" /></figure>

AZ Cyg 2011
1.47µm-1.52µm

<figure><img src="image_55.png" /></figure>

AZ Cyg 2014
1.47µm-1.52µm

<figure><img src="image_56.png" /></figure>

AZ Cyg 2015
1.47µm-1.52µm

<figure><img src="image_57.png" /></figure>

AZ Cyg 2016
1.47µm-1.52µm

<figure><img src="image_58.png" /></figure>

AZ Cyg 2011
1.52µm-1.56µm

<figure><img src="image_59.png" /></figure>

AZ Cyg 2014
1.52µm-1.56µm

<figure><img src="image_60.png" /></figure>

AZ Cyg 2015
1.52µm-1.56µm

<figure><img src="image_61.png" /></figure>

AZ Cyg 2016
1.52µm-1.56µm

<figure><img src="image_62.png" /></figure>

AZ Cyg 2011
1.56µm-1.60µm

<figure><img src="image_63.png" /></figure>

AZ Cyg 2014
1.56µm-1.60µm

<figure><img src="image_64.png" /></figure>

AZ Cyg 2015
1.56µm-1.60µm

<figure><img src="image_65.png" /></figure>

AZ Cyg 2016
1.56µm-1.60µm

<figure><img src="image_66.png" /></figure>

AZ Cyg 2011
1.60µm-1.63µm

<figure><img src="image_67.png" /></figure>

AZ Cyg 2014
1.60µm-1.63µm

<figure><img src="image_68.png" /></figure>

AZ Cyg 2015
1.60µm-1.63µm

<figure><img src="image_69.png" /></figure>

AZ Cyg 2016
1.60µm-1.63µm

**Figure A4.** Channel separated images of AZ Cyg. The beam size in the right-hand corner represents the resolution given by the maximum projected baseline of each dataset.
<figure><img src="image_70.png" /><figcaption>AZ Cyg 2011
1.63µm-1.67µm</figcaption></figure>

<figure><img src="image_71.png" /><figcaption>AZ Cyg 2014
1.63µm-1.67µm</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>AZ Cyg 2015
1.63µm-1.67µm</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>AZ Cyg 2016
1.63µm-1.67µm</figcaption></figure>

<figure><img src="image_74.png" /><figcaption>AZ Cyg 2011
1.67µm-1.70µm</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>AZ Cyg 2014
1.67µm-1.70µm</figcaption></figure>

<figure><img src="image_76.png" /><figcaption>AZ Cyg 2015
1.67µm-1.70µm</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>AZ Cyg 2016
1.67µm-1.70µm</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>AZ Cyg 2011
1.70µm-1.73µm</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>AZ Cyg 2014
1.70µm-1.73µm</figcaption></figure>

<figure><img src="image_80.png" /><figcaption>AZ Cyg 2015
1.70µm-1.73µm</figcaption></figure>

<figure><img src="image_81.png" /><figcaption>AZ Cyg 2016
1.70µm-1.73µm</figcaption></figure>

<figure><img src="image_82.png" /><figcaption>AZ Cyg 2011
1.73µm-1.76µm</figcaption></figure>

<figure><img src="image_83.png" /><figcaption>AZ Cyg 2014
1.73µm-1.76µm</figcaption></figure>

<figure><img src="image_84.png" /><figcaption>AZ Cyg 2015
1.73µm-1.76µm</figcaption></figure>

<figure><img src="image_85.png" /><figcaption>AZ Cyg 2016
1.73µm-1.76µm</figcaption></figure>

Figure A4. (Continued.)

ORCID iDs

Ryan P. Norris https://orcid.org/0000-0002-9120-9728

Fabien R. Baron https://orcid.org/0000-0002-8376-8941

John D. Monnier https://orcid.org/0000-0002-3380-3307

Arturo O. Martinez https://orcid.org/0000-0002-3311-4085

Gail H. Schaefer https://orcid.org/0000-0001-5415-9189

Andrea Chiavassa https://orcid.org/0000-0003-3891-7554

Christopher D. Farrington https://orcid.org/0000-0001-9939-2830

Douglas R. Gies https://orcid.org/0000-0001-8537-3583

László L. Kiss https://orcid.org/0000-0002-3234-1374

Miguel Montargès https://orcid.org/0000-0002-7540-999X

Hilding R. Neilson https://orcid.org/0000-0002-7322-7236

Ettore Pedretti https://orcid.org/0000-0002-6214-0282

Stephen T. Ridgway https://orcid.org/0000-0003-2557-7132

Rachael M. Roettenbacher https://orcid.org/0000-0002-9288-3482

Nicholas J. Scott https://orcid.org/0000-0003-1038-9702

Theo A. ten Brummelaar https://orcid.org/0000-0002-0114-7915

References

Alvarez, R., & Plez, B. 1998, A&A, 330, 1109

Antia, H. M., Chitre, S. M., & Narasimha, D. 1984, ApJ, 282, 574

Arroyo-Torres, B., Wittkowski, M., Chiavassa, A., et al. 2015, A&A, 575, A50

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147

Baron, F., Monnier, J. D., Kiss, L. L., et al. 2014, ApJ, 785, 46

Baron, F., Monnier, J. D., & Kloppenborg, B. 2010, Proc. SPIE, 7734, 77342I

Baron, F., Monnier, J. D., Pedretti, E., et al. 2012, ApJ, 752, 20
Buscher, D. F., Haniff, C. A., Baldwin, J. E., & Warner, P. J. 1990, MNRAS, 245, 7P

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Castelli, F., & Kurucz, R. L. 2004, arXiv:astro-ph/0405087

Chatys, F. W., Bedding, T. R., Murphy, S. J., et al. 2019, MNRAS, 487, 4832

Chiavassa, A., Freytag, B., Masseron, T., & Plez, B. 2011a, A&A, 535, A22

Chiavassa, A., Haubois, X., Young, J. S., et al. 2010, A&A, 515, A12

Chiavassa, A., Pasquato, E., Jorissen, A., et al. 2011b, A&A, 528, A120

Chiavassa, A., Plez, B., Josselin, E., & Freytag, B. 2009, A&A, 506, 1351

Climent, J. B., Wittkowski, M., Chiavassa, A., et al. 2020, A&A, 635, A160

Cotton, D. V., Bailey, J., Horta, A. D., Norris, B. R. M., & Lomax, J. R. 2020, RNAAS, 4, 39

Cushing, M. C., Vacca, W. D., & Rayner, J. T. 2004, PASP, 116, 362

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, yCat, 2246, 0

Davies, B., Kudritzki, R.-P., Plez, B., et al. 2013, ApJ, 767, 3

De Beck, E., Decin, L., de Koter, A., et al. 2010, A&A, 523, A18

Dharmawardena, T. E., Mairs, S., Scicluna, P., et al. 2020, ApJL, 897, L9

Dupree, A. K., Strassmeier, K. G., Matthews, L. D., et al. 2020, ApJ, 899, 68

Freytag, B., Holweger, H., Steffen, M., & Ludwig, H.-G. 1997, in Science with the VLT Interferometer, ed. F. Paresce (Berlin: Springer), 316

Gaia Collaboration 2018, yCat, 1345, 0

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Gomes, N., Garcia, P. J. V., & Thiébaut, É. 2017, MNRAS, 465, 3823

Guinan, E. F., & Wasatonic, R. J. 2020, ATel, 13410, 1

Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951

Harvey, J. 1985, in ESA Special Publication 235, Future Missions in Solar, Heliospheric & Space Plasma Physics, ed. E. Rolfe & B. Battrick (Paris: ESA), 199

Haubois, X., Perrin, G., Lacour, S., et al. 2009, A&A, 508, 923

Hestroffer, D. 1997, A&A, 327, 199

Josselin, E., & Plez, B. 2007, A&A, 469, 671

Joyce, M., Leung, S.-C., Molnár, L., et al. 2020, ApJ, 902, 63

Kiss, L. L., Monnier, J. D., Bedding, T. R., et al. 2010, in ASP Conf. Ser. 425, Hot and Cool: Bridging Gaps in Massive Star Evolution, ed. C. Leitherer et al. (San Francisco, CA: ASP), 140

Kiss, L. L., Szabó, G. M., & Bedding, T. R. 2006, MNRAS, 372, 1721

Kravchenko, K., Chiavassa, A., Van Eck, S., et al. 2019, A&A, 632, A28

Kurucz, R. L. 1993, SYNTHE Spectrum Synthesis Programs and Line Data (Cambridge, MA: Smithsonian Astrophysical Observatory)

Lacour, S., Meimon, S., Thiébaut, E., et al. 2008, A&A, 485, 561

Lançon, A., Hauschildt, P. H., Ladjal, D., & Mouhcine, M. 2007, A&A, 468, 205

Lancon, A., & Rocca-Volmerange, B. 1992, A&AS, 96, 593

Lester, J. B., & Neilson, H. R. 2008, A&A, 491, 633

Levesque, E. M., & Massey, P. 2020, ApJL, 891, L37

Levesque, E. M., Massey, P., Olsen, K. A. G., et al. 2005, ApJ, 628, 973

López Ariste, A., Mathias, P., Tessore, B., et al. 2018, A&A, 620, A199

Massey, P., Plez, B., Levesque, E. M., et al. 2005, ApJ, 634, 1286

Mauron, N., & Josselin, E. 2011, A&A, 526, A156

Meynet, G., Chomienne, V., Ekström, S., et al. 2015, A&A, 575, A60

Monnier, J. D., Berger, J.-P., Millan-Gabet, R., & ten Brummelaar, T. A. 2004, Proc. SPIE, 5491, 1370

Monnier, J. D., Che, X., Zhao, M., et al. 2012, ApJL, 761, L3

Monnier, J. D., Zhao, M., Pedretti, E., et al. 2007, Sci, 317, 342

Montargès, M., Cannon, E., & Lagadec, E. 2021, Natur, 594, 365

Montargès, M., Kervella, P., Perrin, G., et al. 2014, A&A, 572, A17

Montargès, M., Norris, R., Chiavassa, A., et al. 2018, A&A, 614, A12

Mourard, D., Monnier, J. D., Meillard, A., et al. 2015, A&A, 577, A51

Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23

Paladini, C., Baron, F., Jorissen, A., et al. 2018, Natur, 553, 310

Pauls, T. A., Young, J. S., Cotton, W. D., & Monnier, J. D. 2005, PASP, 117, 1255

Percy, J. R., & Sato, H. 2009, JRASC, 103, 11

Pickles, A. J. 1998, PASP, 110, 863

Plez, B. 2012, Turbospectrum: Code for Spectral Synthesis, v19.1, Astrophysics Source Code Library, ascl:1205.004

Rayner, J. T., Cushing, M. C., & Vacca, W. D. 2009, ApJS, 185, 289

Rayner, J. T., Toomey, D. W., Onaka, P. M., et al. 2003, PASP, 115, 362

Ren, Y., & Jiang, B.-W. 2020, ApJ, 898, 24

Renard, S., Thiébaut, E., & Malbet, F. 2011, A&A, 533, A64

Sbordone, L. 2005, MSAIS, 8, 61

Schwarzschild, M. 1975, ApJ, 195, 137

Stothers, R. 1972, A&A, 18, 325

ten Brummelaar, T. A., McAlister, H. A., Ridgway, S. T., et al. 2005, ApJ, 628, 453

Thiebaut, E. 2002, Proc. SPIE, 4847, 174

Trampedach, R., Asplund, M., Collet, R., Nordlund, Å., & Stein, R. F. 2013, ApJ, 769, 18

Tremblay, P. E., Ludwig, H. G., Freytag, B., Steffen, M., & Caffau, E. 2013, A&A, 557, A7

Tuthill, P. G., Haniff, C. A., & Baldwin, J. E. 1997, MNRAS, 285, 529

Vacca, W. D., Cushing, M. C., & Rayner, J. T. 2003, PASP, 115, 389

van Belle, G. T., Creech-Eakman, M. J., & Hart, A. 2009, MNRAS, 394, 1925

van Loon, J. T., Cioni, M.-R. L., Zijlstra, A. A., & Loup, C. 2005, A&A, 438, 273

Villaume, A., Conroy, C., Johnson, B., et al. 2017, ApJS, 230, 23

Wittkowski, M., Abellán, F. J., Arroyo-Torres, B., et al. 2017, A&A, 606, L1

Thiébaut, É., & Young, J. 2017, JOSAA, 34, 904

Young, J. S., Baldwin, J. E., Boysen, R. C., et al. 2000, MNRAS, 315, 635

Zhao, M., Gies, D., Monnier, J. D., et al. 2008, ApJL, 684, L95

Zhao, M., Monnier, J. D., Che, X., et al. 2011, PASP, 123, 964