Measuring the Characteristic Sizes of Convection Structures in AGB Stars with Fourier Decomposition Analyses: the Stellar Intensity Analyzer (SIA) Pipeline.

(1)

Measuring the Characteristic Sizes of Convection Structures in AGB Stars with Fourier Decomposition Analyses.

Author: Miquel Colom i Bernadich,^∗ and Supervisor: Bernd Freytag^† Abstract:

Context. Theoretical studies predict that the length scale of convection in stellar atmospheres is proportional to the pressure scale height, which implies that giant and supergiant stars should have convection granules of sizes comparable to their radii. Numerical simulations and the observation of anisotropies on stellar discs agree well with this prediction.

Aims. To measure the characteristic sizes of convection structures of models simulated with the CO5BOLD code, to look at how they vary between models and to study their limitations due to numerical resolution.

Methods. Fourier analyses are performed to frames from the models to achieve spatial spectral power distributions which are averaged over time. The position of the main peak and the average value of the wavevector are taken as indicators of these sizes. The general shape of the intensity map of the disc in the frame is fitted and subtracted so that it does not contaminate the Fourier analysis.

Results. A general relationship of the convection granule size being more or less ten times larger than the pressure length scale is found. The expected wavevector value of the time-averaged spectral power distributions is higher than the position of the main peak. Loose increasing trends with the characteristic sizes by the pressure scale height increasing against stellar mass, radius, luminosity, temperature and gravity are found, while a decreasing trends are found with the radius and model resolution. Bad resolution subtracts signals on the slope at the side of the main peak towards larger wavevector values and in extreme cases it creates spurious signal towards the end of the spectrum due to artifacts appearing on the frames.

Conclusions. The wavevector position of the absolute maximum in the time-averaged spectral power distribution is the best measure of the most prominent sizes in the stellar surfaces. The proportionality constant between granule size and pressure length scale is of the same order of magnitude as the one in the literature, however, models present sizes larger than the ones expected, likely because the of prominent features do not correspond to convection granules but to larger features hovering above them. Further studies on models with higher resolution will help in drawing more conclusive results.

I. INTRODUCTION

Low mass stars (M < 8 M⊙, M⊙= 2 · 10³⁰ kg) evolve into the Asymptotic Giant Branch (AGB) at the end of their lifetime (τ ∼ 10⁹− 10¹⁰years). At this stage, their radius becomes several hundred times larger than their original one and their surface temperature drops to T ≈ 3000 K. Stars with M > 8 M⊙ go through a similar process, albeit in a much shorter time-span (τ ∼ 10⁶− 10⁸years), with larger radii and stronger gravity on their surfaces, receiving the name of supergiants. [1]

The opaque conditions at their interior make of convection the most efficient form of energy transport. Stel- lar convection was for a time a speculative topic, as it was impossible to observe it on any star besides our Sun. Precisely from observing the phenomena in our Sun, it was proposed in [2] that the convection length scale is proportional to the pressure scale height HP in stellar atmospheres. Therefore, unlike in Sun-like stars (R ∼ R⊙= 7 · 10⁵ km), where convection granules have a typical size of thousands of kilometers (σ ∼ 10³ km),

∗

it was predicted that in the low gravity conditions of the surface of giants and supergiants (R ∼ 10²− 10³ R_⊙, log(g [cm/s]) < −0.1) convection takes over a significant fraction of their radius, so that their surfaces would dis- play only around a dozen of very large convection cells.

This claim is also supported by 2D simulations [3]. Addi- tionally, 3D simulations also predict an irregular clump- ing of dust and the formation of shock waves in the stellar envelopes with similar or larger size scales [4].

Within available observations, it is seen that some AGB stars present time-dependant anisotropies in their disc of sizes comparable to their radius [5–7]. These observations often include variability in wavelength and po- larization, which is interpreted as inhomogeneous dust formation [8, 9]. Of special relevance are the observations of π¹ Gruis, a M7 III AGB [10], 1.5 M⊙ [11] star at a distance of d=163 pc [12]. Its stellar disc of angu- lar diameter 10 mas has been imaged using H band observations from the Very Large Telescope Interferometer with a resolution of 2 mas, revealing a patchy landscape with a characteristic structure size of 5.3 ±0.5 mas, measured with the use of spatial spectral power distributions (SPDs) [13].

The aim of this work is to repeat the same experiment,

(2)

Measuring the Characteristic Sizes of Convection Structures in AGB Stars with Fourier

Decomposition Analyses. Miquel Colom i Bernadich

[4, 14], to test the validity of the relation σ ≈ 10HP

found in [3], and to investigate possible shortcomings of the simulations. In the following section, the method for measuring the characteristic sizes of features on the stellar discs is explained. In the third section, the results are presented, and in the fourth one the limitations of the method and possible remedies are discussed.

II. METHODOLOGY

For the analysis of the structure sizes on the surface of stellar models, time-sequences of 2D frames showing one side of the stars are used. These frames are intensity maps It(~n), where ~n = (nx, ny) are the spatial pixel coor- dinates and t represents the time index. Spatial SPDs are used to determine the characteristic sizes of structures in the stellar surfaces, a technique also used in other works such as [13, 15]. The reader is refereed to the appendix for an extensive and detailed description of the method and the used scripts, as in the following subsections a summary of the implementation is shown instead.

A. The Sample Material.

The analyzed frames are retrieved from radiation- hydrodynamical simulations performed with the CO5BOLD code [14]. This code is specialized in repro- ducing energy transport inside of stars, be it radiative or convective, in a realistic manner. Its main advantage is that it does not need the insertion a free mixing-length parameter to characterize convection, unlike in the more traditional 1D stellar-structure calculations where parameter tuning is usually required. At the same time, the 3D aspect also allows for the characterization of structures in stellar surfaces, which is fundamental for this work.

In particular, this work focuses on a total on 39 simulations of giant and supergiant stars, all of them listed in Table I, some of which were already presented in [4], spanning a range of masses M = 1, 1.5, 5, 12 M⊙, radii R ∼ 400 − 800 R⊙, surface gravities log(g) ∼

−(0.3−1, 2) [cm/s] and temperatures T = 2500−3600 K.

Two other important parameters form the models are the pressure scale height, HP, and the amount of pixels needed to resolve it, NP, refereed from now on as model resolution. All these parameters are specified in Table I for every model.

B. Image Processing.

The It(~n) frames contain a stellar disc at their center, and as such they have a very powerful signal seen as a very strong peak in the SPD around the radius of the star, which outshines the signal corresponding to the smaller features. For this reason, it is compulsory to find a good

fit for the general profile of the star and subtract it, so only the signal from the features remains.

The first step to do this is to find a mask map Mt(~n) with the approximate shape of the disc, which is in general an irregular circle. Mt(~n) is a binary function, taking the value of 0 when ~n is outside of the mask and 1 when

~n is found inside. It is calculated by establishing an intensity threshold to a smoothed version of It(~n). This smoothing is applied because otherwise dark features on the disc would be left outside of the mask. If chosen cor- rectly, the threshold should also leave out the tails of the disc, which correspond to the limb and the inner atmo- sphere, usually presenting interesting features but that are not to be analyzed in this work.

Then, the interior of the disc is fitted with a limb dark- ening function centered in the frame, Bt(~n). This processing can be expressed as

I_t^′(~n) = Mt(~n) × (It(~n) − Bt(~n)) , (1) With this subtraction, most of the signal given by the star itself in the SPD is removed, allowing for the features on the surface to imprint a detectable signal. In a way, it is a form of imposing that the entire map has a total luminosity of 0, so that the only signal comes from variations instead of the total.

Before the processing is complete, a few more details need to be taken into account. Because of the 3D nature of the star, features close to the edge of the disc suffer from distortion and reddening. On top of that, the cut applied by Mt(~n) tends to leave a very sharp jump to 0 at the edge of the disc, which could potentially lead to artificial signal in the SPD distribution. Because of this, instead of having a hard jump from 0 to 1, the edges of Mt(~n) are set to 0 and then to gradually increase inwards until it reaches the value of 1. This way, features close to the edge of the disc are smeared out and a smooth transition from the outside to the inside of the mask is ensured. In Figs. 1a and 1b, the original It(~n) map and the processed map for the frame tframe= 700 of the model st28gm05n024 are shown as an example. For more details on how the processing is actually implemented in the code, the reader is refereed to the appendix. From now on, the expression It(~n) will represent to only processed maps.

C. Reading the SPD Distribution.

Given a sequence of frames with intensity maps It(~n), a Fourier transform is applied to it to get its 2D spectrum Iˆt(~k), where ~k is the wavevector.

Iˆt(~k) = F T [It(~n)], (2) Then, the power of ˆIt(~k) is integrated over rings of constant k = |~k| to get the SPD,

ft(k)dk = Z π/2

0

| ˆIt(k, ϕ)|²kdϕ

!

dk, (3)

(3)

st28gm05n024, t_frame=700

0 100 200 300

x(px,100px=400RO •)

y(px)

(a) Original image of the frame tf rame= 700 from the model st28gm05n024.

t_frame=700, processed

0 100 200 300

x(px,100px=400RO •)

y(px)

(b) The frame of FIG. 1a after going through the processing.

Notice that after subtraction, many pixels of the It(~n) map actually do have negative values, since they are darker than the pixels with 0 value outside of the star. This is to be expected after the subtraction of the solar limb fitting, and necessary if the total luminosity of the frame is to be set to 0.

FIG. 1

being ft(k) the SPD distribution in k = |~k|, an essential quantity for this work. In Fig. 2, the SPD distribution resulting from running the Fourier analysis on the processed image of the tframe = 700 from the model st28gm05n024 is shown. The most powerful peaks correspond to the most prominent structure sizes on the image in Fig. 1b.

From this sequence of distributions, many wavevector magnitudes can be measures, out of which we choose two.

The first one is the position of the peak in the time- averaged spectra.

kmax= max(f (k))k= max (hft(k)i)_k, (4) where max(f (k))k is the wavevector value k that holds

st28gm05n024, t_frame=700

0.0 0.2 0.4 0.6 0.8

k(px^-1) [f700(k)]1/2 (px1/2 )

FIG. 2: Resulting SPD from analyzing the frame tf rame= 700 of the model st28gm05n024 up to where it

has interesting features. Beyond k = 0.8px⁻¹ it just falls asymptotically to 0.

the absolute maximum of the distribution f (k), and the brackets indicate time-averaging. The second relevant measurement is defined as the mean value of the wavevector.

kmean=

* R∞

0 kft(k)dk R∞

0 ft(k)dk , +

(5)

Where the time average is performed as it is done in (4).

As it will be shown later, the distribution f (k) is usually smooth and contains a single peak. Finally, for every model, the characteristic granule sizes are computed as:

2π 2π

(4)

III. RESULTS

In Table 1 the results for σmax and σmean in every model are shown. The mean values are:

σmax

HP

= 25 ± 8, σmean

HP

= 8 ± 1, (7) If the 12 M⊙ stars are excluded then they become:

σmax

HP

= 24 ± 6, σmean

HP

= 8 ± 1, (8) Several plots are shown in Figs. 3, 4, 5 and 6, where σmax/HP and σmean/HP are plotted against stellar mass, radius, surface gravity, stellar luminosity, HP/R and model resolution, each of them with two linear regressions, one taking into account the 12 M⊙ stars (blue) and the other disregarding them (black). The plots show an increasing trend of σmax/HP with mass, radius, gravity and luminosity, and a decreasing one with HP/R and model resolution. The trends in σmean/HP are similar but with larger uncertainties on the slopes.

More subtle results are only seen by studying the resulting power spectra. For instance, in Figs. 7 and 8 it is made apparent that lower resolution models suffer a spurious increase of the signal at the highest value of the wave vectors towards the end of the spectrum, mainly due to artificially bright pixels and other sharp artifacts appearing on the analyzed frames. On the flip-side, models with higher NP do not showcase this artificial behaviour and showcase a finer structure with more power and features on the slope after the peak, or even a displacement of the peak towards smaller sizes. As an extra example, st28gm05n028 and st28gm05n029 strive to model the same star, with st28gm05n028 having larger results for σmean and σmaxand a resolution of NP = 1.8 px, while st28gm05n029 has been computed with NP = 2.5 px.

Looking at Fig. 11, it can be seen that st28gm05n028 has the higher resolution..

IV. INTERPRETATION AND DISCUSSION.

A striking result form the previous section is the great difference between the average values of σmax/HP and σmean/HP. The explanation is that σmax and σmean

are measurements of very different nature. σmean corresponds to the mean value of the time averaged SPD, and therefore it a less robust measurement, as it depends on the resolution of the SPDs, the size of the frames, and is influenced by artificial signal. σmax, by contrast, tells about which wavevector has the strongest signal, and therefore it is a meaningful physical measurement of the most predominant sizes appearing on the frames.

Keeping the previous statements in mind, it is then natural to ask why does σmean/HP come closer to 10 as in [3] than σmax/HP, if the latter is the most meaningful measurement. The issue lies in that it is not always certain that σmax is actually measuring the size

of convection granules. For instance, in the models st35gm03n19, st35gm03n18, st35gm03n17, st35gm03n07 and st35gm03n11, and even in st35gm04n34 and st35gm04n36, presented in Figs 7, 8 and 9, it is very easy to visually check that the sizes of σmax presented in Table I more or less correspond to the size of convection granules seen in the frames. But if one repeats the same experiment for st28gm06n040, st28gm06n041, and st28gm06n042 (Fig. 10), one will see that σmaxcor- responds instead more or less to the size of clouds and shock fronts hovering above the stellar surface, while the one from st28gm05n024 (Fig. 1) corresponds to the large bright and dark contrasts produced by hovering matter.

It is also seen in Figs 11, 12 and 13, which cover most of the model types, that the granulation is not strictly visi- ble, and what σmaxmeasures are most likely larger structures and shocks hovering above the surface. Therefore it makes sense that in all of these models σmax/HP > 10.

Therefore, is is difficult to interpret the behaviour of σmax/HP and σmax/HP against stellar mass, radius, luminosity, temperature, the fraction HP/R and model resolution. A great difference between the two is that many of the trends in σmean/HP are disrupted by the outlying signal given by the 12 M⊙ stars, which at the same time are the ones who present the most exaggerated values of σmax/HP. They are also the only ones that showcase their bare convection granules. It could then be surprising that their σmean/HP deviates so much from 10 too, but this situation is understandable because they are the oldest models and the most impaired by insuf- ficient model resolution, which suppresses fine structure and creates artificial signal towards the end of the SPD distributions, decreasing the measurement of σmean and increasing the one of σmax. It should come as no sur- prise that the highest resolution models from the 12 M⊙

stars are also the ones with the most moderate values for σmax amongst their kind, putting even more weight to the argument of resolution.

Taking these problematic models out of the equation, the trends of σmax/HP and σmean/HP are actually very similar: increasing in mass, radius, surface gravity, temperature and luminosity but decreasing in the fraction HP/R and the resolution. It is hard to see just from this which abscissa axis is the dominant, as mass, gravity, temperature and luminosity increase together, while HP/R is always lower the higher is the mass, forcing larger models to have a lower resolution. One possibility is that the trends are physical and that they speak about phenomena related to the physics of the structures that hover above the surfaces instead of the granules, as it has been said before. Another one is that the trends stem from the improvement of resolution, which tends to bring σmax/HP closer and closer to 10 (in this case, the decrease of σmean/HP with resolution should not con- tradict what has been said in the previous paragraph because the artificial signal at the end of the spectrum is not present anymore). If that is the case, then simulations should strive for a value of at least (NP >

∼ 6) px

(5)

if the regression in Fig. 4c is extrapolated up to fit with σmax/HP = 10. The trends also present large dispersion and have been made with a relatively small sample, mak- ing the possibility of them being spurious also worthy of consideration.

As a closing remark, within the available models, no reliable conclusion can be taken from the behavior of models using non-grey absorption coefficients of light and the grey ones.

V. CONCLUSION

The most reliable measure of the characteristic sizes performed in this work comes from finding which part of the time-averaged spectral power distribution that has the most signal. The average of this quantity divided by the pressure scale height σmax/HP amongst all the models results in a surprisingly high value of σmax/HP =

24 ± 6 in contrast to the value of σ/HP ≃ 10 found [3]. The most likely explanation is that this method is measuring structures hovering above the stellar surface instead of raw convection granules.

It has also been seen that bad resolution erodes and subtracts signal in the slope after the peak because the physics are not well resolved, and in the worst cases it creates artificial signal towards the end of the spectra. These problems are specially prominent in the oldest models of 12 M⊙ It has also been postulated that this might be one of the causes of the observed trends in the behaviour σmax/HP, albeit physical reasons or even the trends being spurious because of an unlucky arrangement in the sample are a possibility too.

In any of the cases, it is remarkable that both σmax/HP = 24 ± 6 and σmean/HP = 8 ± 1 are still of the same order of magnitude as σ/HP ≃ 10, indicating that the pressure scale can indeed be considered a good indicator for the convection length scale and granule sizes on stellar surfaces even in giant and supergiant stars as proposed in [2].

(6)

(a) Plot of the σmax/HP against stellar mass. Symbol legend explained in the larger figure caption.

(b) Plot of the σmax/HP against stellar radius. Symbol legend explained in the larger figure caption.

(c) Plot of the σmax/HP against surface gravity. Symbol legend explained in the larger figure caption.

(d) Plot of the σmax/HP against surface temperature. Symbol legend explained in the larger figure caption.

FIG. 3: Four plots of σmax/HP against stellar parameters such as mass, radius, surface gravity and temperature.

Color/temperature: Red: < 2300 K. Orange: 2300 - 2800 K. Yellow: 2800 - 3200 K. Cyan: 3200 - 4000 K.

Number of spikes/luminosity: 3: 4900 - 5100 L⊙. 4: 6600 - 7600 L⊙. 5: 9000 - 11000 L⊙. 6: 41000 - 43000 L⊙. 7: > 4300 L⊙.

Size/mass: Tiny: 1 M⊙. Small: 1.5 M⊙. Large: 5 M⊙. Very large: 12 M⊙.

(7)

(a) Plot of the σmax/HP against the fraction HP/R. Symbol legend explained in the larger figure caption.

(b) Plot of the σmax/HP against stellar luminosity. Symbol legend explained in the larger figure caption.

(c) Plot of the σmax/HP against simulation resolution, defined in this report as the number of pixels resolving the pressure

scale height. Symbol legend explained in the larger figure caption.

FIG. 4: Three plots of σmax/HP against the stellar parameters of luminosity, and others such as the fraction HP/R and simulation resolution.

(8)

Model M L T log g R HP NP σmax σmean Specifications

(M⊙) (L⊙) (K) [cm/s] (R⊙) (R⊙) (px) (R⊙) (R⊙)

st26gm07n001 1. 6993 2635 -0.77 401 15 3.0 251 111 Dust-free AGB star model.

st28gm05n016 1.5 7517 2952 -0.43 332 8 2.0 178 69 Grey AGB star model.

st28gm05n022 1. 5068 2660 -0.61 335 11 2.8 274 88 Grey dust-free AGB star model.

st28gm05n023 1.5 6997 2955 -0.40 319 8 1.9 217 64 Grey dust-free AGB star model.

st28gm05n024 1.5 6685 2898 -0.41 324 8 1.9 318 74 Non-grey dust-free AGB star model.

st28gm05n030 1.5 7011 2958 -0.40 319 8 1.9 320 61 Non-grey dust-free AGB star model.

st28gm06n038 1. 7090 2417 -0.92 480 20 4.0 435 164 Dusty large non-grey AGB star model.

st28gm06n039 1. 7068 2690 -0.74 387 14 2.9 317 108 Dusty large grey AGB star model.

st28gm06n040 1. 7009 2342 -0.97 509 21 4.3 470 180 Dusty large 2-bin non-grey AGB star model.

st28gm07n001 1. 10086 2506 -1.02 533 25 5.1 486 156 Grey dust-free AGB star model.

st28gm07n002 1. 10043 2251 -1.20 659 34 4.8 685 268 Grey AGB star model.

st28gm07n004 1. 10015 2250 -1.20 659 34 4.9 573 228 Grey AGB star model.

st35gm03n07 12. 92188 3485 -0.34 833 8 0.9 358 47 Low surface gravity Betelgeuse model.

st35gm04n34 5. 42055 3325 -0.45 618 10 1.9 222 103 Low surface gravity Supergiant model.

TABLE I: Table with the name of each model, the stellar parameters, the resolution, the measured characteristic sizes and some informative details.

(9)

(a) Plot of the σmean/HP against stellar mass. Symbol legend explained in the larger figure caption.

(b) Plot of the σmean/HP against stellar radius. Symbol legend explained in the larger figure caption.

(c) Plot of the σmean/HP against surface gravity. Symbol legend explained in the larger figure caption.

(d) Plot of the σmean/HP against surface temperature.

Symbol legend explained in the larger figure caption.

FIG. 5: Three plots of σmean/HP against the stellar parameters of luminosity, and others such as the fraction HP/R and simulation resolution. Asterisks on the uncertainties imply that they are of the order of unity, rendering the

regressions unreliable.

Size/mass: Tiny: 1 M . Small: 1.5 M . Large: 5 M . Very large: 12 M .

(10)

(a) Plot of the σmean/HP against the fraction HP/R. Symbol legend explained in the larger figure caption.

(b) Plot of the σmean/HP against stellar luminosity. Symbol legend explained in the larger figure caption.

(c) Plot of the σmean/HP against surface gravity. Symbol legend explained in the larger figure caption.

FIG. 6: Three plots of σmean/HP against the stellar parameters of luminosity, and others such as the fraction HP/R and simulation resolution. Asterisks on the uncertainties imply that they are of the order of unity, rendering the

regressions unreliable.

Size/mass: Tiny: 1 M⊙. Small: 1.5 M⊙. Large: 5 M⊙. Very large: 12 M⊙.

(11)

(a) The square root of the time-averaged SPD of models st35gm03n19 and st35gm03n18.

(b) Model st35gm03n19, square root of the intensity map.

σmax= 168 R⊙, σmean= 49 R⊙, NP = 1.3 px.

(c) Model st35gm03n18, square root of the intensity map.

σmax= 326 R⊙, σmean= 49 R⊙, NP = 0.9 px. The numerical artifacts are very obvious to plain sight.

FIG. 7: The square root of the time-averaged SPD of models st35gm03n19 and st35gm03n18 plus some frames showcasing their features. The comparison of these two in particular showcases how a slight increase in resolution

highly improves the finer structure and eliminates the very obvious numerical artifacts.

(12)

(a) The square root of the time-averaged SPD of models st35gm03n17, st35gm03n07 and st35gm03n11.

(d) Model st35gm03n11, square root of the intensity map.

FIG. 8: The square root of the time-averaged SPD of models st35gm03n17, st35gm03n11 and st35gm03n07 plus some frames showcasing their features. Similarly as in Fig. 7, this comparison showcases how a slight increase in

resolution highly improves the finer structure and eliminates the very obvious numerical artifacts.

(13)

FIG. 9: The square root of the time-averaged SPD of models st35gm03n17, st35gm04n34 and st35gm04n36 plus some frames showcasing their features.

(14)

(a) The square root of the time-averaged SPD of models st28gm06n042, st28gm06n041 and st28gm06n040.

(d) Model st28gm06n042, square root of the intensity map.

FIG. 10: The square root of the time-averaged SPD of models st28gm06n040, st28gm06n041 and st28gm06n042 plus some frames showcasing their features.

(15)

FIG. 11: The square root of the time-averaged SPD of models st28gm05n28, and st28gm05n29 plus some frames showcasing their features. This is a further example of how a slight increase in resolution highly improves the finer

structure.

(16)

FIG. 12: The square root of the time-averaged SPD of models st28gm05n001 and st28gm05n002 plus some frames showcasing their features.

(17)

FIG. 13: The square root of the time-averaged SPD of models st26gm06n001 and st26gm07n002 plus some frames showcasing their features.

(18)

[1] B. W. Carroll and D. A. Ostlie, An introduction to mod- ern astrophysics (Cambridge University Press, 2017).

[2] M. Schwarzschild, Astrophys. J 195, 137 (1975).

[3] B. Freytag, H. Holweger, M. Steffen, and H.-G. Ludwig, in Science with the VLT Interferometer (Springer, 1997) pp. 316–317.

[4] B. Freytag, S. Liljegren, and S. H¨ofner, A& A 600, A137 (2017).

[5] M. Karovska, P. Nisenson, C. Papaliolios, and R. Boyle, The Astrophysical Journal 374, L51 (1991).

[6] M. Karovska, W. Hack, J. Raymond, and E. Guinan, Astrophys. J Letters 482, L175 (1997).

[7] C. Haniff and D. Buscher, A&A 334, L5 (1998).

[8] K. Ohnaka, G. Weigelt, and K.-H. Hofmann, A&A 589, A91 (2016).

[9] M. Wittkowski, K.-H. Hofmann, S. H¨ofner,

J. Le Bouquin, W. Nowotny, C. Paladini, J. Young, J.-P.

Berger, M. Brunner, I. de Gregorio-Monsalvo, et al., A&A 601, A3 (2017).

[10] M. Feast, Mon. Not. R. Astron. Soc 113, 510 (1953).

[11] A. Mayer, A. Jorissen, C. Paladini, F. Kerschbaum, D. Pourbaix, C. Siopis, R. Ottensamer, M. Meˇcina, N. Cox, M. Groenewegen, et al., A&A 570, A113 (2014).

[12] F. Van Leeuwen, A&A 474, 653 (2007).

[13] C. Paladini, F. Baron, A. Jorissen, J.-B. Le Bouquin, B. Freytag, S. Van Eck, M. Wittkowski, J. Hron, A. Chi- avassa, J.-P. Berger, et al., Nature 553, 310 (2018).

[14] B. Freytag, M. Steffen, H.-G. Ludwig, S. Wedemeyer- B¨ohm, W. Schaffenberger, and O. Steiner, J. Comput.

Phys. 231, 919 (2012).

[15] R. Trampedach, M. Asplund, R. Collet, ˚A. Nordlund, and R. F. Stein, The Astrophysical Journal 769, 18 (2013).

(19)

APPENDIX: the Stellar Intensity Analyzer (SIA) Pipeline.

Abstract

The SIA pipeline takes a set of time-dependent pictures of stellar disks and uses a Fourier Analysis to measure the characteristic sizes of their features and other useful quantities, such as standard deviations or the spatial power distributions of features. The main core of the pipeline consists in identifying the stellar disc in the frames and subtracting their signal from the spatial power distributions through a general fit of the disc intensity.

To analyze a time sequence, the SIA pipeline requires at least two commands from the user. The first command orders the SIA pipeline to read the .sav IDL data structure file where the frame sequence is stored and to produce another .sav file with information on the spectral power distributions, the second command orders the reading of such file to produce two more .sav files, one containing time-averaged size measurements and their deviations while the other breaking down time-dependant information and other arrays used for the calculations.

The SIA pipeline has been entirely written in Interactive Data Language (IDL). Most of the procedures used here are original from the SIA pipeline, but a small handfull like ima3 distancetransform.pro, power2d1d.pro extremum.pro and smooth2d.pro from Bernd Freytag and peaks.pro and compile opt.pro amongst others are actually external.

Contents

1 Introduction 2

1.1 About this Appendix . . . . 2

1.2 General view of the data flow . . . . 2

1.3 The structure of this appendix . . . . 2

2 Data structures 4 2.1 model intens.sav . . . . 4

2.2 model sia data.sav . . . . 4

2.3 model sia statistics.sav . . . . 5

2.4 model sia vectors.sav . . . . 6

2.5 allbasDB.sav . . . . 7

2.6 allseqDB.sav . . . . 8

3 Main routines of the SIA pipeline. 9 3.1 sia subtract back area.pro . . . . 9

3.2 sia analyze timeseq2.pro . . . . 12

3.3 sia data collect.pro . . . . 13

3.4 sia data statistics.pro . . . . 16

4 Optional but major auxiliary routines. 19 4.1 sia models loop . . . . 19

4.2 sia example.pro . . . . 19

4.3 sia prepare postscripts.pro . . . . 20

5 Essential but minor auxiliary routines. 22 5.1 sia padding 0.pro . . . . 22

5.2 sia unpadding.pro . . . . 22

5.3 sia initialize grid.pro . . . . 22

5.4 sia moment calc.pro: . . . . 23

5.5 sia 2dfit gauss n.pro . . . . 23

5.6 sia fourier analyzer.pro . . . . 24

(20)

1 Introduction

1.1 About this Appendix

This document is intended to act both as a user manual and as a blueprint for the Stellar Intensity Analyzer (SIA) routines used to perform a Fourier Analysis on the frames of stars. As such, it should be the first and most important source for any SIA user. It includes: information on the input and output files, an exposition of the roles of every script, a detailed following on the steps done in every major original script and a description of the usage of many other auxiliary ones, including several captions and diagrams to make their usage and functionality more comprehensible.

1.2 General view of the data flow

The SIA pipeline takes a set of time-dependent frames containing a stellar disc of a single star at their center, processes them to leave only the features of the star on them, and then performs a Fourier Analysis to build their spatial spectral power distributions (SPDs). In the end, it measures time-averaged characteristic sizes on the stellar disc and other relevant quantities along with their typical deviations. The frames are to be contained inside of an IDL data structure file named model intens.sav, where the word model identifies the sequence.

To analyze a set of frames, the user is required to run a program two times, identified as the two steps in Fig. 1. The first run computes the necessary elements to build the SPDs out of the Fourier decomposition of the processed frames, which are stored inside of model sia data.sav, another IDL data structure. The second one reads the characteristic sizes from the SPDs, and it produces the files model sia statistics.sav and model sia vectors.sav. The former is another IDL data structure containing a list of measurements of typical spatial wavelengths and wavevectors and their typical and standard deviations. The latter involves a full expansion of the vectors used to calculate the time-averages.

It also contains the time-averaged SPD and smoothed versions of it used in the code. The time-step between two different frames is assumed constant. Since this is hard to achieve with real observations, the pipeline is actually more useful to analyze simulation outputs.

1.3 The structure of this appendix

This document has four main sections besides this introduction. In the following section or Section 2, the four IDL data structures that are a direct input or output of the SIA pipeline and an extra one used by an auxiliary routine are listed along with the meaning of their variables.

Section 3 consists in a listing of all the scripts that are absolutely essential for the processing of information across the SIA pipeline and are recommended to be at least overviewed if the user wants to have an idea on how information is read and transformed. These include: sia subtract back area.pro, sia analyze timeseq2.pro, sia data collect.pro and sia data statistics.pro.

Section 4 follows with the description of routines and functions that, despite not being essential for the SIA pipeline, are of recommended use by the user. These routines usually showcase some functionality of the SIA pipeline, prepare plots or act as a shortcut for calling some parts of the SIA pipeline.

And finally, Section 5 will includes a less exhaustive overview of some essential but minor routines that act as building blocks of the SIA pipeline by for example performing calculations or initializing arrays.

(21)

model intens.sav

Fourier Analysis Pipeline:

sia data collect, sia analyze timeseq2, sia subtract back area, power2d1d

model sia data.sav allbasDB.sav

sia data statistics

model sia statistics.sav model sia vectors.sav

Figure 1: Schematic layout of the entire SIA data flow and the routines acting at every stage.

model intens.sav is the generic name for the input file, which is supposed to contain a time sequence of stellar frames. model sia data.sav is created after sending the frames through the Fourier Analysis Pipeline, and model sia statistics.sav and model sia vectors.sav are intended to be the final output from which useful information can be read, such as time-averaged quantities and their deviations.

allbasDB.sav is an extra database containing some physical information that is required to be able to normalize the units of measurements in the frames.

sia data collect Loops over

sia analyze timeseq2

Runs over

1. sia subtract back area 2. power2d1d

Figure 2: Schematic layout of the routine hierarchy in the Fourier Analysis Section. sia data collect.sav takes an entire model or sequence of frames, while sia analyze timeseq2.pro takes on individual frames.

(22)

2 Data structures

In Figs 1 and 2 an outline of the data flow is shown. Fig. 1 in particular shows the different stages of data stored in structures, going from inputs to outputs. As it can be seen, the IDL structure files described in this section constitute the spinal cord of the data flow in the SIA pipeline, and understanding their variables is equivalent to understanding the shape of the input and output information. The structures are as follows:

2.1 model intens.sav

This is a generic name for the data structure containing the sequence of frames. It is the main input of the SIA pipeline, ad it contains a sequence of frames. The scripts assumes that at least the two following variables come with it:

• intens.i3r: a 3D array with the time dependant intensity maps Itp~nq or frames. The first two indices indicate the pixel position ~n“ pnx, nyq inside of the frame, while the third index is the time index.

• intens.model: a string with the model name.

2.2 model sia data.sav

This is the output of running the Fourier Analysis Pipeline on to a model intens.sav which, as showcased in Fig. 2, it basically means running the mother script of sia data collect.pro. It contains raw information of the SPD, but it has the flaw that all the data is normalized according to the size of the padded box (see the description of sia analyze timeseq2.pro in Section 3.2). For instance, if the variable n padd is set to the default value of 1089, then a wavevector of magnitude k implies a wavelength of λ“ 1088{k px. This happens because the Fourier Analysis Pipeline does not deal with any physical information of the analyzed stars, as it is carried by the allbasDB.sav file instead (see Section 2.5). The variables that it contains are all called sia data.variable, where .variable can take the names of:

• .model: a string with the model name.

• .TotalTimeSteps: an integer containing the total amount of frames.

• .TimeSteps: a vector running from 0 to sia data.T otalT imeSteps´ 1. Each component is the time index t of the corresponding frame, and every vector listed after this one has the same size.

• .TotalIntensity: a vector resulting from adding up the intensity value of every pixel ~n in every frame Itp~nq.

L_total,t“ÿ

~ n

I_tp~nq (1)

• .SecondMomentRad: a vector resulting from calculating the second-moment radius of the distribution in every frame.

r_2,t“ dř

~

n|~n ´ ~ncenter|Itp~nq ř

~

nI_tp~nq (2)

• .MeanMaskRad: a vector containing the mean radius of the disc mask (see the description of the procedure sia subtract back area.pro in section 3.1).

• .MaxMaskDist: a vector containing the largest distance of the edge of the mask from the central pixel ~ncenter.

• .MinMaskDist: a vector containing the shortest distance of the edge of the mask from the central pixel ~ncenter.

• .TotalArea: an integer containing the amount of pixels inside the mask.

• .FittingPreset: a string that tells which fitting preset has been used inside of the mask (see section 3.1).

• .PolyFitCoef: an integer that tells the order of the fitting inside of the mask.

• .WaveDelta: a number telling the wavevector spacing used for the Fourier Analysis (see the external routine power2d1d.pro).

• .n original: a number telling the original size of the frame.