A Comprehensive Method of Estimating Electric Fields from Vector Magnetic Field and Doppler Measurements
Photospheric electric fields, estimated from sequences of vector magnetic field and Doppler measurements, can be used to estimate the flux of magnetic energy (the Poynting flux) into the corona and as time-dependent boundary conditions for dynamic models of the coronal magnetic field. We have modified and extended an existing method to estimate photospheric electric fields that combines a poloidal-toroidal (PTD) decomposition of the evolving magnetic field vector with Doppler and horizontal plasma velocities. Our current, more comprehensive method, which we dub the “PTD-Doppler-FLCT Ideal” (PDFI) technique, can now incorporate Doppler velocities from non-normal viewing angles. It uses the FISHPACK software package to solve several two-dimensional Poisson equations, a faster and more robust approach than our previous implementations. Here, we describe systematic, quantitative tests of the accuracy and robustness of the PDFI technique using synthetic data from anelastic MHD (ANMHD) simulations, which have been used in similar tests in the past. We find that the PDFI method has less than error in the total Poynting flux and a error in the helicity flux rate at a normal viewing angle ) and less than and errors respectively at large viewing angles (). We compare our results with other inversion methods at zero viewing angle, and find that our method’s estimates of the fluxes of magnetic energy and helicity are comparable to or more accurate than other methods. We also discuss the limitations of the PDFI method and its uncertainties.
- 1 INTRODUCTION
2 FINDING ELECTRIC FIELDS
- 2.1 The Inductive PTD Electric Field:
- 2.2 Enforcing the Condition with the potential
- 2.3 Non-inductive Doppler and Fourier Local Correlation Tracking (FLCT) electric fields: and
- 3 USING ELECTRIC FIELDS TO COMPUTE POYNTING AND HELICITY FLUXES
- 4 VALIDATION OF THE ELECTRIC FIELD INVERSION METHOD USING ANMHD SIMULATIONS
- 5 DISCUSSION AND CONCLUSIONS
- A Generalization Of PTD Formalism to Spherical Coordinates
Energy from magnetic fields on the Sun powers nearly all manifestations of solar activity, from heating in active regions and the solar wind, to dramatic events like flares and coronal mass ejections (CMEs). Quantitative studies of the flow of magnetic energy through the solar atmosphere require a knowledge of both the magnetic and electric field vectors. The goal of this paper is to show how time sequences of vector magnetic field maps (vector magnetograms), along with other observational constraints such as Doppler flow measurements, can be used in a practical way to construct maps of the electric field vector. Knowing both the electric and magnetic field, we can find the Poynting flux of electromagnetic energy and the flux of relative magnetic helicity, important quantities that describe how magnetic energy is transported, stored, and released in the solar atmosphere.
Past attempts to determine the electric field distribution on the Sun have followed two approaches: (1) direct spectroscopic measurements of the electric field using the linear Stark effect, (2) indirect determinations of the electric field by using Faraday’s law, relating the temporal derivative of the magnetic field to the curl of the electric field. Wien (1916) was the first to suggest measuring the electric field of solar plasma directly using the Stark effect. In the 1980s, some attempts were made to measure the electric field at the Sun using helium and silicon spectra, finding an electric field of V/cm (Davis, 1977) and V/cm (Jordan et al., 1980). Later, examining Stark-broadened Paschen emission lines in hydrogen, Moran and Foukal (1991) found an upper limit of V/cm. In the same work they pointed out that the direct measurement of the electric field using the Stark effect is difficult because of the low sensitivity of the measurements.
Indirect determinations of the electric field have been considerably more successful. The implementation of local-correlation tracking techniques (LCT, e.g. November and Simon 1988; Berger et al. 1995; Hurlburt et al. 1995; Fisher and Welsch 2008), applied to time-sequences of magnetograms allows one to determine a “pattern motion” velocity of the line-of-sight component of the magnetic field in the plane of the magnetogram (Démoulin and Berger, 2003). Interpreting the pattern motion as a plasma velocity, and assuming an ideal MHD Ohm’s law for the electric field (), allows one to determine some of the information about the velocity field, and hence by assumption, the electric field. These techniques can be improved by explicitly incorporating the normal component of the magnetic induction equation into the solution for the components of the velocity (Kusano et al., 2002; Welsch et al., 2004; Longcope, 2004; Schuck, 2006, 2008; Chae and Sakurai, 2008), resulting in improvements in the velocity inversion (Welsch et al., 2007).
More recently, we have developed inversion techniques for the electric field itself rather than the velocity field, based on a “Poloidal-Toroidal Decomposition” (PTD) of the magnetic field and its time derivative. These techniques use Faraday’s law and other theoretical and observational constraints to determine the electric field (Fisher et al., 2010, 2012). Inversions for the electric field instead of the velocity field offer a few distinct advantages. First, the PTD method incorporates additional information from evolution of the horizontal magnetic field (the field parallel to the photosphere). This is used to make an estimate of the normal electric field, independent of the horizontal electric field inferred by methods that derive from the normal component of the induction equation. Second, in regions where the magnetic field is relatively weak or uncertain, the determination of the velocity field is especially ill-posed: for instance, outside of active regions, fluctuations in the normal magnetic field can be dominated by noise, resulting in a wildly varying and unphysical behavior in the inferred . Complex masking and filtering techniques must be applied post facto to suppress this behavior. The solutions for the electric field, on the other hand, may contain small errors in regions with small magnetic field strength, but the solutions vary smoothly within regions containing both strong and weak magnetic fields. Third, if a model for non-ideal effects is included in the electric field assumed to drive the change in , they can be captured by the PTD solutions, while the velocity formalism may give spurious results.
Fisher et al. (2010) described in detail the PTD techniques necessary to find the “inductive” electric field that satisfies Faraday’s law, given the spatial distribution of the temporal derivative of the magnetic field measured across a region of the solar photosphere. They also noted that the resulting solutions are not unique, in that gradients of arbitrary scalar functions could be added to the PTD solutions without affecting Faraday’s law. They presented several approaches to specify the electric field uniquely, including requiring the electric field to be perpendicular to the magnetic field, consistent with the ideal Ohm’s law. Fisher et al. (2012) demonstrated that by adding to the PTD solutions the gradients of scalar functions that are constructed to match the electric field near polarity inversion lines (PILs), the accuracy of electric field reconstructions could be substantially improved, beyond all of the velocity inversion techniques considered in the comparative study of Welsch et al. (2007). To match the electric field near PILs, one can use measurements of the Doppler velocity and the transverse magnetic field distribution in those regions.
In this paper, we build on Fisher et al. (2012) results in a number of ways, with the goal of carefully describing and validating a practical implementation of the introduced methods that could be used to routinely find the electric field from sequences of vector magnetogram observations, such as those available from the Helioseismic and Magnetic Imager (HMI) (Schou et al., 2012; Scherrer et al., 2012) on NASA’s SDO mission. First, we adopt the use of the FISHPACK software library for Helmholtz equations (Schwarztrauber, 1975), developed at NCAR in the 1970s, for solving the two-dimensional Poisson equations that are at the heart of the PTD inversion technique. This software, based on the cyclic reduction technique, is very efficient, and explicitly allows for the specification of normal-derivative (Neumann) boundary conditions that are consistent with our desired boundary conditions on the electric field at the edges of vector magnetogram images. Second, we develop the ability to incorporate line-of-sight Doppler measurements that are taken at non-normal viewing angles, appropriate for vector magnetograms that are not at disk center, to determine electric fields near PILs. Finally, to validate our inversion methods, we analyze the performance of our techniques, comparing electric field inversions with a test case from an MHD simulation where the true electric field is known (Welsch et al., 2007). Besides the electric field, we compare derived values of the Poynting flux and the relative helicity fluxes with the known values from the test case. Finally, in Appendix A we describe how we adapt our electric field solution techniques to spherical polar coordinates, appropriate for larger fields of view on the Sun.
The remainder of the paper is structured as follows. In § 2 we describe the electric field inversion technique, starting from the basic PTD formalism, and also improvements that account for a non-zero viewing angle. In § 3 we describe how to calculate Poynting and helicity fluxes, using our derived electric field. This includes a decomposition of the Poynting flux into fluxes of potential and free magnetic energies. In § 4 we perform tests of the electric field inversion method using the ANMHD simulation test case, and provide validation metrics that quantify how well the inversion recovers the true solution. Finally, in § 5, we summarize the strengths and weaknesses of our technique and describe how it will be used to analyze HMI vector magnetograms.
2 Finding Electric Fields
2.1 The Inductive PTD Electric Field:
Poloidal-toroidal decomposition (PTD) allows one to estimate the photospheric electric field vector from the evolution of the photospheric magnetic field vector. Here we present a brief synopsis of the PTD method (Fisher et al., 2010), plus improvements we have made since that article was published.
The fundamental idea of PTD is that the magnetic field vector has a solenoidal nature and hence can be specified by two scalar functions, the “poloidal” and “toroidal” potentials (Chandrasekhar, 1961; Moffatt, 1978), which we denote and following the notation used in descriptions of anelastic MHD which employ these potentials (Lantz and Fan, 1999):
Here we consider a locally Cartesian coordinate system , that has its axis oriented perpendicular to the photosphere, with the positive direction away from the Sun’s center (see Figure 1; the spherical case is described in Appendix A). We use subscript to denote vector components or derivatives in the vertical direction, and subscript to denote vector components or derivatives in the horizontal directions, i.e. parallel to the photosphere.
Taking the partial time derivative of Eq. (1), one finds
By examining the -component of Eq. (2), its horizontal divergence, and the z-component of its curl, we find three two-dimensional Poisson equations for the unknown functions , and in terms of known physical quantities:
For future use, we can rewrite this result as:
Uncurling the Eq. (6), we derive the electric field in terms of and :
Here, is the purely inductive contribution to the electric field determined by PTD potentials, and . Within this article, will be referred to as P-solution, or the PTD electric field, where superscript “P” stands for the PTD method. At this point, two comments can provide some insight into the nature of : (1) the horizontal components of only depend upon , which was derived from ; and (2) the vertical component of only depends upon , which is derived from . The component of the total electric field arising from is the non-inductive contribution, for which the PTD solution to Faraday’s law, Eq. (7), reveals no information. For the simulation data analyzed by Welsch et al. (2007), Fisher et al. (2010) demonstrated that when , when non-inductive contributions to the electric field are ignored, does a poor job of representing the actual electric field. The key to a more accurate reconstruction of the electric field is incorporating other physical constraints, and/or additional observational information into the solutions for . Procedures for doing this will be described further below. But we first complete our discussion for obtaining the full (PTD) solution.
To find and and hence the solution, we must solve the three Poisson equations (3)-(5) with well-posed boundary conditions. To permit a net gradient in each potential over the field of view, we choose Neumann boundary conditions
where subscript denotes components or derivatives in the direction of the outward normal to magnetogram’s boundary, and subscript denotes components or derivatives in the counter-clockwise direction along the magnetogram boundary. The first boundary condition, Eq. (9), implies that the tangential component of the electric field around the magnetogram boundary vanishes, implying that the average value of within the magnetogram is zero. If the average value of is not zero (i.e., a change in the flux balance, or monopole term), then we can add a correction term to the electric field post-facto (see Appendix C of Fisher et al. (2010)). The second and third boundary conditions, Eq. (10)-(11), derived from evaluating Eq. (7) at the magnetogram boundary, are degenerate, i.e. there is a family of coupled non-zero solutions, and , which correspond to zero time derivative of the horizontal magnetic field:
Since the solutions to these Cauchy-Riemann equations each satisfy the two-dimensional Laplace equation, they can be added to solutions of equations (4) and (5) without changing the time derivative of the horizontal field on the boundary. This means that there is some freedom to specify the solutions of or at the boundary. In practice, this means that one could choose to set the derivative parallel to the boundary of one of these two functions to zero, while still obeying the coupled boundary conditions (10) and (11). To remove the coupling between the boundary conditions (10) and (11), we choose along the magnetogram boundary, meaning that is assumed uniform along the magnetogram boundary. (This implies no change in the net signed vertical current through the field of view.) This also allows us to set the value of for to zero at the magnetogram boundary by simply subtracting the boundary value from the entire solution for . Choosing the condition also means that the solution of the Poisson equations must be done in a certain order. First, the solution for must be obtained, then the solution for can be obtained, using the variation of along the magnetogram boundary to specify the normal derivative of . The solution for is independent of the other two solutions, and so can be obtained either before or after obtaining the other two solutions.
When solving Poisson equations (3)-(5) and boundary conditions (9-11) numerically, it is important to implement a fast and robust numerical scheme, especially when dealing with large magnetogram datasets. In the past, Fisher et al. (2012, 2010) used the Newton-Krylov technique adapted from the RADMHD code (Abbett, 2007). Its main disadvantage is slow computational speed and poor robustness when applied to this particular problem. In this article, we adopt FISHPACK, a fast and robust collection of Fortran subprograms developed at NCAR (Schwarztrauber, 1975). FISHPACK applies the cyclic reduction technique (Sweet, 1974) and the standard five-point finite difference approximation for the Laplacian to solve the two-dimensional Helmholz equation in Cartesian or spherical (Appendix A) coordinates. This software was designed to solve Helmholz equations (the Poisson equation is a special case of the Helmholz equation) with Neumann, Dirichlet, or periodic boundary conditions, and can use either a centered or staggered grid. To convert our equations into a form that is compatible with the FISHPACK subroutine HWSCRT, we use centered, 2nd-order accurate finite differences to approximate first derivatives in and , and the centered 5-point expression for the Laplacian of a function in a Cartesian, two dimensional geometry:
where in the last equation we assume . Similar to spatial differences, to evaluate time-derivative of function at time , , the source terms (the left-hand sides) of Eq. (3)-(5), we use the centered finite difference approximation
where is the time step between two consecutive frames . Once a solution for has been obtained, the electric field components and are computed by using the finite derivative approximations (Eq. (14-15)) in the expression for ; for we use the solution . We add a single layer of ghost zones around the periphery of each solution domain to ensure that the Neumann boundary conditions are obeyed when the above finite difference expressions are used. Given the known Neumann (normal first derivative) boundary conditions (Eq. (9-11)) and the interior solution points that are returned from HWSCRT, using the expressions in Eq. (14-15), we then determine the values of the solution variables in the ghost zone layers. This procedure leaves the four corners of the ghost zone values undefined; for cosmetic reasons when displaying the solution, we assume each corner value is the average of its two nearest-neighbor ghost zone values.
As Schwarztrauber (1975) points out, it is possible to specify Neumann boundary conditions that are incompatible with some Poisson equations. Their approach to this quandary is to add a constant value to the source term of the Poisson equation, with that value varied until the best possible fit to the boundary values can be obtained. This offset value is returned by HWSCRT subroutine as the PERTRB variable, which we record along with the solution itself. The value of PERTRB is then monitored to ensure that it remains small compared to typical values of the source term of our Poisson equations.
Finite difference approximations for the first derivative (Eq. (14-15)) and the Laplacian (Eq. (16)) introduce incompatibilities, since the -point Laplacian implicitly assumes that first derivatives are evaluated at half-integer grid locations. These incompatibilities could be eliminated by using a staggered grid, in which electric field variables are defined at cell edges, while magnetic field values are defined at cell centers. However, we have found that the resulting errors are usually small, and the introduction of a staggered grid introduces a number of complications, when adding the electric fields derived from the gradients of scalar potentials to the solutions. Here, we have decided to accept these incompatibilities as a cost of using a centered grid, where the electric field and magnetic field variables are assumed co-spatial. We may adopt a staggered grid in a future version of the software.
One case where the incompatibility errors may not be small, is if the solutions are dominated by high-frequency noise near the magnetogram boundaries because of poor signal-to-noise ratios in areas of weak magnetic fields. In this case, we have found it useful to add “zero-padding” to the solution domain for a modest number of zones around the periphery of the observed magnetograms, in which all magnetic field values and their time derivatives are set to zero, forcing the source terms of all of the Poisson equations to be zero within a thin ribbon inward of the computational boundary. The smoothness of the resulting Poisson equation solutions ensures that the incompatibilities of the finite difference approximations described above are not important in the boundary regions, and therefore that the finite difference expressions for the Neumann boundary conditions are good approximations for the normal first derivatives.
2.2 Enforcing the Condition with the potential
One of the properties of the derived inductive PTD electric field (Eq. (8)) is that for most of the cases and even for ideal MHD case, , that has perpendicular to , the derived is not perpendicular to . To resolve this inconsistency, in this section we review the “iterative” method described in §3.2 of Fisher et al. (2010) plus the changes that we have made to it since then. This technique will be used extensively in the discussions of other contributions to the electric field below.
The goal of the iterative method is to construct the gradient of an electric potential that obeys the constraint
for all points within the vector magnetogram domain. Once found, can then be subtracted from to yield an electric field which is normal to , but which has the same curl as . In practice, and as described in §3.2 of Fisher et al. 2010, the procedure is to find two two-dimensional functions, and over the 2-d domain of the vector magnetogram, such that
where and are respectively the horizontal and vertical components of the unit vector that points in the direction of . Once and have been found, then we define the quantity
where then obeys both Faraday’s law and the ideal constraint .
We continue to follow the iterative method procedures described in §3.2 of Fisher et al. (2010) to determine and , but with the following changes: First, in Step 4 of that procedure, in which a horizontal Poisson equation is solved to update the iterative guess for , we now use the FISHPACK routine HWSCRT to solve the Poisson equation, instead of using FFT techniques that implicitly assume periodic boundary conditions. We assume homogenous Neumann boundary conditions for , such that on the outer boundary, meaning that the normal component of the horizontal electric field due to the scalar potential is zero on the outer boundary of the vector magnetogram. We use equations (14-15) to approximate the components of the horizontal gradient of , and equation (16) to approximate the Laplacian. Ghost zone values are defined in a similar fashion as in the previous section for determining . The solution for , in Step 5 of the iterative procedure, is usually negligibly small at the outer boundary of the magnetogram, provided that the horizontal components are also small there. This results in the desired property that is essentially zero at the outer magnetogram boundary. An additional change we have made from Fisher et al. (2010) is that instead of evaluating an error criterion to end the iteration sequence, we have found it more efficient to simply specify a maximum number of iterations, typically chosen to be 25 (see discussion in §4.1).
An important point about the iterative method is that it can be applied not only to the PTD solution itself, but to any solution for the electric field which contains components parallel to , that one might like to minimize. In particular, we can apply the iterative method to solutions, which include both the PTD solution and additional electric field contributions from scalar potentials determined from Doppler measurements or correlation tracking results, as are described in further detail below. Once these interim solutions are obtained, the iterative method can then be applied as a final step to generate the scalar potential (and ) needed to satisfy .
Finally, we comment on the uniqueness of the solutions for found from the iterative method, given a fixed input electric field . As noted in Fisher et al. (2010), mathematically the constraint (18) does not result in a unique solution for and . Empirically, however, we find that the iterative procedure of §3.2 of Fisher et al. (2010) (and also with the modifications described here) consistently results in the same solution for and , even when initialized with substantially different trial functions. We surmise that applying the iterative procedure itself imposes some additional constraints beyond (18) which we do not fully understand, but which then result in a unique solution for the input . In fact, however, given any that is consistent with both Faraday’s and Ohm’s laws, an additional potential functions can be constructed that can be superposed without violating either constraint. So, while our solution method returns a unique answer, Faraday’s and Ohm’s laws do not, by themselves, fully constrain .
2.3 Non-inductive Doppler and Fourier Local Correlation Tracking (FLCT) electric fields: and
When examining the first and second rows of Figure 2 of Fisher et al. (2012), comparing the actual components of the test simulation electric field with the same components of , as well as the actual and inverted vertical Poynting flux, one of the most obvious discrepancies is the far lower range of the Poynting flux values originating from the electric field inversion as compared to the true values (see Figure 3 in Fisher et al. (2012)). The and components of also show significant discrepancies, with the derived horizontal components of being smaller in magnitude than the actual values. Is there some additional information about the horizontal electric field that can be determined from other observational data besides the magnetic field?
Estimating the transport of magnetic energy across the photosphere requires accurately reconstructing photospheric electric fields during flux emergence. From studying ANMHD simulations of an emerging bipolar magnetic configuration, we have found that, by itself, Faraday’s law applied to does not completely capture significant aspects of horizontal electric fields present during magnetic flux emergence. Here we review two possible sources of electric fields which may not be reflected in solutions of Faraday’s law: Electric fields determined from Doppler measurements, and those determined using correlation tracking methods. In the discussion below, we will first consider the computation of these electric fields from the conceptually easier vantage point of zero viewing angle, where line-of-sight is normal to the solar surface (§2.3.1 and §2.3.2). Once we have introduced all of the physical concepts, we will describe the modifications necessary to account for non-normal viewing angles (§2.3.3).
Doppler Electric Fields: Normal Viewing Angle
Certainly in flux-emergence regions where the vertical magnetic field is near zero ( near PILs), a knowledge of the vertical velocity would yield an unambiguous determination of the horizontal electric field:
if we can assume that an ideal MHD Ohm’s law applies in the photosphere. Such measurements could be determined by Doppler velocity measurements viewed from directly above the active region. Figure 1 of Fisher et al. (2012) illustrates with a thought-experiment that such a flux-emergence electric field might even have zero curl near the PIL, and yet, if the Doppler velocity component had a large amplitude along the PIL region, would still result in a very large electric field. One possible approach to incorporating such data would be to simply replace the PTD solution with for those parts of the magnetogram near the PIL location. If one did this, however, the resulting solution near the PIL would not necessarily obey Faraday’s law, and furthermore, could vary discontinuously and/or unphysically, as one crossed the boundary between regions where the two solution methods for were applied.
The solution to this dilemma is to combine the two solutions in such a way that Faraday’s law is obeyed everywhere in the magnetogram region, but that the non-inductive (curl-free) component of is added to , to provide a more realistic electric field estimate in flux emergence regions. To do this, we need two additional ingredients: (1) an estimate for that one can apply when one leaves the vicinity of the PIL, and (2) a technique for removing the inductive contribution from , since will contain all the needed information about Faraday’s law.
Addressing the first of these ingredients, as one moves away from the PIL, Eq. (21) becomes increasingly uncertain for two reasons: first, as the vertical component of the field becomes nonzero, it becomes increasingly possible for any existing flows parallel to the magnetic field to contribute to the Doppler signal without contributing to an electric field. Second, once is no longer zero, there is a contribution to the horizontal electric field from horizontal flows crossed into , which is not reflected in Eq. (21). We adopt the ad-hoc solution of Fisher et al. (2012), namely that is multiplied by an attenuation or confidence factor, , that reduces the amplitude of as the ratio increases above zero. We adopt the same functional form of introduced in Fisher et al. (2012), namely
where is an adjustable parameter that can be tuned by comparing with test case data, where the electric field is known. The determination of the optimal value of for the ANMHD test case is discussed in §4.1.
Addressing the second ingredient, we must find a way of introducing only the curl-free contribution from the electric field estimate . This can be accomplished if we first write the resulting electric field in terms of a scalar potential, , such that
and then set the divergence of to the divergence of , resulting in the Poisson equation
that we solve for using FISHPACK. Since is the gradient of a scalar function, it automatically has zero curl.
When derived from is added to the PTD solution , we denote the resulting solution as the PD electric field. If we then add the ideal correction to the PD solution, by applying the constraint, we denote that resulting electric field as the PDI solution.
FLCT-derived Electric Fields
Local Correlation Tracking (LCT) techniques have been used extensively in recent years to estimate flows in the solar atmosphere, as noted in the introduction. These methods use observed pattern motions of some quantity, such as the specific intensity, or the normal component of the magnetic field, to provide an estimate of the velocity components of the plasma that are parallel to the plane of the image. The input to LCT techniques is the following: two images of the quantity that is being tracked, the time difference between the two images, and an a-priori assumption about the coherence length of the flow-field. The output from LCT is a two- dimensional flow field that warps the first image into the second image over the elapsed time between the images. The fundamental assumption of LCT techniques is that the observed pattern motion corresponds to real physical motion. In the remainder of this paper, we will focus on the use of the Fourier Local Correlation Tracking (FLCT) technique (Welsch et al., 2004; Fisher and Welsch, 2008), but any other optical flow method (Schuck, 2006) could be used instead.
When the quantity being tracked is the normal component of the magnetic field, velocity estimates obtained from LCT methods have been found to provide a good estimate on their own for the “inductive” flows responsible for changes in the normal field (see Figure 6 of Welsch et al. 2007). However, LCT results also contain information about flows that do not contribute directly to the observed changes in the magnetic field tracer, and thus provide additional information about the velocity field (or the electric field, assuming an ideal Ohm’s law). The ILCT method (Welsch et al., 2004) uses this knowledge to explicitly separate the flow field into inductive and non-inductive contributions.
In our case, already contains the “inductive” contribution to the electric field, so a straightforward addition of the electric field from flows derived from FLCT would add redundant information that has already been incorporated, while also adding new information from the non-inductive flows. We therefore want to determine only the non-inductive contribution to the electric field from FLCT-derived flows, and add that result to . The technique for doing this is very similar to how we add the non-inductive electric field determined from the Doppler measurements.
Assuming the ideal Ohm’s law, we can write the electric field derived from FLCT as
where here is the horizontal flow-field from FLCT. This results in a contribution to both (the first term in Eq. (25)), and , the second term in that equation. Considering first the contribution of FLCT to , we note that as one approaches the PIL, this estimate for becomes increasingly suspect since , the quantity being tracked, becomes very small, and the validity of the correlation tracking paradigm to describe the evolution of becomes increasingly questionable. Since the FLCT results get unreliable precisely where we think the Doppler electric field estimate becomes reliable, we make the Ansatz that the FLCT-derived electric field should be multiplied by the complement () of the same confidence factor, , used to modulate the Doppler electric field.
To include only the non-inductive contribution, we write the electric field in terms of the gradient of a scalar potential,
This electric field has no curl, and hence no inductive contribution. To determine , we multiply the right hand side of Eq. (25) by , take its horizontal divergence, and then set that equal to the divergence of Eq. (26):
The result of this solution is an electric field that includes the non-inductive contributions from FLCT which can then be added to . This procedure is very similar to ILCT (Welsch et al., 2004), apart from (1) using electric fields instead of flow-fields and (2) the use of the confidence factor.
Now, we consider the contributions to from the second term in equation (25), and how to correct this for the inductive contribution to from . First, note that the FLCT estimate for is approximately correct, but has a lot of scatter (see the second panel of Figure 11 in (Welsch et al., 2007)). Our own comparison of from FLCT with that from shows that the FLCT contribution contains significant information that overlaps with . Since the inductive contribution is already completely contained within , we need to remove any inductive signature from the FLCT estimate for . The result of doing this is a residual contribution to which satisfies the horizontal Laplace equation, but obeys the Dirichlet boundary conditions for that are returned from FLCT. If the boundary electric fields are near zero, as we generally desire, then the resulting FLCT residual contribution is generally small, since solutions of Laplace’s equation achieve their extremum values at the boundaries. The residual contribution is certainly very small in the validation test case considered later in this paper, as well as for our planned HMI inversions, and we therefore do not include it in our inversion estimates at this time.
When the solution is added to the PTD electric field , we denote the resulting electric field as the PF solution. If the ideal electric field is added to that solution to impose the constraint , we denote the resulting electric field as the PFI solution. If the FLCT electric field and the Doppler electric field are added to , we denote the resulting electric field solution as the PDF solution. Finally, if we impose the ideal electric field constraint to the PDF solution, we obtain what we denote as the PDFI solution. As discussed in greater detail later in this article, we generally find that the PDFI solution provides the best results when compared to the electric fields in our standard test case.
Doppler Electric Fields: Non-normal Viewing Angles
Viewing the solar surface from a non-normal angle ( ) introduces some complications for estimating the electric field. These are due to (1) foreshortening effects, and (2) changes in the direction of the flows derived from Doppler measurements. Because our assumed magnetic field data is from vector magnetograms, we can: correct it for foreshortening; interpolate it to a uniform grid on the solar surface; and reproject the magnetic field directions into the locally normal and horizontal components. This means that the PTD inversion technique for can be applied to the remapped data without any change in the method, as can FLCT, to derive part of the non-inductive the electric field. Using the Doppler velocity to compute non-inductive electric fields along the solar surface, however, becomes more difficult, and some additional steps must be taken.
We first note that if the angle of the observation is written in terms of the unit line-of-sight (LOS) vector , where the direction of this vector is toward the observer as measured from the solar surface, then the velocity vector defined by the Doppler shift is given by
where is the plasma velocity at the source point on the solar surface. The geometry is shown in Figure 1.
The electric field due to the Doppler velocity near LOS PILs is then given by
where here represents the components of in the directions transverse to the LOS direction. Note that this electric field is no longer confined to directions parallel to the solar surface, but instead is confined to the directions normal to . Similar to the normal-viewing angle case, we anticipate that as the ratio increases above zero, our confidence in this electric field to represent the transverse electric field decreases due to possible parallel flows and other contributions to the transverse electric field. Therefore we multiply by an ad-hoc confidence factor
where has the same value as in Eq. (22) for . Note that the confidence factor used for determining is and not , since the FLCT formalism is not affected by non-normal viewing angles.
To determine non-inductive contributions from the Doppler electric field, we borrow the iterative technique from §2.2, and find a scalar potential function that obeys the constraint
where is the unit vector pointing in the direction of . Note the analogy between the roles of and in Eq. (31) and and in Eq. (18). Accordingly, as in Eq. (19), we then decompose Eq. (31) into the directions parallel and perpendicular to the solar surface:
where and represent the horizontal and vertical components of , respectively. We then solve Eq. (32) for and , using exactly the same “iterative” procedure described in §2.2. The result is a curl-free electric field which can then be added to . Henceforth, even if the viewing angle is normal to the solar surface, we use the non-normal formalism for computing electric fields from Doppler shifts, but then set in expression for .
Combinations of Electric Field Contributions: A Summary
Now that we have introduced the different possible electric field contributions, we summarize the naming convention we use henceforth in this article, depending on which contributions are included in an electric-field estimate (see Table 1). To distinguish between different versions of electric field inversions we use combinations of the following four letters: “P” (for PTD, inductive), “D” (for Doppler, non-inductive), “F” (for FLCT, non-inductive), and “I” (for ideal, non-inductive). Solutions with an “I” in their name have the ideal-MHD constraint, the addition of the field from the potential , imposed as the last step. We also consider the full FLCT results, denoted FI, the full Doppler electric field results, denoted DI, and their sum, denoted DFI. For FI, DI and DFI we did not impose , because their total electric field is already ideal. Consequently, we have a notation for the electric field, helicity and Poynting fluxes computed with a variety of different electric field estimation techniques.
|Name||Denoted||Equation for||Input Data / Constraints|
||, (generally, )|
||, , , FLCT output|
||, , , Dopp. data|
|PDFI solution||, , , Dopp. & FLCT|
||Dopp. & FLCT data|
The biggest changes from the solutions described here and those of Fisher et al. (2012) are the adoption of the FISHPACK software to solve the two-dimensional Poisson equations, the ability to compute contributions to Doppler-shift electric fields from non-normal viewing angles, and a much more systematic and quantitative testing of the accuracy and robustness of the technique (described further in §4).
3 Using Electric Fields to Compute Poynting and Helicity Fluxes
3.1 Poynting Flux Components: S and S
The Poynting-flux vector
measures the flow of electromagnetic energy at the photosphere, where magnetic field vector is defined from the observations, and the electric field vector is taken from the techniques described in §2. Since we are interested in the amount of energy flowing into and out of the corona, we focus most of our attention on the vertical component of Poynting flux, . This depends upon the horizontal components of both the electric field and the magnetic field:
Computing is straightforward, given our techniques for finding all three components of as described above, and the availability of the components of from the vector magnetogram data. But we can go further, and decompose into two contributions, the flux of potential-field energy, and the flux of free magnetic energy. The distinction between these two contributions was first made by Welsch (2006), who expressed the Poynting flux in terms of velocities. An expression for the flux of free magnetic energy was given in Fisher et al. (2010) in terms of the electric field, and a decomposition of the horizontal magnetic field in terms of the poloidal and toroidal potentials. The basic idea is that the horizontal magnetic field can be divided into a potential-field contribution, and a contribution due to currents that flow into the atmosphere from the photosphere. Because there are several ways that one can construct a potential field from the photospheric vector magnetogram data, we first discuss what we believe is the best way of performing this decomposition, and then briefly mention another alternative we have considered.
Horizontal electric fields derived from Faraday’s law are derived, in part, from the observed evolution of the vertical magnetic field. Therefore, the most self-consistent description of the potential magnetic field is to find the potential field that best matches the values of the measured normal component of the magnetic field at the photosphere. There are a number of ways to do this using, Green’s function techniques, or FFTs (though in the latter case, one must insure the region is flux balanced, and that sufficient padding is provided outside of the active region to mitigate against artifacts from periodic boundary conditions). Once the potential function has been derived, the horizontal potential magnetic field at the photosphere can be evaluated as
The flux of potential-field energy is then given by
To compute the flux of free magnetic energy, we use the residual horizontal field that is the difference between the measured and potential components:
We have also experimented with using the potential field formulation suggested in Appendix A of Fisher et al. (2010). This has the advantage of being extremely easy to compute in terms of the poloidal-toroidal potential function , which has the same PTD formalism we use to solve for . However, as Welsch and Fisher (2014) have shown, this potential field matches the observed values of at the photosphere, and may not accurately reflect the measured value of . If the actual magnetic field were in fact potential, Welsch and Fisher (2014) show that the two formulations would provide the same result. Since currents are present, however, the two formulations can and do show differences in at the photosphere. Since our derived electric fields are driven primarily by changes in , we therefore choose a potential field that explicitly matches at the photosphere.
For short sequences of vector magnetograms, like e.g. the analysis of the ANMHD test data described in §4, we compute the potential magnetic field by using the Green’s function technique (see Appendix B of Welsch and Fisher 2014), as we have found it to be the most robust and accurate. It is, however, very computationally intensive, and can be impractical for large datasets. For long cadences of large vector magnetograms, we use a padded FFT solution for the potential magnetic field, with the size of the needed padding area calibrated by solutions from the Green’s function technique.
3.2 Helicity Flux Rate:
The magnetic helicity of magnetic field in a volume
measures the linkage of magnetic field lines and is useful for describing their topology. Magnetic helicity is physically meaningful only for magnetic fields that are fully contained in a volume , a condition which for active regions with magnetic field lines extending far above and below the photosphere, is not satisfied. In this case, subtracting the helicity of the potential field that matches the normal component of on , one can define topologically meaningful and gauge-invariant relative magnetic helicity, i.e. helicity relative to the potential field (Berger and Field, 1984; Finn and Antonsen, 1985)
Due to absence of accurate magnetic field measurements throughout the corona, it is hard to derive relative magnetic helicity directly from the observations. Instead, rate of change of relative helicity, is used to characterize the coronal magnetic field complexity. If we know the electric field, then the time-rate-of-change of relative helicity in volume — equivalent to the flux of helicity into — is given by the surface integral (Eq. (62) in Berger and Field (1984))
where is the vector potential that generates the potential field in , which matches the photospheric normal field . We require that obey the Coulomb gauge condition (divergence-free in ) and be tangential on the bounding surface. For a closed magnetic field rooted in the photosphere, the helicity flow across the upper boundaries is zero and only photospheric surface term remains.
The first term corresponds to magnetic helicity generated by flux tubes moving horizontally in the photosphere (braiding), and the second term corresponds to helicity injection due to emergence from the solar interior into the corona. Note, that adopting the naming convention from § 2.3.4, braiding and emergence terms would correspond to the FI and DI electric-field solutions respectively, and their sum to the DFI solution. When using the electric field form of the relative helicity flux (Eq. (41)), the braiding and emergence terms are not neatly separated.
4 Validation of the Electric Field Inversion Method Using Anmhd Simulations
4.1 ANMHD Run: Input, Output, Parameters
In this paper, we use a specific set of anelastic pseudo-spectral ANMHD simulations (Fan et al., 1999; Abbett et al., 2000, 2004) of an emerging magnetic bipole in a convecting box (Welsch et al., 2007) to test our improved electric field inversion technique. From ANMHD magnetic fields and plasma velocities, we know the actual electric fields, which we can compare to electric fields derived using the simulation’s evolving magnetic field. In the past, this ANMHD simulation has been used for several studies of velocity field inversions (Welsch et al., 2007; Schuck, 2008), as well as in test-cases for the first electric field inversions (Fisher et al., 2010, 2012). Fisher et al. (2012) showed that the PDFI-solution significantly improves the accuracy of the derived PI-solution, and the improvement from the knowledge of the Doppler velocity (PDI) is significantly more important than that of the horizontal velocity (PFI), at least in this example of magnetic flux emergence. In this paper, we perform a series of improvements to PDFI method beyond Fisher et al. (2012): we expand the derivation of the non-inductive contribution to non-normal viewing angles, introduce spherical coordinates for the PTD solution (Appendix A) and significantly speed up the Poisson equation solutions by using FISHPACK. With the above upgrades, we find that the PDFI method yields a good estimate of the electric field, Poynting and helicity fluxes and is ready to be applied routinely to the observed vector magnetograms.
We use a pair of ANMHD vector magnetograms, separated by
s with a pixel size of km, and a LOS
velocity map, observed at a specific viewing angle, to derive the
electric field and the vertical Poynting flux for a given set of
parameters (see Table 2). To estimate the horizontal
velocity field, , we use the Fourier local-correlation
tracking (FLCT) technique (Welsch et al., 2004; Fisher and Welsch, 2008) (
denote the Gaussian window size scale (a parameter in FLCT) by
. To suppress noisy behaviour, we only calculate
velocities where G, or of maximum
(Welsch et al., 2007).
In Table 2 we summarize all the input and output variables of the PDFI run and their typical ranges. Further, we vary the observed viewing angle within -range to estimate the accuracy of the PDFI method at non-zero viewing angles. We also vary the other free parameters: the number of iterations in the perpendicularization technique (, see § 2.2), the width of the Gaussian window used in the FLCT estimate for the horizontal velocity (, see below) and the width of polarity inversion line (, see § 2.3.1) to find the best values of the parameter set that yields the most accurate electric field solution (Right column).
For the remainder of this paper, to assess the performance of our methods, i.e. of the reconstruction of the ANMHD variable , we use the following metrics: (1) a fraction of the integrated total , (2) the slope or linear coefficient in the least-squares, polynomial fit, : , (3) the linear Pearson correlation coefficient , where is a standard deviation of , and (4) the normalized error of , , means that the error of the reconstruction is relative to the characteristic range of . Note that to exclude the weak-field background of , where the observed magnetic field (e.g. in HMI/SDO) tends to be noisy, we estimated and in locations where G. In the HMI vector magnetogram case, a similar threshold will be determined from estimated errors in the magnetic field values. To summarize, the ideal reconstruction of variable satisfies: , , and .
|FLCT velocity field||km/s|
|PDFI electric field||V/cm|
|Poynting flux||ergs/(s cm)|
|Parameters||Range (Best value)|
|No. of iterations in =0||iterations|
|Gaussian window width||pixels|
Figure 2 shows three panels that quantitavely justify selection of the best set of PDFI parameters (vertical dotted lines): number of iterations (, left panel), width of the Gaussian window (, middle panel), and the PIL width (, right panel). The left panel, , shows dependence of the RMS angle between the electric and magnetic field vectors on the number of iterations using the perpendicularization technique (see § 2.2). When looking for the optimal number of iterations , we try to keep as low as possible to achieve a high-speed performance, while still aiming for an angle close to . For the ANMHD case, without perpendiculatization (), the angle between the magnetic and electric field vectors . After only one iteration (), . The angle slowly increases to by . The convergence rate slows down, and reaches by =50. We chose as the optimal number of iterations, since above the convergence of toward is too slow to justify the additional computational effort. We also use when applying the non-normal viewing angle technique (§ 2.3.3) for finding the needed scalar potential .
The Middle Panel of Figure 2 shows how the quality of the horizontal velocity reconstruction, using the FLCT technique in the ANMHD test case, depends on the Gaussian window width . Since velocities parallel to the magnetic field do not affect the time evolution of the magnetic field and hence the fluxes of magnetic energy and helicity, in this plot we only compare components of the flow field that are perpendicular to . We find that pixels yield the worst agreement with the actual horizontal ANMHD plasma speed; pixels yields the best agreement, we therefore adopt it as the default value. For comparison, in Welsch et al. (2007) the optimal Gaussian window size was also chosen to be and in Fisher et al. (2012) .
Finally, the Right Panel of Figure 2, “Finding the best ”, shows the quality of the PDFI electric field components (red, blue, orange) and Poynting flux (black) reconstruction for different ’s, a free parameter that reflects the lack of confidence in the accuracy of the horizontal Doppler electric field away from PILs (see § 2.3). The panel shows that has a small effect on the quality of , and , above . We find that yields the best ratio between the total reconstructed and ANMHD Poynting fluxes: .
Using the best set of parameters (, and ) found above, in § 4.2, we estimate the quality of the PDFI reconstructed electric fields, helicity and Poynting fluxes and also estimate the uncertainties in the results at a zero viewing angle. In § 4.3 we describe how observing at the non-zero viewing angles affects the quality of the reconstruction. Finally in § 4.4 we estimate the roles that different non-inductive Doppler- and FLCT contributions play in the reconstruction and compare current results to Fisher et al. (2012).
4.2 Results: Electric Field, Poynting and Helicity Fluxes
Figures 3 and 4 show validation plots that compare electric field components and vertical Poynting fluxes derived from the PDFI-method with the actual ANMHD quantities. The PDFI-method reconstructs the ANMHD electric-field components quite well: the top and middle rows of Figures 3 look almost identical. The slope of the linear fit to the reconstructed versus the actual electric-field component ranges from to . The correlation coefficient, describing the quality of linear fit, in all cases is close to one. Using these electric-field components, we also find a good agreement for the vertical Poynting flux (Figure 4): the correlation coefficient , the slope and the fraction, . For comparison, Fisher et al. (2012) found slightly worse results: (Right Panel). It should be noted that the MEF method (Longcope, 2004) also accurately reconstructed the total Poynting flux in the tests by Welsch et al. (2007), but the spatial correlation between the true and reconstructed fluxes was significantly worse, at (see their Figure 14). Decomposing the total Poynting flux into potential and free components (see § 3.1), we find that the PDFI reconstructs of both the potential and free components () and the slope between the reconstructed and PDFI is one (). The free and potential components comprise and of the total unsigned Poynting flux respectively.
We also test how sensitive the Poynting fluxes are to errors in the vertical Doppler velocity. We find that if there is a random Doppler velocity noise of -km/s amplitude (), i.e. around to of the signal, then it does not significantly affect the Poynting flux: the error increases slightly from to , the slope remains close to one () and the fraction . However, if there is a bias Doppler velocity that increases all the velocities by km/s (towards the viewer) (), then the slope and the fraction increase to and the error is . Similarly, if all velocities decrease by km/s (), then the slope decreases to and the fraction to . This test demonstrates how important it is to remove the Doppler velocity bias (Welsch et al., 2013), when inferring electric fields from the observations.
In Figure 5 we compare actual ANMHD and the PDFI helicity flux rates calculated from . We find a good agreement between the two: the correlation coefficient , the slope , the fraction and the error . For comparison, using Fisher et al. (2012) electric fields we find a very similar helicity flux rate, but with a slightly larger scatter (): , , . The differences between our approach here and that of Fisher et al. (2012) are the adoption of the FISHPACK software to solve the two-dimensional Poisson equations, the ability to compute contributions to Doppler-shift electric fields from non-normal viewing angles, and a much more systematic and quantitative testing of the accuracy and robustness of the technique and its parameters.
4.3 Quality Of Electric Field and Poynting Flux Reconstructions At Non-Zero Viewing Angles
Figure 6 shows that at the error between reconstructed and the actual ANMHD electric fields is the smallest, and as the angle increases, the quality of the reconstruction decreases. For for all E-components, PDFI correctly identifies the slope between reconstructed and the actual ANMHD electric fields, within a difference. At the largest angle, , the slope is for and for , and the error in these variables reaches up to . In contrast to the horizontal electric field, the vertical component is relatively insensitive to the viewing angle: the slope varies within . The latter is not surprising, since the Doppler contribution, which has most of the angular dependence, primarily constrains the horizontal field. For all -components the correlation coefficient at is quite high, , implying that the estimates of the slope shown on the plot adequately describe the quality of the reconstruction.
Figure 7 describes the quality of the vertical Poynting flux reconstruction at different viewing angles. Two left panels show a point-to-point comparison between the ANMHD and PDFI Poynting fluxes at (left panel) and (middle panel). At the agreement between the ANMHD and PDFI is very good: the slope . Increasing the angle, at the scatter increases and the slope decreases to , the error and PDFI method recovers of the total actual flux. The right panel summarizes the quality of the Poynting flux reconstruction for viewing angles within -range. At close-to-normal angles, , PDFI recovers more than of the total energy flux and the error between reconstructed and the actual Poynting fluxes is less than 15% (). At larger angles, , of the total flux is recovered and the error . For the helicity flux rate, at , more than of the total flux is recovered correctly and the error .
4.4 Comparison Of Electric Field Reconstruction Techniques
In this section we analyze roles that the inductive and non-inductive components play in reconstructed vertical Poynting flux (Figure 8), electric field and helicity fluxes (Figures 9 and 10, Table 3).
Figure 8 shows scatter plots comparing the actual vertical Poynting flux with the the Poynting fluxes derived with different reconstruction methods: (1) P, (2) PI, (3) PFI, (4) PDI, (5) PDFI, (6) PDFI at non-normal angle , (7) FI, (8) DI, (9) DFI. The nomenclature that we use here is described in § 2.3.4 and Table 1. Using just the inductive part of the electric field (P), we reconstruct only of the total Poynting flux and the scatter from the linear dependence is large (). Adding the ideal MHD assumption (PI) adds more Poynting flux, leading to much less scatter (larger ) and a smaller error (). Inclusion of the non-ideal contribution due to horizontal plasma velocities, inferred from the FLCT (PFI), does not improve the solution. However, when we include the Doppler non-inductive component alone (PDI), we reconstruct of the flux, i.e. the role of the Doppler contribution is much higher than that from the FLCT horizontal velocities. Finally, when we add both the FLCT and the Doppler contributions (PDFI), the final Poynting flux is the closest to the actual , with a slightly smaller error () than in the PDI case (). On a separate note, if we use a non-PTD ideal inversion technique instead of the PTD, , (FI), we reconstruct only of the flux, i.e. the agreement between the reconstruction and the actual is poor (). If we know the vertical Doppler velocity (DI), then of the actual flux is reconstructed and with a good agreement (). With both vertical and horizontal velocities (DFI), we extract of the total Poynting flux, the slope and the error . How is this different from from the PDFI? The PDFI electric field that includes both inductive and non-inductive contributions yields a factor of two smaller error in () than the DFI field, and shows less scatter, especially in the regions of strong Poynting fluxes, thus it better represents the vertical Poynting flux. For completeness, in Figure 9 we also compare helicity flux rates, , from different reconstruction methods with the actual helicity flux rate. We remark that the crucial role the Doppler signal plays in reconstructing the Poynting and helicity fluxes in this case might be due to the process being modeled in the ANMHD simulation: an emerging magnetic bipole. It is possible that Doppler inputs to electric field estimates are less important when flux is not actively emerging.
In Figure 10 and Table 3 we summarize the quality of different electric field inversion techniques in the same way we did for and for all the variables, , where are the three components of the electric field, is the vertical Poynting flux, is the helicity flux rate (see Eq. (41-42)). Besides the methods shown previously (see Figure 8), we also add the method number 6, “PDFI (Fisher)”, that compares variables reconstructed in Fisher et al. (2012) with the actual ANMHD values.
The P-solution yields the worst reconstructions – the slope between the actual and derived and is