Concept of the generalized reduced-order particle-in-cell scheme and verification in an axial-azimuthal Hall thruster configuration

The reduced-order particle-in-cell (PIC) scheme is a novel modeling approach that enables computationally efficient electrostatic kinetic simulations of plasma. In our previous publications, we demonstrated that a proof-of-concept implementation of this novel PIC scheme resolves the multi-dimensional plasma processes and their interactions in a Hall thruster in a manner close to traditional electrostatic PIC codes. In this work, we extend our efforts on this topic and present a mathematically mature formulation for the dimensionality reduction of Poisson’s equation in the Vlasov–Poisson system, which enables the generalized reduced-order ‘quasi-multi-dimensional’ PIC scheme. The applicability of the dimensionality-reduction approach to solve general 2D Poisson problems is numerically verified. Next, we present several reduced-order quasi-2D (Q2D) simulations of a well-defined axial-azimuthal simulation case from the literature using approximation orders of the 2D problem whose computational costs are 2%–15% of a full-2D simulation. It is shown that these reduced-order simulations allow us to recover the same characteristics, behaviors and effects reported in the literature regarding the azimuthal instabilities in Hall thrusters. Moreover, in terms of the time-averaged plasma properties, it was found that, when increasing the approximation order, the error associated with the Q2D simulations’ predictions decreases from 15% to 4% for the electric field and from 20% to 2% for the ion number density. We have additionally discussed a series of sensitivity analysis results, including the influence of the initial number of macroparticles per cell on the predictions of the Q2D simulations. According to the detailed results and analyses presented, we conclude that the generalized reduced-order PIC scheme serves as a rigorous foundation for eventual cost-effective and comprehensive three-dimensional kinetic studies of the physics in Hall thrusters and similar electrostatic plasma technologies.


Introduction
Hall thrusters are plasma-based in-space propulsion devices that can be used in a variety of mission scenarios, ranging from on-orbit servicing to near-Earth transportation and interplanetary exploration. They rely on electromagnetic fields to ionize the propellant and to accelerate the created plasma to generate thrust. Although Hall thrusters are relatively simple from an engineering standpoint and their basic operating principle is quite straightforward to understand, their underlying physics of operation is rather intricate. The complexity arises, in part, from the fact that, in the presence of a perpendicular configuration of electric and magnetic field in these thrusters, the plasma is subject to strong gradients and anisotropies, which are sources of a variety of instabilities and turbulence across a vast range of spatial and temporal scales. From an operational perspective, these instabilities significantly affect the thruster's performance and lifetime, which ties the understanding of physics to the development of efficient devices and their reliable in-space application.
Accordingly, Hall thrusters have gathered notable scientific attention in the past decades. A considerable portion of research on Hall thrusters has been focused on deriving insights from the investigation of the plasma processes and interactions using high-fidelity kinetic simulations to realize predictive models of Hall thrusters in a bid to enable their cost-effective numerically aided design and testing [1,2]. However, despite numerous efforts to this end [3][4][5][6][7], reliable prediction of the plasma behavior in Hall thrusters is still an elusive target. This is partly because the dynamics of the charged particles and, in particular, the electrons' cross-field mobility, in the presence of the plasma instabilities and turbulence is not yet fully understood.
In this regard, it is becoming increasingly evident that a comprehensive understanding of the three-dimensional kinetic nature of the plasma dynamics in a Hall thruster is essential to fully comprehend the interplay among various plasma phenomena in the thruster and to establish links between these phenomena and their effect on the global discharge behavior [8]. Achieving this goal requires 3D kinetic simulations since important aspects of the underlying physics can be lost by resorting to lower-dimensional 2D kinetic codes. However, 3D kinetic simulations are currently computationally unfeasible for real-world thruster geometries. Thus, in the absence of a comprehensive knowledge regarding the multi-scale, multi-dimensional plasma phenomena and their complex interactions in Hall thrusters, cost-effective predictive numerical simulations are yet to be achieved.
Noting the above-mentioned necessity of 3D kinetic codes and to address the issue of enormous computational cost associated with conventional multi-dimensional kinetic simulations that renders them impractical for real-size Hall devices [9], we introduced the 'pseudo-2D' electrostatic particle-incell (PIC) scheme in a prior publication [10]. Underlying this PIC scheme was a preliminary formulation to split the 2D Poisson's equation into a set of decoupled 1D ordinary differential equations (ODEs). We demonstrated in [10,11] that this splitting of Poisson's equation allows us to remarkably reduce the required number of computational cells and the total number of macroparticles, hence, gaining significant computational advantage over traditional PIC codes. We performed proof-of-concept studies on the pseudo-2D kinetic simulation along various coordinates of a Hall thruster, i.e. axial-azimuthal, azimuthal-radial, and axial-radial, comparing the results against the observations in full-2D reference cases [10,11]. In this respect, despite the rather simple dimensionality-reduction formulation that were used in the 'pseudo-2D' simulations, the essential capability of the approach to resolve the main physics and interactions among plasma phenomena in different 2D configurations was demonstrated [10,11].
According to the remarkable potentials of the pseudo-2D PIC, we present in this article the generalized reduced-order or 'quasi-multi-dimensional' PIC scheme. The rigorous and flexible implementation of our PIC scheme, on which we will elaborate in the next section, is achieved through the development of a mathematically mature, generalizable formulation for the decomposition of Poisson's equation in the Vlasov-Poisson system. We believe that the current implementation of the reduced-order PIC allows us to ultimately realize computationally affordable high-fidelity 3D kinetic simulations of Hall thrusters.
Nevertheless, before proceeding to 3D, the novelty of our approach necessitates performing meticulous verifications. In this regard, our step-by-step, ground-up plan to verify the reduced-order PIC scheme involves carrying out quasi-2D (Q2D) simulations in various 2D configurations of a Hall thruster and comparing our results against those from full-2D codes. Accordingly, in this paper, following the introduction of the concept of the generalized reduced-order PIC, we report the results of its in-depth verification in the axial-azimuthal coordinates of a Hall thruster using the so-called 'high-order' approximations of the 2D problem. This is to demonstrate that, in this approximation limit, the results from our reduced-order PIC simulations practically converge to those from the full-2D ones, capturing the same physics but at a fraction of the computational cost.

The concept of the reduced-order PIC scheme
The reduced-order PIC scheme is first and foremost enabled through a dimensionality-reduction technique to decompose the multi-dimensional Poisson's equation into several onedimensional components along the involved simulation directions. Accordingly, we first present in section 2.1 the mature formulation of the dimensionality-reduction technique. Second, in section 2.2, we discuss the main algorithmic modifications that are needed in a conventional PIC code to accommodate the dimensionally reduced Poisson's equation and, hence, to achieve the reduced-order PIC scheme. We also present the computational advantage of the reduced-order PIC compared to the conventional schemes.
Without losing generality, and in line with the objective of this work to verify the reduced-order PIC scheme in a 2D axial-azimuthal configuration, we focus from this point forward on a two-dimensional problem and, thus, the 2D Poisson's equation.

Formulation of the dimensionality-reduction approach for Poisson's equation
To introduce the basic concepts behind the dimensionalityreduction approach, we first refer to figure 1 in which a 2D domain with the extents L x and L y along the x and y directions, respectively, is represented by two perpendicular computation grids, one along x and one along y. The computational cells along the x-direction have the width ∆x, determined according to the physics of the problem, and are extended over the entire L y extent of the domain. Similarly, the computational cells along y have the width ∆y and a length equal to L x . We term this representation of the domain the 'single-region' decomposition, since we have a single set of perpendicular grids representing the entire domain.
We define two potential functions ϕ x and ϕ y , which are, by definition, only a function of the x-and y-coordinate, respectively, within a region. We approximate the 2D potential function ϕ (x, y) as the superimposition of these two potential functions, i.e., Before moving on with the derivations, we note that the above approximation of the 2D potential function associated with the single-region domain decomposition has a limited applicability in terms of the solutions it can reproduce. However, the approximation given by equation (1) enables the lowest order quasi-multi-dimensional simulations (i.e. the single region) that offer the maximum computational advantage over traditional PIC codes. Of course, a single-region simulation may not provide reasonable predictions in cases of highly multi-dimensional geometries and/or plasma processes.
Nonetheless, the generalization of the single-region domain decomposition to a 'multi-region' one is demonstrated numerically in the subsequent sections of this article to allow recovering the full-2D solutions, both for general Poisson problems and in PIC simulations. It is also important to note that, in the limit where the lengths of regions become as small as 2D cells' size, equation (1) will no longer be an approximation and becomes exact. We present a detailed mathematical proof of convergence for the dimensionality-reduction approach in section A.1 of the appendix.
Following the above clarification, substituting equation (1) in the 2D Poisson's equation and integrating once along the horizontal (x) and once along the vertical (y) computation grid yields a coupled system of 1D ODEs for ϕ x and ϕ y , respectively, In the above equations, ϵ 0 is the permittivity of free space, and x g and y g are the locations of the y and x computation grid, respectively, with respect to the origin. ρ (x, y) is the 2D charge density distribution.
In order for the dimensionality-reduction approach and, hence, the reduced-order PIC scheme to be readily generalizable, we extend the single-region decomposition formulation described above to the multi-region decomposition (figure 2), where the 2D domain is now decomposed into a grid of N × M regions, denoted by Ω n,m . The regions' boundary is defined such that it is equidistant from the two computational nodes adjacent to it. Within each region, there is a set of perpendicular grids along which the 1D potential functions ϕ x n and ϕ y m are resolved. Moreover, the 2D potential function is now approximated as ϕ (x, y) = ϕ x n + ϕ y m , with x, yϵΩ n,m .
Accordingly, equations (3) and (4) are written as in equations (6) and (7) for the multi-region decomposition, yielding the general formulation of the dimensionalityreduction approach,  ρ (x, y) dx.
In equations (6) and (7), L x,m and L y,n are the horizontal and vertical extents of each region Ω n,m , x i is the x-coordinate of each computational node i along a horizontal grid and y j is similarly the y-coordinate of each node j along a vertical grid.
Having solved the coupled system of equations given by equations (6) and (7) for ϕ x n and ϕ y m with the appropriate boundary conditions, to be described shortly, the electric field is determined within each region using ⃗ E = − ⃗ ∇ϕ (x, y) = − ⃗ ∇ (ϕ x n (x) + ϕ y m (y)). In equations (6) and (7), it is important to note the appearance of the so-called 'cross-derivative' terms, ∂ϕ x n, ∂y and ∂ϕ y m ∂x , that take into account the variations in ϕ x n and ϕ y m when transiting from one region to another along the y-and x-direction, respectively.
Moreover, recalling that, inside each region Ω n,m , ϕ x n is, by definition, constant along y, and ϕ y m is likewise constant along x, the integral terms within the second parentheses on the lefthand side of equations (6) and (7) are zero for the computational cells that reside entirely inside each region. In other words, the integrand terms  (7) are non-zero only when one side of a computational cell falls on a region's boundary.
The boundary conditions accompanying equations (6) and (7) are of the following general form The formulation presented above for the dimensionalityreduction approach provides a generalizable 'quasi-multidimensional' approximation of the Poisson's problem, which is included in the PIC scheme to yield the reduced-order PIC.

Numerical implementation of the reduced-order PIC scheme and its computational benefit
The reduced-order PIC scheme is currently implemented as an electrostatic, explicit code, which is named IPPL (Imperial Plasma Propulsion Laboratory) Q2D code. The IPPL-Q2D code is written entirely in Julia language [12]. It, overall, follows the standard algorithm of the PIC method [13,14]. The main difference between the reduced-order PIC code and a conventional one pertains to the module solving the electric potential. This module, which we term the 'Reduced-Dimension' Poisson Solver, or RDPS, is explained in the following. The other only difference is in the module that scatters the particle-based data (such as the electric charge) onto the grids and gathers the grid-based data (such as the electric field) onto the particles' position. The particles' information in the reduced-order PIC is scattered onto the computation grids corresponding to the dimensionality-reduction approach, i.e. in the most general case, N horizontal grids along y and M vertical grids along x (horizontal and vertical green lines in figure 2). Similarly, the field information used to update the particles' velocity and position are gathered from these computation grids.
The RDPS numerically solves the system of equation given by equations (6) and (7) using Julia's built-in direct matrix solve algorithm based on the LU decomposition method. However, for the solver's implementation in the PIC scheme, instead of directly solving the system of coupled equations (6) and (7), the terms corresponding to the gradients of the potential function ϕ y m in equation (6) and, similarly, the terms representing the gradients of ϕ x n in equation (7) are moved to the right-hand side and are, thus, calculated from the previous timestep of the simulation. In this way, the system of equation is kind of decoupled, which is computationally beneficial as the size of the matrix to be solved at each time step is reduced by half. Moreover, we observed an increased robustness and stability of the solver when following this approach without any noticeable issues concerning the numerical artifacts and/or instabilities.
It is also important to mention that, before calculating the right-hand-side derivative terms specified in the preceding paragraph and solely for the purpose of these calculations, a moving average filter with a pre-fixed window size is applied to the solutions of potential functions ϕ x n and ϕ y m . The analysis of the sensitivity of the reduced-order simulation results to the size of the filter's window is presented in section A.2 of the appendix.
Following the above description of the numerical implementation of the reduced-order PIC scheme, we demonstrate quantitatively the computational advantage that our novel PIC scheme has over traditional simulations in a 2D case. In this regard, table 1 provides a comparison between the full-2D and the Q2D PIC schemes in terms of two simulation parameters, the number of cells (N cells ) and the total number of macroparticles at the beginning of the simulation (N p,init ).
For the comparison in table 1, we have considered a domain with N i = 500 cells along the x-direction and N j = 256 cells along the y-direction. These are the same number of computational cells used for the full-2D axial-azimuthal benchmark case in [15,16], with which we will later compare the results of our reduced-order Q2D simulations. Looking at the last column of table 1, the single-region Q2D PIC results in 256 times reduction in the computational cost. Moreover, decomposing the domain into multiple regions, e.g. N = 20, M = 40, which as will be seen in section 4 practically provides the full-2D results in all respects, still translates into a computational gain by a factor of about 6. Table 1. Comparison between a full-2D and quasi-2D PIC simulation in terms of the grid size (N cells ) and the total number of macroparticles at the simulation start (N p,init ); N PPC is the number of macroparticles per cell.

Numerical verification of the generalizability of the dimensionality-reduction formulation and the resulting RDPS
The discussions in this section are aimed at verifying the generalizability of the dimensionality-reduction approach itself. We complement the numerical verifications presented here by a mathematical proof in section A.1 of the appendix to demonstrate that the formulation developed for the dimensionality reduction, discussed in section 2.1, mathematically converges to the 2D Poisson's equation in the limit where the number of regions is the same as the number of 2D computational cells, i.e. N = N j and M = N i . Concerning the numerical standalone verifications in the following, we intend to demonstrate that low-number-ofregion approximations of a 2D Poisson problem are sufficiently representative of the problem's full-2D solution and that, by increasing the number of regions or, equivalently, the approximation order, the approximate solutions from the RDPS converge to the 2D solution.
To this end, three canonical two-dimensional Poisson problems are chosen. The definitions of these problems together with their corresponding boundary conditions are presented in table 2. In all cases, the domain is a Cartesian (x − y) plane of unit length along the x and y directions. For simplicity, the permittivity of free space (ϵ 0 ) is assumed to be one for all problems. The number of computational cells for all cases are the same and equal to 200, i.e. N i = N j = 200.
Referring to table 2, it is worth pointing out the reasons for the selection of these three specific Poisson problems. Regarding Case 1, the problem has a rather complex analytical solution ϕ (x, y) = exp (π x) sin (π y) + 1 2 (xy) 2 . As a result, the ability of the RDPS to provide accurate numerical approximations of this solution proves the versatility of the approach. Moreover, the importance of Case 2 lies in its periodic boundary condition along the y direction. In this respect, as this boundary condition is a typical one applied along the azimuth in Hall thruster simulations that involve this coordinate, it was deemed necessary to demonstrate the applicability of the reduced-dimension solver to the problems with periodic boundary conditions as well. Finally, Case 3 represents a problem with an oscillatory charge distribution of relatively small wavelength. This is a common situation occurring in a Hall   figure 3 illustrates the analytical solution for Case 1, namely, ϕ (x, y) = exp (π x) sin (π y) + 1 2 (xy) 2 , whereas, in figure 4, plot (d) depicts the full-2D solution of Case 2 obtained using Julia's direct Poisson solver function.
Looking at the plots in figure 3, it is observed that, despite the complex nature of the analytical solution, the RDPS is able to capture all of the solution's main features (e.g. the symmetry in the y-direction and the convex curvature along x) even at a Furthermore, in both Cases 1 and 2, as the number of regions N and M are increased, the approximate solution is seen to evolve toward the 'ground-truth' solution such that, at M = N = 50, the reconstructed 2D solution from the RDPS has essentially converged to the 'true' one. The same converging behavior of the solution is observed for the Poisson problem in Case 3 (figure 5). Indeed, from the right-hand-side plot of figure 5, which shows the distribution of the error, as defined by equation (11), between the approximate solution at M = N = 100 (figure 5 (left)) and the full-2D solution, the maximum error is found to be ±1.2%, which is in the same order of the error corresponding to the full-2D solver itself,  The convergence characteristic of the RDPS and the decrease in error between the approximate and the groundtruth solutions when increasing the number of regions can be better understood by looking at the plots of figure 6, which present the evolution of the root-mean-square (RMS) error (equation (12)) for different Poisson problems, It is observed that, in Case 1 and 2, the RMS error falls rapidly when increasing the number of regions, approaching error values in the same order of a full-2D Poisson solver. For Case 3 (figure 6(c)), even though the error falls below 10 −5 % from about M = N = 62, the decrease in the error is rather linear in contrast to the other two cases. However, this is due to the extreme nature of charge density gradient in Case 3, which implies that higher-number-of-region approximations are necessary to provide solutions that are closely representative of the full-2D. Hence, it is possible to conclude that the error convergence rate of the RDPS and, in turn, the reducedorder PIC vs the number of regions is problem dependent and is primarily determined by the significance of the gradients in ϕ x and ϕ y from one region to the adjacent ones, i.e. the importance of the 'cross-derivative' terms in equations (6) and (7).
As the final point, the reason why the number of regions in the test cases of this section, which all corresponded to a 2D 200 × 200 grid, were not increased beyond 100 is that our numerical implementation of the dimensionality-reduction formulation requires that at least two nodes fall within each region.

Verification of the reduced-order PIC scheme in Q2D axial-azimuthal Hall thruster simulations
In this section, we present the verification of the reduced-order PIC scheme in an axial-azimuthal Hall thruster configuration using 'high-order' approximations of the 2D problem. For this purpose, we have adopted the well-defined 2D axial-azimuthal simulation case of [15]. The simulation setup and conditions are explained in section 4.1.
The 'high-order' Q2D simulations are performed using various multi-region decompositions of the 2D domain comprising the following cases:

Overview of the simulation setup and conditions
The setup of the simulations is largely similar to that described in [15,16]. The simulation domain is a 2D Cartesian (x − y) plane, with x along the axial direction and y representing the azimuthal coordinate. The adopted computational and physical parameters, including the domain's extents, cell size, time step, total simulation time, and the initial conditions, are also the same as those for the full-2D reference case [15,16], which are summarized in table 3 as well.
Nonetheless, there are some aspects of the setup that are specific to the Q2D simulations. We specify these in the following paragraphs.
Concerning the initial total number of macroparticles and, hence, the initial macroparticles per cell count in the Q2D simulations, to ensure that a minimum of about 100 macroparticles per cell along the azimuthal direction is maintained throughout the simulation in all cases, we chose the initial total number of macroparticles to be 250 000 for case i, 500 000 for case ii, 1000 000 for case iii, and 2000 000 for case iv.
For the injection of the electrons into the domain to sustain the discharge, we applied the quasi-neutrality (QN) condition [16] on the cathode plane at x c = 2.5 cm for all simulations presented in this article. Therefore, the results from our Q2D simulations are compared against those results reported in [16] that correspond to the QN electrons' injection scheme from the cathode.
The plasma potential in the Q2D simulations is solved using the RDPS explained in section 2.2. The associated boundary conditions, whose general form were given by equation (8) in section 2.1, are specified for the PIC simulations as follows: along the azimuthal computation grids, a periodic boundary condition is considered, which is represented as a zero-volt reference potential applied to both ends, i.e. ϕ y m (y b ) = 0, where y b denotes the domain's boundary locations along the azimuth. This implementation of the periodic boundary condition is in line with the 2D reference simulation case [15,16].
Along the axial computation grids, Dirichlet boundary conditions of 200 V at the anode side (x a = 0) and 0 V at the cathode side of the domain are implemented using the relations below Now, similar to the reference 2D case [15,16], the Q2D simulations feature a Gaussian magnetic field profile and a cosine-shaped, time-invariant imposed ionization source. The axial profiles of the magnetic field intensity and the ionization source are shown in figure 7. The peak magnetic field intensity is 100 G, whereas the maximum current density corresponding to the peak of the ionization source is 400 A m −2 .

Results and discussion
The results from our Q2D axial-azimuthal simulations are presented in this section in terms of the evolution of the normalized electron and ion currents, the time-averaged axial profiles of the plasma properties, the 2D reconstruction of the plasma properties' distribution, and the dispersion characteristics of the azimuthal instabilities. Figure 8 shows the temporal variation of the normalized electron and ion currents from the Q2D simulations. The ion current (I i , figure 8(b)) corresponds to the flux of the ions leaving the cathode boundary of the domain. Moreover, the difference between the electron and ion flux reaching the anode boundary gives the 'discharge' current, whose difference with respect to I i , in turn, yields the electron current (I e , figure 8(a)). In the plots of figure 8, the y-axis is normalized with respect to the total ion current, i.e. J · A, where J is 400 A m − 2 and A is the domain's cross-sectional area. In addition, the dashed horizontal lines indicate the mean value of the currents, calculated from t = 5 µs.
Looking at plots in figure 8, it is observed that, after an initial transient, all simulations have arrived at the quasi-steady state from about 5 µs. Furthermore, increasing the number of regions from N = 5, M = 10 (case i) to N = 40, M = 80 (case iv), the mean values of the normalized electron and ion currents converge to about 3 and 0.8, respectively. These mean current values are very similar to those from the reference 2D simulation [16]. Finally, it is interesting to note that, despite a rather different currents' trace between case iii (N = 20, M = 40) and case iv, their mean currents are almost the same, implying that these two simulations are expected to provide quite similar predictions in terms of the time-averaged macroscopic plasma properties. We will see in the following that this is indeed the case.
The axial distributions of the plasma properties from the Q2D simulations with different numbers of regions N and M are shown in figure 9. The profiles are averaged over the last 10 µs of the simulation time to be consistent with the 2D reference case [16].
In the plots of figure 9, the axial and azimuthal components of the temperature are calculated using the relations T y (x, y) = m en In the above expressions, v x and v y are the axial and azimuthal velocities, V d, x and V d, y are the axial and azimuthal drift velocities, f (v) is the particles' distribution function, n is number density, e is the elementary charge, and m is the particle's mass.
The plots in figure 9 clearly indicate that increasing the number of regions lead to the convergence of the plasma profiles to those from the 2D simulation. To provide quantitative accuracy metrics for the 'highorder' Q2D simulations, we refer to figure 10(a), which shows the variation of the errors in the axial electric field and ion number density with respect to the number of regions used along the axial direction (M). The error in the electric field is calculated based on the difference in the field's peak intensity between each Q2D profile and the 2D one. The error in the ion number density is the axially averaged node-by-node difference between the Q2D and full-2D profiles.
We observe from plot (a) in figure 10 that the error in terms of both plasma properties decreases from case i to case iv. The error values for case iv are almost within the margin of error observed between various 2D simulations performed by different participating groups in the Hall thruster axial-azimuthal benchmarking activity [15]. This is while it is evident from figure 10(b) that the Q2D simulation of case iv takes only about 15% of the time required for a full-2D simulation. Furthermore, it is interesting to note that the Q2D simulation of case i (M = 10) only requires about 2% of the 2D simulation's time to complete whereas its error figures of about 20% in ion number density and about 15% in peak electric field are fairly reasonable.
Before proceeding further with the discussions and looking into the 2D distributions of the plasma properties and the characteristics of the resolved azimuthal fluctuations, it is noted that one of the main instabilities exciting along the azimuthal direction in plasmas with E × B configuration, such as that in a Hall thruster, is the electron cyclotron drift instability (ECDI) [19,20]. This unstable mode features a rapid nonlinear growth, followed by a transition to ion acoustic instability which has a smaller growth rate [21,22]. The details of the  initial non-linear evolution phase of the ECDI are not detectable in the present simulations and, hence, we only focus in the following on the pseudo-saturated state of the instabilities, represented by the modified ion acoustic wave modes [23]. From the literature, the wavevector of the ion acoustic instability in Hall thrusters is known to be mostly along the azimuthal direction with a finite component along the axial direction as well [24].
Having the above remarks in mind, we refer to figure 11, which provides the reconstructed 2D snapshots of the distributions of the plasma properties, the ion number density and the azimuthal electric field, in the x − y plane from the Q2D simulations at time t = 30 µs.
Looking at the reconstructed 2D maps, it is observed that, although in case i (N = 5, M = 10), the simulation predicts fluctuations that have only an azimuthal wavevector  figure 11. This is the evidence that, in these number-of-region limits, the reduced-order PIC is able to resolve the axial wavevector component (k x ) of the azimuthal instabilities. This observation, together with the error and computational cost figures presented in figure 10, points to the remarkable capability of the reduced-order scheme in the high-order approximation limit to capture the underlying physics in the very same manner as a full-2D simulation but at a notably reduced computational cost.
To explain why in cases iii and iv the Q2D simulations can resolve the axial wavenumber of the azimuthal instabilities, it is underlined that, should the regions' axial extent be smaller than the waves' axial wavelength, i.e. l x,m ≪ λ x , the simulations are indeed expected to be able to capture the axial component of the wavevector. This criterion is met in both cases iii and iv. For instance, focusing on case iii, as it can be noticed from the x − y fluctuation map of the azimuthal electric field in figure 11, the azimuthal wavelength (λ y ) of the waves is about 1 mm. Now, we recall that, according to the literature [25], the instabilities' wavelength along the axial direction is much larger than that along the azimuthal direction (λ x ≫ λ y ). Therefore, considering that the axial length of the vertical regions (l x,m ) in case iii is 0.625 mm, we have l x,m < λ y and λ y ≪ λ x , thus, l x,m < λ x . This same analysis is valid for case iv, where the length of the vertical regions is smaller than that in case iii but the waves' azimuthal wavelength is almost the same.
Accordingly, the reconstructed 2D fluctuation maps for cases iii and iv (the third and fourth row in figure 11) were found to be in remarkable agreement with the similar maps from the 2D simulation in [16] in terms of the azimuthal wavelength and the tilt in the waves before the channel exit. Moreover, in the 2D simulation, the instabilities were observed to become almost purely azimuthal downstream of the channel exit plane, which is also seen in the fluctuation maps of cases iii and iv.
Finally, concerning case ii (N = 10, M = 20), we notice that the axial wavenumber of the azimuthal waves is partially captured but the axial extent of the vertical regions has not been small enough in this simulation case to allow for an accurate resolution of k x .
To derive the wavenumber and frequency contents of the azimuthal fluctuations resolved by the Q2D simulations, we have used MATLAB's [26] built-in 2D fast Fourier transform (fft2) function. The 2D FFT is applied to the last 10 µs (i.e. from 20 to 30 µs) of the spatiotemporal signals of the azimuthal electric field. We have plotted the signals' frequency versus the azimuthal wavenumber to illustrate the numerical dispersion of the waves which is then compared against the theoretical dispersion relation of the ion acoustic instabilities in the laboratory's reference frame [27] where, ω R is the real frequency, ⃗ k is the wavevector, and k y is the azimuthal wavenumber. Also, λ D , C s and ⃗ V di are, respectively, the Debye length, the ion sound speed, and the ions' drift velocity. The first term on the right-hand side of equation (16) ( ⃗ k · ⃗ V di ) describes the convection of the waves with the ion flow [27,28]. In the dispersion plots shown below, however, we have ignored the convection term in the theoretical dispersion relation and have additionally considered the onesided dispersion comprising the positive frequencies only, i.e. Figure 12 illustrates the dispersion plots of the azimuthal electric field fluctuations at three axial locations from the 'high-order' Q2D simulations. In the dispersion plots of figure 12, the frequencies and wavenumbers are normalized, respectively, with respect to the ion plasma frequency (ω pi ) and the Debye length. The yellow curves in the plots of figure 12 correspond the local theoretical ion acoustic dispersion relation at each axial location, or more accurately, within the vertical region M whose mid-plane falls at the specified axial location. The green curves in the middle-and rightcolumn plots of figure 12 represent the ion acoustic dispersion relation corresponding to x = 0.4 cm.
Looking at the dispersion plots, we notice that, across all 'high-order' simulation cases, there is only a minor variation in the loci of the dominant frequencies and wavenumbers from the axial location of 0.4 cm to 2 cm. The variation in the wavenumbers is mostly from larger to smaller values when moving from inside the channel (x = 0.4 cm) toward the nearplume. In addition, the waves show a broader dispersion in the near-plume region.
The increasing deviation observed between the local ion acoustic dispersion relations and the one corresponding to the upstream is the evidence of the downstream convection of the azimuthal instabilities by the ion beam, as also reported in [16,28]. This deviation between the local dispersion relation and the upstream one is essentially due to the axial component of the term ⃗ k · ⃗ V di , i.e. k x V di,x , in the theoretical modified ion acoustic dispersion relation (equation (16)). In this respect, the increasing deviation occurs due to the accelerating ion beam.

Analysis of the sensitivity of the reduced-order axial-azimuthal simulations
In this section, we present the results from a series of Q2D simulations that we carried out to evaluate the sensitivity of the reduced-order scheme to the initial number of macroparticles per cell, and to assess the influence that decomposing the domain only along the axial direction has on the simulations' predictions.

Effect of the number of macroparticles per cell
The number of macroparticles per cell is a measure of the statistical quality of sampling the local distribution function in each computational cell. To assess the effect of this parameter on the reduced-order simulation predictions, we have selected the multi-region simulation case with N = 20 and M = 40 (case iii) because its results are closely similar to the full-2D ones whereas it is computationally less expensive than case iv.
We carried out several simulations with the initial number of macroparticles per cell along the azimuthal direction (N ppcy ) varying from 100 (i.e. the nominal value) to 50, and from 100 to 200 and 400. Figure 13 shows the axial profiles of the timeaveraged plasma properties from the simulations with various N ppcy values. Looking at the plots, it is noticed that the results are only marginally sensitive to the value of the N ppcy . The observed sensitivity of the Q2D simulations to N ppcy value was found to be the same as the sensitivity of the 2D reference simulation [15,16].

Effect of a single vs. multiple regions along the azimuthal direction
The current understanding of the plasma behavior in the axial-azimuthal coordinates of Hall thrusters indicates that the axial electric field does not have a notable gradient along the azimuthal direction. As a result, one may wonder if it might be actually necessary to divide the azimuthal coordinate into multiple regions in the reduced-order axial-azimuthal simulations, since such a division would be essentially intended to resolve the azimuthal gradients of the axial electric field. The results and analysis here are aimed at answering this question.
To this end, we have carried out dedicated 'high-order' multi-region Q2D simulations with only a single horizontal region along the azimuth.  simulation capturing all main features of the reconstructed distributions from the N = 20 simulation.
The above results underline the fact that, as initially claimed, the gradient of the axial electric field along the azimuthal direction is rather negligible, and a simulation that does not resolve this gradient provides predictions that compare very well with the simulations that do resolve this aspect. Consequently, we can state that, in general, an a-priori understanding of the physics of a specific problem, if available, can inform the decision concerning the degree of domain decomposition (i.e. the number of regions) such that the computational gain from the reduced-order scheme could be maximized.

Conclusions
In this article, we presented a mathematically rigorous dimensionality-reduction technique for Poisson's equation that enables a generalized reduced-order PIC scheme. It was underlined that the main difference between this novel scheme and a conventional PIC is related to the calculation of the electric potential using the RDPS. In this respect, as the generalizability of the RDPS was demonstrated in several standalone Poisson problems, the resulting reduced-order PIC scheme has, in principle, a broad applicability as well. Therefore, this effort significantly advances the idea of the 'pseudo-2D' PIC, introduced in our previous publication [10], by addressing the shortcomings that we had identified for the preliminary pseudo-2D scheme, namely, its lack of a well-founded formulation for the splitting of Poisson's equation, and the very limited generalizability either with respect to the order of approximation of the 2D problem or applicability to various simulation cases and physical configurations. Consequently, the concept of the reduced-order PIC scheme and its numerical implementation is demonstrated in this work to be mature enough to allow progress toward a computationally efficient, high-fidelity quasi-3D simulation tool for the study of the physics in Hall thrusters and potentially other low-temperature plasma systems.
Indeed, we showed in this paper that, in a 2D axialazimuthal simulation configuration and compared to the results from a corresponding reference 2D case, as we increase the order of approximation of the 2D problem by increasing the number of regions, the predictions of the Q2D simulations converge to those from the 2D reference simulation. We demonstrated that Q2D simulations with computational costs in the range of 5%-15% of a 2D simulation can reproduce the time-averaged plasma properties with an error of about 10%-2% while resolving the same features and characteristics observed for the azimuthal waves and their 2D distribution in full-2D PIC simulations.
Moreover, it was highlighted that even the sensitivity of the 'high-order' Q2D simulations to the number of macroparticles per cell is similar to that observed for the 2D codes. In addition, we showed that the azimuthal gradients of the axial electric field are rather negligible in the axial-azimuthal configuration, which can obviate the need for the decomposition of the domain along the azimuthal direction. This, in turn, yields additional computational advantage compared to Q2D simulations with the domain decomposition along both the axial and azimuthal coordinates.
Even though not presented in this article, the 'low-order' Q2D simulations, which offer about two orders of magnitude reduction in the computational cost, have been also verified in the 2D axial-azimuthal configuration, similar to the effort reported in [10] using the preliminary pseudo-2D PIC. We have observed that not only the same results reported in that article can be faithfully reproduced but also the improvements in the formulation of the reduced-order PIC scheme, discussed in this work, have notably augmented the predictive capability of the scheme in this low-computational-cost limit, enabling the approach to recover more accurately the results from the reference full-2D PIC simulations in several conditions. Interested readers are kindly referred to [29] where the detailed results of the verification of the generalized reduced-order PIC scheme at the low-order approximation limit are available.
The present effort serves as a conclusion for the verification of the reduced-order PIC scheme in the axial-azimuthal configuration. The verifications of the scheme in the radialazimuthal and axial-radial configurations of Hall thrusters are being pursued, on the path to the realization of a reliable quasi-3D kinetic simulation code. Also, as future work, we intend to verify the capabilities of the reduced-order scheme in a wider range of physical setups and applications by carrying out reduced-order simulations in cases representative of other low-temperature plasma configurations.

Data availability statement
The data cannot be made publicly available upon publication because the cost of preparing, depositing and hosting the data would be prohibitive within the terms of this research project. The data that support the findings of this study are available upon reasonable request from the authors. in equations (21) and (22) for each 2D horizontal and vertical computational element, respectively, In equations (21) and (22), Q x and Q y are the total charge averaged over the y-and the x-direction, respectively, and ρ x and ρ y are the corresponding charge densities.
As we increase the number of regions, their lengths L y,n and L x,m decrease. In the limit of having as many regions as there are 2D computational cells, the length L y,n becomes equal to ∆y and the length L x,m becomes equal to ∆x. In this situation, the two computational elements shown in figure 16 effectively collapse on each other for each computational node with positions x and y within the 2D plane. Consequently, ρ x and ρ y become equal to ρ (x, y), which is the charge density at the computational node with coordinates (x, y). Accordingly, in the limit, the Gauss' law for any computational element is the summation of equations (21) and (22), which can be written as below by substituting the definition of each electric field component from figure 16, Dividing both sides of equation (23) by 2 and by ∆x ∆y, it is seen that the resulting equation, equation (24), becomes identical to the 2D Poisson's equation as written in equation (19). This proves the convergence and, thus, the mathematical consistency of the dimensionality-reduction formulation, A.2. Sensitivity of the reduced-order PIC scheme to the size of the moving-average window in the RDPS As described in section 2.2, the implementation of the RDPS in the PIC code involves modifying slightly the coupled system of 1D ODEs given by equations (6) and (7) such that some of the cross-derivative terms are obtained from the previous timestep of the simulation. It was mentioned in section 2.2 that, before calculating these derivative terms at each timestep, we applied a moving-average filter with a certain window size on the potential solutions from the previous timestep. Accordingly, it is important to assess the impact of the size of the adopted window on the Q2D simulations to ensure that this numerical solver parameter does not have a major influence on the results. In this regard, figure 17 shows the axial profiles of the time-averaged plasma properties from the multi-region (N = 5, M = 10) Q2D simulations for various moving-average window sizes. It is observed that the plasma distributions are very closely comparable for all window sizes, with the plasma profiles from the simulations with the window size of 10 (the nominal value), 20 and 40 being almost identical. Of course, the simulation case with the window size of 5 predicts a minimally higher electric field peak and a slightly lower ion number density peak, but its predicted profiles are overall in line with the other cases. As a result, the reduced-order simulations are seen to be almost insensitive to the moving-average filter's window size, particularly when the size of the window is equal to or greater than 10. In a more general sense, the results from the reduced-order simulations are not affected by the window size once the application of the movingaverage filter effectively smooths out the noise on the potential solutions.