Sparse reconstruction of electric fields from radial magnetic data

TL;DR

This paper introduces a sparse reconstruction method for estimating localized horizontal electric fields on the Sun's surface from radial magnetic data, improving modeling of solar corona energy transfer.

Contribution

A novel sparse basis pursuit algorithm is developed to produce localized electric field estimates from magnetic maps, addressing the non-uniqueness of the inverse problem.

Findings

Performs well on flux-balanced magnetic maps in Cartesian and spherical geometries.

Localization fails with flux-imbalanced maps common in data assimilation.

Main obstacle is flux imbalance in input data, affecting realistic electric field estimation.

Abstract

Accurate estimates of the horizontal electric field on the Sun's visible surface are important not only for estimating the Poynting flux of magnetic energy into the corona but also for driving time-dependent magnetohydrodynamic models of the corona. In this paper, a method is developed for estimating the horizontal electric field from a sequence of radial-component magnetic field maps. This problem of inverting Faraday's law has no unique solution. Unfortunately, the simplest solution (a divergence-free electric field) is not realistically localized in regions of non-zero magnetic field, as would be expected from Ohm's law. Our new method generates instead a localized solution, using a basis pursuit algorithm to find a sparse solution for the electric field. The method is shown to perform well on test cases where the input magnetic maps are flux balanced, in both Cartesian and spherical geometries. However, we show that if the input maps have a significant imbalance of flux - usually arising from data assimilation - then it is not possible to find a localized, realistic, electric field solution. This is the main obstacle to driving coronal models from time sequences of solar surface magnetic maps.

Equations

$${\partial }_t\mathbf{B}=-\nabla \times \mathbf{E},$$

$${\partial }_t{B}_r=-{\mathbf{e}}_r\cdot \nabla \times {\mathbf{E}}_{\perp }.$$

$${\mathbf{E}}_{\perp}=-\nabla\times\big{(}\Phi{\mathbf{e}}_{r}\big{)}-\nabla_{\perp}\Psi,$$

$$-{\nabla }^2\Phi ={\partial }_t{B}_r,$$

$$\Phi (x,y,t)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\partial }_t{B}_z(x^{\prime },y^{\prime },t)$$

$$\displaystyle\qquad\times\log\big{|}(x-x^{\prime})^{2}+(y-y^{\prime})^{2}\big{|}^{1/2}\,\mathrm{d}x^{\prime}\,\mathrm{d}y^{\prime}.$$

$${\partial }_t{B}_z=\exp \left(-\frac{(x+\rho {)}^2+y^2}{{\delta }^2}\right)-\exp \left(-\frac{(x-\rho {)}^2+y^2}{{\delta }^2}\right).$$

$$\Phi (x,y)=\frac{{\delta }^2}{2}\log \left(\frac{r_{+}}{r_{-}}\right)+\frac{{\delta }^2}{4}{E}_1\left(\frac{r_{+}^2}{{\delta }^2}\right)-\frac{{\delta }^2}{4}{E}_1\left(\frac{r_{-}^2}{{\delta }^2}\right),$$

$$\Vert {\mathbf{E}}_{\perp }{\Vert }_2^2=\Vert {\mathbf{E}}_{\perp 0}{\Vert }_2^2+\Vert {\nabla }_{\perp }\Psi {\Vert }_2^2+2{\oint }_{\partial D}\Psi {\mathbf{E}}_{\perp 0}\cdot \mathbf{n}\mathrm{d}S.$$

$$\Delta x\Delta y\frac{\partial B_z^{i,j}}{\partial t}=\Delta x{E}_x^{i,j+1/2}-\Delta x{E}_x^{i,j-1/2}$$

$$+\Delta y{E}_y^{i-1/2,j}-\Delta y{E}_y^{i+1/2,j}$$

$$E_{x0}^{i,j+1/2}$$

$$E_{y0}^{i+1/2,j}$$

$$\displaystyle\frac{\partial B_{z}^{i,j}}{\partial t}=\frac{1}{\Delta y^{2}}\Big{(}2\Phi^{i,j}-\Phi^{i,j+1}-\Phi^{i,j-1}\Big{)}$$

$$\displaystyle+\frac{1}{\Delta x^{2}}\Big{(}2\Phi^{i,j}-\Phi^{i-1,j}-\Phi^{i+1,j}\Big{)},$$

$$\Vert \mathbf{x}{\Vert }_1:=\sum _i|x_i|\text{subject to}A\mathbf{x}=\mathbf{b}.$$

$${\partial }_t{B}_z=x\exp \left(-\frac{x^2+y^2}{{\delta }^2}\right).$$

$${\mathbf{E}}_{\perp }=\frac{{\delta }^2}{2}\exp \left(-\frac{x^2+y^2}{{\delta }^2}\right).$$

$${\beta }^{i_0,j_0}(i-i_0)\exp \left(-\frac{(i-i_0{)}^2+(j-j_0{)}^2}{{\delta }^{i_0,j_0}}\right)$$

$$B_z(x,y,t)=\frac{a^2}{a^2+4\eta t}\exp \left(-\frac{x^2+y^2}{a^2+4\eta t}\right).$$

$$E_x=-\frac{2a^{2}y}{(a^2+4\eta t{)}^2}\exp \left(-\frac{x^2+y^2}{a^2+4\eta t}\right),$$

$$E_y=\frac{2a^{2}x}{(a^2+4\eta t{)}^2}\exp \left(-\frac{x^2+y^2}{a^2+4\eta t}\right).$$

$${\partial }_t{B}_z=4a^2\eta \left[\frac{x^2+y^2}{(a^2+4\eta t{)}^3}-\frac{1}{(a^2+4\eta t{)}^2}\right]$$

$$\times \exp \left(-\frac{x^2+y^2}{a^2+4\eta t}\right).$$

Full text

magnetohydrodynamics, Sun: activity, Sun: evolution, Sun: magnetic fields, Sun: photosphere

\AuthorCallLimit

=1 \fullcollaborationNameThe Friends of AASTeX Collaboration

1 Introduction

Magneto-hydrodynamic (MHD) simulations of the magnetic field in the Sun’s corona are typically driven by an imposed evolution on the solar surface. Since the magnetic field ${\mathbf{B}}$ evolves according to Faraday’s law

[TABLE]

the required boundary condition is the horizontal electric field ${\mathbf{E}}_{\perp}=E_{\theta}(\theta,\phi,t){\mathbf{e}}_{\theta}+E_{\phi}(\theta,\phi,t){\mathbf{e}}_{\phi}$, written here in spherical polar coordinates. Unfortunately, this electric field cannot be observed directly, but must be reconstructed from other observations. If we assume that the plasma obeys the ideal Ohm’s law ${\mathbf{E}}=-{\mathbf{v}}\times{\mathbf{B}}$ then ${\mathbf{E}}_{\perp}$ can, in principle, be computed from vector observations of ${\mathbf{B}}$ and of the plasma velocity ${\mathbf{v}}$.

In practice, however, such observations are not routinely available for the full solar surface. Many authors have therefore opted to estimate ${\mathbf{E}}_{\perp}$ purely by inverting Equation (1), typically using only a time sequence of $B_{r}(\theta,\phi,t)$ data (Mikić et al., 1999; Wu et al., 2006; Mackay et al., 2011; Cheung & DeRosa, 2012; Yang et al., 2012). In this case, only the radial component of (1) is used,

[TABLE]

The present paper seeks to address one of the potential pitfalls of this technique.

The fundamental problem of electric field inversion from Faraday’s law (1) is lack of uniqueness. If we write ${\mathbf{E}}_{\perp}$ as a Helmholtz decomposition

[TABLE]

then $\nabla\times{\mathbf{E}}_{\perp}$ depends only on the potential $\Phi(\theta,\phi,t)$, and is independent of the choice of the second potential $\Psi(\theta,\phi,t)$. Note that our notation interchanges $\Phi$ and $\Psi$ compared to Mikić et al. (1999). Given that $\partial_{t}B_{r}$ is independent of $\Psi$, the simplest way to obtain a solution consistent with the given $B_{r}$ data is to set $\Psi\equiv 0$. Then (2) becomes the Poisson equation

[TABLE]

which may be solved uniquely for $\Phi$ given appropriate boundary conditions. We call the resulting ${\mathbf{E}}_{\perp}$ the inductive solution for the electric field, and denote it by ${\mathbf{E}}_{\perp 0}$. This solution has been used to drive coronal models in numerous studies (e.g., Mikić et al., 1999; Amari et al., 2003; Mackay et al., 2011; Cheung & DeRosa, 2012; Gibb et al., 2014; Yang et al., 2012; Feng et al., 2015). Unfortunately, there is good reason to believe that the non-inductive potential $\Psi$ is non-negligible on the real solar surface.

One problem with the inductive solution is its lack of localization. If the real ${\mathbf{E}}_{\perp}$ satisfies Ohm’s law, then it ought to vanish outside of patches of strong $B_{r}$. As discussed in Section 2, this is not the case with the inductive solution; the main focus of this paper will be to choose $\Psi$ in order to correct this deficiency.

A second problem with the inductive solution is that it cannot detect electric fields associated with plasma flows along contours of $B_{r}$, since $\partial_{t}B_{r}=0$ in that case. The corresponding contribution to $\Psi$ may be added if the flows are known. For example, Kazachenko et al. (2014, 2015) have made this correction using observations of plasma velocities in active regions, and the resulting ${\mathbf{E}}_{\perp}$ has been used to drive magneto-frictional simulations (Fisher et al., 2015). On the other hand, Weinzierl et al. (2016) have included the contribution to ${\mathbf{E}}_{\perp}$ from a prescribed large-scale differential rotation in their global magneto-frictional simulations. In this paper, however, we will put this issue aside and focus only on the localization problem.

In Section 2, we describe the localization problem with the inductive solution, before characterizing the inductive solution in a variational framework. This makes the link to our proposed sparse electric field solution, described in Section 3, that is designed to restore localization. Numerical tests in both Cartesian and spherical geometry are presented in Section 4, while Section 5 considers the limitations of the sparse solution.

2 Inductive Electric Field

2.1 The localization problem

Since the inductive potential $\Phi$ solves the Poisson equation (4), the solution may be expressed in terms of the Green’s function as

[TABLE]

Here we have taken Cartesian coordinates and an infinite plane for simplicity, but similar solutions hold on a spherical surface. For a localized source $\partial_{t}B_{z}$, we therefore have that $\Phi\sim\log(r)$ for large distance $r$ from the source. The corresponding inductive electric field ${\mathbf{E}}_{\perp 0}$ therefore decays only as $r^{-1}$. In particular, it may be non-zero well outside the regions of non-zero $B_{z}$.

As an explicit example, consider a bipolar distribution

[TABLE]

In this case, one may solve Equation (4) exactly to obtain the closed-form solution

[TABLE]

where $r_{\pm}=\sqrt{(x\pm\rho)^{2}+y^{2}}$ and $E_{1}(x):=\int_{x}^{\infty}\mathrm{e}^{-t}/t\,\mathrm{d}t$ is the exponential integral. The logarithmic behaviour is clear from this expression, and the corresponding electric field components are plotted in Figure 1.

2.2 Variational formulation

Another way to characterise the inductive solution ${\mathbf{E}}_{\perp 0}$ is by the fact that it minimizes the $L_{2}$-norm $\|{\mathbf{E}}_{\perp}\|_{2}:=\left(\int_{S}|{\mathbf{E}}_{\perp}|^{2}\,\mathrm{d}S\right)^{1/2}$ among all possible solutions to (2). This is straightforward to see by inserting the expression (3) into the norm, which leads to

[TABLE]

On a periodic or infinite domain, the boundary term vanishes, and it follows that $\|{\mathbf{E}}_{\perp}\|_{2}$ is minimized by choosing $\Psi\equiv 0$. In this sense, the inductive electric field ${\mathbf{E}}_{\perp 0}$ is the smoothest possible ${\mathbf{E}}_{\perp}$ satisfying Faraday’s law for a given $\partial_{t}B_{z}$.

2.3 Discrete formulation

In practice, we work with ${\mathbf{E}}_{\perp}$ defined on a discrete numerical grid, so it is useful to formulate the analogous discrete problem. In this paper, we discretise ${\mathbf{E}}_{\perp}$ on a staggered grid (Yee, 1966), where $E_{x}$ is defined on the horizontal cell edges, $E_{y}$ on the vertical cell edges, and $\partial_{t}B_{z}$ at the cell centers (Figure 2). Such grids are commonly used in MHD simulations. We discretise Faraday’s law (2) using Stokes’ theorem, so that the equation

[TABLE]

must hold for each grid cell $i,j=1,\ldots,n$. (In more general curvilinear coordinates, the edge lengths and cell areas would be different for each cell.) This constitutes a system of linear equations for the unknowns $E_{x}$ and $E_{y}$ on each edge. However, the system is highly under-determined, since there are $2n(n+1)$ unknowns but only $n^{2}$ equations. This non-uniqueness of solution reflects our freedom in choosing the non-inductive component.

In the discrete case, the inductive solution may be computed by writing

[TABLE]

where $\Phi^{i,j}$ is defined at the cell centres. Then (9) becomes

[TABLE]

which is simply the standard 5-point stencil for the Poisson equation. In the examples below, we solve this using a standard fast-Poisson solver (Press et al., 1992).

However, it is instructive to think of the system of equations (9) in the more abstract form $A{\mathbf{x}}={\mathbf{b}}$, where ${\mathbf{x}}$ is the vector of unknowns ($E_{x}$, $E_{y}$), ${\mathbf{b}}$ is the vector of $\partial_{t}B_{z}$ values, and $A$ is the $n^{2}\times 2n(n+1)$ matrix corresponding to equations (9). The inductive solution is then the solution to $A{\mathbf{x}}={\mathbf{b}}$ that minimizes the discrete $\ell_{2}$-norm $\|{\mathbf{x}}\|^{2}_{2}:=\sum_{i}x_{i}^{2}$. In other words, it is the (unique) least-squares solution to this under-determined system of equations. As such, the solution may be written in terms of the Moore-Penrose pseudo-inverse as ${\mathbf{x}}=A^{\top}(AA^{\top})^{-1}{\mathbf{b}}$, although in practice it is much more efficient to use the fast-Poisson solver. Nevertheless, viewing the inductive solution as the $\ell_{2}$-minimum makes the connection with the sparse solution that we will describe in Section 3. The sparse solution is found by minimizing a different discrete norm of ${\mathbf{x}}$.

3 Sparse Electric Field

Our idea is to find a sparse solution for ${\mathbf{x}}$ (i.e., $E_{x}$, $E_{y}$) that minimizes the number of non-zero values. This should be more localized than the inductive solution ${\mathbf{E}}_{\perp 0}$. Since it will differ from ${\mathbf{E}}_{\perp 0}$, this new solution will have a non-zero potential $\Psi$. However, we will work directly with the system $A{\mathbf{x}}={\mathbf{b}}$ (Equation 9), rather than solving for $\Psi$ (or $\Phi$) explicitly.

Instead of minimizing the $\ell_{2}$-norm of ${\mathbf{x}}$, as in the inductive solution, we propose to minimize the $\ell_{1}$-norm. In other words, minimize

[TABLE]

In the optimization literature, problem (13) is known as basis pursuit (Chen et al., 2001). It is used in a wide variety of fields, and has recently been applied to the determination of differential emission measures in the solar corona (Cheung et al., 2015). Although minimizing the $\ell_{1}$-norm is not strictly equivalent to minimizing the number of non-zero components of ${\mathbf{x}}$, the minimum-$\ell_{1}$ solution to an under-determined system of linear equations is often the sparsest solution to that system (Candes et al., 2006; Donoho & Tsaig, 2008). Numerically, it is preferable as it is a convex optimisation problem that can be efficiently solved using linear programming (Gill et al., 2011). For the examples in Section 4, we used the numerical implementation of basis pursuit in the SparseLab111SparseLab is freely available for download from http://sparselab.stanford.edu. library. This implementation of basis pursuit uses a primal dual method (Chen et al., 2001) and has a for the optimization.

For all of the test cases in Sections 4 and 5, with tolerance $10^{-12}$, the optimization took only 10 to 20 seconds on a desktop workstation, not including the time to construct the matrix $A$ (which is only needed once for each grid). This is certainly fast enough for practical application to a time sequence of $\partial_{t}B_{r}$ maps, as would be needed for driving a coronal MHD simulation.

4 Numerical Tests

In this section, we compare the inductive and sparse ${\mathbf{E}}_{\perp}$ solutions for two test cases. In the first case (Section 4.1), the true ${\mathbf{E}}_{\perp}$ is known through Ohm’s law, while in the second (Section 4.2) it is not. Moreover, the first case uses Cartesian geometry and the second case uses spherical geometry. The additional problems that arise when applying the technique to real magnetic data are discussed in Section 5.

4.1 Bipolar distribution

Our first test case is based on the simplest configuration that may arise from Ohm’s law. In fact, due to the linearity of both Ohm’s law and Equation (2), more general test configurations may readily be obtained by superimposing multiple copies of this basic configuration. For simplicity, the computation is done in coordinates, in a square domain $|x|<3$, $|y|<3$, using the discretization described in Section 2.3. Periodic boundary conditions are applied in both directions.

Our basic solution is the bipolar distribution

[TABLE]

Physically, this could arise from uniform advection of a single magnetic polarity $B_{z}=\exp[-(x^{2}+y^{2})/\delta^{2}]$ under an ideal Ohm’s law ${\mathbf{E}}_{\perp}=-{\mathbf{v}}_{\perp}\times(B_{z}{\mathbf{e}}_{z})$ with velocity ${\mathbf{v}}_{\perp}=\delta^{2}/2{\mathbf{e}}_{x}$, in which case

[TABLE]

(In this case, we are taking a snapshot at $t=0$, as the polarity is centred on the origin.) Alternatively, the same $\partial_{t}B_{z}$ (with different $B_{z}$) might correspond to the emergence of a bipolar magnetic field from the solar interior. The target ${\mathbf{E}}_{\perp}$ given by (15), along with the $\partial_{t}B_{z}$ in (14), are illustrated in the top row of Figure 3.

The middle row of Figure 3 shows the inductive solution ${\mathbf{E}}_{\perp 0}$ for this $\partial_{t}B_{z}$, which is qualitatively similar to the illustrative solution in Section 2. Not only is ${\mathbf{E}}_{\perp 0}$ not localized within the region of non-zero $\partial_{t}B_{z}$, but the topology is wrong: there is a substantial $E_{x0}$ component, despite the fact that $E_{x}=0$ for the target solution. The sparse solution using basis pursuit, shown in the bottom row of Figure 3, is substantially more accurate. To demonstrate convergence, Figure 4 shows the absolute errors for the sparse solution, as a function of both grid resolution and of the tolerance parameter in the basis pursuit algorithm. Firstly, the error in $\partial_{t}B_{z}$ is independent of tolerance or grid resolution, indicating that the constraint $A{\mathbf{x}}={\mathbf{b}}$ is preserved to high accuracy in all cases. Secondly, there is convergence in both $E_{x}$ and $E_{y}$ as the tolerance is reduced. The error in $E_{x}$ is independent of resolution, whereas that in $E_{y}$ saturates at a (higher) level depending on grid resolution. This simply reflects the truncation error in the numerical approximation (9) for the curl, which is quadratic in $\Delta x$. The saturation is not seen in $E_{x}$ because of our particular target solution $E_{x}=0$.

4.2 Spherical distribution

To demonstrate how the sparse reconstruction performs on a more realistic example, we have taken a snapshot from a flux-transport simulation of $B_{r}(\theta,\phi,t)$ in spherical coordinates, covering the full solar surface. The reason for using a simulated map rather than an observed synoptic magnetogram is to ensure perfect flux balance; the consequences of not having flux balance will become evident in Section 5.

We used the flux-transport model described by Yeates et al. (2015), for which the numerical code is freely available (https://github.com/antyeates1983/sft_data). As input data, we used synoptic maps of $B_{r}$ from the Global Oscillation Network Group (GONG, gong.nso.edu/data/magmap/), starting in Carrington rotation 2073 and generating the output snapshot in rotation 2109. As described by Yeates et al. (2015), the input maps are used (a) to initialize $B_{r}$ at the beginning of the simulation, and (b) to extract strong flux regions used to update $B_{r}$ during the evolution. The computation used a $360\times 180$ grid, equally spaced in longitude ($\phi$) and sine latitude ($\cos\theta$).

For context, Figure 5(a) shows the $B_{r}$ map from the flux-transport simulation, while Figure 5(b) shows the $\partial_{t}B_{r}$ map from the same time. Since the simulation uses a supergranular diffusion rather than imposing convective velocities, the map appears smoother than would a real observed map. To make the $\partial_{t}B_{r}$ map look more realistic, we have added random noise to produce the map in Figure 5(c). The noise is carefully constructed to preserve local flux balance, and is generated by adding a bipolar distribution

[TABLE]

centred at each pixel $(i_{0},j_{0})$ on the $(s,\phi)$ grid. The magnitude $\beta^{i_{0},j_{0}}$ at each pixel is chosen from a normal distribution with mean zero and $\sigma=0.005\max_{i_{0},j_{0}}|\partial_{t}B_{r}|$, and the size $\delta^{i_{0},j_{0}}$ is either 1 or 2 pixels, with equal probability. In addition, we rotate the pattern by 90 degrees at each particular pixel with equal probability. This generates a more realistic map as shown in Figure 5(c).

Figure 6 shows the inductive and sparse reconstructions of ${\mathbf{E}}_{\perp}$ from the map in Figure 5(c). To account for the spherical coordinates, Equations (10)–(12) have been modified to include the necessary coordinate factors. The right-hand column of Figure 6 shows that both methods preserve $\partial_{t}B_{r}$ to high accuracy, as in the Cartesian test. Unlike in the Cartesian test, however, we no longer compare to a target ${\mathbf{E}}_{\perp}$, since the ${\mathbf{E}}_{\perp}$ corresponding to our added noise is unknown. Nevertheless, that the sparse reconstruction is the inductive reconstruction. the $\partial_{t}B_{r}$ distribution is more complex than in the previous test, the sparse ${\mathbf{E}}_{\perp}$ is still much localized within regions of strong $\partial_{t}B_{r}$ . Moreover, it has a weaker $E_{\phi}$ than $E_{\theta}$, consistent with the fact that the dominant Ohm’s law contribution in the flux-transport simulation was from differential rotation ($v_{\phi}$). This is not the case in the inductive solution. Thus we are confident that the sparse solution not just for simple test cases, but also for more complex $\partial_{t}B_{r}$ distributions that are qualitatively similar to those on the Sun.

5 Limitations

As with any inverse problem, care is needed in the application of our sparse reconstruction technique. In this section, we consider two limitations which can be important in the practical application to the driving of coronal MHD models.

5.1 Diffuse distributions of $\partial_{t}B_{r}$

The sparse solution makes the physical assumption that ${\mathbf{E}}_{\perp}$ should be localized because it satisfies Ohm’s law for a localized ${\mathbf{B}}$. If the magnetic field distribution is too diffuse, then we can not expect the assumption of localization to recover the correct (diffuse) ${\mathbf{E}}_{\perp}$.

To illustrate what happens if we nevertheless try to apply the sparse reconstruction in such a situation, consider the horizontal diffusive spreading of an initial gaussian magnetic polarity

[TABLE]

From the resistive Ohm’s law ${\mathbf{E}}_{\perp}=\eta\nabla\times(B_{z}{\mathbf{e}}_{z})$, we get

[TABLE]

Faraday’s law then gives

[TABLE]

For this exercise, we fix $a=0.5$ and $\eta t=0.1$. This solution is shown in the top row of Figure 7.

It is important to note that, in this case, Ohm’s law is compatible with a purely inductive solution, because ${\mathbf{E}}_{\perp}$ from Ohm’s law has the inductive form with $\Phi(x,y,t)=\eta B_{z}(x,y,t)$. Indeed, the middle row of Figure 7 shows the inductive solution to reproduce the target ${\mathbf{E}}_{\perp}$ in this case. However, the bottom row shows that the sparse solution using basis pursuit does not converge to the target. Rather, it favors a concentration of $E_{x}$ along vertical lines, and $E_{y}$ along horizontal lines. This behavior is typical when the target ${\mathbf{E}}_{\perp}$ is too diffuse. It is not possible for the total spatial extent of ${\mathbf{E}}_{\perp}$ to be more localized than $\partial_{t}B_{r}$, owing to the constraint of satisfying (2). But it is possible for more of the electric field to be concentrated in thinner regions, and this is what minimizing the $\ell_{1}$-norm does.

This limitation is an obvious one, and limits the sparse solution technique to input maps where ${\mathbf{B}}$ is sufficiently localized. Fortunately, this situation is typical on the Sun, given high enough resolution. The flux-transport test case in Section 4.2 showed the method to work for realistic solar magnetic maps, at the resolutions of present-day simulations. However, there is another more subtle problem with real data, that we describe in the following section.

5.2 Flux imbalance in $\partial_{t}B_{r}$ maps

For a localized ${\mathbf{E}}_{\perp}$ to exist, it is necessary not only that $\partial_{t}B_{r}$ is localized, but that the net flux $\int_{S}\partial_{t}B_{r}\,\mathrm{d}A$ vanishes over the region $S$ of localization. To see this, apply Stokes’ theorem on some closed curve that encircles $S$. If the net flux in $S$ is non-zero, then there must be non-zero ${\mathbf{E}}_{\perp}$ somewhere on the bounding curve, and indeed on any curve enclosing non-zero net flux.

This condition of vanishing net flux is satisfied by both of our simple examples (14) and (20), and by every local flux concentration in the flux-transport test case (Section 4.2). But it need not be satisfied in an observationally-derived map of $\partial_{t}B_{r}$, even if a flux correction has been applied. We will demonstrate this problem firstly in a controlled test and secondly in a “real” data-assimilative model.

For the controlled test, we modify the input map from Section 4.2 to simulate the typical situation where new magnetogram observations are assimilated only within a finite region on the visible face of the Sun. The most severe problems of flux imbalance occur when a new active region has emerged on the far-side of the Sun, and is gradually assimilated into the magnetic map as it rotates into view.

Figure 8 shows our test. Here the original $\partial_{t}B_{r}$ map, from Figure 5(c), has been modified to simulate the gradual assimilation of a new bipolar active region. In the right-hand column of Figure 8, the assimilation region is indicated on the modified $\partial_{t}B_{r}$ maps. Each row corresponds to a different time as the new active region rotates into view, with the times differing by one day of solar rotation ($27.2753^{\circ}$ in these Carrington maps).

When the new region is fully contained within the assimilation window (bottom row of Figure 8), it makes only a localized contribution to the sparse ${\mathbf{E}}_{\perp}$, correctly reproducing ${\mathbf{E}}_{\perp}$ as in Section 4.1. But when the region is only partially observed, there is an unbalanced flux of $\partial_{t}B_{r}$, both in the assimilation window and in the entire map. Correspondingly, it is impossible to find a localized electric field, and indeed the sparse solution fails. In the area of the new region, we see similar linear structures to the diffusive example in Figure 7, where the flux was also not locally balanced. But notice that the electric field is modified not only near to the new region, but also at farther distances.

Of course, this is a rather unrealistic test. If this ${\mathbf{E}}_{\perp}$ was being used to drive a coronal MHD model, for example, one would correct the $\partial_{t}B_{r}$ map for flux balance, before trying to compute ${\mathbf{E}}_{\perp}$. In Figure 9 we show the effect on $E_{\theta}$ of first making a flux-balance correction to the $\partial_{t}B_{r}$ map (from the top row of Figure 8). The middle and right plots show the results with two different methods of flux correction. In the “additive” method, an equal amount is subtracted from each pixel in the $360\times 180$ grid to remove the imbalance. In the “multiplicative” method, all of the pixels where $\partial_{t}B_{r}>0$ are multiplied by a constant factor $f_{+}$, and all pixels where $\partial_{t}B_{r}<0$ are multiplied by another constant factor $f_{-}$. The factors $f_{+}$, $f_{-}$ are chosen so that the net flux vanishes after scaling, and is equal to the mean of the positive and negative fluxes before scaling. Comparing the results in Figure 8, it is evident that neither method of flux correction gives much improvement in the reconstructed $E_{\theta}$ (or $E_{\phi}$). Only by limiting the flux correction to the active region itself could the solution be improved.

Similar behavior is found when we apply the sparse reconstruction technique to a time sequence of $B_{r}$ maps generated by the Air Force Data-Assimilative Photospheric Flux Transport (ADAPT) model (Arge et al., 2010; Henney et al., 2012). This model assimilates observed magnetograms from the visible face of the Sun into a surface flux-transport simulation, so as to approximate the global distribution of $B_{r}$ on the solar surface, as a function of time. The data assimilation means that the electric field ${\mathbf{E}}_{\perp}$ is not known everywhere in the map, so must be reconstructed, if these maps are used to drive coronal MHD simulations.

As an illustration, we choose an ADAPT run driven by GONG magnetograms, and extract $\partial_{t}B_{r}$ for 2014 November 16, 00:00 UT. This particular date has been deliberately chosen since Weinzierl et al. (2016) found a significant non-localized ${\mathbf{E}}_{\perp 0}$ caused by re-assimilation of a large active region (see their Figure 11). Here, we process the ADAPT data by (i) applying a multiplicative correction for flux balance; (ii) applying a spatial smoothing to the maps; and (iii) remapping to a $360\times 180$ grid in $\phi$ and $\cos\theta$. The time derivative $\partial_{t}B_{r}$ is then estimated using a Savitzky-Golay smoothing filter (of total width 18 hrs), generating the map shown in Figure 10. The position of the data-assimilation window is evident in the map of $\partial_{t}B_{r}$, with the largest contribution coming from the main active region AR12209. (This is not newly emerging, but has significantly changed its structure since it was last observed on the previous Carrington rotation.)

In Figure 11, we show the inductive and sparse solutions for ${\mathbf{E}}_{\perp}$ in this case. Since there is such a dominant source of $\partial_{t}B_{r}$, there are significant non-local electric fields in the inductive solution, in contravention of Ohm’s law. The sparse $E_{\phi}$ is rather better, but the sparse $E_{\theta}$ is still not correctly localized. The cause of this is the flux imbalance in the original ADAPT map for $\partial_{t}B_{r}$, as in the previous test (Figure 8). Again the global flux correction has not really helped the problem. We conclude from these tests that lack of local flux-balance is the main obstacle to driving coronal MHD models with data-assimilative magnetic maps.

6 Conclusion

In this paper, we have shown how a sparse reconstruction technique based on the idea of searching for a localized solution allows one to recover accurate Ohm’s law ${\mathbf{E}}_{\perp}$ based on seemingly insufficient data of just $\partial_{t}B_{r}$. The technique is potentially useful for driving coronal MHD simulations from time-sequences of photospheric $B_{r}$ maps.

However, we have also identified a difficulty that can arise if one tries to compute the sparse electric field from maps where observational magnetogram data have been assimilated (e.g., Schrijver & De Rosa, 2003; Upton & Hathaway, 2014; Hickmann et al., 2015). Namely, errors in $\partial_{t}B_{r}$ that lead to local flux imbalance can prevent the sparse solution from being a good approximation to the real ${\mathbf{E}}_{\perp}$ (meaning the ${\mathbf{E}}_{\perp}$ that satisfies both Faraday’s law and Ohm’s law). This implies that the $\partial_{t}B_{r}$ map must first be corrected if ${\mathbf{E}}_{\perp}$ is to be reconstructed successfully. We have shown in Figure 9 that straightforward “global” methods of correcting the flux – whether additive or multiplicative – are insufficient to remove the problem. But this does not mean that a more sophisticated pre-processing method could not be successful; this is a possible direction for future research.

The problem of flux balance is probably best addressed at source, when the $B_{r}$ maps themselves are produced. This is easy to achieve in an idealized flux-transport model, but not in one with direct assimilation of observed magnetogram data. Yeates et al. (2015) have introduced a flux-transport model where entire, flux-balanced, bipolar regions are assimilated based on synoptic magnetograms (as used for our test case in Section 4.2). But it remains an open problem to assimilate higher-resolution magnetogram data in a manner that enforces, as far as possible, local flux-balance in $B_{r}$.

The author is grateful to Lockheed Martin and UC Berkeley for hosting my research visit in 2015, and in particular to Mark Cheung for suggesting the idea of looking for sparse solutions. I also thank Zoran Mikić and Cooper Downs (Predictive Science Inc.), and Roger Fletcher (University of Dundee, now sadly deceased) for useful discussions. The work was supported by a grant from the US Air Force Office of Scientific Research (AFOSR) in the AFOSR Basic Research Initiative “Understanding the Interaction of CMEs with the Solar-Terrestrial Environment.” I thank Carl Henney for supplying the ADAPT maps. This work utilises data obtained by the Global Oscillation Network Group (GONG) program, managed by the National Solar Observatory, which is operated by AURA, Inc. under a cooperative agreement with the National Science Foundation. The data were acquired by instruments operated by the Big Bear Solar Observatory, High Altitude Observatory, Learmonth Solar Observatory, Udaipur Solar Observatory, Instituto de Astrofísica de Canarias, and Cerro Tololo Interamerican Observatory.