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

The M-type asymptotic giant branch (AGB) star R Doradus (see also the first section of the Methods ) was observed with the Atacama Large Millimeter/submillimeter Array (ALMA) in five epochs spread over 4 weeks between 2 July and 2 August 2023. The observations around 338 GHz were made using the longest available ALMA baselines. Using superuniform weighting, this allowed us to reconstruct images of the stellar surface at an angular resolution between 8 and 25 mas (see the first section of the Methods ). At submillimetre wavelengths, the main opacity source is the electron-neutral free–free interactions ¹⁵ . As a result, the observations are sensitive to the temperature and density structure in an extended atmosphere above the stellar photosphere. In particular, submillimetre 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 ¹¹ . The surface maps for the three observational epochs with the highest angular resolution are shown in Fig. 1 ; the remaining epochs are shown in Extended Data Fig. 1 . The maps show a predominantly circular stellar disk with a radius of R _338GHz = 1.64 ± 0.09 astronomical units ( au ) and a brightness temperature of T _b = 2,270 ± 130 K (see the first section of the Methods ). Considering the correspondence between several epochs, including observations taken at 225 GHz shown in the first section of the Methods , we conclude that the structures that we observe are intrinsic to the star and are probably 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 3 weeks.

To estimate the size of the dominant structures on the submillimetre surface, we use a spatial power spectral density (PSD) analysis (see the description in the first section of the Methods ). The PSD analysis was carried out using the calibrated interferometric visibilities after subtracting the best fit model for the stellar disk 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.6 mas. At a distance of R Doradus of 55 ± 3 pc (ref. ¹⁶ ), 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 corresponds 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 about 25% of that in the dominant structure.

The structures inside the stellar disk 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 Extended Data Fig. 4 , also the 225 GHz observations show similar features, confirming that they are intrinsic to the star. By determining the difference in the 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 inwards to the star. Although 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 ⁻¹ . These velocities can be compared with the local sound speed V _s ≈ 6 km s ⁻¹ in the part of the extended atmosphere probed by the observations and are consistent with supersonic shocks induced by convection (for example, ref. ¹⁹ ). Although the observed motions depend on a combination of changes in gas density, temperature, ionization 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 submillimetre radius for a star with the parameters adopted for R Doradus is V _esc = 30 ± 3 km s ⁻¹ . Hence, because 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 ¹⁴ . 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 (for example, ref. ²⁰ ). This timescale is given by t _surf ≈ x /∆ V = 33 ± 3 days, which is independent of the assumed distance to the star. The measurement of this timescale corresponds well with the timescale of more than 3 weeks estimated from the visual inspection of the images and is consistent with model predictions ²¹ but is only about 30% of the convective decay time extrapolated using the empirical formula (equation (16)) from ref. ¹⁹ . As this empirical 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 ± 1.2, 3.2 ± 0.2 and 1.5 ± 0.2% during the last three epochs (Fig. 2 ). This is notably less than the contrast of about 12% observed for the star π ¹ Gruis in the near-infrared. As the contrast relates to the imprint of the convective motions on the shock structure in the extended atmosphere, 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 about 700 and 1,500 K. As the brightness temperature is closely related to the real temperature of the extended atmosphere (for example, refs. ^{11 , 15} ), this is sufficient to affect the chemistry in the stellar atmosphere (for example, ref. ²² ). Previous observations have also revealed brighter and more compact hotspots on the submillimetre surface of AGB stars, including R Doradus (for example, refs. ^{11 , 23 , 24} ). These hotspots had a contrast in excess of 5% that, because of their compact size (<0.3 au ), corresponded to an increase in brightness temperature that can reach ∆ T _b > 50,000 K. These hotspots have different characteristics from the convective shock structures described here and are rarer. Further monitoring observations are needed to determine the timescales related to the bright hotspots and characterize their origin. Although in the future, our observations can serve as a benchmark for the models of convection in M-type AGB stars, at present it is possible only to compare the results with parametric formulae that predict the granular size scale on the basis of 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 first section of the Methods ) with the three relations that were also compared to the measurements of π ¹ Gruis ¹⁰ . Although 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 ⁴ , with x = αH _p . For AGB stars, the scaling parameter α is assumed to be about 10 (ref. ²⁵ ). For the more massive red supergiant stars, α ≈ 30–50 seems to better match the hydrodynamical simulations and observations ^{8 , 26} . Using H _p = k _B T / µg , in which 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:

giving a result in astronomical units. Using the two AGB measurements, we can estimate a value of \(\alpha ={17}_{-3}^{+5}\) 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 formulae 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 in which convection is generally captured in a free parameter that has a large effect on galaxy evolution models (for example, ref. ²⁷ ).

Source properties

R Doradus is an M-type AGB star that belongs to the class of semi-regular pulsators. On a timescale of about 1,000 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 ⁶ . There is no usable Gaia parallax. From CO observations, it was determined that R Doradus has a relatively low mass loss rate (about 10 ⁻⁷ M _⊙ yr ⁻¹ ) and wind expansion velocity (about 5.7 km s ⁻¹ ) ³⁰ . Previous ALMA observations also indicated that R Doradus rotates fast for a giant star, with a rotation velocity at the surface of about 1.0 ± 0.1 km s ⁻¹ compared to a rotation velocity of a few tens of metres per second expected for solitary AGB stars ³¹ . It has been suggested that the apparent rotation could be the result of a chance alignment of convective cells ²⁰ . However, the rotation has been observed in several molecular lines at four different epochs that, including the observations presented here, span more than 6 years (for example, refs. ^{12 , 31} ). 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 = 2,710 ± 70 K. For the surface gravity, we use log[ g ] = −0.6 ± 0.1, based on models that indicate that 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 (ref. ²⁸ ). 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 disk 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 = 2,270 ± 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 the ALMA project 2022.1.01071.S (principal investigator: T.K.). The band 7 observations were taken between 5 July and 2 August 2023 using four spectral windows centred at 331.2, 333.0, 342.1 and 345.1 GHz. Each spectral window had a bandwidth of 1.875 GHz and 1,920 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 Extended Data Table 1 . The calibration of the last three epochs was carried out using the ALMA pipeline in CASA v6.4.1.12 ³² . The first two epochs were labelled as semi-pass in the ALMA quality assurance and for these the calibration was carried out 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 were averaged to an integration time of 6.06 s and to 50 channels per spectral window. Subsequently, two steps of phase-only self-calibration were carried out on the stellar continuum. The self-calibration improved the signal-to-noise ratio by a factor of about 2.5 on the continuum. Finally, images, using all four spectral windows, were produced for the five epochs using superuniform visibility weighting ³³ . This method increases the relative weight of the visibilities at the longer baselines, which minimizes the beam size at the expense of signal-to-noise ratio. The superuniform beam characteristics and continuum root mean square (r.m.s.) noise are also given in Extended Data Table 1 . The increase of r.m.s. noise compared to more regular Briggs weighting (with a robust parameter of 0.5) depends on the telescope distribution and was a factor of about 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 about 2% in the final epoch to about 15% in the third epoch. The superuniform-weighted images of the final three epochs are presented in Fig. 1 and those of the first two epochs are shown 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 disk and the spatial PSD analysis were carried out on the calibrated visibilities. Although we focus our analysis on the higher-resolution band 7 observations, we also include the continuum result from the band 6 observations in the first section of the Methods . The observational details for these observations are also included in Extended Data Table 1 , and the calibration, self-calibration and imaging steps were identical to those carried out for the band 7 observations. The four spectral windows are centred at 218.9, 220.8, 230.0 and 232.9 GHz, and the increase in the continuum r.m.s. noise between Briggs weighting and superuniform weighting is a factor of 1.5.

Spatial PSD analysis

The spatial PSD is regularly used to derive information about, for example, the turbulent structure of the interstellar medium (for example, refs. ^{34 , 35} ) as well as the convective structure of the solar photosphere (for example, ref. ³⁶ ). The spatial PSD is given by the two-dimensional Fourier transform of an image. However, 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 ³⁴ . We can thus calculate the PSD directly from our interferometric visibilities without introducing potential artefacts during the imaging process. As the PSD would be dominated by the power at the scales of the stellar disk, we first use the uv -fitting code uvmultifit ³⁷ to fit the stellar disk in uv -coordinates (for a discussion on the disk profile, see the first section of the Methods ). Subsequently, we subtract the disk from the visibilities, after which, we calculate, for a phase centre towards the centre 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 as to avoid a contribution to the off-source PSD from the source signal. As the first two epochs have worse spatial resolution, we produce the PSD only 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 disk profile

To check how well a top-hat-shaped stellar disk profile fits the observations, we have also investigated the stellar disk 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 disk model convolved with the interferometric beam. The results of this comparison are shown in Extended Data Fig. 2 . The stellar disk model can accurately describe the observations with residuals at a level of about 2% of the peak emission. This means that the 338 GHz optical depth τ _338GHz increases steeply over only a small change in the radius. The radial motions we observe thus reflect the physical motion of the 338-GHz optical depth surface. As the optical depth is a strong function of the density ¹¹ , 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 Figs. 3 and 4 , respectively. Fitting the interferometric visibilities yields a completely circular stellar disk with F _225GHz = 221.4 ± 0.1 mJy and D _225GHz = 61.8 ± 0.1 mas. This means a brightness temperature T _b = 2,006 ± 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. 4 , it is clear that there is a very good correspondence. As 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.

1. Lecoanet, D. & Edelmann, P. V. F. Multidimensional simulations of core convection. Galaxies 11 , 89 (2023).

   Article ADS Google Scholar

2. Kupka, F. & Muthsam, H. J. Modelling of stellar convection. Living Rev. Comput. Astrophys. 3 , 1 (2017).

   Article ADS PubMed PubMed Central Google Scholar

3. Nordlund, Å., Stein, R. F. & Asplund, M. Solar surface convection. Living Rev. Sol. Phys. 6 , 2 (2009).

   Article ADS PubMed PubMed Central Google Scholar

4. Schwarzschild, M. On the scale of photospheric convection in red giants and supergiants. Astrophys. J. 195 , 137–144 (1975).

   Article ADS CAS Google Scholar

5. Ravi, V. et al. The non-uniform, dynamic atmosphere of Betelgeuse observed at mid-infrared wavelengths. Astrophys. J. 740 , 24 (2011).

   Article ADS Google Scholar

6. Montargès, M. et al. The convective surface of the red supergiant Antares. VLTI/PIONIER interferometry in the near infrared. Astron. Astrophys. 605 , A108 (2017).

   Article Google Scholar

7. Climent, J. B. et al. VLTI-PIONIER imaging of the red supergiant V602 Carinae. Astron. Astrophys. 635 , A160 (2020).

   Article Google Scholar

8. Norris, R. P. et al. Long term evolution of surface features on the red supergiant AZ Cyg. Astrophys. J. 919 , 124 (2021).

   Article ADS CAS Google Scholar

9. Wittkowski, M. et al. Near-infrared spectro-interferometry of Mira variables and comparisons to 1D dynamic model atmospheres and 3D convection simulations. Astron. Astrophys. 587 , A12 (2016).

   Article Google Scholar

10. Paladini, C. et al. Large granulation cells on the surface of the giant star π ¹ Gruis. Nature 553 , 310–312 (2018).

    Article ADS CAS PubMed Google Scholar

11. Vlemmings, W. H. T., Khouri, T. & Olofsson, H. Resolving the extended stellar atmospheres of asymptotic giant branch stars at (sub)millimetre wavelengths. Astron. Astrophys. 626 , A81 (2019).

    Article ADS CAS Google Scholar

12. Khouri, T. et al. An empirical view of the extended atmosphere and inner envelope of the AGB star R Doradus I. Physical model based on CO lines. Astron. Astrophys. 685, A11 (2024).

13. Velilla-Prieto, L. et al. Atmospheric molecular blobs shape up circumstellar envelopes of AGB stars. Nature 617 , 696–700 (2023).

    Article ADS CAS PubMed Google Scholar

14. Freytag, B. & Höfner, S. Global 3D radiation-hydrodynamical models of AGB stars with dust-driven winds. Astron. Astrophys. 669 , A155 (2023).

    Article ADS CAS Google Scholar

15. Reid, M. J. & Menten, K. M. Radio photospheres of long-period variable stars. Astrophys. J. 476 , 327–346 (1997).

    Article ADS CAS Google Scholar

16. van Leeuwen, F. Validation of the new Hipparcos reduction. Astron. Astrophys. 474 , 653–664 (2007).

    Article ADS Google Scholar

17. Rast, M. P. The scales of granulation, mesogranulation, and supergranulation. Astrophys. J. 597 , 1200–1210 (2003).

    Article ADS CAS Google Scholar

18. Rincon, F. & Rieutord, M. The Sun’s supergranulation. Living Rev. Sol. Phys. 15 , 6 (2018).

    Article ADS Google Scholar

19. Tremblay, P. E., Ludwig, H. G., Freytag, B., Steffen, M. & Caffau, E. Granulation properties of giants, dwarfs, and white dwarfs from the CIFIST 3D model atmosphere grid. Astron. Astrophys. 557 , A7 (2013).

    Article Google Scholar

20. Ma, J.-Z. et al. Is Betelgeuse really rotating? Synthetic ALMA observations of large-scale convection in 3D simulations of red supergiants. Astrophys. J. Lett. 962 , L36 (2024).

    Article ADS Google Scholar

21. Freytag, B., Liljegren, S. & Höfner, S. Global 3D radiation-hydrodynamics models of AGB stars. Effects of convection and radial pulsations on atmospheric structures. Astron. Astrophys. 600 , A137 (2017).

    Article ADS Google Scholar

22. Cherchneff, I. A chemical study of the inner winds of asymptotic giant branch stars. Astron. Astrophys. 456 , 1001–1012 (2006).

    Article ADS CAS Google Scholar

23. Vlemmings, W. H. T. et al. Resolving the stellar activity of the Mira AB binary with ALMA. Astron. Astrophys. 577 , L4 (2015).

    Article ADS Google Scholar

24. Vlemmings, W. et al. The shock-heated atmosphere of an asymptotic giant branch star resolved by ALMA. Nat. Astron. 1 , 848–853 (2017).

    Article ADS Google Scholar

25. Freytag, B., Holweger, H., Steffen, M. & Ludwig, H. G. On the scale of photospheric convection. In Science with the VLT Interferometer: ESO Astrophysics Symposia (ed. Paresce, F.) 316–317 (Springer, Berlin, Heidelberg, 1997).

26. Chiavassa, A., Plez, B., Josselin, E. & Freytag, B. Radiative hydrodynamics simulations of red supergiant stars. I. Interpretation of interferometric observations. Astron. Astrophys. 506 , 1351–1365 (2009).

    Article ADS Google Scholar

27. Lu, S. et al. Strong asymptotic giant branch stars’ spectral features in distant quiescent galaxies: impact on galaxy evolution. Preprint at https://arxiv.org/abs/2403.07414 (2024).

28. Ohnaka, K., Weigelt, G. & Hofmann, K.-H. Infrared interferometric three-dimensional diagnosis of the atmospheric dynamics of the AGB star R Dor with VLTI/AMBER. Astrophys. J. 883 , 89 (2019).

    Article ADS CAS Google Scholar

29. Trampedach, R., Asplund, M., Collet, R., Nordlund, Å. & Stein, R. F. A grid of three-dimensional stellar atmosphere models of solar metallicity. I. General properties, granulation, and atmospheric expansion. Astrophys. J. 769 , 18 (2013).

    Article ADS Google Scholar

30. Ramstedt, S. & Olofsson, H. The ¹² CO/ ¹³ CO ratio in AGB stars of different chemical type. Connection to the ¹² C/ ¹³ C ratio and the evolution along the AGB. Astron. Astrophys. 566 , 145 (2014).

    Article ADS Google Scholar

31. Vlemmings, W. H. T. et al. Rotation of the asymptotic giant branch star R Doradus. Astron. Astrophys. 613 , L4 (2018).

    Article ADS Google Scholar

32. Hunter, T. R. et al. The ALMA interferometric pipeline heuristics. Publ. Astron. Soc. Pac. 135 , 074501 (2023).

    Article ADS Google Scholar

33. CASA Team. et al. CASA, the Common Astronomy Software Applications for radio astronomy. Publ. Astron. Soc. Pac. 134 , 114501 (2022).

    Article ADS Google Scholar

34. Crovisier, J. & Dickey, J. M. The spatial power spectrum of galactic neutral hydrogen from observations of the 21-cm emission line. Astron. Astrophys. 122 , 282–296 (1983).

    ADS CAS Google Scholar

35. Green, D. A. A power spectrum analysis of the angular scale of Galactic neutral hydrogen emission towards L=140 deg, B=0 deg. Mon. Not. R. Astron. Soc. 262 , 327–342 (1993).

    Article ADS CAS Google Scholar

36. Wedemeyer-Böhm, S. & Rouppe van der Voort, L. On the continuum intensity distribution of the solar photosphere. Astron. Astrophys. 503 , 225–239 (2009).

    Article ADS Google Scholar

37. Martí-Vidal, I., Vlemmings, W. H. T., Muller, S. & Casey, S. UVMULTIFIT: a versatile tool for fitting astronomical radio interferometric data. Astron. Astrophys. 563 , 136 (2014).

    Article ADS Google Scholar

Download references

Authors and Affiliations

1. Department of Space, Earth and Environment, Chalmers University of Technology, Gothenburg, Sweden

   Wouter Vlemmings, Theo Khouri, Behzad Bojnordi Arbab, Elvire De Beck & Matthias Maercker

Contributions

W.V. wrote the paper and carried out the analysis of the data. T.K. wrote the observing proposal. W.V., T.K., B.B., E.D.B. and M.M. contributed to the interpretation of the results.

Corresponding author

Correspondence to Wouter Vlemmings .

Competing interests

The authors declare no competing interests.

Peer review information

Nature thanks Claudia Paladini and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/ .

Reprints and permissions