Laboratory evidence of dynamo amplification of magnetic fields in a turbulent plasma

Magnetic fields are ubiquitous in the Universe. The energy density of these fields is typically comparable to the energy density of the fluid motions of the plasma in which they are embedded, making magnetic fields essential players in the dynamics of the luminous matter. The standard theoretical model for the origin of these strong magnetic fields is through the amplification of tiny seed fields via turbulent dynamo to the level consistent with current observations. However, experimental demonstration of the turbulent dynamo mechanism has remained elusive, since it requires plasma conditions that are extremely hard to re-create in terrestrial laboratories. Here we demonstrate, using laser-produced colliding plasma flows, that turbulence is indeed capable of rapidly amplifying seed fields to near equipartition with the turbulent fluid motions. These results support the notion that turbulent dynamo is a viable mechanism responsible for the observed present-day magnetization.

A Thomson scattering diagnostic was used to characterize various properties of the colliding jets. The Thomson scattering beam is a 30 J, 1 ns frequency doubled (526.5 nm wavelength) laser that probes the plasma within a 50 µm 2 region. The scattered light is collected with a 63 o scattering angle and the geometry is such that the scattering wavenumber is parallel to the axis of the flow. When fitting the Thomson scattering data, a 0.025 nm instrument function was used and an electron density of 10 20 cm −3 was assumed, as inferred from absolute calibration of the total scattered signal (see below) and predicted by FLASH simulations. Within a factor of 10 from this value for the electron density, the fitted Thomson scattering spectra showed no appreciable change. The resolution of the streak camera was ∼50 ps and so the Thomson scattering data was fitted every 100 ps.
We have investigated the potential plasma-heating effect resulting from the Thomson scattering probe beam. A simple estimate of this heating is given by [1] ∆T e T e 5.32 × 10 −7 Zn e T 3/2 e λ 3 1 − e − ω/eTe I 0 τ, where T e is the electron temperature (in eV), n e the electron density (in cm −3 ), Z the mean ion charge, λ the laser wavelength (in cm), ω the laser frequency, I 0 the probe laser intensity (in W cm −2 ) and τ the laser pulse duration (in s). The above equation should be regarded as an upper bound, because it does not take into account the rapid transport of heat away from the laser focus by both convection and conduction. Indeed, the high electron thermal conductivity plays an important role in making the temperature uniform in the interaction region (see below). Taking an average value of n e ∼ 10 20 cm −3 and assuming T e ∼ 450 eV, we get ∆T e /T e 20%. More quantitatively, the probe-heating effect was fully included in our FLASH simulations which predicted similar values to these simple estimates.
In order to observe the flow properties before jet collision, a single-jet experiment was carried out. The laser beams were fired on one side of the target only, and the Thomson scattering probe beam started 28 ns, after a 10 ns long laser drive. For this shot, the Thomson scattering 1 Supplementary Figure 1: Thomson scattering from a 10 ns laser drive. Flow velocity towards grid B (full blue circles), turbulent velocity (full green squares) and electron temperature (full red diamonds) as measured by Thomson scattering for the case of a 10 ns laser drive on both a single-jet and double-jet interaction. The error bars are determined in the same way as in Figure  3b of the main text.
laser beam probed the plasma on the axis of the flow, at the midpoint between the two targets. As shown in Supplementary Figure 1, before the collision, the jets had a velocity of ∼200 km s −1 and electron temperature ∼220 eV.
For direct comparison with the single-jet experiment, we also used the Thomson scattering diagnostic to determine flow properties from a double-grid shot, but with a 10 ns pulse shape. In this case, the probe beam was initiated with the same timing and positioning as was used with the 5 ns laser drive (described in the main text). The results are also shown in Supplementary Figure 1: the turbulent velocity for the 10 ns drive is lower than the 5 ns one, as is the associated electron temperature. Both are consistent with expected properties of the respective laser drives.
We have performed a full photometric calibration of the Thomson scattering system, recovering the electron density from each shot via a comparison between the incident and scattered power. Namely, where P scatt is the scattered power, P inc is the incident power of the probe beam, dΩ is the collecting solid angle, r 0 = 2.818 × 10 −13 cm is the classical electron radius, TS = 50 µm is the interaction length,ŝ is the unit Poynting vector,Ê i0 is the probe beam's electric field unit vector, S(k) is the spectral density function, Z is the mean ion charge and n i is the ion density. The incident power is provided by the on-shot calorimetry performed at the Omega laser. The effective f-number of the optics was 9.1, giving a solid angle of 10 −2 sr. The spectral density function is [1] S(k) = Zα 4 (1 + α 2 ) [1 + α 2 (1 + Z)] , which can be calculated using the previous Thomson scattering fits to the data. Here α is 1/kλ D , with λ D the Debye length. In our experiment we have α 1; in this limit Landau damping vanishes [1]. We can therefore interpret the Thomson scattering signal using the Salpeter approximation. The shape of the ion feature (i.e., the low-frequency part of the scattering spectrum due to the ion motions) is then a function of the parameter β 2 = Z α 2 1+α 2 Te T i ∼ Z [1]. For carbon ions and nearly equal electron and ion temperatures, β ∼ 2.4 < 3.45 (see also p. 223 in Ref. 37). In this case, the real part of dielectric function has no real zeros and the spectrum consists of a wide (Doppler broadened) central feature and two peaks at the wavelengths corresponding to the minima in the dielectric function. The scattered power was determined via full photometric calibration of all the optics to convert the number of detected counts (on the CCD camera) to energy. We find that for the collided jets with a 5 ns drive, the electron density was n e = Zn i ≈ 7 × 10 19 cm −3 , in agreement with FLASH predictions (see Supplementary . The red square corresponds to the estimated electron density from Thomson scattering photometric calibration. The vertical error is determined by the standard deviation from different calibration shots. The horizontal error bar is determined by the fact the probe laser has a 1 ns pulse length.

Supplementary Method 2: X-ray Imaging and Power Spectrum of Density Fluctuations
The effects of turbulence on the X-ray emissivity were investigated using images of the plasma obtained from a framing camera. These images have a pixel size of δx ≈ 9 µm, with an estimated resolution l ≈ 50 − 100 µm from the pinhole size and detector sensitivity. We characterized relative fluctuations in detected intensity by performing spectral analysis of the jetinteraction region, visible approximately half-way between the grids (see Figure 2 of the main text). In order to distinguish small-scale density variation from large-scale inhomogeneities of the interaction region, we constructed relative-intensity maps based on 'coarse-grained' mean fields calculated via a 61×61 pixel smoothing filter -corresponding to an outer length scale L ≈ 600 µm -combined with a two-dimensional window function to remove edge effects. Power spectra were then calculated by taking average values from a binned annular histogram of a two-dimensional fast Fourier transform (FFT) applied to the relative-intensity image. To reduce the effect of spectral distortion due to image defects [2] such as neutron hits or striations, various filtering methods were employed. Extreme outliers were systematically detected by testing deviation from surrounding points, then removed using a median filter. Localized defects tend to appear as anisotropic features in Fourier space. Their impact on the assigned power-spectrum values was mitigated by using the median, rather than the mean, as a measure of central tendency for each bin: the former is more robust against the outlying values typically associated with such features. The efficacy of these methods was successfully tested by applying such defects on simulated Gaussian fields with prescribed spectra. To check the consistency of the spectra obtained, the above analysis was performed on multiple regions in the image, and various mean-field estimation methods were employed, including a variety of windowed regions. The jump observed in the experimental spectrum at high wavenumbers ( Figure 2c of the main text) is most likely due to Poisson noise. This follows because, in the two-dimensional spectrum, Poisson noise manifests itself as a constant feature: the distortion to the predicted 1D spectrum is proportional to the wavenumber.
Fluctuations in the detected X-ray intensity are typically related to density fluctuations [3], which in turn can be related to velocity perturbations under certain conditions. The radiative energy flux (emissivity) from an optically thin plasma with T = T e = T i is given by a cooling function [4] = ρκ P σT where σ is the Stefan-Boltzmann constant, ρ the mass density, and the functional form of the Planck opacity has been assumed to be κ P = κ 0 ρ α T β , with κ 0 a suitable constant. For Bremsstrahlung-dominated radiation, β = −7/2 and α = 1, resulting in weak dependence on temperature [4]. Since ρ ∝ n, we conclude that the emissivity can be written as where C is a constant. Decomposing the number density and temperature into mean and fluc-tuating parts, viz., and assuming n 0 δn, T 0 δT , we find that the emissivity becomes where the mean emissivity is 0 = Cn α+1 The intensity detected on the screen is then given by the line-of-sight integral with mean-field intensity I 0 defined by I 0 ≡ n 0 0 (x, y, z) dz, where n is the size of the plasma. Assuming n to be larger than the scale of the dominant fluctuations l, and homogeneity at small scales (an assumption discussed at the end of this subsection), we can separate the various terms with a WKB-type (Wentzel -Kramers -Brillouin) approximation to give δI I 0 (x, y) ≈ (α + 1) Note that we cannot necessarily neglect the quadratic terms in Supplementary Equation 9, because the linear integrals are reduced in magnitude by a factor (l/ n ) 1/2 due to the fluctuating integrands; this reduction is not seen in the quadratic terms. In order to ignore the latter, we would require δn/n 0 < (l/ n ) 1/2 and δT /T 0 < (l/ n ) 1/2 . Yet, for 1D power spectra shallower than k −2 this condition is not satisfied.
To proceed further, we need to use an equation of state to relate density and temperature fluctuations. For example, in the case of subsonic, ideal gas motions, it can be shown from pressure balance that This gives a quadratic relation between intensity and density fluctuations of the form provided c 1 and c 2 do not vanish. Alternatively, a similar result can be obtained by assuming that the temperature remains uniform in the plasma, thus neglecting temperature fluctuations with respect to density fluctuations. Since thermal conductivity is large in our experiment, such an assumption would be reasonable; furthermore, a relatively uniform temperature across the interaction region is predicted by FLASH simulations (see Supplementary Figure 2). The above considerations allow the (assumed homogeneous) intensity and density correlation functions, defined by to be related via where we have used typical decay conditions on the density correlation function to extend the integrals formally over the entire real line, and have neglected the cross-correlation between δn/n 0 and δn 2 /n 2 0 on the grounds that such a term contains the same reduction factor as the linear term, and so is negligible under the small fluctuation assumption δn/n 0 1. The 2D intensity spectrum P 2D (k ⊥ ), 3D density spectrum P 3D (k), and 3D squared-density spectrum P (2) 3D (k), given by the Fourier transforms of Supplementary Equations 12, 13, and 14, respectively, combined with the Wiener-Khinchin Theorem [5], are therefore related by Isotropy of the statistics of the power spectrum was invoked for the last equality, and is justified subsequently. We simplify Supplementary Equation 16 by noting that P 3D (k) can be written (17) On account of the largest density fluctuations occurring at the smallest wavenumbers, the sum on the RHS of Supplementary Equation 17 is dominated by wavemodes where k , k k. This is equivalent to assuming that the 1D power spectrum is decreasing, so that |δn k 1 | > |δn k 2 | if |k 1 | < |k 2 |; and as presented in the main text, for k 2π/ n the 1D power spectrum is indeed 6 decreasing. Under this assumption, Supplementary Equation 17 becomes (18) In short, for a decaying spectrum with small density fluctuations at k 2π/ n the contribution of P (2) 3D (k) to P 2D (k) is simply to modify its amplitude. The desired relation follows from Supplementary Equation 18 applied to Supplementary Equation 16.
In the subsonic regime, this relationship should hold more generally than for an ideal gas, provided any linearized relationship deduced between density and temperature fluctuations from the equation of state does not lead to the coefficients c 1 and c 2 vanishing in the quadraticintensity result, Supplementary Equation 11.
The approach described here was validated against numerical results. We have created a synthetic X-ray emission image (Supplementary Figure 3a) using FLASH simulation results, postprocessed with the multi-dimensional collisional-radiative spectral analysis code SPECT3D (http://www.prism-cs.com). Supplementary Equation 18 can be tested with simulated stochastic density fields, and numerical results support it. For example, Supplementary Figure 3c shows a comparison of P 3D (k) to P (2) 3D (k) for the density field resulting from FLASH simulations. At large wavenumbers, the slopes of the density and squared-density spectra are similar, in line with expectations from Supplementary Equation 18. We note that at such wavenumbers, P 3D (k) is greater in magnitude than P 3D (k), confirming the importance of keeping quadratic terms in the expansion in Supplementary Equation 9. The discrepancy at smaller wavenumbers is likely indicative of the breakdown of the decreasing-spectrum approximation inherent in Supplementary Equation 18.
Supplementary Figure 3b also provides a numerical test of Supplementary Equation 16, by comparing for FLASH simulations the predicted 2D relative intensity spectrum to the true intensity spectrum, assuming radiation is Bremsstrahlung dominated. With the exception of the lowest wavenumbers, reasonable agreement is again seen. The difference between the two spectra is presumed to arise from the fact that in the FLASH simulations the size of relative density and intensity fluctuations on the largest scales are not asymptotically small: (δn/n 0 ) rms ≈ 0.5, and (δI/I 0 ) rms ≈ 0.3. We believe that the disagreement for small wavenumbers is not a source of concern, on the grounds that the calculated spectrum at such scales is already somewhat uncertain due to k n ∼ 1. To summarize, the numerical tests using FLASH simulations show that inferring the true density spectrum from the experimental intensity spectrum using Supplementary Equation 16 seems to be a sensible approach.
In deriving Supplementary Equation 19, we assumed that the quantities δn/n 0 and δI/I 0 possess homogeneous and isotropic statistics. These are non-trivial assumptions, as the experiment does not produce an infinite plasma; and while the target is close to being symmetric Supplementary Figure 3: Testing intensity and density spectral relations in FLASH. a) Synthetic x-ray emission image for the 5 ns laser drive at t = 35 ns. The image has been generated with the code SPECT3D using FLASH simulation results. b) Comparison of the exact intensity spectrum resulting from FLASH simulations in the case of Bremsstrahlung-dominated radiation (β = −7/2 and α = 1) to the intensity spectrum predicted by Supplementary Equation 16. c) Comparison of the 3D density spectrum to the 3D squared-density spectrum, as calculated from FLASH simulations. According to Supplementary Equation 18, at sufficiently large wavenumbers these spectra should be proportional to each other. The region from the FLASH simulations used to calculate these spectra is the same as the region used to calculate the 1D density spectrum shown in Figure 2c of the main text.
around the direction parallel to the initial jet flow direction (hereon referred to as the line of centers), the directions parallel and perpendicular to the line of centers could potentially lead to anisotropic plasma statistics. However, these assumptions can be tested using the experimental images, whose horizontal axis is parallel to the line of centers, and vertical axis perpendicular. First, if emission statistics were inhomogeneous, we would expect to obtain different spectral slopes from neighboring regions. Supplementary Figure 4a shows this not to be the case, with close agreement in the spectra calculated from three adjacent regions composing the rectangular region denoted in Figure 2b of the main text. Second, for anisotropic intensity statistics, the spectral slopes obtained by averaging sectors in k-space predominantly parallel and perpendicular to the line of centers respectively would be distinct. Supplementary Figure 4b refutes this by showing that the density spectra calculated using such a procedure have the same slope.

Supplementary Method 3: Faraday Rotation
We have implemented Faraday rotation measurements in the plasma using the Thomson scattering diagnostic system already available at the Omega laser facility. Using a Wollaston prism, the Thomson scattered light beam was separated into two orthogonal components, S and P. The magnitude of the detected S and P polarizations I S and I P depends on three factors: the initial polarization of the linearly-polarized probe beam, the relative transmission of the S and P com-Supplementary Figure 4: Testing homogeneity and isotropy assumptions for density spectrum. (a) Test of homogeneity assumption for the density spectrum, performed by dividing the region of the X-ray image selected for analysis (denoted in Figure 2b of the main text by the red rectangle) into three equally sized regions, then calculating the predicted 3D density spectra for each. Blue saltires denote the upper region, squares the middle region, and circles the lower region. For comparison, the spectra are plotted with a 3D Kolmogorov power law. The upper wavenumber cutoff denotes the nominal resolution limit of the X-ray imaging diagnostic. (b) Test of isotropy assumption for the density spectrum. The test consists of calculating the two-dimensional FFT applied to the full region denoted in Figure 2b of the main text, but calculating partial rather than full binned annular histograms from the FFT. Wavenumbers are separated into those whose predominant direction is horizontal, and vertical, corresponding to directions parallel and perpendicular to the initial direction of the plasma flows (the line of centers). The parallel and perpendicular spectra (denoted by crosses and saltires respectively) are then calculated using the average values from each of the partial annular regions.
ponents through the diagnostic optics subsequent to their separation by the prism, and changes to polarization due to Faraday rotation. Namely, where I 0 is the probe laser intensity, A S and A P the (distinct) transmission factors for each polarization, θ c the angle at which the S and P components are split by the prism, and ∆θ the rotation due to the Faraday effect. Taking the ratio of Supplementary Equations 20 and 21, we find If θ c and the ratio A S /A P are known, the degree of Faraday rotation ∆θ can be calculated from In this experiment, the Faraday rotation diagnostic was set up so that the S and P components are equally split by the prism, that is θ c = 45 o . To calculate the transmission ratio, we consider the magnitude of the two components for a single jet case. The single-jet case can be used for calibration, because the proton radiography diagnostic indicates that the magnetic fields are small enough that the Faraday rotation angle is below detector sensitivity. In Supplementary Figure 5, we show the raw image of the split polarizations and plot the summed intensity of the signal for its entire duration. We determine the transmission ratio to be A S /A P = 0.43. The use of this data for calibration is validated by the fact that the measured rotation angle remains constant in time, as shown in Figure 3c of the main text. Since the resolution of the streak camera is ∼50 ps, the Faraday rotation angle for other shots is then calculated at 100 ps intervals.
When the Faraday-rotation angle is small, it is given (in Gaussian units) by [5] It is appropriate to use twice the size of the plasma n for the path length of the integral because the geometry of the jet-interaction region is such that the Thomson scattering volume lies on the opposing side to that from which the probe beam originates. The incident laser's wavelength is λ = 5.27 × 10 −5 cm, so this implies If the mean magnetic field is small or zero, the mean of the Faraday-rotation measure should vanish and so measured values of it correspond to the standard deviation of the line-of-sight integral in Supplementary Equation 24. This can be estimated by a random-walk argument: assuming a random field has correlation scale B , the typical deviation is equal to the deviation acquired across one correlated structure multiplied by the square root of the number of such structures encountered. Supplementary Equation 24 then implies where n is the length scale of the density profile and B the scale of the magnetic structures that have the largest amplitude. This estimate can be re-arranged to give the result stated in the main text.
There are a range of effects to which the probe beam is subject but that can be reasonably neglected in our analysis. The probe laser changes frequency -and hence wavelength -after scattering, due to collective plasma effects in the probe-interaction region; however, the largest such change will be of the order of the plasma frequency, and this corresponds to a maximum wavelength shift ∆λ/λ 30%. The error in Faraday rotation associated to this wavelength modification is 20%, which is less than the noise in the streak camera data. Due to the small size of the electron Larmor frequency (Ω e ≈ 2 × 10 12 s −1 ) compared to the laser frequency (ω ≈ 4 × 10 15 s −1 ), it can also be shown that other magnetic effects on the polarization -such as the Cotton-Mouton effect -are much smaller in magnitude than the experimental error [6] and can be safely omitted.

Supplementary Method 4: Proton Radiography
Proton deflections in the proton radiography diagnostic are primarily due to magnetic fields rather than electric fields or effects such as Coulomb collisions and kinetic beam instabilities. This can be seen from the energy map on the CR-39 plates. Supplementary Figure 6 shows that, in the jet-interaction region, the mean proton energy has a uniform distribution without pronounced density structures.
For arbitrary magnetic-field configurations and imaging beam parameters, relating analytical distributions of flux and the magnetic fields from which they arise is non-trivial. However, for radiographic setups similar to (and including) this experiment, a number of simplifications can be made to improve tractability of the "reconstruction" problem of extracting path-integrated magnetic fields directly from a given proton flux distribution. Based on such an approach, here we discuss the formulas used to derive the images of path-integrated fields, shown in Figure 4 of the main text, as well as the assumptions under which these formulas are valid. A full derivation of these results is presented elsewhere [7].
We calculate path-integrated fields from the experimentally obtained proton-flux images, assuming that the distribution of flux Ψ is related to the path-integrated field by a Monge- In this expression, ∇ ≡ ∂/∂x ⊥0 is a gradient operator defined on a two-dimensional coordinate system x ⊥0 co-planar with the imaging screen (CR-39 plate), but with a scale reduced by the imaging magnification factor where r s ≈ 27 cm is the distance from the plasma to the screen and r i ≈ 1 cm is the distance from the proton source to the plasma. Ψ 0 is the initial uniform flux distribution, but with its magnitude multiplied by a scaling factor of M −2 . ∇φ is a two-dimensional gradient map related to the path-integrated magnetic field by where e is the proton charge, m p the proton mass, v p the proton speed,ẑ a unit vector perpendicular to the imaging screen, l i the size of the plasma, s the path length, and x(s) the trajectory of a proton through the plasma with initial position x 0 = (x ⊥0 , 0) and velocity v 0 . The gradient map defined by Supplementary Equation 28 does not include explicitly the imaging magnification factor M, which is instead incorporated into the scheme via the re-scaled initial flux Ψ 0 . Although Supplementary Equation 26 is non-linear, it has a unique solution for ∇φ given any (positive) flux distribution Ψ with a boundary condition where δS is the boundary of the region containing the flux distribution of interest, andn the normal to that boundary [8].
The flux relation of Supplementary Equation 26 with boundary condition the Supplementary  Equation 29 was inverted numerically using a finite-difference scheme applied to the parabolic Monge-Ampère equation, the steady-state solution of which is the solution of the conventional Monge-Ampère equation [9]. In implementing this field-reconstruction algorithm for actual data, two filtering procedures were utilized. First, a Lucy-Richardson deconvolution algorithm was applied to flux images [10,11]. This is because the strength of path-integrated fields predicted by solving the Monge-Ampère equation is reduced in the case where strong, narrow flux features are affected by finite image resolution ("smearing"). Some of these effects, such as stochastic magnetic diffusion due to small-scale fields, are difficult to eliminate, due to unknown parameters associated with smearing. However, the effect on image resolution due to the finite size of the proton source has been shown to be a convolution of the "unsmeared" flux image resulting from a point proton source with a point-spread function whose precise form is determined by the velocity profile of the protons [7]. In the case of a fusion capsule, this pointspread function has been well-described as Gaussian, with full-width-half-maximum (FWHM) of 50M µm ≈ 0.14 cm when including the imaging-magnification factor [12,13]. Assuming such a point-spread function for the experimental data, the Lucy-Richardson deconvolution scheme was then implemented. For the particular flux images shown in Figure 4 of the main text and Supplementary Figure 8, ten iterations were found to balance optimally the recovery of the flux distribution with the undesired side effect of Poisson noise enhancement introduced by the deconvolution process. The efficacy of the deconvolution scheme is further enhanced by the general robustness of the non-linear reconstruction algorithm to small-scale noise [7]. An illustration of successful application of Lucy-Richardson deconvolution for simulated FLASH proton radiography is given later (Supplementary Figure 10).
In addition to the deconvolution, a Gaussian high-pass filter (FWHM 2.2Ml i ≈ 6.16 cm) was applied to a selected flux region, in order to remove long-scale variation in the flux distribution. This procedure is required because, when implemented practically, the initial flux produced by fusion capsules can vary by as much as 50% across a single CR-39 sample, though variations over solid angles ≤ 1.1 • are typically small compared to the mean flux [14]. Since the dominant structures remaining in the reconstruction are of the order l B ∼ 300 µm, and order-unity flux features are still retained in the filtered flux image, we conclude that our fieldreconstruction algorithm captures the energetically dominant magnetic structures.
The flux regions of interest (shown in Figure 4 of the main text between the grids and rods) typically have irregular shapes; to make them compatible with a rectangular finite-difference grid, we embed each section of flux in a larger rectangular region, with the edge filled with uniform mean flux (as calculated from the chosen sample). The initial flux is then chosen to be entirely uniform, with mean value set equal to that of the experimental flux region. The reconstruction algorithm is then applied, with Supplementary Equation 29 specified on the boundary of the larger rectangular region. With the reconstruction completed, an image of the magnitude of the calculated path-integrated field is then re-oriented to the original position of its associated flux region.
The applicability of this scheme to our radiography setup depends on various assumptions: paraxiality, small deflections, point-projection, and radiographic injectivity. Since the distance r i from the proton source to the plasma is much greater than the size n of the plasma, we can approximate the section of the beam passing through the plasma as planar, despite the fact that proton beams generated by fusion reactions in a D 2 capsule implosion generally take the form of a uniformly expanding spherical shell [13]. For a given proton, this paraxial approximation is effectively an expansion of the position and velocity of the particle in terms of the paraxial parameter the ratio of the size n of the region being imaged to the radius r i of curvature. For our experiment, we have n ≈ 0.06 cm, r i ≈ 1 cm, giving ϕ 0 ≈ 0.03 1, as required. This means that the initial velocity v 0 of a proton with initial position x 0 before interaction with magnetic fields can be written as Next, typical proton-velocity deflections are assumed small compared to the initial proton velocity. More precisely, it can be shown by integrating the equations of motion under this assumption that a proton with initial position x 0 and velocity v 0 acquires a velocity perturbation w due to magnetic forces given by The second approximation is obtained by substituting the initial velocity (Supplementary Equation 31) and neglecting O(ϕ 0 ) terms. Supplementary Equation 32 implies that the deflection velocity of a proton is perpendicular to its initial direction -that is, w · v 0 ≈ w ·ẑ ≈ 0 -and so the deflection angle of the proton is given by For our experiment, Figure 4f of the main text shows that the maximum predicted path-integrated magnetic field takes a value of 6 kG cm, which in turn gives maximum predicted deflection velocity w ≈ 5.7 × 10 7 cm s −1 . For the slowest species of proton used for imaging (3.3 MeV), we have v p ≈ 2.5 × 10 9 cm s −1 , so ∆ϕ ≤ 0.02 for all imaging protons. For magnetic fields varying over shorter length scales, this value is further reduced. The point-projection assumption -that the distance r s from the plasma to the screen is much greater than n -also holds, meaning that displacements from undeflected proton trajectories acquired inside the plasma due to magnetic forces are negligible compared to the displacements resulting from the free motion of particles after they have exited the plasma, with deflected velocity.
The mapping Supplementary Equation 26 implicitly assumes that no trajectories cross before the protons reach the screen, i.e., that the gradient map (Supplementary Equation 28) is injective. It is not necessarily clear from any particular image whether or not it is injective, an issue discussed at length in [7]. Non-injective images can occur in radiographic setups, and are typically characterized by strong flux features known as caustics, which can be hard to distinguish from strong yet injective flux structures in the presence of effects limiting spatial resolution (smearing effects), such as finite proton-source size. Applying a field-reconstruction algorithm to an image formed by a non-injective gradient map can lead to inaccurately inferred magnetic field morphologies, and underestimates of field strengths [7]. Indeed, the solution of the Monge-Ampère equation is known to provide a lower bound on the RMS of all pathintegrated field strengths that can generate a particular flux distribution [8]. However, we claim that for the 15.0 MeV proton images from our experiment, non-injectivity of the gradient map would not significantly alter the reconstructed path-integrated field. The justification of this claim is given in subsequent paragraphs.
An alternative test of the field reconstruction algorithm can be carried out by generating an artificial proton flux image, with the deflections of simulated protons determined by the reconstructed path-integrated field combined with Supplementary Equation 32 (see Supplementary  Figure 7a). The test consists in comparing these synthetic images to the experimental ones. The result is that they agree fairly well, although the strength of the strong, narrow flux features in the predicted flux image is reduced compared to the experimental one. Enhancing the path-integrated field by 20% recovers the observed strength of flux features (see Supplementary  Figure 7b), suggesting that the reconstructed path-integrated field as shown in Figure 4f of the main text is marginally underestimated.
The imploding capsule also produces 3.3 MeV protons from DD reactions, with associated smaller velocity v ≈ 2.5 × 10 9 cm s −1 . These slower protons are recorded independently from the faster, 15.0 MeV protons; the DD equivalents of the D 3 He images (displayed in Figure 4 of the main text) are shown in Supplementary Figure 8 Despite the presence of caustics, the 3.3 MeV proton radiographs can be used to verify the assumption of injectivity of the 15.0 MeV images. The turnover time of the plasma motions (L/u L ∼ 6 ns) is much greater than the difference in transit time of the two species of protons used for imaging (≤ 0.2 ns). Therefore, to a good approximation, the 3.3 MeV and 15.0 MeV protons are deflected by the same magnetic fields. As a result, any path-integrated field reconstructed from the flux image of one proton species should be consistent with flux images obtained for the other species. This can be tested by generating a predicted image using a reconstructed field from one proton species combined with mapping (Supplementary Equation 28), but varying the initial proton speed to match that of the other proton species. When such a test is carried out for the 15.0 MeV proton reconstructions, the morphology of predicted 3.3 MeV proton images agrees well with the actual 3.3 MeV proton images -see Supplementary Figure 7c for example. Further, numerical experiments investigating the use of our reconstruction algorithm in the regime where caustic appearance cannot be detected by the test given above suggest that the predicted path-integrated field morphology remains relatively robust, though field strength can be somewhat underestimated [7]. This is consistent with the observation (see Supplementary Figures 7b and 7d) that a 20% increase in field strength still leads to similar flux features.
The spectrum of magnetic field energy E B (k) can be calculated directly from the pathintegrated magnetic field if the field's statistics are assumed isotropic [7]. Specifically, the experimental spectrum shown in Figure 4g Figure 4f of the main text. The image is created by seeding uniformly test protons in the region over which the path-integrated field is defined (marked by a dashed red line), before assigning velocities to those particles. These velocities are given by the combination of a particle's initial velocity (Supplementary Equation 31) before interacting with fields, and velocity perturbation (Supplementary Equation 32) due to magnetic fields. The particles are then allowed to propagate to the screen, with dimensions set equal to those in the experiment. The mean flux in each region is set equal to the mean flux of the experimental image. To account for finite source size, a Gaussian point spread function with FWHM 50 M µm ≈ 0.14 cm is subsequently applied to the images. (b) Same as (a), but with the strength of the path-integrated field increased by 20%. (c) 3.3-MeV flux image predicted using reconstructed path-integrated magnetic field shown in Figure 4f of the main text (cf. Supplementary Figure 8d). The path-integrated field experienced by the slower proton species is the same, but their initial perpendicular velocity is reduced, leading to larger deflections. (d) Same as (a), but with the strength of the path-integrated field increased by 20%.
where E path (k) is the one-dimensional spectrum of the path-integrated magnetic field. The RMS value of the magnetic field can then be evaluated by integrating over all wavenumbers to give In practice, the values of the RMS field strength given in the main paper are calculated by summing over pixel values of a Fourier-transformed sample of zero-mean path-integrated magnetic field (the region denoted by a dashed curve in Figure 4f of the main text). Prior to the jet collision, the electron-density scale length is larger, of order the grid's lateral dimensions n ≥ 3 mm (we know this from FLASH simulations). However, as the collision occurs, and the turbulent jet-interaction region forms, the density scale increases to n ∼ 0.6 mm, as stated in the main text.
We have also noted in the main text that the magnetic energy spectrum derived from the experimental proton radiography data using Supplementary Equation 34 is suppressed at small wavelengths ( 100 µm). We emphasize that this may not be indicative of the true spectrum, which could follow a flatter power law -yet this would not be reproduced using the reconstruction method. In essence, this is because the injectivity assumption required above for the reconstruction to be valid is scale-dependent. Small-scale magnetic structures are more likely to lead to the crossing of trajectories of neighboring protons than larger structures of the same field strength, even though the deflections associated with the small structures are smaller. For multi-scale stochastic magnetic fields, this means that the reconstruction algorithm applied to small-scale flux structures will produce an underestimate of small-scale path-integrated fieldswhich in turns leads to the suppression of small wavelengths in the magnetic energy spectrum. In the limit of very small structures, this lack of injectivity manifests itself diffusively, effectively leading to an additional smearing effect. This phenomenon is demonstrated below for the stochastic magnetic fields produced in the FLASH simulation.

Supplementary Method 5: The Collisional Approximation
Starting with the measured values for flow velocities, temperatures, magnetic fields, and densities as given by various experimental diagnostics, we can fully characterize the plasma state. The calculated Debye length is λ D ≈ 6.8 × 10 −7 cm, which is larger than the ion separation length n −1/3 i ≈ 3.5 × 10 −7 cm: this indicates that the average number of particles in the Debye sphere is large, and hence the plasma can be described classically. The Coulomb logarithm is found to be log Λ ≈ 7.
Assuming an underlying Maxwellian distribution, the ion-ion mean free path is λ ii ≈ 1 µm, which (as would be expected for a classical plasma) is larger than the Debye length, but smaller  than all relevant hydrodynamic length scales. This would seem to imply that the plasma is collisional. However, we also need to account for the fact that the plasma formation of the central interaction region is due to two counter-propagating streams. This configuration [16] can lead to kinetic effects such as the Weibel instability that can generate and amplify magnetic fields, as has been observed in laboratory experiments [17,18]. However, the Weibel instability can only occur in a collisionless or weakly collisional plasma, since collisions strongly suppress its growth rate [19].
To assess the role of the Weibel instability in our experiments, it is helpful to consider three temporal phases: (i) what happens in the two plasma flows before they collide; (ii) what happens when the two plasma flows initially collide, creating a region where counter-propagating flows exist; (iii) and what happens after the interaction region forms.
After their formation, the two plasma flows are strongly collisional and, before they collide, there is no region where counter-propagating flows exist. Therefore the Weibel instability does not occur during phase (i).
Counter-propagating flows exist during phase (ii). The characteristic length scale of the Weibel instability is EM ∼ c/ω pi , the ion inertial length, where ω pi is the ion plasma frequency. If the mean free path for carbon ions (which are the majority of the high-Z ions in the plasma) is smaller than EM , carbon ion-ion collisions will suppress the Weibel instability.
To assess the role of the Weibel instability during phase (ii), we use the plasma parameters for the single jet case (T e ∼ T i ∼ 220 eV, U ∼ 2 × 10 7 cm s −1 ) measured using Thomson scattering for each of the two counter-propagating plasma flows (see earlier discussion on Thomson scattering). From the FLASH simulations, we estimate a typical ion density for the flows n i ∼ 1.1 × 10 18 cm −3 . Given the plasma composition, the carbon ion number density is then n Z ∼ 5.5 × 10 17 cm −3 . Using these parameters, the characteristic length of the Weibel instability for carbon ions is: By comparison, the ion-ion mean free path for carbon collisions when the temperature of the colliding flows is much smaller than the ion energy due to bulk flow velocity [20][21][22][23] is: ∼ 160 µm Z 6 −4 M 12 a.m.u. Since mfp < EM we conclude that the Weibel instability is not important during phase (ii).
For phase (iii), we consider two spatial regions: (a) the surface layers of the interaction region that the plasma flows penetrate and (b) the remaining volume of the interaction region.
To assess the role of the Weibel instability in the surface layers of the interaction region that the plasma flows penetrate, we use the plasma parameters for the interaction region given in Supplementary Table 1 (T e ∼ T i ∼ 450 eV, n i ∼ 2.3 × 10 19 cm −3 ). The number density of the carbon ions in the interaction region is then n Z ∼ 1.2 × 10 19 cm −3 and each plasma flow collides with the interaction region with velocity U ∼ 2 × 10 7 cm s −1 . The characteristic length scale of the Weibel instability is then EM ∼ 38 µm. By comparison, the mean free path for ion-ion collisions is mfp ∼ 7 µm. Since mfp < EM the Weibel instability is not important in region (a). The plasma in the remaining volume of the interaction region is strongly collisional, so the Weibel instability is not important in region (b), either.
Thus the Weibel instability is not important during any of the three temporal phases of the experiment.
Since kinetic effects are not important, it is reasonable to describe the developed interactionregion plasma as collisional, with approximately Maxwellian distributions for both electrons and ions. Under this assumption, the temperature equilibration time between ions and electrons is calculated to be τ ie ≈ 8.4 ns, which is consistent with the assumption of a plasma in thermal equilibrium. Following on from this conclusion, we can estimate various transport coefficients. The appropriate form of these depends on the magnetization of the component species in the plasma; the Larmor frequencies for electrons and ions (assuming B ≈ 120 kG) multiplied by their respective collision times are This implies that the electrons are weakly magnetized, and ions are entirely unmagnetized. The magnetized state of the electrons has various consequences: anisotropization of the thermal and electrical conductivities and the appearance of various terms in the generalized Ohm's law associated with a finite electron Larmor radius. However, since ion flows dominate bulk motion, it might be inferred that the plasma dynamics would be relatively unaffected by this. Thus, we conclude that conventional resistive MHD provides a reasonable description of the plasma dynamics. Since T e ∼ T i , the electrical conductivity is essentially determined by the electrons, and the viscosity by ions [24]. This gives diffusivities A summary of relevant parameters related to the experiment is given in Supplementary Table 1.

Supplementary Method 6: FLASH 3D Simulations
The experimental platform was designed using FLASH radiation-MHD simulations. FLASH is a parallel, multi-physics, adaptive-mesh-refinement (AMR), finite-volume Eulerian code [25,26]. The code scales well to over a 100,000 processors, and uses a variety of parallelization techniques including domain decomposition, mesh replication, and threading to utilize hardware resources optimally. The code is publicly available (http://flash.uchicago.edu) and allows one to model problems in a wide range of disciplines including astrophysics, cosmology, combustion, fluid dynamics, turbulence, and high-energy-density laboratory plasmas (HEDLP).
To design the experiments, we performed an extensive series of moderate-fidelity 2D cylindrical FLASH radiation-MHD simulations on the Beagle 2 cluster at the University of Chicago followed by a smaller set of high-fidelity 3D FLASH radiation-MHD simulations on the Mira supercomputer at the Argonne National Laboratory. This simulation campaign is described in detail in a companion paper [27]. The numerical modeling of the experimental platform employed the entire suite of HEDLP capabilities [27,28] of the FLASH code: three-temperature (3T) MHD solvers [29], non-ideal MHD effects such as magnetic resistivity [28] and Biermann battery [30,31], heat exchange between ions and electrons, implicit (using the HYPRE library, https://computation-rnd.llnl.gov/linear solvers/software.php) thermal conduction and radiation transport in the multi-group diffusion approximation, 3T tabulated equations of state and material opacities, and laser beams that are modeled with geometric-optics ray-tracing and deposit energy via inverse Bremsstrahlung. The FLASH code has a well-established record in modeling turbulent flows, magnetized turbulence, and dynamo action. A number of articles have used the code to simulate compressible and weakly compressible hydrodynamic turbulence [32][33][34][35] and turbulent amplification of magnetic fields [36][37][38].
We have used the FLASH code to simulate the experiment and reproduce the diagnostics results as shown in Figure 4 of the main text. This is illustrated in Supplementary Figure 9. The computational domain spans 0.625 cm in the X and Y directions, and 1.25 cm along Zthe line of centers between the two targets. The spatial resolution is ∼ 25 µm. Supplementary  Figure 9 clearly shows the same trend that is seen in the experimental data: a significant increase of the magnetic field following the collision of the two plasma flows. This is reflected in both the change in morphology of the proton radiographs, and in the path-integrated magnetic fields.
The FLASH code can also be used to investigate the efficacy of some of the techniques employed for analyzing the proton radiographs: in particular, the Lucy-Richardson deconvolution scheme used to correct for a finite proton source size, and the calculation of B rms . Synthetic proton radiographs can be generated using the FLASH simulation with both a point source (Supplementary Figure 10c) and a more realistic finite 40 µm source (Supplementary Figure 10a). This allows the reconstructed path-integrated field from the point source to be compared with the reconstructed field from the finite source, both with and without the use of the deconvolution algorithm (Supplementary Figure 10). Supplementary Figure 10a shows that without the deconvolution algorithm, the reconstruction process leads to lower field strengths. However, Supplementary Figure 10b shows that with the deconvolution algorithm, a much closer match is achieved. This is also reflected in the magnetic-field spectra (Supplementary Figure 10d) predicted from the proton radiography using Supplementary Equation 34.
As mentioned in the main text, FLASH also gives the magnitude of the RMS field generated initially by the laser interaction as 4 kG. The simulated time history of the magnetic field and its amplification (peak values B max and RMS values B rms in a 500 µm control volume that tracks the plasma front and the interaction region, see also [27]) is given in Supplementary Figure 11b. In order to assess the impact of the Biermann battery mechanism in the amplification process, we repeat our calculations and turn the Biermann battery term off in the interaction region ( Supplementary Figure 11a), for times t > 23 ns. The resulting time histories (dashed lines in Supplementary Figure 11b) do not show significant difference from those done including the Biermann effect at all times (solid lines in Supplementary Figure 11b). This applies to both peak and RMS values of the magnetic field in the interaction region.
The FLASH magnetic fields can also be used to illustrate the phenomenon described above of suppression of small-wavelength structures when reconstructing magnetic spectra from proton radiographs. Supplementary Figure 11c shows the magnetic energy spectrum derived from the path-integrated field as measured by test protons used for imaging, as well as from the reconstructed path-integrated magnetic field from synthetic proton images of FLASH magnetic fields. The latter in particular shows a much steeper tail than the true FLASH spectrum, with a slope that matches the experimental results of Figure 4g of the main text.
The FLASH simulation results have also been used to recover estimates of the Faraday rotation angle, as in Figure 3c of the main text. To compute the angle, we integrate Supplementary Equation 23 along the complete path of the optical Thomson scattering beam. The latter is discretized into multiple two-segment lines, each of which connect the entry port of the diagnostic laser with a computational cell in the scattering volume, and the computational cell in the scattering volume with the collection lens. The line integrals are subsequently averaged and weighed by the electron density of the cell in which they "scattered," to yield a time-resolved Faraday rotation angle. Supplementary Figure 11d shows the resulting angle (solid blue line) for a 3 ns interval, starting at 32.5 ns as the blue line in Figure 3c of the main text. The simulated rotation angle is qualitatively similar to the experimental result and undergoes sign reversals on comparable timescales. As a consistency check, the simulated RMS angle can also be used in the scaling equation for B ,rms in the main text (assuming n e = 10 20 cm −3 and n B = 0.2 mm 2 as in the experiment). The resulting RMS magnetic field is in good agreement with the in situ measurement of B rms shown in Supplementary Figure 11b, for the same time interval, corroborating the estimate of n B . Conversely, the simulated Faraday rotation angle for the single flow case (green solid line in Supplementary Figure 11d) has negligible values, as expected: the seed magnetic fields carried by the plasma flow are too weak to incur a significant change in the polarization angle (see also the green line in Figure 3c of the main text).
Supplementary Figure 11: Magnetic field generation in the jet interaction region. (a) Rendering of the FLASH magnetic field strength at t = 40 ns for the 5 ns laser drive. The region where the Biermann battery is turned off is outlined by white boundaries. (b) Evolution of maximum field strength, B max , and the RMS field, B rms , in the control volume indicated by a square box in panel (a). The plot includes both the results of the full simulation (solid lines), and a simulation where the Biermann battery term was switched off after 23 ns (dashed lines). (c) Magnetic energy spectrum calculated inside the 500 µm control volume using three methods: using FLASH magnetic field values (red squares, see also Figure 2d of the main text); integrating the simulated magnetic field in the control volume along the proton probe path (blue diamonds); and from the reconstructed field from FLASH synthetic proton radiography images (purple circles). (d) Synthetic Faraday rotation angle from the FLASH simulations: the green solid line corresponds to the single flow case when the diagnostic is probing the seed magnetic fields that are advected into the turbulent region. The line is flat and near zero, consistent with weak magnetic fields. The blue solid line corresponds to the colliding flows case, when the diagnostic is probing the magnetic fields after their amplification (see also Figure 3c of the main text). The line exhibits time-variability and Faraday rotation angles of a few degrees, consistent with fluctuating strong magnetic fields.