One month convection timescale on the surface of a giant evolved star

The transport of energy through convection is important during many stages of stellar evolution1,2, and is best studied in our Sun3 or giant evolved stars4. Features that are attributed to convection are found on the surface of massive red supergiant stars5–8. Also for lower-mass evolved stars, indications of convection are found9–13, but convective timescales and sizes remain poorly constrained. Models indicate that convective motions are crucial to produce strong winds that return the products of stellar nucleosynthesis into the interstellar medium14. Here we report a series of reconstructed interferometric images of the surface of the evolved giant star R Doradus. The images reveal a stellar disk with prominent small-scale features that provide the structure and motions of convection on the stellar surface. We find that the dominant structure size of the features on the stellar disk is 0.72 ± 0.05 astronomical units. We measure the velocity of the surface motions to vary between −18 and +20 km s−1, which means that the convective timescale is approximately one month. This indicates a possible difference between the convection properties of low-mass and high-mass evolved stars.

The transport of energy through convection is important during many stages of stellar evolution [1,2], and is best studied in our Sun [3] or giant evolved stars [4].Features that are attributed to convection are found on the surface of massive red supergiant stars [5][6][7][8].Also for lower mass evolved stars, indications of convection are found [9][10][11][12][13], but convective timescales and sizes remain poorly constrained.Models indicate that convective motions are crucial for the production of strong winds that return the products of stellar nucleosynthesis into the interstellar medium [14].
Here we report a series of reconstructed interferometric images of the surface of the evolved giant star R Doradus.The images reveal a stellar disc with prominent small scale features that provide the structure and motions of convection on the stellar surface.We find that the dominant structure size of the features on the stellar disc is 0.72 ± 0.05 astronomical units (au).We measure the velocity of the surface motions to vary between −18 and +20 km s −1 , which means the convective timescale is approximately one month.This indicates a possible difference between the convection properties of low-mass and high-mass evolved stars.The M-type asymptotic giant branch (AGB) star R Doradus (see also the Methods section 1) was observed with the Atacama Large Millimeter/submillimeter Array (ALMA) in five epochs spread over four weeks between 2023 July 2 and August 2.The observations around 338 GHz were done using the longest available ALMA baselines.Using superuniform weighthing, this allowed us to reconstruct images of the stellar surface at an angular resolution between 8 − 25 mas (see Method section 1).At submillimeter wavelengths, the main opacity source are the electron-neutral free-free interactions [15].As a result, the observations are sensitive to the temperature and density structure in an extended atmosphere above the stellar photosphere.In particular, submillimeter observations probe the dynamics and properties of the shocks excited by large convective cells on the stellar photosphere, while not being hindered by dust or molecular opacity sources that dominate at other wavelengths [11].The surface maps for the three observational epochs with the highest angular resolution are shown in Fig. 1, the remaining epochs are shown in the supplementary material.The maps show a predominantly circular stellar disc with a radius of R 338 GHz = 1.64 ± 0.09 au and a brightness temperature of T b = 2270 ± 130 K (see Methods section 1).Considering the correspondence between several epochs, including observations taken at 225 GHz shown in the Methods section 1, we conclude that the structures that we observe are intrinsic to the star and are very likely induced by surface convection.From the visual inspection, we find that the structures we observe at the surface have a typical lifetime of at least three weeks.
In order to estimate the size of the dominant structures on the submillimeter surface, we employ a spatial power spectral density (PSD) analysis (see the description in the Method section 1).The PSD analysis was performed using the calibrated interferometric visibilities after subtracting the best fit model for the stellar disc in the three last, highest angular resolution, epochs.We show the results in Fig. 2. We find that structure exists at several scales on the surface.The structure varies between the three epochs, particularly in amplitude, and most of the power is concentrated in small scale structure with a size of x = 13.1 ± 0.  The spatial power spectrum density for three epochs of observations of R Doradus.The spatial power spectrum density (PSD), in units of mJy 2 , was determined directly from the interferometric visibilities after subtracting the best fit model for the stellar disc.The three panels represent the final three, highest angular resolution, epochs.The solid circles denote the measurements and the curve represents a cubic spline fit to the observations.The largest blue stars, with corresponding 1σ s.d.error bars, are the PSD determined on-source, while the smaller orange dots represent the off-source measurements.In the last two epochs, the second bin at ∼ 12 mas contains only limited visibilities and hence does not present a significant detection.As the peak at the smallest angular size dominates, this indicates that the smallest granule structure dominates.From an error weighted average of the three epochs, we find a value of x = 13.1 ± 0.6 mas.For a distance of 55 ± 3 pc [16] this corresponds to x = 0.72 ± 0.05 au, compared to a stellar diameter at 338 GHz of 3.28 au.55 ± 3 pc [16], this corresponds to a size of x = 0.72 ± 0.05 au.Other peaks present in the power spectrum are located at 16.0 ± 1.0, 17.9 ± 0.6, 22.1 ± 1.3, and 35.2 ± 3.6 mas.The largest of these correspond to more than half the visible hemisphere.The different scales might be related to the equivalent scales of granulation, mesogranulation and, considering the scales, in particular supergranulation inferred for our Sun [17,18], or could be the result of a superposition of several independent granules.In the absence of a better name we subsequently use the term granules for the smallest structures, although concluding a direct link with Solar granules is not yet possible.On average, the power on each of the four larger scales is ∼ 25% of that in the dominant structure.
The structures inside the stellar disc can also be compared with structures seen on the limb of the star.In Fig. 3a, we show the half-power radius as a function of position angle for the three highest angular resolution epochs.There is a good correspondence between the epochs, with specifically the last two epochs showing very similar structures in the half-power radius in Fig. 3a.As shown in the supplementary information 1, also the 225 GHz observations show similar features confirming that they are intrinsic to the star.By determining the difference in radius between consecutive epochs at every position angle, we can derive the average velocity profile of the 338 GHz optical depth surface as shown in Fig. 3. Negative velocities indicate a radial motion inward to the star.While the average velocity is small, the movement of the optical depth surface as a function of angle varies between V = −18 and +20 km s −1 .These velocities can be compared with the local sound speed V s ∼6 km s −1 in the part of the extended atmosphere probed by the observations and are consistent with supersonic shocks induced by convection [e.g.22].Although the observed motions depend on a combination of changes in gas density, temperature, ionisation and velocity, they are, on short timescales, a good representation of the (radial) motions of the shocks.They thus represent a lower limit to the actual gas velocities.The escape velocity at the measured submillimeter radius for a star with the parameters adopted for R Doradus is V esc = 30 ± 3 km s −1 .Hence, as a consequence of these shocks, only a very small fraction of the gas will be able to escape the gravity field of the star before being accelerated by radiation pressure on dust [14].
Considering the measurement of the size of the granules as well as the velocity of the shocks, we can determine the timescale for the surface structures to readjust [e.g.20].This timescale is given by t surf ∼ x/∆V = 33±3 days, which is independent of the assumed distance to the star.This represents the first measurement of this timescale, which corresponds well with the timescale of more than three weeks estimated from the visual inspection of the images and is consistent with model predictions [21] but is only ∼ 30% of the convective decay time extrapolated using the emperical formula (Equation 16) from [22].As this emperical formula was derived for models of stars with both higher effective temperatures and surface gravity, our result points to a possible difference of convection properties for AGB stars.
The contrast between the average brightness of the stellar disk and the granules induced by the surface convective motions varies between 2.8, 3.2, and 1.5% during the last three epochs (Fig. 2).This is notably less than the ∼ 12% contrast observed for the star π 1 Gruis in the near-infrared.As the contrast relates to the imprint of the convective motions on the shock structure in the extended atmopshere, we can use the contrast and size to determine the average brightness temperature increase ∆T b caused by these shocks.We find ∆T b to be between ∼ 700 and 1500 K.As the brightness temperature is closely related to the real temperature of the extended atmosphere had a contrast in excess of 5% which, because of their compact size (< 0.3 au), corresponded to an increase in brightness temperature that can reach ∆T b > 50000 K.These hotspots have different characteristics from the convective shock structures described here, and are more rare.Further monitoring observations are needed to determine the timescales related to the bright hotspots and characterise their origin.
Although in the future, our observations can serve as a benchmark for the models of convection in M-type AGB stars, it is currently only possible to compare the Gruis.The other lines (dotted [26], dashed [30], and dash-dotted [22]) represent parametric models from the literature, with the thick and thin lines calculated for the stellar temperatures of 2710 and 3200 K, respectively.For the model from [22] we assume solar metallicity.
results with parametric formulas that predict the granular size scale based on fundamental stellar properties such as effective temperature, surface gravity and chemical composition.All relations are based on mixing length theory and extrapolated models of less evolved stars, and there are no models yet that specifically aim to reproduce R Doradus.In Fig. 4, we compare our measurements (using the stellar parameters described in the Method section 1) with the three relations that were also compared to the measurements of π 1 Gruis [10].While there is generally good agreement, none of the relations is a perfect fit.The size of the granulation is generally thought to scale with the pressure scale height H p immediately below the photosphere [4], with x = αH p .For AGB stars, the scaling parameter α is assumed to be ∼ 10 [26].For the more massive red supergiant stars, α ∼ 30 − 50 appears to better match the hydrodynamical simulations and observations [8,27].Using H p = k B T /µg, where k B is the Boltzmann constant, µ is the mean molecular mass and g is the acceleration due to surface gravity, this leads to the relation: Using the two AGB measurements, we can estimate a value of α = 17 +5 −3 that represents the best fit to both observations.
Our results indicate that the convective structure on the surface of AGB stars matches the existing parametric formulas, in size but reveal a difference between AGB stars and the more massive red supergiants.The results match the timescales found in AGB models, but the timescales differ from extrapolations of models based on less evolved stars and thus serve as a unique benchmark for the existing theory.The derived convection properties can also be implemented in stellar evolution and population synthesis models where convection is generally captured in a free parameter that has a large effect on galaxy evolution models [e.g.28].

Source properties
R Doradus is an M-type AGB star that belongs to the class of semi-regular pulsators.On a time-scale of ∼ 1000 days it switches between two pulsation modes with periods of 362 and 175 days.The distance to R Doradus has been determined to be 55 ± 3 pc using revised Hipparcos measurements [16].There is no useable Gaia parallax.From CO observations, it was determined that R Doradus has a relatively low mass loss rate (∼ 10 −7 M ⊙ yr −1 ) and wind expansion velocity (∼ 5.7 km s −1 ) [1]).Previous ALMA observations also indicated that R Doradus rotates fast for a giant star, with a rotation velocity at the surface of ∼ 1.0 ± 0.1 km s −1 compared to a rotation velocity of a few 10s of m s −1 expected for solitary AGB stars [19].It has been suggested that the apparent rotation could be the result of a chance alignment of convective cells [20].However, the rotation has been observed in multiple molecular lines at at least four different epochs which, including the observations presented here, span more than six years [e.g.12,19].This is much longer than the convective timescales found in our analysis, and hence a chance alignment of convective cells can be ruled out.
In our comparison with the convective theory, we adopt the values for the effective temperature of T eff = 2710 ± 70 K.For the surface gravity we use log(g) = −0.6 ± 0.1, based on models that indicate the initial mass was 1 − 1.25 M ⊙ and that the current mass is 0.7 − 1.0 M ⊙ , combined with interferometric measurements that yield a stellar diameter in the infrared of D IR = 51.18± 2.24 mas [29].It is expected that this diameter, which corresponds to a radius of R IR = 1.4±0.1 au = 298±21 R ⊙ , indicates the size of the stellar photosphere.We can compare this with the (τ = 1) size of the star determined with ALMA at 338 GHz obtained by visibility fitting.Using a combination of the last three epochs, we fit a nearly circular stellar disc with F 338GHz = 521 ± 18 mJy, D 338GHz = 59.8 ± 0.4 mas and an axis ratio of 0.99 ± 0.01.This means a brightness temperature T b = 2270 ± 130 K and, taking into account the uncertainty on the distance, a radius R 338GHz = 1.64 ± 0.09 au = 353 ± 19 R ⊙ = 1.18 ± 0.11 R IR .

Observations, data reduction and imaging
The AGB star R Doradus was observed in ALMA Bands 6 and 7 as part of ALMA project 2022.1.01071.S (PI: Khouri).The Band 7 observations were taken between July 5 and August 2 2023 using four spectral windows (spws) centered at 331.2, 333.0, 342.1, and 345.1 GHz.Each spw had a bandwidth of 1.875 GHz and 1920 channels.The integration time of the individual visibilities was set to 2.02 s.The observations were taken in the largest ALMA configurations (C-9 and C-10) with the quasars J0519-4546 and J0516-6207 as bandpass/amplitude and phase calibrator, respectively.Details of the observations are presented in Table 1.The calibration of the last three epochs was performed using the ALMA pipeline in CASA v6.4.1.12[2].The first two epochs were labelled as semi-pass in the ALMA quality assurance and for these the calibration was performed manually by staff from the Nordic ALMA regional center node using CASA v6.5.4.9.In the first epoch, there was an issue with the bandpass calibration that needed to be solved manually.In the second epoch, one of the antennas needed to be flagged, resulting in a loss of some of the longest baselines.For both epochs, the requested angular resolution and sensitivity were not reached.After the initial calibration of each epoch, molecular lines were identified and flagged before the data was averaged to an integration time of 6.06 s and to 50 channels per spw.Subsequently, two steps of phase-only self-calibration was performed on the stellar continuum.The self-calibration improved the signal-to-noise ratio (snr) by a factor of ∼ 2.5 on the continuum.Finally, images, using all four spw, were produced for the five epochs using superuniform visibility weighting [3].This method increases the relative weight of the visibilities at the longer baselines, which minimises the beamsize at the expense of snr.The superuniform beam characteristics and continuum rms noise are also given in Table 1.The increase of rms noise compared to more regular Briggs weighting (with a robust parameter of 0.5) depends on the telescope distribution and was a factor of ∼ 2, 1.5, 1.5, 3, and 7 for the five epochs, respectively.The improvement in angular resolution between superuniform and uniform weighting ranged from ∼ 2% in the final epoch to ∼ 15% in the third epoch.The superuniform weighted images of the final three epochs are presented in Fig. 1 and the first two epochs in Extended Data Fig. 1.The images of the highest angular resolution epochs were used to derive the angular radial profile presented in Fig. 3.We verified that reducing the angular resolution to match that of the epoch with the largest beam smooths out the observed structures and the derived velocities and thus use the highest angular resolution results.The fits of the stellar disc and the spatial power spectral density analysis were performed on the calibrated visibilities.While we focus our analysis on the higher resolution Band 7 observations, we also include the continuum result from the Band 6 observations in the Methods section 1.The observational details for these observations are also included in Table 1 and the calibration, self-calibration, and imaging steps were identical to those performed for the Band 7 observations.The four spws are centered at 218.9, 220.8, 230.0, and 232.9 GHz and the increase in the continuum rms noise between Briggs weighting and superuniform weighting is a factor 1.5.

Spatial Power Spectral Density analysis
The spatial power spectral density (PSD) is regularly used to derive information about e.g. the turbulent structure of the interstellar medium [e.g.4,5] as well as the convective structure of the solar photosphere [e.g.6].The spatial PSD is given by the 2D Fourier transform of an image.But considering interferometric images are themselves the Fourier transform of the interferometric visibilities, the spatial PSD is equal to the modulus squared of the complex visibilities [4].We can thus calculate the PSD directly from our interferometric visibilities without introducing potential artifacts during the imaging process.Because the PSD would be dominated by the power at the scales of the stellar disc, we first use the uv−fitting code uvmultifit [7] to fit the stellar disc (for a discussion on the disc profile see Section 1).Subsequently, we subtract the disc from the visibilities after which we calculate, for a phase center towards the center of the star, the PSD using visibilities annularly averaged in equally spaced bins of uv−distance (in units of kλ).From the inverse of the uv−distance we can directly obtain the angular scale x in milliarcseconds.In addition, we also determine the PSD for a position offset from the star.We chose an offset of 7 arcseconds to avoid a contribution to the off-source PSD from the source signal.Since the first two epochs have worse spatial resolution, we only produce the PSD for the final three epochs.The results are shown in Fig. 2. By comparing the on-and off-source PSD we can determine which structures are significant and which are possibly due to correlated noise in the visibilities.

Stellar disc profile
In order to check how well a top-hat shaped stellar disc profile fits the observations we have also investigated the stellar disc profile in the image plane.We produced a radially averaged profile of R Doradus based on the combined data for the final three epochs.We then compared this profile with a top-hat shaped stellar disc model convolved with the interferometric beam.The results of this comparison are shown in Extended Data Fig. 2. The stellar disc model can accurately describe the observations with residuals at a level of ∼ 2% of the peak emission.This means that the 338 GHz opacity τ 338GHz increases steeply over only a small change in radius.The radial motions we observe thus reflect the physical motion of the 338 GHz optical depth surface.Since the optical depth is a strong function of the density [11], the motions closely reflect the motion of the shocks induced by the convection.

Band 6 observations
The ALMA Band 6 (225 GHz) observations and the radius as a function of position angle are shown in Extended Data Fig. 3 4, it is clear that there is a very good correspondence.Since the different observations are completely independent, this shows that the pattern seen in the radii is intrinsic to the source at the time of the observations.Extended Data Figure 4: As Fig. 3a, the half-power radius of R Doradus as a function of the position angle.The curves indicate the half-power radius measured for the final epochs of the ALMA Band 7 (338 GHz) observations compared with the ALMA Band 6 (225 GHz) observations taken 15 days later.The red bar is the granule size derived from the higher resolution Band 7 observations.The improved signal-tonoise in Band 6 compensates for the lower flux and larger beam so that the typical 1σ s.d.radius uncertainty is also of order 0.2 mas.

Fig. 1 :
Fig. 1: The stellar surface of the AGB star R Doradus.The panels represent the three highest angular resolutions epochs of ALMA observations at 338 GHz.The black ellipse in the panels indicates the average size of the stellar disc at this frequency.The red solid contours and the blue dashed contours indicate the positive and negative 4, 5, and 6σ features with respect to the mean emission of the stellar disc.The size and orientation of the interferometric beam is indicated at the bottom left of each panel.
6 mas.At a distance of R Doradus of

Fig. 2 :
Fig.2: The spatial power spectrum density for three epochs of observations of R Doradus.The spatial power spectrum density (PSD), in units of mJy 2 , was determined directly from the interferometric visibilities after subtracting the best fit model for the stellar disc.The three panels represent the final three, highest angular resolution, epochs.The solid circles denote the measurements and the curve represents a cubic spline fit to the observations.The largest blue stars, with corresponding 1σ s.d.error bars, are the PSD determined on-source, while the smaller orange dots represent the off-source measurements.In the last two epochs, the second bin at ∼ 12 mas contains only limited visibilities and hence does not present a significant detection.As the peak at the smallest angular size dominates, this indicates that the smallest granule structure dominates.From an error weighted average of the three epochs, we find a value of x = 13.1 ± 0.6 mas.For a distance of 55 ± 3 pc[16] this corresponds to x = 0.72 ± 0.05 au, compared to a stellar diameter at 338 GHz of 3.28 au.

Fig. 3 :
Fig. 3: The radius and radial velocity of the stellar surface of R Doradus.Panel a The half-power radius R as a function of the position angle (P.A.) with respect to the north celestial pole.Positive angles point in the direction of the right ascension.The curves indicate R measured for three epochs of ALMA Band 7 (338 GHz) observations.The horizontal lines indicate the average radius for each epoch.The lower figure boundary on the radius is the estimated photospheric radius of 25.6 mas, measured at 2.3µm [29].The red vertical bar indicates the granule size measured from the spatial power spectral density, with the linear size translated to an angular size at R. The 1σ s.d. on the radius (similar for each epoch) is plotted in the bottom right corner (magenta errorbar) and is at most 0.17 mas.Panel b The radial velocity V at the surface of R Doradus.The velocity is determined from the difference of radius between the third and fourth epoch (18 and 27 July 2023; solid line) and fourth and fifth epoch (27 July and 2 Aug 2023; dashed line).The horizontal lines denote the average velocity.The velocity determined in this way corresponds to the movement of the τ 338 GHz = 1 optical depth surface and is an average between the respective epochs.This indicates the movement of the shocks induced by the convective motions.The 1σ s.d. in the radius determination translates in an uncertainty on the velocity of ∼ < 2.6 km s −1 and is indicated in the top right of the panel (magenta errorbar).

Fig. 4 :
Fig.4: The smallest size of the granules on the surface of R Doradus determined using ALMA compared with theoretical predictions and the previous VLTI measurement for π 1 Gruis[10].The observations are plotted vs. log(g) and the 1σ s.d.error bars include the uncertainty on the distance and stellar mass.The solid black lines indicate the parametric scaling presented in this paper.The thick line was calculated for a stellar effective temperature T eff = 2710 K, which is the value for R Doradus.The thin line is the same model for T eff = 3200 K, which is closer to the value for π 1 Gruis.The other lines (dotted[26], dashed[30], and dash-dotted[22]) represent parametric models from the literature, with the thick and thin lines calculated for the stellar temperatures of 2710 and 3200 K, respectively.For the model from[22] we assume solar metallicity.
and Extended Data 4 respectively.Fitting the interferometric visibilities yields a completely circular stellar disc with F 225GHz = 221.4± 0.1 mJy and D 225GHz = 61.8± 0.1 mas.This means a brightness temperature T b = 2006 ± 8 K and, including the distance uncertainty, a radius R 225GHz = 1.70 ± 0.09 au = 365 ± 20 R ⊙ = 1.22 ± 0.11 R IR .In a comparison between the radii of the 225 GHz observations and the last epochs of the 338 GHz observations in Fig.

1 : 2 :
The stellar surface of the AGB star R Doradus.The panels represent the first two epochs of ALMA observations at 338 GHz.The black ellipse in the panels indicates the average size of the stellar disc at this frequency.The size and orientation of the interferometric beam is indicated at the bottom left of each panel.The observed azimuthally averaged radial profile of the stellar disc of R Doradus.The radial profile, in the top panel, is generated from a combination of the three last epochs.The dashed line indicates the uniform disc model that minimizes the residual flux shown in the bottom panel.The error bars on the radial profile are close to the line width, with a largest error of 5 × 10 −5 Jy beam −1 .The error bars on the residual fraction are indicated with the vertical line segments.

3 :
The stellar surface of R Doradus at 225 GHz.The observations were taken using ALMA band 6 on August 17 2023.The black ellipse in the panels indicates the average size of the stellar disc at 225 GHz.

Table 1 :
Extended Data Observational details 1 Superuniform weighting 2 Data marked as semi-pass in ALMA data quality assessment3Observations in ALMA Band 6