Toward Realistic Solar Flare Models: An explicit Particle-In-Cell solver in the DISPATCH framework
Michael Haahr Boris V. Gudiksen Åke Nordlund
 




   

   An explicit Particle-In-Cell solver in the DISPATCH framework

   

    Rosseland Centre for Solar Physics, University of Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo, Norway

    michhaa@uio.no
         Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo, Norway 

         Centre for Star and Planet Formation, Niels Bohr Institute, University of Copenhagen, Øster Voldgade 5-7, 1350 Copenhagen, Denmark

             
             
             
             

   
   


 
  Simulating solar flares, which involve large-scale dynamics and small-scale magnetic reconnection, poses significant computational challenges.
  
   This study aims to develop an explicit Particle-In-Cell (PIC) solver within the DISPATCH framework to model the small-scale kinetic processes in solar corona setting. This study in the first in a series with the ultimate goal to develop a hybrid PIC-MHD solver, to simulate solar flares.
  
   The PIC solver, inspired by the PhotonPlasma code, solves the Vlasov-Maxwell equations in a collisionless regime using explicit time-staggering and spatial-staggering techniques. Validation included unit tests, plasma frequency recovery, two-stream instability, and current sheet dynamics.
  
   Validation tests confirmed the solver's accuracy and robustness in modeling plasma dynamics and electromagnetic fields.
  
   The integration of the explicit PIC solver into the DISPATCH framework is the first step towards bridging the gap between large and small scale dynamics, providing a robust platform for future solar physics research.

   

   Toward Realistic Solar Flare Models
    M. Haahr
          1,2
          B. V. Gudiksen1,2
          Å. Nordlund 1,3
          
    September 9, 2024
===========================================================================================



§ INTRODUCTION

Simulating solar flares, encompassing both the macroscopic evolution and the microscopic details of current sheets that form during magnetic reconnection events, remains an unsolved computational challenge. According to recent work 

<cit.>, no computational tool currently bridges the gap between the large-scale dynamics of solar flares and the intricate processes within magnetic reconnection regions. 


Solar flares are generally accepted to result from the release of magnetic energy in magnetic reconnection regions <cit.>. Within these regions, some particles are accelerated to non-thermal, relativistic speeds. These high-energy particles exchange energy with the surrounding plasma. The non-local effects of these particles influence plasma behavior on large scales, which in turn affects the small scale dynamics. This interplay between large and small scales is critical for understanding the initiation and development of solar flares. 

Kinetic approaches, particularly the Particle-In-Cell (PIC) method  <cit.>, have been extensively applied to study reconnection at microscales. Despite their utility, these studies often resort to simplified, idealized configurations, such as the 2D Harris sheet model <cit.>, due to computational constraints. Given that real-world solar flares can extend over vast distances, employing PIC methods exclusively becomes impractical. Although there have been efforts to combine Magnetohydrodynamics (MHD) and PIC techniques <cit.>, challenges persist in accurately capturing the multiscale nature of solar flares. A hybrid modeling approach that can seamlessly integrate different temporal and spatial scales is thus paramount for realistic simulations.





Handling the wide range of temporal and spatial scales in solar flare simulations, while maintaining computational performance, is a significant challenge. Developing a sophisticated multi-scale, multiphysics solver from scratch would be daunting. Fortunately, the DISPATCH framework <cit.> provides a robust foundation for handling the scales required to simulate solar flares. DISPATCH organizes the simulation domain into "patches," each updated semi-autonomously. This localized data management enables each patch to function as a unique type of solver within a simulation, facilitating simulations with multiphysics solvers. DISPATCH has already proven capable of modeling complex, multiscale environments, such as the entire solar convection zone <cit.>.



This paper marks the initial phase in developing a hybrid solver within the DISPATCH framework by integrating a PIC solver inspired by the PhotonPlasma code <cit.>. We detail the process of embedding the PIC solver into DISPATCH, noting the essential modifications made from the original PhotonPlasma code to ensure compatibility and effectiveness within the DISPATCH architecture. The PIC solver is tailored for future integration with MHD simulations. Through this integration, we aim to bridge the gap between large-scale flare dynamics and the detailed kinetic processes occurring within reconnection zones, thereby contributing to a more comprehensive understanding of solar flare mechanisms.



§ METHODS




 §.§ Governing Equations



The PIC solver implemented in this study solves the Vlasov-Maxwell system of equations. 

For the sake of reproducibility and clarity, we outline only the most important equations here. For an exhaustive description, readers are referred to the original description of the PhotonPlasma code <cit.>. In our notation, we use bold letters to indicate vectors. The Vlasov-Maxwell system is represented as follows:

    ∂𝐟^s/∂ t  + 𝐮·∂𝐟^s/∂𝐱 + q^s/m^s (𝐄 + k_F 𝐮×𝐁) ·∂𝐟^s/∂ (𝐮γ) = C

Here, 𝐟 is the distribution function, which is a function of time, velocity (𝐮), and space (𝐱). s denotes particle species (e.g., electrons, protons), γ = 1/√(1 - (u/c)^2) the relativistic Lorentz factor and C the collision operator. Finally, 𝐄 and 𝐁 are the electric and magnetic fields respectively. For this study, we focus on the collisionless regime, setting C=0. 

The electric and magnetic fields' evolution follow the Maxwell equations:

    ∂𝐄/∂ t = k_E/k_B∇×𝐁 - 4π k_E 𝐉


    ∂𝐁/∂ t = - 1/k_F∇×𝐄,


    ∇·𝐁 = 0,


    ∇·𝐄 = k_E 4 πρ_c,

where 𝐉 and ρ_c represent the current density and charge density, respectively. 

The constants k_F k_E, and k_B ensure the equations' consistency across different unit systems, as detailed by <cit.>. We show the values corresponding to SI, CGS, and HL units in table <ref>. The use of this unit-agnostic formulation is commented on in section <ref>

The PIC solver approximates the Vlasov equation (eq:Vlasov) by sampling the phase space with "macro particles", each representing a collective of real particles, as discussed in works such as <cit.> and <cit.>. Macro particles are assigned a weight that represents the number of real particles they sample. The macro particles are then treated as individual particles, with their velocities influenced by the Lorentz force:



    d (𝐮_p γ)/dt = q_s/m_s (𝐄 + k_F 𝐮_p ×𝐁),

where their positions evolve as per:

    d 𝐱_𝐩/dt = 𝐮_p.

Here 𝐮_p and 𝐱_p is the velocity and position of each individual macro particle. q_s and m_s is the charge and mass of the particle species.

Interpolating electric and magnetic fields to macro particle positions require assigning specific shapes to these particles, affecting field interpolation. Following PhotonPlasma's methodology, we have implemented various shape functions, including Nearest-Grid-Point (NGP), Cloud-in-Cell (CIC), Triangular-Shaped-Cloud (TSC), and Piecewise-Cubic-Spline (PCS), noting that NGP can introduce instabilities <cit.>. We use 'CIC' interpolation and 2nd order derivatives as the default in our solver. 

Determining properties such as charge density from particle distributions uses similar shape function, and we refer to this as 'particle deposition' in this article. The interpolation and deposition shape functions can be defined independently as simulation input parameters.






 §.§ Transparent Scaling



  §.§.§ Scaling

Scaling of simulation parameters plays a critical role in ensuring the stability and optimizing the performance of simulations. Although frequently mentioned in scientific literature, the detailed methodologies for scaling are often not fully elaborated, posing challenges to reproducibility. In response to this, we aim to clearly outline the scaling techniques employed in our study.

Within the DISPATCH framework, we have incorporated a scaling module to handle scaling for SI, CGS, and HL unit systems. This module relies on three primary parameters for scaling: mass density, length, and time. The scaling of all subsequent units stems from these foundational parameters. 

Often, variables of interest such as the proton mass, the elementary charge, and the speed of light, are set to 1 in code units. Along with this, pressure, density, and field strength are usually given in code units as well, making it difficult to recover the 'real' physical units. For ease of physical interpretation and to ensure that input parameters remain intuitive, we have based all scaling on real physical units. As such, the code is designed to accept input parameters in familiar terms, including common length scales, magnetic field strengths, number densities, and temperatures.



  §.§.§ Fudging


A common practice in PIC simulations is the "fudging" of physical constants to adjust time and length scales appropriately. This process, although critical for the fine-tuning of simulations, is frequently glossed over in the literature, leaving a gap in understanding regarding the specific adjustments made. Our goal is to demystify this aspect by ensuring complete transparency in our scaling process. 

In our simulations, we permit the modification of three key constants: the speed of light, electron mass, and elementary charge. The modification of the speed of light is facilitated through changes to the vacuum permeability, μ_0, or  vacuum permittivity, ϵ_0. Specifically, in CGS and HL units where μ_0 is dimensionless, alterations to the speed of light are by adjusting ϵ_0. 





  §.§.§ Unit-Agnostic Equations



The Vlasov-Maxwell system of equations varies depending on the unit system used. Often, the unit system is not explicitly stated and the equations written in a 'natural units' where factors such as 4π are omitted. Gaussian and Heaviside-Lorentz (HL) are examples of 'natural units' and of often both referred to as CGS. Although units may appear identical in Gaussian and HL units their actual values differ. For instance, a magnetic field strength of 50G in Gaussian units corresponds to approximately 14.1G in HL units, differing by a factor √(4π). While an experienced reader may deduce the underlying system, this can be problematic for inexperienced readers. This can make it challenging to reproduce results if the exact unit system is not specified.

'Fudging' the physical constants also alters the equations. When ambiguous fudging is combined with an ambiguous unit system, it presents a complex puzzle for the reader to solve.

To assist the reader, we have written all equations in a unit-agnostic form. Constants are included in the equations, and these constants vary depending on the unit system and the 'fudging' used, as shown in Table <ref>. This approach allows us to use the same equations in both 'real' units and code units, simplifying implementation and improving the clarity of our code.






 §.§ Numerical Approach




  §.§.§ DISPATCH


A fundamental understanding of the DISPATCH framework is needed to understand some of our choices made when implementing the PIC code into DISPATCH. While a detailed description  of DISPATCH can be found in <cit.>, we here summarize the key aspects relevant to our PIC solver integration.

DISPATCH organizes the simulation domain into distinct sections referred to as "patches". These patches are updated semi-autonomously, interacting with adjacent patches primarily to update ghost zones. This architecture offers substantial benefits, specifically the ability for each patch to proceed with its own timestep. This feature significantly boosts the potential for simulation efficiency and enables the simultaneous resolution of diverse timescales within a single experiment.

The framework's design, focusing on localized data management, inherently limits the applicability of global convergence schemes like implicit methods. Furthermore, DISPATCH imposes specific restrictions on the communication between patches, designed to manage the complexity of data interactions efficiently. 



Adapting to DISPATCH's focus on localized computations, our implementation utilizes an explicit PIC solver, staggered in both time and space. This choice is compatible with DISPATCH's decentralized data handling paradigm, ensuring computational efficiency and precision while navigating the framework's structural nuances.



  §.§.§ Discretization






For the numerical stability and reliability of our explicit computational scheme, we incorporate a time-staggering method akin to the leap-frog integration technique, as visually represented in fig:time_stagger. The leap-frog scheme (<cit.> section 9.6, <cit.> p. 56) is renowned for its efficacy in maintaining stability during oscillatory motion— a common occurrence within PIC simulations. The electric field is staggered by a full timestep,  Δ t in the past. Magnetic field, particle velocities, and current density are staggered backwards by 0.5Δ t. Particle position is centered in time.


Spatially, we stagger mesh variables in accordance with a Yee lattice setup (<cit.>), to align with DISPATCH standards. This spatial organization involves positioning magnetic fields on the cell faces, electric fields and current densities are staggered along cell edges. Consistent with DISPATCH conventions, mass density, momentum, and energy metrics are centered within each cell. Uniquely, charge density is allocated to the lower corner of each cell, as shown in fig:space_stagger. 

This arrangement of variables in time and space minimizes the interpolation needed, thereby improving the simulation's computational efficiency and accuracy.



In our simulation, both electric and magnetic fields are interpolated from the mesh to each particle's position. Interpolating fields to particle positions involves considerations of both time and space. As outlined in fig:time_stagger, the magnetic field requires time interpolation between the timestep i, and timestep i+1. However, with varying timesteps this introduces the potential for this time interpolation to become desynchronized. Within the DISPATCH framework, each patch dynamically sets its timestep based on the local maximum velocity to comply with the Courant condition <cit.>. This approach is applied to all patches in the simulation, including the PIC solver patches. Due to the PIC solver's focus on electromagnetic phenomena and the need to accurately simulate the propagation of electromagnetic waves at the speed of light, the "maximum velocity" in PIC patches defaults to the speed of light. This condition results in a uniform timestep across these patches. Consequently, time interpolation for the magnetic field is simplified to averaging B^i and B^i+1, as outlined in fig:time_stagger.

Spatially, after addressing time interpolation, fields are then interpolated to the exact spatial locations of particles. This spatial interpolation relies on the shape functions introduced in Section <ref>.




Particle deposition into current density, mass density, and other mesh variables happens after the particles have been updated, as described in section <ref>. This process ensures that changes in particle distribution are accurately reflected in the mesh. Consequently, electromagnetic field calculations can incorporate the most current state of the plasma. We optionally allow the shape function for deposition to differ from the one used for interpolation, but as default the shape functions are identical.






  §.§.§ Solving Maxwell's Equations

The update cycle for each patch begins with updating electric and magnetic fields. The integration of Maxwell's equations in our solver uses an explicit method, sequentially updating the electric and magnetic fields.
The updating process initiates by updating E^i to E^i+1, (fig:time_stagger), using eq:Ampere_law. This calculation uses the magnetic field, B^i,


alongside the current density, J^i. The magnetic field is then updated in line with eq:Faraday_law, using the just-updated electric field.



Accumulation of errors over time can require the use of 'cleaners' for both the electric and magnetic fields. Within the DISPATCH framework, divergence cleaners for the magnetic field are employed to uphold eq:div_B in the guard zones, where time-interpolation in general does not conserve div(B) exactly. A similar cleaning strategy is applied to the electric field to ensure it meets the requirements of eq:div_E. This cleaning procedure begins by identifying the divergence error, ϵ, through the calculation of the electric field's divergence, from which the charge density is then subtracted:

    ϵ = ∇·𝐄 - k_E 4 πρ_c.

A subsequent Poisson filtering step addresses only the errors local to each patch to compute 
Φ_ϵ from 

    ΔΦ_ϵ = ϵ ,

after which the electric field is updated to a 'cleaned' state:

    𝐄_clean = 𝐄 - ∇Φ_ϵ.


The implementation allows users to decide the frequency of the electric field cleaning step and whether to use it at all. 





  §.§.§ Particle Movement

Next in the update cycle, particle velocities are updated, affecting both particles in the "inner" domain and those within the ghost zones. As highlighted earlier, the electric and magnetic fields are interpolated to the particles' positions before their velocities are updated.

The update of velocities uses the Lorentz force equation (eq:Lorentz_foce). <cit.> compares several different explicit approaches to implement this numerically. We have chosen to use the Vay particle pusher algorithm <cit.> due to its ability to correctly 'cancel' field contributions when both electric and magnetic fields act on particles. Following the velocity updates, we calculate the energy for each macro particle. Then, their positions are updated according to eq:particle_position.

With updated positions and velocities, the simulation is prepared for the next and final phase – deposition of particle data back to the mesh. This stage involves computing the current density, which is essential for updating the electric field in the next timestep. Additionally, other mesh variables, such as mass density and bulk momentum, are computed from particle deposition. 






  §.§.§ Current Density Calculation


Current density calculation is a critical step in ensuring the accuracy of PIC simulations. Often, PIC codes employ a 'charge conserving' method. For example, the PhotonPlasma uses an extension of the charge conserving method developed by <cit.>. This method involves calculating the divergence of the current density by examining how the deposition weights for each macro particle change between timesteps. Following the determination of the divergence, the actual current density is obtained through a prefix sum operation. 

In DISPATCH, such prefix sum operation would require communication across patches to align ghost-values with neighbor patches. The current communication scheme does not allow this, and we therefore opted for a simpler and more direct current density calculation: 


    𝐉 = ∑_s q^s ∑_i w_i 𝐯_𝐢 W(x_p - x_c),

where W(x_p - x_c) represent deposition weight function, w_i the weight of the macro particle, 𝐯_𝐢 the velocity of the macro particle, and q_s the charge of the particle species. This method calculates current density directly from the macro particles' properties.

While not exactly charge conserving, down to each macroparticles gradual contribution, the method is nevertheless on average charge conserving, with errors expected mainly on cell scales, where E-field cleaning (Eqs. <ref>-<ref>) can easily take care of making the electric field exactly consistent with actual charge distribution.

 


 §.§ Particle Sampling



  §.§.§ Initial Condition Sampling




The initial sampling of particles within our simulation is based on mesh variables such as mass density, current density, bulk velocity, and temperature. These mesh variables, critical for defining initial conditions, are aligned with those used in MHD simulations, to streamline PIC-MHD integration efforts.

To distribute particles evenly at the start, we segment each cell into a predefined number of sectors. The number of sectors, n_sec, varies with the dimensionality of the simulation:

    
  * In 3D, n_sec = ⌊ n_pc^1/3⌋,
    
  * In 2D, n_sec = ⌊ n_pc^1/2⌋,
    
  * In 1D, n_sec = n_pc,

where n_pc denotes the number of particles per cell.

For example, in a 2D setup with 10 particles per cell, n_sec = 3, leading to each cell being divided into 3 × 3 = 9 sectors. One particle is placed randomly within each sector, with any additional particles randomly located within the entire cell, as depicted in fig:position_sampling.

After establishing positions, particles' velocities are sampled from a Maxwellian distribution. The initial temperature for this distribution is either provided directly or derived from the internal or total energy, depending on the initial conditions.


Additionally, we incorporate bulk and current velocities into the velocity assignment for particles, similar to the approach described by Daldorf et al. (2014):

    𝐯_𝐞𝐥𝐞𝐜 = 𝐮_𝐛𝐮𝐥𝐤 + m_i/q_e𝐉/ρ,


    𝐯_𝐩𝐫𝐨𝐭𝐨𝐧 = 𝐮_𝐛𝐮𝐥𝐤 + m_e/q_i𝐉/ρ,

where 𝐮_𝐛𝐮𝐥𝐤 represents the bulk velocity, 𝐉 is the current density derived from ∇×𝐁, ρ is the mass density from grid, and m_i, q_i, m_e, q_e are the masses and charges of protons and electrons, respectively. 







  §.§.§ Resampling



Particles within the inner domain of a patch can freely move between cells. Over time, this freedom can lead to significant discrepancies in the number of particles per cell, with some cells having too many particles and others too few. Similarly, particle distribution discrepancies can occur at the patch level. To maintain a roughly constant number of particles per cell and ensure workload balance across patches, we've implemented a resampling strategy that involves "splitting" particles in underpopulated cells and "merging" particles in overcrowded ones.


Splitting, or up-sampling, is straightforward. The heaviest particle in a cell is divided into two particles, each retaining the original velocity but carrying half the weight. To prevent creating identical particles, their positions are slightly offset. This offset is randomly determined in each dimension from a Maxwell distribution with σ = 0.1, with one particle receiving half the offset subtracted and the other half the offset added. This process can then be repeated until the desired number of particles per cell is achieved.


Merging, or down-sampling, presents more complexity, as it can potentially violate conservation laws for mass, energy, and momentum. Drawing on strategies discussed by <cit.>, we've adopted the 'leveling' approach from their work and introduced a 'close-phase' strategy. For the latter, we evaluate particle velocities within the same cell using the Euclidean distance metric, merging the two particles closest in velocity space into one. This new particle's energy is then compared to the total energy of the original pair. If the energy discrepancy falls within an acceptable tolerance, the merge is finalized; otherwise, the particle is discarded. This process repeats until the target particle count is achieved or all potential merges are explored without success due to excessive energy differences.

The frequency and activation of resampling are left to user discretion, offering flexibility in simulation management.
















    
   
  
 















 §.§ Optimization

Optimization is crucial for our goal of developing a comprehensive and realistic solar flare model across multiple scales. We've introduced several optimizations to boost the performance of our simulation framework, preparing for future enhancements.

Our simulation adopts a Particle-in-Patch (PIP) approach, similar to DISPATCH's strategy of updating patches instead of individual cells. This method optimizes cache usage and enables effective vectorization of computations. In our PIC solver, each patch hosts several arrays of particle information, facilitating vectorization. Though these arrays are too large for cache, making them less than ideal for traditional computing, their structure is ideal for GPU execution. GPU acceleration remains a key development goal, given its potential to significantly enhance simulation performance.

A fundamental part of our optimization strategy involves sorting particles to improve memory access patterns. When particles spatially close to each other are also adjacent in memory, the efficiency of interpolation and deposition operations increases markedly. We store particle positions in two arrays: an integer array for cell location and a floating array for position within the cell. This dual-array setup simplifies sorting, allowing for integer-based comparisons that improve performance and numerical precision.

The choice of sorting algorithm impacts performance. We have implemented both insertion sort and quick-sort algorithms <cit.>. Given the minimal expected changes in particle ordering between timesteps, insertion sort seemed the likely choice. However, initial results indicate better performance with a complete re-sort using quick-sort, even when sorting every timestep. 




Interpolation and deposition are, especially for higher-order schemes, the most demanding processes in our PIC solver, due to inefficient memory access patterns when accessing cells in higher dimensions. Sorting particles addresses this issue to some extent. In our effort to optimize the interpolation routine, we found a more efficient method than interpolating for each particle during the particle pusher step. Instead, interpolating field values for all particles prior to the particle pusher step and storing these values in temporary arrays performed better.






§ CODE VALIDATION

To validate the PIC solver integrated within the DISPATCH framework, we have devised a comprehensive suite of test cases that progressively increase in complexity. Initially, we conduct a series of unit tests aimed at examining specific equations governing particle motion and field evolution in isolation, i.e., without any feedback between particles and fields. Subsequently, we introduce complexity by incorporating particle-field feedback mechanisms in both 1D and 2D test scenarios. All tests described herein uses CGS units for consistency, and we did not activate resampling to preserve the integrity of conserved variables. 



 §.§ Unit tests



  §.§.§ Constant Particle Motion




This test, designed to exclusively evaluate eq:particle_position, initiates with 18 particles, each with distinct velocities across every direction. Both electric and magnetic fields are nullified, preventing any particle-induced field alterations. The experiment is conducted within a 3D domain, comprising a grid of 3x3x3 patches, each containing 8x8x8 cells, to examine the solver's three-dimensional computational integrity. The particle trajectories are tracked to detect any discrepancies as they traverse patch boundaries. fig:unit_test illustrates the outcome of this test, affirming that all observed errors remain within the bounds of numerical precision. 
 





  §.§.§ Constant Particle Acceleration




This test specifically aims to validate eq:Lorentz_foce in the absence of a magnetic field, focusing on the effects of an electric field on particle acceleration. We initialize the experiment with seven particles, each moving along one of the Cartesian axes (X, Y, and Z) in positive and negative direction. The seventh particle has initial zero velocity. To evaluate the solver's accuracy in calculating particle acceleration due to an electric field, we conduct three separate runs. In each run, a constant electric field is applied along one of the axes.

The setup for this test mirrors the configuration used in the previous section, utilizing a 3D grid of 3x3x3 patches, where each patch consists of 8x8x8 cells. In each of the three runs, we analyze the change in proper velocity of the particle moving along the axis aligned with the electric field. The primary objective is to confirm that the solver accurately captures the expected velocity change between snapshots, with a focus on the incorporation of relativistic effects. All observed changes in velocity fall within the realm of numerical precision for both electrons and protons as shown in fig:const_acc_unit.



  §.§.§ E cross B drift



The E cross B drift test extends our exploration of eq:Lorentz_foce to scenarios involving both electric and magnetic fields positioned perpendicularly to each other. This orthogonal arrangement results in a drift velocity expressed as v_drift = c E_y/B. Specifically, with an electric field oriented in the y-direction (E_y) and a magnetic field in the z-direction (B_z), the theoretical positions of particles can be determined as follows:

    x = x_0 + v_sign· r_L ·sin(ω· t) + c E_y/B_z t,


    y = y_0 + v_sign· q_sign· r_L · (cos (ω· t) - 1),

where r_L denotes the Larmor radius within the reference frame where E_y = 0 — namely, the frame moving at the drift velocity v_drift = c E_y/B_z. The variable v_sign represents the sign of the initial velocity in this specific frame, while q_sign indicates the charge sign of the particle, and ω denotes the gyro frequency.

fig:ExB_drfit presents the results from this test. Notably, the figure illustrates that the drift velocity error for electrons stays within the bounds of numerical precision up to t=40, equating to around 40,000 timesteps.






  §.§.§ EM-wave in Vacuum




This test assesses the solver's implementation of eq:Ampere_law and eq:Faraday_law in a vacuum scenario with zero current density (𝐉=0). Executed within a 1D setup, the experiment examines the propagation of an electromagnetic wave along a single axis. Separate runs were conducted for each of the three spatial directions for a comprehensive validation. Here, we detail the results for the z-direction. The theoretical model of the wave is described by the following set of equations:

    Bx    =    - B_0 cos(k z - ω t) 
    
        By    =    B_0 cos(k z - ω t) 
     
        Bz    =    0 
    
        Ex    =    B_0 cos(k z - ω t) 
    
        Ey    =    B_0 cos(k z - ω t) 
    
        Ez    =    0

where k represents the wavenumber and ω the frequency of the electromagnetic wave.

We initialize with k = 2 π· m_0/L, with L indicating the domain's extent and m_0 the designated number of wavelengths within this domain. The frequency, ω, is derived using ω = k · c.

The domain's span is set to 1000 cm, and m_0 = 2. Scales are adjusted so that the domain's length corresponds to 1 in code units, and similarly, time is scaled to set the speed of light, c, to 1 in these units. Various cell counts, including 250, 500, and 1000, were tested to gauge the solver's accuracy across different resolutions. fig:EM_wave displays the relative errors encountered in these tests, with derivatives of 2nd, 4th, and 6th orders being examined. Errors predominantly arise from phase shifts and amplitude variations, with phase-shift discrepancies dominating in higher-order cases, while amplitude errors dominate in the 2nd order outcomes.

Error convergence for each order of derivation was critically analyzed. The 2nd order tests showed a 'linear' reduction in errors on a log-log scale, with a slope approximating -2.2, aligning closely with the expected -2.0 slope indicative of its second-order nature. In contrast, 4th and 6th order tests revealed a diminished slope of -0.65, hinting at marginal gains from enhanced resolution. This outcome suggests that the 2nd order time integration of the leapfrog scheme emerges as a significant constraint, limiting the benefits of higher-order spatial resolutions.

To avoid indexing errors in proximity to patch boundaries requires increasing the number of guard zones when using higher-order derivatives, thus incrementally increasing the computational workload. Considering also that a change from 4th to 6th order derivatives does not markedly improve precision, while significantly inflating the computational demands, we conclude that lop-sided increase of only spatial order is not a good idea.



 §.§ Recovery of Plasma frequency



The plasma frequency is pivotal in characterizing plasma dynamics. Accurately capturing this frequency in simulations is crucial for validating the scalability and feedback mechanisms within our solver. We explored this through Fourier analysis in a controlled 2D setup, aiming to recover the plasma frequency.

The experimental setup involved initializing electrons and protons at a uniform density of 10^9  cm^-3 within a 5 cm× 5 cm domain. The domain is segmented into 5x5 patches with 10x10 cells in each for a total of 50x50 grid points. We initialized each species with a thermal velocity corresponding to a temperature of 1 MK. With these parameters the cell width is approximately one Debye length.

We scaled mass density to unity, time such that ω_e = 2 π, and did not scale lengths.

As depicted in fig:plasma_freq, the Fourrier analysis reveal distinct peaks corresponding precisely to the plasma frequency. Time has been scaled such that the theoretical plasma period is exactly 1 'code Hz'. The peak is most notable in the field component orthogonal to the testing plane. This test demonstrates the solver's adeptness at recovering the correct plasma frequency.



 §.§ Two-Stream Instability






Investigating the two-stream instability provides a rigorous assessment of the solver's capability to simulate plasma interactions and instability dynamics. Several different versions of the two stream experiment exists and have been tested  <cit.>. Our setup closely resembles <cit.>, but is notably different due to different scaling parameters. We initialize our setup with two counter-streaming electron populations against a backdrop of stationary ions. Electrons and protons are initialized with Maxwellian velocities at temperature T, and electrons additionally with a drift velocity of ± V_0.

Analyzing the cold case (T → 0) through the dispersion relation:

    0 = 1 - ω_e^21/2(1/(ω - k V_0)^2 + 1/(ω + k V_0)^2) - ω_p^2/ω,

yields the maximum growing wave mode at k_max≈√(3/8)ω_e/V_0, with growth rate γ≈ 0.35355 ω_e, when neglecting the contribution from ions. Delving into the warm case revealed the growth rate and mode size changes by less than 1% for v_e ≤ 0.05V_0, informing our temperature selection to closely model realistic cold-case plasma conditions.

The computational domain span 4.105 d_e where d_e = c/ω_e represents the electron skin depth. This domain size was chosen to accommodate two wavelengths of the maximum growing mode. We conduct a series of simulations with varying cell counts (n=[64,128,256,512]) and particles per cell (ppc=[64,128,256,512,1024]). For n=64 we used 4 patches with 16 cells in each. For the other cell counts we used 32 cells per patch and [4,8,16] patches, respectively. 




Physical variables were kept unscaled, and we normalized length such that d_e=1 in code units and ω_e=1 for time, rendering the speed of light unity in code units as well. We chose a drift velocity of V_0=0.2c, which sets the temperature T≈ 592,990 K for v_e = 0.05 V_0. Additionally, we used a number density of 10^10 cm^-3 and scaled mass, such that the mass density is 1 in code units. 

The simulations proceeded without field cleaners or particle resampling until t=40. fig:two_stream_evo illustrates the evolution of the linear phase and the initial emergence of the nonlinear phase, showing the formation of dual structures in phase space as anticipated from the experiment's design.

fig:two_stream_growth shows the electric field's growth rate for the dominant wavemode, illustrating consistency with the anticipated growth rates derived from linear theory. This alignment is consistent across various experimental configurations, underscoring the solver's robustness in accurately modeling the intricate processes underlying two-stream instabilities under a spectrum of simulation conditions.



 §.§ Current Sheet









The final test is the behavior of the physical effects under different parameter scaling scenarios. We chose a current sheet setup as such a scenario is inherently unstable and should quickly highlight any differences in evolution. Current sheet behavior is an open research topic with many authors investigating current sheet in space and solar contexts, wiht varying degree of realism (e.g. <cit.>, <cit.>, <cit.>, <cit.>). For our initial validation, we have chosen a simpler setup, inspired by <cit.>. The magnetic field configuration is defined by two hyperbolic tangent functions:


    B_x(y) = 
     
    - B_0 tanh( y + 3.6/0.5 · L)    if  y < 0, 
    
    B_0 tanh( y - 3.6/0.5 · L)    if  y ≥ 0,

where B_0 = 50 G represents the background magnetic field strength and L the initial width of the current sheets. We set L=r_p where r_p is the proton gyroradius. The simulation initializes electrons and protons with Maxwellian distributions at T=1 MK and uniform density n=10^9 cm^-3 within a domain sized [25.6,12.8]r_p. Periodic boundary conditions are applied along both the x and y axes. An initial velocity perturbation:

    v_y = v_0 ·sin( 2 π/L_x x ),

with v_0 = 0.001 c, stimulates the onset of reconnection.

Experimentation varied over cell counts ([128x64, 256x128, 512x256, 1024x512, 2048x1024]) and particles per cell ( = [32, 64, 128, 256]). For  = [128x64, 256x128] we used 16x16 cells per patch and [8x4,16x8] patches. For  = [512x256, 1024x512, 2048x1024] we used 32x32 cells per patch and [16x8, 32x16, 64x32] patches.

Four distinct 'fudging' strategies were applied to probe their influence on simulation behavior:




    
  * Pure Runs: Physical variables remained unaltered, providing a control scenario for comparison.
    
  * m400 Runs: Electron mass was scaled to set the mass ratio m_p / m_e = 400.
    
  * c10 Runs: Speed of light was reduced tenfold, achieved by augmenting ϵ_0 by 100.
    
  * c10m400: Combined the adjustments from m400 and c10 runs for cumulative effect analysis.


Modifying the elementary charge primarily reduces the gap between microscopic and macroscopic scales <cit.>. In a pure PIC simulation, changing the elementary charge only affects the length and time scales. When scaled to code units, there is no difference between runs with and without elementary charge scaling. This is unlike scaling the electron mass or the speed of light, which explicitly changes the ratios between parameters of interest, even in code units. Therefore, we did not make any adjustments to the elementary charge.

We kept scaling consistent across scenarios for analytical clarity. We scaled length such that r_p = 1, time such that plasma frequency ω_e = 2 π, and initial mass density so ρ_0 = 1 in code units for the non-fudged run.






The current density evolution depicted in fig:J_comparison illustrates the effects of electron mass and light speed adjustments on plasmoid dynamics. Particularly, the c10 and c10m400 runs showcase a prevalence of smaller-scale plasmoids, contrasting with the more homogeneous current dissipation in pure and m400 scenarios.

Figure <ref> illustrates convergence towards zero relative total energy error with increasing resolution, in alignment with theoretical expectations for a second order solver. A slope of -2.2 was found, similar to our findings from the EM-wave test.

The influence of the particle count per cell on electric field energy stabilization is investigated in fig:electric_energy, affirming the solver's stability under cell sizes below 10 Debye Lengths. Theoretically, Δ x < 3.4 λ_D must be satisfied to avoid self-induced instabilities <cit.>. The theoretical analysis assumes that the number of particles inside the Debye radius is large. In our experiment, we show that the solver is stable for as little as 2 particles per Debye area for 2D simulations.


This comprehensive suite of tests validates the solver's robustness in modeling the intricate dynamics of current sheet formation and magnetic reconnection, setting a solid foundation for expanding investigations into 3D simulations.



§ DISCUSSIONS AND CONCLUSION



 §.§ Speed of light and related issues

The integration of the explicit PIC solver within the DISPATCH framework marks our first step forward towards simulating the complex plasma dynamics of solar flares. This integration also highlights certain limitations inherent to the nature of explicit solvers, and how this plays out relative to the DISPATCH code framework. 

A primary challenge lies in the limitations on the timestep size by the speed of light; a consequence of the need to accurately model electromagnetic field propagation. This constraint, while potentially limiting in other, less extreme contexts, aligns well with the solver's intended application to scenarios involving relativistic particle acceleration, such as observed in solar flares.  One can hardly claim to realistically model such scenarios without embracing and accepting that one must deal with the full range of velocities.

Moreover, relativistic conditions offer other advantages, such as effectively reducing the impact of the difference of mass between electrons and ions, with both species moving at essentially the same velocities, despite large difference in mass.

The explicit solver's reliance on uniform timesteps might seem to diverge from DISPATCH's otherwise innovative approach of employing local time stepping. To some extent, this is an unavoidable consequence of the fundamental impact of the speed of light in relativistic dynamics—when information speed is actually essentially constant there is no advantage to gain from variations in speed, but the framework still allows for coping with large differences in time steps due to differences in spatial resolution, coming from fixed or adaptive mesh refinement.



 §.§ Exposing the consequences of `fudging'

As previously mentioned, the 'fudging' of physical constants is a key component in many PIC codes. However, the descriptions of these adjustments are often unclear. With our 'transparent scaling,' we address this issue by clearly outlining the modifications and their implications.

Some examples of the impact of fudging are examined in Section <ref>. We demonstrate how varying the speed of light and the mass of the electron alters the fundamental plasma scales. Specifically, reducing the speed of light by a factor of 10 increases the allowable timestep by the same factor. However, simultaneously, the relevant scales are reduced by approximately a factor of 10, and thus the benefits offered by this fudging might be partially or fully offset. As illustrated by fig:J_comparison, the results with and without fudging may differ dramatically, and thus the methodology requires very careful vetting of assumptions and consequences, as made possible by the clear distinction made here between fudging on the one hand, and scaling to code units and choice of unit system on the other hand.



 §.§ Conclusion

We have successfully introduced an explicit PIC solver into the DISPATCH framework, adapting and building upon the PhotonPlasma code's foundations. 

Our comprehensive validation process, encompassing both straightforward unit tests and more intricate simulations such as the double current sheet test, has rigorously confirmed the PIC solver's robustness and precision. These efforts underscore its capacity to simulate the nuanced interplay of particle dynamics and electromagnetic fields, offering a robust platform for cutting-edge research in solar and plasma physics.

The integration into the DISPATCH code framework not only broadens the solver's applicability across diverse simulation scenarios but also lays the groundwork for a hybrid PIC-MHD implementation. Such an implementation will greatly extend our ability to model solar flares as well as other complex astrophysical phenomena with unprecedented realism, by using the cheaper MHD representation where possible, switching to the much more expensive PIC representation only where needed. 


We would like to thank Troels Haugbølle for his invaluable discussions and insights regarding the PhotonPlasma code and plasma physics. We also wish to acknowledge Jan Truelsen for his assistance and contributions to the theoretical foundations of this work, particularly insight into the warm case analysis of the two-stream instability test. Additionally, we thank Andrius Popovas and Mikolaj Szydlarski for their guidance and support throughout the code development process.


This research was supported by the Research Council of Norway through its Centres of Excellence scheme, project number 262622.





plainnat  
  


