SSubmitted to Ap. J.
Tracking Vector Magnetograms with the Magnetic Induction Equation
P. W. Schuck † Plasma Physics Division, United States Naval Research Laboratory4555 Overlook Ave., SW, Washington, DC 20375-5346
ABSTRACT
The differential affine velocity estimator (DAVE) developed in Schuck (2006) for estimatingvelocities from line-of-sight magnetograms is modified to directly incorporate horizontal magneticfields to produce a differential affine velocity estimator for vector magnetograms (DAVE4VM).The DAVE4VM’s performance is demonstrated on the synthetic data from the anelastic pseu-dospectral ANMHD simulations that were used in the recent comparison of velocity inversiontechniques by Welsch et al. (2007). The DAVE4VM predicts roughly 95% of the helicity rate and75% of the power transmitted through the simulation slice. Inter-comparison between DAVE4VMand DAVE and further analysis of the DAVE method demonstrates that line-of-sight trackingmethods capture the shearing motion of magnetic footpoints but are insensitive to flux emer-gence — the velocities determined from line-of-sight methods are more consistent with horizontalplasma velocities than with flux transport velocities. These results suggest that previous studiesthat rely on velocities determined from line-of-sight methods such as the DAVE or local cor-relation tracking may substantially misrepresent the total helicity rates and power through thephotosphere.
Subject headings: magnetic fields — Sun: atmospheric motions — methods: data analysis
1. Introduction
Coronal mass ejections (CMEs) are now recognized as the primary solar driver of geomagnetic storms(Gosling 1993). Several theoretical mechanisms have been proposed as drivers of CMEs, including largescale coronal reconnection (Sweet 1958; Parker 1957; Antiochos et al. 1999), emerging flux cancellation ofthe overlying coronal field (Linker et al. 2001), flux injection (Chen 1989, 1996), the kink instability offilaments (Rust & Kumar 1996; T¨or¨ok et al. 2004; Kliem et al. 2004), and photospheric footpoint shearing(Amari et al. 2000, 2003a,b; Schrijver et al. 2005). All of these CME mechanisms are driven by magneticforces. The main differences depend on whether the magnetic helicity and energy are first stored in thecorona and later released by reconnection and instability or whether the helicity and Poynting fluxes areroughly concomitant with the eruption. The timing and magnitude of the transport of magnetic helicity andenergy through the photosphere provides an important discriminator between the mechanisms. In addition,eruption precursors in the photospheric magnetic field might provide reliable forecasting for space weather † [email protected] a r X i v : . [ a s t r o - ph ] J u l ∂ t B z + ∇ h · ( B z V h − V z B h ) = 0 , (1)where the plasma velocity V and the magnetic fields B are decomposed into a local right-handed Cartesiancoordinate system with vertical direction along the z -axis and the horizontal plane, denoted generically bythe subscript “h,” containing the x - and y -axes.D´emoulin & Berger (2003) observed that the geometry of magnetic fields embedded in the photosphereimplied that F = U B z ≡ B z V h − V z B h = (cid:98) z × ( V × B ) = (cid:98) z × ( V ⊥ × B ) , (2a)where F denotes the flux transport vector, U is the horizontal footpoint velocity or flux transport velocity( U · (cid:98) z = 0) and V ⊥ is the plasma velocity perpendicular to the magnetic field V ⊥ · B = 0. The flux transportvectors are composed of two terms B z V h and V z B h representing shearing due to horizontal motion and fluxemergence due to vertical motion respectively. Equation (2a) may be used to transform (1) into a continuityequation for the vertical magnetic field ∂ t B z + ∇ h · ( U B z ) = 0 , (2b)where plasma velocity may be written generally in terms of the flux transport velocity as V = U − ( U · B h ) B | B | + V (cid:107) B | B | , (3a) V ⊥ h = U − ( U · B h ) B h | B | , (3b) V ⊥ z = − ( U · B h ) B z | B | , (3c)and the subscripts “ (cid:107) ” and “ ⊥ ” denote plasma velocities parallel and perpendicular to the magnetic fieldrespectively. Equations (3a-c) are the algebraic decomposition (Welsch et al. 2004) generalized for arbitraryparallel velocity V (cid:107) , but the value of V (cid:107) does not affect the perpendicular plasma velocity (3b)-(3c) or theperpendicular electric field c E ⊥ = − V × B = − E ⊥ h (cid:122) (cid:125)(cid:124) (cid:123) U × (cid:98) z B z − E ⊥ z (cid:122) (cid:125)(cid:124) (cid:123) U × B h , (4)which both depend only on the flux transport velocity U . 3 –Equations (1)-(3) should be formally distinguished from the inverse problem for determining an estimateof the plasma velocity v from vector magnetograms using the normal component of the magnetic inductionequation ∂ t B z + ∇ h · ( B z v h − v z B h ) = 0 , (5a)where f = u B z = B z v h − v z B h = (cid:98) z × ( v × B ) = (cid:98) z × ( v ⊥ × B ) , (5b)and the inverse problem for determining flux transport velocity u from the evolution of the vertical magneticfield or line-of-sight component ∂ t B z + ∇ h · ( ϑ B z ) = 0 . (6)The notation ϑ , denoting an optical flow estimate, emphasizes that ϑ determined from (6) is not necessarilyimmediately identified with the flux transport velocity u . Equations (5)-(6) are ill-posed inverse problemsbecause of two ambiguities:1. The Helmholtz decomposition of the flux transport vectors (Welsch et al. 2004; Longcope 2004) f = u B z = B z v h − v z B h = − ( ∇ h φ + ∇ h ψ × (cid:98) z ) , (7)where φ is the inductive potential and ψ is the electrostatic potential manifestly demonstrates thatonly inductive potential ψ may be unambiguously determined from the local evolution of B z in (5).The electrostatic potential φ must be constrained by additional assumptions. By analogy, (6) is alsoill-posed for the same reason; ∇ h × ( ϑ B z ) is not constrained by the local evolution of ∂ t B z .2. For (5a) and (5b) V (cid:107) is not constrained by the local evolution of ∂ t B z . For (6), there is no a priori relationship between ϑ and u or ϑ and v for the inverse problem. However, if ϑ is identified with theflux transport velocity u then u and v will satisfy the same relationships as U and V in (3).The first ambiguity may be resolved for (5a) by the induction method (IM) (Kusano et al. 2002, 2004),minimum energy fit (MEF) (Longcope 2004), or the differential affine velocity estimator for vector magne-tograms (DAVE4VM) presented in §
2. These methods produce a unique solution for v , but not necessarily the unique solution that corresponds to V . The first ambiguity may be resolved for (6) by local optical flowmethods such as the differential affine velocity estimator (DAVE) (Schuck 2006), its nonlinear generalization(Schuck 2005), global methods (Wildes et al. 2000), the minimum structure reconstruction (MSR) (Geor-goulis & LaBonte 2006) which imposes v ⊥ z = 0 as an assumption, or hybrid local-global methods such asinductive local correlation tracking (ILCT) (Welsch et al. 2004). These methods produce a unique solution for ϑ , but not necessarily the unique solution that corresponds to u .Several assumptions have been used either explicitly or implicitly to resolve the second ambiguity. Chaeet al. (2001) conjecture that local correlation tracking (LCT) (Leese et al. 1970, 1971; November & Simon1988) provides a direct estimate of the horizontal photospheric plasma velocity: ϑ (LCT) = V h . D´emoulin& Berger (2003) conjecture that line-of-sight tracking methods, and in particular LCT, estimate the totalflux transport velocity ϑ (LCT) = U . Schuck (2005) formally demonstrated that LCT is consistent with theadvection equation ∂ t B z + ϑ (LCT) · ∇ h B z = 0 , (8)not the continuity equation in (6), but that LCT could be modified to be consistent with (6) by directintegration along Lagrangian trajectories in an affine velocity profile. Nonetheless, both conjectures may be 4 –considered in the context of (6). Under Chae et al.’s (2001) assumption, the flux transport velocity wouldbe derived from line-of-sight optical flow methods via u B z ≡ ϑ B z − v z B h , (9)where in principle, v z might be approximately determined from Doppler velocities near disk center. UnderD´emoulin & Berger’s (2003) assumption, the total flux transport velocity would be derived from line-of-sightoptical flow methods via u ≡ ϑ for B z (cid:54) = 0 . (10)The Ansatz u = ϑ has important implications for solar observations. This conjecture implies that the totalhelicity and Poynting flux may be estimated by tracking the vertical magnetic field or by tracking the line-of-sight component near disk center as a proxy for the vertical magnetic field. D´emoulin & Berger’s (2003) Ansatz has largely been accepted by the solar community (Welsch et al. 2004; Welsch et al. 2007; Kusanoet al. 2004; Schuck 2005, 2006; LaBonte et al. 2007; Santos & Bchner 2007; Tian & Alexander 2008; Zhanget al. 2008; Wang et al. 2008). However, equivalence between ϑ and u for line-of-sight methods has neverbeen practically established. These two different hypotheses (9) and (10) for the interpretation of ϑ inferredby DAVE will be considered in § v (cid:107) = 0 (Longcope2004). ILCT and the original algebraic decomposition both assume v (cid:107) = 0 (Welsch et al. 2004). Georgoulis& LaBonte (2006) describe a method for inferring v (cid:107) from Doppler measurements for MSR. For DAVE4VMthe second ambiguity is resolved simultaneously with the first. The DAVE4VM method estimates a fieldaligned plasma velocity from only magnetic field observations!Using established computer vision techniques (Lucas & Kanade 1981; Lucas 1984; Baker & Matthews2004), Schuck (2006) developed the DAVE from a short time-expansion of the modified LCT method dis-cussed in Schuck (2005) for estimating velocities from line-of-sight magnetograms. The DAVE locally mini-mizes the square of the continuity equation (2b) subject to an affine velocity profile. Using “moving paint”experiments, Schuck (2006) demonstrated that this technique was faster and more accurate than existingLCT algorithms for data satisfying (2b). The DAVE method has been used to study the apparent motionof active regions (Schuck 2006), flux pile up in the photosphere (Litvinenko et al. 2007), and helicity flux inthe photosphere (Chae 2007). However, nagging questions remain about its performance.Welsch et al. (2007) set an important new standard for evaluating scientific optical flow methods usedfor studying the Sun. For the first time many existing methods for estimating photospheric velocities frommagnetograms were tested on a reasonable approximation to synthetic photospheric data from anelasticpseudospectral ANMHD simulations (Fan et al. 1999; Abbett et al. 2000, 2004). The methods tested wereLockheed Martin’s Solar and Astrophysical Laboratory’s (LMSAL) LCT code (DeRosa 2001), Fourier LCT(FLCT) (Welsch et al. 2004), the DAVE (Schuck 2006), the IM (Kusano et al. 2002, 2004), ILCT (Welsch et al.2004), the MEF (Longcope 2004), and MSR (Georgoulis & LaBonte 2006). Unfortunately the results werenot entirely encouraging. Welsch et al. (2007) treated the velocities estimated from line-of-sight methods asthe flux transport velocities consistent with the hypothesis of D´emoulin & Berger (2003) in (10). Evaluationof the DAVE’s performance on the ANMHD data under this assumption revealed that the DAVE methoddid not estimate the helicity flux or Poynting flux reliably. In fact none of the pure line-of-sight methods:LMSAL’s LCT, FLCT, or the DAVE—estimated these fluxes reliably, reproducing (at best) respectively 5 –11%, 9%, and 23% of the helicity rate, and reproducing respectively 6%, 11%, and 22% of the power injectedthrough the surface.Of course the ANMHD data have limitations. The simulation models the rise of a buoyant magneticflux rope in the convection zone and represents the magnetic structure of granulation or super-granulationrather than the dynamics of an active region (See § they fail to fully explain the poor performance of the tracking methodsto accurately reproduce the quantity they were designed to estimate, namely the helicity flux! This paper has two primary goals:1. Develop a modified DAVE (Schuck 2006) that incorporates horizontal magnetic fields, termed the“differential affine velocity estimator for vector magnetograms” (DAVE4VM), and demonstrate itsperformance on the ANMHD simulation data. DAVE4VM performs much better than the originalDAVE technique and roughly on par with the minimum energy fit (MEF) method developed by Long-cope (2004) which was deemed to have performed the best overall in Welsch et al.’s (2007) comparisonof velocity-inversion techniques.2. Identify the reasons for the poor performance of DAVE in Welsch et al. (2007).The paper attempts to follow, as closely as possible, the presentation of the DAVE in Schuck (2006) and theanalysis of velocity inversion techniques by Welsch et al. (2007). For the remainder of this paper, lower casevariables are used to represent the flux transport vector, flux transport velocity, plasma velocity, electricfield, Poynting flux, and helicity flux estimates from the DAVE4VM and DAVE: f , u , v ⊥ , e ⊥ , s z and h andthe corresponding uppercase variables are used to represent the “ground truth” from ANMHD: F , U , V ⊥ , E ⊥ , S z , and H . The one deviation from this notation involves ϑ which denotes an optical flow estimatebased on (6). Section (2) describes the DAVE4VM model and § ϑ = u . In § ϑ = u for the DAVE is relaxed and compared with an alternative hypotheses that ϑ = v h —that the DAVE produces a biased estimate of the total horizontal plasma velocity.
2. The DAVE4VM Model
The extension of the DAVE for horizontal magnetic fields is straight-forward. The plasma velocity ismodeled with a three-dimensional affine velocity profile: v ( P ; x ) = (cid:98) u (cid:98) v (cid:98) w + (cid:98) u x (cid:98) u y (cid:98) v x (cid:98) v y (cid:98) w x (cid:98) w y (cid:18) xy (cid:19) , (11) 6 –where the hatted variables model the local plasma velocity profile. The coordinate system for the affinevelocity profile is not aligned with the magnetic field. Therefore, the velocities are not guaranteed to beorthogonal to B . However, the parallel ( (cid:107) ) and perpendicular ( ⊥ ) components of the plasma velocity maybe determined from v (cid:107) = ( v · B ) B B , (12a) v ⊥ = v − ( v · B ) B B , (12b) v ⊥ h = v h − ( v · B ) B h B , (12c) v ⊥ z = v z − ( v · B ) B z B . (12d)The error metric C SSD = (cid:90) dtdx w ( x − χ , t − τ ) { ∂ t B z ( x , t ) + ∇ h · [ B z ( x , t ) v h ( P , x − χ ) , (13a) − v z ( P , x − χ ) B h ( x , t )] } , = η · (cid:104) S (cid:105) · η (13b)characterizes how well the local velocity profile satisfies the magnetic induction equation over a subregion ofthe magnetogram sequence defined by the window function w ( x − χ , t − τ ) where P = ( (cid:98) u , (cid:98) v , (cid:98) u x , (cid:98) v y , (cid:98) u y , (cid:98) v x , (cid:98) w , (cid:98) w x , (cid:98) w y )is a vector of parameters and η ≡ ( P , v ( P , x − χ ) in (13a) is referenced from thecenter of the window at x = χ so that (cid:98) u , (cid:98) v , and (cid:98) w represent the plasma velocities at the center ofthe window and the subscripted parameters represent the best fit local shears in the plasma flows, i.e. (cid:98) u x = ∂ x ( (cid:98) x · v ). The matrix elements of (cid:104) S (cid:105) are defined by (cid:104) S (cid:105) = (cid:90) dt dx w ( x − χ , t − τ ) S ( χ ; x, t ) , (13c)where S ( χ ; x, t ) ≡ (cid:20) A bb G (cid:21) = G · · · · · · · · ·G G · · · · · · · ·G G G · · · · · · ·G G G G · · · · · ·G G G G G · · · · ·G G G G G G · · · · s s s s s s s · · · s s s s s s s s · · s s s s s s s s s ·G G G G G G s s s G , (14)is a real symmetric S = S ∗ positive semidefinite structure tensor where a superscript “*” indicates thematrix transpose. The matrix elements of S are provided in Appendix A. The elements G ij correspond tothe original DAVE method (Schuck 2006) and the remainder s ij represent corrections due to the horizontalcomponents of the magnetic field and flows normal to the surface. The least-squares solution is P = − (cid:104) A (cid:105) − · (cid:104) b (cid:105) , (15) 7 –when the aperture problem is completely resolved det ( (cid:104) A (cid:105) ) (cid:54) = 0 and the velocity field is unambiguous. How-ever, there are important new terms in the structure tensor (cid:104) S (cid:105) involving B h . Situations where det ( (cid:104) A (cid:105) ) = 0because B h = 0 or B z = 0 over the region contained within the window must be considered. In general,the Moore-Penrose pseudo-inverse (cid:104) A (cid:105) † provides a numerically stable estimate of the optical flow parameterseven when det ( (cid:104) A (cid:105) ) = 0 P = − (cid:104) A (cid:105) † · (cid:104) b (cid:105) , (16a)where (cid:104) A (cid:105) † ≡ V Σ † U ∗ , (16b)is defined in terms of the singular value decomposition (Golub & Van Loan 1980) (cid:104) A (cid:105) = U Σ V ∗ . (16c)Here U and V are orthonormal matrices corresponding to the nine principle directions, Σ is a diagonalmatrix containing the nine singular values, and Σ † is computed by replacing every nonzero element of Σ byits reciprocal. If B h = 0, the singular values along the vertical direction are zero and (cid:104) A (cid:105) is rank deficient.In this case, the method implemented produces the minimum norm least squares solution resulting in novertical flows: w = w x = w y = 0.
3. Application to ANMHD Simulations
This paper considers the pair of vector magnetograms B separated by the shortest time interval ∆ t ≈
250 s between data dumps of the ANMHD simulation slice archived by Welsch et al. (2007). The “groundtruth” data are derived from the time-averaged velocity and magnetic fields from ANMHD over the shortest time interval. The region of interest in the ANMHD simulations corresponds to a 101 ×
101 pixel regioncentered on a convection cell. Welsch et al. (2007) thresholded on the vertical magnetic field and consideredonly pixels with | B z | >
370 G for all plots and quantities except for the total helicity where a differentmasking of results was used (Welsch et al. 2008). In a departure from the original presentation of Welschet al. (2007) this paper considers pixels with B = (cid:113) B x + B y + B z >
370 G. This corresponds to 7013 pixelsor roughly 70% of the region of interest. This difference in thresholding is important because most of the fluxemergence and helicity flux in the simulation occurs along the neutral line in regions of weak vertical field thatare missed with vertical field thresholding used in the original study. Note that the modified thresholdingmask contains weak vertical field regions and contains substantially more points than the roughly 3800 usedin Welsch et al. (2007). The difference between the helicity flux in this comparison region and the totalsimulation is 0.023%.Since the DAVE4VM is a local optical flow method that determines the plasma velocities within awindowed subregion by constraining the local velocity profile, the choice of window size is a crucial issue forestimating velocities accurately. The window must be large enough to contain enough structure to uniquelydetermine the coefficients of the flow profile and resolve the aperture problem, but not so large as to violatethe affine velocity profile (11) which is only valid locally (See Schuck 2006, for discussion of the apertureproblem in the context of the DAVE). In Welsch et al. (2007), the optimal window size was chosen for theDAVE by examining the Pearson correlation and slope between ∇ h · ( ϑ B z ) and ∆ B z / ∆ t from ANMHD.If the method satisfies the induction equation exactly everywhere, the Pearson correlation and slope wouldboth be equal to −
1. However both the DAVE and DAVE4VM were conceived with the recognition that realmagnetograms contain noise and should not satisfy the magnetic induction equation exactly; these methods 8 –satisfy the induction equation statistically within the window by minimizing the mean squared deviationsfrom the ideal induction equations. Consequently, how well these methods satisfy the magnetic inductionequation globally or over a subset of pixels can be used to assess overall performance. In Welsch et al. (2007),the DAVE was “optimized” over a subset of pixels with | B z | >
370 G that did not include the weak verticalfield regions discussed in this study. Using the Pearson correlation and slope, an asymmetric window of21 ×
39 performed the best on the ANMHD data. For the present study, the DAVE was “re-optimized”over the new criteria | B | >
370 G to provide an “apples-to-apples” comparison. However, I emphasizethat generally this optimization cannot be carried out for the DAVE because this method was proposed forderiving flux-transport velocities from line-of-sight magnetograms. In this situation, the threshold mask canonly be applied on the line-of-sight component as a proxy for the vertical magnetic field. Nonetheless, asa practical matter, understanding the accuracy limitation of the line-of-sight method in comparison to thevector method when synthetic vector magnetograms are available will reveal the relative reliability of helicityflux studies that used optical flow methods that rely exclusively on the line-of-sight magnetic field underthe “best case scenario” for the tracking methods: the true vertical magnetic field is tracked and regionsthat contain interesting physics are known a-priori . By using the “ground truth” vertical magnetic field,the evaluation is biased to favor the performance of the DAVE method over what would be possible underrealistic conditions where only the line-of-sight magnetic field is available.Five different criteria were used to optimize window selection. Only symmetric windows were consideredand some improvement in the results can be achieved by implementing asymmetric windows as in Welschet al. (2007). Figure 1 shows the optimization curves for the DAVE ( left ) and DAVE4VM ( right ). Thetop plots show the Spearman rank order ( ρ ; dashed black ) and Pearson ( C ; solid black curve ) correlationsand slope ( S ; red curve ) between ∇ h · ( ϑ B z ) and ∆ B z / ∆ t for DAVE and ∇ h · ( u B z ) and ∆ B z / ∆ t forDAVE4VM using pixels that satisfy | B | >
370 G. The gradients were computed with 5-point least-squaresoptimized derivatives (J¨ahne 2004). Window sizes between 15 and 30 pixels provide the best balance forachieving performance approaching C = ρ = S = −
1. Increasing the window size beyond 30 pixels continuesto improve the Spearman and Pearson correlations but with diminishing returns while the slopes degradesignificantly. The bottom plots show the power ( black curve ) and helicity rate ( red curve ) as a function ofwindow size for the DAVE and DAVE4VM. For DAVE these calculations require the assumption that ϑ = u .The horizontal black and red dashed lines correspond to the ground truth Poynting and helicity flux fromANMHD. The magnitude of these fluxes are maximum near 20 pixels with roughly uniform performancebetween 15 and 30 pixels. A window size of 23 pixels was chosen for both the DAVE and DAVE4VMas indicated by the vertical dot-dashed lines in Figure 1. These objective metrics for evaluation of globalperformance can be implemented without knowledge of ground truth. Only future tests with more realizationsof synthetic data can reveal whether they are robust metrics for optimizing window choice. Table 1 presentsa summary of the correlation coefficients on the mask | B | >
370 G characterizing the accuracy of the DAVEand DAVE4VM for the quantities discussed in this section.
Determining the flux transport vectors f = u B z is equivalent to determining the perpendicular plasmavelocities v ⊥ . The accuracy of the flux transport vectors is critical for estimating other MHD quantities: Only windows with odd numbers of pixels on each side of the window are possible in these implementations of the DAVEand DAVE4VM. left ) and DAVE4VM ( right ). (top) The Spearman rank order( dashed black curve ) and Pearson ( solid black curve ) correlations and slope ( red curve ) between ∇ h · ( ϑ B z )and ∆ B z / ∆ t for DAVE and ∇ h · ( u B z ) and ∆ B z / ∆ t for DAVE4VM as a function of window size. (bottom)The power (black) and helicity (red) as a function of window size; this assumes ϑ = u for DAVE. Thehorizontal black and red dashed lines correspond to the ground truth power and helicity rate through thesimulation slice from ANMHD. The vertical dot-dashed lines indicate the “optimal” window size of 23 pixelschosen for both DAVE and DAVE4VM. 10 –Table 1. Comparison between the DAVE and DAVE4VM over the 7013 pixels that satisfy | B | >
370 G inFigure 2. This mask contain regions of weak vertical field not considered in Welsch et al. (2007).DAVE (Assuming ϑ = u ) DAVE4VMQuantities Spearman Pearson Slope Spearman Pearson Slope u x B z U x B z u y B z U y B z v ⊥ x V ⊥ x v ⊥ y V ⊥ y v ⊥ z V ⊥ z ∇ h · ( u B z ) ∆ B z / ∆ t -0.85 -0.95 -0.97 -0.79 -0.94 -0.97 e ⊥ x E ⊥ x e ⊥ y E ⊥ y e ⊥ z E ⊥ z s z S z | B | >
370 G in Figure 2. This mask contains regions of weak vertical field notconsidered in Welsch et al. (2007). Here (cid:104)| δ (cid:101) f |(cid:105) is the average fractional error, (cid:104) δ | (cid:101) f |(cid:105) is the average error inmagnitude, C vec is the vector correlation, C CS is the direction correlation, and (cid:104) cos θ (cid:105) W is weighteddirection cosine.DAVE (Assuming ϑ = u ) DAVE4VM f (cid:104)| δ (cid:101) f |(cid:105) (cid:104) δ | (cid:101) f |(cid:105) C vec C CS (cid:104) cos θ (cid:105) W (cid:104)| δ (cid:101) f |(cid:105) (cid:104) δ | (cid:101) f |(cid:105) C vec C CS (cid:104) cos θ (cid:105) W u B z ± ± ± ± v ⊥ ± ± ± ± B z overlaid with the horizontal magneticfield vectors B h in aqua. Blue contours indicate smoothed neutral lines. Top right: Distribution of anglesfor B h in the horizontal plane. Bottom: (red arrows) ϑ B z for the DAVE ( left ) and u B z for the DAVE4VM( right ) plotted over ANMHD’s flux transport vectors U B z (green arrows). Vectors are shown only for pixelsin which | B | >
370 G, and for clarity, only every third vector is displayed. 12 –perpendicular plasma velocities, electric field, helicity flux, Poynting flux, etc, since all of these quantitiesmay be derived directly from flux transport vectors.
The top left of Figure 2 shows the region of interest from the ANMHD simulations with grey scaleimage of vertical magnetic field overlaid with the horizontal magnetic field vectors B h in aqua. The bluecontours indicate smoothed neutral lines. The top right shows the distribution of angles for B h in thehorizontal plane. The horizontal magnetic field is largely aligned with the (cid:98) x -axis as indicated by the aquavectors in the left panel and the strong peak in the histogram near arctan ( B y , B x ) ≈ ◦ in the right panel. There is also significant alignment of the magnetic field with ± ◦ and alignment of weak fields with − ◦ .The bottom panels show ϑ B z from DAVE ( left ) and u B z from DAVE4VM in red ( right ) and the fluxtransport vectors F = U B z from ANMHD in green. The improvement between the DAVE and DAVE4VMis manifest — finding a region where the DAVE4VM performs qualitatively worse than the DAVE is difficult.The DAVE4VM performs the worst in the region 140 − × −
170 where the flux transport vectors runroughly anti-parallel to the horizontal magnetic field and there is little structure in the vertical component.Figure 3 shows scatter plots of the ϑ B z from DAVE ( left ) and u B z from DAVE4VM ( right ) versusthe flux transport vectors U B z from ANMHD. Red and blue are used to distinguish x - and y -components,respectively. The nonparametric Spearman rank-order correlation coefficients ( ρ ), Pearson correlation coeffi-cients ( C ), and slopes ( S ) estimated by the least absolute deviation method are shown for both componentsof the flux transport vectors. Both visually and quantitatively the DAVE4VM’s correlation with ANMHD ismuch higher than the DAVE’s. The correlation coefficients even match or exceed the correlation coefficientsfor the flux transport vectors from the DAVE and MEF reported for the restricted mask | B z | >
370 G inWelsch et al. (2007). In particular, DAVE does not accurately estimate the flux transport vectors in the (cid:98) x -direction. The correlation coefficients for this (cid:98) x -component of the flux transport vectors are ρ x = 0 . C x = 0 .
57 with a slope of S x = 0 .
15. Since that the (cid:98) x -direction is the predominant direction of thehorizontal magnetic for the ANMHD data, the low correlation coefficients suggest that DAVE is insensitiveto flux emergence which is proportional to v z B h . This will be discussed further in § v ⊥ from the DAVE as-suming ϑ = u ( left ) and DAVE4VM ( right ) versus ANMHD’s perpendicular plasma velocities V ⊥ . Red,blue, and black correspond to the x -, y -, and z -components respectively. The nonparametric Spearmanrank-order correlation coefficients ( ρ ), Pearson correlation coefficients ( C ), and slopes ( S ) estimated by theleast absolute deviation method are shown for each component of the perpendicular plasma velocities. TheDAVE4VM’s correlation coefficients match or exceed the correlation coefficients for the perpendicular plasmavelocities from the DAVE. Particularly striking is the DAVE4VM’s relatively higher correlation for the per-pendicular vertical plasma velocity v ⊥ z which exceeds the correlation for the DAVE by roughly 0 . − . explicit inclusion of horizontal magnetic fieldsand vertical flows. The flux transport and perpendicular plasma velocity estimates are further quantified byconsidering the metrics used by Schrijver et al. (2006), Welsch et al. (2007), and Metcalf et al. (2008). The arctan ( y, x ) ≡ arctan ( y/x ).
13 –Fig. 3.— Scatter plots of ϑ B z for the DAVE ( left ) and u B z for the DAVE4VM ( right ) versus ANMHD’sflux transport vectors, F = U B z . Red and blue are used to distinguish x - and y -components, respectively.The nonparametric Spearman rank-order correlation coefficients ( ρ ), Pearson correlation coefficients ( C ),and slopes ( S ) estimated by the least absolute deviation method are shown for both components of the fluxtransport vectors.Fig. 4.— Scatter plots of the estimated perpendicular plasma velocities v ⊥ from DAVE assuming ϑ = u ( left ) and DAVE4VM ( right ) versus the perpendicular plasma velocities V ⊥ from ANMHD. Red, blue,and black correspond to the x -, y - and z -components respectively. The nonparametric Spearman rank-ordercorrelation coefficients ( ρ ), Pearson correlation coefficients ( C ), and slopes ( S ) are shown for each componentof the perpendicular plasma velocities. 14 –fractional error between the estimated vector f and the true vector F at the i th pixel is | δ (cid:101) f i | ≡ | f i − F i || F i | , (17a)whereas the fractional error in magnitude is δ | (cid:101) f i | ≡ | f i | − | F i || F i | . (17b)The moments of these error metrics or any quantity q may be accumulated over the N pixels within themasks (either | B | >
370 G or | B z | >
370 G) producing the average (cid:104) q (cid:105) ≡ N N (cid:88) i =1 q i , (18a)and the variance σ q ≡ N − N (cid:88) i =1 ( q i − (cid:104) q (cid:105) ) . (18b)For perfect agreement between the estimates and the “ground truth” from ANMHD, (cid:104)| δ (cid:101) f |(cid:105) , (cid:104) δ | (cid:101) f |(cid:105) , and theirassociated variances would be zero. Two measures of directional error are considered, the vector correlation C vec = (cid:104) f · F (cid:105) (cid:113)(cid:10) f (cid:11) (cid:10) F (cid:11) , (19a)and the direction correlation C CS = (cid:42) f · F (cid:112) f F (cid:43) ≡ (cid:104) cos θ (cid:105) . (19b)Both metrics range from − | B | >
370 G. The DAVE4VM has fractional errors less than or equal to 0 . . − .
09 and − .
47 for DAVE4VM and DAVE respectivelywhich corresponds to a factor of 5 improvement. For the plasma velocity, the bias error in magnitude is 0.01and 0.09 for DAVE4VM and DAVE respectively which corresponds to a factor of 9 improvement. The vectorcorrelation is larger for DAVE4VM than for DAVE. For DAVE4VM C vec (cid:38) . C vec = 0 .
61 and the perpendicular plasma velocities with C vec = 0 .
81. Finally the directional errors are smaller for the DAVE4VM than for the DAVE. For theDAVE4VM C CS (cid:38) . C CS = 0 . C CS = 0 . C CS is difficult to translate into average angular error because it is a nonlinearfunction of θ and does not indicate whether the estimated vectors “lead” or “lag” the “ground truth” onaverage. For the 2D flux transport vectors the moments of the distribution of angular errors θ = arctan [( u × U ) z , u · U ] , (20) 15 –can be more informative. Figure 5 shows histograms of the angles between u B z and U B z for the DAVE( left ) and DAVE4VM ( right ). This is a quantitative estimate of the errors in directions of the flux transportvectors. The DAVE4VM represents a dramatic improvement over the DAVE. The DAVE4VM produces anearly unimodal distribution peaked near 0 ◦ , whereas the DAVE produces a multi-peaked distribution withthe largest peak at 80 ◦ and a variance that is more than twice as large as DAVE4VM.Metrics such as (17) and (19) weight all estimates equally. To address this Metcalf et al. (2008) suggestedweighting the errors. For example, the weighted direction cosine between an inferred vector and the groundtruth vector may be defined as (cid:104) cos θ (cid:105) W ≡ (cid:104) W cos θ (cid:105)(cid:104) W (cid:105) (21)where W represents weights. For the flux transport velocity, the errors in the orientation of u are moreimportant where | U B z | is large and less important where | U B z | is small which suggests a weighting factor W i = | ( U B z ) | . For perfect agreement (cid:104) cos θ (cid:105) W = 1. The weighted direction cosines for the flux transportvectors and the plasma velocities are reported in Table 2. Comparing the values of C CS and (cid:104) cos θ (cid:105) W demonstrates that weighting the direction cosine improves the apparent performance of DAVE but theresults for DAVE4VM are essentially unchanged. This suggests that DAVE4VM estimates velocities betterthan DAVE in regions of weak flux transport. Under ideal conditions, the magnetic field is only affected by ∇× ( v × B ). Consequently, only theinductive potential φ in (7) may be uniquely determined from the evolution of B alone. The electrostaticpotential ψ must and can be estimated with additional judicious assumptions. These additional assumptionscorrespond to the minimum photospheric velocity consistent with (7) for the global method MEF and tothe prescribed affine form of the local plasma velocity for the local method DAVE4VM. The constraint ofthe affine velocity profile permits DAVE4VM to determine the electrostatic potential ψ from the nonlocal structure of the inductive potential φ . DAVE4VM uses unambiguous “pieces” of the plasma velocity withinthe window aperture to reconstruct the total plasma velocity at the center of the aperture. Within thenotation of the Helmholtz decomposition, DAVE4VM estimates the local electrostatic field from the structureof the nonlocal inductive field within the aperture window by imposing a smoothness constraint on thevelocity (the affine velocity profile). The accuracy of this estimate for the electrostatic field depends onthe validity of the local affine velocity profile and the amount of structure in the local magnetic field; localmethods cannot detect motion in regions of uniform magnetic field.While researchers have widely recognized that estimating the electrostatic potential ψ from (5a) requiresadditional assumptions, they have not generally recognized that the parallel velocity v (cid:107) may be estimatedby the analogous arguments. In the absences of a reference flow, the MEF constrains the velocity to beperpendicular to the magnetic field: v (cid:107) = 0. In contrast to the velocity estimated by DAVE4VM is notconstrained to be perpendicular to the local magnetic field. Instead, DAVE4VM fits an affine velocity modelto the magnetic induction equation in an aperture window. This affine velocity model couples the dynamicsacross pixels within the window aperture. If there is sufficient structuring in the direction of the magneticfield within the aperture, i.e., the perpendicular plasma velocity points in different directions at differentpixels within the aperture, then DAVE4VM can resolve the ambiguity in the field-aligned component of theplasma velocity at the center of the aperture. 16 –Fig. 5.— Histograms of the angles between ϑ B z and U B z for the DAVE ( left ) and between u B z and U B z for the DAVE4VM ( right ). The mean (bias) and standard deviation are reported in the upper left-handcorners.Fig. 6.— Schematic diagram of a uniform plasma flow across a diverging magnetic field above the photosphereat h = 0. The black arrows indicate the strength and direction of the magnetic field, the red arrows indicatethe direction of the spatially uniform total plasma velocity v , and the blue arrows indicate the magnitudeand direction of the perpendicular plasma velocity v ⊥ . The aperture in the photosphere is indicated by thegray box. 17 –Consider the simplified two-dimensional situation illustrated by the schematic diagram in Figure 6 ofspatially uniform plasma flow across a diverging magnetic field above the photosphere at h = 0. The blackarrows indicate the magnitude and direction of the magnetic field, the red arrows indicate the direction ofthe spatially uniform total plasma velocity v , and the blue arrows indicate the magnitude and direction ofthe perpendicular plasma velocity v ⊥ . The aperture in the photosphere is indicated by the gray box. Withinthe aperture, the perpendicular plasma velocity captures a different component of the total plasma velocityat different locations; this is a consequence of the structuring of the magnetic field. Under the smoothnessassumption of a uniform velocity profile, the velocity along the magnetic field in the (cid:98) z -direction at x = 0may be determined from the components of the perpendicular plasma velocity in the (cid:98) z direction at otherlocations within the aperture . Using the N pixels in the window aperture results in an overdetermined systemfor the total plasma velocity : B z ( x ) B ( x ) − B x ( x ) B z ( x ) B ( x ) − B x ( x ) B z ( x ) B ( x ) B x ( x ) B ( x ) ... ... B z ( x N ) B ( x N ) − B x ( x N ) B z ( x N ) B ( x N ) − B x ( x N ) B z ( x N ) B ( x N ) B x ( x N ) B ( x N ) (cid:124) (cid:123)(cid:122) (cid:125) D (cid:18) v x v z (cid:19) = v ⊥ x ( x ) v ⊥ z ( x )... v ⊥ x ( x N ) v ⊥ z ( x N ) (cid:124) (cid:123)(cid:122) (cid:125) d , (22a)which has the solution (Golub & Van Loan 1980) (cid:18) (cid:98) v x (cid:98) v z (cid:19) = ( D ∗ D ) − D ∗ d . (22b)Note that D ∗ D is analogous to (cid:104) A (cid:105) and D ∗ d is analogous to (cid:104) b (cid:105) in (15).This pedagogical example illustrates how DAVE4VM may analogously estimate the field-aligned plasmavelocity for the more general case of a spatially variable plasma flow in an inhomogeneous magnetic fieldfor (15). The accuracy of the estimate of the parallel velocity will be limited by the structuring in directionof the magnetic field within the aperture; if the magnetic field has a uniform orientation in the aperturewindow, no useful estimate of the field-aligned plasma velocity can be made from the magnetic measurementsalone. The quality of the estimate may be assessed with the conditioning of D ∗ D in (22b) for the pedagogicalexample or (cid:104) A (cid:105) in (15) for the full system.Figure (7) shows scatter plots of ( left ) the estimated parallel plasma velocities v (cid:107) from DAVE4VMversus the parallel plasma velocities V (cid:107) from ANMHD and ( right ) the estimated total plasma velocity v fromDAVE4VM versus the total plasma velocities V from ANMHD. The nonparametric Spearman rank-ordercorrelation coefficients ( ρ ), Pearson correlation coefficients ( C ), and slopes ( S ) estimated by the least absolutedeviation method are shown. The comma-separated pairs of numbers in the right plot, corresponding tocorrelations between v and V and v ⊥ and V respectively, represent the relative improvement in total velocityestimate over the simple null hypothesis H : v = v ⊥ that the total plasma velocity is the perpendicularplasma velocity. The correlation of the (cid:98) x -component of total velocity is significantly improved over the nullhypothesis H . This improvement id interesting since the horizontal magnetic field is predominantly alignedwith the x -axis (See Figure 2). The correlation of the (cid:98) y -component of total velocity is slightly worse thanthe null hypothesis. Finally, the correlation of the (cid:98) z -component of total velocity is mixed with the Spearmancorrelation ρ slightly worse than the null hypothesis and the Pearson correlation C slightly better than thenull hypothesis. However, the slopes of all three components are improved over the null hypothesis. 18 –The significance of the correlations in the left plot may be tested against the null hypothesis H : ρ = 0by the Fisher permutation test. Fieller et al. (1957) have demonstrated with analysis backed Monte-Carlosimulation that Fisher’s z -transform of the correlation coefficient z S ( ρ ) = 12 log (cid:12)(cid:12)(cid:12)(cid:12) ρ − ρ (cid:12)(cid:12)(cid:12)(cid:12) , (23)produces approximately normally distributed values. For example, permuting the values of v (cid:107) and V (cid:107) (cid:104) z S (cid:105) = 0 . ± .
01. The Spearman correlation coefficient ρ = 0 .
53 has a z -transform of z S (0 .
53) = 0 .
60 which is roughly 50 standard deviations from the mean ofthe null distribution indicating that the parallel velocity correlation is statistically significant and not dueto sampling error. However, the correlation ρ = 0 .
53 is small and the parallel velocity estimates may not bescientifically significant for accurately predicting the parallel velocity. The plasma velocities may be furtherconstrained by introducing Doppler velocities, but this is beyond the scope of the present discussion. v h and v z Redundant?
DAVE4VM has incorporated an additional component of the velocity over DAVE by introducing threeadditional variables (cid:98) w , (cid:98) w x and (cid:98) w y . Consequently, one may reasonably wonder “are the terms v h B z and v z B h redundant for DAVE4VM?” The answer is a clear “No” for the ANMHD data. Equation (2a) iscomposed of two terms B z V h and V z B h describing shearing motion and emergence respectively. Figure 8shows scatterplots of the two terms for the (cid:98) x ( left ) and (cid:98) y ( right ) components of (2a). The scatterplots indicatethe lack of correlation between the terms describing shearing motions B z V h and emergence V z B h . Redpoints indicate the results for DAVE4VM and blue points indicate the results for ANMHD. The Spearmanrank order ( ρ ) and Pearson ( C ) correlations between the two terms, summarized in Table 3, are very lowfor both components of the flux transport velocities from DAVE4VM or ANMHD. These terms describedifferent physics, that are uncorrelated, and which require independent variables to describe. Figure 9 shows ∇ h · ( ϑ B z ) from the DAVE ( left ) and ∇ h · ( u B z ) from DAVE4VM ( right ) versus ∆ B z / ∆ t from ANMHD. The derivatives for these plots were estimated from 5-point optimized least squares. Theseplots indicate how well the two methods satisfy the MHD induction equation globally. The DAVE has higherTable 3. The Spearman rank order ( ρ ) and Pearson ( C ) correlations between the terms in the fluxtransport velocity describing shearing motion and flux emergence.Correlates DAVE4VM ANMHDSpearman Pearson Spearman Pearson v ⊥ x B z v ⊥ z B x -0.23 -0.08 -0.12 -0.01 v ⊥ y B z v ⊥ z B y -0.07 -0.11 -0.09 -0.02 v ⊥ x v ⊥ z -0.09 -0.10 -0.15 -0.15 v ⊥ y v ⊥ z -0.34 -0.33 -0.30 -0.20 19 –Fig. 7.— Scatter plots of ( left ) the estimated parallel plasma velocities v (cid:107) from DAVE4VM versus the parallelplasma velocities V (cid:107) from ANMHD and ( right ) the estimated total plasma velocity v from DAVE4VM versusthe total plasma velocities V from ANMHD. The nonparametric Spearman rank-order correlation coefficients( ρ ), Pearson correlation coefficients ( C ), and slopes ( S ) estimated by the least absolute deviation methodare shown. The pairs of numbers represent correlations between v and V and v ⊥ and V .Fig. 8.— Scatterplots for the (cid:98) x ( left ) and (cid:98) y ( right ) components of (2a). The scatterplots indicate the lack ofcorrelation between the terms describing shearing motion B z V h and emergence V z B h . Red points indicatethe results for DAVE4VM and blue points indicate the results for ANMHD. The Spearman rank order ( ρ )and Pearson ( C ) correlations between the two terms are very low for both components of the flux transportvelocities. 20 –correlations than the DAVE4VM but the slopes are equivalent. For the DAVE4VM, the most significantdeviations from the MHD induction equation occur near ∆ B z / ∆ t ≈
0. Neither the DAVE nor DAVE4VMsatisfy the induction equation exactly. This is by design, because real magnetogram data are likely to containsignificant noise which will contaminate velocity estimates if the induction equation is satisfied exactly.Furthermore, how well a method satisfies the induction equation will generally depend on the differencingtemplate. Consequently, if the velocity estimates are to be used as boundary values for ideal MHD coronalfield models, then the velocities of any method will have to be adjusted to satisfy the induction equationon the differencing template implemented by the simulation. Using the Helmholtz decomposition (7), theinductive potential may be computed for the simulation directly from the magnetogram sequence (on thesimulation differencing template) ∂ t B z = ∇ h φ, (24a)and the electrostatic potential may be derived from the flux transport vectors determined by the optical flowmethod (Welsch et al. 2004) ∇ h ψ = (cid:98) z · [ ∇× ( u B z )] . (24b)Incorporating photospheric velocity estimate into boundary conditions for a coronal MHD simulation, in aminimally consistent way with the normal component of the magnetic induction equation, requires solvingtwo Poisson equations on the photospheric boundary using the differencing template of the MHD code.Figure 10 shows scatter plots of the estimated perpendicular electric fields e ⊥ from DAVE assuming ϑ = u ( left ) and from DAVE4VM ( right ) versus the electric fields E ⊥ from ANMHD. Red, blue, andblack correspond to the x -, y -, and z -components, respectively. The nonparametric Spearman rank-ordercorrelation coefficients ( ρ ) and Pearson correlation coefficients ( C ) are shown for each component of theelectric field. On the present mask the DAVE4VM estimates improve or essentially match the correlationand slopes of the DAVE’s estimates for all three components of the electric field. Particularly dramatic is theimprovement in the (cid:98) y component of the electric field which the DAVE does not estimate accurately eitheron the present mask | B | >
370 G or the restricted mask | B z | >
370 G (Welsch et al. 2007).
D´emoulin & Berger (2003) show that the Poynting flux can be expressed concisely in terms of the fluxtransport vectors u B z s z ( x ) = − π B h · ( B z v h − v z B h ) = − B h · ( u B z )4 π . (25)Figure 11 shows scatterplots of the estimated Poynting flux s z from the DAVE assuming ϑ = u ( left ) andDAVE4VM ( right ) versus ANMHD’s Poynting flux S z . The correspondence for DAVE4VM, or lack there offor DAVE, indicates the accuracy of the velocity estimates in the direction of the horizontal magnetic field B h . The nonparametric Spearman rank-order correlation coefficients ( ρ ), Pearson correlation coefficients( C ), and slopes ( S ) estimated by the least absolute deviation method are shown, as is the ratio of theintegrated estimated Poynting flux to the integrated ANMHD Poynting flux R s z = (cid:80) s z / (cid:80) S z . TheDAVE4VM’s estimate of Poynting flux is a significant improvement over the DAVE’s. The correlationshave improved by roughly a factor of 4 −
6, the slope has improved by nearly a factor of 18, and theratio of the totals has improved by nearly a factor of 5. Again, DAVE does not reliably estimate the fluxtransport velocity in the direction of the horizontal magnetic field B h suggesting that DAVE is insensitive 21 –Fig. 9.— Scatter plots of ∇ h · ( ϑ B z ) from the DAVE ( left ) and ∇ h · ( u B z ) from the DAVE4VM ( right )versus ∆ B z / ∆ t from ANMHD. The nonparametric Spearman rank-order correlation coefficients ( ρ ), Pearsoncorrelation coefficients ( C ), and slopes ( S ) estimated by the least absolute deviation method are shown.Fig. 10.— Scatter plots of the estimated perpendicular electric field e ⊥ from DAVE assuming ϑ = u ( left )and from DAVE4VM ( right ) versus the electric field E ⊥ from ANMHD. Red, blue, and black correspond tothe x -, y -, and z -components respectively. The nonparametric Spearman rank-order correlation coefficients( ρ ), Pearson correlation coefficients ( C ), and slopes ( S ) estimated by the least absolute deviation methodare shown for each component of the electric field. 22 –Fig. 11.— Scatter plots of the estimated Poynting flux from DAVE assuming ϑ = u ( left ) and DAVE4VM( right ) versus the Poynting flux from ANMHD. The nonparametric Spearman rank-order correlation coef-ficients ( ρ ), Pearson correlation coefficients ( C ), and slopes ( S ) estimated by the least absolute deviationmethod are shown, as is the ratio of the integrated estimated Poynting flux to the integrated ANMHDPoynting flux.Fig. 12.— Scatter plots of the estimated helicity flux DAVE assuming ϑ = u ( left ) and DAVE4VM ( right )versus ANMHD’s helicity flux. The nonparametric Spearman rank-order correlation coefficients ( ρ ), Pearsoncorrelation coefficients ( C ), and slopes ( S ) estimated by the least absolute deviation method are shown, asis the ratio of the integrated estimated helicity flux to the integrated ANMHD helicity flux. 23 –to flux emergence which is proportional to v z B h . The “ground truth” total power through the mask is dP/dt = (cid:80) S z = 7 . × ergs / s.D´emoulin & Berger (2003) show that the gauge-invariant helicity flux (Berger & Field 1984) can beexpressed concisely in terms of the flux transport vectors u B z g A ( x ) = − A p · ( B z v h − v z B h ) = − A p · ( u B z ) (26)where A p = (cid:98) z ×∇ Φ p is the potential reference field (with zero helicity) which satisfies (cid:98) z · ( ∇× A p ) = ∇ h Φ p = B z , (27)and ∇ · A p = (cid:98) z · A p = 0. To estimate the helicity flux density, Φ p was computed on a 257 ×
257 squarecentered on the region of interest with Dirichlet boundary conditions using MUDPACK (Adams 1993).While interpretation of maps of helicity flux g A ( x ) through the photosphere is problematic (Pariat et al.2005, 2007), a comparison of g A ( x ) estimated from the DAVE or DAVE4VM verses G A ( x ) calculatedfrom ANMHD indicates the accuracy of the estimated flux transport vectors in the direction of the vectorpotential. Figure 12 shows scatter plots of the estimated helicity flux from DAVE assuming ϑ = u ( left )and DAVE4VM ( right ) versus ANMHD’s helicity flux. The nonparametric Spearman rank-order correlationcoefficients ( ρ ) and Pearson correlation coefficients ( C ) are shown, as is the ratio of the integrated estimatedhelicity flux to the integrated ANMHD helicity flux R g A = (cid:80) g A / (cid:80) G A . The DAVE4VM’s estimatesrepresent a significant improvement over the DAVE’s, improving the correlation coefficients by roughly 0 . .
2. Furthermore, the ratio of totals has improved by roughly a factor of 4 from 0 . .
94 for the DAVE4VM. The “ground truth” helicity injected through the surface is dH A /dt = (cid:80) G A = − . × Mx / s.
4. Discussion and Conclusions
For completeness, Tables 4 and 5 provide a summary of metrics and correlation coefficients for the DAVEand DAVE4VM on the original mask | B z | >
370 G used by Welsch et al. (2007) in the same format as inTables 1 and 2. MEF performed the best overall in the original study by Welsch et al. (2007) although therewere some metrics where the DAVE outperformed MEF such as in the accuracy of the plasma velocities listedin Table 5 (compare with Figure 8 in Welsch et al. (2007)). The DAVE4VM’s estimates are a substantialimprovement over the results of the DAVE assuming ϑ = u on this mask. Comparing the rank-orderSpearman correlation coefficients for the flux transport vectors, perpendicular plasma velocity, and electricfield in Table 4, the DAVE4VM equals or out-performs MEF. Particularly, the DAVE4VM’s estimate of thevertical perpendicular plasma velocity is substantially better than MEF with a rank-order of 0.76 in theformer and 0.61 in the latter case. Accurate vertical flows are necessary to diagnose flux emergence andaccurately estimate the helicity flux. The one area where MEF exhibits superiority is in the estimate ofthe Poynting flux where the DAVE4VM captures 76% and MEF captures 100% (Welsch et al. 2007). Thefractional errors (cid:104)| δ (cid:101) f |(cid:105) and (cid:104) δ | (cid:101) f |(cid:105) are substantially lower than the DAVE for both the flux transport vectorsand plasma velocities. The DAVE had the largest vector correlation C vec and the direction correlation C CS in the original study and the DAVE4VM improves over this performance exhibiting correlation coefficients of This estimate differs by about 10% from the helicity estimate in Welsch et al. (2007). The discrepancy is caused bythe different methodologies and boundaries used for computing the vector potential A p . MUDPACK was used in this studywith (27) whereas Welsch et al. (2007) used a Green’s function scheme to compute A p .
24 –Table 4. Comparison between the DAVE and DAVE4VM over the 3815 pixels that satisfy | B z | >
370 G inFigure 2. This corresponds roughly to the mask used in Welsch et al. (2007).DAVE (Assuming ϑ = u ) DAVE4VMQuantities Spearman Pearson Slope Spearman Pearson Slope u x B z U x B z u y B z U y B z v ⊥ x V ⊥ x v ⊥ y V ⊥ y v ⊥ z V ⊥ z ∇ h · ( u B z ) ∆ B z / ∆ t -0.92 -0.95 -0.99 -0.84 -0.95 -0.97 e ⊥ x E ⊥ x e ⊥ y E ⊥ y e ⊥ z E ⊥ z s z S z | B z | >
370 G in Figure 2. This corresponds roughly to the mask used in Welschet al. (2007).DAVE (Assuming ϑ = u ) DAVE4VM f (cid:104)| δ (cid:101) f |(cid:105) (cid:104) δ | (cid:101) f |(cid:105) C vec C CS (cid:104) cos θ (cid:105) W (cid:104)| δ (cid:101) f |(cid:105) (cid:104) δ | (cid:101) f |(cid:105) C vec C CS (cid:104) cos θ (cid:105) W u B z ± ± ± ± v ⊥ ± ± ± ± .
9. The improvement for the flux transport vectors is particularly dramatic. The plasma velocitiesare more accurate and exhibit considerably less bias than those reported for MEF in Welsch et al. (2007).The DAVE4VM offers some minor advantages over MEF. The DAVE4VM is somewhat faster than MEF;the DAVE4VM(DAVE) requires 30(10) seconds to process the full 288 ×
288 pixel frame from ANMHDwhereas MEF requires roughly 10 minutes to converge on a reduced mask of the ANMHD data (privatecommunication with Belur Ravindra). The DAVEVM is local and directly estimates velocities across neutrallines and across broader weak field regions whereas MEF is a iterative global method that requires judiciouschoice of boundaries to ensure convergence. In concert, the DAVE4VM’s velocity estimate might be used withMEF either as an initial guess for the electrostatic potential ψ via (24b) or as ancillary inaccurate velocitymeasurements in the MEF variational term that constrains the photospheric plasma velocities (Longcope2004). Since the DAVE4VM is fast and does not require supervision beyond choosing a window size (andeven this could be automated according to the criteria discussed in § SolarDynamics Observatory .What is responsible for the DAVE4VM’s improved performance? The only differences between theDAVE and DAVE4VM are the terms s ij in the structure tensor (14) that describe the local structure of thehorizontal magnetic fields necessary for the description of vertical flows. There are two circumstances whenline-of-sight methods such as the DAVE, LMSAL’s LCT, and FLCT will produce accurate estimates of theflux transport velocity:1. B h = 0: The magnetic field is purely vertical and ϑ = v h = v ⊥ h . If the horizontal magnetic fieldsand their associated derivatives are zeroed, the DAVE and DAVE4VM produce identical flux transportand perpendicular plasma velocity estimates. The DAVE is consistent with the assumption that themagnetic field is purely vertical:lim B h → ∂ t B z + ∇ h · ( B z v h − v z B h ) = ∂ t B z + ∇ h · ( B z v h ) . (28a)2. v z = 0: There are no net upflow/downflows and ϑ = v h (cid:54) = v ⊥ h . In this situation, there must beprojected vertical flows along the magnetic field to cancel any projected vertical flow perpendicularto the magnetic field with v ⊥ z = − v (cid:107) z = − B z ( v h · B h ) /B . Consequently, v (cid:107) (cid:54) = 0. The DAVE isconsistent with the assumption that there are no vertical flows v z = 0:lim v z → ∂ t B z + ∇ h · ( B z v h − v z B h ) = ∂ t B z + ∇ h · ( B z v h ) . (28b)Both limits (28a) and (28b) are isomorphic with (6). By induction, (6) is consistent with the assumptionsleading to (28a) and (28b). Since DAVE does not consider corrections s ij due to the horizontal magneticfield, ϑ should generally be considered a biased estimate of the horizontal plasma velocity ϑ = v h and notthe flux transport velocity! Formally the alternative hypothesis H : ϑ = v h may be tested against the nullhypothesis H : ϑ = u of D´emoulin & Berger (2003). The null hypothesis H : ϑ = u is represented by theleft panel in Figure 3. The alternative hypothesis, represented by H : ϑ = v h , is characterized by the scatter The routines were all coded in Interactive Data Language (IDL 2002) and the computations were performed on a dualprocessor AMD Opteron 240 running at 1.4 GHz with a one megabyte memory cache and ten gigabytes of Random AccessMemory.
26 –plot of the estimated velocities ϑ B z for DAVE versus the horizontal plasma velocities V h B z from ANMHDin Figure 13. The nonparametric Spearman rank-order correlation coefficients ( ρ ), Pearson correlationcoefficients ( C ), and slopes ( S ) estimated by the least absolute deviation method are all significantly betterfor the alternative hypothesis than for the null hypothesis. The null hypothesis that the velocities inferred byDAVE represent the flux transport velocities may be rejected in favor of the alternative hypothesis ϑ = v h .These results explain why V ⊥ z and v ⊥ z are poorly correlated for the DAVE in Figure 4 and the slopebetween them is nearly zero — the DAVE is consistent with the assumption v z = 0 when B h (cid:54) = 0. Generally,in regions of flux emergence, the accuracy v ⊥ z is critical for estimating the flux transport vectors which inturn is critical for estimating the helicity and Poynting fluxes. When horizontal magnetic fields and verticalflows are present, the flux transport vectors estimated from methods that rely exclusively on the line-of-sightor vertical component (DAVE, LMSAL’s LCT, FLCT) cannot be trusted to provide the total fluxes. Thisis particularly true along neutral lines where flux is emerging or submerging! Under the best case scenarios ,only the shearing or “horizontal fluxes” across the photosphere. dpdt (cid:12)(cid:12)(cid:12)(cid:12) h = − π (cid:90) S dx B h · ( v h B z ) , (29a)and dh A dt (cid:12)(cid:12)(cid:12)(cid:12) h = − (cid:90) S dx A p · ( v h B z ) , (29b)may be estimated from the line-of-sight tracking methods. In the best case scenarios only the shearingfluxes are captured by line-of-sight tracking methods in partial agreement with the Ansatz of Chae (2001)and in disagreement with the geometrical arguments of D´emoulin & Berger (2003) who argue that line-of-sight tracking methods capture both the shearing and emergence. Again, for the energy and helicity,the alternative hypothesis H : ϑ = v h (shearing) may be tested against the null hypothesis H : ϑ = u (shearing and emergence). The left panel of Figure 14, representing the null hypothesis is a combination ofthe left-hand panels of Figure 11 and Figure 12. The right panel of Figure 14, representing the alternativehypothesis, combines the shearing term estimated from DAVE with the emergence term from ANMHD.Generally, in the presence of vertical flows and horizontal magnetic fields, line-of-sight tracking methods donot accurately capture the complete footpoint dynamics and the null hypothesis that the velocities inferredby DAVE represent the flux transport velocities may be rejected in favor of the alternative hypothesis ϑ = v h .The implementation of vector magnetograms in optical flow methods presents practical challenges. First,the transverse magnetic field components are known to be noisier than the line-of-sight component and thenoise variance will likely change from pixel to pixel due to variable photon statistics (heteroscedastic errors).Second the line-of-sight component and transverse components are determined from different polarizationsand require inter-calibration. Third, the orientation of the transverse component is ambiguous by 180 ◦ . Thefirst issue may be addressed within the total least squares framework discussed by Schuck (2006), BranhamJr. (1999) and others (See references in Schuck 2006). The main obstacle to resolving the first issue isestimating a covariance matrix for the structure tensor (cid:104) S (cid:105) .The second and third issues both may be interpreted as inter-calibration bias where the estimatedhorizontal magnetic field (cid:98) B h = α B h is proportional to the true horizontal magnetic field B h ; these errors are This does not imply that the alternative hypothesis is correct. D´emoulin & Berger (2003) terms these fluxes the “tangential fluxes” but “horizontal” is more appropriate in the contextof Welsch’s terminology used in this paper.
27 –Fig. 13.— Scatter plots of the estimated velocities ϑ B z from DAVE versus the horizontal plasma veloc-ities V h B z from ANMHD. The nonparametric Spearman rank-order correlation coefficients ( ρ ), Pearsoncorrelation coefficients ( C ), and slopes ( S ) estimated by the least absolute deviation method are shown.Fig. 14.— ( left ) Scatter plots of the estimated Poynting flux (red) and helicity flux (blue) from DAVEassuming ϑ = u versus the Poynting and helicity flux from ANMHD. ( right ) Scatter plots of the estimatedPoynting flux (red) and helicity flux (blue) combining DAVE assuming ϑ = v h with the emergence termfrom ANMHD versus the Poynting and helicity flux from ANMHD. The nonparametric Spearman rank-order correlation coefficients ( ρ ), Pearson correlation coefficients ( C ), and slopes ( S ) estimated by the leastabsolute deviation method are shown. 28 –not random. The flux transport velocities estimated from DAVE4VM are robust to overall inter-calibrationerrors which include the 180 ◦ ambiguity resolution errors. Changing the overall magnitude or sign of B h has no effect on the flux transport velocities because (5a) is invariant with respect to the transformation v z → v z /α and B h → α B h (private communication with Pascal Demoulin). However, the estimated verticalperpendicular plasma velocity and vertical perpendicular electric field will be anti-correlated with the groundtruth when α <
0. Nonetheless, an overall rescaling of the horizontal magnetic field will have no effect onthe helicity flux. However, the Poynting flux will be incorrect by a factor of α including perhaps a signerror because of the rescaling horizontal magnetic field which is inherent in the energy estimate (25). Moretroublesome are the effects of spatially varying bias errors in inter-calibration or ambiguity resolution. Theconsequences of these errors, particularly along the boundaries between proper and improper ambiguityresolution, are presently unknown and should be investigated with future end-to-end analysis of syntheticmagnetograms. However, local methods such as DAVE4VM are probably more robust than global methodsto spatially dependent errors in inter-calibration or ambiguity resolution because local methods inherentlylocalize the effect of bias errors by isolating subregions with the window aperture whereas global methodscouple the entire solution region together permitting bias errors in one subregion to influence the solution inother subregions.In light of the DAVE4VM’s dramatic improvement in performance by simply including horizontal mag-netic fields, speculation that the ANMHD simulation data are not appropriate for testing tracking methodscannot be correct. Rather, aside from issues of image structure, the ANMHD simulation data represent an ideal case for the line-of-sight methods because the vertical magnetic field is known (not simply the line-of-sight component). The results of this study suggest that horizontal magnetic fields and vertical flows willrender velocity estimates from “pure” tracking methods inaccurate if they are treated as the flux transportvelocities u . This conclusion holds equally true for velocity estimates near disk center as the ANMHD simula-tions represents disk-center data! The good agreement between the performance of MEF and the DAVE4VMon the ANMHD data implies that incorporating the right physics is more important for producing accuratevelocity estimates than is the particular method used the solve the equations.Presently, the only way to explore the “image” physics is by testing the “optical flow” methods onsynthetic data from well-designed MHD simulations that attempt to reproduce the physics of the Sun. Naive“moving paint” experiments (Schuck 2006) cannot critically test optical flow methods for magnetogramsbecause the test data are consistent with the two circumstances when pure tracking methods will certainlyperform well: B h = 0 and v z = 0. Consequently, good performance of an optical flow method in naive“moving paint” experiments should not be considered evidence that a method will produce accurate estimatesof plasma physics quantities.In the interest of reproducibility (Joyner & Stein 2007), all of the software used to perform the calcula-tions, create the figures, and draw the conclusions for this paper are archived with the Astrophysical Journal as a tgz file. Updates to the DAVE/DAVE4VM software are also available. I thank the referee for constructive criticism, Graham Barnes for encouraging the publication of thiswork, Bill Abbett for providing the ANMHD data that formed the core of this research, and Pascal Demoulinand Brian Welsch for encouraging the clarification of several issues discussed in this manuscript. I alsogratefully acknowledge useful conversations with George Fisher, Bill Amatucci, and Etienne Pariat. I thankJulie Schuck for editing the manuscript. This work was supported by NASA LWS TR&T grant NNH06AD87I,
29 –LWS TR&T Strategic Capability grant NNH07AG26I, and ONR.
A. Matrix Elements of S G = ( ∂ x B z ) G = ( ∂ x B z ) ( ∂ y B z ) G = ( ∂ y B z ) G = B z ( ∂ x B z ) + ( ∂ x B z ) x (cid:48) G = B z ( ∂ y B z ) + ( ∂ x B z ) ( ∂ y B z ) x (cid:48) G = B z + 2 B z ( ∂ x B z ) x (cid:48) + ( ∂ x B z ) x (cid:48) G = B z ( ∂ x B z ) + ( ∂ x B z ) ( ∂ y B z ) y (cid:48) G = B z ( ∂ y B z ) + ( ∂ y B z ) y (cid:48) G = B z + B z ( ∂ x B z ) x (cid:48) + B z ( ∂ y B z ) y (cid:48) + ( ∂ x B z ) ( ∂ y B z ) x (cid:48) y (cid:48) G = B z + 2 B z ( ∂ y B z ) y (cid:48) + ( ∂ y B z ) y (cid:48) G = ( ∂ x B z ) y (cid:48) G = ( ∂ x B z ) ( ∂ y B z ) y (cid:48) G = B z ( ∂ x B z ) y (cid:48) + ( ∂ x B z ) x (cid:48) y (cid:48) G = B z ( ∂ x B z ) y (cid:48) + ( ∂ x B z ) ( ∂ y B z ) y (cid:48) G = ( ∂ x B z ) y (cid:48) G = ( ∂ x B z ) ( ∂ y B z ) x (cid:48) G = ( ∂ y B z ) x (cid:48) G = B z ( ∂ y B z ) x (cid:48) + ( ∂ x B z ) ( ∂ y B z ) x (cid:48) G = B z ( ∂ y B z ) x (cid:48) + ( ∂ y B z ) x (cid:48) y (cid:48) G = ( ∂ x B z ) ( ∂ y B z ) x (cid:48) y (cid:48) G = ( ∂ y B z ) x (cid:48) s = − ( ∂ x B x ) ( ∂ x B z ) − ( ∂ y B z ) ( ∂ x B z ) s = − ( ∂ x B x ) ( ∂ y B z ) − ( ∂ y B z ) ( ∂ y B z ) s = − ( ∂ x B x ) B z − ( ∂ y B z ) B z − ( ∂ x B x ) ( ∂ x B z ) x (cid:48) − ( ∂ y B z ) ( ∂ x B z ) x (cid:48) s = − ( ∂ x B x ) B z − ( ∂ y B z ) B z − ( ∂ x B x ) ( ∂ y B z ) y (cid:48) − ( ∂ y B z ) ( ∂ y B z ) y (cid:48) s = − ( ∂ x B x ) ( ∂ x B z ) y (cid:48) − ( ∂ y B z ) ( ∂ x B z ) y (cid:48) s = − ( ∂ x B x ) ( ∂ y B z ) x (cid:48) − ( ∂ y B z ) ( ∂ y B z ) x (cid:48) s = ( ∂ x B x ) + 2 ( ∂ x B x ) ( ∂ y B z ) + ( ∂ y B z ) s = − B x ( ∂ x B z ) − ( ∂ x B x ) ( ∂ x B z ) x (cid:48) − ( ∂ y B z ) ( ∂ x B z ) x (cid:48) s = − B x ( ∂ y B z ) − ( ∂ x B x ) ( ∂ y B z ) x (cid:48) − ( ∂ y B z ) ( ∂ y B z ) x (cid:48) s = − B x B z − ( ∂ x B x ) B z x (cid:48) − ( ∂ y B z ) B z x (cid:48) − B x ( ∂ x B z ) x (cid:48) − ( ∂ x B x ) ( ∂ x B z ) x (cid:48) − ( ∂ y B z ) ( ∂ x B z ) x (cid:48) s = − B x B z − ( ∂ x B x ) B z x (cid:48) − ( ∂ y B z ) B z x (cid:48) − B x ( ∂ y B z ) y (cid:48) − ( ∂ x B x ) ( ∂ y B z ) x (cid:48) y (cid:48) − ( ∂ y B z ) ( ∂ y B z ) x (cid:48) y (cid:48)