{'fraction_non_alphanumeric': 0.035495974031448926, 'fraction_numerical': 0.011803886851313183, 'mean_word_length': 4.573778195488722, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Advanced quasistatic approximation

The quasistatic approximation (QSA) is an efficient method of simulating laser- and beam-driven plasma wakefield acceleration, but it becomes imprecise if some plasma particles make long longitudinal excursions in a strongly nonlinear wave, or if waves with non-zero group velocity are present in the plasma, or the plasma density gradients are sharp, or the beam shape changes rapidly. We present an extension to QSA that is free from its limitations and retains its main advantages of speed and reduced dimensionality. The new approach takes into account the exchange of information between adjacent plasma layers. We introduce the physical model, describe its numerical implementation, and compare the simulation results with available analytical solutions and other codes.


Introduction
Particle acceleration in plasma and in particular plasma wakefield acceleration is a promising direction of accelerator development [1]. In this method, a driver (a bunch of charged particles or a short laser pulse) propagates through the plasma and creates a Langmuir wave in it, which accelerates the particle bunch called a witness. Both driver and witness move at almost the speed of light, so the plasma-mediated energy exchange between them can be long-lasting and efficient.
The complex phenomena occuring in plasma wakefield accelerators can only be analysed analytically in simplest approximations, which is why numerical simulations are the main method of theoretical research [2]. Most of the interesting processes take place in a small area of space that moves with the beams. This makes simulations easier, as it allows the use of a short moving window. However, there are widely different temporal and spatial scales in this problem [3,4], which can vary by many orders of magnitude, from the wavelength of the laser pulse (about a micron) to the full acceleration length (hundreds of meters) [5][6][7]. For this reason, it is not always possible to simulate plasma wakefield acceleration with particle-in-cell (PIC) codes based on first-principle equations, and modified (e.g. boosted frame [8,9]) or reduced models have to be used [2]. One popular reduced model is the so-called quasistatic approximation (QSA) [10,11]. It has been implemented in a number of codes [12][13][14][15][16], sometimes in combination with other simplifications such as envelope equation for the laser pulse [11,[17][18][19] or fluid approximation for the plasma [20][21][22]. The QSA relies on the fact that the properties of investigated objects change much more slowly in the co-moving frame than in the laboratory frame (figure 1). This applies to properties of the plasma wave (feature 1 in figure 1), densities and velocities of particle beams (feature 2), and intensity of the laser driver (feature 3). If the beams propagate in the z-direction, the change of variables from the longitudinal coordinate z and time t to s = z, ξ = z − ct, (1) Figure 1: Illustration of objects simulated in the context of plasma wakefield acceleration: plasma wave (1), witness beam (2), laser driver (3), gradient of the plasma density (4). The simulation grids used in conventional PIC codes (5),(7) and in quasistatic codes (6),(8) for modelling beam dynamics (5), (6) and temporal evolution of the plasma wave (7), (8). Duality of interpretation of quasistatic results (9).
where c is the speed of light, results in functions that depend on s (at constant ξ) much slower than on ξ (at constant s) or on z and ct separately. While conventional PIC codes require small grid steps in both time and longitudinal coordinate to resolve rapid changes in simulated quantities (feature 5), QSA requires small steps only in ξ and allows large steps in s (feature 6). This increased grid step makes QSA codes several orders of magnitude faster than conventional PIC codes. The gain in speed is determined by the ratio of beam evolution length to the plasma or laser wavelength, depending on the driver type. Note also that in QSA, the time dependence of simulated quantities at some z can be easily converted into their dependence on z at a fixed t (feature 9), if the beam and plasma properties change insignificantly at the corresponding spatial and temporal intervals. The QSA is exact if the beams of unchanging shape propagate in a longitudinally uniform plasma. In this case, an additional symmetry appears in the problem: identical plasma particles originally located at the same transverse position r ⊥ , but different z, copy the motion of each other with some time delay. The symmetry leads to another advantage of QSA: the dimensionality of the problem reduces by one when computing the plasma response. Namely, one spatial coordinate (z) disappears merging with time t into the time-like coordinate ξ. The reduced dimensionality allows a much smaller number of "macro-particles" to be used for calculating the plasma response. These plasma "macro-particles" are not bunches of real particles as in PIC codes, but "particle jets" grouping real particles that start their motion from a given transverse position with a given initial momentum, but different z.
The third advantage of QSA is its effectiveness in simulating long-term dynamics of the plasma wave. In PIC codes, to calculate the plasma state at the next point in time, we need information about the current state of the neighbouring plasma layers. Therefore, it is only possible to simulate the temporal evolution in long regions covering several wakefield periods λ p to minimize the influence of longitudinal boundaries (feature 7 in figure 1). In QSA, the state of the neighboring layers is not needed for advancing the plasma in ξ (feature 8). The model assumes that these layers copy the considered plasma layer with some time delay or advance. The resulting gain in speed is several times the ratio λ p /∆z. Thanks to this advantage, QSA codes hold the record for simulating long-term evolution of the plasma wave [23,24].
If the beams change shape slowly, or the unperturbed plasma density is weakly dependent on z, then QSA becomes imprecise, but still applicable. The inaccuracy arises because of incorrect information exchange between neighbouring z-layers. The plasma solver assumes that plasma and beam are the same at all z, but they are not. The beam evolves as it propagates in the plasma, and the initial plasma state may be different at different z. Typically, the region of space the state of which affects the plasma behavior at some z is about the plasma wavelength in the longitudinal direction and does not expand with time, since the group velocity of the Langmuir wave in a cold plasma is zero. The relative accuracy of QSA in this case can be estimated as the ratio of the plasma wavelength to the typical distance of beam change or the scale of plasma non-uniformity. However, if some plasma particles make longer longitudinal excursions in a strongly nonlinear wave, or there are waves with non-zero group velocity in the plasma, or the plasma density gradients are sharp, the QSA becomes less accurate and the error accumulates over time.
In this paper, we present an advanced quasistatic approximation (AQSA), which is free from the limitations of conventional QSA and retains its main advantages of speed and reduced dimensionality. The new approach takes into account the exchange of information between adjacent plasma layers and therefore allows correct simulation of longitudinal plasma non-uniformities, fast particles appearing in the plasma, and waves with a non-zero group velocity, if these features are resolved by the simulation grid. In section 2, we introduce the physical model and discuss its numerical implementation using a two-dimensional geometry as an example. In section 3, we compare the simulation results with available analytical solutions and with QSA and PIC codes and show that the new model reproduces effects missing in QSA. In section 4, we discuss directions of development and the applicability area of the new model.

Physical model
We start from the Maxwell equations where E and B are electric and magnetic fields, and j and ρ are current and charge densities. The equations can be combined to yield When passing to the variables (s, ξ), the derivatives change as We assume that all quantities depend slowly on s, so the corresponding derivative is a small parameter: In QSA, all small terms containing ∂ s are omitted. Here we retain the first-order terms: where∇ = (∇ ⊥ , ∂ s + ∂ ξ ), and the subscript ⊥ denotes two-dimensional (transverse) vectors and operators. Similarly, for the magnetic field we obtain The mixed derivatives in equations (6) and (7) are responsible for propagation of free radiation.
In addition to the equations for the electromagnetic field, we must also modify the laws of motion of the plasma macro-particles, that is, the laws of changing their parameters. The parameters include transverse position r ⊥ , three momentum components p and charge q, and they must be calculated as functions of ξ for a given s (figure 2). The longitudinal position of a macro-particle is not a parameter, as it is combined with time t and moved from parameters to arguments. The mass m is proportional to the charge q and does not need a separate consideration. In usual variables, parameters of plasma particles change according to the equations where v = pc/ p 2 + m 2 c 2 is the particle velocity. Denote an arbitrary parameter by χ(s, ξ) and the righthand side of the corresponding equation in (8) by F . In the quasistatic variables (1), we obtain Equation (10) always applies to the particle in the "particle jet" which at the moment has the longitudinal coordinate s. This particle may initially (before the driver comes) be at some distance from this crosssection. Thus, equation (10) describes parameters of different physical particles at different ξ. If we neglect the derivative ∂ s , equation (10) reduces to the usual QSA [25]. To complete the system, equations (6), (7) and (10) must be supplemented by the equations for charge and current, which are the same as in QSA: where j b and ρ b are current and charge densities of the driver, A is a normalization factor, and summation is over "particle jets" crossing the cell in which we calculate j and ρ. The denominators in equation (11) appear because the contribution of a "particle jet" to the density and current depends on the macro-particle velocity in the co-moving simulation window. The equations for the particle driver are as usual (see, e.g., [25]): where ( r b⊥ , ξ b ) are coordinates of the beam macro-particle in the simulation window, and m b , q b and p b are its mass, charge and momentum. The model can be augmented with the envelope equation for the laser driver and the ponderomotive force acting on the plasma particles [11], but we will not use this option in the paper.

Numerical implementation
We have added the new features to the existing quasistatic 2D3V code LCODE [26]. The code can work in Cartesian or axisymmetric geometry. We use both geometries for testing the new model, but here we present the equations for the Cartesian geometry only. Modifications to the cylindrical geometry are straightforward [27]. In the test cases considered, three field components are identically zero (E ϕ = B r = B z = 0 in cylindrical geometry and E y = B x = B z = 0 in Cartesian geometry), so we omit the equations for these components.
We calculate the plasma response in LCODE layer-by-layer in the decreasing ξ direction (from right to left in figure 2). As we need the derivatives of the current with respect to ξ in equation (6) to calculate the field E x , we use the following predictor-corrector scheme: • move the plasma macro-particles from layer a to layer b under the action of the fields of layer a; • calculate the current and charge density in layer b; • calculate all fields in layer b; • move plasma particles from layer a to layer b under the action of the average fields of layers a and b; • again calculate the current and charge density in layer b; • again calculate all fields in layer b; • for the third time, move plasma particles from layer a to layer b under the action of the average fields.
To ensure stability of the algorithm, we solve in finite differences, instead of equation (6), the following equation for the transverse electric field: whereẼ x is some prediction for E x . We present this equation in finite differences as where the indexes l, m, n denote grid layers in x, s, ξ, respectively (figure 3), ∆x, ∆s, ∆ξ are grid steps, and half-integer indices denote the half-sum of values at nearby grid points, for example, When we calculate the fields in the layer with index n for the first time (predictor),Ẽ m x,n,l = E m x,n−1,l , otherwise (corrector)Ẽ m x,n,l = E m x,n− 1 2 ,l , with E m x,n,l found at the predictor step. In choosing a finitedifference approximation of the field derivatives with respect to ξ, we rely on the causality principle. We assume that no perturbation propagates in the opposite direction to the beams, and that all perturbations in the system are driven by the beams and propagate with them or lag behind them. In a window moving at the speed of light, this means that subsequent layers in ξ cannot affect the previous ones, so only the righthand derivatives with respect to ξ should be used. This limits the applicability area, but allows for pipeline parallelization of simulations [13]. For the transverse laplacian, we use a stable Crank-Nicolson scheme [28], which has proven itself in solving the laser envelope equation in plasma wakefield simulations [19,29]. This results in an implicit numerical scheme which, if necessary, can be generalised to the three-dimensional case using the alternating-direction implicit method [30,31].
To avoid calculating the derivative with respect to ξ in the longitudinal component of equation (6), we use the continuity equation ∂ρ ∂t + div j = 0 (16) and obtain ∂ 2 ∂x 2 + 2 which in finite differences yields Systems of equations (14) and (18) are solved using the tridiagonal matrix algorithm. The equation for the magnetic field can be simplified in two-dimensional geometry. Subtracting the z-components of the first of equations (3) and the second of equations (2) yields This equation allows the magnetic field to be calculated to a constant. In cylindrical geometry, B ϕ = E r = 0 on the axis, which defines the constant. In Cartesian geometry, we need the x-component of the first of equations (3): It is sufficient to find the constant for only one value of the transverse index l. We can choose this value (l ) at the periphery of the simulation window, where the currents and fields are low and ∂ s can be neglected.
whence we find B m y,n,l , knowing the other terms. To calculate the particle parameters, we use an implicit scheme, which for equation (10) is The boundary conditions have not changed from the documented version of LCODE [27]. The transverse boundaries are perfectly conducting walls that reflect particles but absorb their energy (change the particle momentum to some small value). There are no fields and no plasma particle motion at the right (front) boundary. Boundary conditions at the left (rear) boundary are not needed in quasistatic codes.
The beam solver is also unchanged. The equations of motion are solved with the modified Euler method (midpoint method). The fields acting on the macro-particle are linearly interpolated to the predicted macroparticle location at half the time step.

Laser pulse propagation in plasma
Let us consider the propagation of a short two-dimensional laser pulse through an underdense plasma of density n 0 in the linear regime. This test shows how perturbations with a nonzero group velocity are described by the new model.
The group velocity v g of the pulse differs from the speed of light c not only because of the plasma, but also because of a finite transverse size of the pulse. Adapting the theory [32] to the flat two-dimensional case, we obtain for the group velocity measured on the axis near the focus point where x 0 is the pulse width at the waist, k 0 is the pulse wave number, k p = 4πn 0 e 2 /(m e c 2 ) is the plasma wave number, e is the elementary charge, and m e is the electron mass. The transverse size of the laser pulse x s grows according to the paraxial theory [32]: where s R = k 0 x 2 0 /2 is the Rayleigh length, and the distance s is measured from the waist. Assume the laser pulse is linearly polarized, and its electric and magnetic fields at s = 0 are   with λ 0 = 2π/k 0 = 810 nm, x 0 = 13 µm, L = 9 µm, E x = 1.26 MV/m (figure 4). The maximum field intensity is initially at ξ c = −30 µm, which is the center of the simulation window (x, ξ) of the size 180 µm × 60 µm. The grid steps are ∆ξ = λ 0 /40, ∆x = λ 0 /4, ∆s = 25λ 0 . Unlike other QSA codes [11,18,19], the propagation of this pulse is provided by the plasma solver. The simulation results quantitatively reproduce the theoretically predicted velocity of the wave packet and its transverse dispersion (figure 5). To find the group velocity from the simulation output, we approximate the pulse envelope on the propagation axis by a Gaussian function and measure the average velocity of the maximum over the time period of 0.13s R /c, starting from the focus point. The pulse width is measured according to the formula where angle brackets denote averaging over the simulation window. Note that the longitudinal grid step ∆s λ 0 , so AQSA retains the main advantage of QSA.

Plasma with longitudinal density gradient
This test shows how the new model describes a strongly nonlinear plasma wave in the presence of density gradients. The wave is driven up to the bubble regime by a short axisymmetric electron beam (figure 6). Plasma electrons at the rear part of the bubble acquire a high longitudinal velocity (close to c), travel far forward from their original location and can transfer information from one plasma layer to another. In a longitudinally uniform plasma, this effect is negated by the same plasma properties in different cross sections, but in the presence of a density gradient it can be significant. The electron beam parameters are: root mean square (rms) transverse and longitudinal dimensions 3 µm and 10 µm, respectively, relativistic factor 4 × 10 4 , and the peak beam current 11 kA. This beam creates a bubble in the plasma, which is well suited to demonstrate the capabilities of the new model and at the same time has no features with too small spatial scale that are difficult to resolve. The longitudinal profile of the plasma density contains a region of linear growth from 1.75 × 10 17 cm −3 at s = 50 µm to 2.5 × 10 17 cm −3 at s = 220 µm, which corresponds to variation of the plasma skin depth k −1 p from 12.7 µm to 10.6 µm. The density is constant before and after this region (figure 7, top row). The simulations are carried out in the cylindrical coordinates (r, ξ) or (r, z).
We compare simulations performed with QSA code LCODE, the axisymmetric PIC code FBPIC [33], and the AQSA solver implemented on the basis of LCODE. The size of the simulation window is 250 µm in r and 170 µm in ξ. In QSA and AQSA simulations, the grid steps are ∆r = ∆ξ = 0.1 µm, ∆s = 10 µm. In FBPIC simulations, ∆r = ∆z = ∆(ct) = 0.1 µm. To compare quasistatic and PIC simulations, we must present the results in the same coordinates. Let these be the quasistatic coordinates (r, ξ). The dependencies on ξ have the meaning of time dependencies at certain z. In a uniform plasma, they can also be interpreted as beam portraits at some time, but not in plasmas with longitudinal density gradients. The standard diagnostics in FBPIC do not allow us to output the time dependencies at a fixed z, so we save the plasma state and fields every fifth time step (to simplify work with large amounts of data), if the coordinate of interest s 0 is in the simulation window. Then, for each saved state, we take the values at the five cross-sections nearest to s 0 and consider them as five consecutive points of the sought time dependence at z = s 0 .
All three approaches are consistent in the regions of uniform plasma density (figures 6 and 7, left and right columns). In the region of density gradient, the classical QSA produces a different plasma response pattern than AQSA and PIC (middle column in figure 7). Unlike QSA, AQSA reproduces the bubble elongation (feature 1 in figure 7) and the appearance of an extended high-density area near the axis (feature 2) observed in PIC simulations. Both effects are caused by longitudinally moving plasma electrons that come from upstream regions where the plasma density is lower and the bubble is longer. The reason for the different shape of the field peak in AQSA and FBPIC (feature 3) is not yet completely clear. The peak amplitude is determined by the spread of electrons in transverse velocity [25], so this difference may appear because of the different numerical heating in the two codes. The shape of the tail wave [34]   a density ridge) contains high energy plasma electrons accelerated in upstream regions. Since the plasma density is lower there, these electrons appear at larger |ξ|, and the tail wave front bends, as seen in AQSA and FBPIC simulations, but not in QSA.
To conclude, the AQSA model reproduces the features arising from the longitudinal plasma nonuniformity similarly to the PIC code FBPIC, but 80 times faster, which is close to the theoretical estimate ∆s/∆ξ = 100.

Discussion
The proposed modification of the quasi-static approximation allows us to extend its applicability area. It opens up the possibility of time-efficient investigation of the laser-plasma interaction, taking into account the physics omitted in the envelope model. The effect of plasma gradients on accelerated particles can be studied more accurately than in the standard QSA. The results obtained are in agreement with both analytical theory and first-principle simulations.
The described model is only the first step towards a more complete and time-efficient simulations of laser and plasma wakefield acceleration based on the quasistatic approximation. Obviously, the next logical step should be to include plasma particle trapping in the model. Trapping plasma particles by the wave followed by their acceleration is a widely used method of injecting particles into the wave [35][36][37]. The present model does not yet include this process. Technically, this is because particle trapping is accompanied by a reduction in the spatial scale of particle parameter change to a level that cannot be resolved. It is for this reason that we have tested the new model with positive density gradients. Negative gradients can lead to particle trapping. However, the particles to be trapped can be treated in a special way, similarly to beam particles. This approach yields qualitatively correct results with the conventional QSA [38,39], and should work even better with AQSA.
The chosen method of numerically solving the equations, although quite reliable, is unlikely the best. Therefore, another direction for improving the model is to optimise the numerical algorithm.
While QSA is very efficient in simulating long-term evolution of the plasma wave, it is not certain that AQSA will perform as well on large timescales. The difference in trajectories and parameters of individual plasma particles at different s increases with the simulation window length that determines the duration of the wave evolution. In order to resolve this difference correctly, the grid step ∆s must be reduced, which depreciates the main advantage of quasistatics.