Sparse reconstruction of electric fields from radial magnetic data

0.96 MB·

by  A. R. Yeates

Draftversion January25,2017
PreprinttypesetusingLATEXstyleAASTeX6v.1.0
SPARSE RECONSTRUCTION OF ELECTRIC FIELDS FROM RADIAL MAGNETIC DATA
Anthony R. Yeates
DepartmentofMathematical Sciences
DurhamUniversity
Durham,DH13LE,UK
7
1
0
ABSTRACT
2
Accurate estimates of the horizontal electric field on the Sun’s visible surface are important not only
n
a forestimatingthePoyntingfluxofmagneticenergyintothecoronabutalsofordrivingtime-dependent
J magnetohydrodynamicmodels of the corona. In this paper, a method is developedfor estimating the

4 horizontal field from a sequence of radial-componentmagnetic field maps. This problem of inverting
2
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
]
R Ohm’slaw. Ournewmethodgeneratesinsteadalocalizedsolution,usingabasispursuitalgorithmto
S find a sparse solution for the electric field. The method is shown to perform well on test cases where
h. the input magnetic maps are flux balanced, in both Cartesianand sphericalgeometries. However,we
p showthatiftheinputmapshaveasignificantimbalanceofflux–usuallyarisingfromdataassimilation
- – then it is not possible to find a localized, realistic, electric field solution. This is the main obstacle
o
r to driving coronal models from time sequences of solar surface magnetic maps.
t
s Keywords: magnetohydrodynamics, Sun: activity, Sun: evolution, Sun: magnetic fields, Sun: photo-
a
[ sphere

1
v 1. INTRODUCTION pitfalls of this technique.
0 The fundamental problem of electric field inversion
Magneto-hydrodynamic (MHD) simulations of the
8
fromFaraday’slaw(1)islackofuniqueness. Ifwewrite
7 magnetic field in the Sun’s corona are typically driven
6 by an imposed evolution on the solar surface. Since the E⊥ as a Helmholtz decomposition
.0 magnetic field B evolves according to Faraday’s law E⊥ =−∇× Φer −∇⊥Ψ, (3)
1
0 ∂tB=−∇×E, (1) then ∇×E⊥ depends only(cid:0) on(cid:1)the potential Φ(θ,φ,t),
7
and is independent of the choice of the secondpotential
1 the required boundary condition is the horizontal elec-
v: tricfieldE⊥ =Eθ(θ,φ,t)eθ+Eφ(θ,φ,t)eφ,writtenhere Ψ(θ,φ,t). Note thatournotationinterchangesΦ andΨ
compared to Miki´c et al. (1999). Given that ∂ B is
i in spherical polar coordinates. Unfortunately, this elec- t r
X independent ofΨ, the simplest wayto obtaina solution
tricfieldcannotbeobserveddirectly,butmustberecon-
ar structedfromother observations. If weassume thatthe consistentwith the givenBr data is to set Ψ≡0. Then
plasmaobeystheidealOhm’slawE=−v×BthenE⊥ (2) becomes the Poisson equation
can, inprinciple, be computed from vectorobservations −∇2Φ=∂ B , (4)
t r
of B and of the plasma velocity v.
In practice, however, such observations are not rou- which may be solved uniquely for Φ given appro-
tinely available for the full solar surface. Many authors priate boundary conditions. We call the resulting
have therefore opted to estimate E⊥ purely by invert- E⊥ the inductive solution for the electric field, and
ing Equation (1), typically using only a time sequence denote it by E⊥0. This solution has been used
of B (θ,φ,t) data (Miki´c et al. 1999; Wu et al. 2006; to drive coronal models in numerous studies (e.g.,
r
Mackay et al.2011;Cheung & DeRosa2012;Yang et al. Miki´c et al.1999;Amari et al.2003;Mackay et al.2011;
2012). In this case, only the radial component of (1) is Cheung & DeRosa 2012; Gibb et al. 2014; Yang et al.
used, 2012; Feng et al. 2015). Unfortunately, there is good
∂tBr =−er·∇×E⊥. (2) reason to believe that the non-inductive potential Ψ is
non-negligible on the real solar surface.
The present paper seeks to address one of the potential One problem with the inductive solution is its lack of
2
∂B Φ E E
t z x0 y0
0.04 0.04
2 0.1 0.05
0.02 0.02
1
y 0 0 0 0 0
-1
-0.02 -0.02
-2 -0.1 -0.05
-0.04 -0.04
-2 0 2 -2 0 2 -2 0 2 -2 0 2
x x x x
Figure 1. Thesimple solution in Equation (7),showing ∂tBz,theinductivepotential Φ, andtheresulting Ex0, Ey0. The color
scales for ∂tBz, Ex0, and Ey0 are clipped at one fifth of their maximum.
localization. If the real E⊥ satisfies Ohm’s law, then source. The corresponding inductive electric field E⊥0
it ought to vanish outside of patches of strong B . As therefore decays only as r−1. In particular, it may be
r
discussed in Section 2, this is not the case with the in- non-zero well outside the regions of non-zero B .
z
ductive solution; the main focus of this paper will be to Asanexplicitexample,considerabipolardistribution
choose Ψ in order to correct this deficiency.
(x+ρ)2+y2 (x−ρ)2+y2
Asecondproblemwiththeinductivesolutionisthatit ∂ B =exp − −exp − .
t z (cid:18) δ2 (cid:19) (cid:18) δ2 (cid:19)
cannotdetectelectricfieldsassociatedwithplasmaflows
(6)
along contours of B , since ∂ B = 0 in that case. The
r t r In this case, one may solve Equation (4) exactly to ob-
corresponding contribution to Ψ may be added if the
tain the closed-form solution
flowsareknown. Forexample,Kazachenko et al.(2014,
2015) have made this correction using observations of Φ(x,y)= δ2 log r+ + δ2E r+2 − δ2E r−2 ,
plasmavelocitiesinactiveregions,andtheresultingE⊥ 2 (cid:18)r−(cid:19) 4 1(cid:18)δ2(cid:19) 4 1(cid:18)δ2(cid:19)
has been used to drive magneto-frictional simulations (7)
(Fisher et al.2015). Ontheotherhand,Weinzierl et al. where r± = (x±ρ)2+y2 and E1(x) := x∞e−t/tdt
(2016)haveincludedthecontributiontoE⊥ fromapre- is the exponepntial integral. The logarithmicRbehaviour
scribed large-scale differential rotation in their global isclearfromthisexpression,andthecorrespondingelec-
magneto-frictional simulations. In this paper, however, tric field components are plotted in Figure 1.
we will put this issue aside and focus only on the local-
ization problem. 2.2. Variational formulation
InSection2,wedescribethelocalizationproblemwith Another way to characterise the inductive solution
the inductive solution, before characterizing the induc- E⊥0 is by the fact that it minimizes the L2-norm
tive solution in a variational framework. This makes kE⊥k2 := S|E⊥|2dS 1/2 among all possible solutions
the link to our proposed sparse electric field solution,
to (2). Th(cid:0)iRs is straigh(cid:1)tforward to see by inserting the
described in Section 3, that is designed to restore local-
expression (3) into the norm, which leads to
ization. Numerical tests in both Cartesian and spheri-
calgeometryarepresentedin Section4, while Section 5 kE⊥k22 =kE⊥0k22+k∇⊥Ψk22+2 ΨE⊥0·ndS.(8)
considers the limitations of the sparse solution. I∂D
On a periodic or infinite domain, the boundary term
2. INDUCTIVE ELECTRIC FIELD
vanishes, and it follows that kE⊥k2 is minimized by
2.1. The localization problem choosing Ψ ≡ 0. In this sense, the inductive electric
Since the inductive potential Φ solves the Poisson field E⊥0 is the smoothest possible E⊥ satisfying Fara-
equation (4), the solution may be expressed in terms day’s law for a given ∂tBz.
of the Green’s function as
∞ ∞ 2.3. Discrete formulation
1
′ ′
Φ(x,y,t)= ∂ B (x,y ,t)
2π Z−∞Z−∞ t z Inpractice,weworkwithE⊥ definedonadiscretenu-
×log (x−x′)2+(y−y′)2 1/2dx′dy′. (5) merical grid, so it is useful to formulate the analogous
discrete problem. In this paper, we discretise E⊥ on a
(cid:12) (cid:12)
Here we have ta(cid:12)ken Cartesian coordi(cid:12)nates and an infi- staggered grid (Yee 1966), where E is defined on the
x
nite plane for simplicity, but similar solutions hold on a horizontal cell edges, E on the vertical cell edges, and
y
sphericalsurface. Fora localizedsource∂ B ,we there- ∂ B at the cell centers (Figure 2). Such grids are com-
t z t z
fore have that Φ ∼ log(r) for large distance r from the monly used in MHD simulations. We discretise Fara-
3
written in terms of the Moore-Penrose pseudo-inverse
as x = A⊤(AA⊤)−1b, although in practice it is much
more efficient to use the fast-Poisson solver. Neverthe-
less, viewing the inductive solution as the ℓ -minimum
2
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 x.
3. SPARSE ELECTRIC FIELD
Our idea is to find a sparse solution for x (i.e., E ,
x
E )thatminimizesthenumberofnon-zerovalues. This
y
should be more localized than the inductive solution
E⊥0. Since it will differ from E⊥0, this new solution
Figure 2. The staggered grid, where Bz and Φ are defined willhaveanon-zeropotentialΨ. However,wewillwork
at cell centres and Ex and Ey on corresponding cell edges. directly with the system Ax = b (Equation 9), rather
than solving for Ψ (or Φ) explicitly.
day’s law (2) using Stokes’ theorem, so that the equa- Instead of minimizing the ℓ -norm of x, as in the in-
2
tion ductive solution, we propose to minimize the ℓ -norm.
1
∂Bi,j In other words, minimize
∆x∆y z =∆xEi,j+1/2−∆xEi,j−1/2
∂t x x kxk := |x | subject to Ax=b. (13)
+∆yEyi−1/2,j −∆yEyi+1/2,j (9) 1 Xi i
must hold for each grid cell i,j = 1,...,n. (In more In the optimization literature, problem (13) is known
generalcurvilinearcoordinates,theedgelengthsandcell as basis pursuit (Chen et al. 2001). It is used in a
areas would be different for each cell.) This constitutes wide variety of fields, and has recently been applied to
a system of linear equations for the unknowns E and the determination of differential emission measures in
x
E on each edge. However, the system is highly under- the solar corona (Cheung et al. 2015). Although min-
y
determined,sincethereare2n(n+1)unknownsbutonly imizing the ℓ -norm is not strictly equivalent to min-
1
n2 equations. This non-uniqueness of solution reflects imizing the number of non-zero components of x, the
our freedom in choosing the non-inductive component. minimum-ℓ solution to an under-determined system of
1
In the discrete case, the inductive solution may be linear equations is often the sparsest solution to that
computed by writing system(Candes et al.2006;Donoho & Tsaig2008). Nu-
Ei,j+1/2=(Φi,j −Φi,j+1)/∆y, (10) merically, it is preferable as it is a convex optimisation
x0 problem that can be efficiently solved using linear pro-
Eyi+01/2,j =(Φi+1,j −Φi,j)/∆x, (11) gramming (Gill et al. 2011). For the examples in Sec-
tion 4, we used the numerical implementation of basis
where Φi,j is defined at the cell centres. Then (9) be-
pursuit in the SparseLab1 library. This implementation
comes
of basis pursuit uses a primal dual method (Chen et al.
∂Bi,j 1
z = 2Φi,j −Φi,j+1−Φi,j−1 2001) and has a for the optimization.
∂t ∆y2
(cid:16) (cid:17) ForallofthetestcasesinSections4and5,withtoler-
+ 1 2Φi,j −Φi−1,j −Φi+1,j , (12) ance10−12,the optimizationtookonly10to20seconds
∆x2(cid:16) (cid:17) onadesktopworkstation,notincludingthetimetocon-
whichissimplythestandard5-pointstencilforthePois- struct the matrix A (which is only neededonce for each
sonequation. Intheexamplesbelow,wesolvethisusing grid). This is certainly fast enough for practical appli-
a standard fast-Poissonsolver (Press et al. 1992). cation to a time sequence of ∂ B maps, as would be
t r
However, it is instructive to think of the system of needed for driving a coronal MHD simulation.
equations (9) in the more abstract form Ax=b, where
x is the vector of unknowns (E , E ), b is the vector 4. NUMERICAL TESTS
x y
of ∂ B values, and A is the n2×2n(n+1) matrix cor-
t z In this section, we compare the inductive and sparse
responding to equations (9). The inductive solution is E⊥ solutionsfortwotestcases. Inthefirstcase(Section
then the solution to Ax = b that minimizes the dis-
crete ℓ -norm kxk2 := x2. In other words, it is the
2 2 i i
(unique)least-squaressoPlutiontothisunder-determined 1 SparseLab is freely available for download from
system of equations. As such, the solution may be http://sparselab.stanford.edu.
4
Target Ex ×10-3 Target Ey ×10-3 Input ∂tBz
1
4 4 0.015
0.5 0.01
2 2
0.005
y 0 0 0 0
-0.005
-2 -2
-0.5 -0.01
-4 -4 -0.015
-1
Inductive Ex0 [3.73e-03] ×10-3 Inductive Ey0 [1.25e-02] ×10-3 Inductive ∂tBz [3.31e-15]
1
4 4 0.015
0.5 0.01
2 2
0.005
y 0 0 0 0
-0.005
-2 -2
-0.5 -0.01
-4 -4 -0.015
-1
Sparse Ex [4.31e-09] ×10-3 Sparse Ey [2.29e-05] ×10-3 Sparse ∂tBz [8.85e-09]
1
4 4 0.015
0.5 0.01
2 2
0.005
y 0 0 0 0
-0.005
-2 -2
-0.5 -0.01
-4 -4 -0.015
-1
-1 0 1 -1 0 1 -1 0 1
x x x
Figure 3. Electric fieldreconstructions from (14)with δ2 =0.05, discretized with grid size n=256 (onlypart of thedomain is
shown). The top row shows the target solution, middle row the inductive solution, and bottom row the sparse solution (with
tolerance 10−12). The numbersin brackets are maximum absolute errors.
4.1), the true E⊥ is knownthroughOhm’s law,while in ary conditions are applied in both directions.
the second (Section 4.2) it is not. Moreover, the first Our basic solution is the bipolar distribution
case uses Cartesian geometry and the second case uses
x2+y2
spherical geometry. The additional problems that arise ∂ B =xexp − . (14)
t z (cid:18) δ2 (cid:19)
when applying the technique to real magnetic data are
discussed in Section 5. Physically, this could arise from uniform advection of a
single magneticpolarityB =exp[−(x2+y2)/δ2]under
z
4.1. Bipolar distribution an ideal Ohm’s law E⊥ = −v⊥×(Bzez) with velocity
Our first test case is based on the simplest configu- v⊥ =δ2/2ex, in which case
ration that may arise from Ohm’s law. In fact, due to
δ2 x2+y2
the linearityofboth Ohm’s law andEquation(2), more E⊥ = exp − . (15)
2 (cid:18) δ2 (cid:19)
general test configurations may readily be obtained by
superimposing multiple copies of this basic configura- (In this case, we are taking a snapshot at t = 0, as
tion. For simplicity, the computation is done in coor- the polarity is centred onthe origin.) Alternatively, the
dinates, in a square domain |x| < 3, |y| < 3, using the same ∂ B (with different B ) might correspond to the
t z z
discretization described in Section 2.3. Periodic bound- emergence of a bipolar magnetic field from the solar in-
5
1
20
0.5
10
θ
os 0 0
max error10−5 EE∂Bxyz/∂t c-0-.51 (a)Br[G] --2100
1 ×10-5
1.5
1
0.5
10−10 θ 0.5
10−10 10−5 os 0 0
tolerance c -0.5
-0.5 -1
Figure 4. Convergence of the sparse solution for (14), as a (b)∂tBr[Gs−1] -1.5
function of the tolerance parameter (x-axis) and grid reso- -1
lution(colors). Themaximumabsoluteerrors areshownfor 1 ×10-5
Ex (solid lines/squares), Ey (dashed lines/diamonds) and 1.5
∂tBr (dot-dashed lines/asterisks). The colors refer to grids 0.5 1
with n=32 (blue), n=64 (red), n=128 (green), n=256 0.5
θ
(magenta), and n=512 (cyan). os 0 0
c -0.5
-0.5
-1
terior. ThetargetE⊥givenby(15),alongwiththe∂tBz -1 (c)∂tB5r0+noi1se00[Gs−1]150 200 250 300 350 -1.5
in (14), are illustrated in the top row of Figure 3. φ
The middle row of Figure 3 shows the inductive so-
Figure 5. Inputdataforthesphericaltestcase,showingthe
lution E⊥0 for this ∂tBz, which is qualitatively similar distributionsof (a)Br and(b)∂tBr from thefluxtransport
to the illustrative solution in Section 2. Not only is model. Panel (c) shows the map of ∂tBr with added noise
E⊥0 not localized within the region of non-zero ∂tBz, tohfamtawxeimuusme ffoorrt(hae),taesntd. ±T5h%e mofampsaxarimeusamtufroarte(bd)aatn±d2(0c%).
but the topology is wrong: there is a substantial E
x0
component, despite the fact that E = 0 for the target
x the numerical code is freely available
solution. Thesparsesolutionusingbasispursuit,shown
(https://github.com/antyeates1983/sft_data).
in the bottom row of Figure 3, is substantially more ac-
As input data, we used synoptic maps of B from
r
curate. Todemonstrateconvergence,Figure4showsthe
the Global Oscillation Network Group (GONG,
absolute errors for the sparse solution, as a function of
gong.nso.edu/data/magmap/), starting in Carrington
both grid resolution and of the tolerance parameter in
rotation 2073 and generating the output snapshot in
the basis pursuit algorithm. Firstly, the error in ∂ B
t z rotation 2109. As described by Yeates et al. (2015),
isindependentoftoleranceorgridresolution,indicating
the input maps are used (a) to initialize B at the
that the constraint Ax = b is preserved to high accu- r
beginning of the simulation, and (b) to extract strong
racy in all cases. Secondly, there is convergencein both
flux regions used to update B during the evolution.
r
E and E as the tolerance is reduced. The error in E
x y x The computation used a 360×180 grid, equally spaced
is independent of resolution, whereas that in E satu-
y in longitude (φ) and sine latitude (cosθ).
rates at a (higher) level depending on grid resolution.
For context, Figure 5(a) shows the B map from the
r
This simply reflects the truncation error in the numer-
flux-transport simulation, while Figure 5(b) shows the
ical approximation (9) for the curl, which is quadratic
∂ B mapfromthesametime. Sincethesimulationuses
t r
in ∆x. The saturation is not seen in E because of our
x a supergranular diffusion rather than imposing convec-
particular target solution E =0.
x tive velocities, the map appears smoother than would a
real observed map. To make the ∂ B map look more
4.2. Spherical distribution t r
realistic, we have added random noise to produce the
To demonstrate how the sparse reconstruction per- map in Figure 5(c). The noise is carefully constructed
forms on a more realistic example, we have taken a topreservelocalfluxbalance,andisgeneratedbyadding
snapshot from a flux-transport simulation of Br(θ,φ,t) a bipolar distribution
in spherical coordinates, covering the full solar surface.
(i−i )2+(j−j )2
Thereasonforusingasimulatedmapratherthananob- βi0,j0(i−i )exp − 0 0 (16)
0 (cid:18) δi0,j0 (cid:19)
served synoptic magnetogram is to ensure perfect flux
balance;theconsequencesofnothavingfluxbalancewill centredateachpixel(i ,j )onthe(s,φ)grid. Themag-
0 0
become evident in Section 5. nitude βi0,j0 at each pixel is chosen from a normal dis-
We used the flux-transport model de- tribution with mean zero and σ =0.005max |∂ B |,
i0,j0 t r
scribed by Yeates et al. (2015), for which and the size δi0,j0 is either 1 or 2 pixels, with equal
6
InductiveEθ0[Vm−1] InductiveEφ0[Vm−1] Inductive∂tBr[Gs−1][3.22e-07] ×10-5
1
1.5
0.05 0.05 1
0.5
0.5
θ
s 0 0 0 0
o
c
-0.5
-0.5
-0.05 -0.05 -1
-1.5
-1
SparseEθ[Vm−1] SparseEφ[Vm−1] Sparse∂tBr[Gs−1][2.97e-08] ×10-5
1
1.5
0.05 0.05 1
0.5
0.5
θ
s 0 0 0 0
o
c
-0.5
-0.5
-0.05 -0.05 -1
-1.5
-1
100 200 300 100 200 300 100 200 300
φ φ φ
Figure 6. Electric field reconstructions for the spherical test case (Figure 5c). The top row shows the inductive solution and
the bottom row the sparse solution (with tolerance 10−12). For comparison, all electric fields are saturated at ±20% of the
maximuminductive|Eθ|,while∂tBr issaturatedat±5%ofthemaximuminputvalue. Thenumbersinbracketsgivemaximum
absolute (discretization) errors.
probability. In addition, we rotate the pattern by 90 5.1. Diffuse distributions of ∂ B
t r
degrees at each particular pixel with equal probability.
The sparse solution makes the physical assumption
This generates a more realistic map as shown in Figure that E⊥ should be localized because it satisfies Ohm’s
5(c). law for a localized B. If the magnetic field distribution
Figure 6 shows the inductive and sparse reconstruc-
is too diffuse, then we cannotexpect the assumptionof
tions of E⊥ from the map in Figure 5(c). To account localization to recover the correct (diffuse) E⊥.
for the spherical coordinates, Equations (10)–(12) have
To illustrate what happens if we nevertheless try to
been modified to include the necessary coordinate fac-
applythesparsereconstructioninsuchasituation,con-
tors. The right-hand column of Figure 6 shows that
siderthehorizontaldiffusivespreadingofaninitialgaus-
both methods preserve ∂ B to high accuracy, as in the
t r sian magnetic polarity
Cartesiantest. UnlikeintheCartesiantest,however,we
no longer compare to a target E⊥, since the E⊥ corre- B (x,y,t)= a2 exp − x2+y2 . (17)
sponding to our added noise is unknown. Nevertheless, z a2+4ηt (cid:18) a2+4ηt(cid:19)
that the sparse reconstruction is the inductive recon-
struction. the ∂ B distribution is more complex than From the resistive Ohm’s law E⊥ = η∇×(Bzez), we
t r get
in the previous test, the sparse E⊥ is still much local-
ized within regions of strong ∂ B . Moreover, it has 2a2y x2+y2
t r
E =− exp − , (18)
a weaker Eφ than Eθ, consistent with the fact that the x (a2+4ηt)2 (cid:18) a2+4ηt(cid:19)
dominant Ohm’s law contribution in the flux-transport 2a2x x2+y2
simulation was from differential rotation (vφ). This is Ey = (a2+4ηt)2 exp(cid:18)−a2+4ηt(cid:19). (19)
notthe casein the inductive solution. Thuswe arecon-
fident that the sparse solution not just for simple test Faraday’s law then gives
cases,butalsoformorecomplex∂ B distributionsthat
t r x2+y2 1
are qualitatively similar to those on the Sun. ∂ B =4a2η −
t z (cid:20)(a2+4ηt)3 (a2+4ηt)2(cid:21)
x2+y2
×exp − . (20)
5. LIMITATIONS (cid:18) a2+4ηt(cid:19)
As with any inverse problem, care is needed in the For this exercise, we fix a = 0.5 and ηt = 0.1. This
application of our sparse reconstruction technique. In solution is shown in the top row of Figure 7.
this section, we consider two limitations which can be It is important to note that, in this case, Ohm’s
important in the practical application to the driving of law is compatible with a purely inductive solution, be-
coronal MHD models. cause E⊥ from Ohm’s law has the inductive form with
7
Target E Target E Input ∂ B
x y t z
3 0.01 0.01
0.04
2
0.005 0.005
0.02
1
y 0 0 0 0
-1
-0.02
-0.005 -0.005
-2
-0.04
-3 -0.01 -0.01
Inductive E [1.68e-05] Inductive E [1.68e-05] Inductive ∂ B [3.14e-06]
x0 y0 t z
3 0.01 0.01
0.04
2
0.005 0.005
0.02
1
y 0 0 0 0
-1
-0.02
-0.005 -0.005
-2
-0.04
-3 -0.01 -0.01
Sparse E [3.79e-02] Sparse E [3.79e-02] Sparse ∂ B [3.17e-06]
x y t z
3 0.01 0.01
0.04
2
0.005 0.005
0.02
1
y 0 0 0 0
-1
-0.02
-0.005 -0.005
-2
-0.04
-3 -0.01 -0.01
-2 0 2 -2 0 2 -2 0 2
x x x
Figure 7. Electric field reconstructions from (20), discretized with grid size n = 256. The top row shows the target solution,
middle row the inductive solution, and bottom row the sparse solution (with tolerance 10−12). The numbers in brackets are
maximum absolute errors.
Φ(x,y,t)=ηB (x,y,t). Indeed, the middle row of Fig- Sun, given high enough resolution. The flux-transport
z
ure 7 showsthe inductive solutionto reproducethe tar- test case in Section 4.2 showed the method to work
get E⊥ in this case. However, the bottom row shows for realistic solar magnetic maps, at the resolutions of
that the sparse solution using basis pursuit does not present-daysimulations. However,thereisanothermore
converge to the target. Rather, it favors a concentra- subtle problem with real data, that we describe in the
tion of E along vertical lines, and E along horizontal following section.
x y
lines. ThisbehavioristypicalwhenthetargetE⊥ istoo
diffuse. It is not possible for the total spatial extent of 5.2. Flux imbalance in ∂ B maps
t r
E⊥ to be more localized than ∂tBr, owing to the con- For a localized E⊥ to exist, it is necessary not only
straint of satisfying (2). But it is possible for more of
that ∂ B is localized, but that the net flux ∂ B dA
the electric field to be concentrated in thinner regions, t r S t r
vanishes over the region S of localization. TRo see this,
and this is what minimizing the ℓ -norm does.
1 apply Stokes’ theorem on some closed curve that encir-
Thislimitationisanobviousone,andlimitsthesparse
cles S. If the net flux in S is non-zero, then there must
solutiontechnique toinput mapswhereB issufficiently be non-zero E⊥ somewhere on the bounding curve, and
localized. Fortunately, this situation is typical on the
indeed on any curve enclosing non-zero net flux.
8
SparseEθ[Vm−1] SparseEφ[Vm−1] ∂Br/∂t[Gs−1] ×10-5
1
0.05 0.05
1
0.5
θ
s 0 0 0 0
o
c
-0.5 -1
-0.05 -0.05
-1
×10-5
1
0.05 0.05
0.5 1
θ
s 0 0 0 0
o
c
-0.5
-1
-0.05 -0.05
-1
×10-5
1
0.05 0.05
0.5 1
θ
s 0 0 0 0
o
c
-0.5
-1
-0.05 -0.05
-1
100 200 300 100 200 300 100 200 300
φ φ φ
Figure 8. Sparse electric field reconstructions from the spherical test case (Figure 5c), modified to simulate assimilation of a
new bipolar region. Each row simulates a different time of observation where the assimilation window (black outline in the
right-handcolumn)hasshiftedleftwardbyonedayofsolarrotation. Herenocorrectionforthefluximbalancein∂tBr hasbeen
applied.
This condition of vanishing net flux is satisfied by ∂ B , both in the assimilationwindow and in the entire
t r
both of our simple examples (14) and (20), and by ev- map. Correspondingly,itisimpossibletofindalocalized
ery local flux concentration in the flux-transport test electricfield,andindeedthesparsesolutionfails. Inthe
case (Section 4.2). But it need not be satisfied in an area of the new region, we see similar linear structures
observationally-derivedmap of ∂ B , even if a flux cor- to the diffusive example in Figure 7, where the flux was
t r
rectionhasbeenapplied. Wewilldemonstratethisprob- also not locally balanced. But notice that the electric
lem firstly in a controlled test and secondly in a “real” field is modified not only near to the new region, but
data-assimilative model. also at farther distances.
Forthecontrolledtest,wemodifytheinputmapfrom Of course, this is a rather unrealistic test. If this E⊥
Section 4.2 to simulate the typical situation where new was being used to drive a coronal MHD model, for ex-
magnetogram observations are assimilated only within ample,onewouldcorrectthe∂ B mapforfluxbalance,
t r
a finite region on the visible face of the Sun. The most before trying to compute E⊥. In Figure 9 we show the
severe problems of flux imbalance occur when a new effect on E of first making a flux-balance correction to
θ
active region has emerged on the far-side of the Sun, the ∂ B map(from the top rowofFigure 8). The mid-
t r
and is gradually assimilated into the magnetic map as dle and right plots show the results with two different
it rotates into view. methods of flux correction. In the “additive” method,
Figure 8 shows our test. Here the original ∂ B map, an equal amount is subtracted from each pixel in the
t r
from Figure 5(c), has been modified to simulate the 360×180 grid to remove the imbalance. In the “mul-
gradual assimilation of a new bipolar active region. In tiplicative” method, all of the pixels where ∂ B > 0
t r
the right-hand column of Figure 8, the assimilation re- are multiplied by a constant factor f , and all pixels
+
gion is indicated on the modified ∂ B maps. Each row where ∂ B <0 are multiplied by another constant fac-
t r t r
corresponds to a different time as the new active region tor f−. The factors f+, f− are chosen so that the net
rotatesinto view, with the times differing by one day of flux vanishes after scaling, and is equal to the mean of
solar rotation (27.2753◦ in these Carrington maps). the positive and negative fluxes before scaling. Com-
When the new region is fully contained within the paring the results in Figure 8, it is evident that neither
assimilationwindow (bottomrowof Figure8), itmakes method of flux correction gives much improvement in
onlyalocalizedcontributiontothesparseE⊥,correctly the reconstructedEθ (or Eφ). Only by limiting the flux
reproducing E⊥ as in Section 4.1. But when the region correction to the active region itself could the solution
isonlypartiallyobserved,thereisanunbalancedfluxof be improved.
9
Eθ[Vm−1]-nocorrection Eθ[Vm−1]-additivecorrection Eθ[Vm−1]-multiplicativecorrection
1
0.06 0.06 0.06
0.5 0.04 0.04 0.04
0.02 0.02 0.02
θ
s 0 0 0 0
o
c -0.02 -0.02 -0.02
-0.5 -0.04 -0.04 -0.04
-0.06 -0.06 -0.06
-1
100 200 300 100 200 300 100 200 300
φ φ φ
Figure 9. Sparse Eθ for different methods of global flux correction.
1 50
Similar behavior is found when we apply the sparse
reconstructiontechnique to atime sequenceofB maps Br[G] 40
r 30
generated by the Air Force Data-Assimilative Photo- 0.5 20
spheric Flux Transport (ADAPT) model (Arge et al. 10
θ
2010; Henney et al. 2012). This model assimilates ob- os 0 0
c -10
served magnetograms from the visible face of the Sun
-20
intoasurfaceflux-transportsimulation,soastoapprox- -0.5
-30
imate the globaldistribution of Br on the solarsurface, -40
asafunctionoftime. The dataassimilationmeansthat -1 -50
the electric field E⊥ is not known everywhere in the 1 ×10-4
map, so must be reconstructed, if these maps are used ∂tBr[Gs−1]
2
to drive coronal MHD simulations. 0.5
As an illustration, we choose an ADAPT run driven 1
by GONG magnetograms, and extract ∂ B for 2014 θ
t r os 0 0
November 16, 00:00 UT. This particular date has been c
-1
deliberatelychosensinceWeinzierl et al.(2016)founda
-0.5
significant non-localized E⊥0 caused by re-assimilation -2
of a large active region (see their Figure 11). Here,
-1 -3
we process the ADAPT data by (i) applying a multi- 50 100 150 200 250 300 350
φ
plicative correctionfor flux balance;(ii) applying a spa-
tial smoothing to the maps; and (iii) remapping to a Figure 10. Input from the ADAPTmodel, showing Br (for
context)andthecomputed∂tBr. Theplotshavebeensatu-
360×180 grid in φ and cosθ. The time derivative ∂tBr rated at ±50G and ±3×10−4Gs−1, respectively.
is then estimated using a Savitzky-Golay smoothing fil-
ter(oftotalwidth18hrs),generatingthemapshownin
Figure10. Thepositionofthedata-assimilationwindow
6. CONCLUSION
isevidentinthemapof∂ B ,withthe largestcontribu-
t r
Inthispaper,wehaveshownhowasparsereconstruc-
tioncomingfromthemainactiveregionAR12209. (This
tiontechniquebasedontheideaofsearchingforalocal-
is not newly emerging, but has significantly changed its
ized solution allows one to recover accurate Ohm’s law
structuresinceitwaslastobservedonthepreviousCar-
rington rotation.) E⊥ based on seemingly insufficient data of just ∂tBr.
The technique is potentially useful for driving coronal
In Figure 11, we show the inductive and sparse solu-
tions for E⊥ in this case. Since there is such a domi- MHD simulations from time-sequences of photospheric
B maps.
nantsourceof∂ B , there aresignificantnon-localelec- r
t r
However, we have also identified a difficulty that
tric fields in the inductive solution, in contravention of
can arise if one tries to compute the sparse elec-
Ohm’s law. The sparse E is rather better, but the
φ
tric field from maps where observational magnetogram
sparse E is still not correctly localized. The cause of
θ
data have been assimilated (e.g., Schrijver & De Rosa
this is the flux imbalance in the original ADAPT map
2003; Upton & Hathaway 2014; Hickmann et al. 2015).
for ∂ B , as in the previous test (Figure 8). Again the
t r
Namely, errors in ∂ B that lead to local flux imbal-
globalfluxcorrectionhasnotreallyhelpedtheproblem. t r
ance can prevent the sparse solution from being a good
We conclude from these tests that lack of local flux-
approximationtotherealE⊥ (meaningtheE⊥ thatsat-
balance is the main obstacle to driving coronal MHD
isfies both Faraday’s law and Ohm’s law). This implies
models with data-assimilative magnetic maps.
thatthe ∂tBr mapmustfirstbecorrectedifE⊥ isto be
reconstructed successfully. We have shown in Figure 9
10
InductiveEθ0[Vm−1] InductiveEφ0[Vm−1] Inductive∂tBr[Gs−1][9.81e-06] ×10-4
1 0.6 0.6
1.5
0.4 0.4
0.5 1
0.2 0.2 0.5
θ
s 0 0 0 0
o
c
-0.2 -0.2 -0.5
-0.5 -1
-0.4 -0.4
-1.5
-1 -0.6 -0.6
SparseEθ[Vm−1] SparseEφ[Vm−1] Sparse∂tBr[Gs−1][4.88e-07] ×10-4
1 0.6 0.6
1.5
0.4 0.4
0.5 1
0.2 0.2 0.5
θ
s 0 0 0 0
o
c
-0.2 -0.2 -0.5
-0.5 -1
-0.4 -0.4
-1.5
-1 -0.6 -0.6
100 200 300 100 200 300 100 200 300
φ φ φ
Figure 11. Inductive(top) and sparse (bottom) electric field solutions, for the∂tBr map derived from ADAPTBr maps. The
electric field components are saturated at ±20% and ∂tBr at ±10%.
that straightforward“global”methods ofcorrectingthe Berkeley for hosting my research visit in 2015, and in
flux – whether additive or multiplicative – are insuffi- particular to Mark Cheung for suggesting the idea of
cient to remove the problem. But this does not mean looking for sparse solutions. I also thank Zoran Miki´c
that a more sophisticated pre-processing method could andCooperDowns(PredictiveScienceInc.), andRoger
not be successful; this is a possible direction for future Fletcher (University of Dundee, now sadly deceased)
research. for useful discussions. The work was supported by a
The problem of flux balance is probably best ad- grant from the US Air Force Office of Scientific Re-
dressed at source, when the B maps themselves are search (AFOSR) in the AFOSR Basic Research Ini-
r
produced. This is easy to achieve in an idealized flux- tiative Understanding the Interaction of CMEs with
transport model, but not in one with direct assimila- the Solar-Terrestrial Environment. I thank Carl Hen-
tionofobservedmagnetogramdata. Yeates et al.(2015) ney for supplying the ADAPT maps. This work utilises
have introduced a flux-transport model where entire, dataobtainedbytheGlobalOscillationNetworkGroup
flux-balanced, bipolar regions are assimilated based on (GONG) program, managed by the National Solar Ob-
synopticmagnetograms(asusedforourtestcaseinSec- servatory, which is operated by AURA, Inc. under a
tion 4.2). But it remains an open problem to assimilate cooperative agreement with the National Science Foun-
higher-resolution magnetogram data in a manner that dation. Thedatawereacquiredbyinstrumentsoperated
enforces, as far as possible, local flux-balance in B . by the Big Bear Solar Observatory, High Altitude Ob-
r
servatory, Learmonth Solar Observatory, Udaipur Solar
Observatory, Instituto de Astrofsica de Canarias, and
The author is grateful to Lockheed Martin and UC Cerro Tololo Interamerican Observatory.
REFERENCES
Amari,T.,Luciani,J.F.,Aly,J.J.,Mikic,Z.,&Linker,J.2003, Feng, X.,Ma,X.,&Xiang,C.2015,JournalofGeophysical
ApJ,595,1231 Research(Space Physics),120,10
Arge,C.N.,Henney,C.J.,Koller,J.,etal.2010,Twelfth Fisher,G.H.,Abbett, W.P.,Bercik,D.J.,etal.2015,Space
International SolarWindConference, 1216,343 Weather, 13,369
Candes,E.J.,Romberg,J.K.,&Tao,T.2006,Communications
Gibb,G.P.S.,Mackay,D.H.,Green,L.M.,&Meyer,K.A.
onpureandappliedmathematics, 59,1207
2014,ApJ,782,71
Chen,S.S.,Donoho,D.L.,&Saunders,M.A.2001,SIAM
Gill,P.R.,Wang,A.,&Molnar,A.2011,IEEETransactionson
review,43,129
SignalProcessing,59,4595
Cheung,M.C.M.,Boerner,P.,Schrijver,C.J.,etal.2015,ApJ,
807,143 Henney, C.J.,Toussaint,W.A.,White,S.M.,&Arge,C.N.
Cheung,M.C.M.,&DeRosa,M.L.2012, ApJ,757,147 2012,SpaceWeather, 10,S02011
Donoho,D.L.,&Tsaig,Y.2008, IEEETransactionson Hickmann,K.S.,Godinez, H.C.,Henney, C.J.,&Arge,C.N.
InformationTheory,54,4789 2015,SoPh,290,1105