{'fraction_non_alphanumeric': 0.029518187011997713, 'fraction_numerical': 0.010029835586872341, 'mean_word_length': 4.697649186256781, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reflection-free finite volume Maxwell's solver for adaptive grids

We present a non-staggered method for the Maxwell equations in adaptively refined grids. The code is based on finite volume central scheme that preserves in a discrete form both divergence-free property of magnetic field and the Gauss law. High spatial accuracy is achieved with help of non-oscillatory extrema preserving piece-wise or piece-wise-quadratic reconstructions. The semi-discrete equations are solved by implicit-explicit Runge-Kutta method. The new adaptive grid Maxwell's solver is examined based on several 1d examples, including the an propagation of a Gaussian pulse through vacuum and partially ionised gas. Two-dimensional extension is tested with a Gaussian pulse incident on dielectric disc. Additionally, we focus on testing computational accuracy and efficiency.


Introduction
Adaptive mesh refinement (AMR) [1] allows to increase grid resolution locally when and where it is required physically relevant shortest length scale of problem under consideration. This makes AMR a highly desirable technique for kinetic simulation of multi-scale relativistic effects in plasma [2]. Our interest to AMR is motivated by recent studies of electron-positron cascades driven by ultra-strong laser fields [3], [4]. Next generation of laser facility will allow to study such processes in the laboratory [5]. Quantum cascades are very localized phenomenon which may cause stiffness of coupled Maxwell's and kinetic plasma equations. The aim of this paper is to develop an adaptive mesh solver for the Maxwell equations in a way which allows to treat a general class of localized and stiff problems.
Often the case that numerical schemes show spurious reflections when applied to adaptive (non-uniform) meshes. These artefacts can be observed as wave packet propagating through change in grid resolution splits into transmitted and reflected parts. Then these reflections can be trapped in regions of fine grid resolution or even lead to numerical instability. Reflections were also observed in Maxwell's solvers [6] and [7] which are based on a popular in plasma simulations Yee's scheme [8]. In central scheme numerical artefacts are often attributed to a non-monotonic character of grid dispersion which supports modes with different sign of the group velocity [9]. Decoupled in uniform mesh these modes can couple across the grid interface. Though dispersion of Yee's scheme is monotonic the grid staggering artificially couples distinct branches at grid interfaces [10]. Effect of spurious reflection can be mitigated to some extent by smoothing the transition between grids [9] or by high-order solution interpolation across the interface [7]. Reflections caused by high frequency modes can be reduced by smoothing of electric and magnetic fields as suggested in [11] for staggered scheme. In this paper we propose to exploit a natural dissipation properties of explicit finite volume method (FV) with complex dispersion law. Such schemes damps out most reflection-affected high-frequency modes automatically. Also conservative cell-centred FV methods [12] can be easily extended to any type of grids including adaptively refinement ones. FV discretization of the Maxwell equation have been presented in [13] where authors proposed a divergence-corrected Godunov's method based on solution of Riemann problem.
In this paper we adopt an alternative Riemann-free semi-discrete central scheme proposed in [14]. This class of methods also allows for un-split integration of stiff source terms using an appropriate time stepping routine [15]. A core part of central FV schemes is a polynomial reconstruction of the solution from cell averages computed at the previous time step. Reconstructed solution is then used to calculate numerical fluxes. Stability of FV schemes heavily relies on monotonic (non-oscillatory) property of the reconstruction enforced with the help of limiters [16]. Roughly speaking the limiting simply reverts the reconstruction to monotonic first order near discontinuities [17]. Limiters are not always good in recognition of smooth extrema where limiting is not needed.
Such unjustified limiting can be prevented with help of additional extrema indicators acting on top of limiters for detection of smooth minima or maxima [18]. An essential ingredient of any Maxwell's solver is maintaining of divergence constraints imposed by the Gauss law ∇E = 4πρ and divergence-free property of magnetic field, i.e. ∇ · B = 0. Violation of divergence conditions causes severe numerical problems like non-conservation of charge in plasma [19].
Accumulated divergence errors can be efficiently cleaned out by introducing a generalized Lagrange multiplier [20]. Application of this technique results in corrected Maxwell equations with additional terms responsible for transport of divergence errors. However for the central finite volume method is more preferable to keep an original form of the Maxwell equations on discrete level. For this reason we adopt an alternative approach based on an appropriate redistribution of numerical fluxes [21] formulated in terms of vertex-based auxiliary potentials [22].
The reminder of this paper will be organized as follows. Section 2 and Section 3 describe the central finite volume method and time stepping routine.
High order solution reconstructions with extension to adaptive meshes is discussed in Section 4. In Section 5 we compare FV-AMR solver with one based on Yee's method and discuss grid dispersion properties. The Maxwell system is extended to nonlinear problem related to a tunnel ionization of plasma in Section 6 Multi-dimension version of the code is tested in Section 7 by simulating the 2d Gaussian pulse incident onto refractive disc. Our conclusions are summarized in Section 8.

Basic equations and numerical scheme
In the following we will numerically integrate evolutionary part of the Maxwell system 1 c where E, B and J are the electric, magnetic fields and current density respectively. Remaining pair of the Maxwell equations serves a role of constraints imposed on electric and magnetic field which provide continuity equation in plasma where ρ is the charge density.
Let us cast the Maxwell equation in form of conservation law system where q is the vector of conserved variables, F (q) = {f x (q), f y (q), f z (q)} is the flux functions and S(q) is the source terms. A detailed transcription on variables reads A usual rout to derive finite volume method [12] is to integrate a conservation law over a cell Ω i,j = ∆x∆y as follows Applying the Gauss theorem in order to replace volume integral by the contour one over cell edges results in equation for the cell average where ∂Ω i,j denotes the surface of control volume, and n is the outward vector.
The cell average is defined as Calculating the contour integral in (9) with mid-point rule we arrive at semidiscrete finite volume scheme where numerical fluxes are approximated with the two-point Lax-Wendroff func- As input parameters this function takes two values q + and q − of reconstructed solution computed at both sides of cell edge (for more details about reconstruction see Section 4) where a is local speed of solution estimated using eigenvalues of the Jacobian In order to satisfy divergence conditions (4) in discrete equations (11) we adopt auxiliary potentials method [22]. Consider a minimal setting for the transverse mode The corresponding conserved variables and flux functions are given by We define vertex-based potentials as follows New fluxes then are given by averaging of two vertex potentials over cell corner Numerical schemes now reads where δ is the difference operator, i.e.
Preservation of ∇ · B = 0 can be easily proven for modified scheme Auxiliary potentials provide a transverse correction over vertexes for electric and magnetic fields which naturally incorporates divergence conditions in the scheme.

Time-stepping with IMEX Runge-Kutta
Semi-discrete formulation of the scheme (11) leaves gives a freedom in selection of suitable time stepping routine to solve an ordinary differential equa- The Runge-Kutta type methods are described by the following set of stage (26) are often implemented to solve semi-discrete equations. In this paper we con- A detailed transcription of scheme (26)(27) reads In absence of the source term this scheme reduces to a fully explicit strongstability preserving Runge-Kutta method [23]. We apply fixed point iteration to solve this non-linear system.

Reconstruction and adaptive grid
As it was noted in the previous Section t flux functions need point-wise values of solution taken from both sides of cell edges. Solution can be obtained using a simple piece-wise reconstruction where D x and D y are the discrete derivatives approximated by centred differences. In smooth regions this reconstruction provides a solution of second order spatial accuracy. In non-smooth regions additional care has to be taken about monotonicity of the reconstruction. This can be done by calculation of discrete derivatives using for example MinMod limiter where the MinMod limiter is defined as is employed with 1 ≤ θ ≤ 2. This limiter is a subject of extrema clipping effects which leads to lower-order solution representation around smooth minima or maxima. The accuracy can be restored with additional indicator for recognition of smooth extrema where limiting is not needed. In our code we have implemented an extrema-preserving limiter proposed in [18] which is shortly An accurate approximation of solution beyond second order can be achieved with higher-reconstructions. One of successful example is to combine both piecewise and quadratic polynomials as suggested in [24] q where θ i is limiter which blends different reconstructions depending on local Where coefficients A, B and C can be uniquely defined from conservation of solution over cell For uniform grid the coefficients are defined by The limiter θ i in equation (35) can be calculated as follows where L i±1/2 = L i (x i±1/2 ). Variables M and m are defined as and Extension to adaptive grid is especially easy in finite volume method since only reconstruction step has to be modified. Communication between cells happens at reconstruction stage when cell averages from neighbour cell are requested for computing discrete derivatives. This communication is organized with help of interface functions which provide information on required resolution level, either coarsened, refined or left unaffected. If refinement ratio between adjacent cell is fixed by factor of at most these transfer functions are especially simple.
For coarsening and refining we apply the same reconstruction polynomial as described above. For example from coarse level to finer one is given by Setting the x i = 0, y j = 0 (cell centre), one can immediately obtain In the opposite case of transfer from finer q l to coarse grid q l−1 the integration of reconstruction polynomial reduces to summation over finer cells

1D Gaussian pulse on adaptive grid
Let us first illustrate how grid reflections look like with example Yee's scheme.
For this reason we consider analytically of a Gaussian pulse propagating in vac- Dispersion law of the Yee scheme is real for all frequencies ω = ω r + i · 0. Thus the Yee scheme propagates all modes without being damped. The FV scheme demonstrate non-monotonic form of dispersion curve peaked at k∆x = 0.6π. 3π/4 Group velocities for both methods v g = ∂ω/∂k are given by Group velocity changes sign for waves with k∆x > 0.6π. These modes are subject or reflection across grid interface due to artificial coupling a dispersion branch which supports backward propagating waves [9] and [10]. These modes are effectively damped out provided by imaginary part of frequency ω = ω r +iγ, γ > 1. The decrement of damping is shown in Figure with grey shaded area.
It can be seen that dissipation gets noticeable already at k > k g , where k g = arccos (3/4) · ∆x −1 . However grid resolved modes in range However we most interested to resolve waves (0, k g ) remains almost unaffected by damping. As one can see in Figure, wave modes below k g are propagated faster then speed of light in FV scheme while Yee's one decelerate all modes. This property of FV method is especially beneficial for relativistic plasma simulation since it prevents spurious resonance between particles in waves.
We proceed with study of numerical accuracy measured using L 2 norm as where ∆x = L x /N x and N x is the number of grid cells. Exact solution q ex is given by (46). For adaptive grid this definition should be slightly modified as follows [25] L amr where Ω i is the volume of cell i (in 1d case Ω i = ∆x i ). We collect results for piece-wise and piece-wise-quadratic reconstructions in Table 1. Each type of reconstruction demonstrates higher order then expected from formal one, i.e.
super-convergence. Refinement is done in form presented in Figure 1, a base grid which contains one or two levels of refinement. One can see that second level of refinement does not always improve accuracy. We attribute this effect to over-resolution when discrete precision errors proportional to a number of grid cell and time steps accumulates faster then decrease of scheme error with grid refinement.

Tunnel ionization by laser pulse
In order to test the ability of our solver to handle non-linear source term we carry out simulation of dynamic ionization by laser field. I The model for initially neutral gas is adopted from [26] This model describes generation of new electrons with density N e from initially neutral gas N g (x). Here we apply the following normalization where N cr is the critical plasma density, E a = m 2 e 5 / 4 = 5.1 × 10 9 V/cm is the field strength at the Bohr radius. Tunnel ionization is governed by probability w(|E|) as given in [27] w(|E|) = 4ω a where I a and I h are the ionization potentials of atom and hydrogen respectively (here I a /I h = 1), ω a = me 4 / 3 = 4.16 × 10 16 s −1 is the atom frequency. Initial profile of neutral gas is defined by where L g = 16c/ω 0 is the characteristic size of gas layer, N 0 = 1. In order to resolve gas layer we make use of a grid with 2 refined levels with finest resolution concentrated at neutral gas profile depicted with shadowed region in Figure 3.

Two-dimensional simulation: Interaction with refractive disc
Two-dimensional tests are carried out for a problem of laser pulse incident on dielectric refractive disc [7]. Electromagnetic effects in medium are taken into account in the following form of the Maxwell equations where solved in a box with 3 refinement levels. It was assumed vacuum magnetic permeability µ = 1, the dielectric susceptibility ε depends on refractive index n = √ µε given by The refractive disc of radius R = 2 is placed at x = 0, y = 0. The laser pulse is initialized with the Gaussian profile E = a 0 · p · cos(k · r) · Env (63) where k, p are the wave and polarization vectors respectively and the envelop is defined by where r = k · r and r ⊥ = |r × k| are the parallel and perpendicular com- gives acceleration about 100times, while finer one-level grid (square markers) gives only factor of 50 but with much better accuracy. Thus a gain in computational performance is considerable however we expect much better results after implementation of adaptive time stepping and dynamic grid adaptivity.

Conclusions
A new Maxwell solver has proven to be very effective for adaptive mesh com- We seek for solution of this equation in form q =q exp(ikx + iωt), where q = {E, B}. Thus where a = 1 + i sin(k∆x)/2 and b = 1 − i sin(k∆x)/2. Thus Fourier transformed fluxes (A.2) take form The complex valued solution reads