Particle-in-cell simulation of plasma-based amplification using a moving window

Current high-power laser amplifiers use chirped-pulse amplification to prevent damage to their solid-state components caused by intense electromagnetic fields. To increase laser power further requires ever larger and more expensive devices. The Raman backscatter instability in plasma facilitates an alternative amplification strategy without the limitations imposed by material damage thresholds. Plasma-based amplification has been experimentally demonstrated, but only with relatively low efficiency. Further progress requires extensive use of numerical simulations, which usually need significant computational resources. Here we present particle-in-cell (PIC) simulation techniques for accurately simulating Raman amplification using a moving window with suitable boundary conditions, reducing computational cost. We show that an analytical model for matched pump propagation in a parabolic plasma channel slightly overestimates amplification as pump laser intensity is increased. However, a method for loading data saved from separate pump-only simulations demonstrates excellent agreement with full PIC simulation. The reduction in required resources will enable parameter scans to be performed to optimize amplification, and stimulate efforts toward developing viable plasma-based laser amplifiers. The methods may also be extended to investigate Brillouin scattering, and for the development of laser wakefield accelerators. Efficient, compact, low-cost amplifiers would have widespread applications in academia and industry.


I. INTRODUCTION
There is significant international effort dedicated to developing ultrahigh-power systems for next-generation laser facilities such as the Extreme Light Infrastructure (ELI) [1,2]. These facilities are providing a platform for developing novel technologies and driving new areas of science, such as strongfield physics. Current high-power lasers are based on the chirped-pulse amplification (CPA) technique [3]. However, CPA lasers have become very large to prevent damage to optical components. One way of reducing their size is to use stimulated Raman backscattering [4,5] as a method of directly amplifying ultrashort-duration pulses in plasma [6][7][8][9][10][11][12][13]. Plasma is a robust and convenient medium as it is already broken down and can be easily replenished. It offers the potential of significantly reducing the size and cost of devices, while opening a route to exawatt powers [14]. This article introduces simulation techniques for accurately simulating amplification using a moving window, which enables parameter scans to identify suitable regimes for sustained and efficient amplification.
Raman scattering is the inelastic scattering of an electromagnetic wave by an oscillating mode excited in matter. In plasma, this is a longitudinal electron oscillation (Langmuir wave). Due to energy conservation, the scattered frequency ω 1 = ω 0 ∓ ω p can be shifted down (Stokes, −) or up (anti-Stokes, +) from the incident pump frequency ω 0 by the frequency of the excitation, in this case the plasma frequency ω p = e √ n e /(m e 0 ) (where n e is the plasma electron density, and −e and m e are the charge and mass of the electron, respectively). This frequency matching condition, which is a consequence of energy conservation and the Manley-Rowe relation, is complemented by k 0 = k 1 ± k p for the respective wave vectors, expressing momentum conservation. The Stokes process is unstable since the beat of the incident and scattered waves reinforces the electron oscillation, which in turn increases the scattering. The growth rate is particularly high for direct backscattering, and an initial counterpropagating, low-intensity seed pulse with a frequency around the downshifted Stokes frequency can be used to stimulate the scattering, resulting in amplification of the seed at the expense of pump energy. Stimulated Raman scattering has been proposed [6][7][8][9][10][11][12][13] as a mechanism for amplifying shortduration intense light pulses using a low-intensity pump laser pulse in plasma, while avoiding breakdown limits that occur in conventional solid-state amplifiers.
Raman experiments are usually quite challenging because of the need to synchronize and overlap several beams within the plasma media, which also need to be carefully tailored and prepared to ensure controlled and sustained amplification. Numerical simulations are usually performed with particle-in-cell (PIC) codes and have an essential role in the preparation of experiments, but also present a nontrivial task in themselves. Amplification usually takes place in plasma media over a length of 1-10 mm, while the beat wavelength of the lasers, λ b = λ 0 λ 1 /(λ 0 + λ 1 ) λ 0 /2, is several orders of magnitude shorter and needs to be adequately resolved. As particle kinetic effects are important, many particles per cell are required to obtain reliable results. Multidimensional simulation of the entire domain therefore requires significant computing resources. While simulations have been performed in 2D and 3D (with very few particles per cell) [24,25,36,38], many investigations are only performed in 1 dimension [22] or modeled with Vlasov or three-wave fluid models [12,32,[39][40][41][42][43][44].
In this paper, we focus on the PIC code FBPIC [45], a quasi-3D GPU-based spectral code that uses azimuthal modal decomposition to represent the 3D volume. This enables going beyond a 2D description to a cylindrical geometry that captures 3-dimensional effects such as diffraction, at a fraction of the computational cost of full 3D simulation. The spectral solver reduces sources of numerical noise, which can become pathological and artificially seed or disturb the amplification process [45,46]. The ability to run on multiple GPUs [47] makes this a fast, efficient, and scalable simulation option. However, we emphasize that the techniques developed and presented here can be implemented in other PIC codes.
A full simulation is extremely expensive (even using a quasi-3D geometry with modal decomposition) which severely restricts the ability to undertake parameter scans necessary to identify suitable regimes for experiments. This paper addresses this significant restriction by introducing methods to enable parameter scans through the use of a moving simulation window and appropriate boundary conditions. Using a moving window is not new [24], but previous efforts have not included the effects of the prior pump-plasma interaction. As will be shown, these effects are necessary to correctly model amplification in realistic systems, particularly for high pump intensities.
Although the methods presented here have been illustrated in the context of Raman-based laser amplifiers, they have wide applicability to any system that is prepared by an extended interaction. This could be the interaction of plasma with a long pump pulse (as in pump-probe investigations of Raman, Compton, or Brillouin scattering) or series of pulses (as in coherent control). In addition to the direct impact on the development of next-generation laser systems (and their wide range of applications), there are many other physical systems where the (possibly time dependent) state of a prepared (and potentially nonlinear) system is probed by a counterpropagating pulse. These powerful methods allow the same system to be probed numerous times using different probe pulses, propagation angles, or timings, without requiring simulation of the entire system each time. Another important application is in the development of plasma-based accelerators, where shaped prepulses (or electron bunches) can be used to condition the plasma prior to the arrival of the high-intensity driver to create a plasma channel, stimulate or control the injection process, or influence beam transport and emittance. These novel techniques could be implemented in other codes to benefit wider communities that rely on intensive numerical simulation to advance their research.
We begin by introducing two different simulation methodologies based on the concept of the moving window, before comparing with results of full simulation.

II. SIMULATION TECHNIQUES FOR RAMAN AMPLIFICATION
While a full simulation of the entire amplification domain is preferable as it self-consistently captures all of the details from the pump-plasma interaction, we explore ways of reducing the simulation domain by using a moving window-a smaller computational domain moving with the seed pulsewhich dramatically reduces the computational cost. The primary challenge in performing reduced simulations using a moving window is how to determine the plasma and laser conditions loaded at the leading edge of this domain, while properly accounting for prior pump-plasma interactions that are not explicitly part of the simulation. We discuss three simulation methods and their corresponding advantages and disadvantages: (i) full simulation; (ii) matched pump, using an analytical model to determine the laser and cold plasma profiles to be loaded; and (iii) data loading, a two-step process where a single full simulation is performed for the pump alone while recording the data in a plane at the location that would correspond to the front of the moving window, which is then loaded at the leading edge of multiple moving window simulations using different seed parameters.
Ideally, one would always perform full simulations to completely capture the pump-plasma-seed interactions in the entire domain, as illustrated in Fig. 1(a). This self-consistently includes the effects of plasma heating, spontaneous Raman backscatter, and pump evolution, but requires substantial computing resources. (For FBPIC, the ability to fit the entire domain on GPUs must also be considered.) However, long propagation times can impose limitations due to numerical heating, which may develop even with spectral solvers.
Since the pump interacts with the plasma over a prolonged period before the pump and seed pulses interact, loading cold, unperturbed plasma would be unphysical and may give rise to enhanced growth. The analytical matched pump model (described in Sec. III) attempts to approximate the electron response for a guided, long-duration pump laser. However, this cannot easily model heating or spontaneous scatter. The analytical estimates are loaded into the leading edge of the moving window, as illustrated in Fig. 1(b). This requires relatively small resources and therefore numerous simulations can be run rapidly. Moreover, the short interaction duration reduces unphysical effects arising from numerical heating.
As with the matched pump model, the data loading method (described in Sec. IV) utilizes a moving window. However, for each pump-plasma configuration one full-scale simulation (without seed) is required to capture the data to be loaded into the leading edge of the moving window, as illustrated in Figs. 1(c) and 1(d). This enables numerous runs with various seed parameters to be performed for those pump and plasma parameters. For GPU-based simulations, the entire simulation must still fit within memory constraints.
Further details for the matched pump model and data loading method will be presented in Secs. III and IV, respectively.

III. ANALYTICAL MATCHED PUMP PROPAGATION
In the standard Raman amplification configuration, the pump travels through the entire length of the plasma before it encounters the seed pulse. This interaction causes both the pump laser and the plasma to evolve from their initial conditions.
Here we consider a long pump laser pulse propagating in a preformed parabolic density channel. In our model, an equilibrium is reached between the ponderomotive force of the laser displacing particles and the electrostatic restoring force acting upon them. The case of a flat density profile has been considered before by Hafizi et al. [48]. In contrast, we consider guiding in an extended preformed parabolic channel produced, for example, in a capillary discharge [49,50].
Assuming a Gaussian transverse amplitude profile, where the peak dimensionless amplitude a(x, t ) is related to the electric field through for laser frequency ω, then the waist w has a stationary value w = w * , where This is calculated from the balance between the ponderomotive force of the laser pulse pushing the electrons away from regions of high laser intensity, and the electrostatic attraction back toward the positively charged ions. The solution depends on the matched waist in an unperturbed parabolic channel of radius r c with electron and ion density profiles Without ponderomotive self-focusing, the matched waistŵ is given byŵ wherer c = r c / 1 − κ 2 a 4 0 4 F 3 (2, 5 2 , 3, 3; 4, 4, 4; −κa 2 0 )/72 is the channel radius corrected for relativistic mass effects, with 4 F 3 (· · · ) a hypergeometric function. The on-axis Lorentz factor represents the averaged mass increase of the electron in the laser field, with κ = 1/2 (1) for linear (circular) polarization. The plasma wave number k p0 = 2π/λ p0 = ω p0 /c corresponds to the plasma frequency ω p0 for n e (0) = n 0 . In addition, w * depends on the matched radiusw for relativistic and ponderomotive self-focusing in homogeneous plasma with density n 0 , If we wish instead to specify the waist, w * , then the corresponding channel radius r c may be calculated fromr c 013227-3 for the desired a 0 and n 0 aŝ The model assumes the background ions to maintain the parabolic density profile (4), whereas the electron density is modified by interaction with the pump, n e (r) → n e (r) = n ion (r) + δn(r), where with σ = r 2 /w 2 * and γ = 1 + κa 2 0 exp(−2σ ).
From the requirement that the forces balance, an expression for the radial electric field can be found as from which it can be shown that with δn(r) given by Eq. (10). The radial dependence of the field is shown in Fig. 2.
Note that close to the axis, the density perturbation (11) corresponds to a reduction in electron density, as shown in Fig. 2. As discussed by Hafizi et al. [48], this is not valid when complete cavitation takes place, so we require n 0 + δn 0 such that δn/n 0 −1. This may be combined with Eqs. (12) and (13) to express the constraint as However, parabolic channel parameters can be chosen such that the constraint (14) is violated for r less than some value. In this case, the model has to be modified to ensure nonnegative electron density [51]. a 0 √ κ that can be used without complete cavitation occurring, as the on-axis density or parabolic channel radius is changed. Panel (c) shows the minimum density required to match a laser pulse with a 0 √ κ, regardless of the channel radius. The divergence of the density (and subsequent unphysical negative value) indicates that the model cannot be used for these parameters. Note that the black dotted line in all parts indicates the value a 0 √ κ = 1.86, an upper limit above which complete cavitation on the laser axis is unavoidable and the model has to be modified to avoid negative electron density. Using the relation this upper limit for model applicability corresponds to an intensity of 1.48 × 10 19 W cm −2 for λ = 800 nm. In current Raman amplification experiments, the pump laser pulse is usually significantly weaker than this, such that complete cavitation is avoided. With these restrictions in mind, as an approximation we load a laser pulse with Gaussian transverse profile characterized by the matched laser waist given by (3), along with a modulated parabolic electron profile according to Eqs. (9)-(11) and the corresponding field (12) at the leading edge of the moving window. We should note that the plasma is loaded with the correct particle momenta for the laser phase, but neglecting any longitudinal motion or detailed heating effects caused by interaction with the pump laser pulse. This method has been used to simulate Raman amplification in plasma and, as an example, Fig. 4 shows the evolution of the seed when using the monochromatic pump detailed in Sec. V. We observe the characteristic steepening as amplification occurs at the front of the pulse.

IV. PUMP DATA LOADING
In contrast to the technique discussed above, which models a guided pump laser pulse in a preformed plasma channel, we have implemented a feature that allows for data from a previous full-scale, pump-only simulation to be saved and subsequently loaded into the leading edge of a moving window simulation.
Before the pump laser reaches the seed for amplification to start, its leading edge travels through almost the entire length of the plasma. The previous model approximates the response of the plasma (and pump laser) to this interaction. In contrast, here, the pump-plasma interaction is simulated once in full for a given set of parameters, and stored to be reused multiple times to study the interaction with various seed pulses.
As illustrated in Fig. 5, a recording plane is set to move through the pump-plasma simulation domain at the location at which the front of a moving window surrounding the seed would be located if the simulation contained a seed. The counterpropagating pump pulse is launched and allowed to interact with the plasma, and timed such that the peak of the seed would meet the leading edge of the flat-top pump at the top of the plasma density ramp. The laser fields and particle positions and momenta are saved at every time step. A moving window simulation can then load the saved particle data into the full simulation at the correct moment, in addition to the evolved pump laser fields. The former allows for recovery of particle kinetic effects that may enhance or saturate amplification, while the latter accounts for energy lost by the pump laser in driving and heating the plasma, as well as to spontaneous Raman scatter.
Approximations have been made for technical reasons. While the pump can lose energy to spontaneous Raman scattering, only the counterpropagating component of the pump laser is loaded into the simulation, and therefore any spontaneous scattering that occurs as the pump travels unseeded through the plasma is removed. Currently, while the electrons and ions can be loaded with different density distributions, the corresponding self-consistent electrostatic fields caused by charge separation and required to satisfy Poisson's equation cannot be loaded by the laser antenna. However, we continue to develop the scheme to include these fields.
To present and compare the 3-dimensional laser pulses, we use the root-mean-square power function defined by where in cylindrical geometry dA ⊥ = rdrdθ . The energy in the domain at time t, or energy passing through a plane at z = z 0 , may then be calculated as Figure 6 compares the two reduced simulation methods, matched pump and data loading, with a full-scale simulation. Both of the reduced simulations are able to capture the amplification process. Figure 6(b) enlarges the moving window region highlighted in panel (a), and shows excellent agreement of the power profiles for the data loading method and the full simulation, with the matched pump overestimating amplification. The time evolution of the seed energy presented in Fig. 6(c) shows amplification taking place. The vertical dotted lines indicate where the peak of the seed pulse reaches the top of the plasma density up (down) ramp and the leading (trailing) edge of the pump laser. The dotted curve shows the total backward-propagating energy in the full simulation, including any contributions from spontaneous scatter produced by the pump, while the window energy contained in the moving window region is shown by the solid line. Excellent agreement is initially observed between the window energy and the reduced simulations; however, the analytical model noticeably overestimates amplification, while the data loading method slightly underestimates it. The spontaneously scattered energy measured in the full domain (dashed) and moving window region (solid) from the pump-only simulation (with no seed pulse) are also indicated in the figure. Final measured seed energies are 947 μJ for the full simulation, of which 801 μJ is within the moving window; 995 μJ for the matched pump; and 782 μJ with data loading. The initial and final seed energy spectra are shown in Fig. 6(d) with all showing amplification peaked slightly below the initial seed central wavelength λ 1 = 810 nm (vertical dotted line). The bandwidth and peak are well reproduced by the reduced simulations.
As noted at the beginning of this section, there are a number of assumptions made by the analytical (matched pump) simulations which may explain the overestimation of amplification. Neglecting the evolution of the pump and plasma as they interact before reaching the seed pulse prevents pump energy loss, in addition to particle kinetic effects such as plasma heating, Landau damping, and saturation processes. However, the low computational cost of the simulations makes them useful for conducting parameter scans to identify regions of potential interest.  The dashed lines in Fig. 7 illustrate the scaling of the absolute growth E 1 (the change in the seed energy) and the relative growth factor [G = E 1 /E 1 (0)] using the analytical model. Full simulation results are indicated by symbols. As can be seen, the analytical model consistently overestimates the amplification, and the discrepancy increases as the pump intensity increases. For comparison, results obtained using the data loading technique are presented by the solid curves. In this case, the agreement with the full simulation results remains excellent, highlighting the importance of the prior pump-plasma interactions as the intensity increases. These simulations required a single full simulation for each pump intensity, and the data are reused for each of the seven seed intensities simulated. We note that for these parameters the full simulations do offer a small correction to data loading re-sults, and for final investigation or comparison to experimental results these should be performed.
The simulations presented here were performed using 4 NVIDIA Tesla P100 16 GB GPUs. For the parameters considered, the full simulations took an average of 23 hours 58 minutes and required 500 GB of storage per simulation. In contrast, the analytical matched pump simulations took on average 4 hours 21 minutes to complete and required 43 GB storage. The pump-only full-scale simulations took on average 21 hours 41 minutes, and the file created to store the data at the loading plane was 14 GB. Finally, the data loading simulations took on average 3 hours 57 minutes each and, as in the analytical matched pump simulations, consumed 43 GB. As soon as there is more than one simulation to perform using the same pump laser configuration, the advantages of the data loading scheme are clear. If storage space is a concern, then again the data loading scheme can help since it only requires a relatively small data file to be stored per pump configuration. Interestingly, the analytical simulations took longer to run, indicating that evaluating the expressions is more costly than loading the data from a disk.

VI. CONCLUSIONS
We have presented analytical and data loading simulation techniques with a moving window that can efficiently be employed to numerically investigate Raman scattering processes, while taking advantage of the immense computing power offered by modern GPUs. The FBPIC particle-in-cell code has a number of features which make it a good candidate for simulating amplification, such as a quasi-3D geometry and a dispersion-free spectral solver for the electromagnetic field. The code has been extended to facilitate the methods described above, but we note that the methods above could be implemented in other PIC codes and for other processes such as Brillouin scattering and laser wakefield acceleration.
While the interaction of the pump and plasma before meeting the seed is important to the amplification process, the results above indicate that reduced simulations are adequate to investigate general properties of Raman amplification. Indeed, their strengths may be used together to produce a work of increasing complexity. Initial investigations where a large number of simulations are required, such as parameter scans, can be performed using the matched pump based on an analytical model for the pump-plasma equilibrium. Following this, full-scale pump-plasma simulations can be performed (at increased computational expense) to support a series of data loading simulations. We note that this technique can be used for arbitrary plasma channel profiles, including periodic modulations or density steps, which could be used to control amplification or to shape amplified pulses. This method remains in active development to self-consistently load the plasma fields. Finally, a small set of full-scale simulations can be performed for the desired parameter combinations for in-depth study and comparison with experimental data.
The role of numerical simulations in studying laser-plasma interactions is steadily increasing, and to make progress toward plasma-based amplifiers they are a necessary tool to develop a better understanding of the amplification process and the conditions for efficient amplification. There is currently a bottleneck in experimental studies where the efficiency of direct amplification is on the order of 10%, and spontaneous scattering occurs with similar efficiency. There is a need to explore the full parameter range so that more directed experiments can be undertaken. Numerical simulations are essential for enabling this approach. One major obstacle is the computational resources required to perform comprehensive simulations, which ultimately compromises progress. There are many remaining challenges in the development of plasmabased laser amplifiers. This paper introduces powerful methods that enable investigations that would otherwise require overwhelming computational resources, and will stimulate progress toward the realization of viable plasma-based amplifiers for next-generation laser systems.
Data are openly available online from the University of Strathclyde KnowledgeBase [52].