Issue 
A&A
Volume 635, March 2020



Article Number  A181  
Number of page(s)  15  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/201937331  
Published online  31 March 2020 
Characterizing the spatial pattern of solar supergranulation using the bispectrum
^{1}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: boening@mps.mpg.de
^{2}
GeorgAugustUniversität, Institut für Astrophysik, FriedrichHundPlatz 1, 37077 Göttingen, Germany
^{3}
Center for Space Science, NYUAD Institute, New York University Abu Dhabi, PO Box 129188, Abu Dhabi, UAE
Received:
17
December
2019
Accepted:
5
February
2020
Context. The spatial power spectrum of supergranulation does not fully characterize the underlying physics of turbulent convection. For example, it does not describe the nonGaussianity in the horizontal flow divergence.
Aims. Our aim is to statistically characterize the spatial pattern of solar supergranulation beyond the power spectrum. The nextorder statistic is the bispectrum. It measures correlations of three Fourier components and is related to the nonlinearities in the underlying physics. It also characterizes how a skewness in the dataset is generated by the coupling of three Fourier components.
Methods. We estimated the bispectrum of supergranular horizontal surface divergence maps that were obtained using local correlation tracking (LCT) and timedistance helioseismology (TD) from one year of data from the helioseismic and magnetic imager onboard the solar dynamics observatory starting in May 2010.
Results. We find significantly nonzero and consistent estimates for the bispectrum using LCT and TD. The strongest nonlinearity is present when the three coupling wave vectors are at the supergranular scale. These are the same wave vectors that are present in regular hexagons, which have been used in analytical studies of solar convection. At these Fourier components, the bispectrum is positive, consistent with the positive skewness in the data and consistent with supergranules preferentially consisting of outflows surrounded by a network of inflows. We use the bispectral estimates to generate synthetic divergence maps that are very similar to the data. This is done by a model that consists of a Gaussian term and a weaker quadratic nonlinear component. Using this method, we estimate the fraction of the variance in the divergence maps from the nonlinear component to be of the order of 4–6%.
Conclusions. We propose that bispectral analysis is useful for understanding the dynamics of solar turbulent convection, for example for comparing observations and numerical models of supergranular flows. This analysis may also be useful to generate synthetic flow fields.
Key words: hydrodynamics / turbulence / convection / Sun: interior / Sun: helioseismology
© Vincent G. A. Böning et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1. Introduction
Supergranular motions are characterized by horizontal velocities of about 300 m s^{−1} and are coherent features with a length scale of about 30 Mm or harmonic degrees of about l ∼ 120 in the power spectrum (e.g., Rincon & Rieutord 2018, and references therein). Close to the surface, supergranules are observed as regions of flow with horizontal divergence (associated with upflows) surrounded by a network of convergent flows (associated with downflows, e.g., Duvall & Gizon 2000).
The origin of supergranulation as a dominant scale of convection remains unclear (see Rincon & Rieutord 2018, for a review). Supergranulation as a particular lengthscale might emerge because of thermal instabilities in the ionization zones (Simon & Leighton 1964), as the largest buoyantly driven scale of convection in the Sun (Cossette & Rast 2016), or the largest convective scale not significantly constrained by rotation (Featherstone & Hindman 2016). Also unexplained is the evolution of the supergranulation pattern, which displays a wavelike behavior (Gizon et al. 2003; Schou 2003) – although this is not the topic of this paper.
Comparisons between theory and observations of supergranularscale flow fields are usually based on the power spectrum of horizontal surface flows (e.g., Rincon et al. 2017). However, it is clear that the power spectrum does not contain much information about the spatial pattern of the flow field, apart from the scaledependence of the amplitude. For example, the power spectrum does not contain any information about the skewness of the underlying field (see e.g., Fig. 5 in Nordlund et al. 1997). Divergence maps of supergranulation have a positive skewness (Duvall & Gizon 2000) which is due to the asymmetric distribution of upflows and downflows.
Such an asymmetry between up and downflows is present in idealized hexagonal convection cells that have been used in analytical studies of convection and supergranulation in particular (e.g., Chandrasekhar 1961; Gough et al. 1975; Busse 2007). These hexagonal cells are characterized by a central upflow region surrounded by a network of downflows (Chandrasekhar 1961; Toomre et al. 1977; Busse 2003), similar to supergranulation. The opposite case of strong downflows surrounded by a weaker network of upflows was also considered (Busse 2004, 2007). For example, theoretical attempts to explain the phase velocity of supergranulation and its contribution to maintaining the nearsurface shear layer (Busse 2004, 2007) are based on decomposing the solution to the NavierStokes equations in terms of planforms that are constructed from such patterns in the horizontal direction.
The main aim of our work is to provide a statistical characterization of the spatial pattern of supergranulation beyond the power spectrum. In other words, we want to find out whether convective patterns, such as hexagons, can be directly related to the supergranulation data. Additionally, our longerterm aim is to construct synthetic supergranular flow fields. The availability of synthetic flow fields could be useful in helioseismic methods for studying subsurface convection (e.g., Hanasoge et al. 2012; Duvall et al. 2014).
Information on the skewness is contained in the bispectrum, which is a nextorder generalization of the power spectrum (e.g., Priestley 1981, Sect. 11.5.1). Just as the power spectrum describes the variance as a function of scale, the bispectrum describes the skewness as a function of two coupling scales (e.g., Watkinson et al. 2017). Measuring the bispectrum of supergranular divergence maps would therefore be useful to characterize the skewness as a function of scale. Furthermore, the bispectrum is directly related to hexagonal patterns because a hexagon has a specific signature in the bispectrum. The bispectrum characterizes the coupling of two wave vectors to a third wave vector, which appears similarly in the nonlinearities in the NavierStokes equations. Such nonlinearities are present in the modal equations of convection studied using horizontal patterns (e.g., Chandrasekhar 1961; Gough et al. 1975; Busse 2007).
To clarify the results obtained in this paper, it is necessary to clarify what we mean by linear and nonlinear. By linear, we mean the hypothesis that the timeevolution of supergranular divergence maps was given by a linear partial differential equation (PDE) and a stochastic forcing term with independent Gaussian Fourier components. In this case, the dynamics at the Fourier components of the solution would also be independent. By nonlinear, we mean the hypothesis that the dynamics of supergranular divergence maps was described by a PDE that in addition includes a nonlinear term. In this case, the Fourier components of the solution would not be independent but would show some dependence that might be measurable. Supergranular dynamics are strongly nonlinear according to Rincon et al. (2017) and divergence maps have a large nonzero skewness (Duvall & Gizon 2000). This should be reflected in the bispectrum.
Introductions to bispectral analysis are given by Brillinger (1965), Kim & Powers (1979), Priestley (1981, Sect. 11.5.1), Nikias & Raghuveer (1987), and Chandran & Elgar (1990). Bispectral analysis has been used in a number of research domains. For example, it was used to measure nonlinear wave–wave coupling (e.g., Kim et al. 1980) of ocean and surface water waves (e.g., Elgar 1989), to study nonGaussianity in the cosmic microwave background (e.g., Komatsu & Spergel 2001; Planck Collaboration XVII 2016), and in laboratory plasma physics (e.g., Itoh et al. 2017). The bispectrum is related to the spectral transfer of kinetic energy in turbulent flows (e.g., Yeh & van Atta 1973). In the context of stellar convective dynamos, measurements of the spectral energy transfer in magnetohydrodynamical simulations have been taken for example by Brandenburg (2001), Moll et al. (2011), and Strugarek et al. (2013).
In this work, we apply bispectral analysis to nearsurface horizontal divergence maps of supergranulation obtained by Langfellner et al. (2018) using local correlation tracking (LCT, November & Simon 1988) and timedistance helioseismology (TD, Duvall et al. 1993) and data from the Helioseismic and Magnatic Imager (HMI, Scherrer et al. 2012; Schou et al. 2012) onboard the Solar Dynamics Observatory (SDO, Pesnell et al. 2012).
2. Divergence maps of supergranular flows
We use maps of the horizontal nearsurface divergence that were obtained by Langfellner et al. (2018) using Dopplergrams and intensity images from HMI/SDO. The data cover 365 days from May 2010 through April 2011.
Regions of size ∼180 × 180 Mm^{2}, which are centered at the equator, were tracked for eleven days from about 70° east to 70° west at the surface rotation rate (Snodgrass 1984). Dopplergrams and intensity maps were remapped and each elevenday timeseries was subdivided into 33 eighthour segments. The seventeenth segment crosses the central meridian.
Using the intensity images and LCT of granules (November & Simon 1988), Langfellner et al. (2018) obtained estimates of the two horizontal flow components for the same spatial and temporal coverage. This results in maps with 97 × 97 pixels, each of size (5 × 0.348 Mm)^{2}. These maps span a slightly smaller region of size ∼170 × 170 Mm^{2}. From these maps, we computed the horizontal divergence by taking derivatives using Savitzky–Golay filters (Savitzky & Golay 1964). In the Savitzky–Golay filters, we used a different filter width (3 instead of 15 pixels) and polynomial order (2 instead of 3) compared to Langfellner et al. (2018) to keep the smoothing as small as possible.
For each segment, Dopplergrams were filtered in the Fourier domain using an fmode ridge filter (Langfellner et al. 2015) which includes power from wave numbers 300 < kR_{⊙} < 2600 with a peak at kR_{⊙} ≈ 1100. Here, R_{⊙} is the solar radius and k is the horizontal wave number. Next, fmode pointtoannulus (Duvall et al. 1996) traveltime maps (Gizon & Birch 2004) were obtained from the filtered Dopplergrams in the inward minus outward sense. This resulted in maps with 512 × 512 pixels, with each pixel being (0.348 Mm)^{2} in size. These traveltime maps are known to be a proxy of nearsurface horizontal divergence of the fluid motion (e.g., Langfellner et al. 2015). To remove the smallestscale noise from the data, we lowpass filter this dataset in the Fourier domain and keep Fourier components with kR_{⊙} < 1200.
In addition, we obtained a smoothed version of the (not lowpass filtered) TD maps in order to better compare to earlier results. We used a twodimensional Gaussian window with full width at half maximum of 5.22 Mm (15 pixels), which is similar to the smoothing done by Duvall & Gizon (2000).
For all data sets discussed here, the horizontal mean of each map has been removed. Figure 1 shows example LCT and TD maps. The most prominent features in these maps are outflows due to supergranules, surrounded by inflows. Figure 1 also shows histograms of the maps. Both LCT and TD histograms are asymmetric about zero and have a nonzero skewness, which is well known (e.g., Duvall & Gizon 2000). The skewness of the ith map X^{(i)}(x) is
Fig. 1. Example divergence maps (top) and corresponding histograms (bottom) for LCT (left), TD (middle), and smoothed TD data (right). The smoothed TD data was apodized with a raised cosine and the histogram was computed for the region inside the dashed rectangle. The histogram and skewness were computed from the entire dataset. 
where the variance is
Here, x = (x, y) denotes the horizontal location, ⟨ ⋅ ⟩ is a horizontal average, ⟨X^{(i)}(x)⟩ = 0 by construction, X^{(i)}(k) denotes the horizontal Fourier transform of X^{(i)}(x) – see Appendix A –, and h_{k} is the resolution in Fourier space. We use the same symbol for a function and its Fourier transform. We use for the (mean) skewness computed from the entire dataset.
The skewness depends slightly on disk position. In the following, we therefore only use the 11 locations that are closest to disk center for each segment, ranging from about 23° east to 23° west, where the skewness is approximately constant. For the entire datasets, this results in a skewness of for the LCT maps, for the TD maps, and for the smoothed TD maps. Duvall & Gizon (2000) reported a skewness of 0.42 for TD maps. As we see in Sect. 5, different filters applied in the spatial domain may result in different values for the skewness.
The fact that the divergence maps have a nonzero skewness suggests that different Fourier components of the divergence field cannot be independent. If they were independent, the data would have zero skewness, just like a random wave field. As the statistical properties of the supergranular divergence field are translation invariant, it is clear that two Fourier components are uncorrelated. As we see in Sect. 4, the nonzero skewness of the data can be related to nonzero correlations between three Fourier components, which leads to the definition of the bispectrum. As such correlations between three Fourier components are also present in simple patterns of convection, we now look at these patterns before introducing the bispectrum.
3. A brief account of convection patterns
We briefly discuss a few patterns used in analytical studies of convection. Related to the specific case of solar supergranulation, Busse (2003, 2004, 2007) studied hexagonal convection cells. This work is based on earlier and more general studies of Boussinesq convection that used convection patterns such as rolls, squares, and hexagons (e.g., Palm 1960; Chandrasekhar 1961; Toomre et al. 1977) and the decomposition of convective motions into planforms, which are made up from such patterns in the horizontal direction (e.g., Chandrasekhar 1961; Gough et al. 1975). Due to the continuity equation and under the assumption of a steady flow, the vertical flows and the horizontal divergence show the same pattern close to the surface.
We now look at a few possible patterns in detail. The top row in Fig. 2 shows simple convective rolls. Such a roll is given by a pattern X(x) for which
Fig. 2. Example patterns X(x) and their Fourier transform X(k). The wave vectors in the right small panels are: k_{0} (red), l_{0} (blue), and m_{0} (orange). The skewness γ and fractional power f_{NL} of the nonlinearly coupled component are given in each panel. Equations (3) and (4) give definitions of the patterns. 
where k_{0} = (k_{x, 0}, k_{y, 0}) is a chosen wave vector and we use the notation +c.c. to indicate that the previous terms are repeated with a complex conjugate. The top row in Fig. 2 shows three examples of rolls for three different wave vectors k_{0}, l_{0}, m_{0}, which are kept fixed for the entire section and are shown in the small subpanels of Fig. 2 as colored arrows. These wave vectors form a triad, that is k_{0} + l_{0} = m_{0}. The wave vectors chosen here have the same length and are at an angle of 60° (up to a small deviation caused by the discretization).
We now consider different patterns of the form
where A is the (possibly complex) amplitude of the third Fourier component. This pattern can be seen as two independent rolls (k_{0} and l_{0}) plus a nonlinear component (m_{0}). This nonlinear component is given by a quadratic coupling (Kim & Powers 1979),
where A is a coupling coefficient. If one thinks of X(m_{0}) as being due to such a coupling, the nonlinear fraction of the variance is
which is given in the panels in the bottom row in Fig. 2.
A hexagonal pattern is given by A = 1; see bottom left panel in Fig. 2. Equation (4) with A = 1 is equivalent to the horizontal dependence in Eq. (249) in Sect. 16c of Chandrasekhar (1961). Hexagonal patterns have the advantage that they can be used to model a convective pattern that is characterized by a central upflow region surrounded by a network of downflows (Chandrasekhar 1961; Toomre et al. 1977; Busse 2003), similar to supergranulation. In addition, hexagonal planforms are not symmetric with respect to a change in sign and they produce a nonzero skewness in the horizontal divergence, similar to the supergranular divergence maps (see Sect. 2); see the skewness of each pattern given in Fig. 2.
The case A = −1 corresponds to a pattern of inverted hexagons (see right panel in the bottom row of Fig. 2). This corresponds to downflows surrounded by a network of upflows as considered by Busse (2004, 2007).
The case A = 0 corresponds to rectangles, and we can vary A continuously between zero and one to obtain patterns with different coupling amplitudes, intermediate between rectangles and hexagons; see the example labeled “small coupling” (A = 0.3) in Fig. 2.
In addition to these patterns, it is possible to tessellate space by a pattern generated from any combination of three wave vectors with k_{0} + l_{0} = m_{0}, each with a separate, possibly complex Fourier amplitude. For a given pattern, the phases of X(k_{0}) and X(l_{0}) may be jointly put to zero by a translation in real space. Statistically, we do not expect any preferred direction for supergranular divergence maps near the equator, and it is therefore reasonable to assume X(k_{0}) = X(l_{0}) for illustrative purposes. The Fourier amplitudes of X may then be normalized to one as done above.
For the patterns considered here, one pattern fills the whole domain. On the Sun, however, every supergranule has its individual shape. It is therefore possible that a coupling is present not only at the Fourier components with equal length (k_{0}=l_{0}=m_{0}, as for regular hexagons), but potentially for any coupling triad of Fourier components. As a consequence, the signature in Fourier space is a more complex combination of the previous characteristics at several Fourier components. This motivates the definition of the bispectrum.
4. Bispectral analysis in two dimensions
Here we introduce and summarize the definition of the bispectrum and related quantities for twodimensional random processes. This part is mainly based on Sect. 11.5.1 in Priestley (1981), Kim & Powers (1979), and Chandran & Elgar (1990). Later in this section, we give a few examples and illustrate the definitions with the help of the patterns introduced in the previous section. In addition, we introduce a method for generating data with a given bispectrum.
4.1. Definition of the bispectrum
Here we consider a random field X(x) that depends on a twodimensional discretely sampled spacial variable x. We assume that all statistical properties of X(x) are translation invariant. This is the generalization of a stationary process to two (spatial) dimensions. Since X(x) is a random variable, it induces the definition of an expectation value, which we denote by 𝔼[⋅]. The spatial Fourier transform of X(x) is again denoted by X(k). By X^{(i)}(x), we denote an individual independent realization of X(x), that is, an individual map.
The bispectrum is the next highestorder statistic beyond the power spectrum, and is defined as
The threepoint correlation function is the Fourier transform of the bispectrum, just as the twopoint (auto)correlation function is the Fourier transform of the power spectrum:
where the expectation value on the lefthand side of Eq. (8) does not depend on x_{0} due to the assumed translation invariance of X(x).
Similarly, the skewness of the random field X(x), γ_{1}, which is the thirdorder moment normalized by the variance, is the sum of the variancenormalized bispectrum,
where the power spectral density P is defined by
Table B.1 provides a summary of key terms defined in this paper for later reference.
In the following, we assume that there is no preferred horizontal direction in the data. In this case, the bispectrum of the divergence field is real. The expectation value 𝔼[X(k)X(l)X^{*}(m)] is zero if m ≠ k + l because of the assumed translation invariance. Therefore, the triadic relation m = k + l is implicitly understood in the following and the third variable is omitted in bispectrumrelated quantities.
4.2. Estimating the bispectrum
In principle, the bispectrum can be estimated using Eq. (7) by simply averaging over independent realizations X^{(i)}(k) of the data (e.g., Kim & Powers 1979). See Appendix B for a formal definition of the estimator , which is obtained this way. We take this approach although it may not be optimal; Nikias & Raghuveer (1987) discuss other possible approaches.
In order to assess the significance of the estimated bispectrum, it is necessary to consider the bispectral coherence (also known as the bicoherence) b, which is a normalized bispectrum (e.g., Kim & Powers 1979),
The definition of b is independent of the convention used to define the Fourier transform because it is normalized and because the bispectrum is real in our case. We see later that the squared bicoherence b^{2} quantifies the fraction of power which is due to quadratic nonlinear coupling at m = k + l. For such a coupling to be present in the data and a bispectral measurement to be significant, it is therefore necessary that the bicoherence be significantly different from zero. The variance of an estimator of the bicoherence that is built by simply averaging over independent realizations for estimating expectation values in Eq. (12) – see Appendix B –, is bounded by (Kim & Powers 1979)
where N_{B} is the number of independent data realizations used in the estimate. The statistical 1σ error in the bicoherence estimate is therefore . Additionally, we average the results azimuthally – see Appendix B.1 – which results in a typical σ_{b} = 0.0035 at mR_{⊙} = 155.5.
4.3. A model with a weak quadratic nonlinearity
We here present a model that helps us to interpret the measured bispectrum signal and to generate synthetic data that have approximatively a given bispectrum. There are several ways to generate data with a given power spectrum and bispectrum; see for example Contaldi & Magueijo (2001) and Wagner et al. (2010). Such a method is not unique (e.g., Wagner et al. 2010). Following Wagner et al. (2010), we assume that the data can be represented as
where X_{G} is a Gaussian field and X_{NL} a small nonlinear component. In the Fourier domain, the Gaussian field follows a complex Gaussian distribution,
where 𝒩_{k} are Gaussian random variables with independent real and imaginary parts, zero mean and unit variance. We set to ensure that X_{G} is real, but otherwise the 𝒩_{k} are independent for different k.
We assume that the nonlinear component is obtained from the Gaussian one by an arbitrary quadratic nonlinearity,
where A_{k, l} denote coupling coefficients of modes k and l = m − k. To guarantee that X_{NL} is real, the coupling coefficients have to fulfill .
Synthetic data with a given bispectrum B_{0}(k, l) can then be generated by setting
Appendix C outlines further details of this model.
4.4. The bispectrum of hexagons
For the patterns introduced in Sect. 3, the bispectrum of the patterns is given by and the coupling coefficient as defined above is given by A_{k0, l0} = A/6. The factor of onesixth appears because in the model introduced in this section it is assumed that a coupling occurs at several Fourier components. In the case of only one triad of coupling Fourier components as for the patterns from Sect. 3, one has A_{k0, l0} = A.
The bicoherence of the patterns with nonzero coupling coefficient A is b(k_{0}, l_{0}) = 1. This is because the signal at Fourier component m_{0} is entirely due to the nonlinear coupling, and there is no noise or any other independent signal at this component.
4.5. The impact of noise on bispectral estimates
If noise is present in the data, this reduces the bispectral estimates. If we assume for example that
where ϵ(x) is the noise, which we assume to have independent Gaussian Fourier components and to be independent of X_{G} and X_{NL}. If the noise is simply added to the data, the bispectrum (as an expectation value) does not change; see Eqs. (7) and (C.6). However, the estimated coupling coefficient A_{k, l} and the bicoherence b are reduced in amplitude because of the increase in power from the noise; see Eqs. (12) and (17). It is therefore not possible to discriminate between a dataset with small noise and a small coupling coefficient, and another one with largeramplitude noise and a larger coupling coefficient.
In other words, for a given power spectrum, which is due to the signal plus noise, stronger noise results in a smaller bispectrum. For the case of possible hexagonal convection cells, it might therefore be that their signal in the bispectrum is reduced by noise.
5. Results
In this section, we show results for data averaged over the azimuthal direction of m. We do not find any directiondependence of the bispectrum.
5.1. The bispectrum of supergranulation
Figure 3 shows the estimated bicoherence from LCT. The bicoherence is significantly nonzero for a large number of combinations of wave vectors. The bicoherence is largest at wave vectors of similar amplitude, k ≈ l. In general, it is strongest if the coupling vectors are close to the supergranular peak in the power spectrum, in our case at mR_{⊙} = 155.5. At these wave vectors, the amount of power due to a nonlinear component in the data is largest.
Fig. 3. Bispectral coherence b(k, m) from LCT for selected values of mR_{⊙} = 77.7, 155.5, 233.2 ± 12.3. The vector mR_{⊙} is fixed for each plot and given by the orange arrow. The bispectral estimates were averaged over the direction of m, which is equivalent to assuming that there is no preferred horizontal direction and which results in a real bispectrum. The direction of m in the plot was chosen simply for convenience. The Fourier vectors for a triad k + l = m that corresponds to hexagons is given by the other two wave vectors (red: k, blue: l). The angles between the displayed Fourier vectors k, m, and l are not exactly 60 degrees because of the discretization. White color denotes locations where no estimate was obtained because one of the three vectors has zero length (as the mean was removed from each map, the bispectral coherence is not defined in this case). See Fig. 5 for a cut through the middle panel including error bars. See Figs. F.2–F.4 for further illustrations of these results and Fig. F.1 for the corresponding TD results. 
Figure 3 displays results as a function of k = (k_{x}, k_{y}) and m for selected values of m. The results are averaged in such a way that m (orange arrow) is rotated to be in the x direction and thereby is fixed; see Appendix B. Each data point is then plotted at the location corresponding to k relative to the direction of m. For each k, the third coupling Fourier component is then given by l = m − k. The strongest coupling is present at a location where the triad k + l = m corresponds to a hexagonal pattern at the supergranular scale; see middle panel in Fig. 3. Figures F.2 and F.4 illustrate these results further.
The results for TD are very similar for mR_{⊙} ≲ 300 (Fig. F.1 and second column in Table 1). For mR_{⊙} > 300, the signal considerably weakens compared to LCT and is dominated by noise and leakage due to the spatial smoothing. This is the case for both TD datasets.
Summary of results for the bispectral analysis of supergranular divergence maps from LCT and TD.
Estimates for the bispectrum are qualitatively very similar to the bicoherence; see Fig. F.3. As the sign of the bispectrum is predominantly positive, the dominant part of the quadratic interaction results in a positive contribution to the skewness; see Sect. 4.1. This is consistent with the positive skewness of the data; see Fig. 1. The peak in the bispectrum coincides with the peak in the bicoherence. The coupling modes that produce the strongest nonlinearity therefore also result in the largest contribution to the positive skewness in the data.
A positive sign in the bispectral coherence at wave vectors similar to hexagons indicates that the divergence field is preferentially given by outflows surrounded by a network of inflows (see the patterns discussed in Sect. 3). The outflows are on average slightly stronger while the inflows occupy a larger area. However, in individual realizations of the data, the phase of X^{(j)}(k)X^{(j)}(l)X^{(j)*}(m) varies substantially. As a consequence, the opposite case of inflows surrounded by outflows is also present in the data, which can be seen in Fig. 1.
In addition, there are regions with negative bispectrum and bicoherence. In those regions, the three vectors k, l, m are approximatively parallel, with one component being nearly twice as long as the other two. Figure F.2 shows one example for this in the central feature in the central panel (mR_{⊙} = 310.9, at k = l = m/2). Due to the symmetry relations B(k, l, m) = B(−k, m, l)^{*} = B(m, −k, l)^{*} and as the bispectrum is real in our case, this location has the same bispectral value as a location near the negative peaks in the middle panel in Fig. 3. The pattern associated with these negative values and the associated Fourier components is a pattern made of asymmetric rolls, where the regions of outflow are slightly larger than the regions of inflow. This is consistent with the observation that hexagons with a positive bicoherence also have regions of outflow that are larger in diameter than the inflow regions.
5.2. Generating synthetic data
We generate synthetic data with approximately the given measured power spectrum and bispectrum using the method described in Sect. 4.3. As there is not a unique way to do so, this is only one possible model for the data. Figure 4 shows an example realization of the synthetic LCT data (second row, left column). For comparison, the panel directly below shows an example realization of the LCT data. The synthetic and the real data are very similar and a difference can hardly be seen by eye, apart from the fact that the synthetic data are periodic. Figure F.5 shows the corresponding results for TD, which are similar.
Fig. 4. Example synthetic data (top two rows) generated from the measured bispectrum and example LCT data (bottom two rows). The different columns are the (synthetic or real) data (left), contours of the estimated quadratic nonlinear component overplotted over the data (right), the residual after subtracting the estimated nonlinear component from the data (middle column, two bottom rows), and the Gaussian field from which the synthetic data were generated (middle column, top two rows). Contour lines in the right column indicate levels of ±1σ_{XNL} (dashed) and ±2σ_{XNL} (solid) of the nonlinear component (positive: red, negative: blue). In addition, the skewness is given in the corresponding histogram; both were computed from the entire dataset. See Fig. F.5 for the corresponding TD results. 
The bicoherence and bispectrum estimated from the synthetic data agree well with the estimates from the real data (Fig. 5) with excellent agreement for the bispectrum. The bicoherence of the synthetic data is slightly lower due to the power in the synthetic data being slightly larger than the real data; see Eq. (C.4), and below for further discussion of this topic.
Fig. 5. Comparison of estimated bispectrum (left) and bicoherence (right) for real and synthetic LCT data. The shaded regions indicate 1σ error bars. The LCT data are cuts at k_{x}R_{⊙} = 77.7 through the middle panels of Fig. F.3 (bispectrum) and Fig. 3 (bicoherence). 
As we know the nonlinear and Gaussian components of the synthetic data for every realization (e.g., middle and right panels in the top two rows in Fig. 4), we determine the average fraction of the total power (or equivalently of the total variance) stemming from the nonlinear component; see Appendix D. We find a value of f_{NL} = 4.0%, see righthand column of Table 1.
In addition, we estimate the nonlinear fraction of power from the bicoherence using Eq. (D.6). For the synthetic LCT data, we estimate a value of ; see third to last column in Table 1. This value is roughly consistent with the true value given the bias introduced by the higherorder term in Eq. (D.5). This bias is due to the actual power spectrum of the synthetic data being larger than the real data by a factor of 1 + f_{NL}(k) (see Eq. (C.4)). This reduces the value of the bicoherence by the same amount and results in a small negative bias; see also Fig. 5. It is possible to compensate for this effect in future applications by an iterative approach in which the power spectrum of the Gaussian component is successively adjusted in amplitude.
Under the assumption that this model is a good model for our data, we similarly estimate the fraction of quadratic nonlinear power in the data. We obtain values of 4.2% for LCT, 3.8% for TD, and 5.6% for the smoothed TD data (see also Table 1). However, this result depends on whether the data are filtered in wavenumber space. This is important because both measurements effectively smooth the flow field. To test this, we estimate f_{NL} for different lowpass filters by limiting the sums over m and k in Eq. (D.6) to k,m,m − k≤k_{max}. Figure 6 shows the resulting dependence of on (kR_{⊙})_{max} = k_{max}R_{⊙}. For LCT and the smoothed TD data, the fractional power of the nonlinear component increases until it levels off, while for TD, decreases after attaining its maximum. This is because at large wave numbers, TD is dominated by noise, which was filtered out in the smoothed data.
Fig. 6. Dependence of on lowpass filter applied to LCT and TD data. 
In addition to this, we estimate the quadratic nonlinear component directly from the divergence field using Eq. (E.1). This estimate is biased and may be influenced by higherorder correlations, for example of six Fourier components, which we do not take into account in our study. It therefore serves as an orientation. Figure 4 (right column) shows how these estimates of the nonlinear component relate to the corresponding realizations of the data. The estimates of the nonlinear component for the real data (right column, bottom panel) and for the synthetic data (right column, top panel) are qualitatively similar. Their relationship to the complete synthetic or real data is similar as well. In both cases, outflows (red contours) are seen in the nonlinear component at locations where an outflow is surrounded by a ring of inflows, or when an inflow is surrounded by a ring of outflows. Inflows (blue contours) are seen mostly in the network between supergranular outflows, either amplifying the network inflow, or suppressing outflows that appear in a networklike location (especially visible in the synthetic case).
Taken together, we conclude that at least 4–6% of the variance in the data is due to quadratic nonlinear interactions, under the assumption that the data are well explained by a random Gaussian field plus a weak quadratic nonlinear component.
5.3. Estimating the skewness from the bispectrum
Finally, we obtain estimates of the skewness from the estimated bispectrum using Eq. (9); see the fourth column in Table 1 () and Appendix B for a definition of the estimator. The results are consistent with the estimates from the data, , shown in the third column of Table 1. Similar to the fractional power of the nonlinear component, the skewness depends on any lowpass filter applied to the data; see Fig. 7. The estimated skewness of the filtered dataset is obtained using the estimated bispectrum and by limiting the summation in Eq. (9) to k,l,k + l ≤ k_{max} = (kR_{⊙})_{max}/R_{⊙} and by estimating σ using the similarly filtered power spectrum.
Fig. 7. Dependence of on lowpass filter applied to LCT and TD data. 
For the synthetic data, the two estimates of the skewness from the data and from the bispectrum agree similarly well; see Table 1. However, the skewness in the synthetic data is slightly reduced compared to the real data. This again is due to the total power in the data being larger than the original power by a factor of 1 + f_{NL}(k) due to secondorder effects in the model. This difference could be compensated by an iterative approach in which the power spectrum of the Gaussian component would be adjusted by 1 − f_{NL}(k). The coupling coefficient A_{k, l} would have to be adjusted in this case as well to guarantee that the resulting bispectrum is unaltered.
6. Conclusions
We measured the bispectrum and bispectral coherence of the horizontal divergence of nearsurface supergranular motions from LCT and TD. The measurements for both datasets agree well and yield significantly nonzero values.
The bispectrum is the next highestorder statistic beyond the power spectrum (e.g., Priestley 1981, Sect. 11.5.1); it characterizes how the skewness in a dataset is generated by the quadratic interaction of wave vectors that form a triad k + l = m. The bispectral coherence, which is a powernormalized bispectrum, characterizes the strength of a quadratic nonlinear component in the underlying field as a function of the coupling Fourier components, relative to the total power at this Fourier component.
The bispectral measurements are most significant at supergranular scales, when the three coupling Fourier vectors are of the same magnitude and thus at an angle of about 60° (middle panel in Fig. 3). These Fourier vectors are the same as those of regular hexagons at the supergranular scale. In other words, supergranular scales are the scales where the horizontal divergence pattern looks closest to hexagonal among the scales probed by LCT and TD. At the same Fourier components, the bispectrum is of maximal amplitude and positive, which means that these Fourier components predominantly generate the positive skewness in the dataset (middle panel in Fig. F.3). The positive skewness is due to a preference for broad outflows (positive divergence) surrounded by a weaker network of inflows (negative divergence), which is consistent with the typical appearance of supergranulation (see Fig. 1). The bispectral measurement is consistent with the skewness of the supergranular divergence maps (see Table 1).
We use the measured bispectrum to generate synthetic data which look very much like supergranulation (Figs. 4 and F.5). The synthetic data have the same bispectrum within error bars and a similar skewness (Fig. 5 and Table 1). To do this, we use a model that consists of a Gaussian field with independent Fourier components plus a weakly nonlinear component, which is generated from the Gaussian one by a quadratic coupling equivalent to the measured bispectrum. This method may be extended in the future to generate synthetic supergranular flow fields. These may be used, for example, for testing helioseismic methods.
Based on this weakly nonlinear model, we find that the nonlinear component is strongest at the supergranular scale. We find that at least 46% of the variance is due to this component. This estimate depends on how the smallscale noise is filtered out and it is biased by possibly a few percent (Fig. 6 and Table 1). However, for all lowpass filters tested, the estimates are below 6%.
We find that it is at the supergranular scale that the nonlinear component contributes the largest fraction to the variance. In this sense, the nonlinearities have the strongest effect at supergranular scales. However, the coupling coefficient A_{k, l} (see Eq. (17)), which is the bispectrum weighted with the inverse of the squared power, increases towards smaller scales.
The bispectral measurements presented in this work may be understood in terms of the example patterns introduced in Sect. 3. The skewness and the strength of the nonlinear component in the supergranular divergence maps is of similar order to the pattern termed “small coupling” (see central panel in the bottom row in Fig. 2). This pattern is essentially the same as regular hexagons, but one Fourier component is lowered by a factor of about three.
The relationship between our finding that 4–6% of the variance is explained by the nonlinear component and the conclusion of Rincon et al. (2017) that supergranular dynamics are strongly nonlinear has yet to be determined. Here we present observations of nearsurface flows. In future work, the depthdependence of the skewness and the nonlinear component may be studied and the evolution in the time domain may be included in the analysis.
By measuring the bispectrum, we provide a statistical characterization both of the spatial pattern of supergranulation and of the nonlinear interaction in supergranular flows. Together with future studies, this may result in important insight into how supergranulation emerges as a particular length scale on the Sun because it lends itself as a test that simulations and theoretical models of supergranulation have to pass. Such a comparison can realistically be made because one needs about 17 days of simulated solar time to measure a bispectral peak as found in this study at a signaltonoise ratio of about three.
Quadratic nonlinear interactions play an important part in turbulent convective dynamics. This is because the underlying triadic interactions, which couple two Fourier modes to a third one, are responsible for kinetic energy transport in turbulent flows (e.g., Alexakis & Biferale 2018). This may be used in future studies to estimate the direction of the kinetic energy cascade in supergranular flows and may shed light on the still mysterious origin of supergranulation.
Acknowledgments
We thank Damien Fournier for pointing to a formula that was helpful for the derivations, Friedrich Kupka and Dan Yang for reading and providing feedback on an earlier version of the manuscript, and Michael Wilczek and Friedrich Kupka for helpful discussions. This work was supported in part by the Max Planck Society through a grant on PLATO Science. The computational resources were provided by the German Data Center for SDO through a grant from the German Aerospace Center (DLR). We acknowledge partial support from the European Research Council Synergy Grant WHOLE SUN #810218.
References
 Alexakis, A., & Biferale, L. 2018, Phys. Rep., 767, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A. 2001, ApJ, 550, 824 [NASA ADS] [CrossRef] [Google Scholar]
 Brillinger, D. R. 1965, Ann. Math. Stat., 36, 1351 [CrossRef] [MathSciNet] [Google Scholar]
 Busse, F. H. 2003, Phys. Rev. Lett., 91, 244501 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H. 2004, Chaos, 14, 803 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H. 2007, Sol. Phys., 245, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Chandran, V., & Elgar, S. 1990, IEEE Trans. Acoust. Speech Signal Process., 38, 2181 [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability (Oxford: Clarendon Press) [Google Scholar]
 Contaldi, C. R., & Magueijo, J. 2001, Phys. Rev. D, 63, 103512 [NASA ADS] [CrossRef] [Google Scholar]
 Cossette, J.F., & Rast, M. P. 2016, ApJ, 829, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, Jr., T. L., & Gizon, L. 2000, Sol. Phys., 192, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, Jr., T. L., Jefferies, S. M., Harvey, J. W., & Pomerantz, M. A. 1993, Nature, 362, 430 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, T. L., D’Silva, S., Jefferies, S. M., Harvey, J. W., & Schou, J. 1996, Nature, 379, 235 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, T. L., Hanasoge, S. M., & Chakraborty, S. 2014, Sol. Phys., 289, 3421 [NASA ADS] [CrossRef] [Google Scholar]
 Elgar, S. 1989, Workshop on HigherOrder Spectral Analysis, 206 [CrossRef] [Google Scholar]
 Featherstone, N. A., & Hindman, B. W. 2016, ApJ, 830, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., & Birch, A. C. 2004, ApJ, 614, 472 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., Duvall, T. L., & Schou, J. 2003, Nature, 421, 43 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gough, D. O., Spiegel, E. A., & Toomre, J. 1975, J. Fluid Mech., 68, 695 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proc. Natl. Acad. Sci., 109, 11928 [Google Scholar]
 Itoh, S.I., Itoh, K., Nagashima, Y., & Kosuga, Y. 2017, Plasma Fusion Res., 12, 1101003 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, Y. C., & Powers, E. J. 1979, IEEE Trans. Plasma Sci., 7, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, Y. C., Beall, J. M., Powers, E. J., & Miksad, R. W. 1980, Phys. Fluids, 23, 258 [NASA ADS] [CrossRef] [Google Scholar]
 Komatsu, E., & Spergel, D. N. 2001, Phys. Rev. D, 63, 063002 [NASA ADS] [CrossRef] [Google Scholar]
 Langfellner, J., Gizon, L., & Birch, A. C. 2015, A&A, 581, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Langfellner, J., Birch, A. C., & Gizon, L. 2018, A&A, 617, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moll, R., Pietarila Graham, J., Pratt, J., et al. 2011, ApJ, 736, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Nikias, C. L., & Raghuveer, M. R. 1987, IEEE Proc., 75, 869 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, A., Spruit, H. C., Ludwig, H. G., & Trampedach, R. 1997, A&A, 328, 229 [NASA ADS] [Google Scholar]
 November, L. J., & Simon, G. W. 1988, ApJ, 333, 427 [Google Scholar]
 Palm, E. 1960, J. Fluid Mech., 8, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration XVII. 2016, A&A, 594, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Priestley, M. B. 1981, Spectral analysis and time series (London: Academic Press) [Google Scholar]
 Rincon, F., & Rieutord, M. 2018, Liv. Rev. Sol. Phys., 15, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Rincon, F., Roudier, T., Schekochihin, A. A., & Rieutord, M. 2017, A&A, 599, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Savitzky, A., & Golay, M. J. E. 1964, Anal. Chem., 36, 1627 [NASA ADS] [CrossRef] [Google Scholar]
 Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207 [Google Scholar]
 Schou, J. 2003, ApJ, 596, L259 [NASA ADS] [CrossRef] [Google Scholar]
 Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229 [Google Scholar]
 Simon, G. W., & Leighton, R. B. 1964, ApJ, 140, 1120 [NASA ADS] [CrossRef] [Google Scholar]
 Snodgrass, H. B. 1984, Sol. Phys., 94, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Strugarek, A., Brun, A. S., Mathis, S., & Sarazin, Y. 2013, ApJ, 764, 189 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, J., Gough, D. O., & Spiegel, E. A. 1977, J. Fluid Mech., 79, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Wagner, C., Verde, L., & Boubekeur, L. 2010, JCAP, 2010, 022 [NASA ADS] [CrossRef] [Google Scholar]
 Watkinson, C. A., Majumdar, S., Pritchard, J. R., & Mondal, R. 2017, MNRAS, 472, 2436 [NASA ADS] [CrossRef] [Google Scholar]
 Yeh, T. T., & van Atta, C. W. 1973, J. Fluid Mech., 58, 233 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Fourier transform convention
Here we use the same convention for the discrete Fourier transform as Gizon & Birch (2004),
where h_{x} = L/N_{x} is the sampling interval in one spatial direction, L the spatial extent of the data and N_{x} the number of samples in one spatial direction. Similarly, h_{k} = 2π/L is the resolution in Fourier space.
Appendix B: Bispectral estimators
We here provide the definitions of the estimators used in this work. Table B.1 summarizes the most important quantities and their definitions.
Key quantities and definitions used in this paper.
B.1. Azimuthal averaging
Assuming that supergranulation at the equator does not have a preferred horizontal direction, it is possible to average the bispectral estimates over circles with the same (or similar) values of m. This reduces the noise in the measurement and can be done as follows. For a given triad m = k + l, all three vectors are rotated clockwise by the angle between m and to obtain rotated vectors . The resulting vector then points in the xdirection and the realization X^{(j)}(k)X^{(j)}(l)X^{(j)}(m)^{*} is taken as a realization of in the averaging procedure.
Equivalently, it is possible to simply average the estimates for the bispectrum in the same way. This reduces the bispectral estimates from fourdimensional to threedimensional variables, which can be expressed either as a function of (k, m) or (k,l,m).
We have to take into account that Fourier components of opposite sign have identical data up to complex conjugation and therefore do not add any new information. The number of relizations of X^{(j)} is N = 4015 and the number of realizations used at m is then N_{B} = N_{m}N, where is half of the number of vectors with the same length as m, that is the number of elements of the set
The azimuthal averaging leads to a different number of realizations being used for different values of m so that the errors in the estimates of the bicoherence and bispectrum depend on m. In addition, this number varies slightly with the value of k due to the discretization and due to rounding. For example, at the supergranular peak in the power spectrum, which in our case is at mR_{⊙} = 155.5, we averaged over a typical number of N = 80 300 realizations. This leads to a statistical error for the bispectral coherence of about σ_{b} = 0.0035 at mR_{⊙} = 155.5.
B.2. Definition of bispectral estimators
To shorten notation, we write
Following Kim & Powers (1979), our estimate of the bispectrum, , can then be written
where the vectors are obtained from k, l = k − m, by a rotation of , see Appendix B.1.
Our estimate for the bicoherence then becomes
and similarly for the powernormalized bicoherence,
As a consequence, one has
Appendix C: Details on the weakly nonlinear model
We here show the basic steps for demonstrating that the bispectrum of the synthetic data generated with the weakly nonlinear model of Sect. 4.3 is the same as the input bispectrum to leading order.
To do that, one first observes that different Fourier components of the nonlinear component are uncorrelated,
where δ_{m, m′} is the Kronecker symbol and the nonlinear component is uncorrelated with the Gaussian component,
Consequently, the total power is proportional to
where f_{NL}(k) is the ratio of then nonlinear power with regard to the background power of the Gaussian component at k,
and where O(X_{NL}^{n}) denotes terms that are of nth order in the nonlinear component.
The bispectrum of X is then (see also Wagner et al. 2010)
Using Eq. (17) and using the symmetry relation B_{0}(k, l, m) = B_{0}(k, −m, −l) = B_{0}(l, −m, −k), which is fulfilled by every bispectrum, we then obtain
where the correction term is given by
If the nonlinearity is weak, the term O(X_{NL}^{3}) can be neglected and this model is suited to generate data with a given bispectrum.
If we generate synthetic data for a bispectrum , which was estimated from the data, we set A_{k, l} to zero if the estimated bicoherence is not significantly different from zero. Using too many nonsignificant coupling terms results in additional random nonlinear terms in the data. In this work, we use a threshold for of about 2σ_{b} (1.7σ_{b} for LCT and 2σ_{b} for TD). This threshold was chosen to balance the removal of noise with keeping the signal. This procedure is in general agreement with the discussion in Kim & Powers (1979).
Appendix D: Estimating the power of the nonlinear component
Given the above model, one can estimate the fraction of power due to the nonlinear interaction based on the bispectral estimates. To do that, one considers the power at a Fourier component m,
and finds
where b_{P} is similar to the bicoherence, but normalized by the power spectra (compare to Eq. (12)),
To first order, the fraction of power of the nonlinear component at m is thus proportional to the squared bicoherence, b^{2}, summed over all triads that add up to m (except at k = l). The factor of 1/18 stems from the multiple permutations that are present when computing expectation values of four multiplied Gaussians in Eqs. (C.6) and (D.2).
The fractional variance from the nonlinear component among the total field can similarly be computed,
The fraction of the total variance that is due to the nonlinear component, f_{NL}, is thus a powerweighted average of the fractional nonlinear power at m, plus a small higherorder correction. We thus estimate the fraction of power due to the nonlinear component by
where is given in Appendix B. From Eqs. (D.5) and (D.6), it is clear that is biased due to the higherorder correction in Eq. (D.5). When computing the sum in Eq. (D.6), we again only keep terms where the bicoherence is significantly nonzero by at least 2σ_{b}. This is consistent with using the same threshold for generating synthetic data as discussed at the end of Appendix C.
The bispectrum B and its variancenormalized version B/σ^{3} are effectively bispectral densities (see Eq. (9), which includes an term). As such, they are in principle invariant under changes in the box size, once the resolution of the Fourier grid is much smaller than the scale on which B varies. The same is not true for the bicoherence (Eqs. (D.2) and (D.5) do not include an term).
Appendix E: Estimating the nonlinear component
In addition, it is possible to estimate the nonlinear component of the divergence field at each realization. One starts from Eq. (16) and simply inserts the data X instead of the Gaussian component,
where A_{k, m − k} is obtained from the estimated bispectrum and Eq. (17) with . Again, we keep only terms where the estimated bicoherence is significantly nonzero. The estimator fulfills
We note that is a biased estimator.
Appendix F: Additional figures
To further illustrate our results from Sect. 5, this appendix includes a few additional figures.
F.1. Results for TD
Figure F.1 is essentially the same as Fig. 3 and shows results for the bicoherence b for the (unsmoothed) TD data for selected values of mR_{⊙} that are closest to those in Fig. 3. Results for the smoothed TD dataset are not shown but very similar with slightly less noise and slightly larger values, which is due to the apodization. Both TD results compare well with LCT. Figure F.5 shows the same comparison of synthetic and real data for the unsmoothed TD dataset that Fig. 4 shows for LCT.
F.2. Further details on LCT results
In Fig. F.2, we show results for the bicoherence for more values of mR_{⊙} as shown in Fig. 3. Figure F.3 shows results for the bispectrum. Finally, in Fig. F.4, we show results for the bicoherence as a function of the amplitude of all coupling Fourier vectors, that is of (k,l,m), for the same selection of m as in Fig. F.4.
Fig. F.2. Same as Fig. 3 for more values of mR_{⊙}. We note that the extent of the axis in the top row is different to the other panels. 
Fig. F.4. Estimates of the bicoherence b from LCT as a function of (k,l) for the same values of m as in Fig. F.2, also noted in the top of each panel. White color denotes locations where either the triadic relation k + l = m cannot be met or where one of the three vectors has zero length (as the mean was removed from each map, the bicoherence is not defined in this case). 
All Tables
Summary of results for the bispectral analysis of supergranular divergence maps from LCT and TD.
All Figures
Fig. 1. Example divergence maps (top) and corresponding histograms (bottom) for LCT (left), TD (middle), and smoothed TD data (right). The smoothed TD data was apodized with a raised cosine and the histogram was computed for the region inside the dashed rectangle. The histogram and skewness were computed from the entire dataset. 

In the text 
Fig. 2. Example patterns X(x) and their Fourier transform X(k). The wave vectors in the right small panels are: k_{0} (red), l_{0} (blue), and m_{0} (orange). The skewness γ and fractional power f_{NL} of the nonlinearly coupled component are given in each panel. Equations (3) and (4) give definitions of the patterns. 

In the text 
Fig. 3. Bispectral coherence b(k, m) from LCT for selected values of mR_{⊙} = 77.7, 155.5, 233.2 ± 12.3. The vector mR_{⊙} is fixed for each plot and given by the orange arrow. The bispectral estimates were averaged over the direction of m, which is equivalent to assuming that there is no preferred horizontal direction and which results in a real bispectrum. The direction of m in the plot was chosen simply for convenience. The Fourier vectors for a triad k + l = m that corresponds to hexagons is given by the other two wave vectors (red: k, blue: l). The angles between the displayed Fourier vectors k, m, and l are not exactly 60 degrees because of the discretization. White color denotes locations where no estimate was obtained because one of the three vectors has zero length (as the mean was removed from each map, the bispectral coherence is not defined in this case). See Fig. 5 for a cut through the middle panel including error bars. See Figs. F.2–F.4 for further illustrations of these results and Fig. F.1 for the corresponding TD results. 

In the text 
Fig. 4. Example synthetic data (top two rows) generated from the measured bispectrum and example LCT data (bottom two rows). The different columns are the (synthetic or real) data (left), contours of the estimated quadratic nonlinear component overplotted over the data (right), the residual after subtracting the estimated nonlinear component from the data (middle column, two bottom rows), and the Gaussian field from which the synthetic data were generated (middle column, top two rows). Contour lines in the right column indicate levels of ±1σ_{XNL} (dashed) and ±2σ_{XNL} (solid) of the nonlinear component (positive: red, negative: blue). In addition, the skewness is given in the corresponding histogram; both were computed from the entire dataset. See Fig. F.5 for the corresponding TD results. 

In the text 
Fig. 5. Comparison of estimated bispectrum (left) and bicoherence (right) for real and synthetic LCT data. The shaded regions indicate 1σ error bars. The LCT data are cuts at k_{x}R_{⊙} = 77.7 through the middle panels of Fig. F.3 (bispectrum) and Fig. 3 (bicoherence). 

In the text 
Fig. 6. Dependence of on lowpass filter applied to LCT and TD data. 

In the text 
Fig. 7. Dependence of on lowpass filter applied to LCT and TD data. 

In the text 
Fig. F.1. Same as Fig. 3, but for TD and mR_{⊙} = 73.6, 147.3, 220.9 ± 13.0. 

In the text 
Fig. F.2. Same as Fig. 3 for more values of mR_{⊙}. We note that the extent of the axis in the top row is different to the other panels. 

In the text 
Fig. F.3. Same as Fig. 3, but for the estimated bispectral density . 

In the text 
Fig. F.4. Estimates of the bicoherence b from LCT as a function of (k,l) for the same values of m as in Fig. F.2, also noted in the top of each panel. White color denotes locations where either the triadic relation k + l = m cannot be met or where one of the three vectors has zero length (as the mean was removed from each map, the bicoherence is not defined in this case). 

In the text 
Fig. F.5. Same as Fig. 4, but for TD data. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.