A MAGNETIC CALIBRATION OF PHOTOSPHERIC DOPPLER VELOCITIES

Magnetic energy passes from the solar interior into the Sun's outer atmosphere as an outward-directed Poynting flux S of magnetic energy across the photosphere,

$$\begin{equation} S_n = \hat{\mathbf {n}}\cdot {\mathbf {S}}= c\, \hat{\mathbf {n}}\cdot ({\mathbf {E}}\times {\mathbf {B}})/4 \pi , \end{equation} \tag{ 1 }$$

where c is the speed of light, E and B are the photospheric electric and magnetic fields, and $\hat{\mathbf {n}}$ is the outward-directed unit vector normal to the photosphere. We will refer to vectors perpendicular to the local normal as horizontal. We use the photosphere as the boundary between the solar interior and outer atmosphere primarily for convenience because the vector magnetic field is only routinely measured at the photosphere. A different atmospheric layer could be used if observations of the vector magnetic field were routinely produced there (see, e.g., Metcalf et al. 1995).

While B at the photosphere can be measured by vector magnetographs, E must be inferred by other means. A popular approach has been to use the observed magnetic evolution in sequences of magnetograms to estimate the electric field by inverting the finite-difference approximation to Faraday's law,

$$\begin{equation} \frac{\Delta {\mathbf {B}}}{\Delta t} = -c (\nabla \times {\mathbf {E}}) , \end{equation} \tag{ 2 }$$

where Δt is the difference between the times t_f and t_i of final and initial magnetograms, respectively, with E representing an average electric field over Δt. Note that for any E that satisfies this expression, (E − ∇ψ), where ψ is an arbitrary scalar potential function, will also satisfy it. Hence, Faraday's law does not fully constrain E.

We assume that the magnetic field is typically frozen to the plasma at the photosphere on length scales observable with the Helioseismic and Magnetic Imager (HMI), which has a pixel size Δx ≃ 05, or ≃ 360 km near disk center. Kubat & Karlicky (1986) estimated the magnetic diffusivity η_Ku from collisions at the photosphere to be 10⁸ cm² s⁻¹ in magnetized regions with B = 100 G, implying that the fluid could slip across the magnetic field with relative velocity u ∼ η_Ku/Δx ≃ 3 cm s⁻¹, which is completely negligible compared to typical photospheric speeds. Ignoring diffusivity, the ideal Ohm's law relates photospheric velocities, v, to the electric field,

$$\begin{equation} c {\mathbf {E}}= -{\mathbf {v}}\times {\mathbf {B}}. \end{equation} \tag{ 3 }$$

This implies that estimates of v can be used to determine the flux of magnetic energy (and magnetic helicity) across the photosphere (e.g., Démoulin & Berger 2003; Schuck 2006). Several techniques have been developed to estimate photospheric flows from ΔB_n/Δt, e.g., Chae (2001), Kusano et al. (2002), Welsch et al. (2004), Schuck (2006, 2008), and Fisher & Welsch (2008). These techniques are, however, imperfect (Rieutord et al. 2001; Welsch et al. 2007; Schuck 2008), so efforts to improve them are ongoing.

Fisher et al. (2010) presented a method to determine E from a sequence of vector magnetograms using a poloidal–toroidal decomposition (PTD) of the magnetic field with Faraday's law. Notably, the PTD method uses evolution of both the normal magnetic field, ΔB_n/Δt, and the normal electric current, ΔJ_n/Δt, to estimate E. While the PTD approach does not rely upon the ideal Ohm's law to estimate E, any Ohm's law with known resistive terms, including the ideal case, can be imposed post facto to constrain the solution for E.

For the ideal case, Fisher et al. (2012) recently presented a method to use Doppler measurements of the line-of-sight (LOS) component of the velocity, v_LOS, to better constrain E. They tested this approach using synthetic magnetograms and Doppler data from an MHD simulation. They noted that flows along B (commonly referred to as parallel or siphon flows) do not contribute to the ideal electric field in Equation (3), but can contribute to v_LOS in regions where the LOS component of the magnetic field is nonzero (B_LOS ≠ 0). Hence, they only incorporated Doppler data from areas (1) near polarity inversion lines (PILs), loci where the component of the magnetic field along their (assumed) LOS changes sign and the LOS field vanishes; and (2) where the field transverse to the LOS is large. Their tests demonstrated that including Doppler data near PILs substantially improves the estimation of both E and the normal Poynting flux S_n. For ideal evolution, this makes sense because, in principle, the Doppler electric field E^D from the LOS velocity v_LOS and transverse magnetic field B_trs along a PIL,

$$\begin{equation} c {\mathbf {E}}^D \equiv -(v_{\rm LOS} \times {\mathbf {B}}_{\rm trs}) \vert _{\rm PIL} , \end{equation} \tag{ 4 }$$

is not uncertain by the gradient of a scalar potential, as are estimates of E from Equation (2) alone.

In real magnetograms, procedures have been developed to automatically identify PILs (also sometimes called neutral lines; e.g., Falconer et al. 2003) of both the LOS components (e.g., Schrijver 2007; Welsch & Li 2008; Welsch et al. 2009) and normal components (Falconer et al. 2003) of magnetogram fields.

One factor hampering studies of Poynting fluxes has been the relative dearth of sequences of vector magnetograms. The launch of the SpectroPolarimeter (SP) instrument with the Solar Optical Telescope (SOT; Tsuneta et al. 2008) aboard the Hinode satellite (Kosugi et al. 2007) provided some seeing-free magnetogram sequences for investigations of PIL dynamics, such as observations suggestive of an emerging flux rope reported by Okamoto et al. (2008). However, the cadence of the SP instrument is relatively slow compared to timescales of photospheric evolution on arcsecond scales, and the field of view (FOV) of SOT is limited. Vector magnetograms of active region (AR) fields at higher cadence and over larger FOVs should be routinely produced by the HMI instrument (Scherrer et al. 2012) aboard the Solar Dynamics Observatory (SDO) and the SOLIS vector spectromagnetograph (VSM; Keller et al. 2003), enabling routine estimates of Poynting flux.

Given the association between solar activity and magnetic fields, a quantitative description of processes responsible for the introduction and removal of magnetic flux into the solar atmosphere is essential to understand solar activity.

New magnetic flux typically appears at the solar photosphere via emergence of Ω-shaped loops from the convection zone (Figure 1, top row, left panel). About 3000 ARs are typically cataloged by NOAA (which requires each to contain a visible sunspot umbra) over one 11-year cycle of solar activity, each with ∼10²² Mx (10¹⁴ Wb) in unsigned flux. Estimates vary for the time required for emergence to replace the small-scale magnetic fields in the quiet-Sun network, with unsigned flux on the order of ∼10²³ Mx over the photospheric surface, but they are generally on the order of a day or less (e.g., Hagenaar et al. 2003). We note that new flux can also appear in photospheric magnetograms via the submergence of the bases of U-shaped loops from above, as seen in simulations by Abbett (2007; Figure 1, top row, right panel), though this process is not often discussed. Hence, the term appearance encompasses both emergence and submergence. This definition of the term differs substantially from that used by Lamb et al. (2008), in which appearance refers to initial identification of a feature by a tracking algorithm, where the definition of "feature" is algorithm-dependent.

Zoom In Zoom Out Reset image size

Much like the appearance of new flux, the removal of magnetic flux from the photosphere also plays a central role in solar activity. Clearly, for the large-scale solar dynamo to operate cyclically, flux that is introduced to the photosphere must eventually be removed or recycled. Similarly, flux emergence in the small-scale dynamo must be statistically balanced by flux removal (Schrijver et al. 1997).

How is flux removed from the photosphere? The short answer is "flux cancellation," which Livi et al. (1985) defined in observational terms as "the mutual apparent loss of magnetic flux in closely spaced features of opposite polarity" in magnetogram sequences. Physically, cancellation could correspond to (1) the emergence of U-shaped magnetic loops (e.g., Lites et al. 1995; van Driel-Gesztelyi et al. 2000), (2) the submergence of Ω-shaped loops (e.g., Rabin et al. 1984; Harvey et al. 1999; Chae et al. 2004; Iida et al. 2010), or (3) reconnection in the magnetogram layer (e.g., Yurchyshyn & Wang 2001; Kubo & Shimizu 2007; see also Welsch 2006). These possibilities are sketched in the bottom row of Figure 1.

Spruit et al. (1987) proposed that opposite polarities within ARs could be connected by U-loop subsurface extensions, and that flux in these extensions might emerge as weak, "sea-serpent" (undulating) fields between strong-field regions, eventually canceling with AR flux via U-loop emergence. Low (2001) proposed a process that also removes most AR flux from the photosphere by the emergence of U-loops. He suggests that the flux rope that forms an AR contains many turns, such that a magnetogram can intersect the flux rope many times after the flux rope has partially emerged, and as the flux rope continues to emerge, much of the flux cancels by U-loop emergence, with successive magnetograms showing less and less flux remaining until little is left. In contrast to these models, van Ballegooijen & Mackay (2007) and van Ballegooijen (2008) constructed numerical models of AR flux systems remaining anchored near the base of the convection zone, with most canceling flux retracting to "repair" the toroidal magnetic field deep in the solar interior.

Spruit et al. (1987) noted that sea-serpent cancellation might occur on small scales. On unobservably small scales, this would produce apparent in situ flux disappearance (Wallenhorst & Howard 1982; Wallenhorst & Topka 1982; Gaizauskas et al. 1983). Further, unlike the "self-cancellation" of AR flux modeled by Spruit et al. (1987), Low (2001), and van Ballegooijen & Mackay (2007), AR fields might also cancel with the ubiquitous, small-scale fields of the quiet Sun (e.g., Lin & Rimmele 1999; Harvey et al. 2007), including unresolved fields (Sánchez Almeida 2009). This process would formally remove equal amounts of AR and quiet-Sun flux from the photosphere, but if the latter were below a given magnetograph's noise threshold, it would be undetectable, and the AR flux would have seemed to disappear. While case studies of cancellation have been made, it is unclear how much AR flux cancels on observable scales over a solar cycle. Hence, it is possible that cancellation with quiet-Sun flux is the dominant method of AR flux removal. Moving magnetic features in moat flows around sunspots can explain rates of sunspot flux loss (e.g., Kubo et al. 2008), but much flux in ARs lies outside sunspots, and its disappearance must still be understood.

Because magnetic fields are divergence-free, all field lines (tangent lines of the vector field B) form closed loops. (These might be infinitely long, or ergodic, but in any case field lines do not end.) Magnetic flux can therefore only emerge or cancel where magnetic fields are tangent to the photosphere. These loci correspond to normal-field PILs, where the normal magnetic field vanishes and regions of positive and negative flux, where emerging (or submerging) field lines thread the photosphere, are nearby. In our terminology, the appearance (or cancellation) of magnetic flux increases (resp., decreases) the total unsigned magnetic flux at the photosphere.

If the photospheric electric field during cancellation is ideal or nearly so, then measurements of time-averaged Doppler shifts along PILs should be able to distinguish between cancellation via either U-loop emergence or Ω-loop submergence. The lack of a clear Doppler signal while LOS flux cancels would be consistent with reconnective cancellation. We note that the magnetic field along PILs away from disk center in LOS magnetograms can have a component that is normal to the photosphere, implying that LOS PILs away from disk center do not exactly correspond to sites of flux appearance (or cancellation). Hence, only Doppler shifts along LOS PILs near disk center can effectively constrain the physical processes at work in cancellation.

Several case studies of Doppler shifts at cancellation sites have been undertaken. Yurchyshyn & Wang (2001) and Bellot Rubio & Beck (2005) reported Doppler shifts consistent with upflows at the cancellation sites they studied, which they interpreted as outflows from reconnective cancellation. Chae et al. (2004) and Iida et al. (2010) found evidence for flux submergence during cancellation. Kubo & Shimizu (2007) studied cancellations along several PILs in ASP magnetograms, and reported mixed Doppler signals, which they interpreted in terms of reconnective cancellation at multiple heights. To determine the rest wavelengths used to compute Doppler velocities, Chae et al. (2004) used quiet-Sun values of line center, while Kubo & Shimizu (2007) and Iida et al. (2010) estimated their rest wavelengths by averaging line centers over their FOVs. As described below, however, these approaches to determining rest wavelengths probably yield biased estimates.

One factor complicating the use of Doppler data both to estimate photospheric electric fields and to constrain dynamics in flux cancellation is inaccuracies in the determination of the rest wavelengths of photospheric lines in ARs. Uncertainties in rest wavelength can arise from both instrumental effects (e.g., from thermal variations in components) and biases in analysis techniques. In addition, the gravitational redshift is present (Takeda & Ueno 2012). Periodicities in magnetic fields estimated by HMI (Liu et al. 2012a) on orbital timescales (12 and 24 hr) suggest that instrumental effects are present in estimated magnetic fields. It is plausible that similar effects should be present in Doppler signals, a point we revisit below.

Because the position of line center is typically computed from sampling predominantly quiet-Sun regions, where line profiles are systematically shifted blueward by the convective blueshift (Dravins et al. 1981; Cavallini et al. 1986; Hathaway 1992; Asplund & Collet 2003; Schuck 2010), Doppler shifts in ARs typically exhibit "pseudo-redshifts." The blueward bias of quiet-Sun line profiles occurs because rising plasma is both (1) brighter than sinking plasma (since rising plasma is hotter) and (2) occupies a greater fraction of an instrument's FOV than sinking plasma (since upwelling convective cells are larger than downflow lanes). Consequently, a determination of line center position based upon the statistical properties (e.g., means or medians) of Doppler images (Dopplergrams) is biased by the upward motion of quiet-Sun plasma. Observations and modeling, by, e.g., Gray (2009) and Asplund & Collet (2003), respectively, suggest that the magnitude of this bias can range from a few hundred m s⁻¹ to nearly 1 km s⁻¹ for various photospheric lines. Because AR magnetic fields inhibit convection (e.g., Welsch et al. 2012), measured line centers in ARs are redshifted relative to any rest wavelength derived from quiet-Sun Doppler measurements. (Helioseismology requires accurate measurement of changes in Doppler shifts, so is insensitive to errors in determination of the rest wavelength.) P. Scherrer (2009, private communication) has stated that one must account for the convective blueshift to accurately determine Doppler shifts in ARs.

By fitting measured Doppler shifts over the disk with profiles that account for differential rotation, meridional flow, and center-to-limb variations, any overall constant Doppler shift (sometimes referred to as the convective limb shift; Hathaway 1992) can be estimated (Snodgrass 1984; Hathaway 1992; Schuck 2010). This approach, however, has a major physical uncertainty: the physics of the center-to-limb variation in average Doppler shift in the particular spectral line used by HMI (or any other spectral line) involves detailed interactions between height of formation, the height of convective turnover, the variation with viewing angle of the average convective flow speed, and the variation with viewing angle of optical depth (Carlsson et al. 2004; Takeda & Ueno 2012). We note that diverging flows tangent to the photosphere in granules can, depending upon optical depth in granules at the formation height of the line, produce a convective blueshift toward the limb, because diverging flows on the near sides of granules approach the viewer, while receding flows on the far sides of the granules are at least partially obscured by the optical depth of the granules. This is sketched in Figure 2. Hence, fitting the observed center-to-limb variation in line-center positions does not imply that all bias from convective motions has been determined. In principle, simulations of convection near the photosphere (e.g., Asplund & Collet 2003; Fleck et al. 2011) could be used to characterize the expected center-to-limb variation in line centers, which could then be used to remove the modeled convective blueshift bias (but not any instrumental biases). We are unaware of published comparisons of observed and modeled center-to-limb variations in the HMI spectral line.

Zoom In Zoom Out Reset image size

The goal of this paper is to describe a "magnetic calibration" technique to estimate any spatially uniform bias in measured Doppler velocities. In basic terms, this method estimates such a bias by requiring statistical consistency between two independent measures of changes in flux near LOS-field PILs, the loci where flux emerges and submerges: (1) ΔΦ_LOS, half the change of total unsigned LOS flux near PILs over a time interval Δt; and (2) ΔΦ_vBL, the amount of transverse flux transported upward or downward across the magnetogram layer over Δt, computed by summing the Doppler electric field—the product of the Doppler velocity and transverse magnetic field strength—along PILs. This consistency relies upon both Faraday's law and the assumption of ideal evolution.

The rest of the paper is organized as follows. In the next section, we describe our calibration method in greater detail. We then apply the method to the initial vector magnetogram sequence released by the HMI Team, starting with a description of our preliminary treatment of the data (Section 3) followed by a step-by-step demonstration of the method (Section 4) and analysis of its results (Section 5). We conclude with a brief summary in Section 6.

To reduce the full data set to a more manageable size, we focus on a subset of the full five-day sequence. To study photospheric magnetic evolution prior to the X flare, and to baseline this evolution against post-flare evolution, we retain about 72 hr of data, from midnight at the start of the 13th until midnight at the end of the 15th, inclusive of endpoints. Given the 12-minute cadence, the resulting time series consists of 361 time steps. The AR was at S19E11 at 00:30UT on 2011 February 13 and S21W27 at 00:30UT on 2011 February 16. Pixels lacking data near the edges of the cutout FOV, from artifacts of the cutout process, were cropped: columns [0–24] and rows [0–5] were removed for all steps.

Next, Doppler velocities were corrected for spacecraft motion. Due to SDO's geosynchronous orbit, its velocity along the radial, Sun-observer line can be large, of the order of ±5 km s⁻¹. Also, because the Earth is orbiting westward about the Sun, and SDO orbits Earth, there is a significant projection of SDO's westward motion onto LOSs to many pixels. In addition, there is a nonzero component of SDO's northward velocity onto the LOS. While the radial component is larger than the W and N components, the latter can be significant—a few tens of m s⁻¹ or more. Accordingly, the spacecraft velocity was projected onto the LOS to each pixel—in the small-angle approximation—and subtracted off the measured Doppler velocities. Then Stonyhurst latitudes and longitudes of each pixel (accounting for the solar B and P angles) were computed to correct Doppler velocities for solar rotation, using the "magnetic/2 day lag" rotation rate found by Snodgrass (1983), for which

$$\begin{eqnarray} &&\qquad\quad \omega _{\rm rot} = A + B \sin (\tilde{\theta })^2 + C \sin (\tilde{\theta })^4 \ \ \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} A &=& 2.902\qquad\ \ \ \ \ \ \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} B &=& -0.464\qquad\ \, \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} C &=& -0.328 ,\qquad \end{eqnarray} \tag{ 28 }$$

where ω_rot is the photospheric rotation rate, $\tilde{\theta }$ is latitude, and A, B, and C are in microrad s⁻¹. Accurate correction for rotation prior to application of our calibration procedure is not necessary; but without removing it, the bias velocity estimates will include (and could be dominated by) the Doppler shift from rotation.

We also note that the Doppler signal from helioseismic p-mode oscillations, which have a peak in power at frequencies near 1/300 Hz, could introduce non-convective signals into the 720 s cadence Dopplergrams. While individual mode amplitudes are on the order of a few tens of cm s⁻¹, the overall motion due to the superposed modes can be 500 m s⁻¹ or more, which is comparable to convective motions. In practice, however, HMI's 720 s Dopplergrams are created by averaging higher-cadence filtergrams with a cosine-apodized, 1215 s boxcar (Liu et al. 2012b) with a FWHM of 720 s (J. Schou 2012, private communication). This effectively acts as a filter, which would reduce a random-phase sinusoidal signal with period of 300 s and amplitude of 500 m s⁻¹ to an averaged value of less than 20 m s⁻¹ (less than the noise level of HMI Dopplergrams).

Accurate estimation of electric fields by PTD (or, really, any method) requires that changes in B between magnetograms arise from physical processes on the Sun, as opposed to artifacts of the measurement process. Hence, any spurious changes in the measured field, where they can be identified, should be removed.

We noticed that in successive images of the transverse-field azimuths, the inferred azimuths in some regions flip by nearly 180° from one frame to the next, as would be expected from errors in resolving the 180 ambiguity. The 7 successive frames (12 minutes apart) in the left column of Figure 5 show an example of this effect, with time increasing downward. The area shown corresponds to that in the solid white box in Figure 4, and the initial frame corresponds to that in Figure 4. Azimuths range over [0, 360] deg. The next column shows maps of the angular difference (the interior angle) between the current and previous azimuths. Differences in azimuths range from [−180, 180] deg (note the rotated color wheel in the top panel). The black patches in rows 4–6 of this column correspond to a region in which azimuths flip by nearly 180° in a region from one frame to the next. The checkerboard pattern in the same region in rows 2–3 of this column implies a spatially alternating pattern of 180° azimuth flips, and probably arises from a lack of convergence in the simulated annealing algorithm used to infer the azimuths (Leka et al. 2009). Also, the red and blue speckled regions toward the left side of each frame in this column exhibit rapid spatial variations in azimuth changes. These speckled regions correspond to areas with relatively weak-field strength (B < 200 G), where azimuth determinations could be more problematic than in strong-field regions.

Zoom In Zoom Out Reset image size

Interpreting frame-to-frame azimuth flips over finite regions as spurious, we seek to identify and remove these artifacts. Our goal is to automate the detection of such flips, and our basic approach is to identify suspicious changes in azimuth, which we envision as "top hats" (or inverted top hats) in the running differences of azimuths in individual pixels: large, positive (or negative) jumps in azimuth for one time step, followed by reversals to the pre-jump level at the next step. The actual patching is simpler: we add 180° to the azimuth, and output the flipped data modulo 360°.

We tried more than one approach to detecting and correcting spurious azimuth changes before finding one that we think works well. To save other researchers from repeating our efforts, we now describe in detail both our failed approach and the approach we found worked best.

In our initial unsuccessful approach, we simply identified all pixels with unsigned changes in azimuth from frame (i − 1) to frame i above a given threshold (we tried, e.g., 120°, 135°, 150°, and 170°), and flipping these in frame i. The updated frame i was then used as a reference for finding large changes in frame (i + 1). Using this approach, once a pixel is flipped, its new state tends to be propagated forward in time. But legitimate magnetic evolution can cause this state to be incorrect, leading to spatially isolated pixels that differ from their neighbors. This contradicts a central tenet of the original ambiguity resolution procedure: currents should be minimized. Further, the number and spatial arrangement of problematic pixels appears to depend strongly on threshold.

The best approach we found is to flip azimuths only in pixels where both (1) the jump in azimuth was larger than 120° between frames; and (2) the jump in azimuth increased the time-averaged, unsigned changes in azimuth over a the set of frames within ±N steps of a given frame. The use of a range of images on either side of the ith image is necessary to capture the essential character of a top hat: a large change, followed by a reversal. (A large azimuth change from one frame to the next, by itself, might be legitimate—e.g., by horizontal advection of a region with a spatial transverse field reversal.) We determined suspicious flips by comparing means of two arrays of 2N azimuth differences: (1) unsigned differences between the azimuths at frame i and azimuths in ±N neighboring frames in time, [i − N: i − 1] and [i + 1: i + N]; and (2) unsigned differences between the flipped azimuths (hypothetical data, with all pixels' azimuths flipped) at frame i with actual azimuths in ±N neighboring frames in time, [i − N: i − 1] and [i + 1: i + N]. Pixels for which flipping would decrease the mean, unsigned, frame-to-frame angular differences are then flipped.

We call frames in which azimuths have been flipped filtered. Filtered azimuths are used for past frames, [i − N: i − 1], while unfiltered azimuths are used for future frames [i + 1: i + N], so this approach is not time symmetric. This approach should be able to deal with top hats that are N steps wide, but not wider. We tried both N = 4 and N = 2, and show results with N = 2 here. Using N = 4 results in 2%–5% fewer flips in a given frame, so it is slightly more restrictive. We did not, however, see much difference between results from N = 2 and N = 4 in strong-field regions, so most differences are probably in weak-field regions where inference of the direction of transverse fields is less reliable anyway.

The third column from left in Figure 5 shows filtered data from the same seven successive frames as in the leftmost column, with suspicious changes in azimuths identified by our approach flipped by 180°. The rightmost column shows maps of the angular difference between the current and previous azimuths (note the rotated color wheel in the top panel) from the third column. The checkerboard and black patches visible in the second column are not present in this column. Also, changes in azimuths in the red and blue speckled regions on the left sides of each frame are not as large in this column as in the second column. This suggests our approach decreases fluctuations of azimuths in weak-field regions, too.

The first step in our approach is to use an automated algorithm to identify PILs in LOS magnetograms. Our procedure creates masks of each polarity for pixels with |B_LOS| above a field strength threshold of B_thr, dilates each mask by one pixel, and finds all areas where the dilated masks overlap (see Schrijver 2007; Welsch & Li 2008). The threshold field is the only input parameter in PIL identification procedure; here, we use B_thr = 50 Mx cm⁻², ensuring that we identify relatively strong fields in close proximity. This approach identifies structures 1–3 pixels wide, which we erode into single-pixel-width lines with IDL's morph_thin.pro procedure. We define the resulting pixels to be PIL pixels. In Figure 6, we show identified PIL pixels, color coded by (possibly biased) Doppler velocity v'_LOS, in a close-up view of the LOS magnetic field at the time step when AR 11158 was closest to disk center.

Zoom In Zoom Out Reset image size

In this subregion of the full AR, N_PIL = 98 PILs were identified. Most of the Doppler signals along PILs appear redshifted. These are almost certainly pseudo-redshifts, arising from the Doppler zero point being defined to be blueward of the true rest wavelength for photospheric plasma with no LOS velocity. The velocity zero point for HMI is derived from the median of Doppler velocities within 90% of the solar radius over 24 hours (Y. Liu & S. Couvidat 2011, private communication). This implies that contributions from quiet-Sun (and therefore strongly convecting) plasma dominate the median, and probably bias the zero point.

Having identified PILs at a given time t₀, our second step is to estimate the set of changes of LOS flux {ΔΦ_LOS} near all PILs over the Δt = 24 minutes between the two LOS magnetograms observed ±12 minutes from the magnetogram at t₀. To do so, for each PIL we create a binary map of that PIL's pixels, then apply IDL's dilate.pro procedure with a (d × d) matrix of 1's to the binary map to define that PIL's "neighborhood mask." The dilation parameter d is the only other free parameter in our method. While granular flows have speeds on the order of ∼1 km s⁻¹, such flows are short lived, with lifetimes of ∼5 minutes. Over Δt = 24 min, therefore, displacements of fluxes tend to be consistent with lower averaged speeds, on the order of ∼0.5 km s⁻¹ or less (see, e.g., Figure 13 of Welsch et al. 2012). Given the HMI pixel size Δx ≃ 05 ≃ 360 km near disk center, magnetic flux should not move by more than about two pixels over Δt = 24 minutes. Here, we use d = 7—that is, dilation by ±3 pixels in all directions—although we have explored other choices (see below). (As Démoulin & Berger 2003 have pointed out, emerging fields that are strongly inclined with respect to the LOS could have very high apparent speeds, implying displacements that would exceed our 3 pixel dilation. Excessive dilation, however, risks incorporating changes in flux that are unrelated to emergence/submergence in PIL neighborhoods.) Next, we multiply the two co-aligned LOS magnetograms from ±12 minutes by this neighborhood mask and sum the unsigned flux in each product, then difference these summed fluxes. Because this counts LOS flux changes in both polarities, we divide the result by two to estimate ΔΦ_LOS over Δt = 24 minutes.

It should be noted that converging or diverging horizontal flows unrelated to emergence/submergence at each PIL can transport "background" flux into the neighborhood of a PIL, and near-PIL flux out of this neighborhood, in some cases dominating our estimates of the change in LOS flux due to emergence or submergence along the PIL. Hence, the computed rate of change in LOS flux along an individual PIL is subject to large errors. Identifying collections of like-polarity pixels as "features" and tracking their evolution (DeForest et al. 2007; Welsch et al. 2011) might be one way to better account for fluctuations in LOS flux due to flux concentration/dispersal (Lamb et al. 2008). Assuming the contributions from converging and diverging flows cancel in the aggregate over many PILs, the statistical properties of changes in LOS fluxes along a set of PILs are more robust to such errors, and therefore more useful for our purposes.

Our third step is to compute, for each PIL, the corresponding expected flux change from the (assumed biased) Doppler signal, ΔΦ'_vBL, and magnetic length, $\overline{BL}$. Since the expressions derived for these quantities (Equations (18) and (23)) in Section 2 refer to the component of the transverse magnetic field perpendicular to the PIL, we investigated the angle of the transverse magnetic field along PILs in our data set. Figure 7 illustrates the geometry of the transverse field along several PILs. The upper left panel shows unit vectors $\hat{\mathbf {x}}_{\perp }$ perpendicular to the PIL, determined from the gradient of B_LOS evaluated at PIL pixels. The upper right panel shows unit vectors $\hat{\mathbf {b}}_{\rm trs}$ along B_trs, evaluated at PIL pixels. Note the predominance of unit vectors $\hat{\mathbf {b}}_{\rm trs}$ pointing toward the bottom of the FOV: this coherence might be spurious, and could arise if the minimization of currents by the ambiguity resolution algorithm introduces artificial, large-scale correlations in azimuths. Color coding in the bottom left panel shows the cosine of the angle χ between these unit vectors—their dot product, $\hat{\mathbf {x}}_{\perp } \cdot \,\hat{\mathbf {b}}_{\rm trs}$—with red positive, and the saturation level set to ±0.866. (The sign of the angle χ is defined such that the angle between $\hat{\mathbf {x}}_\perp$ and $\hat{\mathbf {x}}_\parallel$ is +90°; see Equation (10).) Fields along PILs evidently exhibit components in both the "normal" (positive-to-negative; blue) and "inverse" (negative-to-positive; red) directions (Martens & Zwaan 2001) perpendicular to the PIL. For each PIL, the mean (PIL-averaged) signed and unsigned cos  (χ) and unsigned sin  (χ) can be computed; the bottom right panel shows the distributions of these averages for all PILs in our data set. The peaks at −1 and 1 for the signed cos  (χ) and around 0.5 for |sin  (χ)| imply that the average field of many PILs: (1) is spatially coherent; and (2) tends to point primarily across (not along) the PIL.

Zoom In Zoom Out Reset image size

As with changes in LOS flux near PILs—which can arise from converging or diverging flows unrelated to emergence/submergence—some component of Doppler velocities along identified PILs might be unrelated to emergence/submergence. Slight inaccuracies in PIL identification, for instance, could lead us to improperly include flows along the magnetic field (siphon flows) in our flux transport rates. Also, filling factors—due to scattered light within the HMI instrument, or imaged emission from unresolved, unmagnetized plasma, or both—could bias the estimated flux transport rate in a given pixel (and contributions from either would plausibly tend to be blueshifted). Hence, as with changes in LOS flux near PILs, the inferred rate of transport of transverse flux along an individual PIL is subject to large errors. Again, the statistical properties of flux transport rates along the aggregated collection PILs should be more robust to error, and therefore more useful for quantifying the pseudo-redshift.

We have also investigated making the simplifying replacement B_⊥ → B_trs when computing ΔΦ_vBL and $\overline{BL}$. Physically, one can justify this approximation by noting that the footpoints of fields with components tangent to the PIL must still thread the photosphere somewhere, if not directly across the PIL, as would be the case for no magnetic component tangent to the PIL. We will show results computed both ways, using either B_⊥ or B_trs at each PIL pixel, but unless otherwise stated, results shown below are derived using B_trs. As will be seen, this did not drastically change the estimates of bias velocities v₀, although quantities derived using B_⊥ were substantially noisier. We also note that using B_trs makes the bias estimation technique independent of ambiguity resolution.

Finally, we can estimate any bias velocity present. We first difference ΔΦ_LOS − ΔΦ'_vBL on each PIL to compute the set of bias fluxes {ΔΦ_bias} from all PILs. We must then estimate the coefficient v₀ from the ratios of {ΔΦ_bias} to $\lbrace \Delta t \overline{BL} \rbrace$.

In the data, errors are present in both measured quantities. Before estimates of typical uncertainties were published by the HMI team, we adopted uniform uncertainties of 20 m s⁻¹ for the Doppler velocities, 25 Mx cm⁻² for B_LOS, and 90 Mx cm⁻² for B_trs. This value for the Doppler uncertainty is consistent with the expected near-disk-center noise levels reported by Schou et al. (2012). The value for noise in B_LOS is significantly larger than the HMI Team's since-published estimate the uncertainty of ∼6 Mx cm⁻² in 720 s data (e.g., Liu et al. 2012a), while the value for B_trs is probably closer to the HMI Team's estimate of noise in the transverse component on the order of 10² Mx cm⁻² (Sun et al. 2012).

We estimate v₀ and its uncertainty by analyzing the set of ratios $\lbrace \Delta \Phi _{\rm bias}/(\Delta t \overline{BL}) \rbrace$ over a subset of identified PILs. Based upon both our prior knowledge that the convective blueshift biases the median-derived estimate of the rest wavelength in HMI Doppler data blueward, and the predominance of redshifted PILs in Figure 6, we expect that the biased PIL velocities are more red than the true velocities. Recalling that we use the astrophysical convention that redshifts correspond to positive velocities with respect to the observer, the bias velocity v₀ should be positive.

Not all bias fluxes, however, are consistent with a positive bias velocity. This is not obvious from Equation (24), since it deals with absolute values. Consequently, we now consider the different possibilities for the relative sizes of ΔΦ_LOS, ΔΦ'_vBL, and ΔΦ_bias. A positive v₀ is consistent with PILs that obey either

$$\begin{equation} \Delta \Phi _{vBL}^{\prime } > |\Delta \Phi _{\rm LOS}| \end{equation} \tag{ 29 }$$

or

$$\begin{equation} -|\Delta \Phi _{\rm LOS}| < \Delta \Phi _{vBL}^{\prime } < 0 , \end{equation} \tag{ 30 }$$

since the correction to ΔΦ'_vBL is $-\Delta \Phi _{\rm bias}^{\prime } = -v_0 \Delta t \overline{BL}$. For PILs that obey either

$$\begin{equation} |\Delta \Phi _{\rm LOS}| > \Delta \Phi _{vBL}^{\prime } > 0 \end{equation} \tag{ 31 }$$

or

$$\begin{equation} \Delta \Phi _{vBL}^{\prime } < -|\Delta \Phi _{\rm LOS}| < 0 , \end{equation} \tag{ 32 }$$

however, v₀ would have to be negative to improve agreement between the flux changes.

What fraction of PILs are consistent with a positive bias velocity v₀? Of the ∼50,000 PILs identified in the full data set (with B_thr = 50 Mx cm⁻²), 73.8% and 5.7% were consistent with Equations (29) and (30), respectively, and 13.9% and 6.6% obeyed Equations (31) and (32), respectively. Consequently, ∼80% of PILs are consistent with a positive bias velocity.

In Figure 8, we show the distribution of bias fluxes for all PILs at all time steps. (In any single time step, there are ∼140 PILs in our FOV, too few to form a continuous distribution.) The distribution is strongly skewed toward positive bias fluxes. We consider PILs matching either Equation (31) or Equation (32) to be pathological, which we attribute to "noise," primarily systematic errors in our estimates of the fluxes. (Errors from our assumed uncertainties in the underlying magnetic and Doppler data would be 3 × 10¹⁶ Mx per pixel or less.)

Zoom In Zoom Out Reset image size

This distribution of bias fluxes can be used to quantify uncertainties in the bias flux. Mirroring the negative bias flux values across zero and fitting the result with a Gaussian enables an empirical estimate of noise in the bias flux from the fitted width, 3.75 × 10¹⁷ Mx. Following the same procedure but using B_⊥ instead of B_trs to compute ΔΦ'_vBL results in a fitted width of 5 × 10¹⁷ Mx, reflecting greater uncertainty in estimates of this flux. These define thresholds in the bias flux: estimates of bias fluxes smaller than these levels fall within our (systematic) uncertainties.

Subtracting the fit from the total distribution of bias fluxes leaves the residual distribution of bias fluxes above our uncertainty level. At each step, we computed separate estimates of the bias velocity from both the full set of PILs and just the subset of PILs with bias fluxes more positive than these estimates of the threshold bias fluxes. Not surprisingly, estimates of v₀ using only PILs with bias fluxes above these thresholds—which we refer to as the "high-bias" subset of PILs—were substantially higher (as shown below) than estimates derived using all PILs. We view use of all PILs as more conservative, so results shown below are derived from all PILs, unless explicitly stated otherwise.

To investigate which of our estimates for v₀ is most reasonable, we now consider the physical implications of our results in terms of the upward and downward transport of flux across the photosphere. Recall that processes that increase photospheric flux are Ω-loop emergence and U-loop submergence, while processes that remove flux are Ω-loop submergence and emergence. Without correcting Doppler velocities, essentially all increases and decreases in flux are attributed to the two submergence processes (from U-loops and Ω-loops, respectively), since nearly all PILs show only redshifts. Adjusting the observed Doppler shifts to compensate for the convective blueshift, however, leads to a different apportionment between the four possible processes, depending upon the applied bias velocity: PILs that, on average, show upward transport of transverse flux and increases (decreases) in LOS flux are ascribed to Ω-loop (U-loop) emergence, while PILs that, on average, show downward transport of transverse flux and decreases (increases) in LOS flux are ascribed to Ω-loop (U-loop) submergence.

In Figure 14, we show the frequencies of these processes, as functions of the number of pixels along all identified PILs at all time steps, for three possible corrections: the average estimates of v₀, derived using B_trs with B_thr = 50 Mx cm⁻² and d = 5 (top); the median of Doppler velocities on pixels in identified PILs (middle); and the average high-bias estimates of v₀, again using B_trs with B_thr = 50 Mx cm⁻² and d = 5 (bottom). In all panels, the black curve shows the distribution of sizes of all PILs. Processes involving blueshifts (Ω-loop and U-loop emergence) are plotted in blue, while processes involving redshifts (Ω-loop and U-loop submergence) are plotted in red. Processes that remove flux via cancellation (U-loop emergence and Ω-loop submergence) are plotted with dashed lines. Even with the average-v₀ correction (top panel), submergence processes still dominate, contradicting the expectation that the observed increase in unsigned radial flux of ∼2 × 10²² Mx over the data sequence occurred via emergence from the interior. Removing the median Doppler velocity over all pixels in identified PILs at each time step gives a more reasonable distribution of processes (middle panel): emergence and submergence processes are more closely balanced. But the requirement of a net flux increase implies that emerging processes should predominate, and this expectation is only consistent with the high-bias-v₀ correction (bottom panel).

Zoom In Zoom Out Reset image size

In principle, observations of changes in total unsigned radial flux in the FOV should be consistent, separately, with estimates of changes in LOS flux near PILs and the upward/downward transport of transverse flux. Our estimates of flux increase by Ω-loop emergence from changes in LOS flux near PILs are roughly consistent with the net change in LOS flux over the entire sequence. Unfortunately, other components of our flux budgeting near PILs are inaccurate, in two ways.

First, our estimates of flux removed by cancellation processes (Ω-loop submergence and U-loop emergence) are too large (from ∼10²² Mx to several times this), whether inferred from LOS changes or rates of transport of transverse flux. This might be explained if some flux is dispersing near PILs that we identify as sites of cancellation, instead of submerging/emerging at these PILs. As discussed in Section 4.2, identifying magnetogram features and tracking their evolution (DeForest et al. 2007), or imposing additional constraints (e.g., requiring cancellation to persist for multiple time steps, or taking the flux lost to be the minimum of that among the two polarities; see Welsch et al. 2011), could more accurately quantify canceled flux. We leave investigation of these approaches for future work.

Second, our estimated rates of transport of transverse flux are systematically too large compared to changes in LOS flux near PILs, for all processes. As noted in Section 4.2 above, inaccuracies in PIL identification could incorporate velocities from siphon flows in our flux transport rates. In addition, a filling factor f < 1 could also play a role, since it would affect |ΔΦ_LOS| and |ΔΦ_vBL| differently: while length transverse to the LOS (two factors of which enter computation of |ΔΦ_LOS|) is scaled by $\sqrt{f}$, length along the LOS (which enters computation of |ΔΦ_vBL|) is not. Hence, if f < 1, one would expect |ΔΦ_vBL| > |ΔΦ_LOS|. Given typical discrepancies in fluxes by factors of 4–6, this would imply filling factors of ∼0.05. This would be implausible in umbrae, where the field is expected to be space filling. We analyze PILs, however, which are, almost as a matter of definition, in mixed-polarity regions, where the field is expected to be more intermittent. The plausibility of a filling factor this small along PILs is unclear.

Another possibility is that flux transport rates are overestimated because the transport is not ideal, i.e., that the plasma is moving at a different speed than the flux transport rate. For instance, if the flux emergence process were not ideal (Pariat et al. 2004), then the rate of flux transport derived from Doppler shifts would underestimate the true emergence rate, since flux would be partially decoupled from the plasma, allowing flux to move upward independently of the plasma. While magnetic diffusivity would be manifested in that case by flux slipping through stationary (or much slower-moving) plasma at a rate faster than that suggested by ideal transport, the ideal transport rate is too high in our estimates. Chae et al. (2008) note that turbulent diffusivity can act on canceling LOS magnetic fluxes; but again, we find the rate of loss of LOS field to be too small relative to ideal transport of horizontal fields, and invoking turbulent diffusivity would make this discrepancy worse. One possible explanation is that unmagnetized plasma, which could be systematically faster-moving than magnetized plasma, could be contributing significantly to Doppler signal.

To quantitatively characterize cases of apparently non-ideal evolution, one could employ a simple non-ideal model, assuming that a uniform magnetic diffusivity η is present. This would imply a corresponding non-ideal electric field parallel to the PIL (in addition to the ideal electric field along the PIL already used to estimate ΔΦ_vBL), given by

$$\begin{equation} c E_\parallel ^{\rm NI} = \eta |\partial _{\rm LOS} B_{\perp } - \partial _\perp B_{\rm LOS}| , \end{equation} \tag{ 33 }$$

where, as above, $\hat{\mathbf {x}}_\perp$ locally points across the PIL (see Equation (10)). For PILs where |ΔΦ_LOS| ≫ |ΔΦ_vBL|, we could attribute the excess ΔΦ_LOS to η ∂_⊥B_LOS integrated along the PIL, taking ∂_LOSB_⊥ to be negligible, as it would be if there were a current sheet along the PIL (e.g., Litvinenko 1999). Values of η derived this way could be compared to the assumed values used by Litvinenko (1999), Chae et al. (2002), and Litvinenko et al. (2007). For PILs where |ΔΦ_LOS| ≪ |ΔΦ_vBL|, we could attribute the excess ΔΦ_vBL to overestimation of the transport rate of transverse flux along the LOS, due to slippage of transverse flux through the moving plasma in the presence of diffusivity. Because ∂_LOSB_⊥ is unknown for such cases, however, we could not estimate η.

We note that spurious changes in LOS fields over a magnetogram FOV, due, e.g., to artifacts arising from orbital motion in the case of HMI (or any other cause) will introduce errors into the bias velocity estimate: our approach attributes changes in unsigned LOS flux to emergence/submergence, so changes unrelated to actual emergence/submergence will introduce error. Are the known artifacts in HMI field strengths significant? Probably not, since the orbital variations occur on timescales much slower than the evolution over Δt = 24 minutes that we analyze. Liu et al. (2012a) report changes in mean |B_LOS| of approx. 150 Mx cm⁻² peak to peak, i.e., over 12 hr (see their Figure 5). This rate of change corresponds to a change of 5 Mx cm⁻² over 24 minutes, which is at the noise level of an individual LOS measurement. This implies these slow variations in LOS field strength would be dominated by fluctuations from noise, or faster evolution.

After filtering azimuths and removing Doppler velocities due to solar rotation and convection, we interpolated the data to sub-pixel accuracy for two purposes: (1) image co-registration, to remove frame-to-frame jitter, primarily from quantized changes in the cutout window boundaries (determined by pixel size) imposed on the smooth rotation of the underlying solar features; and (2) Mercator projection of the spherical data onto a uniform grid in a Cartesian frame. While interpolation propagates errors from pathological pixels, and can introduce other artifacts (e.g., smoothing) into the data, it should nonetheless be undertaken, since significant whole-frame shifts were found (some frame-to-frame shifts were on the order of 1 pixel, and cumulative shifts showed both secular evolution and jumps) in the x and y directions by Fourier cross-correlation in the B_LOS image sequence.

Since azimuth filtering and interpolation modify the input data, we characterized uncertainties in the data set empirically instead of relying upon error arrays produced in the HMI processing pipeline. To do so, we constructed histograms of the LOS and transverse fields, B_LOS and B_trs, respectively, at each time step, and statistically characterized these distributions.⁹

In each distribution of B_LOS, we interpret the weak-field core of distribution as a measure of noise (e.g., DeForest et al. 2007), which we quantified by fitting the core of the distribution (±6 Mx cm⁻²) with a Gaussian of width σ_L, and by finding the flux density at half max (HM) in the left and right wings. Statistically, the mean and standard deviation of fitted σ_L's were 4.1 ± 0.3 Mx cm⁻² and combined left width and right width at HM were 4.7 ± 0.5 Mx cm⁻², but fluctuations with a periodicity around 6 hr were observed. In SDO's geosynchronous orbit, (1) Doppler motion toward and away from the Sun and (2) thermal cycling due to Earth's dayside albedo would both exhibit 12 hr periodicities, and would be a quarter cycle out of phase, so could combine to produce variations in instrumental response with such a characteristic 6 hr periodicity. This analysis roughly agrees with the HMI Team's estimate of the uncertainty of 6 Mx cm⁻², and suggests that the uncertainty in the LOS field component assumed above, 20 Mx cm⁻², probably overestimates the true uncertainty.

Since B_trs is a magnitude, it is necessarily positive, and we find that the distribution peaked away from zero. If the distributions of each of the three observed field components were Gaussian, with the same width σ, then the distribution of transverse field strengths would be that of a 3D Maxwell distribution, ∝B²_trsexp (− B_trs²/2σ²). In fact, we find that a Gaussian displaced from zero (fitted to data within ±25 Mx cm⁻² of the tallest peak) consistently fit each distribution of B_trs more closely (as quantified by χ² over the whole distribution) than a 3D (or 2D) Maxwellian. It should be noted, however, that the shape of the distribution varied substantially over the observing sequence—e.g., exhibiting two peaks at some time steps—so the Gaussian fits, although better than the Maxwellians, were still poor in many cases. The means and standard deviations of the Gaussians' fitted centroids B_t0 and widths σ_t were 81 ± 17 Mx cm⁻² and 29 ± 6 Mx cm⁻², respectively. Both also exhibited periodicities, though with more complicated temporal structure than with fits to the distribution of B_LOS, with the clearest period at ∼24 hr. This analysis suggests that our assumed uncertainty in the transverse field component, 90 Mx cm⁻², is reasonable.