An Exactly Energy-conserving Electromagnetic Particle-in-cell Method in Curvilinear Coordinates

In this paper, we introduce and discuss an exactly energy-conserving particle-in-cell method for arbitrary curvilinear coordinates. The flexibility provided by curvilinear coordinates enables the study of plasmas in complex-shaped domains by aligning the grid to the given geometry or by focusing grid resolution on regions of interest without overresolving the surrounding, potentially uninteresting domain. We have achieved this through the introduction of the metric tensor, the Jacobian matrix, and contravariant operators combined with an energy-conserving fully implicit solver. We demonstrate the method’s capabilities using a Python implementation to study several one- and two-dimensional test cases: the electrostatic two-stream instability, the electromagnetic Weibel instability, and the geomagnetic environment modeling reconnection challenge. The test results confirm the capability of our new method to reproduce theoretical expectations (e.g., instability growth rates) and the corresponding results obtained with a Cartesian uniform grid when using curvilinear grids. Simultaneously, we show that the method conserves energy to machine precision in all cases.


Introduction
Advanced plasma simulation methods have become a key tool for studying a wide variety of systems and phenomena throughout the fields of plasma physics and astrophysics.The difficulty of recreating the relevant plasma conditions in a laboratory environment makes simulations irreplaceable for obtaining otherwise unachievable key insight.Among plasma simulations techniques, particle-in-cell (PIC) methods stand out as particularly useful when studying plasmas where kinetic effects are important.The state of the art of PIC methods is constantly evolving with new capabilities and improvements being developed.Recently, the astrophysics and plasma physics communities have steered their interest toward more advanced PIC methods, especially involving flexible (i.e., adaptive) grids and high-stability algorithms, to allow for longterm simulations of systems where Cartesian (uniform) grids might be ill suited (for example tokamak simulations using magnetic coordinates as shown in, e.g., Jolliet et al. 2007, or accretion disks around astrophysical compact objects using spherical Kerr-Schild coordinates as shown in, e.g., Crinquand et al. 2022).A nonuniform, adaptive grid could be particularly beneficial in several cases, e.g., (i) when the physical shape of the domain of interest is nontrivial, such that a significant part of an encompassing Cartesian box would be occupied by empty space or otherwise uninteresting regions where plasma behavior is trivial or quasistatic; and (ii) when a small region of interest is embedded in a much larger domain, and the latter might be significantly overresolved (in terms of spatial and temporal resolution) due to physical requirements only imposed by the small region of interests.In both cases, employing a standard uniform grid can substantially increase the computational cost of the simulation, in the worst case making certain studies entirely unfeasible.The introduction of a nonuniform grid is one way to relax these limitations by adapting the grid to the physical setup in question and thus assigning the computational resources more efficiently (see, e.g., Fichtl et al. 2012;Delzanno et al. 2013;Chacón & Chen 2016;Stanier & Chacón 2022).
Several methods have been developed to use curvilinear grids, but most are implemented with one specific grid in mind and thus lack flexibility (see, e.g., Ringle 2011;Gonzalez-Herrero et al. 2019).Others are fully general, i.e., they work with an arbitrary grid but come with other trade-offs such as: (i) using a reduced set of equations, rather than a full electromagnetic implementation, thus limiting their applicability to the constraints of these approximations (e.g., an electrostatic model in Fichtl et al. 2012;Delzanno et al. 2013 or the Vlasov-Darwin approximation in Chacón & Chen 2016); or (ii) introducing large deviations from energy conservation in long simulations, due to discretization errors that are intrinsic of the PIC implementation with explicit methods.Energy conservation is a universal physical property, of particular interest in setups where energy is converted from one form to another.For example, when modeling instabilities or reconnection events, preserving energy is extremely important to avoid numerical heating or cooling, which could significantly influence the plasma dynamics and potentially invalidate the results (e.g., Markidis & Lapenta 2011).
In this work, we provide a method combining general curvilinear grids through the introduction of coordinate transformations, with a fully electromagnetic-PIC implementation using an exactly energy-conserving, Jacobian-free Newton-Krylov solver.The method was designed with the architecture of the ECSIM code (Lapenta 2017) in mind, to allow for easy implementation in a production-ready Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.infrastructure.To the best of our knowledge, this may represent the first production-scale code that combines the aforementioned properties (curvilinear grids and exact energy conservation).
This paper is organized as follows: the mathematical foundation of the method is described in detail in Section 2; in Section 3, a series of one-and two-dimensional tests are discussed to assess the validity of the new method; the main results, conclusions, potential applications, and next steps in this line of research are discussed in Section 4.

The Curvilinear PIC Method
The method described in this paper combines a fully implicit-PIC implementation (Chen et al. 2011;Markidis & Lapenta 2011) using finite differences and a Jacobian-free Newton-Krylov solver, with a curvilinear grid through the introduction of the metric tensor, Jacobian matrix, and covariant differentiation operators.As we will show in Section 3, the method was successfully implemented and tested in Python.This test code has been made publicly available.4

Governing Equations
The PIC method is a first-principles approach for fully kinetic plasma simulations.It describes the plasma as a coupled system between freely moving particles existing in a positionvelocity phase space and a discrete representation of electric and magnetic fields on a computational grid.
To evolve the fields in time, we use Maxwell's equations, here in centimetre-gram-second units: where E and B are the electric and magnetic fields respectively, ρ is the charge density, c is the speed of light, and I is the current density, which is computed from the particle velocities (see Section 2.3 below).Note that the PIC method only requires Faraday's law (Equation ( 3)) and Ampère's law (Equation ( 4)) to form a closed system of equations.It can easily be shown that Gauss's law for magnetism (Equation (2)) will always be satisfied if the spatial discretization on the computational grid mimics the continuous vector identity ∇ • ∇ × V = 0 for a generic vector V (Lapenta 2017).Gauss's law for E (Equation (1)), however, is not automatically satisfied in fully implicit PIC.While it is possible to construct a method that also satisfies this equation (e.g., Chen et al. 2011;Chen & Tóth 2019), we are currently only concerned with presenting a first implementation of our method, which conserves energy exactly in curvilinear coordinates and which can then be further improved upon in future work.Moreover, in our tests, we did not detect numerical artifacts linked to charge conservation, implying that numerical errors are limited.
Computational particles in the PIC method are updated using the standard equations of motion, where x and v are the particle position and velocity, q is the particle charge, and m is the particle mass.
In the PIC paradigm, the sources for Maxwell's equations (I and ρ) are collected from the particles onto the computational grid via interpolation; electromagnetic fields are interpolated at the particle positions to obtain the Lorentz force needed to update the particles.Because of this, particle and field equations are in general nonlinearly coupled.In the method described in this paper, a fully implicit nonlinear iterative solver is used to solve this set of equations, which are discretized as described in Section 2.3.

Curvilinear PIC: Coordinate Transformation
To work with a curved grid, we consider an arbitrary mapping The space described by the new coordinates ξ will be called the "logical" space, and the space described by the original set of coordinates x will be called the "physical" space.The mapping f is chosen such that the curvilinear grid in physical space becomes a regular orthogonal grid in logical space.This is depicted in Figure 1.Next, we introduce the Jacobian and its inverse, and the associated metric defined as Here, Greek indices are used to denote ξ coordinates in the logical space and Latin indices to denote x coordinates in physical space.In this way, the difficulty of working with a nonuniform curvilinear grid has been circumvented by transforming to the logical grid, at the cost of a nontrivial metric tensor and Jacobian matrix.Calculating the metric and Jacobian only has to be done once at initialization for each grid point and can subsequently be stored in memory.Note that, in our implementation, the grid is static, and therefore, the Jacobian and metric are time invariant.Due to the introduction of a nontrivial metric, the vector operators have to be changed from their usual Cartesian implementation to a covariant implementation.The only differential operator appearing in our PIC equations is the curl operator, e.g., for a generic (covariant) where J = |j| is the Jacobian determinant and ò μ ν κ is the Levi-Civita symbol.Furthermore, the metric tensor can be used to raise and lower indices, e.g., V κ = g κλ V λ .
For the particle equations of motion we keep the Cartesian description from Equations (5)-( 6).We avoid solving these equations in logical space since this would require the Jacobian and metric tensor to be known at each particle position at every time step.Given the computational complexity of deriving these values from the initial mapping function, this was not considered a practical operation to be carried out at every iteration for every particle.This operation is in principle feasible if all mapping quantities are analytically known, but this imposes additional constraints on the mapping function (i.e., invertibility and differentiability), which we want to avoid for generality.Calculating the equations of motion in Cartesian coordinates instead introduces the need to convert particle positions and vector-field components between physical and logical spaces.This strategy is still less computationally demanding than evolving particles in the logical space.Particle positions are transformed from Cartesian physical space to general logical space simply using the original mapping function f in Equation (7).Contravariant vector components can be transformed using the Jacobian and inverse-Jacobian matrices as This is necessary in the interpolation of the current density and the fields, as we show in the following section.For further details or interest in the topic we recommend the book by Liseikin (1999).

Discretization
The notation in this paper uses superscript parentheses to denote discrete time steps and subscript parentheses to denote discrete spatial coordinates.All quantities are known at integer time steps; when the method requires half-integer times (see Equations ( 16)-( 20)), the relevant quantities will be linearly interpolated between subsequent time levels.Field quantities are defined at discrete locations on the computational grid, e.g., in two spatial dimensions, x h V n , ( ) ( ) denotes a quantity V at the grid location (ξ, η) in the logical space and discrete time level n.This allows us to avoid confusion with indices of vector components.When the exact spatial indices are not relevant, they are replaced by subscript g to indicate a quantity known at a generic grid point or subscript p for a quantity known at a particle's position.
To discretize the governing equations, we use a finitedifference scheme.Magnetic-field components are placed on grid nodes, with integer spatial indices, while electric-field components are placed on the grid centers, with half-integer spatial indices.5 Maxwell's Equations (3)-( 4) require spatial derivatives of the magnetic field to be collocated with the electric field and conversely the spatial derivatives of the electric field to be collocated with the magnetic field.This requires the quantities to be averaged along the direction perpendicular to that of the derivative before taking finite differences.This allows for a centered finite-difference spatial derivative, which is second-order accurate; in two dimensions, the discrete derivatives for a generic quantity V along the ξ or η directions are where the spatial indices can be integers or half-integers depending on whether the specific quantity exists on nodes or centers respectively.Time derivatives are discretized using a second-order accurate central-difference scheme.Using these discrete derivatives, we can now write the discrete Maxwell's Equations (3)-( 4): The current density on each grid point must be gathered from the particles as where w pg is the interpolation function linking particle p with the grid point g.w pg is a function of the particle position in logical coordinates; this can be found using the mapping function ξ = f (x).The inverse-Jacobian factor in the equation above transforms the current from physical to logical components as described by Equation (13).
Since particles are evolved in the physical space, the discretization is straightforward: Note that the fields in these equations are located at the particle positions, which requires an interpolation: for a generic field V, where w gp is the interpolation function from grid points to particle positions.By construction we have chosen w gp = w pg , which is a requirement for energy conservation (see the Appendix).As described by Equation (12), we added the Jacobian m j i to convert V μ from curvilinear to Cartesian coordinates, as required by the Cartesian equations of motion.
To solve for Equations ( 16)-( 17) and ( 19)-( 20), a nonlinear Newton-Krylov iterative solver is used.Convergence of this solver is typically case based and is not guaranteed when spatiotemporal resolution is insufficient to capture relevant physical phenomena.As a rule of thumb for the tests in the next sections, the temporal and spatial resolutions were constrained to satisfy cΔt/Δx < 1, with a typical value for this ratio usually set to ∼0.25.This was sufficient in all test cases to achieve absolute and relative errors below 10 −14 , typically within ∼10 iterations.Note that Δx may vary throughout our nonuniform grids, in which case the smallest value of Δx in the grid must be considered for this constraint.

Discrete Energy
To verify the exactly energy-conserving nature of this method, a proper definition for the relevant energies must be established.The electric-and magnetic-field energies are, respectively, defined as which, combined, give the total electromagnetic-field energy The particle kinetic energy is defined as The sum of the kinetic and electromagnetic-field energy is the total energy in the system: In the Appendix we show analytically that, using these definitions, the total energy will be conserved exactly (i.e., to machine precision).

Test Cases
To validate the methods developed in this paper, we considered three different test cases: the one-dimensional, electrostatic two-stream instability; the one-dimensional, electromagnetic Weibel (or filamentation) instability (Weibel 1959;Bret et al. 2010); and the two-dimensional Geospace Environmental Modeling (GEM) challenge, i.e., a paradigmatic magnetic-reconnection setup (Birn et al. 2001).

Two-stream Instability
The two-stream instability represents a classical onedimensional test for PIC codes.The initial setup consists of two counterstreaming electron beams with uniform number density n 0 /2.We initialize the particles in each beam according to a drifting Maxwellian with thermal speed v th /c = 0.001.The drift-velocity component is along the only spatial direction of the system x and is added to both beams with the opposite sign, i.e., v d /c = ± 0.2.The analytic solution of the dispersion relation predicts that the system is unstable if kv d /ω p < 1 (Goldston 1995), with k the wavenumber and ω p the plasma frequency.This results in an unstable, purely electrostatic evolution where charge bunching will generate a self-reinforcing electric-field perturbation.During the linear instability phase, the electric-field energy (Equation ( 22)) in the fastestgrowing unstable mode will grow as To test the capability of our code to handle curvilinear coordinates, we employed two nonuniform, one-dimensional grids alongside a standard Cartesian grid used as a reference.The first curvilinear grid has a small sinusoidal perturbation of one period along the domain, with mapping function where L x the length of the domain and ε parameterizes the strength of the perturbation, which we set to ε = 0.02.The second grid has a perturbation generated by a hyperbolic tangent, producing a grid with a squared hyperbolic secant profile for the cell density.A high-density peak is present near the center of the domain; the mapping function is where w determines the width of the high-density region, which we set to 4 with ε = 0.04.These two perturbed grids are constructed to ensure the total system length L x is maintained.This is necessary to ensure a fair comparison to the Cartesian reference grid.The parameters used in the sinusoidal and hyperbolic grids are chosen such that the ratio of the cell width between the smallest and largest cells is ∼0.75.This provides a meaningful deviation from the Cartesian case while still allowing the same time steps without breaking the constraint cΔt/Δx < 1 mentioned in Section 2.3.Figure 2 depicts the  25) during the linear-growth phase of the instability.The bottom left panel plots the time history of the relative error in the total energy, which confirms energy conservation to machine precision for all three cases.Lastly, the three panels on the right show a snapshot of the position-velocity phase space (x, v x ), where we observe the formation of holes that are typical for the nonlinear phase of this instability.These results give a good indication that the method is capable of faithfully reproducing electrostatic effects in a nonuniform, one-dimensional setup while maintaining exact energy conservation.

Weibel Instability
To test our method in a more general electromagnetic case, we study the Weibel (i.e., filamentation) instability.The setup is almost equivalent to the two-stream instability: N x = 2048, , and ppc = 144 uniformly distributed and divided between the two electron populations.The difference here is the direction of the beam drift velocity, which is now perpendicular to the domain's direction, i.e., v y = v d = ± 0.2c, and similarly, the thermal velocity only has a component along the perpendicular direction as well, with magnitude v th = 0.001c.The Weibel instability is excited by current separation generated by the Lorentz force in a perturbed magnetic field when acting on the two counterstreaming beams.The resulting current filaments will reinforce the initial magnetic-field perturbations.This leads to an exchange of energy between particles and magnetic fields.The magneticfield energy in the fastest-growing unstable mode will evolve as during the linear phase of the instability.
Figure 4 shows the results of the Weibel instability tests.The upper left panel shows the evolution in time of the magneticfield energy for all grids.For all cases, we observe a slight discrepancy between the simulated and theoretical values of the growth rate during the linear instability phase.However, the results of the different grids are in strong agreement with each other.The slight discrepancy with the theoretical growth rate could be attributed to different factors: for example, the finite spatial resolution prohibits resolving large wavenumbers, which would otherwise increase the observed growth rate, potentially agreeing better with the theoretical maximum growth rate that corresponds to k → ∞; in addition, noise induced by the motion of the limited number of particles drowns out small waves that might otherwise contribute to the instability.At any rate, the discrepancy is very small, and the results appear robust in terms of agreement between the different cases we run.The lower left panel shows again the evolution of the relative energy error, indicating conservation of energy to machine precision for all three grid setups.The panels on the right-hand side of the figure show the current filaments forming and growing throughout the domain as time increases.Considering these results, the method appears capable of correctly reproducing fully electromagnetic effects while conserving energy to machine precision, at least in this simple one-dimensional test case.

GEM Reconnection Challenge
The last test case presented in this paper is the GEM reconnection challenge.This is a well-known two-(or three-) dimensional setup to study magnetic-reconnection events (Birn et al. 2001).We adopt a two-dimensional xy geometry where the B x field undergoes a sign reversal along the y-axis, creating a current sheet in the inversion region.For simplicity, we adopt periodic boundary conditions, initializing two successive inversions to respect periodicity.The initial setup of the magnetic field is To respect the c∇ × B = 4πI condition, we initialize a drifting electron and positron population in the current sheets with a number density of Here B 0 and n 0 are the peak magnetic-field strength and density in the current sheet, respectively, and δ is the width of the current sheet, which we choose equal to 0.5c/ω p .At initialization, the thermal pressure of the particles in the current sheet must be able to withstand the magnetic pressure to ensure the structure neither collapses nor expands, i.e., for pressure equilibrium with k b the Boltzmann constant and T the plasma temperature.
The factor 2 is introduced to include both electron and positron contribution to the density.By equating k b T with mv th 2 , a constraint for the thermal speed in the current sheet can be found, Note that we do not manually perturb the initial equilibrium so that tearing modes eventually causing reconnection are only excited by random particle noise.
As in the one-dimensional cases, a Cartesian grid is used to establish a baseline to which we compare nonuniform-grid results.The Cartesian grid has dimensions L x = 16c/ω p by L y = 32c/ω p with N x = 160 by N y = 320 cells and ppc = 16 particles per cell, which provides a sufficient signal-to-noise ratio for this setup.The number of particles required in a setup will always depend on the specific subject and requirements of the simulation.The time step was set to w D = t 0.025 p 1 for nt = 30,000 steps, for a total time of w - 750 p 1 .This satisfies the convergence constraint cΔt/Δx < 1.
To test the behavior of our code with curvilinear coordinates, we construct a nonuniform grid such that we achieve the same density of computational cells (i.e., the same numerical resolution) near the current sheets while significantly reducing the grid resolution in the upstream regions.This naturally results in fewer total grid cells than in the uniform Cartesian case.This grid setup enables us to resolve the reconnection region with the same precision of the Cartesian setup while greatly reducing the total number of cells in regions of no physical interest.The mapping function for this grid was created using piecewise linear functions that were connected using hyperbolic tangents as an approximation for the Heaviside step function, )).Here, r n is the cell-density ratio with respect to the low-density upstream regions.This is set to 1 in the upstream regions and to a desired larger ratio for the current-sheet regions.s determines the sharpness of the transition, with larger values more closely resembling the Heaviside function and lower values creating a smoother transition region.This value was set to s = 5 for all grids.The b n terms determine the break points, i.e., the ycoordinate where the piecewise function switches from one linear function to the next.These break points therefore determine the points of transition between high-and lowdensity regions.These are set at L y /8, 3L y /8, 5L y /8, and 7L y /8 respectively.The p n factors are the offsets of the individual linear functions, which are fully constrained by the previous parameters.Finally, p max is the maximum value reached by the piecewise function and is introduced such that η ä [0, L y ].

g y r y p r y p s y b tanh
We employed two nonuniform grids, the first with a density ratio r = 5 and number of cells in the y-direction N y = 192, and the second with ratio r = 10 and N y = 176.The number of cells in the x-direction remains unchanged from the Cartesian case (N x = 160).By design, the curvilinear grid has the same physical size L x = 16c/ω p and L y = 32c/ω p of the Cartesian grid.Furthermore, the time step and number of iterations were deliberately set to be identical in both cases to again satisfy the convergence constrained since the smallest Δx is also unchanged.The number of particles was also unchanged to maintain the same amount of particles per cell and thus a similar signal-to-noise ratio, in the high-grid-density region.Since the Python test code is computationally limited by the number of particles, this means the observed run time was very similar between the three cases (Cartesian, curvilinear with r = 5, and curvilinear with r = 10).However, it is at the moment not possible to assess the computational efficiency of the method based on our test implementation since it lacks lowlevel optimization and parallelization commonly found in production-ready codes.
Our results are shown in Figure 5.The left panel shows the Cartesian grid in the top half of the domain and the curvilinear grid with density ratio r = 10 in the bottom half.In the middle panel, we show a comparison of the out-of-plane current J log z 10 | | and the in-plane magnetic-field lines between the Cartesian grid (top) and the curvilinear grid with density ratio r = 10 (bottom) at time w = t 750 p 1 .Although it is clear that the upstream regions have considerably lower resolution in the curvilinear case, the reconnection regions of interest look very similar.Since the exact position of the formed X-points and magnetic islands is due to the random initialization of the particle properties, they appear in different places in the top and bottom current sheets in the domain.In the top right panel, we plot total magnetic and kinetic energies for all cases to show the energy exchange between magnetic fields and particles over time.The initial magnetic field contains the majority of the energy, and as the magnetic fields reconfigure into a more relaxed state via reconnection, the particles gain energy and get accelerated.The bottom right panel shows the evolution in time of the relative error in the total energy, which is again conserved to machine precision, supporting our findings from the one-dimensional tests and the theory in the Appendix.

Conclusions
In this paper, a method is presented to adapt a fully implicit-PIC method to work with arbitrary coordinate systems, including curvilinear coordinates, while maintaining exact (i.e., to machine precision) energy conservation.While previous methods have been published for reduced equations, like the Darwin approximation (Chacón & Chen 2016), this is the first time to our knowledge when this is achieved for a PIC with the full Maxwell equations.This is achieved by changing the coordinates to a system where the grid becomes logically rectangular.To change coordinates, we introduce the metric tensor, Jacobian matrix, and contravariant curl operator into Maxwell's equations for the PIC algorithm.The particle equations of motion are solved in Cartesian coordinates to avoid the complexity of calculating the metric and Jacobian at each particle position and at each iteration.
To validate the method, we have presented three different test cases.The two-stream instability and Weibel instability were tested in one spatial dimension, and the GEM reconnection challenge was tested in two dimensions.The tests run with nonuniform (curvilinear) grids were compared to the results obtained with a uniform Cartesian grid, taken as a baseline.Two nonuniform grids were tested in the one-dimensional cases, both generated by perturbing the baseline Cartesian grid: the first with a sinusoidal perturbation and the second with a hyperbolic-tangent perturbation.The nonuniform grid in the GEM challenge was designed to greatly reduce the total number of grid cells while maintaining the same spatial resolution of the uniform grid around the current sheets to properly resolve reconnection dynamics while underresolving the unimportant upstream regions.Comparisons were performed on global quantities such as individual energy components and the total energy of the system.
We demonstrated that our new method achieves energy conservation to machine precision in all three test cases with both Cartesian and curvilinear grids.The two-stream instability modeled with curvilinear grids shows a growth rate of the electric-field energy, which closely matches both the (Cartesian uniform) baseline and the theoretical growth rate.The Weibel instability also shows excellent alignment of the magnetic-field energy to the Cartesian baseline.In our nonuniform-grid runs, we observed that the two-dimensional GEM challenge produced reconnection events inside the current-sheet region as expected of this setup, despite the nonuniform grid having only 55%-60% the total number of cells with respect to the uniform Cartesian run.Reconnection can be observed in the inplane magnetic-field line rearrangement, which shows the formation of magnetic islands, and in the exchange between magnetic and particle energy, showing a conversion of magnetic-field energy into kinetic energy as particles get accelerated during the reconnection process.
Although our method possesses desirable properties (exact energy conservation and grid adaptability to complex geometries), ample ground is left for further improvements, which we will pursue in future work.Unlike the energy, the total charge is not conserved since ∇ • E = 4πρ is not enforced in our implementation.However, as shown in other works (e.g., Chen & Tóth 2019), a charge-conserving adaptation of fully implicit-PIC methods is possible (Chen et al. 2011) and could be implemented in our code in the future.Another shortcoming is represented by the convergence of the nonlinear iteration, on which our method is based, which is not guaranteed; from our experiments, we observed that a sufficiently high spatiotemporal resolution must be enforced to avoid nonconvergence issues, which could break the temporal iteration.In many practical cases, however, this will not pose a problem since a fine spatiotemporal resolution is required anyway to resolve the physics of interest.In our tests, convergence (i.e., an absolute error below a tolerance ∼10 −14 ) was always attained with ∼10 nonlinear iterations; hence, the computational cost in our example simulations was approximately a factor 10 larger than a standard explicit PIC method since each nonlinear iteration roughly corresponds to one explicit time cycle.The greater computational cost of our method is, however, counterbalanced by numerous advantages, as discussed below.In addition, speeding up convergence is, in principle, possible with more advanced strategies: for example, Chen et al. (2014) have shown that by preconditioning the nonlinear solver, the average number of iterations can be drastically reduced and the spatiotemporal resolution constraints can be relaxed.
By combining a curvilinear grid with an exactly energyconserving PIC method, we have achieved a unique set of features with numerous potential applications.In general, nonuniform grids are especially well suited for setups with a small region of interest within a much larger domain or setups where the bulk of the plasma is contained within a spatial region with a nontrivial shape.Meanwhile, energy conservation is important in many problems where energy conversion from one type to another occurs and has significant impact on the plasma dynamics, e.g., reconnection events or instabilities (see Markidis & Lapenta 2011 and references therein).In several cases, the combination of energy conservation and custom grids can be of high interest: in fusion plasma physics, our method could be used in simulations of the scrape-off layer (i.e., the plasma layer close to the separatrix) of a tokamak plasma (e.g., Ohtsuka et al. 1978;Xu & Cohen 1998;Krasheninnikov 2001;Fundamenski et al. 2007).The nonuniform grid can be tailored to match the geometry of the boundary layer, greatly reducing the total computational cost.Another example is the expanding solar wind (e.g., Innocenti et al. 2019Innocenti et al. , 2020;;Bott et al. 2021;Micera et al. 2021): here, a large domain, which expands in time as the solar wind travels, makes it difficult to employ a standard Cartesian box while maintaining sufficient resolution throughout the domain.By using a grid that mimics the shape of the expanding solar wind, with a high resolution closer to the Sun and decreasing grid density farther out (where less resolution is required), the computational cost of this problem can be significantly reduced.Conversely, when considering compression-driven dynamics in plasmas, a compressing box has been used in the past in the context of simulations of the intergalactic medium (e.g., Sironi & Narayan 2015), which could similarly be replaced with our nonuniform grids.Other solar wind-related phenomena whose modeling could benefit from a custom grid are switchbacks.Switchbacks are magnetic structures in the solar wind characterized by a typical S shape of the plasma stream and may be related to interchange reconnection events and particle scattering and acceleration (e.g., Drake et al. 2021;Wyper et al. 2022).With our method, a curved grid can be fitted to this S-shaped stream, reducing the total grid resolution and cutting out the surrounding plasma that might be of little interest.This might open up the possibility to a fully kinetic study of the driving phenomena, which was thus far unfeasible.In higher-energy plasma scenarios, our method could also find applicability, e.g., for general-relativistic plasma simulations.Recent works have simulated plasma accretion around compact objects such as black holes (e.g., Parfrey et al. 2019;Bransgrove et al. 2021;Crinquand et al. 2022;El Mellah et al. 2022;Galishnikova et al. 2023), where plasma flows in structures such as accretion disks and jets.These works present large-scale simulations utilizing explicit codes, usually implemented with one specific (four-dimensional) metric corresponding to the particular physical case in question.With the method presented in this work (modified by adding a time component of the metric and relativistic effects; see, e.g., Bacchini et al. 2019;Bacchini 2023), we can simultaneously avoid the aforementioned downsides of explicit codes (i.e., numerical instability and lack of energy conservation) and generalize such simulations to any arbitrary metric tensor that might be of interest, opening the possibility for more physical cases to be studied.
In all potential applications, it is important to keep in mind that the curvilinear grids can potentially influence the transmission of electromagnetic waves and can do so nonuniformly and anisotropically if the grid resolution changes dramatically across the domain.High-frequency waves will be dampened in regions with insufficient resolution along the wave direction of propagation.Therefore, if such waves have a notable impact in the simulation, one should ensure sufficient grid resolution in the appropriate regions in order to capture their propagation.
We foresee several potential developments for the future.The logical next step in this line of work is to implement the method in a production-ready, parallel infrastructure such as those currently employed for implicit-PIC simulations (e.g., IPIC3D, ECSIM).This will open the door to larger-scale simulations outside the capabilities of our current test implementation in Python.The most obvious candidate would be the ECSIM code (Lapenta 2017) since it shares (by design) a very similar discretization of our governing PIC equations.This step will be carried out in the future.It would also be of great interest to consider novel architectures such as GPUs.The traditional bottlenecks in PIC codes are the particle gathering and deposition steps.With the increasing popularity of GPUaccelerated codes, these bottlenecks could be mitigated since the gathering and deposition steps are well suited to parallelization (e.g., Joseph et al. 2011;Decyk & Singh 2014;Vasileska et al. 2021).The field solver, on the other hand, whether implemented on CPUs or GPUs, requires global communication of grid quantities, which hinders scaling behavior to a large number of cores, so it becomes imperative to further reduce the cost of the field-solution step through other means.With the method described in this paper, this can be achieved by optimizing the grid and thus reducing the total number of cells, which is the main driving factor of the field solver's computational cost.
Subtracting Equation (A1) from Equation (A2) and rearranging terms gives On the left-hand side, the time derivatives are discretized as described in Section 2.3, and half-integer time step values are substituted with a temporal linear interpolation, producing Summing over all grid points and multiplying with J g ΔξΔη/4π while working out the products on the left-hand side results in On the left-hand side, we can now recognize the definition of the field energies as described by Equations ( 22)-( 23).The first term on the right-hand side is the divergence of the Poyntingflux vector; this discrete integral (expressed by the sum over all grid points g) vanishes under periodic boundary conditions.This can be easily shown for the one-dimensional case by expanding the sum, where every term can indeed be canceled out if the zeroth and Nth grid points are identical.This argument can easily be extended to two and three dimensions as well.The right-hand side in Equation (A5) thus reduces to the second term and can be used to denote the change in field energy between two subsequent time levels, where the remaining term represents energy exchange between fields and particles.
We now have to verify that this term is identical to the energy-exchange term obtained from the particle equations.To do so, we start from Equation (20) and contract both sides with v i p In the last term on the right-hand side, we can recognize v i ò ijk v j B k = 0 since this is the dot product of v with the cross product between itself and B. Like before, we rewrite the halfinteger time steps as temporal averages in the left-hand side, and we take the sum over all particles p to get A 9 We can recognize the definition of the kinetic energy as shown in Equation (24) and use the expression for interpolated quantities (Equation ( 21)) in the right-hand side: Since the interpolation weight from grid to particles w gp is the same as the weight from particles to grid w pg , it is possible to switch the order of the summations, and by adding the missing factors we can recover the righthand side of Equation (18) for the current density, å å Since this vanishes exactly between any two time steps, the method conserves the total energy U exactly (to machine precision) independently of the time step, grid spacing, and any other numerical parameter.

Figure 1 .
Figure1.Diagram of the coordinate transformation between physical and logical spaces.The Jacobian and its inverse are used to transform contravariant vector components.In physical space the Euclidean metric is used, while a nontrivial metric is used in logical space.

Figure 2 .
Figure 2. The three one-dimensional grids used in the study of the two-stream and Weibel instabilities: from top to bottom, Cartesian (i.e., the usual grid employed in standard PIC), sinusoidal, and hyperbolic-tangent mapping functions.They are shown here for a 64-cell grid with the gray scale indicating the local grid density.

Figure 3 .
Figure 3. Top left: evolution of the electric-field energy for the different grids as well as the theoretical growth rate associated with the linear-growth phase of the twostream instability.Bottom left: evolution of the relative error in the total energy for the different grids.Right: snapshots of the (x, v x ) phase space during the nonlinear stage for the different grids.

Figure 4 .
Figure 4. Top left: evolution of the magnetic-field energy for the different grids as well as the theoretical growth rate associated with the linear-growth phase of the Weibel instability.Bottom left: evolution of the relative error in the total energy for the different grids.Right: formation and growth of current filaments in the domain with evolving time for the different grids.

Figure 5 .
Figure 5. Left: comparison of the uniform Cartesian grid (top) and curvilinear grid with r = 10 (bottom), which has a reduced resolution in the upstream regions, while maintaining the same resolution in the current sheets that was used in the Cartesian grid.The magnetic islands formed in the top and bottom halves are not collocated since their positions are randomized by the starting values of the particle properties.Middle: comparison of a snapshot of the out-of-plane current and the in-plane magnetic-field lines for the aforementioned grids.Top right: evolution of the magnetic-field energy and kinetic energy for the different grids.Bottom right: evolution of the relative error in the total energy for the different grids.
-exchange term with the field energy is recovered.Combining the change in kinetic energy (Equation (A14)) and field energy (Equation (A7)) gives the change in total energy of the system between time steps,