Using Bright Point Shapes to Constrain Wave Heating of the Solar Corona: Predictions for DKIST

Magnetic bright points on the solar photosphere mark the footpoints of kilogauss magnetic flux tubes extending toward the corona. Convective buffeting of these tubes is believed to excite magnetohydrodynamic waves, which can propagate to the corona and deposit heat there. Measuring wave excitation via bright point motion can thus constrain coronal and heliospheric models, and this has been done extensively with centroid tracking, which can estimate kink-mode wave excitation. DKIST is the first telescope to provide well-resolved observations of bright points, allowing shape and size measurements to probe the excitation of other wave modes that have been difficult, if not impossible, to study to date. In this work, we demonstrate a method of automatic bright point tracking that robustly identifies the shapes of bright points, and we develop a technique for interpreting measured bright point shape changes as the driving of a range of thin-tube wave modes. We demonstrate these techniques on a MURaM simulation of DKIST-like resolution. These initial results suggest that modes other than the long-studied kink mode could increase the total available energy budget for wave heating by 50%. Pending observational verification as well as modeling of the propagation and dissipation of these additional wave modes, this could represent a significant increase in the potency of wave-turbulence heating models.


INTRODUCTION
Photospheric bright points are intensity enhancements of order 100 km in size and with lifetimes of order 5 min (Berger et al. 1995;Goode et al. 2010;Liu et al. 2018a), seen in the intergranular lanes of the photosphere, with greater densities along the supergranular network and in plage regions.Bright points usually correspond with kilogauss concentrations of vertical magnetic flux (Utz et al. 2013), sometimes called flux elements, which are believed to be the bases of flux tubes that may rise to the corona, making them the source of much of the open magnetic flux in coronal holes (Hofmeister et al. 2019).The rise of these flux tubes to the low chromosphere has been observed as brightenings in the chromospheric line Ca II H coinciding with photospheric bright points (Xiong et al. 2017;Liu et al. 2018b;Narang et al. 2019).As they rise, these flux tubes expand rapidly (Kuckein 2019) due to the decreasing plasma pressure (Spruit 1976).These magnetic concentrations are responsible for the appearance of bright points, as the magnetic pressure reduces the plasma density, allowing emission from slightly deeper, hotter layers.This produces continuum enhancement and also much greater enhancement in the G-band and other spectral lines (e.g.Uitenbroek & Tritschler 2006).
Bright points are seen to be in nearly-continuous motion, buffeted by the constantly-evolving granular pattern (van Ballegooijen et al. 1998;Nisenson et al. 2003;Utz et al. 2010;Chitta et al. 2012), and this motion is believed to excite flux tube waves which propagate up the tube to the corona, where they can deposit heat through turbulent dissipation (Cranmer & van Balle-gooijen 2005;Soler et al. 2019;Cranmer & Winebarger 2019).Evidence of oscillations in the chromosphere has been observed above magnetic elements, e.g. by Jafarzadeh et al. 2017;Stangalini et al. 2017, supporting this model.This mechanism of wave-driving is proposed as one potential explanation for the driving of high temperatures in the solar corona.A number of studies have used bright-point motions as a proxy for the driving of transverse, kink-mode waves in observations (e.g.Nisenson et al. 2003;Yang et al. 2014) and simulations (e.g.Keys et al. 2011) and have modeled the resulting waves (e.g.Cranmer & van Ballegooijen 2005;De Moortel & Howson 2022), but each has been limited by the barelyresolved nature of these small bright points.Therefore, while kink-mode waves, associated with the overall centroid motion of the bright point, have been studied extensively, other wave modes associated with the shape of the flux-tube cross section are yet under-explored.While efforts have been made to probe wave modes associated with bright-point area changes (e.g.Gao et al. 2021), observations have been very limited in the wave modes than can be studied.
The new Daniel K. Inouye Solar Telescope (DKIST; Rimmele et al. 2020;Rast et al. 2021) is now changing this situation.Its Visible Broadband Imager (VBI; Wöger et al. 2021) achieves the highest spatial resolution ever for solar observations: a diffraction-limited, Gband resolution of 15 km on the photosphere, which can reveal the shapes and sub-structure of the ∼ 100 km bright points.Paired with the instrument's high cadence (a few seconds), it is also possible to probe brightpoint motions at very high frequencies-the frequencies most relevant to wave-turbulence heating models (e.g.Pelouze et al. 2023).Additionally, DKIST's spectropolarimetric instruments, including the Diffraction-Limited Near-Infrared Spectro-Polarimeter (DL-NIRSP; Jaeggli et al. 2022) and the Visible SpectroPolarimeter (ViSP; de Wijn et al. 2022) provide unprecedented resolution and sensitivity, allowing detailed, small-scale probes into the magnetic structure of bright points.
We propose a new technique for connecting observations of bright-point shape evolution to the driving of thin-tube waves in the overlying flux tubes.Our overarching goals in the present work are to produce an initial estimate of the energy flux which could not be measured before DKIST, in order to determine whether further study is warranted (we find that it is), and to propose ideas which may inspire further development.We demonstrate this technique, as well as a bright-point identification algorithm that robustly identifies brightpoint shapes, on a MURaM simulation of DKIST-like resolution, and we use this demonstration to motivate future work.

MURaM Simulation
To demonstrate and analyze the techniques presented in the following sections, we use a MURaM simulation of DKIST-like resolution.Rempel (2014) used an extensively-modified version of MURaM (Vögler et al. 2005;Rempel et al. 2009), a radiative magnetohydrodynamic (MHD) simulation code.We use the outputs of a MURaM run similar to the O16b run of Rempel (2014) but with an expanded vertical domain above the photosphere and with radiative transfer computed in four opacity bins.This run was also analyzed by Agrawal et al. (2018) and the Appendix of Van Kooten & Cranmer (2017).The simulation domain is 6×6×4 Mm 3 at a 16 km grid spacing, with the photosphere approximately 1.7 Mm below the top of the domain.The snapshot cadence is two seconds over about one hour of simulated time.The simulation includes convective flows at well-resolved granular scales, with bright points arising naturally from a magnetic field produced by a smallscale dynamo.The upward-directed, white-light intensity is computed through radiative transfer with four opacity bins and produces an analog of observational images which we use in our analysis.
We show sample data from the simulation in Figure 1.It can be seen that small concentrations of vertical flux abound in the granular downflow lanes, and the strongest ones are associated with intensity enhancements.The flux tube associated with a bright point is also shown.The strength and verticality of the flux tube is greatest near the optical surface.Deeper down, convective forces bend and twist the flux tube out of the plane of the plot.Above the photosphere, the flux tube expands rapidly as the plasma pressure drops, reducing the field strength.

Bright-Point Identification
To identify bright points, we employ a version of the tracking code of Van Kooten & Cranmer (2017), modified by Van Kooten (2021) to achieve robust identification of the edges of bright points (rather than focusing only on centroid location).It consists of the following steps.(Values for all parameters are given in Table 1.) 1.We calculate the intensity of each pixel relative to the mean of its eight neighboring pixels (a discrete analog of the Laplacian), and we identify those pixels for which this value lies n σ standard deviations above the mean of this value (calculated  across that frame).This provides a collection of "seed" pixels which are significantly brighter than their immediate surroundings.This is ideal for bright points, which stand out strongly against the surrounding dark lanes but which are not necessarily the brightest features in the image.
2. Each contiguous set of seed pixels is expanded through n exp iterations to include eligible neighboring pixels.For each set, we identify the maximum brightness B max of the constituent pixels as well as the minimum brightness B min of the pixels in a surrounding box (the minimal rectangle containing the seed pixels, widened in all directions by n exp pixels).In each round of expansion, pixels are eligible to be added if their brightness exceeds B min + f contour (B max − B min ).Effectively, expansion is constrained to be within a contour drawn at that threshold value, which is locally adapted to each bright point's own contrast with its surroundings.We also define a parameter ∆B max : when B max − B min > ∆B max , the contour is instead set at B min + f contour ∆B max .This allows f contour to be kept at an appropriate level for dimmer bright points without constricting the shape of the brightest bright points by placing their contours too high.A fixed number n exp of expansion rounds are conducted, and the neighboring pixels added to the feature in each round are added to the set of seed pixels for the next round of expansion.The number of rounds of expansion is tuned to the typical pixel and bright-point sizes.
3. Following expansion, a "false-positive" rejection step eliminates features that would expand significantly in an additional expansion round.These are typically features which are actually especiallybright portions of a granule, or bright points without sufficiently-clear separation from surrounding granular pixels.This is as opposed to true and clearly-identifiable bright points, which are typically small and compact.We calculate for each feature the fraction of its immediately surrounding pixels which would be added to the feature if another round of expansion were conducted.If that fraction falls above a set threshold f fp , the feature is rejected.
4. Features are also rejected if they touch the edge of the frame, if they are very large or very small (having an area A < A min or > A max , or a minimal bounding rectangle with a diagonal length exceeding d max ), or if they are extremely close to another bright point (within a distance d prox ) and thus possibly part of a larger, complex feature that was mis-identified as multiple fragments.
5. Remaining features are connected from frame to frame based on mutual overlap.Merging and splitting events (i.e., a feature in one frame overlapping multiple features in the preceding or following frame) are treated as the end of the feature(s) going into the event and the beginning of the new feature(s) arising from the event.If, between two frames, a feature's size changes by more than ∆s max,px pixels and that change in area, measured relative to the smaller of the two sizes, is greater than ∆s max,% , the feature is not connected across those two frames.That is, the first feature is ended in the earlier frame, and the second, very differently-sized feature is regarded as beginning in the second frame.This avoids very large impulsive velocities inferred when a sudden, large change in bright-point area occurs, which is usually due to a large number of pixels being very close to a relevant threshold rather than a large change in the underlying intensity distribution.Manual inspection of a sample of such events supports this characterization.A final lifetime criterion eliminates features that are tracked for fewer than t min seconds.
In the present work, we tune many of the thresholds and cutoff values (listed in Table 1) from the values of our past work to better match the tracking to our data set and to improve the subjective quality of the tracking.When considering "contiguous" or "neighboring" pixels, the present work uses only the four above/below/left/right pixels and does not include the four diagonal connections (which requires a large increase in the value of n exp ), to avoid cases in which an identified bright point consists of two separate regions with only a very tenuous, diagonal connection between them..Our most significant change over our past work is the new expansion criterion of step 2. This change significantly improves detection of the edges of each bright point, supporting our current analysis of shapes and sizes.In comparing bright point boundaries drawn with our previous and current approaches, we find from manual inspection that in the majority of cases the two criteria produce comparable boundaries, with each approach at times performing better in the subjective judgement of the viewer.However, we also find a minority of cases in which our new approach produces significantly more plausible boundaries, whereas we did not identify any cases in which the new approach performed significantly worse.

The Identified Boundaries
In Figure 2 we show a sample of identified brightpoint boundaries.Across these samples it can be seen that the boundaries typically align quite well with the boundary that one would draw by eye.In the included animation-which shows a variety of conditions, including splits, mergers, tracking imperfections, and a variety of both accepted and rejected features-it can be seen that the identified outlines very closely follow the evolving bright-point shape and tend to be very stable from frame to frame.Additionally, comparison with the maps of vertical magnetic field strength shows that brightpoint boundaries are typically, but not always, wellaligned with the boundaries of the magnetic enhancement, though it can also be seen that the flux concentration sometimes extends in a weakened form to other pixels not identified as part of the bright point.It can also be seen that weak-field regions not associated with a white-light enhancement are abundant.(These relations between intensity and magnetic flux at the smallest scales are relations that should be verified observationally with DKIST, and not necessarily accepted as fact based on these simulations.It is also important to note that the intensity maps represent a line-of-sight integration, while the field-strength maps show only the value Figure 2. Sample of identified bright points (continued on the next page).These bright points were selected from a sample of 30 randomly-chosen bright points to illustrate different shapes and scenarios.The left column provides an unobstructed view of the intensity pattern, and the center and right columns show bright point boundaries overlaid on intensity and vertical magnetic field strength at τ = 1, respectively.Green lines mark the boundaries of accepted features.Magenta lines (with crosshatches) mark regions rejected by our false-positive rejection step.White lines mark regions rejected due to the minimum-lifetime constraint.Features rejected for other reasons are not present in these samples.Color maps range from 0.7 (black) to 1.75 (white) times the mean intensity for the intensity maps, and -2 (red) to 2 (blue) kG for magnetic field strength.An animated version of the first row is available in the online journal.The 38-second animation shows the evolution of multiple bright points over 382 s of simulation time, at a cadence of 2 s per frame.along the τ = 1 surface, and so in any case differences between the two maps are expected.)

Bright Point Statistics
We identify 1,064 bright points in the MURaM data cube.In Figure 3 we show the distribution of brightpoint sizes, lifetimes, and centroid velocities.We include the distributions using our current tracking approach and, for comparison, distributions produced by applying the tracking algorithm of our past work (Van Kooten & Cranmer 2017).It can be seen that larger bright points are produced by our new approach (driven directly by the increase in A max , d max , and n exp , though the effect of increasing n exp is partially to compensate for disallowing diagonal pixel connectivity).Lifetimes with the new tracking are very slightly skewed toward longer lifetimes, suggesting that the updated feature identification makes the inter-frame linking more robust.Centroid velocities are significantly reduced by the new tracking, which we interpret as being caused by the increased frame-to-frame consistency of the identified bright-point boundaries and therefore a reduction in random centroid "jitter"-this was a strong motivation for the contourbased expansion step.The distribution of per-brightpoint mean vertical magnetic field strengths is shifted toward higher values with our updated tracking.This may be because our new approach is more selective in the features it identifies, and the more marginal features that are now rejected may be less likely to be true or strong-field bright points.Mean values of these measured parameters with the current tracking are 95 km for equivalent diameter, 42 s for lifetime, 855 G for absolute vertical magnetic field strength, and 1.8 km s −1 for centroid velocity.
It can be seen in Figure 3 that the diameter distribution peaks at the smallest allowed size and tapers off with a rather consistent slope (in log space) to ∼140 km, after which the slope flattens until reaching the largest sizes in the detected population.This broken slope suggests our population of bright points may be described well by a two-component model (indeed, we produce a plausible fit with a two-Gaussian model), consistent with the findings of Berrios Saavedra et al. (2022), who present observational evidence for the possibility of a two-population nature for bright points.
Figure 4 shows the distribution of per-pixel values within our identified bright points.It can be seen that these values skew heavily toward strong vertical magnetic field strengths, high intensities, and low depths of continuum formation.(Recall that bright points appear bright because the high magnetic pressure offsets gas pressure and reduces gas density, lowering the τ = 1 surface to deeper, hotter plasma.)Additionally, the clear majority of strong-field pixels are included in identified bright points (though not the majority of flux: only 3.8% of unsigned vertical magnetic flux is within a bright point).Of note, however, these distributions do still include significant numbers of pixels with values that would not be expected for the "ideal" bright point, namely weak fields, low intensity, and high continuum formation.While imperfect tracking no doubt contributes some of these values, others are due to the fact that, while the intensity enhancement and the magnetic flux concentration are related and coupled phenomena, there is no exact relation between the two quantities.With this caveat in mind, the strong-field nature of most detected bright-point pixels suggests that the intensity enhancements do serve as a useful proxy of the fine structure of flux concentrations, even at these small scales.(Observational validation of this claim at the highest resolutions must await DKIST analysis, whereas this has been well demonstrated at medium scales, e.g.Ishikawa et al. 2007.)In Van Kooten (2021) we explore in greater depth the properties of the identified bright points and the impacts our algorithmic changes have on our results.We show that our bright-point boundaries, on an individual basis, are at locations where the optical depth, intensity, and vertical magnetic field strength vary in ways characteristic of the edge of a flux tube.We show that we get largely similar results from applying our tracking code (with few modifications) to photospheric slices of vertical magnetic field strength to track the magnetic flux enhancements that correspond to flux tubes.We also show that our algorithmic changes significantly reduce the amount of centroid "jitter," resulting in power spectra of centroid motion that are much more consistent in slope across a range of frequencies.

Overview
In this section, we present our technique for connecting thin-tube wave modes to bright-point observations.In this method, we "unroll" the outlines of bright points in polar coordinates, which we then fit with a sum of sinusoids, each representing a shape perturbation of a specific order n.Treating changes in these fitted outlines as describing changes in the cross-sectional shape of a flux tube, we connect them to a range of MHD thin-tube wave modes, and we estimate an upward energy flux associated with each.We assume that changes in the bright-point boundary are strictly caused by, and directly connectable to, motions of the flux-tube wall, and that these changes occur at one single height along the strictly-vertical flux tube.These assumptions are not guaranteed to hold (and in fact, it would be astonishing if they were universally true), but they are reasonable extensions of the assumptions in traditional centroid tracking.
We focus solely on the edges of the bright point, which can be measured in broadband observations at the highest possible resolution and cadence.The edges are relatively sharp and easy to track, and we take them as a proxy for the edge of the overlying flux tube.We disregard the intensity pattern within the bright point, which is more difficult to interpret physically.It is very difficult in such an analysis, if not impossible, to fully incorporate the complex dynamics of bright points and the mechanisms by which the visible-light enhancement connected to plasma flows, and we do not claim to do so in the present work.Instead, we emphasize that our approach is a proof-of-concept, to be developed further in future work and with the benefit of spectropolarimetric DKIST observations (which will allow an understand-ing of the degree of correspondence between flux-tubes and bright-point shapes and dynamics at the smallest scales).
We are not the first to consider incorporating the shapes of bright points to more accurately treat their behavior and evolution (see, e.g., Keys et al. 2020, who proposed modeling bright points as ellipses), but we believe that our approach of connecting the shapes of resolved bright points to estimated flux-tube wave energy fluxes provides useful new insights.

Summary of Linear Oscillation Properties
The properties of small first-order wave-like perturbations have been studied extensively, and we will use much of the notation of Spruit (1982) and Edwin & Roberts (1983).We assume a zeroth-order (unperturbed) flux tube that is cylindrical, with radius r 0 , and oriented vertically in the solar photosphere.We will use cylindrical coordinates (r, ϕ, z) to describe po- sitional variations inside (r < r 0 ) and outside (r > r 0 ) the tube.We assume the background magnetic field is always pointed along the z-axis of the cylinder.The pressure, mass density, and field strength are P 0 , ρ 0 , and B 0 inside the tube, and P e , ρ e , and B e outside.We also define the interior and exterior (adiabatic) sound speeds and Alfvén speeds in the standard way, and we further assume that the exterior region is fieldfree, such that B e ≈ V Ae ≈ 0.
The standard derivation for compressive MHD waves begins by assuming that the divergence of the velocity fluctuations behaves as where R(r) is a yet-undetermined function, and the azimuthal mode number n is the primary identifier of the mode of the oscillation.We will solve a dispersion relation for the frequency ω as a function of the axial wavenumber k.Two other key quantities derived by Spruit (1982) and Edwin & Roberts (1983) are the radial wavenumbers for the interior and exterior of the tube: and where the tube speed c T is defined inside and outside the tube as Note that m 2 e is assumed to always be positive, while m 2 0 can be positive or negative.When m 2 0 < 0, we define µ 2 0 = −m 2 0 > 0. The solutions for R(r) take the form of Bessel functions.Both the sign of m 2 0 and the boundary conditions determine which type of Bessel functions to use: where A is a normalization constant with units of s −1 .Spruit (1982) derived the expressions for all three components of the first-order velocity perturbation.Leaving off the exponential term in Equation (3), these are and where the general Bessel function B n must be replaced with a function and argument as in Equation ( 7) above.Note also that the interior or exterior versions of m and c must be substituted as well.
The last key quantity is the dispersion relation for ω(k).It is given by ensuring that both v r and the total pressure are continuous across the boundary at r 0 .Edwin & Roberts (1983) specifies two versions for two distinct cases.For "surface waves" (m 2 0 > 0), and for "body waves" (m 2 0 < 0), where

The Thin Tube Approximation
We will assume that the arguments of the Bessel functions are all ≪ 1 (the thin-tube limit), and we will evaluate the appropriateness of this assumption shortly.
In this limit, the surface-wave and body-wave dispersion relations produce the same expressions for ω 2 /k 2 , and these expressions depend on the value of n.For n = 0 ("sausage-mode" waves), For n ≥ 1 (general "kink-mode" waves), where defines the kink speed c k .The phase velocity is constant, indicating these are dispersionless waves with V gr = V ph .With these phase velocities in hand, we now calculate values for many of these quantities and evaluate the appropriateness of the thin-tube approximation.Drawing from the model of Cranmer & van Ballegooijen (2005), at the photosphere of a typical intergranular bright point surrounded by a field-free granule, we have B 0 = 1430 G, ρ e = 3 × 10 −7 g cm −3 , and a density contrast of ρ 0 /ρ e = 0.316.If the temperature is the same inside and outside the tube, then where we also assume T = 5770 K and the mean atomic weight µ = 1.2 for the mostly neutral photosphere.The gas pressures (P = ρkT /µm H ) inside and outside the tube are roughly 3.765×10 4 and 1.190×10 5 dyne cm −2 , respectively.Thus, c T0 = 6.91 km s −1 , c Te = 0 , c k = 6.42 km s −1 .
The field-free limit for the external medium simplifies the evaluation of m e , since in this case which is indeed always positive for the above values.For n = 0, m e ≈ 0.528k, and for n ≥ 1, m e ≈ 0.614k.The interior radial wavenumber m 0 can be written as and we examine two cases: • For n = 0, with V ph = c T0 , the denominator of Equation ( 17) is zero.However, Edwin & Roberts (1983) show numerical solutions where c 2 T0 < V 2 ph < c 2 0 (with V ph very close to c T0 for small values of the argument to the Bessel functions) producing m 2 0 < 0 and requiring J n Bessel functions.Later, we will show that m 0 (or rather, µ 0 ) is not present in the final equations for the n = 0 modes, so that the resulting very large magnitude for m 0 is not physically relevant.
To evaluate the suitability of the thin-tube limit, we take r 0 = 50 km (for a flux tube of order 100 km across) and V ph = 6.5 km s −1 for the n ≥ 1 modes.Taking m ≈ k, the Bessel-function argument x ≈ kr 0 = ωr 0 /V ph , which depends on the frequency.Representative timescales for bright-point evolution range from about 2 to 20 minutes.This produces a span from x ≈ 0.04 (for the longer period) to x ≈ 0.4 (for the shorter period).For the n = 0 mode, taking V ph = 1.1cT0 , those same representative timescales produce values from x ≈ 0.02 (for the longer period) to x ≈ 0.2 (for the shorter period).We see that the thin-tube approximation is not suitable for any faster flux-tube motions, but it is certainly appropriate for representative motions over the bulk of a bright point's life.

Kink Modes: Velocity Amplitudes and Energy Fluxes
Here we analyze the n ≥ 1 (kink) waves, with n = 0 waves deferred for later.First, we derive the velocity field both inside and outside the tube.In the thin-tube limit, I ′ n (x) ≈ nI n (x)/x, so that Equations ( 9) and ( 10) give v ϕ = iv r .Additionally, the axial (i.e., vertical) velocity amplitude is negligible for a kink-mode wave, since where N 0 = m 0 /k is a dimensionless number.For the example numbers given above, N 0 ≈ 1.448 and |v z /v r | ≈ 5.57kr/n.We are treating n values of order one, and the thin-tube limit sets kr ≪ 1.Thus, we can ignore the contributions of v z motions to the total energy fluxes of modes with n ≥ 1.We aim to derive the full eigenfunctions for v r (r) and v ϕ (r) for different values of n.The constants A given in Equations ( 8) through (10) are not continuous across the boundary of the tube, but v r itself must be continuous.Thus, we define A 0 inside the tube and A e outside the tube, with where, outside the tube, we use K ′ n (x) ≈ −nK n (x)/x.Equating the two velocities at r 0 , and assuming for now that c 0 = c e , we find However, if we choose to normalize everything by the value of v r at the tube boundary (which we label V 0 ), the expression simplifies to Note that the complete expression for v r (r, ϕ, z, t) is the above function multiplied by exp(iωt + inϕ + ikz).It can be seen that the n = 1 mode displays uniform radial motion throughout the body of the flux tube (consistent with its nature of offsetting the entire tube cross section).Higher-n modes are increasingly concentrated at the tube boundary.
As an aside, we showed earlier that inside the flux tube, v ϕ = iv r .However, outside the flux tube, v ϕ = −iv r .Thus, if v r is forced to be continuous at the boundary, then v ϕ must be discontinuous (by changing sign).This was discussed by Goossens et al. (2009).
The instantaneous kinetic energy density carried by a wave with complex amplitudes is defined as where ρ is the zero-order density (ρ 0 inside the tube and ρ e outside the tube), Re(v) is the real part of vector v, and angle brackets denote an average over one wave period.For sinusoidal time dependence, this is equivalent to (see, e.g., Mihalas & Mihalas 1984;Walker 2005).Note that the above expressions for ε are fully spatially dependent, but since we have used Equation (3), the only remaining spatial variation is in the r direction.Thus, if where we define f (r) as the dimensionless quantity to the right of the braces in Equation ( 22), then Traditionally, one multiplies ε by the group velocity to obtain a kinetic energy flux (with units W m −2 ) of However, because each isolated flux tube contains both internal and external fluctuations, we integrate over the horizontal plane to obtain the power (in W) associated with that one tube.Later, we can sum this power over all detected bright points and divide by the total area under study to produce a spatially-averaged energy flux across a patch of the solar surface.We write the power (associated with the kinetic energy) as Note that the two relevant integrals over f 2 are identical: so the power is given by Usually, for linear MHD waves, the kinetic energy is exactly half of the total energy carried by all forms of variability (with the other half being magnetic and thermal energy; see Section 2.5 of Kulsrud 2005), so Ėtot = 2 ĖK .
In any case, if we can measure V 0 for each mode n ≥ 1, we can compute each mode's contribution to the total power (see also Goossens et al. 2013;Van Doorsselaere et al. 2014).

Fluxes
We now consider the separate case of n = 0 modes.The velocity amplitudes are given by Equations ( 8) through ( 10) with appropriate choices of the Bessel functions.Note that v ϕ = 0 for these circularly-symmetric oscillations.In the thin-tube limit, the vertical velocity amplitudes are and the radial velocity amplitudes are Note also that because the explicit dependence of v r on µ 0 (in the r < r 0 case) cancels out in the thin-tube limit, so its value need not be computed.
As with the n ≥ 1 case, we can write the radial velocity in terms of its value at the tube boundary: By requiring v r be continuous at r = r 0 , we can find expressions for A 0 and A e (assuming c 0 = c e as before): Note that in the thin-tube limit, the dominant velocity is the vertical v z inside the flux tube, and this quantity is also spatially constant inside the tube.
To compute the upward kinetic energy flux, we assume that v•v * ≈ |v 2 z | in the interior of the tube, and v•v * ≈ 0 outside, so that and which agrees with Equation (24) of Moreels et al. (2015).
(Recall that ĖK must be doubled to account for all forms of energy.)We discuss next how we infer v z from observations.

Connecting Shape Changes to Waves
To connect measured shape changes to this traditional framework of thin-tube waves, we take our bright-point shapes and "unroll" their outlines into polar coordinates relative to their centroid, producing a function r(ϕ) that describes the outline.We then fit r(ϕ) with a sum of sinusoids where Re indicates the real part of the expression.This represents an unperturbed circle of radius r 0 with radial perturbations for wave modes n = 2 through n max , a chosen parameter.(The n = 1 mode, representing centroid motion, is treated separately, as r(ϕ) is defined relative to the centroid location.)A n and B n represent the fitted parameters (together representing an amplitude and phase).
After fitting each observed shape across a bright point's lifetime, we have time series of r 0 (t), A n (t) and B n (t).Differentiating these time series (and thereby differentiating Equation ( 41)), we obtain measured values

Re
Ḃn + i Ȧn e inϕ . (42) Each term in the summation represents a perturbing wave of the form of Equation ( 25), at some fixed height z (introducing an arbitrary phase offset) and with r fixed at the bright-point radius r 0 and therefore f (r) = f (r 0 ) = 1.The expression Ḃn + i Ȧn therefore fills the role of V 0 , and can be carried through to Equation ( 30).
(As shown earlier, v z is negligible for these modes, and v ϕ = iv r is accounted for in the progression to Equation (30).)The previous derivations presented V 0 as a constant, while here it is a time series of values due to the multi-frequency nature of real oscillations.An RMS of that V 0 series can provide a mean effective amplitude and energy flux, while a power spectrum can probe the frequency content of the observations.Attempting this fitting exposes an immediate problem, that many bright-point shapes do not unroll into single-valued functions r(ϕ) (i.e., the bright-point edge crosses a ϕ value multiple times).To address this, we replace these problematic shapes with approximations that are single-valued functions.We term our chosen solution "shrinkwrapping," in which we draw a circle of points, uniformly-spaced in ϕ and at a large radial distance centered on the bright-point centroid, we draw each point radially in until it first touches the actual bright-point edge, and we treat these points as our approximating shape.Some bright points are problematic in that their centroid lies outside their boundary, and so a suitable stopping condition must be defined as each point is drawn radially in.For points that reach the centroid before touching the bright-point border, we stop the points at r = 0.
Our approach measures perturbations to the assumed equilibrium, circular shape of the bright point and treats the measured centroid of the bright point in any one image as the true, unperturbed centroid.Since n = 1 perturbations (in the small-amplitude limit) represent offsets of the instantaneous centroid relative to the unperturbed centroid, they do not fit in this framework, and so we measure them in the traditional way (see, e.g., Van Kooten & Cranmer 2017) of measuring the velocity of the measured centroid location between successive frames, which serves as V 0 for this mode.(The n = 1 mode is the only mode that perturbs the centroid location.) The n = 0 mode also requires a different approach, as this mode has its dominant velocity in v z , while only v r can be directly inferred from changes in the shape of the flux tube.At the tube boundary (where |v r | = V 0 ), we can use Equations ( 31) and (33) to write (43) The above expression reproduces Equation ( 12) of Defouw (1976) and confirms that in the thin-tube limit (kr 0 ≪ 1), v z > v r .
We write the time series of measured bright-point areas as The constant term πr 2 0 is the mean value of the area A(t).We take the RMS after subtracting the mean (equivalent to taking the standard deviation, which we denote as σ( • )), (using the fact that the standard deviation of a sinusoid is 1/ √ 2 times its amplitude).With this expression, an effective value of V 0 can be computed from a time series of area measurements.V 0 is the magnitude of the radial velocity oscillations at the flux-tube boundary, but we require a value for the dominant vertical velocities.With Equation 43, we find where we have used ω/k = c T0 , a constant.This expression serves as an estimate of the amplitude of v z .An observed sequence of A(t) and v r values will, in fact, be a sum of these sinusoidal oscillations at varying frequencies.Taking advantage of this sinusoidal nature, |v z |/ √ 2 can be taken to be an effective RMS of v z , and it can be used directly in Equation ( 40) to compute a mean upward energy transfer rate over time for the sum of the individual wave components.
If we desire a power spectrum of the energy flux, however, a different approach must be taken.Each observed v r value will be due to the contributions of those multiple sinusoids of varying frequencies, and it derives not from the amplitudes of those sinusoids but from individual values taken from any point along those sinusoidal curves.Additionally, because of the presence of multiple frequency components, the frequency-dependent Equation ( 43) cannot be used to convert these instantaneous v r values into v z values.(It is this requirement to convert the measurable v r values into the dominant v z values that sets the n = 0 mode apart from the n ≥ 1 modes.)Instead, we re-write Equation (44) as which expresses the area of the flux-tube cross section as the area of the unperturbed cross section and an area perturbation due to a small perturbation r 1 (t) in radius.
Taking the time derivative yields where we have used v r (r = r 0 ) = dr 1 (t)/dt, which follows from the radial position of the tube edge, r(t) = r 0 +r 1 (t).Equation ( 48) directly connects the derivative of a time series of A(t) measurements to a time series of v r measurements (assuming a constant value of r 0 , which we take to be the average of r(t) = A(t)/π).
We then take the Fourier transform of the v r time series, which then allows us to re-scale each frequency bin from v r to v z using Equation ( 43), with k = ω/c T0 .We then take an inverse-Fourier transform to recover a time series of estimated v z values.Repeating for each bright point then allows us to produce an average v z (or Ė) power spectrum for the n = 0 mode.

Temporal smoothing
One potential hazard in this approach is that of overdiscretization of the bright point's shape changes by a very rapid cadence.Consider, for example, the case of uniform rotation of an ellipsoidal bright point.While the underlying motion is smooth, pixels must enter or leave the identified bright point at discrete points in time, and this causes a degree of "jumpiness" in the identified outlines of the bright point.To address this, we apply a temporal smoothing to our outline fitting with a window size of 5 frames.For each time step in a bright point's lifetime, we fit Equation ( 42) to a set of r(ϕ) points containing the outline at that time step, with a weight of 1 for each point, and the outlines identified in the two frames before and two frames after, with weights of 2/3 for the nearer frames and 1/3 for the further frames.This ensures that pixels' influence gradually enters and leaves the feature, that the fitted perturbations vary smoothly with time, and that highfrequency noise in the fits is suppressed.Our choice of this smoothing window is discussed further in Section 4.5.

RESULTS
We present here results from applying this technique to bright points identified in simulated observations from the MURaM simulation described in Section 2.1.An archive of our analysis code and data is available (Van Kooten & Cranmer 2023).

Fitted Outlines
In Figure 5 we show a sample of bright-point shapes fit using the method of 3.2, including an illustration of the shrinkwrapping process and the amplitudes of each fitted term.It can be seen that these three samples, which were chosen somewhat arbitrarily, are reasonablywell approximated by the fitted shapes.While we show fitted amplitudes for one point in time for three bright points, it is the time derivative of these amplitudes that determines the wave energy flux.
To understand the impact of the shrinkwrapping process, which is necessary to produce single-valued r(ϕ) fits, we show this procedure's impact in Figure 6.We compute the area contained within our identified brightpoint edges before and after shrinkwrapping, and find that this area change is typically very small (a few percent), and it also follows a more narrow distribution than the distribution of area changes between two consecutive frames, indicating that the impact of shrinkwrapping is smaller than the actual bright-point shape changes that we seek to measure.One could imagine a situation in which shrinkwrapping both adds area to the bright point on one side and subtracts area on another side in roughly equal amounts, conspiring to produce significant shape changes but little total area change.We thus also separately measure the total amount of area added and subtracted from each bright point due to shrinkwrapping.We do this by overlaying the two shape outlines, both before and after wrapping, on a very fine pixel grid and summing the pixels that are contained in one but not the other of the two shapes.We find that area subtractions tend to be very small or non-existent, with area changes almost exclusively being area additions.This is consistent with the general behavior of shrinkwrapping being to remove concave portions of the bright-point shape.
In Figure 7, we show the distribution of the amplitudes of our fitted perturbations.Bright points are confined within the long, thin lanes between granules, which tends to produce bright points with long, thin shapes.This is clearly reflected in the distribution of  41)).The n = 1 mode is not represented on these plots, as it must be defined relative to an unperturbed centroid location, and so it is not a single-frame quantity.The third column shows a particularly interesting case, in which the n = 3 component is dominant, rather than the n = 2 component that dominates most bright points.This example also illustrates well how some shape fidelity can be lost to shrinkwrapping.42)) expressed relative to r0.The n = 1 amplitude is measured as the offset in centroid position between two time steps.On the right are the radial velocity amplitudes (the time derivatives of the radial perturbations).2).The gray bar indicates the total flux of all modes n ̸ = 1, for comparison to the longstudied n = 1 flux.n = 2 perturbation amplitudes (and, to a lesser extent, n = 4), which displays a strong offset from zero as well as much larger typical values compared to the other modes.When taking the time derivative of these perturbations to find the radial velocity of the brightpoint edge (which is the quantity we connect to energy flux), the dominance of the n = 2 mode is diminished and we see the largest values in the n = 1 mode.This leads directly to the dominance of that mode in energy flux, which we show next, which is amplified by a 1/n scaling between velocity amplitudes and energy fluxes.

Energy Fluxes
We described methods for inferring energy fluxes in Section 3.2.3 for n ≥ 1 modes and Section 3.2.4 for n = 0 modes.For each bright point, we use these techniques to produce time series of Ėtot values, indicating the rate of vertical energy transfer through the flux tube, integrated over the area of the bright point (and surrounding areas, for the n ≥ 1 modes where external perturbations are significant).These values can be integrated across the lifetime of each bright point, and then summed across all bright points to produce a total amount of energy transferred across the simulation run-time.We then divide by the simulation duration and area to produce a spatially-and temporally-averaged vertical energy flux across the simulation domain.The mean flux values we compute, as well as mean values for Ėtot , are shown in Table 2 for each wave modes, and the fluxes are also represented in Figure 8.These fluxes have all been scaled by a correction factor as described at the end of this sub-section.
The n = 1 mode is clearly dominant in both measures, as was evident in Figure 7 where the n = 1 mode was dominant in perturbation velocities.But the n = 1 mode is only dominant by a factor of 4 compared to the next-strongest mode (n = 0), and its flux is only about 50% larger than the sum of the n ̸ = 1 fluxes.This indicates that n ̸ = 1 modes (primarily the n = 0, 2 modes) may make a significant contribution to the energy flux budget.We believe these wave modes are therefore a very compelling target for future study, as they have the potential to significantly increase the potency of waveturbulence models of coronal heating.It is important to note, though, that the upward propagation of these modes through the rapidly-varying plasma properties of the chromosphere and transition region must be modeled before any implications to the coronal energy budget become clear.
In computing these fluxes, we have first applied a filtering step to remove bright points that may be problematic.We remove any bright points whose fitted outline r(ϕ) is negative at any ϕ, which is approximately 1% of bright points.These negative values typically indicate a severe violation of the assumption of small perturbations on a cylindrical flux tube.We remove about 8% of bright points for which, at any point in the identified brightpoint lifespan, the shrinkwrapping process changes the total area by more that 15% (i.e.cases where the postshrinkwrapping shape is not a good approximation for the true shape).We finally remove the 22% of bright points which show an area change of more than 35% at any point in their lifetime.These tend to be small bright points, where an otherwise-small area change represents a large percentage difference.Nonetheless, such a large relative change may cause us to infer a large, impulsive wave flux that is more likely a result of the limited resolution than a true effect.28% of all bright points fall into these categories.Since these removed bright points are identified bright points, just with shape histories that are more difficult to treat with this approach, ignoring them can cause an under-estimation of the total waveenergy flux.To account for this, we have divided our reported fluxes by the fraction of bright points that survive these filters-effectively assuming that each ignored bright point has an energy flux equal to the average of the accepted bright points.As these filters primarily remove bright points with extremely large, impulsive changes in area, the removed bright points tend to have correspondingly-large inferred energy fluxes, and so this imputation may represent a conservative underestimate of the energy flux for any incorrectly-rejected bright points.

Flux Comparisons
These flux values can be put in context by comparing them to the upward Poynting flux in the MURaM simulation.We calculated this flux, averaged horizontally and across 19 data cubes evenly spaced through the simulation.We found that this flux is negative immediately below the average τ = 1 surface, rises rapidly to a peak value of 33 kW m −2 at a height of 200 km above the photosphere, and drops to a value of zero at the simulation's upper boundary.We interpret the negative fluxes below the photosphere as being associated with the strong convective downflows in the lanes where bright points exist, as well as downward-propagating waves excited near the photosphere.We interpret the drop-off toward the top of the simulation as damping by the upper boundary condition.We therefore take the peak value of 33 kW m −2 as the most proper comparison point, where the influence of both convective downflows and boundary-condition damping is minimized.When comparing the Poynting flux to the wave energy fluxes above, it is important to remember that the wave fluxes were arrived at by doubling the kinetic energy flux, since MHD-wave energy transport is half kinetic energy and half "other" or potential energy-the latter half includes magnetic and thermal energy (see Kulsrud 2005).Appropriate values for the magnetic energy flux that can be compared with the Poynting flux are therefore at most half of the wave fluxes above (e.g.3.5, 13, and 2.3 kW m −2 for the n = 0, 1, 2 modes, respectively).These values fit together within the Poynting flux budget (which is not expected to consist solely of flux tube waves above bright points), indicating both that the order of magnitude of our fluxes is reasonable, and that the Poynting flux does not show our fluxes to be over-estimates.
In our past work (Van Kooten & Cranmer 2017), we computed an n = 1 flux value from the observational bright-point centroid velocities of Chitta et al. (2012).Disregarding for now the factor used in that work to account for expected wave reflection, and using the ρ and B values of Section 3, these observational velocities produce a comparison flux of 10.7 kW m −2 for the n = 1 mode.Our value is approximately 2.5 times larger, an effect that may be due to increased spatial resolution allowing finer motions to be resolved and more bright points to be identified (see Section 4.6).It is worth comparing these estimated energy fluxes to those required to maintain coronal temperatures.As a simple approach to estimating how all of these wave modes propagate to the corona, we can multiply each estimated flux by (1 − R)/(1 + R), where R is a reflection coefficient.Cranmer & van Ballegooijen (2005) modeled the propagation of the n = 1 mode from the photosphere to the corona, and found typical values of R of ∼ 0.9 below the transition region (with dramatically-reduced reflection above the transition region).Simply applying this factor to all wave modes as a zeroth-order estimation produces energy fluxes at the base of the corona of 360, 1,400, and 240 W m −2 for the n = 0, 1, 2 modes, respectively.These fluxes compare favorably to the classical fluxes required to maintain coronal temperatures of Withbroe & Noyes (1977), i.e., 300 and 800 W m −2 for quiet-Sun and coronal hole regions, emphasizing that these wave modes, as part of a wave-turbulence-based heating model, are plausible as a source of coronal heating.

Spectra
In Figure 9 we plot averaged power spectra of the Ė time series (choosing to show spectra for Ė rather than velocities to account for the mode-dependent relation between the two quantities).We restrict this analysis to only the 58 bright points which we tracked for at least 40 frames, to ensure adequate frequency content in each Figure 10.Power spectra of n = 1 and 2 waves for different temporal smoothing windows.The temporal smoothing window width is given for each spectrum, with weights dropping linearly from 1 at the window center to 0 at the first snapshot outside the window.One time step is 2 s.
sequence as we produce our average spectrum.We employ Welch's method (Welch 1967), dividing these time series into a total of 283 overlapping, Hann-windowed, 40-frame segments (using zero padding at the ends) and we average the spectra of each of these segments.The n ≥ 2 spectra generally have very similar and featureless spectra, while the n = 0, 1 spectra show steepening at higher frequencies.In general, all the spectra show broad frequency content, with no one clearly-dominant frequency range.

Effect of Smoothing
We consider here the effect of the temporal smoothing we employ in this work, implemented via the pixel weights described in Section 3.2.6.In Figure 10 we show power spectra for the n = 1 and 2 modes with various sizes of the temporal smoothing window.The un-smoothed curve for both modes changes to an upward trend at the lowest frequencies, indicative of contamination from high-frequency "jitter" in the brightpoint shape outlines (i.e.jitter due to the impulsive way pixels are added to or removed from identified features, as well occasional pixels that are added to and then rapidly removed (or vice versa) from an identified feature, if those pixel values remain close to a relevant feature-identification threshold).This jitter was a concern discussed in the Appendix of Van Kooten & Cranmer (2017).A smoothing window of 3 time steps (6 s) can be seen to make a large reduction in this effect, while 5 time steps (10 s) is the smallest window that tends to ensure monotonicity.This window size does tend to produce a significant drop in power at the highest frequencies, which may indicate excessive damping.At and above this 5 time-step window size, spectral power begins to decrease in the middle of the spectrum, suggesting that any larger window may over-smooth the bright points at all frequencies.We therefore select a temporal smoothing window of 5 time steps, which seems to be a middle choice that removes high-frequency jitter while erring only modestly toward over-smoothing.

Effect of Resolution and Blurring the Data
A relevant question in this "preparing for DKIST" context is how our analyses change if the MURaM images are degraded to a resolution more comparable to pre-DKIST observations.We investigate this by convolving the MURaM images with an Airy disk of radius 100 km (or 6.25 MURaM pixels).We run our brightpoint tracking on these blurred images, finding through manual inspection of a sample of bright points that no modifications to our algorithm or its parameters are required for satisfactory tracking of the blurred bright points.The number of tracked and accepted bright points is reduced by a factor of 3.3, from 1,064 to 324 bright points.The bright points that are not detected in the blurred data tend to be difficult or impossible to detect by eye in the blurred data.The average brightpoint lifetime increases by a factor of 1.2, from 42 s to 52 s, and the average bright-point area increases by a factor of 1.6, from 36 px to 57 px, suggesting that it is predominantly small and short-lived bright points that are lost.This interpretation is confirmed by the size ratios of individual bright points.The 324 tracked bright points in the blurred data account for 8,261 individual, single-frame outlines.Of these outlines, 7,246 can be unambiguously linked to a single outline in the unblurred tracked results-that is, the blurred feature overlaps ex- actly one unblurred feature, and that unblurred feature overlaps exactly one blurred feature.Among these oneto-one matches, the geometric mean of the ratios of the size of a bright point as detected in the blurred and unblurred tracking is 1.03, indicating that blurring does not typically affect the identified size of a bright point.The increase in average bright-point size can therefore be interpreted as a loss of small bright points.(Only a third of the 324 bright points in the blurred data show a one-to-one match that is consistent across their entire lifetime, so we do not attempt a similar analysis for the changes in lifetimes.) Following the methods described earlier, we compute energy fluxes using our blurred tracking results.This reveals a mode-dependent reduction in flux, shown in Figure 11.The n = 0 flux is reduced from 6.9 to 4.2 kW m −2 , the n = 1 flux from 26 to 6.4 kW m −2 , and the n = 2 flux from 4.5 to 0.67 kW m −2 , with the higher modes all reduced to under 0.2 kW m −2 .This blurring produces an n = 1 flux much more comparable to the value of 10.7 kW m −2 we computed from the observations of Chitta et al. (2012) (discussed in Section 4.3), suggesting that much of the increased n = 1 flux we detect in the MURaM simulations may be due to the presence of smaller bright points that are only detectable in the higher-resolution data, and/or the ability to resolve finer motions.
The flux estimated from the blurred data is approximately 60% of the unblurred n = 0 flux, 25% of the unblurred n = 1 mode, 15% of the unblurred n = 2 mode, and 10% for the higher modes, suggesting that the higher modes are most sensitive to shape information being lost at low resolution, which is to be expected.

DISCUSSION
In addition to applying these techniques to DKIST observations (which we intend to do in future work), At least four factors must be considered before moving forward with the estimated energy fluxes.
First is the "true" filling factor of bright points.Our tracking algorithm does not capture every bright point, and not every strong-field region corresponds to a visible bright point-this can be seen, for example, in Figure 2.This means that we may be underestimating the total wave energy flux.This can be resolved, in part, by further work refining the bright point tracking, but that still will not account for all vertical magnetic flux.Our identified bright points account for 42% of pixels with |B z | > 1000 G, 14% of pixels with |B z | > 500 G, and 4% of all unsigned vertical flux.(This distribution is shown in Figure 4.) Thus, an ad-hoc correction factor of 1/0.42 might be considered to approximately account for the energy flux in strong-field regions not included in the identified bright-point ensemble-regions which might be expected to experience wave driving very similar to the bright points we measure.An alternative way to approach this problem is to track magnetic flux concentrations directly in magnetograms, to apply our techniques to flux-tube cross sections inferred from those measurements, and to consider differences in the fluxes computed from bright-point tracking and flux-element tracking as a rough estimate of the methodological uncertainty present in each flux value, including, in part, the variability due to the difficulty in identifying all flux tubes.
Second, one must assess the degree to which flux tubes reach the corona, as opposed to bending over and reentering the photosphere (e.g., through an oppositepolarity bright point).In regions with mixed magnetic polarity, the latter scenario may be expected to occur with some frequency.Flux tubes with a large vertical extent may therefore be more commonly rooted in large, unipolar network regions.(The MURaM simulation we use does not include the supergranular network, and so this is an effect we are unable to test.)An estimate of the fraction of bright points reaching the corona might be obtained by using a potential field extrapolation from measurements of the average magnetic field within individual bright points (see, e.g., Close et al. 2003;van Ballegooijen & Hasan 2003;Wiegelmann & Solanki 2004).
Third, the degree of propagation of these waves (accounting for both reflection and damping) must be understood before our estimated energy fluxes can be considered in a coronal context.For n = 1 flux-tube waves, Cranmer & van Ballegooijen (2005) modeled the wave propagation from the photosphere, through the heights where the flux tubes are expected to merge into a "canopy," through the rapidly-changing plasma properties of the chromosphere and transition region, and into the corona and heliosphere.They found that approximately 5% of the wave flux reaches beyond the transition region.(Soler et al. 2019 find a similarly strong reduction in flux for the torsional Alfvén wave, a mode we do not consider in the present work.)This factor cannot be directly transferred to the n ̸ = 1 modes, but it is to be expected that these other modes will also have a reduced fraction of their flux succeed in reaching the corona.An increase in the wave flux may also affect the rate at which waves experience reflection or mode conversion at the sharp transition region boundary (see, e.g., Cranmer & Molnar 2023).Thus, for a variety of reasons, further modeling work must be undertaken before the implications of our inferred n ̸ = 1 wave-mode driving can be understood, including the degree to which the energy budget for wave dissipation in the corona is changed by the addition of these wave modes.We hope that the present work provides new motivation for a thorough modeling effort of these modes.
Fourth, this study should be repeated on highresolution observations (which we intend to do with DKIST).While this work provides strong motivation for doing so, the present study only considers simulated images.While observations will of course provide stronger results, they will bring their own difficulties, such as accounting for viewing-angle and projection effects in the observations if far from disk-center, which may complicate both the identification of bright points and the interpretation of their foreshortened shapes; and accounting for seeing effects and ensuring corrections do not significantly damp observed bright-point motions.Such factors will be addressed in future work.
In addition, future work might reassess our theoretical interpretation of the measured shape changes to account for the non-circular shape of bright points.Our approach assumed perturbations on a cylindrical fluxtube cross section, but many bright points have more elliptical shapes due to their confinement in granular downflow lanes.MHD waves in elliptical flux tubes have been considered in some scenarios (Aldhafeeri et al. 2021), and such a treatment could increase the strength of our analysis.Further, not every bright-point shape is close to a circular (or elliptical) shape, meaning the "perturbations" required to describe its shape may not be low-amplitude.(See Figure 7 for the distribution of amplitudes.)This challenges our framework for interpreting those perturbations as waves, and future work should consider whether large-amplitude perturbations should be excluded from the analysis.
One interesting avenue for future exploration is to track flux-tube footpoints by tracking the corresponding magnetic flux enhancement as seen in photospheric magnetograms (i.e., to track the magnetic manifestation of the footpoint rather than the white-light manifestation).This may be expected to more directly probe the cross-sectional shape of the magnetic flux tube.However, magnetograms require spectropolarimetric observations that are much more demanding than the broadband observations that support bright-point measurements.Tracking flux concentrations therefore requires a compromise in cadence, limiting the wave frequencies than can be probed, or in field-of-view, limiting the number of flux concentrations that can be seen.The utility of flux-concentration tracking-even with DKISTmay therefore be limited to more specific case studies where the evolution of a few flux concentrations are compared to that of their corresponding bright points, giving opportunities to constrain and validate the conclusions drawn from bright-point tracking.

CONCLUSIONS
In this paper, we have focused on preparing for DKIST observations of bright points, which can resolve the shapes of bright points at a level never before achieved.These resolved shapes can provide new insights on wave-turbulence heating models; however, doing so requires analysis techniques beyond the centroid tracking that has traditionally been applied to unresolved bright points to study n = 1 waves.We introduced what we believe to be a novel method for analyzing the shape changes of bright points and connecting them to the modeled excitation of flux-tube waves of different modes.We decompose the shapes of bright points into a sum of sinusoids, and we connect variations in these sinusoids to wave-driven perturbations in the cross-section of a thin flux tube.We considered n = 0-10 modes, and our framework can be extended to higher-n modes (though higher modes appear to be unimportant).
Our estimated energy fluxes, when applying this technique to simulated observations from a MURaM simulation of DKIST-like resolution, show that the n = 1 mode is dominant, but the n = 0 and 2 modes make significant contributions to the wave energy budget, together being approximately half the n = 1 flux.This suggests that these wave modes may be significant contributors to wave heating of the coronal in quiet-sun regions, beyond the n = 1 energy flux that has long been the focus of bright-point studies.
These estimated fluxes, though, are to be understood as existing at and immediately above the photosphere.Further modeling must be done to understand if and how these wave modes propagate upward, where they dissipate, and thus whether they make significant contributions to the heating budget in the corona.The results we present now provide a strong motivation for such modeling efforts, and also strongly motivate the application of these techniques to DKIST observations, which we intend to do soon.
DKIST's VBI instrument offers a diffraction limit matching the grid spacing of the MURaM simulation we analyze presently and an even finer pixel scale, which will reduce the impact of the discretization of brightpoint shapes (potentially reducing the necessity of our temporal smoothing step).We thus anticipate that our technique will transfer well to observations and perhaps allow even cleaner inference of wave fluxes.

Figure 1 .
Figure 1.A sample bright point seen in the MURaM simulation.Panel (a) shows the vertical magnetic field strength in a vertical slice through a bright point, with the τ = 1 surface marked by the black line.A strong concentration of vertical flux is seen, coinciding with a lowering of the optical surface by ∼ 100 km.The concentration drops rapidly with height as the flux tube expands.Panels (b) and (c) show the vertical field strength at the average τ = 1 level and the upward, white-light intensity, both for a subsection of the full domain.The granular pattern is seen in the intensity image, and vertical flux concentrations are confined to the downflow lanes.Horizontal lines in these plots indicate the location of the vertical slice in panel (a).

Figure 3 .
Figure3.Distributions of bright point properties using our 2017 tracking and our current tracking.Equivalent diameter is the diameter of a circle of the same area as an identified bright point, and we plot the mean equivalent diameter across the lifetime of each bright point.Mean absolute vertical magnetic field strength is the mean of the absolute value of each pixel across the lifetime of the bright point and is measured at the τ = 1 surface.Centroid velocities are measured and plotted in between each pair of subsequent time steps.

Figure 4 .
Figure 4. Distributions of pixel values for bright points and full images.For each quantity, we show the distributions of the values of each pixel included in a bright point across all bright points and all snapshots ("BP pixels"), the values of each pixel not included in a bright point, and the values of all pixels.Vertical magnetic field strength is measured at the τ = 1 surface.The height of the τ = 1 surface, computed independently at each pixel for a vertical line of sight, is measured relative to the bottom of the simulation box (which has a total height of ∼ 4 Mm).

Figure 5 .
Figure5.A sample of three identified bright points (each of the three columns), illustrating the fitting process.In each column, (a) shows the identified bright point pixels in white, the shrinkwrapped outline as rainbow dots, and the centroid as an "x"; (b) shows the unrolled outline in polar coordinates; (c) shows the fitted outline in orange and the source outlines being fit, with solid blue lines indicating the outline with weight 1, gray dashes indicating the neighboring frames with weights 2/3, and gray dots indicating the more distant neighboring frames with weights 1/3; (d) shows the fitted and source outlines in polar coordinates; and (e) shows the amplitudes of each fitted component (|An + iBn|, using the terminology of Equation (41)).The n = 1 mode is not represented on these plots, as it must be defined relative to an unperturbed centroid location, and so it is not a single-frame quantity.The third column shows a particularly interesting case, in which the n = 3 component is dominant, rather than the n = 2 component that dominates most bright points.This example also illustrates well how some shape fidelity can be lost to shrinkwrapping.

Figure 6 .Figure 7 .
Figure 6.Impact of shrinkwrapping on bright point areas.On the left, we show a 2D histogram showing, for each bright point and each timestep, the change in that bright point's area due to shrinkwrapping and the change in its (un-shrinkwrapped) area over the following timestep.It can be seen that in almost all cases, the area change due to shrinkwrapping is smaller than that which occurs naturally over time.On the right, we show the area of regions added to the bright point and of regions removed from the bright point by the shrinkwrapping process.

Figure 8 .
Figure8.Per-mode energy fluxes integrated over the entire spatial and temporal domain of the simulation (that is, the "mean flux" values of Table2).The gray bar indicates the total flux of all modes n ̸ = 1, for comparison to the longstudied n = 1 flux.

Figure 9 .
Figure 9. Power spectra of the Ė time series for each wave mode.

Figure 11 .
Figure 11.Energy fluxes computed from blurred and unblurred images.

Table 1 .
Van Kooten & Cranmer (2017)Van Kooten & Cranmer (2017)and in the present work.The meaning of each quantity is explained in the text.Dashes indicate quantities not used in the previous work.

Table 2 .
Mean vertical flux and Ėtot values for each wave mode.The fluxes are averaged over the full simulation domain, and the mean Ėtot values are across all bright points and time steps.The large difference in the magnitudes of these quantities is due to the very small filling factor of bright points.The fluxes are corrected as described in Section 4.2.