Numerical simulations of laser-driven experiments of ion acceleration in stochastic magnetic fields

Cosmic rays have been studied for well over a century, and it was only in the 1910s that we realized they came not from the Earth's crust, but instead bombarded the Earth's upper atmosphere.^1,2 Despite more than a century of investigation, both the sites of acceleration³ and plasma processes⁴ that generate the range of energies measured in cosmic rays remain open questions and, thus, topics of active research.

The energization mechanisms that are expected to operate in astrophysical systems rely on collisionless interaction with the electric and magnetic fields supported by the tenuous plasmas typical of these environments. The most widely invoked acceleration mechanism in astrophysics is diffusive shock acceleration,^5–8 an example of a first-order Fermi mechanism. Two key requirements for such mechanisms to occur are the ability to lift a subset of particles from the thermal gas to energies where they decouple from the bulk flow, and a scattering process that can sustain the acceleration and lead to the formation of a non-thermal component, i.e., not simply a heating process. Both processes rely on the interplay between particle transport and plasma turbulence. Large-scale kinetic simulations of shocks suggest that the diffusive shock acceleration process can be efficient under a broad range of plasma environmental conditions.⁹ Moreover, satellite measurements present compelling evidence for energization due to interactions with whistler¹⁰ as well as lower-hybrid¹¹ waves in the case of electrons, or ion-cyclotron modes¹² in the case of ions. These provide possible routes to injecting particles into the shock acceleration process. More recently, laboratory experiments and simulations of collisionless shocks have shown electron and ion spectra consistent with shock acceleration predictions.^13–15 

That said, particle acceleration does not occur only at shocks. The acceleration mechanism originally proposed by Fermi¹⁶ posited that particles are accelerated through repeated interaction with moving magnetic field structures. In the original picture, charged particles propagating through a magnetized medium can interact with magnetic fields, the latter acting as mirrors changing the particle energy proportional to the mirror velocity. If a mirror is moving toward a particle, it will result in an energy increase. Conversely, if receding, the particle gets decelerated. Unlike at shocks, where there exists a clear preferred direction, here, the mirrors are assumed to be randomly oriented and the resulting energy gain averaged over many interactions is thus second order in the mirrors' velocity.^17,18 While less efficient than first-order Fermi, the large extent of interstellar and intracluster mediums can increase the likelihood of meaningful contributions from second-order Fermi acceleration. Hence, it is expected to modify, rather than dominate, the cosmic ray spectrum.¹⁹ Despite detailed theoretical studies,^20,21 laboratory evidence of second-order Fermi ion acceleration has remained elusive.

Inhomogeneous environments, which are common in astrophysical contexts, can also lead to instabilities. The lower-hybrid drift instability, which develops from misaligned gradients of density and magnetic field, is one such example.²² When these gradients are present, lower-hybrid waves can develop, generating electric fields that in turn accelerate charged particles. Theoretical models have been proposed to explain the hard x-ray tail emission in Cassiopeia, arguing the lower-hybrid drift instability provides the necessary acceleration for electrons.²³ Simulations have expanded upon the quasi-linear theoretical framework,²² and have demonstrated that electrons can be self-consistently accelerated.²⁴ Measurements of accelerated auroral electrons have been analyzed, showing lower-hybrid waves as the cause of acceleration.²⁵ Space probes have measured the lower-hybrid drift instability in the magnetotail current sheet,²⁶ showing that lower-hybrid waves are generated at the separatrix. In the laboratory, high-power laser experiments have demonstrated electron acceleration from reflected ions in analogy to the solar wind interacting with a comet.²⁷ 

While many acceleration mechanisms may be operating in astrophysical environments, laboratory experiments offer the advantage of tailoring plasma conditions to probe individual mechanisms. The objective of this paper is to model experiments that measure ion acceleration in a stochastic magnetized plasma at the GSI Helmholtz Centre for Heavy Ion Research facility.²⁸ We use a platform similar to that fielded previously to study fluctuation dynamo,^29–34 coupled with the UNIversal Linear ACcelerator (UNILAC) of the GSI to launch well-characterized ions through a stochastic magnetic field. Utilizing the radiation-magneto-hydrodynamic (MHD) code FLASH, and the particle-in-cell (PIC) code OSIRIS, we calculate the plasma conditions along the calcium–ion flight path. From the plasma conditions achieved, we can estimate the efficiency of the various acceleration mechanisms and identify the dominant one. Experimental results from a campaign were used to validate the predictive capabilities of the FLASH code.

In Sec. II, we explain the platform employed at the GSI Helmholtz Centre as well as the initialization of the FLASH simulation. Section III describes the simulation results, the temporal evolution of the plasma, plasma conditions achieved, and compares the simulated results to the experiment results. In Sec. IV, we estimate energization of the UNILAC ions from the plasma conditions obtained.

The experimental campaign was carried out at the GSI Helmholtz Institute for Heavy Ion research in Darmstadt, Germany, to observe ion acceleration in a stochastic magnetic field. The experimental configuration is shown in Fig. 1. Two 20  $μ$m thick, 1 mm diameter CH foils are separated by 1.95 mm. These foils are machined to have a grooved pattern on their inner surface, facing the opposing targets. The grooves consist of ridges 10  $μ$m in height, which are 20  $μ$m wide, with a 60  $μ$m spacing between them.

Two lasers, NHELIX (Nanosecond High Energy Laser for Ion beam eXperiments) and PHELIX (Petawatt High-Energy Laser for Heavy Ion EXperiments), are used to ablate two targets and form a colliding plasma. The UNILAC generates Ca⁺¹⁸ ions with an energy of 242 MeV (5.052 MeV per nucleon) and spread in energy of 0.1%. The ions are directed through the center of the experimental chamber, between the two targets, forming a 100  $μ$m diameter, traveling into a diamond detector on the other side of the chamber (see also Fig. 1). An optical spectroscopy diagnostic collects light from over a path integrated area with a 4.0 mm diameter, encompassing the two target, to measure electron density, electron temperature, and the RMS (root mean square) velocity.

To model this experiment, we utilized the FLASH code. FLASH^35,36 is a publicly available, high-performance computing, multi-physics, adaptive mesh refinement (AMR), finite-volume Eulerian hydrodynamics and MHD code.³⁷ FLASH has extensive high-energy-density physics (HEDP) capabilities³⁶ incorporated in recent years, which enable it to model laser–plasma experiments. HEDP experiments typically cannot be described by a thermal equilibrium between electrons, ions, and radiation. Therefore, we employ a three-temperature (3T) description to solve the extended radiation–MHD equations that describe our system. We also leverage other capabilities in FLASH, such as laser ray-tracing, radiation diffusion, and multi-material equation of state (EOS) and opacities. For a detailed explanation, see Tzeferacos et al.^36,38

The 3D FLASH simulations were configured to replicate the experiments. The simulation domain extends 1.28 mm in X and Y and 3.84 mm in Z. The domain is resolved with $∼$ 6.3  $×$ 10⁶ cells, with a spacing of 10  $μ$m, which is sufficient to fully resolve the grooves, recesses, and target thickness. Outflow boundary conditions are employed to emulate the large experimental campaign. These boundary conditions fix the gradients of the flow properties to zero, except for the normal component of the magnetic field, which is set through the solenoidality condition. The multi-group, flux-limited radiation diffusion scheme uses six energy groups that range from 0.1 eV to 100 keV. For heat exchange, magnetic resistivity, and electron thermal conductivity, the Spitzer formula are used.³⁹ The generation of the magnetic field is modeled using a generalized Ohm's law, which includes the resistive and Biermann-battery terms.^40,41

Two CH foils are initialized in the domain, using the parameters from the experimental platform. To accurately model the EOS and opacity of CH, a PROPACEOS⁴² table is used. The density of the material is set to 1.04 g cm⁻³ and the initial temperature is set to room temperature. The remaining computational domain is filled with He at $1.0×10−6$ g cm⁻³ and room temperature.

The NHELIX and PHELIX beams, each simulated using 10 000 rays, have a Gaussian radial profile with 100  $μ$m Gaussian FWHM (full-width at half-maximum) radius of the focal spot. The PHELIX's laser wavelength is set to 1053 nm, while NHELIX's is set to 1064 nm. The beam power is described by a 3 ns ramp up, 4 ns plateau, and 3 ns ramp down, with a plateau power of $7.85×109$ W.

Here, we describe the temporal evolution of the experiment as seen in the 3D FLASH simulations. The evolution consists of three distinct phases. Phase I describes the laser ablation of the targets, which generates two fully ionized, counterpropagating jets. Phase II begins when the jets collide between the targets. During this phase, 1D OSIRIS and FLASH simulations are employed to assess the impact of kinetic effects during the formation of the interaction region, which is the region bounded by two reverse shocks. Finally, Phase III describes the temporal evolution and relaxation of the interaction region.

The temporal evolution in the simulation is shown in Fig. 2. The inner face of each target is illuminated by a single beam, absorbing laser energy through electron free-free absorption, otherwise known as inverse bremsstrahlung. This results in ablation of the targets, generating two expanding, supersonic, fully ionized, asymmetric jets. The laser–target interaction generates Biermann-battery magnetic fields,^43,44 which are flux-frozen into the expanding jets and advected with the jets into the center of the domain. The jets retain a density variation, as a result of the grooved pattern present on the target surface. The properties of the jets are summarized in Table I. Approximately, 3 ns after the laser is turned on, the jets collide in the center of the domain. Phase I concludes with the collision of the two jets.

To that end, we execute a 1D OSIRIS^46–48 PIC simulation (which can capture counterstreaming and kinetic effects). These 1D simulations will result in a longer counterstreaming time than a 3D case, thus giving an upper bound on the counterstreaming phase. The OSIRIS simulation was initialized at the time when the two jets began to merge. Two counterstreaming plasma jets, extending 12.4 mm, were initialized with spatially varying density and velocity profiles in the 28 mm simulation domain. This setup ensured that the particle and electromagnetic-field boundaries were transparent for the physical region of interest (the central 2 mm, or −1 mm < z < 1 mm, region). By fine-tuning the density and velocity spatial profiles in the OSIRIS simulation, the time-varying density and velocity profiles of the inflowing plasma jets at z =  $±$ 1 mm were matched to those in the FLASH simulation. Within the 10 ns simulation time, the inflowing jet velocity decreases from $6×107$ to $6×106$ cm s⁻¹; the jet density increases from $3×1018$ to $6×1018$ cm⁻³. The grid size of the OSIRIS simulation was set to $3.9×10−5$mm, resolving the Debye length for the initial temperature. The numbers of particle-per-cell (PPC) for the three species, namely, electron, carbon, and hydrogen, were 28, 24, and 4, respectively. The choice of PPC numbers ensures the same numerical weight of the charge for all particles, and thus conserves energy and momentum for every binary collision event that is modeled by the collision package⁴⁹ in OSIRIS. In Fig. 3, we show a composite image of the jet properties from the OSIRIS PIC simulation. The interpenetration is captured in Fig. 3(a), evidenced by the net positive velocity of all particles that constitute the left jet. However, at 4 ns post-collision [Fig. 4(b)], the interaction region has stagnated at the center of the domain. A comparison of the density to the FLASH code [Fig. 3(c)] shows similar interaction regions at 3 ns. Note that the particle interpenetration is confined to the 1.5 mm region in the PIC simulations, indicating that the particles are bound between the 2.0 mm region between the two targets.

To build confidence in the 3D FLASH results, which cannot capture the short-lived counterstreaming phase, we compare 1D FLASH simulations with 1D OSIRIS simulations to assess how discrepant the plasma properties are in the interaction region. The 1D FLASH run is initialized in an identical manner as is done in the OSIRIS simulations. In Fig. 4, we compare the temporal evolution of the average electron density over a 100  $μ$m region around the center of the computational domains. Initially the electron density is higher in the FLASH calculations, as the jet is immediately stopped and cannot counterstream. However, by 2 ns FLASH and OSIRIS predict similar densities. This result suggests that while the two codes compute different plasma conditions when the flows are counterstreaming, they quickly converge to similar values. Comparing the size of the interaction region shows a smaller interaction region at first in the FLASH results, as expected. At 3 ns, the two simulations yield interaction regions of similar size as well.

At this stage, the density has increased in the center of the domain resulting in a developed interaction region whose properties are summarized in Table I. The interaction region is subsonic and expands radially. Because the inflowing plasma jets retain the variations from the grooved pattern, their interaction results in shear and stochastic fluctuations in the interaction region. While the plasma jets continue to flow toward the interaction region, at 7 ns, the laser power ramps down and turns off at 10 ns, leading to a reduction in the density, temperature, and velocity. As the simulation progresses, the plasma temperature in the center reduces from hundreds to tens of electron volts, while the electron density drops from 5  $×1019$ to 2  $×1019$ cm⁻³.

The laser–target interactions generate large temperature and density gradients, whose cross product gives rise to a Biermann-battery magnetic field. There is no initial or background magnetic field. The magnetic field is characterized by an RMS value of 165 kG and is flux-frozen into the plasma, as seen by the large magnetic Reynolds number (>1000) in Table I. The magnetic field is advected toward the center of the domain and undergoes some flux compression, which increases its value to ∼200 kG.

We leverage the ongoing experimental campaign to validate the FLASH simulation with optical spectroscopy data. To provide electron density, electron temperature, and a RMS velocity, an optical spectroscopy diagnostic was employed, gathering data over a path integrated area with a 4.0 mm diameter, encompassing the two targets. From this large collection region, a feature in the emitted light, near 470  $μm$⁠, was fit using PrismSPECT,⁵⁵ providing temperature, electron density, and RMS velocity measurements. This peak in the spectra was the only consistent feature across the experimental data and was therefore the one used for the fitting procedure. The error bars on the values were calculated using the 100 best data fits. The uncertainty of the timing of the optical spectroscopy and the fitting procedure only gives approximate values. We therefore present individual values for each plasma property from the experiment, which are representative of the conditions from 15 to 25 ns, shown in Table II under the “Experiment” column. To directly compare experimental and simulation data, we average the FLASH results over a 1.28  $×$ 1.28  $×$ 3.84 mm³ region (the entire simulation domain), masking for the unablated target. These values are then averaged from 15 to 25 ns, which corresponds to the experimental temporal window when the optical spectroscopy measurement occurred.

The experimental plasma properties and values calculated in the FLASH simulation (under the “Large Volume” column) are shown in Table II. The simulated electron density is within the error bars of the experiment, in good agreement with the spectroscopic data. The ion temperature is somewhat higher than in the experiment, which is due to the lack of the initial counterstreaming in the FLASH simulations, resulting in higher temperature when the two jets collide. This is also why the RMS velocity is somewhat lower in the FLASH simulations, as more kinetic energy is converted into thermal energy than in the experiment.

The comparison of the FLASH calculations and the experimental results shows agreement over a large volume, but the values of interest for energization are those along the calcium–ion flight path. Employing the FLASH simulations allows us to obtain the plasma conditions and magnetic field of the calcium–ion flight path. Both are required to estimate the ion acceleration and both are different between the large volume and the flight path. Since the FLASH calculations include the magnetic field generation terms and, as shown in Sec. III, the Weibel instability is not expected to affect the magnetic field in the interaction region, we can utilize the properties in the flight path (Table II) to estimate the energization of the calcium ions.

The plasma properties, through a 100  $μ$m diameter and a 1280  $μ$m long cylinder, were calculated, i.e., the calcium–ion flight path. This cylindrical region is in the center of the domain, parallel to the two targets, and the results of averaging over this volume are shown in Fig. 5 (dashed line). As expected, the ions travel through higher densities and temperatures than those measured in the large volume (solid line) region. Conversely, the RMS velocity is lower in the calcium–ion flight path because it does not include the inflowing material. The reduced RMS velocities and increased electron temperature and electron density highlight the need for the FLASH simulations, which show that there is a substantial difference between the large volume and the calcium–ion flight path and will allow a more accurate estimation of particle energization. We calculate the magnetic field in the calcium–ion flight path, finding that it peaks at the 10 ns mark, with a value of 200 kG, dropping to ∼80 kG by the 15 ns mark, and reduces ∼50 kG by the 30 ns mark. The calculated plasma properties in the calcium–ion flight path are shown in Table II.

The ion beam is not simulated in this work. The reason is twofold, FLASH cannot accurately capture the plasma and ion–beam interaction and the computational cost associated with running this form of simulation in OSIRIS is prohibitive. Instead scaling laws are used to estimate the particle energization.

From the late-time plasma conditions that FLASH predicts in the calcium–ion path, we estimate the expected acceleration from first-order and second-order Fermi mechanism, as well as the lower-hybrid drift instability. The time frame of interest is 15–25 ns, i.e., after the interaction region has relaxed to the plasma conditions, also seen in the experiment. These times were chosen as this corresponds to when the ion bunches were launched through the domain in the experiment and is also when the FLASH results match the experiment. To mirror the experimental results, we average the plasma parameters from 15 to 25 ns in the calcium–ion flight path and display them in Table II under the “Calcium–Ion Flight Path” column.

The first-order Fermi mechanism relies on shocks. Since the interaction region is subsonic, and the reverse shocks are outside the calcium–ion path, this mechanism is not expected to contribute to ion acceleration. However, second-order Fermi and lower-hybrid drift can contribute, so we examine each separately.

Using the plasma conditions and magnetic field obtained from the FLASH simulations, we estimate the energization the calcium ions would experience to be $ΔE∼$ 9.7 keV. Comparing this value to the 0.1% (or 240 keV) energy spread in the generated ions shows that this mechanism will not result in an observable change to the energy spectrum.

In this paper, we have shown results of FLASH simulations of a laser-driven experimental campaign at the GSI Helmholtz Institute for Heavy Ion Research, which aims to observe ion acceleration in stochastic magnetic fields. The laser ablation of two targets generates weakly magnetized supersonic jets, which form an interaction region at the center of the chamber. The magnetic field originates from the propagation of the Biermann-battery generated magnetic field from the laser–target interaction. While the Weibel instability may appear in the interaction region during the counterstreaming phase, we found that its growth rate is too small to have an appreciable effect on the magnetic field.

We validated the simulation results with optical spectroscopy from the experiments and then used the validated FLASH simulations to infer the plasma properties on the calcium–ion flight path. Using the plasma properties calculated using FLASH, we estimated the energization due to second-order Fermi acceleration and lower-hybrid drift instability that would be expected in this experiment. While second-order Fermi acceleration might be acting on the ions, the effect will not be observable. Conversely, the electric field from the lower-hybrid drift instability was found to have a more appreciable effect on the ion energy. The latter is expected to result in energy gains of $ΔE∼$ 264 keV, which can be resolved by the time-of-flight diagnostics at the GSI.

The results presented here are based on the experiment P-22-00089, which was performed at the target station Z6 at the GSI Helmholtzzentrum für Schwerionenforschung in Darmstadt (Germany) in the frame of FAIR Phase-0. The TDYNO team would like to thank the operators and technicians of GSI, specifically those who assisted on the UNILAC accelerator, the Z6 experimental hall, and the PHELIX and NHELIX lasers. This material was based upon work supported by the Department of Energy National Nuclear Security Administration under Award No. DE-NA0004144 through the Horton Fellowship Program at the Laboratory for Laser Energetics of the University of Rochester. The Flash Center acknowledges support by the U.S. DOE NNSA under Award Nos. DE-NA0002724, DE-NA0003605, DE-NA0003934, DE-NA0004144, DE-NA0004147, and Subcontract 630138 with LANL; the U.S. DOE Office of Science under Award No. DE-SC0021990; and the NSF under Award Nos. PHY-2033925 and PHY-2308844. A.F.A.B. received support from UKRI (Grant No. MR/W006723/1). The simulations were performed using HPC resources at the Laboratory for Laser Energetics, with support from the Laboratory for Laser Energetics HPC Group. The software used in this work was developed in part by the Flash Center for Computational Science at the University of Chicago and the University of Rochester, with support from both the U.S. DOE NNSA and the U.S. DOE Office of Science.

The authors have no conflicts to disclose.

K. Moczulski: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). H. Wen: Conceptualization (equal); Investigation (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). T. Campbell: Formal analysis (equal); Investigation (equal). A. Scopatz: Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). C. A. J. Palmer: Investigation (equal); Writing – review & editing (supporting). A. F. A. Bott: Conceptualization (equal); Investigation (supporting); Methodology (supporting); Writing – review & editing (supporting). C. D. Arrowsmith: Investigation (supporting); Writing – review & editing (supporting). K. A. Beyer: Writing – review & editing (equal). A. Blazevic: Investigation (equal). V. Bagnoud: Investigation (equal). S. Feister: Conceptualization (supporting); Writing – review & editing (supporting). J. Halliday: Methodology (supporting); Writing – review & editing (supporting). O. Karnbach: Writing – review & editing (supporting). M. Metternich: Investigation (supporting). H. Nazary: Investigation (equal); Methodology (equal). P. Neumayer: Investigation (supporting). A. Reyes: Investigation (supporting); Methodology (supporting). E. C. Hansen: Investigation (supporting); Methodology (supporting). D. Schumacher: Investigation (equal). C. Spindloe: Conceptualization (supporting); Investigation (supporting). S. Sarkar: Conceptualization (supporting); Writing – review & editing (supporting). A. R. Bell: Conceptualization (equal); Methodology (equal). R. Bingham: Conceptualization (supporting); Formal analysis (equal); Investigation (equal); Writing – review & editing (supporting). F. Miniati: Conceptualization (supporting). A. A. Schekochihin: Conceptualization (supporting); Writing – review & editing (supporting). B. Reville: Writing – review & editing (supporting). D. Q. Lamb: Conceptualization (supporting). G. Gregori: Conceptualization (equal); Investigation (supporting); Methodology (supporting); Writing – review & editing (supporting). P. Tzeferacos: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (equal); Investigation (supporting); Methodology (supporting); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (supporting); Visualization (supporting); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.