Low and Lou analytical NLFFF solution

The non-linear force-free equations cannot be solved exactly in general, so we have no simple formula for the coronal field. This is a problem when we want to test an extrapolation code: if we feed it a magnetogram and it returns a 3D field, how do we know the answer is right? The usual metrics only tell us the field is self-consistent (small force, small divergence), not that it matches the true field.

Low & Lou (1990) got around this by finding one special case that can be solved. They assumed the field is axisymmetric (the same all the way around an axis), which is enough to collapse the equations down to a single ordinary differential equation that we can solve on a computer to high accuracy. The trick at the end is to bury the solution below the surface and tilt it, so that on a flat photosphere it looks like an ordinary bipolar active region rather than something obviously symmetric.

NF2 ships this as get_analytic_b_field(...) in nf2/data/analytical_field.py. It is the synthetic partner to the real-data AR 11158 benchmark.

---

Setting up an axisymmetric field

We work in spherical coordinates $(r,\theta ,\phi )$ and assume nothing depends on the angle around the axis, so ${\partial }_{\phi }=0$. Any such field can be built from a single scalar function $A(r,\theta )$, called the flux function, plus a separate sideways (azimuthal) part $B_{\phi }$:

$$\mathbf{B}=\nabla A\times \nabla \phi +B_{\phi }\hat{\mathbf{\phi }}.$$

Writing this out gives the three components

$$B_r=\frac{1}{r^2\sin \theta }\frac{\partial A}{\partial \theta },B_{\theta }=-\frac{1}{r\sin \theta }\frac{\partial A}{\partial r},B_{\phi }=\frac{Q}{r\sin \theta },$$

where in the last term we have just renamed the azimuthal field, defining $Q\equiv r\sin \theta B_{\phi }$ because it turns out to be the more convenient variable. The reason for building the field from $A$ this way is that it is automatically divergence-free: with these definitions $\nabla \cdot \mathbf{B}=0$ holds exactly, so one of Maxwell’s conditions is satisfied before we do any work. These are the same formulas used in analytical_field.py:

Br     = (r**2 * np.sin(theta))**-1 * dA_dtheta
Btheta = -(r * np.sin(theta))**-1 * dA_dr
Bphi   = (r * np.sin(theta))**-1 * Q

---

Applying the force-free condition

A force-free field obeys $\nabla \times \mathbf{B}=\alpha \mathbf{B}$, where $\alpha$ is the torsion function and is constant along each field line. We now put our axisymmetric field into this condition. Because taking the curl swaps the in-plane and sideways parts of the field, the single vector equation splits into two pieces.

The sideways part gives one equation linking $A$ and $Q$:

$$\frac{{\partial }^{2}A}{\partial r^2}+\frac{\sin \theta }{r^2}\frac{\partial }{\partial \theta }\left(\frac{1}{\sin \theta }\frac{\partial A}{\partial \theta }\right)=-\alpha Q.$$

The in-plane part gives $\nabla Q=\alpha \nabla A$. This says $Q$ changes in exactly the same direction as $A$ everywhere, which can only happen if $Q$ depends on position only through $A$. In other words

$$Q=Q(A),\alpha =\frac{dQ}{dA}.$$

This is just the axisymmetric version of the familiar result that $\alpha$ is the same all along a field line. Substituting $\alpha =dQ/dA$ back into the first equation leaves a single equation for $A$ alone (this is a force-free Grad–Shafranov equation, the standard equation for axisymmetric magnetic equilibria):

$$\frac{{\partial }^{2}A}{\partial r^2}+\frac{\sin \theta }{r^2}\frac{\partial }{\partial \theta }\left(\frac{1}{\sin \theta }\frac{\partial A}{\partial \theta }\right)+Q\frac{dQ}{dA}=0.$$

So we have gone from a system of vector equations to one scalar equation, as long as we pick the function $Q(A)$.

---

Simplifying the angle

It helps to swap the angle $\theta$ for $\mu =\cos \theta$. The messy angular term then tidies up to $(1-{\mu }^2){\partial }^{2}A/\partial {\mu }^2$, and the equation becomes

$$\frac{{\partial }^{2}A}{\partial r^2}+\frac{1-{\mu }^2}{r^2}\frac{{\partial }^{2}A}{\partial {\mu }^2}+Q\frac{dQ}{dA}=0.$$

---

The Low & Lou guess

Low & Lou looked for a solution where the $r$ and $\mu$ dependence separate, of the form

$$A(r,\mu )=\frac{P(\mu )}{r^n},$$

for some number $n$ and some shape function $P(\mu )$. With this guess, the first two terms in the equation both fall off as $r^{-(n+2)}$. For the third term to fall off at the same rate (so the $r$‘s cancel and we are left with an equation in $\mu$ only), we have to choose

$$Q(A)=a{A}^{1+1/n},$$

for a constant $a$. Note this makes $\alpha =dQ/dA\propto A^{1/n}$, which varies in space — so this really is a non-linear force-free field, not the constant-$\alpha$ kind. In the code $Q=\sqrt{a^2}A|A{|}^{1/n}$, where a2 is $a^2$ and the absolute value keeps the fractional power well-defined when $A$ changes sign.

Putting the guess in and cancelling the $r$ dependence leaves a single ordinary differential equation for $P$:

$$(1-{\mu }^2)\frac{d^{2}P}{d{\mu }^2}+n(n+1)P+a^2\frac{n+1}{n}{P}^{1+2/n}=0.$$

This is the Low & Lou equation, and it is exactly what _differential_equation solves:

d2P_dmu2 = -(n*(n+1)*P + a2*(1+n)/n * P**(1+2/n)) / (1 - mu**2 + 1e-8)

(the +1e-8 just avoids dividing by zero at the poles, where $\mu =\pm 1$).

Boundary conditions and why $a$ is fixed

At the two poles ($\mu =\pm 1$) the field has to stay finite, which forces $P(-1)=P(+1)=0$. That is two conditions on a second-order equation, which is normally one too many. The catch is that $a$ is not free: only for certain special values of $a$ does a solution exist that is zero at both poles. Finding those values is the whole job. Because the equation is non-linear, the size of $P$ also matters, so we fix it by also setting the slope at one pole. NF2 uses

$$P(-1)=0,P(+1)=0,P^{\prime }(-1)=10,$$

and solves the problem numerically with scipy.integrate.solve_bvp, scanning for the allowed values of $a$. The integer $m$ counts how many times $P$ crosses zero in between, and selects which solution we want. The standard choice (and the NF2 default) is $n=1,m=1$: this gives a clean, dipole-like field with a single sign change.

---

Burying and tilting it

As derived, the solution sits at the centre of its coordinate system and is perfectly symmetric about the axis. That is not yet useful: the centre is a singular point, and a symmetric field is a poor test. Low & Lou fix both by moving the whole solution. Two numbers control this:

• $l$ - how far below the photosphere the centre sits. Burying it (NF2 default $l=0.3$) keeps the singular point out of the box, so the field is smooth everywhere we look.
• $\Phi$ (called psi in the code) — the tilt of the symmetry axis away from vertical. A tilt (NF2 default $\Phi =\pi /4$) is what makes the field look asymmetric on the flat surface, turning the symmetric solution into a realistic, sheared bipole.

In the code, each point $(x,y,z)$ in the box is shifted down by $l$ and rotated by $\Phi$ to find where it sits relative to the buried, axis-aligned solution:

X = x*np.cos(psi) - (z + l)*np.sin(psi)
Y = y
Z = x*np.sin(psi) + (z + l)*np.cos(psi)

We then read off $r$, $\theta$, $\phi$ in that frame, evaluate the field, and rotate the result back into the box.

---

Building the final field

With $P(\mu )$ and $a^2$ known, everything is now an explicit formula. We compute

$$A=\frac{P(\mu )}{r^n},\frac{\partial A}{\partial \theta }=-\frac{\sin \theta }{r^n}\frac{dP}{d\mu },\frac{\partial A}{\partial r}=-\frac{nP(\mu )}{r^{n+1}},Q=\sqrt{a^2}A|A{|}^{1/n},$$

get $(B_r,B_{\theta },B_{\phi })$ from the formulas at the top, convert these to Cartesian components, and finally rotate back through $\Phi$ into the box frame. The output is a 3D array $\mathbf{B}(x,y,z)$ that, to numerical precision, is an exact non-linear force-free field everywhere.

---

How NF2 uses it

The test is simple: build the full cube, hand the code only the bottom slice $\mathbf{B}(x,y,0)$ as a fake magnetogram, let it extrapolate, and compare its 3D answer against the cube we already know (using $C_{\text{vec}}$, energy ratios, and the other metrics).

This is a stronger check than the force-free and divergence residuals NF2 minimises during training. A small residual only says the field is self-consistent; Low & Lou lets us confirm it is the correct field. For example, a current-free potential field can have a small residual yet completely miss the twist and free energy that the true field carries.

Its weakness as a test is that it is very clean: one smooth buried source, no noise, and perfectly force-free boundary data. Passing it shows the solver’s maths is right, but not that it can cope with the messy, partly non-force-free photospheric data of a real active region like AR 11158.

---

Related

• From MHD to the force-free equations — where $\nabla \times \mathbf{B}=\alpha \mathbf{B}$ comes from
• Force-free and NLFFF coronal fields — the equations this solution satisfies exactly
• Classical NLFFF extrapolation methods — the codes the benchmark tests
• NLFFF quality metrics — the metrics used in the comparison
• NOAA AR 11158 as an NF2 benchmark — the real-data partner benchmark
• Potential field extrapolation — the simpler field Low & Lou is not
• NF2 codebase pipeline — where get_analytic_b_field sits in the data path