© ESO 2018

1. Introduction

The Sun exhibits convective motions at its surface on a range of scales, reaching from granulation to giant cells. Supergranulation (Rieutord & Rincon 2010) is a pattern on intermediate scales (∼30 Mm) with horizontal velocities of a few hundred meters per second. Even half a century after the discovery of supergranulation (Hart 1954), its origin is still not understood, although many attempts have been made (see,e.g., Rieutord & Rincon 2010, for an extensive list). Recently Cossette & Rast (2016) suggested, based on results from anelastic simulations, that supergranulation might be the largest spatial scale of convection that is driven by buoyancy. Featherstone & Hindman (2016), on the other hand, argued that solar rotation might suppress convection on scales larger than supergranules. Other scenarios involve a prominent role of the magnetic field. For example, Crouch et al. (2007) and Thibault et al. (2012) consider that the network field could arise from random clustering of magnetic elements, without the need for supergranular flows. Crouch et al. (2007) then hypothesize that the network might drive the supergranular flows.

Another prominent feature of the supergranular pattern is its super-rotation, meaning that the pattern rotation rate is higher than that of granules and magnetic elements (e.g., Duvall 1980; Snodgrass & Ulrich 1990; Meunier & Roudier 2007). Using time-distance helioseismology and direct Doppler imaging, respectively, Gizon et al. (2003) and Schou (2003) found signatures in the supergranular power spectrum that are characteristic of waves with periods of 6–9 days. At the equator, these waves were found to travel predominantly in the east-west direction, with more power in the prograde direction. This excess westward power leads to the observed super-rotation of the pattern. In the current paper we present new and independent evidence in support of these findings.We note that the reality of the wavelike properties of supergranulation was questioned by Rast et al. (2004) and Hathaway et al. (2006).

To further characterize the supergranulation phenomena observed by Gizon et al. (2003), Schou (2003), here we have analyzed Dopplergrams and intensity images from the Helioseismic and Magnetic Imager (HMI, Schou et al. 2012) onboard the Solar Dynamics Observatory (SDO, Pesnell et al. 2012) using time-distance helioseismology (TD) and local correlation tracking (LCT), two independent techniques that are known to deliver consistent results for horizontal flow measurements near the solar surface (e.g., Langfellner et al. 2015b).We then related the flow evolution to the evolution of the magnetic field to study their interplay. To that end, we compared magnetograms measured by HMI with simple simulations, where corks represent the magnetic elements.

After introducing the observations, in this paper we show the signature of supergranulation in power spectra of the travel-time and LCT divergence maps on an extended wavenumber interval and study how this signature is influenced by the tracking rate. Then we present the horizontal flow evolution of the average supergranule over a range of spatial scales. In addition, we compare and discuss the results of the observed and simulated magnetic field. Finally, we draw conclusions for the dynamics of the supergranulation.

2. Observations and data processing

We used one year of SDO/HMI Dopplergrams and intensity maps in the period from May 2010 through April 2011. We tracked regions of size ∼180 × 180 Mm² at the equator for 11 days, following those regions from about 70° east to 70° west. For the tracking rate, we used two options: a rate consistent with the rotation rate at the equator from Snodgrass (1984), and a rate faster than the Snodgrass rate by 60ms⁻¹ (in the rest of the paper dubbed “fast tracking”). The fast rate corresponds to the supergranulation pattern rotation rate at the equator as measured by Meunier & Roudier (2007), but is slower than the values obtained by Duvall (1980) and Snodgrass & Ulrich (1990). We remapped the regions using Postel’s projection and a pixel size of 0:348Mm at a cadence of 45 s.

We then divided the 11-day tracked and mapped datacubes into 33 non-overlapping segments of eight hours each. Segment 17 thus crosses the central meridian. For each Dopplergram segment, we selected the f modes by applying the ridge filter described in Langfellner et al. (2015b).We then computed point-to-annulus wave travel times (Duvall et al. 1996) with an annulus radius of 10 Mm, using the method of Gizon & Birch (2004). These travel times, measured as the difference between the inward and outward travel times, are sensitive to the divergence of the horizontal flows near the solar surface. For comparison, we also computed the horizontal flow divergence from local correlation tracking (LCT) of granules (November & Simon 1988). We first obtained the horizontal flowcomponents using the code described in Langfellner et al. (2015b) and then computed the divergence using Savitzky-Golay filters (Savitzky & Golay 1964).

Following the procedures from our previous work (Langfellner et al. 2015b), we constructed the average supergranule. To minimize edge effects, we cut away the outer 10 Mm in the travel-time maps and then removed the map-averaged travel time. We also removed the noise on scales kR_⊙ > 300 by applying a low-pass filter in Fourier space. In these processed travel-time maps, we identified the positions of both supergranular outflow and inflow centers, only using locations that are at least 18 Mm away from the map edges. We then co-aligned both the processed travel-time maps and the LCT horizontal divergence maps at these positions, always using the coordinates obtained for the segment in which the supergranules cross the central meridian (segment 17). This will allow us to study the evolution of the flows in the average supergranule near the solar surface.

To study the magnetic field in the context of the supergranular flows, we also obtained HMI line-of-sight magnetograms for the same times and regions as the Dopplergrams and intensity images. We averaged the magnetograms over eight hours to reduce fluctuations on short timescales and took the absolute value to avoid cancellation of opposite polarities in the further processing. We then co-aligned these processed magnetograms at the supergranular positions.

3. Supergranular power spectrum

3.1. Dependence on tracking rate

Using time-distance helioseismology and direct Doppler imaging, Gizon et al. (2003) and Schou (2003) found that supergranulation exhibits a wave-like behavior. If studied at the equator, the wave power shows an excess power forwaves traveling in the west (prograde) direction. Rast et al. (2004) claimed that these observations do not necessarily imply the existence ofwaves but could rather be artifacts arising from a tracking rate that is too slow.

To check this, we studied the power spectra of the divergence signal for two different tracking rates, one that is 0.028930 μrad s⁻¹ (about 20 m s⁻¹) slower (the Snodgrass 1984, tracking rate) and the other 0.057139 μ rad s⁻¹ (about 40 m s⁻¹) faster (the “fast” tracking rate) than the Carrington rotation rate applied by Gizon et al. (2003). For each 11-day datacube, we collated the 33 travel time and LCT divergence maps in three-dimensional arrays and Fourier-transformed these arrays to obtain the power spectra.

Figure 1 shows various cuts through the three-dimensional power spectra for the travel-time cubes at the solar equator, using the two tracking rates and TD. The power is concentrated in a ring at supergranular scales (kR_⊙ ∼ 100), with most power in the prograde direction. In the east-west direction, cutting at k_x ≈ 120/R_⊙ and k_y = 0 (see Fig. 2, upper-left panel), two peaks are visible for the Snodgrass tracking rate. The stronger peak is located at about 3 μHz, the weaker one at ∼0 μHz. For the faster tracking rate, the power appears shifted to the lower frequencies, so that the two peaks have equal distance from zero frequency (∼ ± 2 μHz), but otherwise the power distribution is very similar to the Snodgrass tracking rate. For LCT (upper-right panel), the power distribution is almost indistinguishable from the time-distance case. The frequency shift of the peaks in the east-west direction can be explained by the Doppler shift Δω/2π = ku_x ≈ 1.7 μHz due to the tracking rate difference u_x ≈ 60 m s⁻¹. In the north-south direction (k_x = 0, k_y ≈ 120/R_⊙, lower panels), the two peaks are roughly equal in height and are located at ∼ ± 2 μHz, irrespective of the tracking rate. Again, TD and LCT agree well.

Fig. 1. Supergranular power spectrum of the divergence signal at the solar equator using time-distance helioseismology. Top row: at ky = 0. The left panel shows the spectrum for a tracking rate set to the Snodgrass (1984) rotation rate (about 20 m s−1 slower than Carrington). The right panel shows the spectrum for a tracking rate that is faster by about 60 m s−1. Bottom row: cut through k plane at frequency 3.2 μHz.

Fig. 2. Supergranular power spectrum of the divergence signal at the solar equator for both tracking rates, cut at kR⊙ ≈ 123, using one year of HMI data: time-distance helioseismology (left) and local correlation tracking right). Top row: East-west direction. Bottom row: North-south direction.

3.2. Fitting a parametric wave model

The power distribution and amplitude ratios of the peaks are in agreement with Gizon et al. (2003, cf. Fig. 1 therein). To extract wave properties like the dispersion relation, we fit the power, separately for each wavenumber k, using a parametric wave model that builds on the model presented in Gizon & Duvall (2004). The model consists of a sum of two Lorentzians, centered at frequencies ±ω₀, whose amplitudes depend on the azimuth ψ. The wavevector k = (k_x, k_y) is parameterized as k = (k cos ψ, k sin ψ). In addition, the model takes into account the Doppler shift k · u from a horizontal flow u = (u_x, u_y), a finite lifetime t_life = 1/γ, and a background term that is independent of frequency. The full model is given by

$$P(k,\omega ) =\frac{F(k,\omega )+F(-k,-\omega )}{2}+B_0(k)+B_2(k)\hspace{0.17em}\cos (2\psi -\alpha ),$$(1)

$$\begin{array}{ll}F(k,\omega )\hfill & =\frac{A_0(k)+A_1(k)\hspace{0.17em}\cos (\psi -{\psi }_{\max })+A_2(k)\hspace{0.17em}\cos (2\psi -\alpha )}{1+{[\omega -{\omega }_0(k)-k\cdot u(k)]}^2/\gamma {(k)}^2},\hfill \end{array}$$(2)

where ψ_max is the direction of maximum power, A₁/A₀ gives the power anisotropy, and the A₂ term describes spatial variations with respect to an angle α that can arise from, for example, instrumental astigmatism. The coefficients B₀ and B₂ describe the background term.

The observed and fitted power as a function of azimuth and frequency is shown for kR_⊙ ≈ 123 in Figs. A.1 and A.2. For the Snodgrass tracking rate, the bands of strong power are located at frequencies that depend on azimuth, consistent with a Doppler shift Δω = k · u = ku_xcosψ due to a flow u_x ≈ 60 m s⁻¹. Adjusting the tracking rate, as we did in the fast tracking case, removes this Doppler shift and thus the cosψ dependence of the frequencies where the power has a maximum, but leaves the power spectrum otherwise largely unaffected (see Fig. A.2).

Figure 3 shows the dispersion relation (with f₀ = ω₀/2π) that we obtained through the fits of the power spectra. The curves for TD and LCT agree well and do not show any dependence on the tracking rate. Moreover, the dispersion relation that we measured is consistent with the empirical power law estimate obtained by Gizon et al. (2003) and Schou (2003) for kR_⊙ < 150. At larger wavenumbers, our measured frequencies are systematically below the power law estimates. However, in the interval 150 < kR_⊙ < 200, the values obtained by Gizon et al. (2003) and Schou (2003) are not well described by their power laws either. Indeed, their frequency values are ∼2.0 μHz for these wavenumbers, which is consistent with our measurements. The regime at kR_⊙ > 200 was not explored previously. Our data show a continuation of the flat dispersion relation until the fits start to diverge at kR_⊙ ∼ 270.

Fig. 3. Parameters from Lorentzian fits of time-distance and LCT divergence power spectra for Snodgrass (1984) and fast tracking rates. Top row: dispersion relation and half-width at half maximum (HWHM). Bottom row: horizontal velocity components of background flow. The tracking rate difference of 60 m s−1 has been subtracted from the ux values for Snodgrass (1984) tracking to facilitate the comparison with fast tracking. Values for kR⊙ > 230 (Snodgrass tracking) and kR⊙ > 271 (fast tracking) are omitted, as the fits are not stable in that regime.

The half-width at half maximum (HWHM), γ/2π, of the Lorentzians is independent of the tracking rate as well. It is roughly constant for kR_⊙ < 100 and then starts to increase at an accelerating rate for larger wavenumbers. At kR_⊙ ≈ 200, the HWHM becomes larger than f₀. This behavior agrees well with Gizon et al. (2003) and Schou (2003), confirming their findings.

The bottom row of Fig. 3 shows that the flow speed u_x is not constant, but decreases with kR_⊙ for kR_⊙ < 240. This could be a sign that different scales of supergranules are affected differently by the radial differential rotation in the near-surface shear layer, and thus probe different depths. The north-south flow component u_y, on the other hand, is consistent with zero for all scales and both tracking rates, as expected from meridional circulation at the solar equator.

For the fast tracking, Fig. 4 shows some selected parameters derived from the dispersion relation and ridge width of the power spectra, together with values derived from the analysis of the average supergranule that is described in Sect. 4.2. The period P = 1/f₀ varies from ten days at kR_⊙ ≈ 40 to 6 days at smaller wavenumbers (kR_⊙ > 120). The lifetime, t_life = 1/γ, decreases with wavenumber, from about three days at the largest scale to half a day at kR_⊙ ≈ 250. The phase speed, c_ph = ω₀/k, drops quickly from about 150 m s⁻¹ at the largest scale to 60 m s⁻¹ at kR_⊙ = 120, and then decreases slowly to 30 m s⁻¹ at kR_⊙ = 250. The group speed c_gr = ∂ω₀/∂k, on the other hand, shows a peculiar decline from about 40 m s⁻¹ at kR_⊙ ≈ 40 to zero at kR_⊙ = 200. This coincides with a drop of the quality factor (f₀/HWHM) below one. The highest quality factor is reached at kR_⊙ ∼ 90. Finally, the e-folding length c_pht_life is monotonically decreasing with wavenumber. There is a conspicuous change of slope in the log-log plot at kR_⊙ ≈ 120, the scale of strongest power for supergranulation. At larger scales, the e-folding length increases more slowly than the wavelength. We might speculate that this is connected to a convective regime change at this scale, as proposed by, for example, Cossette & Rast (2016) and Rincon et al. (2017).

Fig. 4. Parameters extracted from Lorentzian fits of divergence power spectra and the average convection cell evolution (see Sect. 4.2) for fast tracking rate at the solar equator, as a function of spatial scale. Values for kR⊙ > 271 are omitted as the fits are not stable in that regime. For f0 the dispersion relations from Gizon et al. (2003) and Schou (2003) are overplotted for comparison.

3.3. Distinguishing wavenumber regimes

When considering the observations of the supergranular power spectrum at a slightly higher level of abstraction, the following wavenumber intervals can be distinguished:

1. kR_⊙ ≲ 120: The dispersion relation can be described by a squareroot function and there are two clearly distinguishable peaks in frequency cuts of the power spectrum (see, e.g., Fig. 2).

2. 120 ≲ kR_⊙ ≲ 200: The dispersion relation is flat (f₀ ≈ 2 μHz) and there is no more visible minimum between the two peaks in the east-west direction – this is mostly because the power in the east direction is much smaller than in the west direction.

3. 200 ≲ kR_⊙ ≲ 270: The quality factor has fallen below one (the HWHM of the ridge exceeds its central frequency f₀).

4. kR_⊙ ≳ 270: Even in the north-south direction, where the two peaks are of equal strength, there is no more visible minimum between the two peaks. Analogous to the Rayleigh criterion in optics, this provides an estimate for the resolution limit. This is supported by the fact that the Lorentzian fits are diverging or completely failing in this regime, so no meaningful dispersion relation is extractable beyond kR_⊙ ∼ 270.

4. Supergranular flow evolution

4.1. The average supergranule

We now analyze the evolution of the intermediate-scale convection in real space by measuring the near-surface flows of the average supergranule at a reference time and following it up to five days into the future and past.

Figure 5 shows the evolution of inward–outward travel times for the average supergranular outflows and inflows¹. The travel-time maps show a complex evolution with an oscillatory component. At the reference time, Δt = 0, the central outflow region is surrounded by a ring of inflows. The outflow region drifts westward for Δt > 0, while the eastern part of the inflow ring follows in the same direction and gains strength. After a short dipolar stage, where both components are roughly equally strong (Δt ∼ 1.7 d), the inflow region moves further west into the center. After about three days, the pattern seen at the reference time has reversed: a central inflow region is now surrounded by a ring of outflows. However, the size of the structure has increased; the ring diameter is roughly 50 Mm compared to about 30 Mm at the reference time. For larger time lags (Δt > 4 d), the signal mostly vanishes beneath the noise background. Yet it can be seen that at Δt = 5.3 d the sign at the origin reappears and switches back to the sign at reference time. Before the reference time (Δt < 0), a similar behavior is observed as for positive Δt, albeit with the time axis reversed and the east-west direction flipped. At times Δt ∼ −3 d and Δt ∼ +3 d, the patterns are very similar – the supergranulation pattern thus appears to undergo an oscillation with a period of approximately six days.

Fig. 5. Evolution of f-mode-filtered inward–outward travel times for the average supergranule. Top: outflow. Bottom: inflow. The frames, from top left to bottom right, span a total of eleven days in non-overlapping steps of eight hours. The central frame (red outline) is for the reference time, Δt = 0. We note that the color scale is different for each individual frame; saturation is reached at the maximum absolute value (given by the lower label boxes) of the respective frame. The color map is symmetrized around zero; outflows are red and inflows are blue.

For the average inflow at reference time, the evolution is very similar to the outflow but with the sign of the travel times reversed. Small differences to the average outflow can be noticed, which manifest as a generally smaller travel-time amplitude and a slightly shorter oscillation period.

In LCT divergence maps (Fig. 6, using the same supergranule positions as before), the evolution of the average supergranule shows the same characteristics as for the travel-time maps. Outflows change into inflows and vice versa in about three days; the ring diameter increases away from the reference time. The patterns in the LCT divergence maps are equal to the travel-time maps down to the smallest scales present in the (low-pass-filtered) inward–outward travel-time maps. There is, however, a slight asymmetry in time. For the average outflow, the divergence is stronger before the reference time, for the average inflow (not shown), it is stronger after the reference time. The reason for this temporal asymmetry is unclear; among the possible causes are the different depth sensitivities of f-mode TD and granule LCT or a systematic influence of the magnetic field on the divergence measurements.

Fig. 6. Evolution of LCT divh for the average supergranular outflow (cf. top of Fig. 5).

To summarize, the supergranular evolution at the equator is both symmetric in time and the east-west direction; inflows and outflows turn into each other and can essentially be regarded as phase-shifted versions of one another. In other words, new supergranules preferentially form in old inflows. This is consistent with the findings by Shine et al. (2000) who observed several individual supergranules forming in LCT divergence maps obtained from Michelson Doppler Imager (Scherrer et al. 1995) intensity images. The sign reversal of the horizontal divergence after about two days was also observed by Greer et al. (2016) using helioseismic ring-diagram analysis (see Figs. 7 and 8a therein). The authors measured a super-rotation velocity of the supergranular pattern of about 50 m s⁻¹ compared to the Carrington rotation rate, which is slightly faster than our fast tracking rate (40 m s⁻¹).

Fig. 7. Evolution for different scales: space-time images. Average over y band (10 Mm) is given for the average outflow. Top row: TD. Bottom row: LCT. The dashed lines give the phase speed from the dispersion relation as obtained from the power spectrum. Outflows are red and inflows are blue. To increase visibility, the color maps are saturated at half the maximum value.

Fig. 8. Evolution of supergranular cells on different scales: Average over circular area with radius (2πR⊙/k)/4 around the origin. Top row: TD outflow and inflow (sign reversed). Bottom row: LCT outflow and inflow.

The oscillatory behavior does not depend on the tracking rate. Figure A.5 in the Appendix shows the evolution of the average supergranular outflow for the Snodgrass (1984) tracking rate, using time-distance helioseismology. The visible patterns are almost indistinguishable from the 60 m s⁻¹ faster supergranule pattern rotation rate (Fig. 5), except for an east-west drift of the pattern in the Snodgrass case that is consistent with the difference in the tracking rates. This clearly contradicts the claims by Rast et al. (2004) and Hathaway et al. (2006) that the measured oscillation is an artifact due to tracking rate or projection effect issues.

In the remainder of this paper, we consider only the faster tracking rate. This will facilitate the analysis of the supergranular pattern evolution over the course of days, as the envelope maximum of the pattern does not drift in the east-west direction.

4.2. Dependence on spatial scale

In order to study the change in size of the supergranular evolution pattern and to make a connection to the findings from the power spectra, we investigated the evolution of the average cell separately for a set of narrow scale intervals, spanning the overall scale range from about 15 Mm to 80 Mm. To do so, we bandpass-filtered the travel-time maps in Fourier space to only select supergranular cells in a narrow scale range before applying the image segmentation algorithm to identify the positions of the cells. For each bandpass filter, we then co-aligned the TD and LCT maps (low-pass-filtered and not filtered, respectively) at these positions as before. The results for convective outflows at three selected scales are shown as space-time plots in Fig. 7.

The oscillatory behavior (outflows turning into inflows and vice versa) can be observed at all the scales we analyzed, with good agreement between TD and LCT. When comparing the scales by eye, two things stand out: The oscillation period decreases with higher wavenumber (roughly from 8 d to 4 d) and the lifetime of the pattern decreases (at kR_⊙ = 52, the pattern is visible in all panels, whereas at kR_⊙ = 222, it disappears by |Δt| ≈ 3 d).

These observations are in line with the dispersion relation and linewidth that we obtained from the power spectra as visualized in Sect. 3, Fig. 4. Whereas at small wavenumbers the travel times fall off slowly at larger time lags, at larger wavenumbers the travel times are increasingly concentrated in thin strips around Δt = 0. The decreasing phase speed with larger wavenumbers manifests itself as a steeper slope of the red and blue diagonal stripes.

For a quantitative comparison of the supergranular evolution pattern and the power spectra, we studied the time dependence of the convective signal at the origin (x = 0, y = 0) and average over a small circular region to increase the signal-to-noise ratio. This is similar to studying the autocorrelation of travel time and LCT divergence maps on a given spatial scale. The results are shown in Fig. 8. Since the autocorrelation is the inverse Fourier transform of the power spectrum, we are motivated to fit the inverse temporal Fourier transform of a Lorentzian to the curves:

$$f(k,\Delta t)=a(k)\hspace{0.17em}\cos ({\omega }_0(k)\Delta t)\hspace{0.17em}\text{exp}(-\gamma (k)|\Delta t|),$$(3)

where ω₀ and γ are as defined in Sect. 3.

The fit results are overplotted on Fig. 4 in Sect. 3. The parameters agree well with the fit parameters obtained from the power spectra. This means that we can identify the structures in the power spectra with the observed supergranular evolution pattern.

The spatial scale kR_⊙ = 271 is the limiting case where the oscillations (corresponding to the cosine term in Eq. (3)) are barely visible, so it is hard to extract a meaningful period beyond this scale. This is reflected by the low quality factor in Fig. 4. On larger scales though, oscillations are clearly present; lifetimes and periods can thus be measured simultaneously even if the periods are longer than the lifetimes (e-folding times).

Figure 8 also shows that LCT is more sensitive to smaller scales than the TD travel times. For example, the LCT divergence amplitudes at zero time lag are larger at kR_⊙ = 271 than at kR_⊙ = 124, whereas the travel-time amplitudes at kR_⊙ = 271 are only about half the amplitude at kR_⊙ = 124.

Figures 5 and 6 can now be better understood in light of the scale dependence: The average supergranule that was found without bandpass filtering is a power-weighted average over a range of spatial scales. The larger-scale flows have longer lifetimes, thus the shorter-lived smaller-scale flows decay more rapidly and the average scale increases with larger |Δt|. The average supergranular inflow apparently has a slightly shorter oscillation period (and a weaker divergence signal) than the outflow. A possible explanation is that slightly smaller scales (which might be concealed due to smoothing of the measurements) are selected on average for the inflows.

5. Magnetic field evolution

5.1. Observations

In addition to the horizontal flows, we studied the evolution of the line-of-sight magnetic field for the average supergranule. The data processing steps are outlined in Sect. 2. Figure 9 shows that the magnetic field follows the oscillations of the supergranular flows. When there is an average inflow, the magnetic field is concentrated and stronger than average. When there is an outflow, on the other hand, the magnetic field is diluted. However, the evolution is not entirely symmetric in time with respect to the reference time (Δt = 0, time of strongest outflow or inflow). For example, the cleared zone in the average outflow is still clearly present at Δt = 1.3 d, but is only starting to form at Δt = −1.3 d.

Fig. 9. Evolution of magnetic field (absolute line-of-sight field, mean subtracted) for the average supergranule. Top: outflow. Bottom: inflow. Blue is weaker than the average field strength and red is stronger than average.

For the average inflow, the maximum magnetic field strength is not reached at the reference time, but a few hours later. Assuming a smooth evolution, a temporal Fourier interpolation yields a time lag of about six hours (see right panel in Fig. 11).

5.2. Cork simulation

The observed time lag in the magnetograms raises the question if the magnetic field is merely passively advected by the supergranular flows or if it also plays an active part in shaping those flows. The latter was proposed by, for example, Crouch et al. (2007), who suggests that random clustering of magnetic elements could trigger downflows on supergranulation scales. Recently, Chatterjee et al. (2017) found evidence for a solar-cycle dependence of supergranular parameters (e.g., the cell size) in century-long Ca II K observation records.

To study how much the magnetic field of the average supergranule is influenced by the supergranular flows, we conducted a cork simulation, ignoring magnetic polarity. Magnetic elements were modeled as non-interacting floating corks with a finite lifetime t_life = 16 h that corresponds to the flux replacement time (Hagenaar et al. 2003). The corks are advected by supergranulation-scale flows and simultaneously perform a random walk owing to the statistical influence of smaller-scale convective flows (as seen by Orozco Suárez et al. 2012; Roudier et al. 2016; Agrawal et al. 2018, for the internetwork field) with a turbulent diffusivity η = 250 km² s⁻¹ (Simon & Weiss 1997). Further correlations in the smaller-scale flows, for example, the coherent flows on mesogranular scales that build up in trees of fragmenting granules, as described by Roudier et al. (2016), were neglected. Details of the simulation are given in Appendix A.4.

We then computed the cork density by counting the number of corks in each pixel (on the observation grid) for a given time step. To facilitate the comparison with the observed magnetic field, we smoothed the cork density temporally by applying an eight-hour boxcar moving average and spatially by convolution with a Gaussian of 0.7 Mm FWHM, similar to the HMI point-spread function (Yeo et al. 2014).

Despite the limitations of the simulation, the evolution of the magnetic field distribution (Fig. 9) is well reproduced by the corks (see Fig. 10)². In particular, the time lag of the strongest cork concentration in the average supergranular inflow matches the time lag of the strongest magnetic field concentration (see Fig. 11). The spatial peak size is reduced though (FWHM is about 5 Mm compared to 7 Mm).

Fig. 10. Simulated evolution of cork density (mean subtracted) for the average supergranule. Top: outflow. Bottom: inflow. Blue is weaker than the average cork density, red is stronger than average. The cork density is taken from a simulation run with η = 250 km2 s−1 and tlife = 16 h (run #5).

Fig. 11. Maximum and minimum magnetic field (observations) and cork density (simulation) over each average supergranule map as a function of time lag. Left: outflow. Right: inflow. The values are shown after subtraction of mean and division by the maximum absolute value to facilitate the comparison. The cork density is taken from a simulation run with η = 250 km2 s−1 and tlife = 16 h (run #5).

To analyze how sensitive the time lag is to the cork simulation parameters, we performed a parameter study. We varied the turbulent diffusivity and the cork lifetime, running nine simulations for various combinations of these parameters (see Table A.2). The original simulation run with η = 250 km² s⁻¹ and t_life = 16 h serves as our reference (run #5 in the table).

For each simulation run and the observed magnetic field, we measured the time of highest cork density and strongest magnetic field with respect to the time of strongest horizontal inflow. To this end, we constructed the time series of maximum cork density and magnetic field values in each frame (as shown by the red curves in Fig. 11) and used two different fitting methods to obtain the time lag from these time series: a parabolic fit using only the three temporal pixels around the maximum, and a Lorentzian fit using the entire time series. The resulting time lags for the different parameter combinations are shown in Table A.2. All simulation runs produce positive time lags in the range between 3.6 and 7.7 h, with the majority lying within one hour of the observed time lag for the magnetic field. Smaller differences are likely not to have a strong meaning, judging from the differences between the two fitting methods, which are typically 0.5–1 h. However, there are (weak) trends of decreasing time lags for increasing diffusivity, and increasing time lags for increasing lifetimes. This results in an underestimation of the time lag for run #7 and an overestimation for runs #2 and #3.

For the average supergranular outflow, we were able to qualitatively reproduce the clearing of the magnetic field in the outflow center, however with a shallower minimum compared to the observed field (see Fig. 12, for simulation run #5). Around the outflow, a ring of increased cork density forms, corresponding to the magnetic network field. This ring of higher cork density stands out more clearly than for the observed magnetic field.

Fig. 12. Simulated cork density (left) and observed magnetic field (right) for the average supergranular outflow at reference time (Δt = 0), after subtraction of mean. The cork density is taken from a simulation run with η = 250 km2 s−1 and tlife = 16 h (run #5). The circles have radii 11 Mm and 18 Mm with respect to the origin.

The observed network field is stronger on the west side of the average supergranule than on the east side, as we detected in our previous study (Langfellner et al. 2015a) and as was later also found by Roudier et al. (2016). This observation is successfully reproduced by the cork simulation (run #5). To study the dependence of this anisotropy on the simulation parameters, we analyzed the azimuthal distribution of the network field by averaging the cork density in a ring (see circles in Fig. 12) around the outflow center over azimuthal bins of 15°. The results for all simulation runs are shown in Fig. 13. The anisotropy is similar to the observed magnetic field for all runs with t_life = 16 h and 24 h. For shorter lifetimes (t_life = 8 h), the anisotropy is much less pronounced, if present at all. In these cases, the corks apparently vanish before the supergranular flows had enough time to aggregate the corks on the western side of the supergranules. This is consistent with the observation by Gošić et al. (2014) that the net contribution of the internetwork field to the network field is positive (amplification stronger than cancellation) and the most important contribution factor, with a flux replacement time of 18–24 h. Also, a higher than observed diffusivity (η = 500 km² s⁻¹) results in a weaker anisotropy. Overall, these results support our previous speculation that the supergranular oscillatory flows acting on the magnetic field are able to produce the network anisotropy.

Fig. 13. Simulated cork density in a ring surrounding the outflow center (11 Mm < r < 18 Mm, as shown in Fig. 12) as a function of azimuthal angle, ψ. Results are shown for nine separate simulation runs with different combinations of η and tlife parameter values (see Table A.2). The angle ψ = 0 points west. All results are at reference time (Δt = 0), after subtraction of the mean over the image. We show the observed magnetic field in all three panels for reference.

6. Summary and discussion

6.1. Supergranular wave power and evolution

We have confirmed the supergranular power spectrum of the divergence signal measured by Gizon et al. (2003) and Schou (2003) at the solar equator and extended it to higher wavenumbers, with a good agreement between TD and LCT of granules. As was found in the previous studies, the wave power is strongest in the prograde (westward) direction. Lorentzian fits to the power spectrum yield at kR_⊙ > 120 a strong deviation from the empirical square-root-of-k dispersion relation in form of a flat profile with an oscillation frequency of ∼2 μHz. For kR_⊙ > 270, the line width is too broad to extract a meaningful dispersion relation.

For the average supergranule in real space, we detected a complex evolution pattern with an oscillatory component, both with TD and LCT. Average supergranular outflows and inflows turn into each other with a mean period of six to seven days, viewed at a tracking rate corresponding to the supergranulation pattern rotation rate (Carrington rate plus 40 m s⁻¹). This means that new supergranules are preferentially born in inflow regions west of existing supergranules.

Indeed, the supergranular evolution pattern is scale-dependent; restricting the selection of supergranules to narrow size ranges in the averaging process reveals an increase of the oscillation period and lifetime for larger cells. These trends are consistent with the parameter values extracted from the Lorentzian fits of the supergranular power spectra. We were thus able to relate the signature of supergranular waves in Fourier space to the evolution pattern of the average supergranule in real space.

Changing the tracking rate to a lower value (Carrington rate minus 20 m s⁻¹) introduces a westward drift of the evolution pattern corresponding to the tracking speed difference (60 m s⁻¹), but does not alter the oscillation pattern. In Fourier space, this drift corresponds to a Doppler frequency shift ku_x.

6.2. Depth of supergranules

As we showed in Fig. 3, the background flow u_x shows a slight dependence on spatial scale. A potential explanation for this might be that supergranules of different sizes are rooted at different depths and are sensitive to the change of the rotation rate in the near-surface shear layer.

To get a rough idea if this might be the case, we relate u_x(k) to measurements of the solar rotation rate from global-mode helioseismology. Let us assume that supergranules have a depth z that depends linearly on their horizontal extent L, with a scaling factor α, i.e., L = αz. With L = 2π/k and z = R_⊙ − r_base, where r_base is the radial distance from the solar center to the base of the supergranule, this gives a relationship between kR_⊙ and r_base that depends on α

$$k{R}_{\odot }=\frac{2\pi }{\alpha \left(1-r_{\text{base}}/R_{\odot }\right)}\cdot$$(4)

We further speculate that the velocity u_x(k) that is “felt” by a supergranule is some weighted average of the rotation rate Ω(r) in the range of r between the top and the base of the supergranule. Let us consider two scenarios for the weights.

In the first scenario, we suppose that supergranules feel an average of the linear velocity Ω(r)r between r_base and the solar surface (R_⊙):

$$u_1(r_{\mathrm{base}})=\frac{{\displaystyle {\int }_{r_{\mathrm{base}}}^{R_{\odot }}\Omega }(r)r\mathrm{d}r}{R_{\odot }-r_{\text{base}}}\cdot$$(5)

In a second scenario, motivated by the idea that the motion of supergranules at the surface reflects the angular velocity at their rooting depths, we suppose that supergranules feel the average of the angular velocity Ω(r). This gives a higher weight to greater depths compared to the first scenario, as Ω(r) is multiplied by the solar radius R_⊙ at all depths instead of the local radius r ≤ R_⊙:

$$u_2(r_{\mathrm{base}})=\frac{{\displaystyle {\int }_{r_{\text{base}}}^{R_{\odot }}\Omega }(r)R_{\odot }\mathrm{d}r}{R_{\odot }-r_{\text{base}}}\cdot$$(6)

The left panel of Fig. 14 shows the linear rotational velocity Ω(r)r and the rotation rate Ω(r) (multiplied by R_⊙ for ease of comparison), where Ω(r) is the equatorial rotation rate from a 2D global-mode RLS inversion of HMI Dopplergrams (360-day average in 2010–11, Larson & Schou 2018). As the rotational velocity never exceeds 1980 m s⁻¹ for any depth, the weighted average u₁ is lower than the supergranular velocity (see right panel) by at least 50 m s⁻¹, irrespective of α.

Fig. 14. Left: linear rotational velocity Ω(r)r and rotation rate Ω(r)R⊙ vs. radial distance from solar center at solar equator obtained from a global-mode inversion (Larson & Schou 2018). Right: vertical average of rotation rate Ω(r)R⊙ between base of supergranule and solar surface vs. kR⊙, for various α values (see text for definitions of u2 and α). Also shown is ux plus tracking speed (2034 m s−1) from supergranular power spectrum fits (TD and LCT fast tracking). We note that the y-axis range is different in both panels.

On the other hand, the quantity Ω(r)R_⊙ reaches 2040 m s⁻¹. With Eq. (4), we can express the weighted average u₂ as a function of horizontal scale kR_⊙. It is shown in Fig. 14 for different values of α, together with u_x from the supergranular power fits where the fast tracking speed (2034 m s⁻¹) was added. While only considering the angular velocity Ω(r_base)R_⊙ at the base of the supergranule would yield projected velocity values that are 10 m s⁻¹ too high, the average u₂ between r_base and R_⊙ matches the supergranular velocity profile quite closely for α ∼ 0.4. In our simple picture, this corresponds to vertically elongated supergranules, for example, L = 30 Mm would yield a depth z ∼ 75 Mm, which is deeper than the near-surface shear layer.

While this match might be pure coincidence, it motivates further study of the supergranular wave properties, for example, using depth inversions, to explore the possibility of deep-rooted supergranules, as originally suggested by Parker (1973). We note that other recent works typically reported shallow supergranules (e.g., Rieutord et al. 2010; Duvall et al. 2014).

6.3. Passive magnetic field

We found that the magnetic field undergoes similar oscillations as the horizontal divergent flows on supergranular scales, but with a time lag of about six hours. The magnetic field evolution, including the time lag, can be well reproduced with a simple simulation, where magnetic elements are modeled as corks. In this simulation, the corks are advected by the average supergranular flows, possess a finite lifetime and are affected by smaller-scale motions through the turbulent diffusivity. The cork simulation is also able to reproduce the previously observed east-west anisotropy of the network field (Langfellner et al. 2015a; Roudier et al. 2016) surrounding the average supergranule.

The occurrence of the six-hour time lag both in the magnetic field observations and the cork simulation is indicative of the passive nature of the magnetic field in supergranules and confirms the findings of Simon & Weiss (1989), Orozco Suárez et al. (2012) and Roudier et al. (2016). In our measurements and simulations, there is no sign of magnetic fields shaping the supergranular flows as proposed by Crouch et al. (2007) or network field arising from random walks (Thibault et al. 2012). Instead, the higher cork densities in our simulations, which correspond to the network field, arise solely due to the advection by the supergranular flows.

However, an influence of the magnetic fields on the (quiet-Sun) supergranular flows in certain cases cannot be completely excluded. For example, Sangeetha & Rajaguru (2016) found that the vorticity in supergranular inflows is reduced if the magnetic field is strong. According to Kobel et al. (2012), in the strongest parts of the magnetic network the convection can be inhibited, causing magnetic elements to appear darker than in regions with a weaker magnetic field. It is unclear though how far the influence of the magnetic field goes.

Attempts to relate the magnetic field strength to supergranulation parameters are difficult; for example, Chatterjee et al. (2017) reported that the supergranule size increases in active regions with stronger magnetic field, whereas a strong quiet-Sun network field is associated with smaller supergranules. Meunier et al. (2007), on the other hand, found that a strong network field is associated with bigger cells and a strong internetwork field is responsible for smaller cells. Thus, the results depend on the exact definition and measurement place of the magnetic field strength and are difficult to compare and interpret.

Movies

Movie of Fig. B1 Access here

Movie of Fig. B1 Access here

Movie of Fig. B2 Access here

Movie of Fig. B2 Access here

Acknowledgments

The HMI data used are courtesy of NASA/SDO and the HMI science team. The data were processed at the German Data Center for SDO (GDC-SDO), funded by the German Aerospace Center (DLR). We thank J. Schou, R. Cameron and T. Duvall, Jr. for the useful discussions. We are grateful to R. Burston and H. Schunker for providing help with the data processing, especially the tracking and mapping. We used the workflow management system Pegasus funded by The National Science Foundation under OCI SI2-SSI program grant #1148515 and the OCI SDCI program grant #0722019 as well as the distributed computing software HTCondor.

References

Appendix A: Additional plots and tables

A.1. Azimuthal dependence of power

This Appendix consists of Figs. A.1 through A.4.

A.2. Fit parameters

This Appendix consists of Table A.1.

A.3. Evolution: dependence on tracking rate and spatial scale

This Appendix consists of Fig. A.5.

A.4. Cork simulation details

The cork simulations that we ran to model the magnetic field on supergranular scales (see Sect. 5.2) span 11 days in time (as long as the observation datasets) with a time step Δt_sim = 15 min on a domain of size 150 × 150 Mm². The simulations were initialized with N_init ≈ 1.8 × 10⁵ corks at random positions of uniform distribution, corresponding to one cork per square pixel on average. This number effectively sets the noise level; the more corks, the lower the noise. At each time step, the simulations were advanced by applying the following actions:

1. Advection of existing corks by horizontal supergranular background flow. We used the time-dependent u = (u_x, u_y) measured with LCT for the average supergranule that was also used to compute div_hu shown in Fig. 6.

2. Random walk of existing corks. The average distance was determined by the turbulent diffusivity η = 250 km² s⁻¹ (Simon & Weiss 1997) for the standard run. For other runs, we also used 125 and 500 km² s⁻¹ (see Table A.2).

3. Corks that leave the simulation domain are removed.

4. The remaining corks have a probability p_dec = 1 − exp(−Δt_sim/t_life) to be removed from the simulation (exponential decay). The cork lifetime t_life was chosen to be equal to the flux replacement time t_repl ∼ 16 h (Hagenaar et al. 2003) for the standard run. For other runs, we also used 8 and 24 h.

5. To compensate for the removed corks and to simulate the emergence of new flux, we spawn N_spawn = p_decN_init new corks at random positions with uniform probability.

Fig. A.1. Supergranular power spectrum at the solar equator: Azimuthal dependence at kR⊙ ≈ 123 for the Snodgrass tracking rate using TD (f modes) and LCT. Top row: observed power. Bottom row: Lorentzian fit. Frequencies beyond ±7 μHz are not taken into account.

Fig. A.2. As Fig. A.1, but for the fast tracking rate.

Fig. A.3. As Fig. A.2, but for kR⊙ = 74.

Fig. A.4. As Fig. A.2, but for kR⊙ = 270.

Fig. A.5. Evolution of f-mode-filtered inward–outward point-to-annulus travel times for the average supergranular outflow, with the Snodgrass (1984) tracking rate (cf. Fig. 5: tracking at supergranulation pattern rotation rate, which is faster by about 60 m s−1).

Appendix B: Evolution movies

This paper is accompanied by four movies that show the temporal evolution of the average supergranule. Figure B.1 shows a still from the evolution movie of the average supergranular outflow, which juxtaposes TD and LCT (see caption for details). The corresponding movie for the average supergranular inflow (no still shown) is composed analogously.

Figure B.2 shows a still from the movie that directly compares the observed evolution of the line-of-sight magnetic field and the cork simulation (run #5) for the average supergranular outflow. In this case, too, the corresponding movie for the average supergranular inflow has the same structure, and thus no still is shown.

Fig. B.1. Still image from supergranular temporal evolution movie at time lag Δt = 30 h. Top row: TD (inward–outward travel times) and LCT divh for the average supergranular outflow, as shown in Figs. 5 and 6. Frames at time lags that are not multiples of eight hours were calculated using temporal Fourier interpolation. Bottom row: maximum and minimum travel time and LCT divh over each average supergranule map as a function of time lag. The highest absolute value at the current time lag is marked with the black circle. This value is reflected by the limits of the colorbar in the top row panels. The full temporal evolution for both, the inflow and the outflow, can be found in the two movies attached to this figure.

Fig. B.2. Still image from supergranular evolution movie. Top row: simulated cork density (run #5) and observed line-of-sight magnetic field for the average supergranular outflow, as shown in Figs. 10 and 9. Bottom row: maximum and minimum cork density and magnetic field strength as in Fig. 11. In addition, the highest absolute value at the current time is marked with the black circle. The full temporal evolution for both, the inflow and the outflow, can be found in the two movies attached to this figure.

All Tables

All Figures

Fig. 1. Supergranular power spectrum of the divergence signal at the solar equator using time-distance helioseismology. Top row: at ky = 0. The left panel shows the spectrum for a tracking rate set to the Snodgrass (1984) rotation rate (about 20 m s−1 slower than Carrington). The right panel shows the spectrum for a tracking rate that is faster by about 60 m s−1. Bottom row: cut through k plane at frequency 3.2 μHz.

In the text

	Fig. 2. Supergranular power spectrum of the divergence signal at the solar equator for both tracking rates, cut at kR⊙ ≈ 123, using one year of HMI data: time-distance helioseismology (left) and local correlation tracking right). Top row: East-west direction. Bottom row: North-south direction.
In the text

Fig. 3. Parameters from Lorentzian fits of time-distance and LCT divergence power spectra for Snodgrass (1984) and fast tracking rates. Top row: dispersion relation and half-width at half maximum (HWHM). Bottom row: horizontal velocity components of background flow. The tracking rate difference of 60 m s−1 has been subtracted from the ux values for Snodgrass (1984) tracking to facilitate the comparison with fast tracking. Values for kR⊙ > 230 (Snodgrass tracking) and kR⊙ > 271 (fast tracking) are omitted, as the fits are not stable in that regime.

In the text

	Fig. 4. Parameters extracted from Lorentzian fits of divergence power spectra and the average convection cell evolution (see Sect. 4.2) for fast tracking rate at the solar equator, as a function of spatial scale. Values for kR⊙ > 271 are omitted as the fits are not stable in that regime. For f0 the dispersion relations from Gizon et al. (2003) and Schou (2003) are overplotted for comparison.
In the text

Fig. 5. Evolution of f-mode-filtered inward–outward travel times for the average supergranule. Top: outflow. Bottom: inflow. The frames, from top left to bottom right, span a total of eleven days in non-overlapping steps of eight hours. The central frame (red outline) is for the reference time, Δt = 0. We note that the color scale is different for each individual frame; saturation is reached at the maximum absolute value (given by the lower label boxes) of the respective frame. The color map is symmetrized around zero; outflows are red and inflows are blue.

In the text

	Fig. 6. Evolution of LCT divh for the average supergranular outflow (cf. top of Fig. 5).
In the text

	Fig. 7. Evolution for different scales: space-time images. Average over y band (10 Mm) is given for the average outflow. Top row: TD. Bottom row: LCT. The dashed lines give the phase speed from the dispersion relation as obtained from the power spectrum. Outflows are red and inflows are blue. To increase visibility, the color maps are saturated at half the maximum value.
In the text

	Fig. 8. Evolution of supergranular cells on different scales: Average over circular area with radius (2πR⊙/k)/4 around the origin. Top row: TD outflow and inflow (sign reversed). Bottom row: LCT outflow and inflow.
In the text

	Fig. 9. Evolution of magnetic field (absolute line-of-sight field, mean subtracted) for the average supergranule. Top: outflow. Bottom: inflow. Blue is weaker than the average field strength and red is stronger than average.
In the text

	Fig. 10. Simulated evolution of cork density (mean subtracted) for the average supergranule. Top: outflow. Bottom: inflow. Blue is weaker than the average cork density, red is stronger than average. The cork density is taken from a simulation run with η = 250 km2 s−1 and tlife = 16 h (run #5).
In the text

	Fig. 11. Maximum and minimum magnetic field (observations) and cork density (simulation) over each average supergranule map as a function of time lag. Left: outflow. Right: inflow. The values are shown after subtraction of mean and division by the maximum absolute value to facilitate the comparison. The cork density is taken from a simulation run with η = 250 km2 s−1 and tlife = 16 h (run #5).
In the text

	Fig. 12. Simulated cork density (left) and observed magnetic field (right) for the average supergranular outflow at reference time (Δt = 0), after subtraction of mean. The cork density is taken from a simulation run with η = 250 km2 s−1 and tlife = 16 h (run #5). The circles have radii 11 Mm and 18 Mm with respect to the origin.
In the text

Fig. 13. Simulated cork density in a ring surrounding the outflow center (11 Mm < r < 18 Mm, as shown in Fig. 12) as a function of azimuthal angle, ψ. Results are shown for nine separate simulation runs with different combinations of η and tlife parameter values (see Table A.2). The angle ψ = 0 points west. All results are at reference time (Δt = 0), after subtraction of the mean over the image. We show the observed magnetic field in all three panels for reference.

In the text

Fig. 14. Left: linear rotational velocity Ω(r)r and rotation rate Ω(r)R⊙ vs. radial distance from solar center at solar equator obtained from a global-mode inversion (Larson & Schou 2018). Right: vertical average of rotation rate Ω(r)R⊙ between base of supergranule and solar surface vs. kR⊙, for various α values (see text for definitions of u2 and α). Also shown is ux plus tracking speed (2034 m s−1) from supergranular power spectrum fits (TD and LCT fast tracking). We note that the y-axis range is different in both panels.

In the text

	Fig. A.1. Supergranular power spectrum at the solar equator: Azimuthal dependence at kR⊙ ≈ 123 for the Snodgrass tracking rate using TD (f modes) and LCT. Top row: observed power. Bottom row: Lorentzian fit. Frequencies beyond ±7 μHz are not taken into account.
In the text

	Fig. A.2. As Fig. A.1, but for the fast tracking rate.
In the text

	Fig. A.3. As Fig. A.2, but for kR⊙ = 74.
In the text

	Fig. A.4. As Fig. A.2, but for kR⊙ = 270.
In the text

	Fig. A.5. Evolution of f-mode-filtered inward–outward point-to-annulus travel times for the average supergranular outflow, with the Snodgrass (1984) tracking rate (cf. Fig. 5: tracking at supergranulation pattern rotation rate, which is faster by about 60 m s−1).
In the text

Fig. B.1. Still image from supergranular temporal evolution movie at time lag Δt = 30 h. Top row: TD (inward–outward travel times) and LCT divh for the average supergranular outflow, as shown in Figs. 5 and 6. Frames at time lags that are not multiples of eight hours were calculated using temporal Fourier interpolation. Bottom row: maximum and minimum travel time and LCT divh over each average supergranule map as a function of time lag. The highest absolute value at the current time lag is marked with the black circle. This value is reflected by the limits of the colorbar in the top row panels. The full temporal evolution for both, the inflow and the outflow, can be found in the two movies attached to this figure.

In the text

Fig. B.2. Still image from supergranular evolution movie. Top row: simulated cork density (run #5) and observed line-of-sight magnetic field for the average supergranular outflow, as shown in Figs. 10 and 9. Bottom row: maximum and minimum cork density and magnetic field strength as in Fig. 11. In addition, the highest absolute value at the current time is marked with the black circle. The full temporal evolution for both, the inflow and the outflow, can be found in the two movies attached to this figure.

In the text