**Abstract**

**|SOURCE|** Proton radiography is an important diagnostic method for laser plasma experiments and is particularly important in the analysis of magnetized plasmas. The theory of radiographic image analysis has heretofore only permitted somewhat limited analysis of the radiographs of such plasmas. We furnish here a theory that remedies this deficiency. We show that to linear order in magnetic field gradients, proton radiographs are projection images of the MHD current along the proton trajectories. We demonstrate that in the linear regime (i.e., the small image contrast regime), the full structure of the projected perpendicular magnetic field can be reconstructed by solving a steady-state inhomogeneous 2-dimensional diffusion equation sourced by the radiograph fluence contrast data. We explore the validity of the scheme with increasing image contrast, as well as limitations of the inversion method due to the Poisson noise, discretization errors, radiograph edge effects, and obstruction by laser target structures. We also provide a separate analysis that is suited to the inference of isotropic-homogeneous magnetic turbulence spectra. Finally, we discuss extension of these results to the nonlinear regime (i.e., the order unity image contrast regime).

**I. INTRODUCTION**

Proton radiography is an experimental technique wherein the electric and magnetic fields in a plasma are imaged using beams of protons. The deflection by electric and magnetic fields of a proton’s path bears information about the morphology and strengths of electric and magnetic fields. This kind of imaging is in use at laser facilities around the world and has advanced our understanding of flows, instabilities, and electric and magnetic fields in high energy density plasmas.

Two distinct types of proton sources have been developed for use in proton imaging. The first type of source uses thermally produced protons accelerated by the extreme electric field gradients generated when a target is illuminated by an intense laser.1–4 The second type of source uses mono-energetic, fusion-produced protons created by the implosion of a D3He capsule.5–8

The ability of proton imaging to provide information about the strength and morphology of electric and magnetic fields is exceptional. Such fields are ubiquitous in high energy density physics (HEDP) experiments, since laser-driven flows create strong electric fields and the mis-aligned electron density and temperature gradients in such flows create strong magnetic fields via the Biermann battery effect.9 Thermally produced protons have been used, for example, to detect electric fields in laser-plasma interaction experiments1 and to image solitons,10 laser-driven implosions,2 and collisionless electrostatic shocks.11 Mono-energetic protons have been used, for example, to image, quantify, and distinguish between electric and magnetic fields in laser-driven plasmas from foils5 and inside hohlraums;12 to study changes in magnetic field topology due to reconnection,13 the morphology of the magnetic fields generated in laser-driven implosions,14–16 the evolution of filamentation and self-generated fields in the coronae of directly driven inertial-confinement-fusion capsules,17 and the Rayleigh-Taylor induced magnetic fields;18 and to observe magnetic field generation via the Weibel instability in collisionless shocks.19

In a few laser plasma experiments, the strength of the projected perpendicular magnetic field was estimated by placing a grid between the proton source and the detector and using the distorted image of the grid to measure the deflection angle of the protons.5 Prior to the publication of the paper by Kugland et al.20 (hereafter K2012), essentially no quantitative information could be extracted from proton radiographs. Thus, it was about a decade before any theory was developed regarding how proton radiography images could be analyzed to recover the electromagnetic structure of the imaged plasma.

The K2012 paper furnished the first systematic quantitative discussion of the analysis of proton radiographs, identifying the principal physical parameters, and identifying the “linear” (small image contrast) and “nonlinear” (large image contrast) regimes. K2012 noted the frequent occurrence of high-definition structures in proton images and highlighted caustics of the proton optical system in the nonlinear regime as the key to interpreting such structures. Starting from a selection of simplified field configurations, K2012 explored the caustic image structures to be expected in radiographs of those configurations and advocated for an approach that relies on recognizing the typology of such structures and using them to infer the fields that produce them. This approach was made more systematic in the forward-modeling code of Levy et al.,21 who model laser-driven point-like proton sources with three-dimensional electromagnetic fields of arbitrary strength and structure, to create synthetic proton radiographs.

In this paper, we provide a first-principles description of the nature of the images of magnetized plasmas produced by proton imaging and discuss the implications for gaining information about the morphology and the strength of magnetic fields in high energy density laboratory experiments. Our concern is strictly with the imaging of magnetic fields; we do not treat the imaging of electric fields.

In this work, we mainly treat the linear regime (i.e., the small image contrast regime). We do venture into the nonlinear regime (order unity image contrast) to explore the validity of the linear method when applied to order unity image contrasts and furnish a path forward for formally addressing the nonlinear regime. We expand on the latter in a subsequent paper.22 We show that in the linear regime, image contrasts in proton radiographs are images of projected MHD current. It follows from this observation that it is perfectly possible for high-definition structures to occur even in the linear regime, where currents with sharp distributions occur, sharp features must occur in their radiographic images, irrespective of whether focusing is strong or weak. We follow up the brief discussion in K2012, which showed that inversion is, in principle, possible in the linear regime by furnishing and verifying a practical field reconstruction algorithm and, separately, an algorithm for recovering turbulent spectra.

In Sec. II we describe the experimental setup, including the proton imaging diagnostic, that we consider in this paper. In Sec. III we provide a theory of proton radiography in the limit of small image contrasts, a situation that applies in HEDP experiments. In Sec. IV we show how the proton radiographic image can be used to reconstruct the average transverse magnetic field using this theory. We present verification tests of the theory in Sec. V and investigate the validity of the method with increasing image contrast. In Sec. VI, we explore and characterize some of the errors that affect the reconstruction of the magnetic field using proton radiographs. In Sec. VII we apply the theory to reconstruction of the spectrum of the turbulent magnetic field produced by a turbulent flow. In Sec. VIII, we discuss an approach that can be used to extend the field reconstruction algorithm to the nonlinear regime. In Sec. IX we discuss potential limitations of caustic analysis, which is often employed to decipher proton radiographs. We discuss our results and draw some conclusions in Sec. X.

**II. PROTON RADIOGRAPHY: DEFINITIONS AND MAGNITUDES**

In this work, we concern ourselves exclusively with proton-radiographic imaging of magnetic structures, neglecting all electric forces on proton beams, unlike the more general investigation of K2012.20

An illustration of the experimental setup for a proton radiographic image system is shown in Fig. 1. Protons emitted isotropically from the implosion source at the left of the figure pass through an interaction region of the non-zero magnetic field, where they suffer small deflections due to the Lorentz force. They then continue on to strike a screen, where their positions are recorded.

If the transverse size of the interaction region di is assumed small compared to the distance ri of the interaction region from the implosion source, then the proton trajectories are nearly parallel. That is, if we assume that the angular spread ω ≡ di/ri ≪

1—the “paraxial” assumption also made in K2012—we may use planar geometry at the screen. We will adopt this assumption in what follows.

We further assume that the longitudinal size of the interaction region, li, is small compared with ri so that the parameter λ ≡ li/ri ≪

1. This assumption is convenient because it allows certain geometric factors arising in integrals to be treated as constants. This assumption was considered equivalent to the previous assumption ω ≪

1 in K2012, since in that work the longitudinal and transverse extents of the interaction region were assumed equal. Since there are two separate uses made here of the “small interaction region” approximation, one (paraxiality) related to the transverse extent and the other related to the longitudinal extent, we prefer to distinguish between the two extensions, even if they are usually physically similar.

A third simplifying assumption made here in common with K2012 is that the angular deflections α induced by the magnetic field are very small. We may quantify this conservatively by assuming a uniform field strength B, which gives rise to a Larmor radius rB=mvceB

for protons of mass m moving at velocity v. When rB ≫li, we may estimate α≈lirB=eBliβmc2, where β≡vc. If the protons have kinetic energy T, then β≈2T/mc2‾‾‾‾‾‾‾√. Therefore,

α=e2mc2√T−1/2Bli=1.80×10−2rad×(T14.7MeV)−1/2×(B105G)×(li0.1cm).

(1)

With proton kinetic energies of 3 MeV or 14.7 MeV characteristic of D3He capsules and magnetic field strengths of even 106 G, even a coherent field produces small angular deflections. A random walk due to a disordered field would produce even smaller deflections.

K2012 notes the importance of a further “small” experimental parameter, which measures the strength of the variations of angular deflection induced by varying fields across the interaction region. Assuming a perturbation of length scale a in the transverse direction inducing a deflection angle of characteristic size α, the definition of the dimensionless number μ is

μ≡liαa.

(2)

In effect, μ measures the amount of deflection per unit (transverse) length in the interaction region. It is clearly closely related to the derivatives of the magnetic field. The particular importance of μ is that it controls the amplitude of density contrast on the radiographic image: small values of μ ≪

1 correspond to small-amplitude contrasts; larger values of μ ≲

1 result in order unity (i.e., nonlinear) contrasts; whereas μ > 1 can lead to a strongly nonlinear regime where proton trajectories cross and the proton focusing may lead to the formation of caustics. Thus, the parameter μ can be used to demarcate different contrast regimes, which we further explore in a separate paper.22

It is noted in K2012 that μ and α are independent parameters, which may be large or small without reference to each other. In particular, one may have fields that are intense (large α) but relatively uniform (small μ) yielding small contrasts or conversely weak fields (small α) with large gradients (large μ) producing large contrasts. The latter regime was of interest in K2012, since it is the regime in which caustic surfaces may occur.

To summarize the research reported in this paper with respect to the parameters ω, λ, α, and μ, we assume that α ≪

1, ω ≪ 1, λ ≪ 1, and μ ≪ 1, which is commonly the case in many HEDP-related proton-radiographic experimental setups. We will see later that in this “linear” (in μ) regime, μ is in fact directly related to MHD current density ∇

× B so that in effect contrast amplitude is directly related (in fact, proportional) to current strength.

As we will see, the identification of μ with MHD current is sufficient to reconstruct the average transverse magnetic field in detail and, separately, to infer the spectra of isotropic-homogeneous magnetic turbulent fields. It also serves as a caution with respect to the use of the caustic structure in the interpretation of radiographs of magnetized plasmas: high-definition structure on a radiograph cannot, in general, be attributed to caustic surfaces of the proton optical system, since they can be produced by sharp current distributions in regimes where magnetic field gradients are too weak to produce caustics. Misidentification of high-definition radiographic “blobs” with caustics can therefore lead to overestimation of magnetic field strengths.

**III. SMALL FLUCTUATION THEORY OF PROTON RADIOGRAPHY**

We now consider the small contrast regime, μ ≪

1. As we will now see, the theory of such contrasts can be constructed in terms of a 2-dimensional deflection vector field in the focal plane. The divergence of this field yields the proton fluence contrast field directly and may be related to derivatives of the magnetic field in the interaction region—specifically to the MHD current.

**A. The lateral displacement field**

The fluence distribution of protons at the screen is represented by an area density function, Ψ(x⊥

) d2x⊥, where x⊥ represents position along the screen. In the absence of a magnetic field, the density function is uniform by assumption, Ψ0(x⊥) = ψ0, a constant. The effect of magnetization in the interaction region is to introduce a field of small lateral displacements w⊥(x⊥) that distort the uniform field Ψ0 to produce Ψ(x⊥). The w⊥

are “small” in the sense that they correspond to the small angular deflections as discussed earlier.

We may relate the change in Ψ due to the magnetic field to the displacement field w⊥

. If we let the undisturbed coordinates be x(0)⊥, so that x⊥=x(0)⊥+w⊥, we have by conservation of protons,

Ψ(x⊥)d2x⊥==≈Ψ0(x(0)⊥)d2x(0)⊥Ψ0(x⊥−w⊥)∣∣∣∂x(0)⊥∂x⊥∣∣∣d2x⊥[Ψ0(x⊥)−w⊥⋅∇⊥Ψ0(x⊥)+ (α2)]×[1−∇⊥⋅w⊥+ (μ2)]d2x⊥,

(3)

to the first order in α and μ. It follows that

Ψ(x⊥)≈Ψ0(x⊥)−∇⊥⋅[Ψ0(x⊥)w⊥],

(4)

which expresses the conservation of the “flux” Ψ0(x⊥)w⊥ in the plane of the detector screen. Since the unperturbed proton flux is uniform by assumption, so that Ψ0(x⊥) = ψ0, a constant, we finally obtain

δΨ(x⊥)/ψ0≡≈Ψ(x⊥)−ψ0ψ0−∇⊥⋅w⊥.

(5)

Equation (5) specifies our chief observable, the contrast field

Λ(x⊥)≡δΨ(x⊥)/ψ0,

(6)

and relates it to the model quantity w⊥

.

It should be clear from the manipulations above that we are making essential use of the conditions α ≪

1 and μ ≪ 1. In particular, |w⊥|/rs∼ (α) and Λ∼ (μ)

, and we drop all terms of higher order in α and μ.

**B. The Lorentz force**

We may relate the lateral displacement field w⊥

to the magnetic field in the interaction region. Between its emission and a later time t, a proton experiences a lateral motion δx⊥(t) given by

δx⊥(t)=∫t0dt1v⊥(t1).

(7)

Here, v⊥(t1) is the lateral velocity experienced at time t1 due to the Lorentz force accelerations received at prior times. It is given by

v⊥(t1)=∫t10dt2emcvn×B(x(t2)),

(8)

where n is the tangent vector to the proton trajectory, which is approximately constant along the trajectory due to the smallness of the angular deflection α.

We treat the implosion source as approximately point-like and use the near-constancy of n so that the argument x(t) of B in the integrand of Eq. (8) is given by x(t) = nv

t. The Lorentz force is perpendicular to n, so the longitudinal velocity v is constant and t = z/v, where z denotes distance along the n-direction. We change variables from t to z, obtaining

δx⊥(z)=1v∫z0dz1v⊥(z1),

(9)

v⊥(z1)=emc∫z10dz2n×B(nz2).

(10)

Combining Eqs. (9) and (10), we get

δx⊥(z)=emcvn×∫z0dz1∫z10dz2B(nz2).

(11)

We may switch the order of integration, to obtain

δx⊥(z)==emcvn×∫z0dz2B(nz2)∫zz2dz1emcvn×∫z0dz2B(nz2)(z−z2).

(12)

The intuitive meaning of this equation is that the lateral displacement δx⊥(z) of a proton at z is a superposition of displacements resulting from additive velocity pulses emc(nv)×Bd(z2/v), each pulse operating for a time (z − z2)/v

.

The relation between n and the transverse location at the screen x⊥

is

n=[zˆ+x⊥/rs]+ (ω2),

(13)

where, again, ω = di/ri is the small spread angle of protons crossing the interaction region. Defining

r(z,x⊥)≡(zˆ+x⊥/rs)z,

(14)

we may express the lateral deflection at the screen, w⊥ = δx⊥(rs) to leading order in ω and α,

w⊥(x⊥)=emcvn×∫rs0dzB(r(z,x⊥))(rs−z).

(15)

As a further simplification, we may take advantage of the fact that in the integrand of Eq. (15), the variation of the magnetic field is much more rapid than that of the factor (rs − z), which may therefore be taken outside of the integral and replaced with (rs − ri) at the cost of a relative error of order ∣∣z−ri∣∣/ri∼li/ri=λ

(this is the “small longitudinal extent” approximation discussed earlier). The result is

w⊥(x⊥)=e(rs−ri)mcvn×∫rs0dzB(r(z,x⊥)).

(16)

**C. Divergence of deviation field**

From Eq. (5), we know that the proton density field δΨ(x⊥

)/ψ0 is connected to the quantity ∇⊥· w⊥, rather that to w⊥ directly. To calculate this quantity, it is convenient to extend w⊥(x⊥) to a vector field u(x) defined in a 3-D neighborhood of the detector screen. The definition of u is

u(x)≡e(rs−ri)mcvx∣∣x∣∣×∫rs0dzB(x∣∣x∣∣z).

(17)

When x=x⊥+rszˆ, it follows that ∣∣u(x)−w⊥(x⊥)∣∣/∣∣w(x⊥)∣∣∼ (ω). We also have ∇⋅u=(∇⊥⋅w⊥)(1+ (ω)), since the component of u along zˆ is (ω)

. This is convenient because it is easier to take unconstrained spatial derivatives of the expression in Eq. (17) than to take “perpendicular” derivatives of the expression in Eq. (16).

We therefore have

Λ(x⊥)==−∇⊥⋅w⊥(x⊥)=−∇⋅u(x)e(rs−ri)mcvx∣∣x∣∣⋅∇×∫rs0dzB(x∣∣x∣∣z),

(18)

where, to get the third line, we have used the identity ∇·(P × Q) = Q·∇×P − P·∇×Q and the fact that ∇×x∣∣x∣∣=0

.

We pass over to tensor notation, denoting a vector q by its components qi, i = 1, 2, 3, and using the Einstein convention that repeated indices imply summation. Letting r(z)≡x∣∣x∣∣z

and taking the derivative inside the integral, we find

Λ(x⊥)=====≈e(rs−ri)mcvxi∣∣x∣∣ϵijk∫rs0dz∂∂xjBk(r)e(rs−ri)mcvϵijkxi∣∣x∣∣∫rs0dz∂∂rlBk(r)∂rl∂xje(rs−ri)mcvϵijkxi∣∣x∣∣∫rs0dz∂∂rlBk(r)z(δjl∣∣x∣∣−xjxl∣∣x∣∣3)e(rs−ri)mcvϵijkxi∣∣x∣∣2∫rs0dzz∂∂rjBk(r)e(rs−ri)mcv1∣∣x∣∣n⋅∫rs0dzz∇×B(r(z))eri(rs−ri)mcvrsn⋅∫rs0dz∇×B(r(z))≡χ(x(0)⊥).

(19)

Here, in the first line, we have used the totally antisymmetric Levi-Civita tensor ϵijk to express the cross product in the tensor form; to get from the third to the fourth equality, we have used the fact that ϵijkxixj = 0 to eliminate the second term in the parentheses of the integrand; and in the final line, we have appealed to the rapid variation of B(r(z)) in the interaction region to replace the factor z in the integrand by ri(1+ (λ)) and take it outside the integral, while also setting |x|=rs(1+ (ω)). In the final line, we have defined the current projection function χ(x(0)⊥), a dimensionless function of the undeflected coordinates x(0)⊥

.

Equation (19) is amenable to an interesting interpretation. The MHD current J is given by c4π∇×B

so that the contrast map Λ(x⊥

) is in fact proportional to the line integral of J along the flight path of the protons. That is, the image on the screen is basically the projection of the z-component of the MHD current field. We can now also interpret the small contrast regime as the regime of small MHD currents.

It follows that in small-contrast theory, proton radiographs of MHD plasmas are, for all intents and purposes, projective images of the MHD current component along the proton beam. This is a very useful observation, as it furnishes an attractive alternative to the forward modeling advocated in K2012, where substructures of field topologies are assembled to interpret the radiographs.

**IV. AVERAGE TRANSVERSE MAGNETIC FIELD RECONSTRUCTION**

An interesting and useful further consequence of the theory developed thus far is that it is possible to reconstruct the proton deflections w⊥

due to a magnetized plasma. By Eq. (16), this is equivalent to reconstructing the z-integrated transverse magnetic field. This possibility was recognized but not exploited in K2012, and the treatment in that work did not recognize an important issue that must be addressed in order for such a reconstruction to be practical, as we now discuss.

We begin by writing Eq. (5) as follows:

−∇0⊥⋅w(x(0)⊥)=Λ(x⊥),

(20)

where w(x(0)⊥) is the lateral deflection as a function of undeflected lateral coordinate x(0)⊥ and the notation ∇0⊥ refers to differentiation with respect to the undeflected coordinates x(0)⊥. The function w(x(0)⊥)

is written in this way to emphasize that it is a function of the unperturbed coordinates (of which it is also a perturbation).

The point about the coordinate dependence of w is worth belaboring because it is the origin of a subtlety in Eq. (20): the source term on the right-hand side is a function of the perturbed coordinates x⊥=x(0)⊥+w(x(0)⊥)

. We therefore have

−∇0⊥⋅w(x(0)⊥)=Λ(x(0)⊥+w(x(0)⊥)).

(21)

The subtlety worth emphasizing here is that we may not, in general, neglect the correction to the argument of Λ on the RHS of Eq. (21), despite the small-deflection approximation α ≪

1. The reason is that while the angular deflections may be a small fraction of a radian, they may still result in a substantial change in Λ, in regions where gradients of Λ are large.

Now, as pointed out in Eq. (78) of K2012, the deflection field w(x(0)⊥)

is irrotational, as a consequence of the solenoidal character of B. That is, there exists a scalar function ϕ(x(0)⊥) such that

w(x(0)⊥)=−∇0⊥ϕ(x(0)⊥).

(22)

Combining Eqs. (21) and (22), we obtain the field reconstruction equation

∇2ϕ(y)=Λ(y−∇ϕ(y)),

(23)

where we have unburdened the expressions of their superscripts and subscripts to reduce notational complexity by introducing y≡x(0)⊥

and stipulating that all gradients and Laplacians are two-dimensional operators in the plane of the detector and differentiate with respect to y.

Equation (23) differs from the reconstruction equation of K2012 [Eq. (20) of K2012, basically the Poisson equation for the projected potential, and its extension to the magnetic case given in Eq. (79)] by the argument of the source term on the RHS. We have found that, as a practical matter, the correction to the argument may not be neglected; doing so results in a very inaccurate algorithm because Λ can change substantially on this scale. The correction may be treated approximately, however, as we now demonstrate.

Equation (23) is a second-order, elliptical, nonlinear partial differential equation for the scalar function ϕ, which must be solved numerically. It is not straightforwardly solvable in its full nonlinear form. However, we may expand the right-hand side to first order in ∇

ϕ, obtaining

∇2ϕ(y)=Λ(y)−∇Λ(y)⋅∇ϕ(y).

(24)

Multiplying Eq. (24) through by the integrating factor exp(Λ(y)), we obtain

∇⋅(eΛ(y)∇ϕ(y))=Λ(y)eΛ(y).

(25)

Equation (25) is a steady-state diffusion equation, with an inhomogeneous diffusion coefficient exp(Λ) and a source term Λ exp(Λ). It may be solved by standard numerical methods, as we demonstrate in Sec. V.

Once a solution has been obtained, it is straightforward to recover the deflection field w = −∇

ϕ. The projection integral of the perpendicular magnetic field may then be obtained using Eq. (16),

∫rs0dzB⊥(r(z,y))=mcve(rs−ri)∇ϕ(y)×z,

(26)

where we have made our now usual approximation n(0) ≈ z at the cost of a small error (ω). The vector product with z in effect rotates the vector ∇

ϕ clockwise by 90°.

We may use Eq. (26) to estimate the longitudinally averaged transverse field,

B⎯⎯⎯⊥≡1li∫rs0dzB⊥(r(z,y)),

(27)

which is obviously interesting information that is well-worth extracting. In order to obtain this average, we clearly need some way of estimating the longitudinal extent of the interaction region, li. Such an estimate may be available from a separate diagnostic measurement of the plasma or from a simulation. Alternatively, one may assume a rough isotropy, implying that the transverse extent of the interaction region is approximately the same as the longitudinal extent. Under this assumption, one may use the radiographic data to estimate the transverse extent and parlay this into the required estimate of li.

**V. NUMERICAL VERIFICATION OF THE THEORY**

In order to verify the theory developed thus far, we have written a small ray-tracing code in Python and used it to simulate large ensembles of 14.7 MeV protons that are fired at a magnetized region from a central source and recorded at a screen. The code uses the ordinary differential equation (ODE) integration routines from the Scipy23 library to integrate the trajectory equations,

dnds=(αli)n×b,

(28)

and

dxds=n,

(29)

as an ODE system. Here, s denotes the physical length along each integral curve of n, and we measure the magnetic field in units of some typical field strength Bt in the domain so that b ≡ B/Bt and the deflection angle α is the “typical” deflection α=eBtlimvc

.

The ODE system is further augmented by the tracer equations

dχds=eri(rs−ri)mcvrsn⋅∇×B,

(30)

and

dds(liB⎯⎯⎯⊥)=B⊥,

(31)

which allow us to trace the current projection function χ and the average perpendicular field B⎯⎯⎯⊥

experienced by each proton, so that we can map them to the screen and compare them with the fluence contrast Λ and the reconstructed perpendicular field, respectively, thus verifying the theory.

Protons are fired in a cone about the z axis, bounded by a circular target in the x–y-plane at the location of the interaction region. The cone is distorted by the Lorentz force, and the protons’ locations and tracer data are recorded at z = rs. The cone is then “collimated” to a square region, with boundaries aligned with the x- and y-axes, such that the square is entirely contained by the circular aperture of the cone.

We simulate point sources only, located at a distance ri = 10 cm from the interaction region and rs = 100 cm from the screen. The magnetic fields that we simulate have characteristic field strengths of 2 × 104 G. The field configuration that we simulate is a magnetic “ellipsoidal blob,” described in K2012. The field is purely azimuthal about an arbitrary axis, with field strength

Bϕ(r′,z′)=B0r′aexp(−r′2a2−z′2b2),

(32)

where r′ is the perpendicular distance from the azimuthal symmetry axis and z′ is the distance along the axis. The axis may be rotated arbitrarily with respect to the arrival direction of the protons.

All images are binned into 128 × 128 pixels. The protons are summed in each bin, while the current projection function χ and the average transverse field liB⎯⎯⎯⊥

, integrated for each proton according to Eqs. (30) and (31), are averaged in each bin.

The pixelated proton counts are converted to fluence contrast Λ using Eq. (6). The field reconstruction is carried out by solving Eq. (25) for ϕ numerically on the 128 × 128 image grid, by the slow but effective method of the unaccelerated Gauss-Seidel iteration (see Sec. 19.5 in the study of Press24), which gives serviceable results for present purposes. We assume the Dirichlet boundary conditions, that is, ϕ(y) = 0 at the image boundary. This scheme, which in effect solves a linear problem of the form Ax = b, iteratively reduces the residual |Ax − b|, where || refers to the 2-norm. We monitor the normalized L2-norm of the residual, which is |Ax − b|/|b|, for convergence.

The initial guess is a solution of the Poisson equation obtained from Eq. (23) by setting the deflection ∇

ϕ = 0 in the argument of Λ on the RHS. This solution can be obtained very rapidly by fast Fourier transform (FFT) techniques but yields a poor match to the true field, and we therefore only use it to initialize the iterative scheme.

The deflections, and the field projection, are then recovered using Eqs. (22) and (26).

The reconstruction code, named PRaLine (Proton Radiography Linear reconstruction), is written in Python and uses routines from the Scipy23 library. We have made the PRaLine code publicly available through GitHub (Publicly available at https://github.com/flash-center/PRaLine) and the FLASH code release (publicly available at https://flash.uchicago.edu). PRaLine is a simple, easy-to-use, and efficient code for analyzing synthetic and experimental proton radiographs.

**A. An isolated blob**

We begin with a simple, clean-room case: an ellipsoidal blob, whose image is well-separated from the boundaries of the image. The field strength B0 = 2 × 104 G, and the azimuthal symmetry axis of the blob is aligned with the z axis. The blob has geometric parameters a = 0.03 cm, b = 0.02 cm and is centered on the z axis (that is, down the middle of the beam).

In this test, we fire 10 × 106 protons at the circular target, and after collimation to the square image we have a yield of 6.1 × 106 protons.

In the top row of Fig. 2, we show the comparison of the nonlinear fluence contrast Λ computed from the pixelated data (left panel) to the current projection function χ, integrated along proton trajectories (middle panel). The agreement is visually very good. To perform a more quantitative comparison, in the right hand panel, we show the normalized difference (Λ − χ)/σΛ, where σΛ is the per-pixel Gaussian error. This is a reasonable procedure, since the mean counts per pixel in the simulation is 373 so that the pixel count distribution is well-described by Gaussian statistics. The error σΛ for a pixel enclosing n protons is computed starting from the Poisson error for the counts, σP=n‾√

, and using the standard propagation of errors formula the functional relationship between Λ and n. In terms of the average counts per bin, n0, we have

σΛ=n‾√n0.

(33)

We can see from the top-right panel of Fig. 2 that the current projection function predicts the fluence contrast perfectly. The fluctuations shown in the figure may further be summed in quadrature to supply a χ2 value of 16 585.1 for 16 384 degrees of freedom (DOF). This corresponds to a P-value of 13.4%, indicating a perfectly acceptable (parameter-free) “fit” of the model to the data.

In the lower panels of Fig. 2, we show analogous results for the same simulations, using raw counts and the predicted counts per pixel. Once again, the comparison is visually and quantitatively excellent.

To explore the range of validity of the linear theory, we repeat the exercise and gradually increase the field strength to reach fluence contrasts with peak values of order unity and larger (i.e., in the nonlinear regime). The results are reported in Fig. 3, where we compare the nonlinear fluence contrasts Λ [panels (a), (d), and (g)] to the corresponding current projection functions χ[panels (b), (e), and (h)]. The first row corresponds to χpeak = 1, the second row corresponds to χpeak = 1.4, and the third row corresponds to χpeak = 1.7. The agreement between panels [(a) and (b)] and [(d) and (e)] is—visually—fairly good, despite the large peak field-strength at the blob’s center, which brings the fluence contrast into the nonlinear regime locally (see also the first four columns of Table I). Conversely, we see deviations from agreement on the first two panels of the third row, where Λ and χ values no longer match. This is well-captured in the normalized difference (Λ − χ)/σΛ [panels (c), (f), and (i)], where we note a gradual increase in the error at the center of the domain.

The reconstruction of the field configuration for the nominal case (Fig. 2) is displayed in the upper panels of Fig. 4. We ran the iterative solver for 4000 iterations, at which point the residual error was 0.51%. The top-left panel shows the actual field projection liB⎯⎯⎯⊥ from the integration of Eq. (31), while the top-right panel shows the reconstructed field. The colormap displays field strength on a logarithmic scale. We see that the field reconstruction in the center of the image, near the peak field-strength location, is very good, while a good deal of unstructured noise appears in portions of the image where the field strength is more than three orders of magnitude less than the peak value. At the location of peak field strength, the difference between actual and reconstructed field strengths is 2.0%. A quantitative comparison of the two images may be performed using the L2 norm,

L2≡⎡⎣⎢⎢∑p∈pixels∣∣B(Reconstructed)p−B(True)p∣∣2∑p∈pixels∣∣B(True)p∣∣2⎤⎦⎥⎥1/2. (34)

For this reconstruction, we find L2 = 0.13. The discrepancy is presumably a composite of the Poisson noise (itself around 5% per pixel), discretization error, and algorithmic error (the algorithm terminated with a 0.53% residual).

The lower panels in Fig. 4 show the concomitant proton deflections, as reconstructed by the algorithm. These figures illustrate the expected focusing effect of the azimuthal field. Again, the vector displacements maintain good accuracy in the range from peak displacement length to displacements that are about three orders of magnitude down from peak, after which unstructured noise appears.

To visualize the accuracy of the magnetic field reconstruction, we compare lineouts of the actual field projection liB⎯⎯⎯⊥

with the reconstructed field, for gradually increasing field strengths. The lineouts shown in Fig. 5 are taken along the x axis at Y = 0 cm and correspond to peak values of the current projection function χpeak = 0.7 [panel (a)], χpeak = 1 [panel (b)], χpeak = 1.4 [panel (c)], and χpeak = 1.7 [panel (d)]. The discrepancy between “true” and reconstructed fields is again quantified using the L2 norm and is reported on Table I. This exercise reveals that the linear reconstruction algorithm performs adequately well even when the radiograph has localized nonlinear contrasts [e.g., Figs. 2(a), 3(a), and 3(d)], provided that the root mean square (RMS) contrast fluctuations remain small (see second column of Table I). This seems to indicate that the RMS contrast fluctuations are a better indicator of the overall “linearity” of a given radiograph. On the other hand, when non-linearity affects a substantial part of the image, as is the case for χpeak = 1.7 [Fig. 3(g)], the reconstruction fails to correctly capture the actual field projection [Fig. 5(d)] and the L2 norm becomes large (L2 = 0.5 for χpeak = 1.7).

**B. Two superposed blobs**

We also simulate a more complex and less symmetric situation, by superposing two blobs of unequal strength—1 at y = +0.025 cm, with B0 = 2 × 104 G, and 1 at y = −0.025 cm, with B0 = 104 G. The geometric parameters of the blobs are both identical to those in the previous example—a = 0.03 cm, b = 0.02 cm. We again fire 10 × 106 protons at the system, recovering 6.1 × 106 protons after collimation to the square image.

Comparisons of fluence contrast to current projection, and of counts per pixel to predicted counts per pixel, are displayed in Fig. 6. The comparison is again both visually and quantitatively satisfactory. The difference of the fluence contrast and current projection maps yields a χ2 = 16 448.4 for 16 384 DOF, a P-value of 36%. The model is quite clearly in good agreement with the data.

The reconstruction is illustrated in Fig. 7. The iterative solution residual after 4000 iterations is 0.60%. Again, the visual impression of the reconstructed field is quite good, with good structural fidelity within three orders of magnitude peak, for both magnetic field and for deflections. At the peak field strength location, the difference between the reconstructed and actual field strengths is 2.1%. The L2 norm of the difference between the reconstructed and actual fields is L2 = 0.13.

**VI. EXPLORATION AND CHARACTERIZATION OF RECONSTRUCTION ERRORS**

In this section, we explore the character of some errors that affect the analysis of radiographs. We are concerned here not with the general subject of systematic errors in proton radiography (a subject that received a very thorough treatment in K2012) but rather with errors that limit the accuracy of the magnetic field reconstruction.

We can identify several individual (although not necessarily independent) sources of error that are relevant to reconstruction:

Poisson noise: In the diffusion equation that we use for reconstruction, Eq. (25), both the source term on the right-hand side and the diffusion coefficient are estimated noisily from binned proton counts. This Poisson noise is a source of error whose principal effect is uncertainty that scales as n−1/2, where n is the average number of protons per pixel. However, we shall see that the Poisson noise can also be expected to produce other indirect effects.

Discretization noise: The numerical solution is carried out on a uniform grid in the image plane. Obviously, the level of refinement of the grid has an impact on the accuracy of the solution. In principle, the mesh should be fine enough to resolve the smallest magnetic structures in the image. However, one may not simply refine the mesh indefinitely, as the reduced error due to discretization will eventually be balanced by the Poisson noise as the average number of protons per pixel drops and the form of the fluence contrast function Λ becomes noisily determined. Thus, there is a trade-off between the discretization noise and Poisson noise that should be optimized somehow.

Edge effects: The solution of an elliptical partial differential equation (PDE) such as Eq. (25) requires assumptions about boundary conditions. A natural choice is Dirichlet boundary conditions (BC), with the function ϕ(y) → 0 at the edges of the radiograph. If the active region is well-separated from the edges, then this choice may be expected to give good results. If the active region is not well-separated from the edges, this choice is not ideal and may at best be regarded as an approximation of convenience. In this case, it seems clear that the ability to recover detailed field structure is compromised. It is interesting to inquire whether one may still infer the correct order of magnitude of |B⊥

| in these circumstances. It should also be noted that even for an active region that is isolated from the boundaries, the Poisson noise can furnish spurious activity near the boundary. Thus, edge effects and Poisson noise effects are also intertwined.

Obstruction effects: It is not uncommon in laser experiments for some part of the physical structure of the target to be illuminated by the beam of protons and thus to cast a shadow on the image. Such a shadow, treated naively, would appear to the present analysis as an enormous, proton-clearing field of large strength and strong gradients. Evidently it is necessary to characterize the importance of this kind of error. One very crude fix is to replace the proton fluence in the shadowed regions by the mean, zero-deflection proton fluence. If the shadow encroaches on an image region of magnetic activity, this procedure will clearly spoil the reconstruction of the detailed field structure. We may still hope to at least recover the typical field strengths, at least approximately.

Later, we explore numerically some of the characteristics of the Poisson noise, discretization noise, edge effects, and obstruction effects.

**A. Poisson noise effects**

We explore the influence of the Poisson noise using the first field simulated in Sec. V, which was displayed in Fig. 4. We launch 108 protons at this field and produce datasets of different fluences by truncating this simulation to 6 × 105, 106, 3 × 106, 6 × 106, 107, 3 × 107, 6 × 107, and 108 protons. In each case, after collimation to the square field-of-view (FOV), about 60% of the protons remain. We perform the reconstruction procedure on each dataset, using the same 128 × 128 grid that we used in Sec. V.

Three sample reconstructions are displayed in Fig. 8. They correspond to the datasets with 6 × 105 protons (3.7 × 105 in the FOV), 6 × 107 protons (3.7 × 107 in the FOV), and 108 protons (6.1 × 107 in the FOV). The progressive degradation of the reconstruction is evident and proceeds from the regions of the weaker field (where the current integral is more poorly determined) inwards to the regions of the stronger field.

The lower-right panel shows the L2-norm of the difference between the reconstructed field and the actual field, normalized to the L2-norm of the actual field, as a function of number of protons in the FOV. The error shows the expected Poisson scaling of n−1/2 for n ≲

107. Evidently, for these simulations, other types of error begin to dominate the error budget above this proton fluence.

**B. Discretization effects**

As discussed earlier, the effect of grid spacing on reconstruction accuracy is inextricably entwined with the effect of the Poisson noise. For a fixed radiograph, the resolution can only be increased up to the point where the per-pixel proton counts begin to be so low that the Poisson noise begins to degrade the reconstruction.

To explore the connection between the resolution and Poisson noise, we conducted three series of reconstructions. Each series comprised a set of six grids—32 × 32, 45 × 45, 64 × 64, 90 × 90, 128 × 128, and 181 × 181 grid points, respectively (each contains twice as many points as the previous one). We performed this series of reconstructions on each of three simulations of a blob. The three simulations differed only in the number of protons—they contained 6.1 × 105, 2.4 × 106, and 6.1 × 106 protons in the square FOV, respectively. The selected blob had parameters that differed from the one in Sec. V A, in the single respect that a = 0.05 cm instead of a = 0.03 cm. We expanded the blob laterally so as to give the reconstructions a better chance to capture the field at the low resolution end.

The results are displayed in Fig. 9. The blue circles show the relative L2 norm of the difference between the reconstructed field and the “true” field. The red crosses show the relative L2 norm of the difference between two “adjacent” resolutions (the crosses are therefore located between their corresponding pairs of circles). The value of the latter measure is that it is available in experimental situations, when the true field is not known.

The figure shows the very strong effect of proton fluence on reconstruction accuracy. Depending on the proton fluence, the reconstruction accuracy may be compromised even at low resolutions, or alternatively it may continue to improve up to the limit of our patience with the reconstruction algorithm (which has a time-to-solution that grows at least as fast as the number of grid points). The difference between adjacent resolutions appears to furnish a reasonable rough indicator of optimality, in the sense that in the low-fluence case (left panel) the crosses begin to rise at about the time when the departure from the true field begins to be significant.

**C. Edge effects**

To illustrate the influence of boundary effects on field reconstruction, we performed a simulation using a blob with the same parameters as the ones used in Sec. V A but displaced so that a substantial part of the previously reconstructed field now overlaps the right edge of the image. To do this, we displace the interaction region to the right by 1.1li. We shoot 107 protons at the target as before and again recover about 6 × 106 protons in the collimated square image.

The resulting reconstruction is displayed in Fig. 10. The comparison with the single isolated blob of Sec. V A is instructive. In this case, the detailed reconstruction accuracy is degraded—the L2 norm of the difference between actual and reconstructed fields has increased to 29% (from 13%), and the difference between actual and reconstructed field strengths at the location of peak field strength has increased to 13.5% (from 2.0%). On the other hand, the structure of the field is still recognizable: the reconstructed contours are still clearly concentric circles to a radius of 0.5 cm from the center point of the field, which was also the case in Fig. 4. It appears, then, that when the image of the interaction region encroaches on the edge of the radiograph, one may still hope to infer field strength peak magnitudes and field structure to some useful degree, despite the inexactly obeyed boundary conditions.

**D. Obstruction effects**

We study the effect of obstruction on field reconstruction using a variation of the two superposed blob simulation described in Sec. V B, where we summarily remove the protons in a vertical strip near the blob and replace the removed fluences by the mean, zero-deflection fluence for the image. We do this for obstruction strips of two different widths.

The results are shown in Fig. 11. The left panel of this figure shows the result of reconstruction using unobstructed data and is therefore the same as the reconstruction in Fig. 7. The middle panel shows the reconstruction that is obtained after obstructing the data in the vertical strip 0.1 < x < 0.2, while the right panel shows the reconstruction obtained after obstructing data in the strip 0.1 < x < 0.5.

What emerges from this study is that the field reconstruction can be remarkably robust to obstruction effects. In both obstructed cases, the typical field intensities are largely unaffected by the obstruction, and the morphology is not wildly distorted. This is, we believe, due to the non-local nature of the reconstruction, which can partly “heal” the data corruption, interpolating to preserve the assumption of continuity and differentiability of the field that is built in to the reconstruction procedure. Naturally, one would expect the fidelity of the reconstruction to be more seriously compromised as the obstruction becomes more serious. It is also likely that the unrealistic simplicity of the two-blob test case is helpful in aiding reconstruction. Nonetheless, this test case offers some reason to place some limited reliance on field strengths and morphologies inferred through the reconstruction procedure, even in the case of obstructed data.

**VII. PROTON RADIOGRAPHY OF ISOTROPIC-HOMOGENEOUS TURBULENT FIELDS**

In experimental applications that result in turbulent flows, the direct reconstruction of the transverse magnetic field is not necessarily desirable. Of greater interest is the spectrum of the turbulent magnetic field. While the spectrum is in principle recoverable from the reconstructed field, such a procedure would unnecessarily entangle the inferred spectrum with some of the reconstruction errors as discussed earlier (i.e., edge and obstruction effects). As we will now demonstrate, under the assumptions that the field has an approximately uniform and isotropic spectrum, the spectrum itself is directly and robustly recoverable from the radiographic image without performing the field reconstruction.

Our starting point is Eq. (19), which we rewrite slightly as

Λ(x⊥)=eri(rs−ri)mcvrszˆ⋅∫ri+li/2ri−li/2dz∇×B(r(z,x⊥)),

(35)

where we accept to incur an error of order (ω) by replacing n(0) with zˆ in the inner product and one of order (λ) by separating the dependence on z from that on x⊥ in the argument of the integrand, setting

r(z,x⊥)≡zzˆ+rirsx⊥.

(36)

In Eq. (35), we have also changed the limits of integration to include only the region in which the magnetic field is appreciable. The reason for this change will be made apparent later.

Note also that in Eq. (35) we are neglecting the small shift in the argument of Λ that was so important to the reconstruction algorithm [compare Eq. (21)]. This is permissible at the cost of stipulating that the spectrum obtained thereby cannot be trusted at scales close to, or smaller than, the angular deflection scale αli. At longer scales, the spectrum should be unaffected by this approximation.

**A. The two-point correlation function and the spectrum**

We define the two-point correlation function for the radiographic image as follows:

η(∣∣x′⊥−x⊥∣∣)≡⟨Λ(x⊥)Λ(x′⊥)⟩,

(37)

where the expectation value ⟨⟩ is an ensemble average. The fact that η(∣∣x′⊥−x⊥∣∣) depends only on the distance between x′⊥ and x⊥

is a consequence of the assumed isotropy and homogeneity of the stochastic field B, as will be shortly evident.

Inserting Eq. (35) into Eq. (37), we get

⟨Λ(x⊥)Λ(x′⊥)⟩=e2r2i(rs−ri)2m2c2v2r2sϵlmnϵl′m′n′zˆlzˆl′×∫ri+li+1/2ri−li/2dz∫r1+li/2ri−li/2dz′∂∂rm∂∂r′m′⟨Bn(r(z,x⊥))Bn′(r(z′,x′⊥)⟩.

(38)

The correlation expression in the integrand defines an isotropic tensor

⟨Bn(x)Bn′(x′)⟩≡Qnn′(x−x′),

(39)

that is analogous to the classic incompressible velocity correlation tensor, since ∇·B = 0 is analogous to the incompressibility condition for velocity. Following the development in Secs. 6.2 and 8.1 of the study of Davidson,25 we write

Qnn′(a)=∫d3keik⋅aΦnn′(k),

(40)

where the isotropic tensor Φnn′(k) has the form

Φnn′(k)=EB(k)k2(δnn′−knkn′/k2).

(41)

The quantity EB(k) is the energy spectrum of the magnetic field and satisfies the relation

∫∞0dkEB(k)=18π⟨B2⟩,

(42)

which is obtainable by setting a = 0 and contracting the indices in Eq. (40).

Inserting Eqs. (39)–(41) into Eq. (38), we obtain

⟨Λ(x⊥)Λ(x′⊥)⟩=e2r2i(rs−ri)2m2c2v2r2sϵlmnϵl′m′n′zˆlzˆl′×∫ri+li+1/2ri−li/2dz∫z1+li/2ri−li/2dz′∫d3keik⋅(r(z,x⊥)−r(z′,x′⊥))×kmkm′EB(k)k2(δnn′−knkn′/k2).

(43)

At this point, it is possible to clarify why it was important to adopt limits of z-integration that are confined to the interaction region. The magnetic field is confined to the interaction region, and we were to Fourier-transform B instead of Q, the various complex phases of the transform of B would contain the information necessary to respect this confinement. However, the spectrum EB(k) that appears in the Fourier transform of Q knows nothing of those phases. It knows that there is an integral scale cutoff k0 ≳

2π/li, but it represents a model of uniform, isotropic turbulence in which all space is pervaded by a turbulent B-field with integral-scale cutoff k0. The restriction on the limits of integration in z in effect prevents that model from incorporating proton deflections from regions where physically the field is zero.

The term knkn′

in the integrand of Eq. (43) vanishes, since ϵlmnkmkn = 0. We may then use the decomposition k=k⊥+kzzˆ and the Levi-Civita tensor’s contraction identity

ϵlmnϵl′m′n′δnn′=δll′δmm′−δlm′δml′

(44)

to obtain

⟨Λ(x⊥)Λ(x′⊥)⟩=e2r2i(rs−ri)2m2c2v2r2s∫∞−∞dkz4sin2(kzli/2)k2z×∫d2k⊥ei(ri/rs)k⊥⋅(x⊥−x′⊥)EB([k2⊥+k2z]1/2)k2⊥k2⊥+k2z.

(45)

As a function of kz. The integrand term 4sin2(kzli/2)/k2z

is sharply peaked near kz = 0 and has a width 1/li. Compared to the k scales present in EB(k), it is well-approximated by a δ-function,

4sin2(kzli/2)k2z≈2πliδ(kz).

(46)

Using this in Eq. (45), we obtain

⟨Λ(x⊥)Λ(x′⊥)⟩==2πe2r2i(rs−ri)2lim2c2v2r2s×∫d2k⊥ei(ri/rs)k⊥⋅(x⊥−x′⊥)EB(k⊥)2πe2(rs−ri)2lim2c2v2×∫d2q⊥eiq⊥⋅(x⊥−x′⊥)EB(rsriq⊥),

(47)

where we have defined the screen-plane wave vector q⊥≡rirsk⊥. Since EB(k⊥) depends only on the magnitude of k⊥, we see that the autocorrelation indeed only depends on the magnitude ∣∣x′⊥−x⊥∣∣

, as previously asserted.

Note that we had not restricted the z integration range, rather leaving it as 0 < z < zs, we would have obtained a similar result but scaled up by rs/li [because instead of the RHS of Eq. (46), we would have found an expression that is approximately 2πrsδ(kz)]. Thus, the predicted contrast spectrum would have too large an amplitude by a factor rs/li, which would correspond to deflections suffered in consequence of a turbulent field pervading the entire region between the implosion capsule and the screen, rather than one confined to a small interaction region.

We Fourier-transform Eq. (47) with respect to both arguments, obtaining

⟨Λ˜(p⊥)Λ˜(p′⊥)∗⟩=2πe2(rs−ri)2lim2c2v2δ2(p⊥−p′⊥)EB(rsrip⊥),

(48)

where we have defined the Fourier transform of the contrast fluctuation map,

Λ˜(p⊥)≡(2π)−2∫d2x⊥e−ip⊥⋅x⊥Λ(x⊥).

(49)

The δ-function in Eq. (48) may be disposed of by noting that the contrast fluctuation map Λ(x⊥

) is available as a discretized array on a mesh of cells, rather than as a continuous function, and its Fourier transform is necessarily a Discrete Fourier Transform (DFT). If the mesh on the screen is an N × N evenly spaced array over a square of lateral dimension L, the DFT will produce a complex array of transform values Λˆ(n), where n are 2-dimensional mode vectors whose components are integers, n = [nx, ny]. n are related to the wave vectors p⊥ by p⊥=2πLn. Λˆ(n) are then related to Λ˜(p⊥) by averaging the latter over the cell in p⊥-space corresponding to the mode n. That is to say,

Λˆ(n)=(2πL)−2∫Cellnd2p⊥Λ˜(p⊥),

(50)

where the integral is over the p⊥-space square cell of side 2π/L centered at 2πLn

.

We therefore average the dependence of Eq. (49) on p⊥

, p′⊥ over cells n, n′, respectively. This process duly removes the δ-function, leaving the result

⟨Λˆ(n)Λˆ(n′)∗⟩=e2(rs−ri)2liL22πm2c2v2δnn′EB(rsri2πL∣∣∣n∣∣∣).

(51)

The final result for the turbulence spectrum is then

EB(rsri2πL∣∣∣n∣∣∣)=2πm2c2v2e2(rs−ri)2liL2⟨Λˆ(n)Λˆ(n)∗⟩.

(52)

In a practical computation, the ensemble average in Eq. (52) may be replaced by a rotational average over cells in n-space, binning |n| to suit one’s convenience.

It is a noteworthy feature of Eq. (52) that the factor rsrip⊥

in the argument of EB expresses the expected stretching of wavelengths by the divergence of the proton rays on their way from the interaction region to the screen.

Note also that according to Eq. (52), we again require an estimate of the longitudinal extent of the interaction region li in order to obtain the amplitude of the spectrum EB. The discussion of the availability of such an estimate of li was given at the end of Sec. IV.

**B. Numerical verification with a Gaussian field**

We verify Eq. (52) by means of a simulation in which protons are fired at a domain containing a Gaussian random magnetic field of a known spectrum. Such a field does not really represent a realistic magnetic turbulence distribution, of course, since turbulence is well-known to be non-Gaussian (see Ref. 25, pp. 99-103, for example). This is not an issue here, since we merely wish to confirm that we can successfully recover the two-point spectrum of a homogeneous-isotropic distribution, irrespective of what the other moments of the distribution may be, so we are free to choose a Gaussian random field on grounds of simplicity.

The optical setup used in this test is ri = 20 cm, rs = 100 cm, and li = 0.1 cm. Protons energies are monochromatic, at 14.7 MeV, and the protons are emitted from a point source at the origin.

The simulated spectrum is a truncated power-law with the following form:

EB(k)=⎧⎩⎨⎪⎪⟨B2⟩8πα+1kα+12−kα+11kα:0:k1<k<k2otherwise.

(53)

This form allows us to specify the normalization in terms of the mean magnetic energy B2/8π, according to Eq. (42). We parametrize the spectral cutoffs in terms of discrete modes n1 and n2, with k1=2πlin1, k2=2πlin2

. In our simulations, we chose 〈B2〉1/2 = 5 × 102 G, α = −3.5, n1 = 2, and n2 = 40. These parameters, together with the optical system parameters rs, ri above, result in expected contrast fluctuations at the screen of about 0.1.

We generate the field in a 128 × 128 × 128 mesh of cells in a cubic domain. We start out in the Fourier space, randomly drawing the (complex) components of the Fourier-transformed vector potential a(k) in each k-space cell independently from a zero-mean normal distribution. The variance of the distribution is chosen to yield the desired spectrum EB. The reality constraint a(k)* = a(−k) is enforced by drawing samples only in the hemisphere kx > 0 and copying the complex conjugate of the field in this hemisphere to the reflected point in the opposite hemisphere. The resulting vector field is then Fourier-transformed to the real domain to obtain A(x), and the discrete curl is taken to yield the final magnetic field, B = ∇

×A.

We simulate firing 5 × 106 protons through this setup, in a cone of radius 0.1 cm at r = ri. The radiographic results are displayed in Fig. 12, which shows a square region of data entirely contained within the image of the interaction region.

To obtain the spectrum, we proceed as indicated in the discussion following Eq. (52). We bin the region imaged in Fig. 12 into a 128 × 128 array of cells, computing the contrast map in each cell, and perform an FFT on the result to estimate Λˆ(n)

. We bin the radial wavenumber n = |n| into equally spaced logarithmic bins and average the contributions in each bin to perform a circular average. Then, using Eq. (52), we estimate EB.

The result of this computation is displayed in Fig. 13. The left panel of this figure shows the result of the analysis just described. The right panel shows the result of repeating the analysis on the current projection function χ, which, as we have seen, predicts the fluence contrast. We note that χ is not an experimental observable (as is the fluence contrast Λ) and is recovered via the tracer Eq. (30) by the protons as they traverse the domain. As such, the spectral analysis on χ is expected to always reproduce the correct spectrum. From Fig. 13(a), it is clear that at low frequency, the method returns the correct spectrum—both the low-frequency cutoff k1 and the spectral slope are correctly recovered. At high frequency, we see a leveling-off of the spectrum in the contrast data that has no counterpart in the spectrum inferred from χ [Fig. 13(b)]. The white-noise spectrum due to the Poisson statistics dominates at large k, producing the leveling-off of the magnetic spectrum.

This observation points to a practical necessity in implementing this algorithm: it is necessary to optimize the trade-off between spectral resolution and the Poisson noise, by choosing the size of the image binning mesh carefully. Too high resolution can lead to low counts-per-bin, and therefore excessive noise, while too low resolution leads to the inability to probe high frequencies. This trade-off must be resolved empirically in the analysis of radiographs, on a case-by-case basis. These considerations are analogous to the ones described in Sec. V B.

Another high-frequency effect that is to be expected in real data, but that we do not simulate here, is that the finite size of the implosion capsule must necessarily produce a high-frequency cutoff in the measured spectrum. The reason is that a finite capsule size s blurs the image of every point in the interaction region by an angular blurring s/ri and hence by a linear blurring length s(rs − ri)/ri in the image. This corresponds to a length scale in the interaction region of λs = s(rs − ri)/ri × ri/rs = s(rs − ri)/rs. This is the limit to which we may expect to resolve any structure in the spectrum, which should therefore be expected to contain a high-frequency cutoff at wavenumber ks = 2π/λs. Naturally, for this cutoff to be visible, it must occur at wavenumbers k below the advent of the white noise plateau as discussed earlier.

In order to gauge the regime of validity of this spectral analysis, we repeat the test problem three more times, scaling up the Gaussian random magnetic field by factors of 2, 4, and 6. It should be noted that we do not consider a new random realization of the field but, rather, a scaled version of the same realization. The resulting fluence contrast, current projection contrast, and normalized difference for the scaled field cases are shown in Fig. 14, equivalent to the first row of Fig. 12. In the new radiographs, the peak predicted values of χ increase from χpeak = 0.3 [Fig. 12(b)] to χpeak = 0.6 [Fig. 14(b)], χpeak = 1.1 [Fig. 14(e)], and χpeak = 1.7 [Fig. 14(h)], thus bringing the fluence contrast well into the nonlinear regime. The RMS fluctuations of the image contrast increase from 0.11 to 0.25, 0.40, and 0.77, respectively. As in the isolated blob problem, we see that χ becomes a progressively worse predictor of Λ and the normalized error values increase considerably, more so where the contrasts are large. Nonetheless, the current projection contrast still captures the general morphology and features of the fluence contrast.

For each case, we construct the magnetic energy spectra utilizing the algorithm described previously and Eq. (52), applied to either the fluence contrast or the current projection function. The spectra are displayed in Fig. 15. Panels (a) and (d) correspond to the χpeak = 0.6 case, panels (b) and (e) correspond to the χpeak = 1.1 case, and panels (c) and (f) correspond to the χpeak = 1.7 case. The top panels are the result of the spectral analysis performed on the fluence contrast, while the bottom panels on the current projection function. The fluence contrast spectral plots [Figs. 15(a)–15(c)] show that the fluctuation spectrum departs from the model spectrum starting at high wavenumbers, eventually peeling off entirely. Unlike the direct reconstruction results for the isolated blob [e.g., Fig. 5(d)], spectral reconstruction appears to overestimate the field strength. On the other hand, the current projection spectra [Figs. 15(d)–15(f)] continue to predict the model spectrum well. This is not surprising since they are constructed from the actual field variables. We conclude that the failure of the spectral analysis in the nonlinear regime is evidently a direct consequence of the failure of the current projection to predict the contrast fluctuations.

To recover the path integrated magnetic field, we solve once more Eq. (25) assuming, as in the isolated blob test problem, the Dirichlet boundary conditions at the image boundary. This choice is arguably expedient but poor, since the magnetic field now fills the 128 × 128 × 128 cubic computational domain and therefore the fluence contrast is not negligible at the edges of the radiograph [Fig. 12(a)]. In Sec. VI C we saw how this can have deleterious effects in the accuracy of the reconstruction. Artificial jumps in the fluence contrast at the boundaries must be interpreted by the algorithm as the result of large deflections, i.e., artificially strong MHD currents. The problem is further compounded by the non-local nature of the reconstruction method; the built-in assumption of continuity and differentiability of the magnetic field that “healed” the data corruption in Sec. VI D will also result in the propagation of boundary condition errors into the domain solution. These effects increase as the RMS fluctuations of the fluence contrast increase.

The lineouts of the reconstructed field are given in Fig. 16, where we compare the actual field projection liB⎯⎯⎯⊥

(blue solid line) with the reconstructed field (red dashed line), for increasing field strengths. The lineouts are taken along the X axis at Y = 0 cm and correspond to peak values of the current projection function χpeak = 0.3 [panel (a)], χpeak = 0.6 [panel (b)], and χpeak = 1.1 [panel (c)]. The reconstruction algorithm failed to converge for the strongest field case with χpeak = 1.7 and is therefore not shown. For the nominal case, the reconstruction is able to capture the field’s characteristic length scales and features [Fig. 16(a)]. The difference between true and reconstructed fields rapidly increases as the RMS values of the fluence contrast become larger [Figs. 16(b) and 16(c)]. This is a result of the current projection function no longer being a good predictor of the fluence contrast in the nonlinear regime plus the use of the Dirichlet boundary conditions in conjunction with the non-local nature of the reconstruction method.

**VIII. ON THE NONLINEAR THEORY OF PROTON RADIOGRAPHY**

An important missing capability that remains in the analysis of proton radiographs is the ability to reconstruct the transverse magnetic field in the nonlinear regime (order unity contrast) in a manner analogous to the linear reconstruction as discussed earlier. From an experimental point of view, this is an important topic to address, since radiographs of laser plasmas can show evidence of nonlinear contrasts (see, for example, Fig. 6 of Ref. 26). In these cases, the linear theory can provide at best rough estimates of field strength and morphology.

Presentation of a complete theory of nonlinear radiography is beyond the scope of this work. However, we wish to outline an approach that can be used to reconstruct the projected perpendicular magnetic field in the nonlinear regime (i.e., the large image contrast regime). We successfully pursue this method—including a complete discussion and numerical verification of it—in a subsequent paper.22

The starting point for the approach we wish to describe is the nonlinear equation of intensity at the screen. This equation, which replaces the linear regime Eq. (4), is given by Eq. (6) of K2012, which we write here as

Ψ(x(0)+w)×det∣∣∣∣∂(x(0)+w)∂x(0)∣∣∣∣=ψ0.

(54)

Here, as before, Ψ is the measured fluence on the radiograph, ψ0 is the fluence that would be measured in the absence of any deflection, x(0) is the location on the screen where a proton would land in the absence of deflection, and w is the lateral deflection.

Now, even in the nonlinear regime, it remains the case that w = −∇

ϕ, as was discussed in Sec. IV—this is essentially a consequence of the irrotationality of the magnetic field. If we define the function

ζ≡12∣∣∣x(0)∣∣∣2−ϕ,

(55)

then we may write Eq. (54) in the following form:

Ψ(∇ζ)×det∣∣D2ζ∣∣=ψ0,

(56)

where D2ζ is the matrix of second derivatives of ζ, that is, [D2ζ]ij=∂2ζ/∂x(0)i∂x(0)j

.

As intractable as this equation appears at first glance, it turns out to be an example of a very famous class of equations: it is the Monge-Ampère equation, first discussed by the French mathematicians Gaspard Monge (in 1784) and André-Marie Ampère (in 1820). It is an equation of central importance in the subject of Monge-Kantorovich optimal mass transportation,27 and it crops up in fields as diverse as differential geometry,28 meteorology,29 cosmology,30 medical imaging,31 economics,32 and many others. For this reason, the properties of this equation have been studied in depth.

The equation arises in its archetypal form in optimal transportation problems, wherein a fixed given distribution of mass in some space is to be mapped to a second fixed distribution of mass, by means of a mapping on the space that minimizes some cost function. In symbols (if not in formal mathematics), the initial distribution μ0, a measure on a space Ω, is to be mapped to a new, given distribution μ1 on Ω by means of a “transport plan,” a diffeomorphism y : Ω → Ω so that μ1 and μ2 are related by the Jacobian of y(x). Obviously, there are many possible transportation plans that can effect the required mapping, but the one of interest is the one that minimizes a certain cost function.

Suppose that the distributions μ0, μ1 are described by densities P1(x)dx, P2(x)dx, respectively. P1 and P2 are related by

P2(y)det∣∣∣∂y∂x∣∣∣=P1(x).

(57)

If the cost function to be minimized is the Kantorovich-Wasserstein L2 distance

d2(P1,P2,y)≡∫ΩdxP1(x)(y(x)−x)2,

(58)

then it has been shown33,34 that the cost-minimizing diffeomorphism y(x) is the gradient of a convex function ζ, that is,

y(x)=∇ζ(x).

(59)

Inserting Eq. (59) into Eq. (57), it then follows that ζ(x) is the solution of a Monge-Ampère equation. It has been further shown that this solution is unique.

There is therefore a curious connection between the nonlinear radiography reconstruction problem and the optimal transportation problem. Whereas in the latter, it is the case that the L2-cost minimizing optimal diffeomorphism is the gradient of the unique solution of the Monge-Ampere equation, in the radiography problem, we know that the mapping between the no-deflection fluence distribution and the observed distribution is necessarily the gradient of some function (by the irrotationality of the magnetic field) and that therefore the unique solution of the Monge-Ampère equation (56) is the optimal—by the L2 cost metric—transport plan between the uniform fluence density ψ0 and the observed density Ψ.

In other words, in our current notation, given the cost function

d2(Ψ,ψ0,w)=ψ0∫d2x∣∣∣w(x)∣∣∣2,

(60)

optimal transportation theory tells us that minimizing Eq. (60) with respect to the set of deflections w(x) that transform ψ0 to Ψ(x) leads to an irrotational vector field w(x) = −∇

ϕ(x) and hence to a solution of Eq. (56).

This connection creates an opportunity. There are a number of numerical algorithms in the literature31,35–37 that have been successfully used to solve the optimal transportation problem with L2 cost in various contexts. The successful application of any such algorithm to the radiographic inversion problem would produce the lateral deviations w and hence, by Eq. (16), the average transverse magnetic field. In a subsequent paper, we successfully pursue this method and provide a detailed and complete discussion that includes applications and numerical verification.22

**IX. LIMITATIONS OF CAUSTIC ANALYSIS**

The theory developed in this article assumes linear contrast departures, which means small field gradients and, particularly, small currents. As μ approaches order unity values and, consequently, the image contrast becomes of order unity, we saw in Sec. VIII that a different approach is warranted that involves the solution of the Monge-Ampère equation (56). In these regimes, there is no strong focusing, and the caustic structures examined and described in K2012 cannot develop in radiographic images. These are only expected for μ > 1, when the proton trajectories cross.22 Nonetheless, it is noteworthy that high-definition structures can develop in images, even in the absence of caustic structure. In fact it is apparent that the appearance of sharp features in a proton radiograph cannot be taken as evidence of caustic formation. Furthermore, as we have seen the essential tool for understanding radiographic structure is not caustic structure but rather the MHD current.

As an aid to the discussion, consider the “wire-in-a-beer-can” field configuration, wherein a uniform current flows along a central wire of finite radius, and the return current flows back along the finite-thickness walls of the can. The wire diameter is r0, the can diameter is r1, the thickness of the can walls is δ, and the height of the can is z0. The net current in the central wire is I. We assume the current density is uniform in the wire and in the can’s vertical side walls.

The magnetic field B is purely azimuthal, B = Bϕeϕ, and we assume azimuthal symmetry so that ∇

· B = r−1∂Bϕ/∂

ϕ = 0.

We choose a factorized form for Bϕ,

Bϕ=f(r)g(z),

(61)

where r is the cylindrical radius coordinate and z is cylindrical height coordinate.

The z-dependence of the field simply limits the field to the can vertically,

g(z)=⎧⎩⎨⎪⎪1,1−|z|−z0/2δ,0,|z|<z0/2,z0/2≤|z|<z0/2+δ,|z|≥z0/2+δ.

(62)

The function f(r) is set by the can/wire geometry described earlier and by Ampere’s law,

Bϕ=2i(r)cr,

(63)

where i(r) is the total net current along the z-direction inside radius r. By considering the region |z| < z0/2, we therefore have

f(r)=⎧⎩⎨⎪⎪⎪⎪2Ircr20,2Icr,2Icri(1−r−riδ),0,r<r0;r0≤r<r1;r1≤r<r1+δ;r≥ri+δ.

(64)

This field configuration yields a uniform current density in the wire and in the vertical walls of the can. The current J bends to flow around the end-caps from the wire to the side walls and back. The azimuthal magnetic field strength rises linearly from the center to the edge of the wire, then declines as r−1 in the interior of the can, and finally drops linearly to zero in the outer wall of the can.

Figure 17 gives two simulated (with 5 × 106 protons) radiographic images of “wire-in-a-beer-can” simulations. The current projection is negative in the central wire, positive in the walls. In the left panel, the parameters are chosen so that the fluctuations are linear—the contrast is everywhere less than about 0.3. In the right panel, the field intensity has been increased to push the contrasts well into the nonlinear regime so that the central region is nearly cleared of protons and the walls have more than double the mean count rate.

There is no question of any caustic structure in the first image, since the contrasts are too small for the necessary nonlinear structure to be present. In the second image, it is likely that a caustic has developed somewhere in the outer ring, given the high contrast there. Nonetheless, we see very similar morphology in the two cases. This illustrates the point that caustic formation is in no sense necessary for the formation of sharp radiographic structure. We conclude that it is not in general correct to identify high-definition features in radiographs with caustic structures, in an effort to infer magnetic structure using the formulae in K2012 to relate such structure to image caustics. At a minimum, one should verify the necessary condition for caustic formation—large contrast fluctuations.

The importance of using MHD current to interpret radiographs can be illustrated if we hypothesize interpretations of the field configuration that resulted in the radiographs of Fig. 17, which in prima facie may appear intuitive but lead to incorrect inferences. In particular, it would be incorrect to associate the peak field intensity with the peak count rate, since according to Eq. (64) the field strength peaks at the surface of the wire, rather than at the side walls of the can. A somewhat less naive interpretation would be to associate field gradients only with the interiors of the images of the wire (the central hole) and the wall (the outer annulus). This would also be incorrect; however, since according to Eq. (64), |df/dr| is also nonzero in the interior of the can (r0 ≤ r < r1), where there is no image contrast at all! Furthermore, on the basis of either image, one might naively ascribe the same (zero) field to the exterior of the can (r > r1 + δ) as to the interior (r0 ≤ r < r1), since the images of both regions are uniform-no-contrast. In fact, the interior of the can contains a non-zero irrotational field, while the exterior field is zero.

All these features are perfectly intelligible when the current is considered. The current along the proton beam direction is negative in the central wire, positive in the walls, and zero elsewhere. That statement gives a completely satisfactory and concise description of both the linear and the nonlinear image.

A further limitation of a caustic-geared analysis is that it would lead the analyst to emphasize the bright annulus of the wall, where the strong focusing occurs, over the hole in the center, where strong de-focusing takes place. But the two are complementary, as can be seen by considering how the image would change where the current simply reversed everywhere—the central hole would become a bright spot (and presumably the locus of a caustic in the nonlinear image), the wall a trough. The complementarity between the two features, which is obvious when their interpretation is in terms of current, is lost when one adopts only caustics as an interpretive key.

This leads to a final point. It is current, not the caustic structure, that makes these images intelligible. Even in the non-linear regime, despite the invalidation of the linear theory, we see that current continues to be a valuable guide to the image structure. This consideration is important for image analysis but vital for experimental design: the decision on where to position and how to orient a radiographic setup should be made, to the maximum extent possible, on the basis of the morphology of the MHD currents that one expects to be produced in the experiment. Setting up the proton beam largely perpendicularly to expected current flows would result in avoidably weak signal; contrariwise, the signal may be optimized by aligning beam and currents to the greatest extent possible. This highlights the importance of high-fidelity simulations of laser plasma experiments, such as the ones described in the study of Tzeferacos et al.38,39 Such simulations allow one to discern the main current flows in the plasma and hence optimize one’s radiographic setup with respect to the experimental setup.

**X. DISCUSSION**

We have shown in this work several new and interesting results in the theory of proton radiography. First, we have shown that in the linear regime (i.e., the small image contrast regime), proton radiographic images of magnetized (non-electrified) plasmas are simply projective images of MHD current. This fact, which was not previously recognized, lends considerably more diagnostic power to the analysis of such radiographs, which may be interpreted directly in terms of the currents that they image. This is a much more direct and less cumbersome analysis than that proposed in K2012,20 wherein radiographs are interpreted using assemblies of field configurations to mimic structures observed in the image.

Second, we have shown that, in the linear regime, it is possible to reconstruct the line-integrated transverse magnetic field by solving a steady-state diffusion equation with a source and an inhomogeneous diffusion coefficient that are both functions of the proton contrast image on the radiograph. Given an independent estimate of the interaction region size li, this translates directly into estimates of the transverse magnetic field, averaged over proton trajectories.

Third, we have shown that proton-radiographic images of isotropic-homogeneous turbulent magnetic fields can be straightforwardly analyzed to infer the field spectrum, which is obtainable from the Fourier transform of the two-point autocorrelation function of the image. The average magnetic field energy density is also easily obtainable from this kind of analysis.

Fourth, by assessing the range of validity of the linear reconstruction method, we provide guidance to experimental groups for obtaining path integrated magnetic fields from proton radiographs. While the applicability of the method depends on the particular physical situation, we provide insights on when the reconstructed field is expected to match the true field in the experiment by examining two representative cases, the isolated blob and the Gaussian random magnetic field. These test cases show the method can be expected to work well when the RMS fluctuations of the fluence contrast are small and when large contrast regions lie well within the radiographic image, as was the case for the isolated blob test problem. For the stochastic field problem, on the other hand, rapidly increasing differences between the reconstructed field and the true field occur with increasing RMS values of the fluence contrast. This is a result of the current projection function no longer being a good predictor of the fluence contrast in the nonlinear regime plus the use of the Dirichlet boundary conditions combined with the non-local nature of the reconstruction method.

Fifth, and perhaps most importantly, we have propounded a view of proton radiograph analysis that is rather different from what has become a standard view, to wit that structures on radiographs are explainable by caustic formation, and that identification of such structures and their analysis in terms of caustics is necessary to estimate magnetic structure from images. That view became prevalent, we believe, because the direct relation between radiographic structure and MHD current structure was not yet understood, and it was felt that the explanation of sharply defined structures on radiographs could not be carried out merely in terms of presumably smooth magnetic fields. For this reason, appeal was made to caustic properties of optical systems to describe the observed sharply defined structures.

In the view described in this work, while caustics are undoubtedly a possible feature of proton radiographic images, it is not generally necessary to focus on caustic structure to interpret such images. As we have shown, sharp current features such as filaments and sheets automatically result in analogous sharp structures on the radiograph, whose presence may therefore not be the result of the caustic structure. In the linear and nonlinear contrast regimes, the transverse magnetic field is directly recoverable from the image22 (as is the energy spectrum, in the case of homogenous-isotropic turbulence). From an aesthetic point of view, it seems to us much more satisfactory to be able to interpret radiographic images not as entangled products of the plasma structure and the focusing properties of the proton optical system but rather as the direct result of the plasma structure alone.

Finally, we have outlined an approach based on the Monge-Ampère equation that is capable of reconstructing the projected perpendicular magnetic field in the nonlinear regime (i.e., the large image contrast regime). We present a complete discussion of this approach in a subsequent paper.

**ACKNOWLEDGMENTS**

This work was supported in part at the University of Chicago by the U.S. Department of Energy (DOE) under Contract No. B523820 to the NNSA ASC/Alliances Center for Astrophysical Thermonuclear Flashes; the Office of Advanced Scientific Computing Research, Office of Science, U.S. DOE, under Contract No. DE-AC02-06CH11357; the Office of Fusion Energy Sciences, Office of Science, U.S. DOE, through Grant No. DE-SC0016566; the U.S. DOE NNSA ASC through the Argonne Institute for Computing in Science under field work proposal 57789; and the U.S. National Science Foundation under Grant Nos. PHY-0903997 and PHY-1619573. Carlo Graziani would like to thank Fausto Cattaneo for very useful advice on numerical matters. The authors are grateful to Alemayehu Solomon Bogale for making the reconstruction code suitable for release.