Characterizing the Motion of Solar Magnetic Bright Points at High Resolution

Rempel (2014) uses MURaM, software for performing radiative MHD simulations, to simulate the upper convective zone and photosphere. The version of MURaM used is extensively modified from the original of Vögler et al. (2005). In this paper, we analyze the results of a run with a domain size of $24.576\times 24.576\times 7.680\,{\mathrm{Mm}}^{3}$, at a horizontal resolution of 16 km per pixel. The data set covers 8280 s, or 2.3 hr, of time at a cadence of one snapshot every 20.7 s (for 401 total snapshots). The simulation run begins with an imposed mean vertical magnetic field of 30 G (as discussed in Section 4.3 of Rempel 2014), but is otherwise identical to the "O16bM" described in the paper's Table 1. The imposed flux is quickly swept into the downflow lanes, resulting in an enhancement of kilogauss features. The top boundary is positioned 1.5 Mm above the $\tau =1$ surface of the photosphere, allowing a convective region depth of 6.2 Mm. From this simulation, we take the $\tau =1$ surface as our solar-observation analog. Within this slice, we use the upward-pointing white-light intensity and the three components of both plasma velocity and magnetic field. A sample patch of the simulation is shown in Figure 1.

Zoom In Zoom Out Reset image size

Table 1.  Properties of the Spectra in Figure 7

Spectrum	Integrated Power (km2 s−2)		Slope at
	over all frequencies	0.005 Hz ≤ f ≤ 0.02 Hz	0.001 Hz	0.01 Hz
Chitta	1.87	0.56	⋯	−1.25
MURaM BPs	3.46	1.07	⋯	−1.27
MURaM Corks	5.42a	1.21	−0.37	−1.59
ROUGH Corks	5.42a	0.57	−0.45	−3.45

Note.

^aROUGH's $V\max$ value is set so its corks' integrated spectrum matches that of the MURaM corks.

Download table as:  ASCIITypeset image

From these MURaM simulations, we extracted granule size distributions, the motion paths of passive tracers, and the motion paths of automatically identified bright points. This section describes the algorithms used for collecting this information, while the data products themselves are analyzed in Section 4.

For the identification of granules, observers such as Abramenko et al. (2012) use intensity from nonspectrographic data to segment granules. With the MURaM data, we segment granules based on vertical velocity instead, using just a positive–negative threshold. While deviating from observational approaches, this produces very similar results with a simple, more direct technique.

We track passive, point-like test particles, called "corks," as one way to analyze the plasma flows in the MURaM data. We begin with a uniformly spaced grid of 4900 corks. (We find the power spectrum of cork motion to be very insensitive to the initial arrangement.) Once initially placed and released, corks move according the horizontal portion of the velocity vector. We linearly interpolate in time to a time step of 4 s when propagating corks, to ensure corks in small, high-plasma velocity regions (e.g., turbulent whirls) do not maintain that high velocity after quickly exiting the region. However, when computing power spectra, we sample the cork velocity every 20 s in order to avoid exceeding the Nyquist rate of the data. With a time step of 4 s, a pixel size of 16 km, and typical plasma velocities of ∼2 km s⁻¹, a cork stays near its initial location during a time step. Cork locations are stored to floating-point precision, and the velocity field is interpolated spatially from the four nearest grid points. Corks are allowed to travel through the periodic boundaries of the simulation. Cork velocities used for computing power spectra are those read out of the interpolated velocity field at cork locations.

For tracking bright points themselves, we use an algorithm that first identifies features within a given frame, following Feng et al. (2012), and then links identified features between frames, following Yang et al. (2014). The intra-frame algorithm keys in on the fact that bright points are small and are characterized by high contrast (they are "bright" relative to their immediate surroundings, but not always in absolute terms—especially in continuum images). In each step, we err toward preventing false positive detections at the cost of more false negatives, since there is such a wealth of bright points available to be detected.

The single-frame stage begins by identifying "seed" pixels in a white-light intensity map. A discrete analog of the Laplacian is calculated by subtracting from each pixel's intensity the mean intensity of the eight surrounding pixels. A local maximum in intensity will produce a significantly positive value, and this indicates a likely bright-point center. We use a threshold of three standard deviations from the mean in the Laplacian to identify a set of likely centers.

Next, the set of candidate pixels is dilated to include pixels that both exceed a lower threshold and neighbor an already-selected pixel (including neighbors along a diagonal). We repeat this dilation process three times, and the threshold for inclusion in this stage is a positive Laplacian value. This process expands the set of selected candidate pixels from the "seed" pixels, which are quite likely to be in bright points, to a nearby cloud of pixels whose likelihood of being in a bright point is enhanced by virtue of being near a pixel already considered a candidate. Since the structure of a bright point tends to be a decrease in intensity from a local maximum to low values at the edge, the use of the Laplacian >0 threshold aims to ensure that the dilation includes pixels within the bright point but stops expanding where it reaches the edge of the bright point. At this point, each contiguous block of selected pixels can be identified and labeled as a unique candidate bright point. An example of such a bright point is shown in Figure 2.

Zoom In Zoom Out Reset image size

The next step is to remove as false positives local intensity maxima inside granules. After the previous three rounds of dilation, these local maxima will produce identified features surrounded by further bright pixels. True bright points, however, tend to be mostly or entirely included within the identified set of pixels at this stage, so the identified feature will be surrounded by dark intergranular lane pixels. (A bright-point diameter is of the order of 100 km, and a pixel in this simulation is 16 km across. Three rounds of dilation allow a maximum expansion of $16\times 3=48$ km or one bright-point radius, meaning a seed pixel at the center can be expanded to cover a typical bright point. Often more than one seed pixel will be identified within a feature, allowing features on the large end of the spectrum to easily be captured.) We compute the fraction of a feature's neighboring pixels, which would be added to the feature if a fourth round of dilation occurred (still using the Laplacian >0 threshold). We reject as false positives any features for which more than 20% of these new neighboring pixels have positive Laplacians.

We apply two more criteria to ensure a high-quality selection of bright points. First, we reject features within four pixels of another feature. This ensures the analysis is of bright points that are cleanly and unambiguously separated from other bright points. We reject as short-lived transients features below four total pixels in area (36 km equivalent diameter). We also reject features above 110 total pixels in area (189 km equivalent diameter) or above 20 pixels (320 km) across (as determined by finding the diagonal of the smallest rectangle to contain the feature), to ensure identified features are not so large or extended that their motion cannot be represented by the motion of their centroid. (This necessarily eliminates from consideration the long chains of connected bright points that can occur in regions of strong background magnetic flux; see Section 5).

This completes the intra-frame identification step. Each frame now has a set of uniquely identified features, but these features must be linked in the temporal dimension for any analysis. This linking step, following Yang et al. (2014), is short: if a feature in frame n overlaps, by at least one pixel, the location of a feature in frame $n-1$, the two features are considered to be the same. This is supported by observations (e.g., Nisenson et al. 2003; Chitta et al. 2012), which find typical bright-point sizes of 100–150 km and typical velocity distributions with maxima around 7 km s⁻¹. In the simulation's 20 s time steps, a bright point would move a maximum of 140 km (with most bright points moving much less), ensuring that most bright points will overlap with their prior position. In situations where a feature in frame n overlaps multiple features in frame $n-1$ (e.g., two bright points merge), or where multiple bright points in frame n overlap one feature in frame $n-1$ (e.g., a bright point splits in two), linkage is ambiguous and is not attempted—instead, the feature(s) in frame n are considered newly found features. Once linking is completed across all frames, we reject bright points with lifetimes under 5 frames (1.6 minutes) as a further guard against transient features.

The result of this algorithm is a set of bright-point identification numbers, with each number corresponding to a set of coordinates within each frame of the bright point's lifespan.

In an attempt to accurately and efficiently capture the laminar portion of the solar granulation pattern, we developed a new Monte Carlo–based model called ROUGH (Random, Observationally motivated, Unphysical, Granulation-based Heliophysics). ROUGH simulations are highly phenomenological, and the properties of the simulation are drawn from observations of solar granulation. The model aims to produce a reasonable representation of granulation at a single slice of constant height while offering very direct control over the properties of that granulation pattern.

In a ROUGH simulation, a set number of granule centers are randomly distributed over a region of space and throughout a given time span. Granule centers are assigned a horizontal drift speed, drawn from a Gaussian distribution of zero mean and width 0.25 km s⁻¹ (Roudier et al. 2012), as well as a drift direction drawn from a uniform distribution. Each center is also assigned a lifetime randomly drawn from the lifetime distribution of Del Moro (2004). We use the exponential fit to lifetimes >5 minutes (${e}^{-t/2.7\min },t\gt 5$ minutes), as this produces excellent agreement between simulated and observed granulation movies. This produces a mean lifetime of 7.7 minutes. At the end of a granule's lifetime, a granule can fade away or it can fragment. Lemmerer et al. (2017) show that the majority of granules fragment into two or more child granules. We give each simulated granule a 70% chance of fragmenting into exactly two children. The child granules each have new lifetimes drawn again from the Del Moro (2004) distribution.

Each granule center is taken as a point source of radially diverging, horizontal flow with a time-dependent amplitude of

$$\begin{eqnarray}&&{V}_{r}(t)=V\max {e}^{-{((t-{t}_{\mathrm{mid}})/\delta )}^{4}},\end{eqnarray} \tag{ 1 }$$

where ${t}_{\mathrm{mid}}$ is the midpoint of the granule center's lifetime and δ, set to a third of the granule's lifetime, determines the falloff rate with time. $V\max$ is set to 5.1 km s⁻¹ to match the integrated ROUGH power spectrum to the integrated MURaM cork-tracking power spectrum. The radial flow is given a spatial dependence of

$$\begin{eqnarray}&&{V}_{r}(r,t)={V}_{r}(t)\,\displaystyle \frac{r}{{r}_{0}}\,\displaystyle \frac{1}{{(1+{(r/{r}_{0})}^{a})}^{b}},\end{eqnarray} \tag{ 2 }$$

where ${r}_{0}=500\,\mathrm{km},a=4$ and $b=1/a$ determine the shape of the function, and r is the horizontal distance from a location to the center of the granule. This functional form is chosen to match the observed granulation pattern (see also Simon & Weiss 1989). At each pixel in the simulation, the horizontal velocity is taken to be the velocity corresponding to the largest value of ${V}_{r}(r,t)/r$ produced at that pixel by any one granule center. In other words, each granule center has a "sphere of influence" in which it is both the sole producer and the strongest possible producer of horizontal flow. ${V}_{r}/r$ is used for this step rather than ${V}_{r}$ itself to flatten the initial ramp-up from ${V}_{r}(r=0)=0$ to the function's maximum value.

At this point, the modeled granulation pattern will have discontinuous transitions at granule boundaries. Gaussian smoothing with a $1/e$ full-width of 317 km is applied to the horizontal velocities to produce realistic inter-granular lanes of finite size. An example of the resulting granulation pattern is shown in Figure 3.

Zoom In Zoom Out Reset image size

As a final step, vertical velocities are calculated by mass conservation in an assumed hydrostatic density profile, using

$$\begin{eqnarray}&&{V}_{z}={H}_{\mathrm{eff}}\left(\displaystyle \frac{{{dV}}_{x}}{{dx}}+\displaystyle \frac{{{dV}}_{y}}{{dy}}\right)\end{eqnarray} \tag{ 3 }$$

(see Simon & Weiss 1989), where the effective scale height ${H}_{\mathrm{eff}}$, defined in terms of the pressure and vertical velocity scale heights by $1/{H}_{\mathrm{eff}}=1/{H}_{P}+1/{H}_{{V}_{z}}$, is set to 75 km computed from values of ${H}_{{V}_{z}}\,\approx$ 150 km (Oba et al. 2017) and ${H}_{P}\,\approx$ 150 km.

ROUGH cannot contain explicit bright points, so instead we analyze the horizontal flow fields using passive tracers ("corks"). A cork is placed initially at the center of the simulation domain. At each time step, the horizontal velocity at the cork's location is calculated by a two-dimensional interpolation from the four nearest pixels, and the cork moves with that velocity until the next time step. The cork's location is not constrained to discrete pixel locations, but is calculated and saved to full numeric precision. When calculating power spectra, we remove the initial portion of each cork's path, up until the cork first reaches a pixel with a downward vertical velocity. This is more representative of actual bright points, though we find that doing so produces little change in the final power spectrum. An example cork motion track is shown in Figure 3.

The ROUGH model contains a number of adjustable parameters that control the granulation pattern, and varying those parameters yields physical insight. We find that the cork power spectra are insensitive to the granule drift velocity distribution, consistent with the fact that granule drift is slow compared to granular flow velocities and to the timescale of change in the granulation pattern. We find increased power in cork motion when the horizontal plasma velocities are increased, with the integrated power spectrum following a $\propto {V}_{\max }^{2}$ scaling for moderate changes in velocity (reflecting the nature of kinetic energy). We also find that total power is decreased when the granule emergence rate (per unit area, per unit time) is increased or when the granule scale size r₀ (Equation (2)) is increased. Increasing the emergence rate increases the number of granules and thus the fraction of area occupied by stationary inter-granular lanes. Lanes are areas of converging horizontal flow, so a cork in a lane is in a relatively stable position and will experience less motion. Furthermore, having more granules causes each individual granule to be smaller. When a new granule emerges and rearranges the local intergranular lanes, then, the lanes will be shifted by a smaller distance and any displaced corks will be dragged along for a shorter distance before being in a stable location again. An increase in the granule scale size pushes the maximum of the radial velocity function into the region affected by the Gaussian smoothing, thus extending the low-horizontal-velocity region from the center of granules to the edges. This means that when the locations of the lanes changes, the granule-edge, horizontal flows that push corks along with the moving lane are weaker, and so corks move more slowly.

The analysis in this paper is of 3072 independent ROUGH simulation runs, with one cork per run. We use a spatial resolution of 50 km/pixel, a time step of 2 s, a simulation duration of 167 minutes, and a granule emergence rate of $5.9\times {10}^{-9}$ km⁻² s⁻¹.

Here we present a comparison of the granulation patterns in ROUGH and MURaM as a validation of the ROUGH approach. In Figure 4, we compare the distribution of granule sizes between the two simulations and the observational values. The observational values, from Abramenko et al. (2012), draw on a field of view comparable to the MURaM simulation domain and a diffraction-limited resolution of 77 km. Abramenko et al. (2012) use intensity thresholding to separate granules from intergranular lanes. As described in Section 2.2, we instead use a ${V}_{z}=0$ threshold, which provides more straightforward granule segmentation.

Zoom In Zoom Out Reset image size

Observed granule sizes are fit by the sum of a Gaussian distribution of large granules and a power-law distribution of small granules (Abramenko et al. 2012), thought to reveal the presence of two distinct populations of granules (as in Hirzberger et al. 1997). We find that the MURaM distribution reproduces very well the observational distribution. However, the ROUGH distribution contains only the Gaussian portion. This suggests that the small granules with a power-law size distribution have their origins in the turbulent phenomena not included in the ROUGH model, whereas the Gaussian-distributed large granules are produced by the large-scale, convective structures of the photosphere.

Figure 4 also compares the distribution of horizontal velocities in the two simulations. While the shape of the ROUGH distribution is set by the radial velocity functions (Equations (1) and (2)), the horizontal scaling is set by the matching of integrated power spectra. This produces good agreement between the two horizontal velocity distributions. Notable, however, is the hard cutoff at $V\max$ in the ROUGH plot, whereas the MURaM plot has a high-velocity tail due to turbulent regions. Figure 4 finally compares two granulation images, mapped in terms of vertical velocity, in which the qualitative agreement is good (aside from the presence of turbulence).

Applying the bright-point identification algorithm described in Section 2.2, we find a total of 3185 bright points meeting all our size and lifetime thresholds in the MURaM simulation data over the simulation time of 2.3 hr (401 snapshots). Of these thresholds, the required lifetime of at least 1.6 minutes has the largest effect. 35 out of every 36 identified bright points are removed by this criterion. (See the steep lifetime distribution in Figure 6.) Of the remaining 3185 bright points, nine were the result of mergers, 34 ended in a merger, 37 were the product of a fragmenting bright point, and 56 ended by fragmenting.

In Figure 5, we show the distribution of observed horizontal velocities for bright points, corks, and the plasma flow. Cork velocities are lower than plasma velocities. This is consistent with the fact that corks in a high-velocity region will travel at high velocity and quickly leave, while corks in a low-velocity region will linger; low-velocity regions are therefore preferentially sampled by corks. Bright-point velocities may be susceptible to contributions from jitter in the centroid location caused by, e.g., bright-point shape changes or variation in the identified edge of the bright point. We note here that bright-point velocities are typically lower than plasma or cork velocities, suggesting that this is not a strong effect. In the Appendix, we show more evidence that this effect appears to be insignificant in our data.

Zoom In Zoom Out Reset image size

In Figure 6, we show selected distributions of bright-point lifetimes, intensities, sizes, and vertical magnetic field strengths. The distributions of bright-point size, vertical flux, and intensity all cluster near a central value (respectively, 109 km, 1.25 kG, and 2.7 I₀). The bright-point lifetime, however, shows a sharp decay from the minimum lifetime threshold (1.6 minutes), with nearly all bright points falling below a 5-minute lifetime. The mean bright-point lifetime (for bright points living longer than 1.6 minutes) is 2.7 minutes, while the maximum is 14.5 minutes. This mean is comparable to the mean lifetime of ∼3 minutes for the longer-lived granule population in Del Moro (2004), indicating that bright points evolve on the same timescale as the granulation pattern. Our mean bright-point lifetime is notably shorter than the 9.33 minutes of Berger et al. (1998). However, they take steps to track bright points through merging and splitting as well as through short (∼40 s) disappearances of bright points. We do not attempt to match these capabilities, as doing so would complicate the calculation of power spectra.

Zoom In Zoom Out Reset image size

Neither bright-point lifetime nor bright-point velocity show a strong correlation with any other measured value. This may be expected for lifetime due to our inability to track bright points through splitting and merging events, meaning that our reported lifetimes do not all correspond to the true lifetimes of bright points, and due to the fact that bright points themselves are only one stage in the constant evolution of magnetic flux. The noncorrelation may also be expected for velocity, since bright-point motion is controlled by external forces. However, strong correlations emerge between bright-point intensity, size, and magnetic field strength. We use for each bright point the maximum size and the 95th percentile pixel value value for intensity and field strength in order to capture the value when the bright point is most fully emerged and vertical but to avoid outlier values. We find a correlation coefficient of 0.59 between intensity and $| {B}_{z}|$, 0.44 between size and $| {B}_{z}|$, and 0.62 between size and intensity. These correlations are all statistically significant and reasonably suggest that a stronger flux concentration will produce a larger, brighter bright point (a conclusion supported observationally, e.g., Ji et al. 2016).

The correlations we find between bright-point lifetime and both diameter and velocity differ from those found in Yang et al. (2014). These Hinode/SOT observations yielded a correlation coefficient between lifetime and size to be +0.83 (contrasted with our value of +0.01 for mean values and +0.1 for peak values) and between lifetime and velocity to be −0.49 (contrasted with our −0.09 for mean values and +0.14 for peak values). The large differences here are likely due to a difference in approach. Bright-point lifetimes are discrete multiples of the imaging cadence, and for each lifetime value, Yang et al. (2014) calculate a mean value of diameter and velocity among all bright points of that lifetime and use these mean values to perform a linear regression. In contrast, we use our full data set with all its variation. Using average values collapses the data to a mean trend line, which risks artificially enhancing the resulting correlation coefficient. We believe that using the full set of data more accurately represents the level of variation in the data (and, therefore, how well one measurable can predict another).

Following the literature (van Ballegooijen et al. 1998; Cranmer & van Ballegooijen 2005; Chitta et al. 2012), we compute power spectra from velocity sequences via the Wiener–Khinchin theorem, which provides a power spectrum as the Fourier transform of the autocorrelation of a data sequence. While perhaps a roundabout approach,¹ this allows individual bright-point velocity sequences to be more straightforwardly averaged together. Since velocity sequences of different lengths are produced, their direct Fourier transforms (and those transforms' absolute-square) have different ranges and spacings in frequency, making a more direct averaging nontrivial. However, as long as different velocity sequences have the same temporal spacing, their autocorrelations will have the same spacing. This allows averaging to be performed on the autocorrelations (analytically equivalent to averaging the final power spectra), which is computationally simpler.

For the MURaM bright points, we calculate power spectra of the finite-difference velocities of the intensity-weighted centroids of bright points. Our observational comparison, Chitta et al. (2012), present their power spectrum and autocorrelation function as fitted, analytical functions. In this section, we present the power spectrum produced by discretely sampling the Chitta et al. (2012) autocorrelation function at its observational cadence of 5 s out to ± 105 s (the range used for their fit to the autocorrelation) and Fourier transforming, to provide a more direct comparison to our own power spectra.

Power spectra for observations, MURaM corks, MURaM bright points, and ROUGH corks are shown in Figure 7. ROUGH is configured so that the integrated power matches the integrated MURaM cork power.

Zoom In Zoom Out Reset image size

At low frequencies, the ROUGH spectrum is similar to the MURaM cork spectrum, whereas at high frequencies it drops off precipitously. This reflects the similarity between the two simulations of the large-scale granulation patterns, which cause low-frequency movement, and the absence in ROUGH of small-scale turbulence (e.g., vortex flows, see Giagkiozis et al. 2017), which drive high-frequency motion. This comparison also suggests that these two phenomena are the dominant drivers of motion in their respective frequency regimes.

The MURaM cork spectrum shows more low-frequency power than the MURaM bright points, suggesting that perhaps corks are more free to drift in steady flows than are bright points, but the two spectra converge at high frequencies, suggesting that they respond similarly to rapid variation in plasma flow. Finally, both MURaM analyses show a consistent enhancement in power over the Chitta data for the frequencies of overlap, quantified later in this section.

The integral of the power spectrum is equal to the variance of the velocity data (via Parseval's theorem). These integrated values are shown in Table 1. We report both the integral of each spectrum over its full extent and the integral over the frequency range covered by all four spectra (0.005 Hz ≤ f ≤ 0.02 Hz). The former represents a spectrum's lower bound on the "true" total power, unconstrained by limited frequency coverage, while the latter allows spectra to be compared over a comparable domain. These numbers reveal more motion and variation for MURaM corks than for bright points, and as noted above, we find a 91% increase in total power measured in the frequencies of overlap for MURaM bright points versus observed bright points. This increase in motion, in both the power spectra and the integrated spectra, might be explained in part by the calibration offset applied by Chitta et al. (2012), a step that could tend to damp out some real motion in addition to spurious motion and that is not needed in simulation analysis. It might also suggest that there is, in fact, more power to be detected with higher-quality observations. Improvements in spatial resolution and sensitivity will make bright-point detection and centroid identification more reliable and make smaller bright points visible, both enabling more robust statistics. Verification of the presence of unmeasured power must wait until DKIST comes online, currently scheduled for 2019. The Visible Broadband Imager on DKIST will have a diffraction-limited, G-band resolution of ∼16 km (comparable to these MURaM simulations) and a 3 s cadence (Elmore et al. 2014), improving on both the ∼160 km resolution of the Hinode SOT observations of Yang et al. (2014, 2015) and the ∼100 km resolution and 5 s cadence of the SST observations of Chitta et al. (2012).

Both the MURaM and Chitta et al. (2012) power spectra shown in Figure 7 indicate power-law behavior with slopes of the order of ${f}^{-1}$ between 10⁻³ and 10⁻² Hz (see Table 1). This is a much flatter spectrum than the presumed exponential drop-off associated with random-walk-like motions such as those found from earlier measurements of bright points (e.g., van Ballegooijen et al. 1998; Cranmer & van Ballegooijen 2005). However, whether this points to the presence of an active turbulent cascade (Petrovay 2001) in the photosphere is still not known. Reduced MHD turbulence simulations for solar flux tubes (van Ballegooijen et al. 2011; Woolsey & Cranmer 2015) showed that low-frequency random-walk motions rapidly produce a steep power-law spectrum ($P\sim {f}^{-4.5}$) in the photosphere. This may be modified for the flatter bright-point driving spectrum found in this work. There is a long history of interplanetary ${f}^{-1}$ fluctuations being interpreted as fossil remnants of oscillations that map down to the solar surface (Matthaeus & Goldstein 1986; Velli et al. 1989; Dmitruk & Matthaeus 2007; Verdini et al. 2012). However, those fluctuations occur at even lower frequencies than those studied in this paper (e.g., $f\lt {10}^{-4}$ Hz) and their physical origin is still being debated.