Moment tracking and their coordinate transformations for macroparticles with an application to plasmas around black holes

Particle-in-cell (PIC) codes usually represent large groups of particles as a single macroparticle. These codes are computationally efficient but lose information about the internal structure of the macroparticle. To improve the accuracy of these codes, this work presents a method in which, as well as tracking the macroparticle, the moments of the macroparticle are also tracked. Although the equations needed to track these moments are known, the coordinate transformations for moments where the space and time coordinates are mixed cannot be calculated using the standard method for representing moments. These coordinate transformations are important in astrophysical plasma, where there is no preferred coordinate system. This work uses the language of Schwartz distributions to calculate the coordinate transformations of moments. Both the moment tracking and coordinate transformation equations are tested by modelling the motion of uncharged particles in a circular orbit around a black hole in both Schwarzschild and Kruskal–Szekeres coordinates. Numerical testing shows that the error in tracking moments is small, and scales quadratically. This error can be improved by including higher order moments. By choosing an appropriate method for using these moments to deposit the charge back onto the grid, a full PIC code can be developed.


I. INTRODUCTION
In numerical simulations involving the dynamics of a large number of particles, for example, in a plasma, it is impossible to track the dynamics of every ion and electron.In particle-in-cell (PIC) codes, the plasma is modelled using macroparticles, where each simulated macroparticle represents a large number of actual particles.To improve the accuracy of a PIC code, there are two options: use more macroparticles, or give each macroparticle more information.A macroparticle typically has only a position and a velocity.This paper presents a different method, in which a macroparticle represents the moments of a group of particles [1][2][3][4][5] (figure 1).Such a model may be more efficient in cases where a large number of particles can be accurately modelled by only a small number of macroparticles and their moments, and the electromagnetic field does not vary much across the extent of the macroparticle.Additionally, in cases where the dominating interactive force can be calculated from the Liénard-Wiechert potential, the electromagnetic fields can be calculated directly from the moments of the macroparticle [6].
There are several existing methods for tracking moments: the code MERLIN implements a transition matrix approach for particle accelerators [4] and a continuous model has been developed through a Hamiltonian approach [5].Moment tracking can also be done continuously by differentiating the definition of a moment and using the Vlasov equation [3,7].In plasma physics, the concept of moment tracking is often used to describe hybrid-Vlasov approaches, where individual species of the plasma may be modelled through their bulk properties [8][9][10].Moments in the context of hybrid-Vlasov approaches are constructed by integrating over velocity space only and can be interpreted as physical quantities such as temperature and pressure.In this work, moments shall be FIG.1: Tracking several individual particles compared to a macroparticle with moments.The 5 particles (the black lines) in the left diagram are replaced by a single macroparticle (the orange line) in the right diagram.By tracking the moments, quantities such as the difference between the centre of charge of the 5 particles and the position of the macroparticle (the first order moment, represented by the horizontal blue arrows), and the variance in position of the particles (the second order moment, represented by the green error bars) can be tracked.
constructed by integrating over both velocity and position space, giving related, but different quantities.The transport of moments is also used in wider fields where the Liouville equation holds, such as particle nucleation [11], crystal growth [12], and nuclear collisions [13].
To perfectly track the moments, an infinite set of moments is required.This is because, as is shown in this paper, higher order moments generate lower order moments.Tracking more moments is more computationally expensive, as these moments come with additional computational work (both more differential equations to solve each time step and more memory usage per macroparticle).This increased computational work is balanced by reducing the total number of macroparticles in the simulation.At the quadrupole

Order of moments tracked
Information known

Number of differential equations
to solve relative to a standard macro-particle Monopole (zeroth order) X µ ,U µ 6 1 Dipole (first order) X µ ,U µ , V a 12 2 Quadrupole (second order) X µ ,U µ , V a , V ab 33 5.5 Octopole (third order) X µ ,U µ , V a , V ab , V abc 89 14.83 Hexadecapole (fourth order) X µ ,U µ , V a , V ab , V abc , V abcd 215 35.83 TABLE I: The amount of information in a macroparticle that also tracks moments compared to a standard macroparticle.This increased information results in more memory usage and a larger number of first order differential equations to solve.Since there are 6 differential equations per time step for a standard macroparticle, the fourth column is the third column divided by six.
(second) order there are 33 equations to solve each time step (21 quadrupole moments, 6 dipole (first order) moments, 3 components of velocity and 3 components of position).This expands to 215 equations if the expansion is carried out to the hexadecapole level (fourth order) (table I).To make moment tracking computationally feasible, a truncation is required.This truncation is the highest order of moments considered, above which the contribution from higher order moments are neglected.Care must be taken to ensure the truncation is of as low an order as possible to minimise computational load, whilst also ensuring that neglecting the higher order moments does not significantly impact accuracy.This article calculates all quantities to quadrupole order, although all results presented can be generalised to arbitrary order.
Recently there has been a focus on using PIC codes to model the dynamics of plasma around black holes [14][15][16][17].Such plasmas may be important in active galactic nuclei, pulsars and gamma-ray bursts [17][18][19].Ref. [20] contains a full review of these studies.In such systems, there is not a preferred choice of coordinate system i.e. when modelling a static uncharged black hole, there is a choice to work in Schwarzschild coordinates, or Kruskal-Szekeres coordinates, amongst others.Because there is a choice in coordinate systems, it is useful to be able to transform between different coordinate systems, especially where the time and space coordinates are mixed together (such as the coordinate transform between Kruskal-Szekeres and Schwarzschild coordinates, shown in figure 2).Coordinate transformations that mix space and time coordinates also appear in particle accelerators.When simulating the motion of a linearly accelerating bunch in a particle accelerator, the transformation into the instantaneous rest frame of an accelerating bunch mixes space and time coordinates (this transformation is similar to the one shown in figure 2), so the spacetime coordinate transformations presented in this work are necessary.
The moments of a macroparticle depend on the choice of time slicing (figure 3).In all coordinate systems considered in this article, the global coordinate time will be used as the time slicing.This time slicing is a foliation given by spatial hypersurfaces of constant coordinate time.In general, different coordinate systems will give different time slicings.This means when transforming between coordinate systems that mix temporal and spatial coordinates, for example, transforming between Schwarzschild and Kruskal-Szekeres coordinates, the time slicing will change (figure 2).This time slicing may be the global Killing timelike vector in relativistic scenarios, or the lab time in particle accelerators.Another possible time slicing is the backward light cone of an observer, which is the frame used when making astrophysical observations.A choice of time slicing commonly used in modelling plasma around black holes is found in the fiducial observer (FIDO) scheme, where the time slicing is given relative to local moving observers [21,22].A fourth possible time slicing is to take all the vectors orthogonal to the velocity of the world line.By using geodesics to propagate these to the world line, these can be used for the Dixon representation of a multipole [23].
One of the key results of this article is the formula for the coordinate transformation of moments between coordinate systems that mix space and time coordinates.The standard integral representation of moments cannot transform moments between coordinate systems where the time slicing changes.Simply transforming the moments into the new coordinate system does not take into account the change in time slicing.If one were to use the standard representation of moments, an 'effective' coordinate transformation can be performed.This is the process of reconstructing the distribution function from the moments, transforming the distribution function into the new coordinate system and retaking the moments in the new coordinate system.In order to find the coordinate transformations that change the time slicing, a different representation of moments must be used, in terms of Schwartz distributions [24,25].In this representation of moments, the moments are transformed into the new coordinate system, then projected onto the new time slicing to find the full coordinate transformation.The dashed black line is the event horizon.The moments in a given coordinate system have no components in the direction of the time slicing in that coordinate system.FIG.3: An example of two different time slicings.When taking moments, the time slicing is the hypersurface the moments are integrated over.By using a different time slicing (e.g.either the horizontal dashed purple lines t or the curved dashed blue lines t), the distance between the centre of the macroparticle (the orange line) and the nearby particles the macroparticle is representing (the green contours) is different, so different moments will be found.
In a standard particle-in-cell code there are three stages: the particle pusher, charge and current deposition, and the field updater.This work will focus only on the particle pusher stage of the PIC code, showing how tracking the moments of the group of particles represented by the macroparticle gives extra structure to the information within the macroparticle.There will be no consideration of particle-particle interactions.Particles will only be affected by external fields.This approach means it is not meaningful to talk about the cells at this stage.Future work will consider how to deposit the charge and current associated with the extra structure of the moments onto the grid.This future work will examine the appropriate cell size needed for this deposition.This article calculates the moment tracking equations and coordinate transformations for arbitrary coordinate systems.The structure of this article is as follows: Section II introduces the Vlasov equation, which describes the dynamics of collisionless charged particles.Section III introduces moments, and the existing methods for modelling them.Section IV introduces a distributional representation of moments in terms of derivatives of Dirac delta functions and uses this representation to show the time evolution of moments.Section V finds the coordinate transformations for the moments, which is the main result of this work.Section VI shows a computational model of uncharged particles circling a black hole to validate the theory at a proof-of-concept level.Section VII uses the language of differential geometry and de Rham currents to present a geometric interpretation of the moment differential equations and finds the coordinate transformations of the moments through this language.Section VIII discusses the considerations that will need to be made in the development of a method that uses the moments to deposit charge and current onto the grid, in order to develop a full PIC code.We conclude in section IX by discussing other interesting direction of future investigation.The appendix contains the more technical proofs needed to establish the theory.

A. Notation of indices
To study the dynamics of plasma, calculations are performed in 7 dimensional phase space.Because of this, several different summation conventions for indices are needed: summations over just space, just velocity, space and velocity, time and space but not velocity, and space, time and velocity.In order to account for all these different summations, this article uses the convention that Latin indices a, b, c represent summations over 0, . . ., 6; Greek indices µ, ν, ρ represent summations from 0, 1, 2, 3; and an underlined index means there is no summation over the 0 index, i.e. a = 1, . . ., 6 and µ = 1, 2, 3.The coordinates used are (t, x µ , u µ ), where t is the global time, x µ is a space coordinate, and u µ is a spatial component of the 4-velocity.The notation (ξ a ) will be used to represent a general coordinate, such that i.e. ξ 4 = u 1 .Dimensions are used such that c = G = 1.The notation f | p will be used to represent the evaluation of a function at a point, i.e. f | p is f evaluated at p.After their introduction, the arguments of functions will not be written, unless needed for emphasis.

II. THE VLASOV EQUATION A. 7 Dimensional phase space
To calculated the complex dynamics of plasmas, calculations are performed in 7 dimensional phase space.These dimensions are time, 3 spatial dimensions and 3 proper velocity (4-velocity) dimensions, where t is the global coordinate time.As previously stated, these will be represented in coordinates as (t, x µ , u µ ).The time component of 4-velocity, u 0 (t, x µ , u µ ), is a function in phase space, defined as the positive solution to the quadratic equation where g(t, x µ ) is the metric with signature (−, +, +, +).In the specific example of Minkowski spacetime with Cartesian coordinates, u 0 is the Lorentz factor in special relativity, given by

B. The Vlasov equation
Consider the dynamics of a group of charged particles, each with charge q and mass m with distribution function f (t, x µ , u µ ).By assuming these particles are collisionless, they can be modelled by the Vlasov equation.To find the Vlasov equation, consider an arbitrary spacetime with metric g µν (t, x µ ), Christoffel symbols Γ µ νρ (t, x µ ) and electromagnetic 2-form F µν (t, x µ ).Let C(t) be a world line parameterised by t, the same parameter as the global time ξ 0 .Let η(t) be the prolongation of C(t), that is, a curve such that For a given value of t, η(t) is a point in 7-dimensional phase space.In terms of coordinate functions (t, x µ , u µ ), this can be represented as In general, η(t) is a curve in position-velocity phase space.In this work η(t) will be the curve moments are taken around.In particle accelerators, the ideal orbit is a natural choice for η(t).This choice does not exist in plasmas; in these cases choices for η(t) could be the position of the macroparticle, or a trajectory based on the initial centre of charge of the macroparticle.
The Vlasov vector field, W , is the integral curves of η, defined as Finding these derivatives gives 7 ODEs, The acceleration can be found using the pregeodesic equation combined with the (pre-)Lorentz force, where the 1/u 0 term in the Lorentz force arises when coordinate time time is the parameterisation, rather than proper time, and κ(t) arises from the parameterisation.κ(t) can be found by noting C 0 = t, hence d 2 C 0 /dt 2 = 0, solving this gives From this, solving equation ( 7) for all η gives The Vlasov equation describes the motion of the whole particle distribution, where This article will work exclusively in coordinate time frames: these are frames where the parameterisation of η is the same as the time coordinate.In any coordinate time frame, the Vlasov vector field has the form of equation (10).This formulation of the Vlasov vector field is distinct from the 3 + 1 formalism used in other approaches for simulating plasma in general relativity [21,22], although the moment tracking and coordinate transformations presented in this work can be used in both formalisms.

III. MOMENTS
The collective properties of a group of charged particles with distribution function f can be modelled by their moments.This article makes the choice that f is a scalar density of weight 1.This choice means that the moments are affected by the measure from a curved spacetime through a change in f .The zeroth order moment, known as the monopole, is defined as and corresponds to the total charge of the macroparticle.The first order moment, known as the dipole moment, is given by These can be combined into a single equation, where ξ a and η a are defined by ( 1) and ( 5) respectively.The dipole moment corresponds to the deviation of the centre of charge from the centre of the macroparticle.If the initial centre of the macroparticle is chosen such that the dipole moments are initially zero, then whilst the centre of the macroparticle will obey the Lorentz force equation, the centre of charge will not [26].Hence, the two paths will diverge.The dipole moment represents the difference between these.
This formalism can be generalised to the nth order moment, where V is totally symmetric and n! is required for counting due to the symmetry.Note the sums are only over space and velocity coordinates, and there are no corresponding 'time' moments.The moments in this work are integrated over both position and momentum space, in contrast to the conventional moment equations for plasmas, which are just integrated over velocity space [8].The naming convention for multipoles scales as 2 k , so the zeroth order moment is the monopole, then the dipole, quadrupole, octopole etc. (table I).The quadrupole moment corresponds to the variance of the macroparticle; the octopole moment corresponds to the skew of the macroparticle, and the hexadecapole moment corresponds to the kurtosis of the macroparticle.

A. Dynamics of the moments
When parameterising with a global time, the dynamics of the moments can be calculated by differentiating (15) with respect to time [1,2], ) By using the Vlasov equation and Taylor expanding around η, the dynamics can be found.As an example, the differential equation for the dipole is This differential equation is an infinite series of higher order moments, and a choice must be made to truncate this series at some point.Truncating at quadrupole (second) order, the differential equations for the dipole and quadrupole are given by The dipole generated from the quadrupole can clearly be seen.
All moment tracking methods suffer from this error -higher order moments generate lower order moments.The importance of the choice of truncation on the accuracy of the model is discussed in section VI.

B. Coordinate transformations of moments where the time slicing is preserved
As moments are defined through integrals, they are highly coordinate dependent objects.For coordinate transformations where the time coordinate is unchanged, this definition of the moments can be used to find the coordinate transformation.Consider a new coordinate system denoted by hatted coordinates (t, ξ â) where t = t (since the time coordinate stays the same) and use hatted indices to represent a summation in the new coordinate system.The distribution function f transforms as a scalar density [27] of weight 1, i.e.
where a hatted variable represents that variable in the new coordinate system, and d 7 ξ /d 7 ξ is the Jacobian.By expanding ξ in terms of ξ around η, and noting η(t) = η(t), the coordinate transformations can be found.These are infinite sums in terms of all higher order moments.Truncating the expansion at quadrupole order, the coordinate transformations up to the quadrupole order for the dipole and quadrupole are given by For coordinate transformations that mix the space and time coordinates, this method for coordinate transformations will not work.The difficulty in transforming coordinates using this method arises from the important dependence on the choice of time.In order to find the moments, the coordinate systems defines a time slicing, which is a hypersurface of constant t (figure 3).The integrals in the definition for the moments in equation ( 15) are over these hypersurfaces.When performing a coordinate transformation that mixes the time and space coordinates, then the time slicing will also change.Changing the time slicing will give different moments.The integral representation of moments in equation ( 15) cannot be transformed between coordinate systems that mix time and space coordinates.In this article we present a new representation of moments using Schwartz distributions.These distributions are introduced in section IV, and the use of this language to find the coordinate transformations for the moments is presented in section V.

IV. SCHWARTZ DISTRIBUTIONAL MULTIPOLES AND THEIR TIME EVOLUTION A. Ellis multipoles
Rather than defining moments explicitly through spatial integrals (equation ( 15)), an alternative representation of moments is through derivatives of Dirac delta functions.This means that rather than modelling a macroparticle as just a Dirac delta function and using a shape function to deposit charge (the cloud-in-cell approach), the macroparticle is represented as a set of derivatives of Dirac delta functions.By using the language of the Schwartz distributions presented in this section, any coordinate transformation of the moments can be calculated, including those mixing space and time coordinates, which, as previously stated in section III B, cannot be done through equation (15).By defining a moment through the method presented below, there is no dependence on the time slicing in the definition of the moment, and as such, the spacetime coordinate transformation can be found.Representing a multipole expansion using derivatives of a Dirac delta function is known as the Ellis representation of a multipole [6,25,28].It may also be known as the Schwartz distribution representation of a multipole.In terms of Dirac delta functions the expansion to the second order (the quadrupole) is where and ) where d 6 ξ = d 3 xd 3 u.Note q has no dependence on time due to the conservation of charge.
The terms V a , X a ,V ab and X ab are called the components of a multipole.As this definition involves derivatives of Dirac delta functions, they are defined by their action on test functions (φ 0 , . . . ,φ 6 ) which are the components of a covector, and have compact support.Equation ( 22) acting on φ a gives The evaluation of the test form at η will not be explicitly written in future representations of J a , but is implicitly present.
The components in (24) are only over (1, . . ., 6) (there are no terms of the form V 0 or V 01 ).By writing a multipole in this way the components of a multipole are unique.The components can be extracted by acting J a on specific test forms.The details of these test forms and how they isolate the components of J a are outlined in appendix A 1. A discussion on why this means they are unique can be found in ref. [25].

B. Relating the components of Ellis multipoles and moments
To show the components of ( 25) are the moments of f (equation ( 24)), the distribution function f is squeezed.This is the process of representing a function using a Dirac delta function and the derivative of Dirac delta functions (figure 4).The moments of f are the coefficients of this expansion.By doing this, it can be shown that the moments of the distribution function f naturally appear in the components of the Ellis representation of a multipole.
Squeezing a distribution can be shown using the language of differential geometry, the process of which can be found in appendix A 2.

C. Time evolution of moments
The evolution of the moments defined by (22) are governed by two conditions: the conservation of charge, described by FIG.4: By squeezing the width of a function down whilst increasing its height, such that its area stays the same, a function can be represented by a series of derivatives of Dirac delta functions.The coefficients of this series are the moments of the function.
the equation and the Vlasov equation, The combination of the equations ( 26) and ( 27) are called the transport equations, and are used to find the evolution of the moments in the Ellis representation.The transport equations can also be defined by their actions on test forms λ and α ab , where α ab is antisymmetric, and both λ and α ab have compact support.These two conditions give Equation ( 18) can be found as a trivial corollary of this, so this method can also be used to find the differential equations for the moments.These differential equations are the same as directly differentiating (15).
Proof.The proof of this is in appendix A 3.
Whilst both methods achieve the same result, there is a different philosophy to each approach.The method presented in section III A gives an infinite Taylor expansion, and picks the truncation point at the end of the method.Alternatively, the Ellis representation makes the choice of truncation when first performing the multipole expansion to a specific order.This order of expansion is then kept throughout.Whilst this work only carries out this expansion to quadrupole order, it is trivial to extend this result to higher orders.

V. COORDINATE TRANSFORMATIONS OF MULTIPOLES WHERE THE TIME SLICING CHANGES
Having introduced the Schwartz distributional multipole we are now in a position to calculate the coordinate transformations where the time slicing changes.Consider a new coordinate system denoted by hatted coordinates (t, ξ â), where the time coordinate is changed, and use hatted indices to represent a summation in the new coordinate system.To find the transformation rules of a multipole, it is assumed J a transforms as a density, i.e.
such that is an invariant quantity.
Since ηa are functions on a world line, they transform as where Recall φ a and ∂ a are tensorial so the transformation rules for these are Note that if summations just run over (1, . . .6) in the original coordinate system, in general the indices in the new coordinate system will run over (0, . . .6).Using this, it can be shown that using these transformation rules gives where and a hatted index represents a sum over coordinates in the new coordinate system.
Proof.The proof of this is in appendix A 4.
Observe that U a ,U ab ,Y a and Y ab have indices ranging over (0, . . ., 6), not just (1, . . ., 6).Multipoles with components that have indices running over (0 . . .6) do not have unique components.As shown in section IV A, for the components of a multipole to be unique, the components must only range over (1 . . .6).The reason the multipole contains terms of the form U 0 etc. is because it is still adapted to the time slicing t from the original coordinate system.To adapt the components to the new time slicing, and thus find the full coordinate transformation, they are projected onto the new time coordinate t.Consider differentiating along a world line, By rearranging this for ∂ 0 , these terms are projected onto com- ponents along ∂ a , and components along the world line, Since U a are functions along a world line, with similar relations for U ab ,Y a ,Y ab , and φ | η .By using this projection, terms like the U 00 term get projected into terms in V ab , Xab , V a , and Xa .
Using this projection allows multipoles to be adapted to the time slicing t in the new coordinate system, i.e. adapted such that the components indices only range over (1 . . .6), giving the full coordinate transformation.The new moments are where Similar terms exist for X ab and X a such that ( 29) is still satisfied, and are shown in equation (A67) in appendix A 5. They can also be found using (24) in the new coordinate system.
Proof.The proof of this is in appendix A 5.
Equation (40) and equation (41) are a new result, and are numerically tested in the next section to show their validity.In the case there is no change in time coordinate, A 0 b = 0, this reduces to equation (21).This work has transformed between two coordinate time frames, but can be generalised to any frame, not just one where η0 = 1, with a more general projection VI. COMPUTATIONAL VALIDATION To test the accuracy of the moment tracking and coordinate transformation equations, a computational model was developed.The code tests if the truncation to quadrupole order is acceptable, or if a higher order multipole expansion should be used for practical cases.The results presented here are an example of a coordinate transformation that mixes time and space coordinates, focusing on the particularly challenging case of black holes.This work is not solely limited to plasma around black holes.In particular, the moment tracking can be applied to any plasma.
To develop the code, the derivatives of the Vlasov equation were calculated using the symbolic algebra software MAPLE.The simulation itself was written in C++.There are 36 first derivatives of the Vlasov equation and 126 second derivatives.Each of these derivatives is very complex, hence the need to calculate them using symbolic algebra software.To give an example of this complexity, in Schwarzschild spacetime, one derivative of the Vlasov field is The model uses the forward Euler method to integrate both the particle motion and the moment equations.Particles are updated using the equations where u 0 is defined by equation (2) and ∆t is the time step size.Moments are updated using the equations Despite the high numerical error associated with the forward Euler method, it is acceptable for use as a test to assess the validity of the equations, as the dominating error is not numerical.
In all cases, only uncharged particles will be tracked.The analytical results calculated in the previous sections can be applied to charged particles.To calculate the inter-macroparticle forces needed for a full PIC code, a method to use the moments to deposit the charge from a macroparticle onto the grid must be developed.This deposition process is beyond the scope of this paper.An outline of a method that can be used to deposit the charge and current onto the grid is given in section VIII.

A. The spacetimes modelled
To validate the model, tests were performed in both Schwarzschild and Kruskal-Szekeres coordinates.The coordinates in Schwarzschild coordinates are denoted with lowercase letters ξ µ (Sw) = (t, r, θ , φ , u r , u θ , u φ ), with metric  µν is the metric in Schwarzschild coordinates.
Coordinates in Kruskal-Szekeres will be denoted with a capital letter, such that the coordinates in Kruskal-Szekeres are ξ µ (KS) = (T, R, Θ, Φ,U R ,U Θ ,U Φ ).The transformation from Schwarzschild coordinates to Kruskal-Szekeres coordinates (shown in figure 2) is given by where the subscript (KS) indicates the coordinates in Kruskal-Szekeres.The metric in Kruskal-Szekeres coordinates is given by where r, the radial coordinate in Schwarzschild coordinates, can be found by the inverse transformation where W 0 is the principal branch of the Lambert W function.
The numerical testing performed in this article uses moments over a size that may be considered small on astrophysical scales.This is because for numerical simulations in Kruskal-Szekeres, it is impossible to run the model with large moments, or for large amounts of time.As R is updated, the particle will eventually cross the event horizon (the point T 2 − R 2 = 1) due to numerical errors from updating the position of the particle.This will happen with any numerical differential equation solver that overestimates the true value.This is another reason it is important to transform coordinates, so a full PIC simulation can be performed in Schwarzschild coordinates, then transformed to Kruskal-Szekeres at the end, avoiding these numerical issues.Although the variation in the metric over the domain represented by the macroparticle and its moments may be small, this does not mean the derivatives of the metric, which the moments couple to, are small.

B. Computational results
To test the model, the motion of 200 particles that began normally distributed at r = 30000 in Schwarzschild coordinates with Schwarzschild radius r s = 3000, and a time step size of ∆t = 0.01 were modelled.These particles were also transformed into Kruskal-Szekeres coordinates.In both spacetimes, the particles were tracked using equation (45) and the moments taken at t = 10, and the moments were taken at t = 0 and the moments tracked using equation (46).When taking the moments, recall that f is a density, so the effect of curved spacetime adding a measure to the integrals is included in the definition of f .To convert f into a scalar field, f can be divided by |det(g)|, where the lack of square root is because the integral is over six dimensional phase space.
By modelling particles at a radius of 10r s , the accretion disc of the black hole can be studied.Once the accuracy of the moment tracking model at this distance from the black hole has been established, the accuracy of our model in the extreme environments close to the Schwarzschild radius can be examined.
The results of the tracking are shown in figure 5, with the moments used to calculate confidence ellipses in Schwarzschild spacetime.Figure 5 shows two things: firstly, whilst there is some deviation between the moment tracking and particle tracking ellipses, neither accurately reflect the underlying distribution of particles.This is because the data develops a large skew, as the faster particles orbits get thrown radially outwards.This means the particle distribution is no longer normally distributed, and as such, cannot be modelled FIG.5: The (r, u φ ) phase space portraits in Schwarzschild coordinates for the individual particles, the centre of orbit with the path η, and the ellipses used to visualise the actual moments and tracked moments of these particles.Also shown are the ellipses generated from tracking the same particles and their moments in Kruskal-Szekeres coordinates, then the moments coordinate transformed back into Schwarzschild coordinates.These ellipses show the range that 95% of particles will be within, if the particles are normally distributed.The reason that 95% of particles are not within the ellipses in 5c and 5b is because the data is no longer normally distributed.Note that if just a standard macroparticle was tracked, only the centre of orbit would be known.[Associated dataset available at http://dx.doi.org/10.5281/zenodo.8082181](Ref.[29]).
accurately with just the first and second order moments.To improve this, higher order moments will also need to be tracked.All models correctly predict the spread in the radial coordinate, such that if the r and u φ axes were equally scaled, all models would correctly predict a long, horizontally thin ellipse.The second feature is that whilst the data is initially uncoupled, the velocity and position very quickly develop a covariance (a V 16 moment), coupling the u φ and r motion together.In this particular case, this is an expected result, as particles with a faster magnitude of u φ than the one required for a circular orbit are spiralling outwards to a higher radius.
There is a noticeable displacement in the radial coordinate in the Kruskal-Szekeres coordinates tracking (figure 5b) from 30000 to 30000.08.This is because in Kruskal-Szekeres coordinates, small errors from the numerical integration will compound, and result in the ideal orbit drifting from its expected position.This compounding of errors is because Kruskal-Szekeres coordinates involves exponential functions, so small deviations can result in substantial offsets, which cannot be reduced by decreasing step size.This effect is small, and can be avoided by prescribing η beforehand, if it is known.If η is not known beforehand it is still not a major issue, as the small offset from these numerical factors is likely to be small compared to the other effects within a plasma.
To quantify the error in the model, there are six different errors that are analyzed: ε sp−sm , ε sp−kp , ε sp−km , ε kp−km , ε kp−sp , and ε kp−sm .These errors are defined in table II.The subscript sp represents that particles were tracked, then the moments taken and the end of the simulation, all in Schwarzschild coordinates.Whilst a subscript km represents moments being tracked in Kruskal-Szekeres coordinates.These subscripts are defined pictorially in figure 6.A hatted V in table II represents a coordinate transformation.As an example, ε sp−sm represents the difference between tracking the group of particles and taking their moments at the end (the black ellipse in figure 5b), compared to tracking the moments using equation ( 18) (the blue ellipse in figure 5b), all in Schwarzschild coordinates.These errors assess either the error in the moment tracking model, the error in the coordinate transformations, or the combined error of both.These six errors allow both the errors in moment tracking and coordinate transformations to be quantified and measured over time.
The error ε km−sm could also be calculated, to show the error in the coordinate transformations between two different of moments that were tracked.The origin of this error would not be discernible, as it would be impossible to distinguish between errors caused by moment tracking in either coordinate system, or the error from the coordinate transformation.
To examine the error, another simulation was performed, again around a black hole with Schwarzschild radius r s = 3000, and an ideal circular orbit at r = 30000 in Schwarzschild coordinates.For these simulations the number of particles was decreased to 20, and the time step used was increased to ∆t = 0.1, with 10 6 total iterations.This adjustment was made to allow the simulation to run for longer, to obtain better information about the long term behaviour of the model.Note the time step size is in the respective frame, so t = 10000 in Schwarzschild coordinates corresponds to T = 1200 in Kruskal-Szekeres coordinates.Figure 7 shows the three different total errors as functions of time, where the particles are normally distributed around the ideal orbit with variance 10 −20 in all dimensions.

Schwarzschild spacetime
Kruskal-Szekeres coordinates The error in all methods is small, although the error from the coordinate transformations is more substantial than the error from the moment tracking.It is postulated that this increased error is because the higher order moments affect the coordinate transformations twice, once during the coordinate transformation (during equation ( 35)), and once in the projection (during equation ( 38)).This suggests that for a moment tracking code that also incorporates a coordinate transformation, a higher order of moments will be needed.The bumps and discontinuities in the results are due to particles passing the macroparticle centre.If a particle is travelling faster than the macroparticle centre, the particle tracking moments (i.e.V sp ) will decrease, then increase once the particle passes the macroparticle centre.The moment tracking code will not see this behaviour, and will track the moments as always either decreasing or increasing, rather than the true mixture of both.These discontinuities are an intrinsic part of modelling moments.They can be avoided by sorting the particles before the modelling starts, so higher speed particles are ahead of the macroparticle centre, but this is no longer realistic.
To study the magnitude of the error, rather than just the shape, 3 main sources of error can be identified: floating point errors, numerical errors (the error arising from finite step size in numerical integration), and truncation errors (the error arising from truncating the multipole expansion at second order).These errors will dominate in different ways depending on the size of the total initial moments µ, where Floating point errors can arise from a sufficiently small time step and small number of iterations in any numerical differential equation solver, but in the case of this work they also dominate if the moments are very small i.e. µ ≈ 10 −15 .Numerical errors arise from the choice of integrator used, and in the case of forward Euler, are linear in ∆t.The truncation errors arise from only running the moment tracking code up to Error and description of the error The error between tracking particles then taking moments, compared to tracking moments, both in Schwarzschild coordinates.
The error between tracking particles and taking moments in Kruskal-Szekeres coordinates then transforming these into moments in Schwarzschild coordinates, compared to tracking particles and taking moments in Schwarzschild coordinates.
The error between tracking moments in Kruskal-Szekeres coordinates then transforming these into moments in Schwarzschild coordinates, compared to tracking moments in Schwarzschild coordinates.
The error between tracking particles then taking moments, compared to tracking moments, both in Kruskal-Szekeres coordinates.
The error between tracking particles and taking moments in Schwarzschild coordinates then transforming these into moments in Kruskal-Szekeres coordinates, compared to tracking particles and taking moments in Kruskal-Szekeres coordinates.
The error between tracking moments in Schwarzschild coordinates then transforming these into moments in Kruskal-Szekeres coordinates, compared to tracking moments in Kruskal-Szekeres coordinates.

TABLE II:
The types of error the numerical testing generates.These errors show the accuracy of both the moment tracking model and the coordinate transformations by comparing the results to a particle tracking model.Figure 6 shows these errors diagrammatically.
quadrupole order.There is an infinite expansion of moments, which is truncated to quadrupole order in this paper.Including more moments will decrease the total error.At the quadrupole level, the truncation error is quadratic.This means total error ε(µ) ≈ µ 2 .This is verified in figure 8.The error is linear in the low µ regime, where numerical errors dominate, and as µ increases, the total error increases at about ε(µ) ≈ µ 1.7 .This is a combination of the predicted quadratic increase, and the numerical error.This quadratic behaviour suggests that if the moments are half the size, the total error will be quartered.[29]).

VII. GEOMETRIC INTERPRETATION OF THE MULTIPOLE TRANSPORT EQUATIONS A. Introducing de Rham current distributions
In this section we present the transport equations and coordinate transformations of multipoles in the language of differential geometry and de Rham currents [24,25].By using this method there is a more obvious split between the V a and X ab components.This split means it is simpler to isolate each term when doing complicated calculations that mix the two terms, such as during coordinate transformations.It also means the evolution of the moments can be described in a coordinate free way.This is in contrast to working directly from equation (15), which is highly dependent on the coordinate system, in particular the choice of time slicing.The Ellis method also requires a coordinate system to define the action on a test form.
On a manifold M with tangent bundle T M the coordinate time frame hypersurface E is defined such that the Z 0 component of any vector on M is 1, The bundle of p-forms is written Λ p E such that a specific p-form (field) is denoted as α ∈ ΓΛ p E. A vector field is denoted as V ∈ ΓT E.
A distributional p-form is defined by its action on a test The definition of the wedge product, Lie derivatives, internal contractions and exterior derivatives for an arbitrary distribution Ψ are defined as ) where β ∈ ΓΛ q E and V ∈ ΓT E.
The order of a p-form distribution is defined as follows.If Ψ[λ k+1 φ ] = 0 for all φ ∈ ΓΛ q E with compact support and λ ∈ ΓΛ 0 E such that η * (λ ) = 0 (57) where η * is the pullback, then the order of Ψ is at most k.Note that a quadrupole (a distribution of order 2) also includes the dipole and the monopole terms.Lie derivatives and exterior derivatives both increase order by one.Internal contractions do not affect the order of a distribution.
Given η : R → E is a closed embedding parameterised by t, the de Rham pushforward with respect to η of a p-form α ∈ ΓΛ p R is given by the distribution The degree of the distribution η ς (α) is 6 + deg(α).Since R is a curve, the degree of α is either 0 or 1, and the degree of η ς (α) is either 6 or 7.An internal contraction decreases the degree by one and an exterior derivative increases the degree by one.Lie derivatives do not affect the degree of a distribution.
A distribution Ψ of degree 6 is a semi-multipole of order at most l if The integral curves of the Vlasov vector field are tangent to the de Rham pushforward, such that This work concerns the dynamics of a semi-quadrupole, which is a semi-multipole J of order 2 and degree 6.In coordinates, this is denoted as where L a is the Lie derivative with respect to ∂ a , and i a is the internal contraction with respect to ∂ a .Note in this representation of the semi-quadrupole, the X ab and V a term are easier to distinguish by the additional internal contraction, as opposed to the Ellis representation (equation ( 22)), in which the separation between the terms was less clear.
The de Rham current representation of a multipole can be related to the Ellis representation through the relationship and conversely In J , there is no term of the form i a L b L c η ς (X abc dt) even though this contains two Lie derivatives.Similarly, there is no X abc term in the Ellis representation.This term is included if J is a full quadrupole, as opposed to a semi-quadrupole.The advantage of the coordinate free approach used in this section is that it can be shown that the X abc term vanishes in a coordinate system adapted to η, and hence J is a semi-quadrupole in this adapted coordinate system.Since the definition of a semi-multipole is coordinate free, J is a semi-quadrupole in every coordinate system.The coordinate free semi-multipoles written in this form correspond to the semi-multipoles of ref. [25], and the electric multipoles of ref. [24].

B. The distributional transport equations
By using the language of differential geometry, a more geometric approach to understanding the origin of the transport equations can be found.To find the dynamics of multipoles defined through distributions, the equivalent of equations ( 26) and ( 27) for distributions are used, where W = W a ∂ a .The dJ = 0 condition corresponds to the conservation of charge, and the i W J condition says that the flow lines of J are the integral curves of the Vlasov field.These equations are the same as equations ( 26) and ( 27), i.e. it is equivalent to the Ellis representation.
Proof.Using equation ( 62) This is zero if and only if ∂ a J a = 0, so dJ = 0 is equivalent to (26).For the i W J term, Since i b i a is antisymmetric, we can take the symmetric part of J a W b , This is zero if and only if J a W b + J b W a = 0, so i W J = 0 is equivalent to equation (27).
Since the conditions in the geometric approach are the same as the Ellis representation, both languages give the same differential equations for the moments.

C. Coordinate transformations
Coordinate transformations can also be found using the language of distributions.The coordinate transformations for internal contractions are given by and for Lie derivatives the coordinate transformations are Note that similarly to the transformation of ∂ a , the indices in the transformed coordinate system run from (0 . . .6), whilst the original indices only ran from (1 . . .6).From these the coordinate transformations for the semi-quadrupole can be found.This is equivalent to the coordinate transformations found through the Ellis representation (equation ( 35)).
Proof.The proof of this is in the appendix A 6.
As before, the transformed quadrupole is still based on the original time slicing.In this case the projections to the new time slicing are based on the internal contraction and Lie derivatives along η.The projections are given by Pushing these through the distributions, ( 40) and ( 41) are found.These coordinate transformations are also the same as the ones found through the Ellis representation.
Proof.The proof of this is in appendix A 7.

VIII. CONSIDERATIONS FOR IMPLEMENTING MACROPARTICLES WITH MOMENTS IN A FULL PARTICLE-IN-CELL CODE
With the work presented, it is possible to calculate all the dynamics for a single macroparticle in an arbitrary spacetime and external electromagnetic field.To convert this into a full PIC code, a process that uses the moments to reconstruct the distribution function and deposit the charge over multiple grid points must be developed.This will allow inter-macroparticle effects to be modelled.This is work currently being undertaken.Once this deposition process has been developed, a standard finite-difference time-domain (FDTD) algorithm can be used to update the fields, so a full PIC code can be developed.
One approach to use moments to deposit the charge and current onto the grid is to use a model function [11,30].A model function is a function that is assumed to be of a similar shape to the actual distribution of particles, then shaped to have to have the required moments.This method works well in cases where the distribution of particles is similar to a top-hat function i.e. the distribution of particles in a plasma.In particle accelerators, the distribution of particles around the ideal orbit is not a top-hat, but rather a Gaussian.The model function depositing method does not work well for modelling these.An alternative approach in this case is to model the particle distribution function as a series of Hermite polynomials [31].This can be related to the moments, to approximate the distribution function using the moments and Hermite polynomials.It is possible to approximate the conventional plasma moments (the pressure tensor, the energy flux density etc.) used in magnetohydrodynamics by using the higher order moments to approximate the underlying distribution.This is through the same process as depositing charge and current onto the grid, and requires tracking moments of a higher order than the plasma moment e.g.finding the pressure tensor would need the hexadecapole moments to be tracked.The reconstructed model function is normalised to have the same monopole moment as the macroparticle, this ensures no extra charge is added into the system through deposition.As this method does not use shape functions that conserve charge exactly [32,33] it may be necessary to perform divergence cleaning to stop accumulation of charge through the FDTD Maxwell solver (See ref. [34] for more detail on this procedure).It should be noted that the exactly charge conserving shape functions can only be used in Cartesian coordinates, and so cannot be applied to the curved spacetimes presented in this work.
As mentioned in the introduction, this work only considered particles interacting with external fields, and as such, it was not meaningful to talk about cell sizes for the system.In a full PIC code, further work will need to be done to find the appropriate cell size for a macroparticle with moments.As shown in this work, the dominating source of error in the moment tracking is when the field cannot be accurately modelled by a small number of derivatives at the macroparticle centre.This means that one limit for the cell size will be based on the density of grid points required to accurately find the derivatives of the fields.Because the macroparticle has an internal structure, it needs to span over several cells to allow the extra structure from tracking the moments to be reflected in the deposited charge and current.This is in contrast to existing methods, which deposit over several cells to reduce numerical instability [34,35].It may be that in cases where the fields can be approximated using a low resolution grid, that it is possible for a single macroparticle to be larger than the Debye length of the plasma.By tracking the moments, information about the structure within the macroparticle is known, and this affects how the charge and current will be deposited onto the grid points.Because of this, it may be that the effects of instabilities generated from using a finite grid (aliasing instabilities) [34,[36][37][38] are not as significant.

IX. CONCLUSION AND DISCUSSION
This paper found the dynamics of the moments of a macroparticle obeying the Vlasov equation (equation ( 18)), and the coordinate transformations of the quadrupole moments (equations ( 40) and ( 41)).By using the Ellis representation or the de Rham current representation of the moments, coordinate transformations can be found between frames that mix the space and time coordinates.By representing a group of particles as a macroparticle and its moments, this can be used in PIC codes.These results were validated numerically for the case of particles orbiting a black hole.This was done by transporting particles and moments in both Schwarzschild and Kruskal-Szekeres coordinates and comparing results.The numerical results show that a large number of particles can be successfully modelled by a single macroparticle with moments, although for a full PIC code it is likely moments will need to be calculated to an order larger than the quadrupole.In addition, more macroparticles will need to be modelled, and a method to use the structure of the moments to deposit the charge and current onto the grid will need to be developed.By doing these, a full PIC code can be developed, after which it is possible to assess how the full process of moment tracking benefits PIC codes.
The dynamics of the moments depend on the Vlasov field and its derivatives.The focus of this paper was on how the number of moments taken affects the accuracy of the moment tracking.For a full PIC code, the derivatives of the Vlasov field will also be important.If the electromagnetic or gravitational fields quickly vary in space (i.e. the higher order derivatives of the fields are large), then this will mean more moments need to be tracked.The moment tracking model is likely to work well in situations where both the distribution of particles represented by a macroparticle can be described by only a small number of moments, and the variation in electromagnetic and gravitational fields across the volume the macroparticle represents is small, such that the Vlasov field across the extent of the macroparticle can be modelled by just the Vlasov field and a small number of its derivatives at the macroparticle centre.Since numerically calculating the derivatives of the electromagnetic field requires a high density grid, the moment tracking method may also work in cases where a high resolution grid is already needed, such as laser-solid interactions [39,40].In such cases, the moment tracking method will be able to model much larger macroparticles, with less particles per cell, even if more cells are needed to compensate.
The modelling in this work was done at 10 Schwarzschild radii away from the black hole.This was to model a stable accretion disc.In the development of the theory, there was no assumption that the macroparticle is far away from the black hole, so there is nothing in the theory that stops a macroparticle being close to the singularity.In practice, close to the Schwarzschild radius, the error from only using a finite number of moments will become significant.There are two reasons for this, firstly, the particles being represented by the macroparticle and its moments will begin to undergo spaghettification, and this will cause the higher order moments to become large.Secondly, close to the singularity, the Schwarzschild Christoffel symbols become very large, and their derivatives will become even larger.This means the rate of change of the moments will be dominated by the higher order moments which are being neglected in the code.Because of this, the accuracy of the model is yet to be determined in the extreme environments close to the Schwarzschild radius.
Another potential application of the moment tracking method is to model inter-bunch forces within particle accelerators.It is possible to calculate the Liénard-Wiechert fields directly from the moments of a moving quadrupole [6].By using this method the electromagnetic field, and its derivatives, can be calculated without the need to deposit the charge onto a grid.This is particularly useful for modelling coherent synchrotron radiation in particle accelerators, where macroparticles are close together compared to the radius of the beam pipe, such that the effects of boundary conditions on the electromagnetic fields can be ignored.
It may be possible to add internal dynamics into the moment equations.Whilst this has been previously done in refs.[41] and [2] for specific applications of muon cooling and space-charge respectively, it may be possible to come up with general internal dynamics equations by modifying the transport equations (equation ( 18)).This would allow the moment tracking method to model intra-bunch effects within particle accelerators.
The applications of the coordinate transformations in astrophysical scenarios are wide.A particularly useful case is the transformation from the global time to the backwards light cone frame.This is particularly useful in cases where black holes are modelled in the fiducial observer (FIDO) frame, where the global time coordinate is significantly different to the backwards light cone frame global time.By doing this the difficulty of calculating the backwards light cone through ray tracing only needs to be done once, rather than each time step.This transformation would allow the moments observed by the observer a finite distance away from the black hole to be calculated.
There are also applications of the coordinate transforma-tion formulae in circular particle accelerators.Results from accelerators are often presented in Frenet-Serret coordinates, where the parameter is the position along the beamline, rather than time.By finding either the Vlasov equation in Frenet-Serret coordinates, or finding the coordinate transformation between Cartesian and Frenet-Serret coordinate systems for a given beamline, the moment tracking can be applied to circular accelerators.Additionally by using the spacetime coordinate transformations presented in this article, the coordinate transformation into the frame of an accelerating bunch can be found (this is a similar transformation to the one presented in figure 3).
Lastly the use of the Vlasov equation to model dynamics may be extended to modelling stress-energy-momentum quadrupoles as a source for linearised gravity.In both refs.[25] and [23] it was shown that the dynamics of stressenergy-momentum quadrupoles contain a number of free components, known as constitutive relations.In the case of a plasma these constitutive relations may be determined by the Vlasov equation.In this case the dynamics will be governed by the divergenceless of the stress-energy-momentum tensor combined with the Vlasov equation.
This section shows that by acting an Ellis distribution J a on specific test forms, the components of a multipole can be extracted.If it is possible to do this, the components of a multipole are unique.The components are isolated using specific tests forms: where ψ : R → R is a test function such that ψ 1 (0) = 1, it is flat about zero and ψ 1 (t)dt = 1, t 0 is the point at which the moments are evaluated, and E, the coordinate time frame hypersurface, is defined by equation (54) in section VII.
Proof.Only the V bc term and X bc + V b ηc term will be shown as the other terms follow trivially.
Consider J a acting on the V bc equation of (A1).The only non-zero derivatives are the ξ b and ξ c terms.There are three possibilities, each ξ b or ξ c can either not be differentiated, differentiated once, or differentiated twice.If it is not differentiated, then the evaluation at η gives (ξ b − η b )| η = η b − η b = 0.If it is differentiated exactly once, then ∂ a (ξ b − η b ) = δ b a , a Kronecker delta.If this is differentiated twice, then the derivative of a Kronecker delta will vanish.Thus the only non zero term when acting J a on the V bc equation of (A1) is the term where the number of derivatives matches the number of ξ b terms.In this case this happens when there are exactly two partial derivatives.This gives Noting ξ a | η = η a , η0 = 1, V bc = V bc (t) and introducing the substitution t = t 0 + εt ′ gives lim .

(A17)
Taylor expanding φ 0 , φ a and W a around η, noting there are only spatial derivatives as ξ 0 | η − η 0 = t − t = 0,  Using the definitions of q, V a , and V ab from equation (A12), Next, expand the partial derivatives, and note W a | η = ηa , For the V a term, Summing all these terms together gives the transformed quadrupole, This gives the coordinate transformations for the quadrupole components, where the components of the quadrupole in the new coordinate system are no longer unique.

Proof of the Ellis multipole projection
This section performs the projection ∂ 0 = ∂ η − ηa ∂ a into the quadrupole with non unique components defined by equation (35).By performing this projection U a ,U ab ,Y a and Y ab can be written in a form where the components are unique, giving the full coordinate transformation for the quadrupole.
Inserting (38) term by term into the non-unique semi-quadrupole

FIG. 2 :
FIG.2:A spacetime diagram in Kruskal-Szekeres coordinates of a particle travelling at a constant r in Schwarzschild coordinates (the orange hyperbola).The diagonal blue lines are time slicings in Schwarzschild coordinates (constant t), and the horizontal purple lines are time slicings in Kruskal-Szekeres coordinates (constant T ).The dashed black line is the event horizon.The moments in a given coordinate system have no components in the direction of the time slicing in that coordinate system.
(a) t = 0, all ellipses overlap.(b) t = 10 showing particles and moments tracked in Schwarzschild coordinates.(c) t = 10 showing the same particles transformed to Kruskal-Szekeres coordinates, then the moments and particles tracked, then transformed back to Schwarzschild coordinates.

FIG. 7 :
FIG.6:The model used to test the moment tracking and coordinate transformation theories.The error in the moment tracking model is the difference between transporting the moments and transporting the particles then taking moments.The error in coordinate transforming is the difference between transporting particles then taking moments in each frame.The combined error is the difference between transporting the moments in one frame, compared to tracking particle in the other.These errors are defined algebraically in table II.
FIG.8:The total error as a function of µ, the initial total moment, for the three different kinds of errors considered.The gradient of the lines is approximately linear in low µ, and approximately 1.7 for larger µ. [Associated dataset available at http://dx.doi.org/10.5281/zenodo.8082181](Ref.[29]).

J a φ a d 6 ξ = 1 2 R
ηa U bc ∂ b ∂ c φ a dt + R Y ab ∂ b φ a dt + R ηa U b ∂ b φ a dt + R Y a φ a dt + R ηa qφ a dt.(A50)