We introduce Astrophysical HybridKinetic simulations with the FLASH code inlineformulatexmath idTM0001 notationLaTeXtt AHKASHtexmathinlineformula  a new Hybrid particleincell PIC code developed within the framework of the multiphysics code FLASH. The new code uses a secondorder accurate Boris integrator and a predictorpredictorcorrector algorithm for advancing the Hybridkinetic equations using the constraint transport method to ensure that magnetic fields are divergencefree. The code supports various interpolation schemes between the particles and grid cells with postinterpolation smoothing to reduce finite particle noise. We further implement a inlineformulatexmath idTM0002 notationLaTeXdelta ftexmathinlineformula method to study instabilities in weakly collisional plasmas. The new code is tested on standard physical problems such as the motion of charged particles in uniform and spatially varying magnetic fields the propagation of Alfvn and whistler waves and Landau damping of ion acoustic waves. We test different interpolation kernels and demonstrate the necessity of performing postinterpolation smoothing. We couple the TURBGEN turbulence driving module to the new Hybrid PIC code allowing us to test the code on the highly complex physical problem of the turbulent dynamo. To investigate steadystate turbulence with a fixed sonic Mach number it is important to maintain isothermal plasma conditions. Therefore we introduce a novel cooling method for Hybrid PIC codes and provide tests and calibrations of this method to keep the plasma isothermal. We describe and test the hybrid precision method which significantly reduces by a factor inlineformulatexmath idTM0003 notationLaTeXsim 1.5texmathinlineformula the computational cost without compromising the accuracy of the numerical solutions. Finally we test the parallel scalability of the new code showing excellent scaling up to 10000 cores.
['Astrophysics - Instrumentation and Methods for Astrophysics', 'Astrophysics - Earth and Planetary Astrophysics', 'Physics - Computational Physics']
AHKASH a new Hybrid particleincell code for simulations of astrophysical collisionless plasma
{'AHKASH : a new Hybrid particle-in-cell code for simulations of astrophysical collisionless plasma': 'Radhika Achikanath Chirakkara,  \n1  \n★  \n1  \n,  \n2  \n<!-- image -->  \n<!-- image -->  \n1 Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia 2 Australian Research Council Centre of Excellence in All Sky Astrophysics (ASTRO3D), Canberra, ACT 2611, Australia  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': "Weintroduce A strophysical H ybridK inetic simulations with the flASH code ( AHKASH ) - a new Hybrid particle-in-cell (PIC) code developed within the framework of the multi-physics code FLASH . The new code uses a second-order accurate Boris integrator and a predictor-predictor-corrector algorithm for advancing the Hybrid-kinetic equations, using the constraint transport method to ensure that magnetic fields are divergence-free. The code supports various interpolation schemes between the particles and grid cells, with post-interpolation smoothing to reduce finite particle noise. We further implement a 𝛿 𝑓 method to study instabilities in weakly collisional plasmas. The new code is tested on standard physical problems such as the motion of charged particles in uniform and spatially varying magnetic fields, the propagation of Alfvén and whistler waves, and Landau damping of ion acoustic waves. We test different interpolation kernels and demonstrate the necessity of performing post-interpolation smoothing. We couple the TurbGen turbulence driving module to the new Hybrid PIC code, allowing us to test the code on the highly complex physical problem of the turbulent dynamo. To investigate steady-state turbulence with a fixed sonic Mach number, it is important to maintain isothermal plasma conditions. Therefore, we introduce a novel cooling method for Hybrid PIC codes and provide tests and calibrations of this method to keep the plasma isothermal. We describe and test the 'hybrid precision' method, which significantly reduces (by a factor ∼ 1 . 5) the computational cost, without compromising the accuracy of the numerical solutions. Finally, we test the parallel scalability of the new code, showing excellent scaling up to 10,000 cores.  \nKey words: methods: numerical - plasmas - turbulence", '1 INTRODUCTION': "Manyastrophysical plasmas, such as the solar wind (Verscharen et al. 2019), hot accretion flows in the galactic center (Event Horizon Telescope Collaboration et al. 2019), and the hot intracluster medium of galaxy clusters (Schekochihin & Cowley 2006; Kunz et al. 2022) are weakly collisional or collisionless plasmas, and their evolution is coupled with magnetic fields of dynamically important strength. The collisionality of a plasma can be determined from the ratio of the Coulomb mean free path ( 𝜆 mfp ) to the characteristic length scale of the system ( 𝐿 ). For example, the hot interstellar medium is characterised by 𝜆 mfp / 𝐿 ∼ 10 -3 -10 -2 , which can be suitably modelled as a collisional plasma (Ferrière 2020), where the magnetohydrodynamical (MHD) approximation holds reasonably well. For the solar wind, 𝜆 mfp / 𝐿 ∼ 1, for hot accretion flows, 𝜆 mfp / 𝐿 ∼ 1, and the hot intracluster medium, 𝜆 mfp / 𝐿 ∼ 0 . 1, making these plasmas weakly collisional or collisionless. Such astrophysical systems are interesting laboratories to study the nature of magnetized collisionless and weakly collisional plasma and it is important to develop numerical schemes that can suitably model such systems, which is the aim of this work.  \nMagnetohydrodynamic (MHD) codes have been used extensively to model magnetized astrophysical gas (Fryxell et al. 2000; Dubey  \nChristoph Federrath,  \net al. 2008a; Mignone et al. 2007; Stone et al. 2008; Brandenburg & Dobler 2010; Stone et al. 2020; Pencil Code Collaboration et al. 2021). Assuming gas as a conducting fluid, these codes solve the continuity equation, Navier Stokes equation, induction equation for magnetic field evolution and energy equation for a given equation of state. MHD numerical simulations have become an established method to model astrophysical plasma in the last few decades. However, when the plasma is weakly collisional or collisionless, the continuum limit of the fluid equations breaks down and MHD is no longer valid for such plasmas. In such situations, we resort to a kinetic treatment of the plasma (Kulsrud 2005), which is our approach for this work.  \nKinetic codes solve the Vlasov or the Vlasov-Landau equation (Birdsall & Langdon 1991), to study the evolution of the ion distribution functions, along with Maxwell's equations. Numerically this can be achieved in many ways, including solving the six-dimensional (three spatial + three velocities) equations for the distribution function (Ganse et al. 2023) or using a particle-in-cell (PIC) approach to sample the moments of the distribution function while using a threedimensional spatial computational grid (Spitkovsky 2005; Cerutti et al. 2013; Shalaby et al. 2017; Derouillat et al. 2018; Nättilä 2022). Here we focus on the PIC method. The PIC method can simulate the physics on the ion and electron time and length scales, but doing so makes this approach computationally expensive and challenging. Moreover, it can be very difficult to run these simulations for larger physically interesting time and length scales, as the scale separa-  \non between physical processes is often large. For example, in the intracluster medium, the dynamic length scale ∼ 10 kpc, the ion Larmor scale ∼ 10 -12 kpc and the electron Larmor scale ∼ 10 -14 kpc, making PIC codes infeasible for certain physical problems.  \nAnother numerical complication is modelling turbulence which is prevalent in most astrophysics systems. Turbulence is also present in weakly collisional astrophysical environments like the hot intracluster medium (ICM) of galaxy clusters and solar wind. Turbulence in such environments may be driven by merger events, feedback from active galactic nuclei, or by in-falling galaxies stirring the ICM (Simionescu et al. 2019). Observational signatures of turbulence in such environments have been studied by Hitomi Collaboration et al. (2016); Gatuzz et al. (2022a,b, 2023). Regions of the ICM can also host subsonic or supersonic turbulence of rotational or irrotational nature (Zinger et al. 2018). This is also true of collisionless astrophysical plasma like the solar wind (Howes et al. 2014; Bruno & Carbone 2013; Klein & Howes 2015; Verscharen et al. 2019).  \nWe adopt the middle ground approach of solving the Hybridkinetic equations using so-called HybridPIC in our new code AHKASH 1 , which builds on the architecture of the FLASH code (Fryxell et al. 2000; Dubey et al. 2008a,b). In Hybrid PIC, positively charged ions are modelled as particles, which sample the moments of the ion distribution function, while the electrons are modelled as a mass-less fluid 2 (Horowitz et al. 1989; Lipatov 2002; Bagdonat & Motschmann 2002; Gargaté et al. 2007; Müller et al. 2011; Kunz et al. 2014b; Le et al. 2023). We note that we solve the Hybrid-kinetic equations in the non-relativistic limit (the ratio of typical speeds to the speed of light in the hot ICM ≲ 10 -3 ). With this approach, we capture the physics at the ion scales, like the Larmor gyration of ions and the damping of ion acoustic waves, without having to encapsulate electron scale physics, such as Debye oscillations. This also significantly reduces the computational cost of performing these simulations relative to the full PIC approach, while still being more computationally expensive than MHD simulations, but allowing us to model the physics of weakly collisional plasmas.  \nThe viscous heating of plasma due to turbulence, and the thermodynamics of weakly collisional plasma are different from collisional MHD plasma (Howes 2010; Kunz et al. 2011; Squire et al. 2023a). It is therefore important to study plasma turbulence and heating of collisionless plasmas like the solar wind, by using PIC methods. For example, recent studies use Hybrid PIC to study the heating of ions and electrons in the solar wind (Arzamasskiy et al. 2019; Squire et al. 2023b). The Hybrid PIC approach also enables the study of magnetic fields in a turbulent, magnetized, weakly collisional plasma (St-Onge & Kunz 2018; Achikanath Chirakkara et al. 2024). While Rincon et al. (2016) solve the Vlasov equation in their numerical experiments, St-Onge & Kunz (2018) and Achikanath Chirakkara et al. (2024) use Hybrid PIC simulations. We note that the numerical experiments described in Achikanath Chirakkara et al. (2024) have been performed using AHKASH .  \nTurbulence affects the thermodynamics of the plasma and controls the evolution of magnetic fields. Therefore, to properly understand the evolution of weakly collisional astrophysical systems, it is important to numerically model turbulence. Several kinetic codes have been used to study turbulence in weakly collisional or collisionless plasma  \n1 AHKASH is an acronym for Astrophysical Hybrid-Kinetic simulations with the flASH code. AHKASH (pronounced as ¯ak¯aś) means 'the sky' in many Indian languages such as Hindi, Malayalam, and Gujarati.  \n(Spitkovsky 2005; Kunz et al. 2014b). Turbulence dissipates energy and heats the plasma which makes it difficult to maintain isothermal conditions, especially for controlled numerical experiments. Here we aim to accurately model the effects of turbulence in collisionless plasma with our code using an Ornstein-Uhlenbeck process (Federrath et al. 2010, 2022), especially to maintain the isothermal nature of plasma via the novel cooling in the presence of plasma heating due to turbulence.  \nWeuse the second-order accurate predictor-predictor-corrector integration scheme to update the particles and electromagnetic fields in our code, similar to Kunz et al. (2014b). We use the Boris integrator (Boris 1970) and the constrained transport method to update the magnetic field while keeping it divergence-free (Yee 1966). We also implemented a 𝛿 𝑓 method to enable the study of instabilities in collisonless plasma. We study the effect of using different interpolation kernels, performing smoothing operations and varying the number of particles per grid cell in our numerical simulations. Finally, we have implemented a new plasma cooling method, which can be used to maintain a steady sonic Mach number in driven fully-developed turbulence in weakly collisional plasmas, even when such plasmas experience strong local compression and/or rarefaction. This method is important for studying turbulence in collisionless plasma, particularly in the supersonic regime.  \nThe rest of this paper is organised as follows. We describe the Hybrid-kinetic equations and the related numerical methods we use in Sec. 2. We test the accuracy of our new Hybrid PIC code with various standard test problems in Sec. 3. We discuss the interpolation tests in Sec. 4 and describe the ion cooling method we developed alongside numerical tests in Sec. 5. We discuss the 'hybrid precision' method, performance of the code and present scaling tests in Sec. 6. Finally, we summarize our study in Sec. 7.", '2 METHODS': 'In this section, we present the equations and methods used in AHKASH . Sec. 2.1 introduces the equations that govern the motion and coupling of positively-charged ions (particles) with an electron fluid in Hybrid-kinetic plasmas. Sec. 2.2 describes the operations that the code performs in each time step. To track ion trajectories, we use the Boris integrator, described in Sec. 2.3. The particle-grid interpolation and smoothing operations are explained in Sec. 2.4. The 𝛿 𝑓 method is introduced in Sec. 2.5, and the magnetic field update is described in Sec. 2.6. We discuss the predictor-predictor-corrector integration scheme in Sec. 2.7, with the corrections to the interpolated electromagnetic fields explained in Sec. 2.8, and higher-order numerical hyper-resistivity discussed in Sec. 2.9. Finally, we describe the turbulent driving in Sec. 2.10, and the novel cooling method, introduced to overcome plasma heating, in Sec. 2.11. Time-step constraints are discussed in Sec. 2.12.', '2.1 Hybrid-kinetic equations': "The evolution of a weakly collisional plasma is governed by the Vlasov equation, which ensures the conservation of phase-space den-  \nsity,  \nd 𝑓 i ( 𝑡, r , v ) d 𝑡 = 0 , 𝜕 𝑓 i 𝜕𝑡 + 𝜕 r 𝜕𝑡 · 𝜕 𝑓 i 𝜕 r + 𝜕 v 𝜕𝑡 · 𝜕 𝑓 i 𝜕 v = 0 , 𝜕 𝑓 i 𝜕𝑡 + v · 𝜕 𝑓 i 𝜕 r + 𝑞 i 𝑚 i ( E + v × B ) · 𝜕 𝑓 i 𝜕 v = 0 , (1)  \nwhich describes the evolution of the ion distribution function, 𝑓 i ( 𝑡, r , v ) , where 𝑡 is time, r is the position vector, and v is the velocity vector 3 . The electromagnetic force on the ions is described by the Lorentz force, ( 𝑞 i / 𝑚 i )( E + v × B ) , where 𝑞 i and 𝑚 i are the charge and mass of the ions, respectively. E and B are the electric and magnetic fields experienced by the ions. We solve Eq. 1 using the PIC approach, for which we can write the equations of motion for the ions,  \nd r 𝑝 d 𝑡 = v 𝑝 , (2)  \n𝑝  \nd  \nv  \n𝑡  \nd  \n=  \n𝑞  \ni  \n𝑚  \ni  \nE  \n+  \nv  \n𝑝  \n×  \nB  \n,  \n(3)  \nwhere r 𝑝 and v 𝑝 are the position and velocity of the 𝑝 th particle, respectively, and 𝑝 = 1 , 2 , . . . , 𝑃 , where 𝑃 is the total number of particles. We note that this model can be easily expanded for multiple positively-charged ion species, however, here we only consider protons.  \nThe evolution of the electrons is described by the Vlasov-Landau equation (similar to Eq. 1, but now with a collision operator on the right-hand side),  \n𝜕 𝑓 e 𝜕𝑡 + v · 𝜕 𝑓 e 𝜕 r + 𝑞 e 𝑚 e ( E + v × B ) · 𝜕 𝑓 e 𝜕 v = 𝐶 [ 𝑓 e ] , (4)  \nwhere 𝑓 e ( 𝑡, r , v ) , 𝑞 e, and 𝑚 e, are the electron distribution function, charge, and mass of the electrons, respectively, and 𝐶 [ 𝑓 e ] is the electron collision operator.  \nPerforming an expansion of the electron distribution function in powers of the electron-to-ion mass ratio ( 𝑚 e / 𝑚 i ) 1 / 2 , we can simplify the Vlasov-Landau equation (Appendix A1.1 of Rosin et al. (2011)). Up to the lowest orders in this expansion, one finds that the zeroth-order distribution function is Maxwellian as shown in Appendix A1.3 of Rosin et al. (2011). Multiplying the expression obtained from Eq. 4 up to order ( 𝑚 e / 𝑚 i ) 1 / 2 by 𝑚 e 𝑣 , integrating over velocity space and using that the zeroth-order electron distribution function is Maxwellian, as detailed in Appendix A1.5 of Rosin et al. (2011), we obtain the following the Ohm's law for the electric field with resistivity (index 'nores'),  \nE nores = -( J e × B ) 𝜌 e + ∇ 𝑝 e 𝜌 e , (5)  \nwhere J e, 𝜌 e, and 𝑝 e, are the electron current, electron charge density, and electron pressure, respectively. The electron collision operator is modelled using Ohmic resistivity, which gives us  \nE = -( J e × B ) 𝜌 e + 𝑘 B 𝑇 e ∇ 𝑛 e 𝜌 e + 𝜂𝜇 0 J , (6)  \nwhere J is the total current and 𝜂 is the magnetic diffusivity. This additional term acts as a sink for magnetic energy. 𝜇 0 is the magnetic permeability constant. We have also used the ideal gas equation of state in Eq. 6 with a constant temperature for the electrons,  \n3 Note that 𝜕𝑓 𝜕 r = GLYPH<16> 𝜕𝑓 𝜕𝑥 , 𝜕𝑓 𝜕𝑦 , 𝜕𝑓 𝜕𝑧 GLYPH<17> , and 𝜕𝑓 𝜕 v = GLYPH<16> 𝜕𝑓 𝜕𝑣𝑥 , 𝜕𝑓 𝜕𝑣𝑦 , 𝜕𝑓 𝜕𝑣𝑧 GLYPH<17> .  \nGLYPH<0>  \nGLYPH<1>  \n𝑝 e = 𝑛 e 𝑘 B 𝑇 e, where 𝑛 e, 𝑇 e, and 𝑘 B , are the electron number density, constant electron temperature, and the Boltzmann constant, respectively. The total current is the sum of the ion current ( J i ) and the electron current ( J e),  \nJ = J i + J e . (7)  \nThe charge density of the ions, 𝜌 i , can be obtained from the zeroth moment of the distribution function and can be written as  \n𝜌 i = 𝑞 i ∫ 𝑓 i 𝑑 3 v , (8)  \nand assuming quasi-neutrality of the plasma, we have  \n𝜌 i + 𝜌 e = 0 , (9)  \nwhich implies that the charge density of the ions and electrons is the same, 𝜌 i = -𝜌 e.  \nThe first moment of the distribution function corresponds to the mean (or bulk) velocity of the ions, which defines the ion current,  \nJ i = 𝑞 i ∫ v 𝑓 i 𝑑 3 v . (10)  \nUsing this and 𝜌 e = 𝑞 e 𝑛 e in Eq. 6, the Ohm's law can be re-written as  \nE = ( J -J i ) × B 𝜌 i -𝑘 B 𝑇 e 𝑞 e ∇ 𝜌 i 𝜌 i + 𝜇 0 𝜂 J . (11)  \nFinally, the evolution of the positively-charged ions and electrons and the Ohm's law (Eq. 11) are coupled to the following Maxwell's equations,  \nJ = ∇ × B 𝜇 0 , (12)  \nwith the magnetic field calculated from Faraday's law,  \n𝜕 B 𝜕𝑡 = -∇ × E , (13)  \nand the solenoidal condition for the magnetic field,  \n∇ · B = 0 . (14)", '2.2 Hybrid PIC algorithm': 'We solve the Hybrid-kinetic equations described in Sec. 2.1 using the PIC method. Particles are used to sample the moments of the ion distribution function, 𝑓 i . In this section, we go through a simplified description of what the code does in one computational time step, Δ 𝑡 . A schematic diagram of this simplified algorithm is shown in Fig. 1 and the computational loop is described below.  \n- (i) Particle evolution: In the presence of electric and magnetic fields, the ions are accelerated by the Lorentz force, 𝑞 i / 𝑚 i ( E nores + v × B ) , where E nores is the electric field experienced by the ions, excluding resistive effects. The resistive term is excluded here to ensure the conservation of momentum of the overall system (see Sec. 4.5.2 of Lipatov (2002) and Bagdonat (2004)). The particle velocities and positions are updated based on this acceleration and velocity, respectively. More details on the particle evolution are provided in Sec. 2.3.\n- (ii) Depositing moments of the distribution function onto the computational grid: The ion charge density and current are obtained from the zeroth and first moments of the distribution function, using Eq. 8 and Eq. 10. These moments are sampled by the particles and deposited onto the grid, using an interpolation method (see Sec. 2.4 and Sec. 2.5).', '(i) Update particle positions and velocities': "Figure 1. Simplified schematic of the Hybrid particle-in-cell (Hybrid PIC) algorithm for one time step Δ 𝑡 (see Sec. 2.2). (i) The particle positions and velocities are updated based on the velocity and acceleration from the Lorentz force, respectively (see Sec. 2.1 and Sec. 2.3). (ii) Following this, the updated charge density, 𝜌 i, and the ion current, J i, are deposited onto the computational grid using an interpolation operation (see Sec. 2.4 and Sec. 2.5). (iii) The charge density and ion current are source terms in the Ohm's law, Eq. 34, which is used to compute the electric field, E , on the grid. (iv) The magnetic field, B , is then calculated using Faraday's law, Eq. 13, on the grid (see Sec. 2.6). Finally, the updated electric and magnetic fields are interpolated from the grid to the particles for the next particle evolution step.  \n<!-- image -->  \n(iii) Computing the electric field (Ohm's law): The charge density and current on the grid are used in the Ohm's law, Eq. 11. The total current is calculated at this stage using Ampère's law, Eq. 12 (see Sec. 2.1).  \n(iv) Computing the magnetic field (Faraday's law): The magnetic field update on the grid is described in more detail in Sec. 2.6. The E and B fields are now interpolated from the grid to the particles for the next cycle.  \nOur code uses the predictor-predictor-corrector integration algorithm (Kunz et al. 2014b), which is more sophisticated and involved than the simplified description above, and is second-order accurate, capable of stably propagating Alfvén and whistler waves. We provide details on the full scheme in Sec. 2.7. The time step, Δ 𝑡 , is computed as the minimum of several particle and wave time steps, as discussed in detail in Sec. 2.12.", '2.3 Particle integrator': "To update the particle positions and velocities, we solve the equations of motion, Eq. 2 and 3, using the Boris integrator. The Boris integrator is commonly used in PIC codes for its ability to accurately evolve particle dynamics for a large number of time steps (Boris 1970; Qin et al. 2013). The velocity update in one time step, Δ 𝑡 , can be denoted as v 𝑡 → v 𝑡 + Δ 𝑡 , where v 𝑡 is the velocity at time 𝑡 and v 𝑡 + Δ 𝑡 is the velocity at time 𝑡 + Δ 𝑡 . Then, v 𝑡 + Δ 𝑡 is calculated via  \nv -= v 𝑡 + 𝑎 e E 𝑡 + Δ 𝑡 / 2 nores (15)  \nv + = v -+ 𝑎𝑏 v ' × e B 𝑡 + Δ 𝑡 / 2 (17)  \nv ' = v -+ 𝑎 v -× e B 𝑡 + Δ 𝑡 / 2 (16)  \nv 𝑡 + Δ 𝑡 = v + + 𝑎 e E 𝑡 + Δ 𝑡 / 2 nores , (18)  \nwhere v -, v ' and v + are intermediate time-step solutions for the velocity. The coefficients 𝑎 and 𝑏 are calculated each time-step and defined as 𝑎 = 𝑞 i /( 2 𝑚 i ) Δ 𝑡 and 𝑏 = 2 /( 1 + 𝑎 2 e B 2 ) . e E 𝑡 + Δ 𝑡 / 2 nores and e B 𝑡 + Δ 𝑡 / 2 are the electric and magnetic fields interpolated to the particles at the intermediate time-step, 𝑡 + Δ 𝑡 / 2 (how these are obtained is described in Sec. 2.7). As shown above, the scheme splits the electric field evolution into two steps. Eq. 15 describes the first electric field acceleration step at the beginning of the scheme. The magnetic field evolution is done in between the electric field updates, in the form of two rotations as shown by Eq. 16 and 17 (see fig. 4 in Boris 1970). The final electric field advancement is done in Eq. 18, and the velocity update is complete.  \nThe update of the particle positions, r 𝑡 → r 𝑡 + Δ 𝑡 , is carried out as  \nr 𝑡 + Δ 𝑡 = r 𝑡 + Δ 𝑡 2 ( v 𝑡 + v 𝑡 + Δ 𝑡 ) . (19)", '2.4.1 Interpolation': "The PIC method requires particle-to-grid and grid-to-particle interpolations (see Fig. 1). The grid is used to evolve the electromagnetic fields, while the particles represent the ions, which exhibit Lagrangian motion governed by the electric and magnetic fields through the equations of motion, Eq. 2 and 3.  \nIn general, the interpolation operation provides us with the interpolated scalar quantity Q at position r , based on a quantity 𝑄 at positions r 𝑙 , via  \nQ( r ) = 𝐿 ∑︁ 𝑙 = 1 𝑄 ( r 𝑙 ) 𝑊 ( r -r 𝑙 ) , (20)  \nwhere 𝑊 ( r -r 𝑙 ) is the weight (interpolation kernel) function, and  \n𝑙 = 1 . . . 𝐿 denotes the index of all 𝐿 particles (cells) within the kernel, for interpolations from the particles (grid cells) to the grid (particles). We describe the nearest-grid-point (NGP), cloud-in-cell (CIC), and the triangular-shaped-cloud (TSC) weight functions in Appendix A. For vector quantities, Eq. (20) is performed componentwise.  \nFor grid-to-particle (G → P) interpolations, r 𝑙 = r 𝑔 and r = r 𝑝 in Eq. 20. In this case, 𝑄 ( r 𝑙 ) = 𝑄 ( r 𝑔 ) is the grid quantity to be interpolated to the particles. For example, the electromagnetic fields are interpolated from the centre of the computational grid cells to each particle's position using the same weight function used to perform the P → G interpolations. This is described by Fig. 2b. We note that the default interpolation kernel we use is the CIC weight function.  \nFor particle-to-grid (P → G) interpolations, r 𝑙 = r 𝑝 are the particle positions and r = r 𝑔 is the centre of each computational grid cell in Eq. 20. This is shown by Fig. 2a. The quantity 𝑄 ( r 𝑙 ) = 𝑄 ( r 𝑝 ) may be further expressed as a product, 𝑄 ( r 𝑝 ) = 𝑞 𝑝 v 𝑚 𝑝 , where 𝑞 𝑝 and v 𝑝 are the charge and the velocity of the 𝑝 th particle, respectively, and 𝑚 = 0 , 1 , 2 , ... corresponds to the 𝑚 th moment of the ion distribution function. The most common application of P → G interpolations is to compute the charge density, which is the zeroth moment ( 𝑚 = 0), 𝜌 i = Q( r )/ Δ 𝑉 , where Δ 𝑉 = Δ 𝑥 Δ 𝑦 Δ 𝑧 , based on the charge of each particle, 𝑄 ( r 𝑝 ) = 𝑞 𝑝 . Likewise, the first moment ( 𝑚 = 1) represents the ion current, J i = Q( r )/ Δ 𝑉 , by setting 𝑄 ( r 𝑝 ) = 𝑞 𝑝 v 𝑝 .  \nThe influence of using different kernels (NGP, CIC, TSC) is quantified in Sec. 4.", '2.4.2 Post-interpolation smoothing to reduce particle noise': 'As described above, we sample the moments of the ion distribution function using averages over a population of particles. However, since one can only use a finite number of particles, the distribution function can only be sampled to a certain accuracy, and therefore finite particle counts always introduce noise in any of the particleto-grid interpolated quantities. This can be particularly problematic for the charge density, as the particle noise may introduce spurious electric fields, due to the fact that the thermo-electric term in Eq. (11) generates electric fields based on gradients in the electron pressure (which is proportional to the charge density). Thus, noise in the charge density may lead to large spurious gradients. To address this issue, we perform a smoothing operation on the grid for the ion charge density and current, every time these moments are deposited from the particles to the grid.  \n𝐷 ( 𝑖, 𝑗, 𝑘 ) → 𝐷 ( 𝑖 -1 , 𝑗 , 𝑘 ) 4 + 𝐷 ( 𝑖, 𝑗, 𝑘 ) 2 + 𝐷 ( 𝑖 + 1 , 𝑗 , 𝑘 ) 4 , (21)  \nThe smoothing (filtering) operation is defined such that a P → G interpolated quantity, 𝐷 , at cell index ( 𝑖, 𝑗, 𝑘 ) is transformed as  \nin the 𝑖 th spatial direction. In case of a 3D simulation, this operation is performed in the 𝑖 th , 𝑗 th and 𝑘 th direction, consecutively. This means that during a single smoothing operation, cell ( 𝑖, 𝑗, 𝑘 ) passes 1/2 of its value to its left and right neighbours, and receives 1/4 of the value of its left and right neighbour, respectively, therefore conserving the overall sum.  \nThe smoothing operation (Eq. 21) can be repeated 𝑁 smooth times to achieve a smoother and smoother representation of the interpolated quantity, 𝐷 ( 𝑖, 𝑗, 𝑘 ) . We investigate the appropriateness of the smoothing operation in Sec. 4 below. We find that using 2 smoothing passes ( 𝑁 smooth = 2) avoids excessive smoothing, retaining smallscale detail, and at the same time provides an appropriate reduction of particle noise, at least in applications of turbulent dynamo amplification. We note that the details of required smoothing passes  \nmay depend on the specific problem and on the affordable number of particles per cell (with tests provided in Sec. 4).', '2.5 Studying perturbations with Hybrid particle-in-cell simulations': 'The 𝛿 𝑓 method can be used to study perturbations and instabilities in collisionless plasmas (Kunz et al. 2014a,b). In this method, the particles are used to sample the difference between the full distribution function, 𝑓 , and the zeroth-order equilibrium distribution function, 𝑓 0 ,  \n𝛿 𝑓 = 𝑓 -𝑓 0 . (22)  \nFirst, the ion charge density is deposited onto the grid as  \n𝜌 i ( r 𝑔 ) = 𝜌 0 ( r 𝑔 ) + ˜ 𝜌 i ( r 𝑔 ) , (23)  \nwhere 𝜌 0 ( r 𝑔 ) is the equilibrium charge density, and ˜ 𝜌 i ( r 𝑔 ) is obtained by a P → G interpolation with 𝑄 ( r 𝑝 ) = 𝑞 𝑝 𝑤 𝑝 in Eq. 20, where 𝑤 𝑝 is the factor the particle moments are weighted with in the 𝛿 𝑓 method (Parker & Lee 1993; Denton & Kotschenreuther 1995). The ion current is deposited as  \nJ i ( r 𝑔 ) = 𝜌 0 ( r 𝑔 ) 𝑣 0 ( r 𝑔 ) + ˜ J i ( r 𝑔 ) , (24)  \nwhere 𝑣 0 ( r 𝑔 ) is the equilibrium mean velocity, obtained from 𝑓 0 . ˜ J i ( r 𝑔 ) is obtained by a P → Ginterpolation with 𝑄 ( r 𝑝 ) = 𝑞 𝑝 v 𝑝 𝑤 𝑝 in Eq. 20. A common choice for the equilibrium distribution function is the Maxwell distribution function, for which 𝜌 0 ( r 𝑔 ) = 𝜌 0 = 𝑞 i 𝑛 0 , where 𝑛 0 is the mean number density, and 𝑣 0 ( r 𝑔 ) = 0. The 𝛿 𝑓 method is useful to study the plasma when 𝛿 𝑓 ≪ 𝑓 . However, as the perturbations grow and 𝛿 𝑓 ∼ 𝑓 , the 𝛿 𝑓 method is no longer suitable to model the plasma, which is usually the case when turbulence is present.', '2.6 Constrained transport and solenoidality of magnetic fields': "To update the magnetic fields we use the constrained transport (CT) method on the three-dimensional computational grid to solve Faraday's law, Eq. 13 (Yee 1966). Fig. 3 shows a schematic of a single grid cell with the location of the ion charge density, ion current, and electric and magnetic fields. The charge and current density are sampled from the particles and deposited at the centre of each grid cell as 𝜌 i ( 𝑖, 𝑗, 𝑘 ) and J i ( 𝑖, 𝑗, 𝑘 ) , respectively. We note that ( 𝑖, 𝑗, 𝑘 ) represents the coordinate of the cell centre denoted as r 𝑔 , which is described in Sec. 2.4.1.  \nThe magnetic fields are calculated on the face-center of the grid cell and represent the face-averaged value of the magnetic fields. The 𝑥 -component of the magnetic field is calculated on the face of the grid cell parallel to the 𝑦 -𝑧 plane, 𝐵 𝑥 ( 𝑖 -1 / 2 , 𝑗 , 𝑘 ) . The magnetic field's 𝑦 and 𝑧 -components are defined similarly to the 𝑥 -component. These points are depicted by the grey circles in Fig. 3.  \nIn the following, we discretize the electromagnetic fields on the grid, where 𝑖 → 𝑖 + 1 represents 𝑥 → 𝑥 + Δ 𝑥 , where Δ 𝑥 is the grid cell length along the 𝑥 direction, and analogously for 𝑗 and 𝑘 in the 𝑦 and 𝑧 direction, respectively. The electric fields, E , required to update the magnetic fields, B , are calculated on the edge centres of the grid cell using the Ohm's law, Eq. 34. The 𝑥 -component of the electric field is calculated on the edge centre parallel to the 𝑥 -axis, 𝐸 𝑥 ( 𝑖, 𝑗 -1 / 2 , 𝑘 -1 / 2 ) . The 𝑦 and 𝑧 -components of the electric field are defined in analogy to the 𝑥 -component. These are shown by the blue squares on the edge centre of the grid cell in Fig. 3.  \nNext, we discretize Faraday's law and update the magnetic fields in  \n<!-- image -->  \n/uni0394 /u1D465 (a) Particle-to-grid (P → G) interpolation  \n/uni0394  \n/u1D465  \n<!-- image -->  \n(b) Grid-to-particle (G → P) interpolation  \n/uni0394 /u1D466 x x This paper has been typeset from a T E X/L A T E X file prepared by the author. Figure 2. Schematic diagram describing (a) particle-to-grid (P → G) and (b) grid-to-particle (G → P) interpolations with the cloud-in-cell weight function in two-dimensions. The shaded region shows the cloud-in-cell interpolation kernel ( 𝑊 ), which has the same dimensions as the grid cell, centred on the particle position, r 𝑝 . The cell centres, r 𝑔 , are depicted by the cross symbols. The arrows show the direction of the interpolation (a) from the particles to the grid cell centres and (b) from the grid cells to the particle (see Sec. 2.4.1 for further details).  \n/uni0394 /u1D466 /u1D44A ( r /u1D45D -r /u1D454 ) time. In one time step, i.e., 𝑡 → 𝑡 + Δ 𝑡 , the evolution of the magnetic field components is given by  \nr  \n/u1D45D  \n/uni0394 /u1D465 r /u1D454 x x 𝐵 𝑡 + Δ 𝑡 𝑥 ( 𝑖 -1 / 2 , 𝑗 , 𝑘 ) = 𝐵 𝑡 𝑥 ( 𝑖 -1 / 2 , 𝑗 , 𝑘 ) + Δ 𝑡 Δ 𝑧 GLYPH<2> e 𝐸 𝑡 + Δ 𝑡 / 2 𝑦 ( 𝑖 -1 / 2 , 𝑗 , 𝑘 + 1 / 2 ) -e 𝐸 𝑡 + Δ 𝑡 / 2 𝑦 ( 𝑖 -1 / 2 , 𝑗 , 𝑘 -1 / 2 ) GLYPH<3> -Δ 𝑡 Δ 𝑦 GLYPH<2> e 𝐸 𝑡 + Δ 𝑡 / 2 𝑧 ( 𝑖 -1 / 2 , 𝑗 + 1 / 2 , 𝑘 ) -e 𝐸 𝑡 + Δ 𝑡 / 2 𝑧 ( 𝑖 -1 / 2 , 𝑗 -1 / 2 , 𝑘 ) GLYPH<3> , (25)  \npaper has been typeset from a T  \nE  \nX/L  \nA  \nT  \nE  \nlated at the cell centre is thus  \n∇ · B 𝑡 + Δ 𝑡 ( 𝑖, 𝑗, 𝑘 ) = 1 Δ 𝑥 GLYPH<2> 𝐵 𝑡 + Δ 𝑡 𝑥 ( 𝑖 + 1 / 2 , 𝑗 , 𝑘 ) -𝐵 𝑡 + Δ 𝑡 𝑥 ( 𝑖 -1 / 2 , 𝑗 , 𝑘 ) GLYPH<3> + 1 Δ 𝑦 GLYPH<2> 𝐵 𝑡 + Δ 𝑡 𝑦 ( 𝑖, 𝑗 + 1 / 2 , 𝑘 ) -𝐵 𝑡 + Δ 𝑡 𝑦 ( 𝑖, 𝑗 -1 / 2 , 𝑘 ) GLYPH<3> + 1 Δ 𝑧 GLYPH<2> 𝐵 𝑡 + Δ 𝑡 𝑧 ( 𝑖, 𝑗, 𝑘 + 1 / 2 ) -𝐵 𝑡 + Δ 𝑡 𝑧 ( 𝑖, 𝑗, 𝑘 -1 / 2 ) GLYPH<3> . (28)  \nX file prepared by the author.  \n/uni2605 𝐵 𝑡 + Δ 𝑡 𝑦 ( 𝑖, 𝑗 -1 / 2 , 𝑘 ) = 𝐵 𝑡 𝑦 ( 𝑖, 𝑗 -1 / 2 , 𝑘 ) + Δ 𝑡 Δ 𝑥 GLYPH<2> e 𝐸 𝑡 + Δ 𝑡 / 2 𝑧 ( 𝑖 + 1 / 2 , 𝑗 -1 / 2 , 𝑘 ) -e 𝐸 𝑡 + Δ 𝑡 / 2 𝑧 ( 𝑖 -1 / 2 , 𝑗 -1 / 2 , 𝑘 ) GLYPH<3> -Δ 𝑡 Δ 𝑧 GLYPH<2> e 𝐸 𝑡 + Δ 𝑡 / 2 𝑥 ( 𝑖, 𝑗 -1 / 2 , 𝑘 + 1 / 2 ) -e 𝐸 𝑡 + Δ 𝑡 / 2 𝑥 ( 𝑖, 𝑗 -1 / 2 , 𝑘 -1 / 2 ) GLYPH<3> , (26)  \n©  \n𝐵 𝑡 + Δ 𝑡 𝑧 ( 𝑖, 𝑗, 𝑘 -1 / 2 ) = 𝐵 𝑡 𝑧 ( 𝑖, 𝑗, 𝑘 -1 / 2 ) + Δ 𝑡 Δ 𝑦 GLYPH<2> e 𝐸 𝑡 + Δ 𝑡 / 2 𝑥 ( 𝑖, 𝑗 + 1 / 2 , 𝑘 -1 / 2 ) -e 𝐸 𝑡 + Δ 𝑡 / 2 𝑥 ( 𝑖, 𝑗 -1 / 2 , 𝑘 -1 / 2 ) GLYPH<3> -Δ 𝑡 Δ 𝑥 GLYPH<2> e 𝐸 𝑡 + Δ 𝑡 / 2 𝑦 ( 𝑖 + 1 / 2 , 𝑗 , 𝑘 -1 / 2 ) -e 𝐸 𝑡 + Δ 𝑡 / 2 𝑦 ( 𝑖 -1 / 2 , 𝑗 , 𝑘 -1 / 2 ) GLYPH<3> . (27)  \nE-mail: publications@ras.ac.uk (KTS) We also require the cell-centred values of the magnetic fields to interpolate these fields to the particle positions. The cell-centered values of the magnetic field at 𝑡 + Δ 𝑡 are obtained by averaging the face-centred values.  \n2024 The Authors The divergence of the magnetic field from the CT method calcu-  \nDue to the staggered positions of the electric and magnetic fields in the construction of the stencil, we find ∇ · B 𝑡 + Δ 𝑡 ( 𝑖, 𝑗, 𝑘 ) = 0 when Eqs. 25-27 are plugged into Eq. 28. The construction of the CT method ensures that the analytical expression, ∇ · (∇ × C ) = 0, is satisfied for any vector C , and therefore it can be used to ensure the solenoidality of the magnetic field, to machine precision.  \nE-mail: publications@ras.ac.uk (KTS) 2024 The Authors For comparison with the CT method, we also describe the cellcentred finite difference method to update magnetic fields in Appendix B, which we show also ensures that the magnetic fields are divergence-free. Similar to the CT method, the cell-centred finite difference method ensures the solenoidality of magnetic fields (although on a larger stencil, which has exactly twice the linear extent of the CT method) using the construction of the numerical stencil. Our code includes both the CT and the cell-centred finite difference method to update magnetic fields, but we use the CT method as the default.", '2.7 Predictor-predictor-corrector integration scheme': "Hybrid PIC codes use different integration schemes to advance the particle positions, velocities, and electromagnetic fields in time. Examples of these are the predictor-corrector method and the cyclic leapfrog (CL) with the current-advance method (CAM-CL) (Holmstrom 2009; Holmström 2011; Winske et al. 2022). Here we use the predictor-predictor-corrector algorithm (Kunz et al. 2014a), which requires two predictions for the electric and magnetic fields at an intermediate time-step, 𝑡 + Δ 𝑡 / 2. The predictor-predictor-corrector method is second-order accurate and can stably propagate Alfvén  \nx  \nx  \nr  \n/u1D454  \n/uni0394  \n/u1D465  \nr  \n/uni0394  \n/u1D465  \nr  \n/u1D45D  \n/u1D45D  \n/u1D446  \n/u1D446  \n/uni0394  \n/u1D465  \nx  \nr  \n(  \n/u1D454  \nx  \n/uni0394  \n/u1D465  \nx  \nr  \n/u1D45D  \nx  \nr  \n-  \nr  \n/u1D454  \n-  \n/u1D454  \nREFERENCES Figure 3. Staggered construction of the electric and magnetic fields for the constrained transport method to solve Faraday's law, Eq. 13. The electric fields are depicted by the blue squares on the edge centres and the magnetic fields are shown by the grey circles on the face centres of the grid cell. The charge density and ion current are always deposited on the grid-cell centre, shown by the yellow circle.  \n<!-- image -->  \nThis paper has been typeset from a T  \nAs described in Sec. 2.3, to update the particle velocity, v 𝑡 → v 𝑡 + Δ 𝑡 , we require the electromagnetic fields e E 𝑡 + Δ 𝑡 / 2 and e B 𝑡 + Δ 𝑡 / 2 at the intermediate time-step, 𝑡 + Δ 𝑡 / 2. Furthermore, the constrained transport method also requires the electric field, e E 𝑡 + Δ 𝑡 / 2 , to advance the magnetic fields, B 𝑡 → B 𝑡 + Δ 𝑡 , as discussed in Sec. 2.6. Therefore, we need estimates of e E 𝑡 + Δ 𝑡 / 2 and e B 𝑡 + Δ 𝑡 / 2 .  \nand whistler waves, as demonstrated in Sec. 3.2. A complete description of the predictor-predictor-corrector integration scheme is provided in Kunz et al. (2014b).  \n(KTS)  \ne E 𝑡 + Δ 𝑡 / 2 = 1 2 ( E 𝑡 + E 𝑡 + Δ 𝑡 pred , 1 ) , (29)  \nThe electric field is calculated from the Ohm's law, Eq. 34, and is a function of the charge density, ion current and magnetic field, E ( 𝜌 i , J i , B ) . In the first predictor step of the predictor-predictorcorrector algorithm, B 𝑡 + Δ 𝑡 pred , 1 is calculated using E 𝑡 . The first prediction for the electric field, E 𝑡 + Δ 𝑡 pred , 1 , is calculated using 𝜌 𝑡 i , J 𝑡 i and B 𝑡 + Δ 𝑡 pred , 1 . Once E 𝑡 + Δ 𝑡 pred , 1 and B 𝑡 + Δ 𝑡 pred , 1 are estimated, the electromagnetic fields at 𝑡 + Δ 𝑡 / 2 are calculated as  \ne B 𝑡 + Δ 𝑡 / 2 = 1 2 ( B 𝑡 + B 𝑡 + Δ 𝑡 pred , 1 ) . (30)  \nThesecondprediction for the magnetic field, B 𝑡 + Δ 𝑡 pred , 2 , is calculated using e E 𝑡 + Δ 𝑡 / 2 . We advance the particles from 𝑡 → 𝑡 + Δ 𝑡 using e E 𝑡 + Δ 𝑡 / 2 and e B 𝑡 + Δ 𝑡 / 2 to obtain 𝜌 𝑡 + Δ 𝑡 i and J 𝑡 + Δ 𝑡 i , which we use along with B 𝑡 + Δ 𝑡 pred , 2 to predict E 𝑡 + Δ 𝑡 pred , 2 . Advancing the particles is computationally the most expensive step in the algorithm, but using an additional particle-advancement step helps us better estimate the electromagnetic fields directly from the motion of the particles. The quantities e E 𝑡 + Δ 𝑡 / 2 and e B 𝑡 + Δ 𝑡 / 2 are then re-calculated using E 𝑡 + Δ 𝑡 pred , 2 and B 𝑡 + Δ 𝑡 pred , 2 ,  \nX/L  \nA T E similar to Eq. 29-30,  \nX file prepared by the author.  \ne E 𝑡 + Δ 𝑡 / 2 = 1 2 ( E 𝑡 + E 𝑡 + Δ 𝑡 pred , 2 ) , (31)  \ne B 𝑡 + Δ 𝑡 / 2 = 1 2 ( B 𝑡 + B 𝑡 + Δ 𝑡 pred , 2 ) . (32)  \nFinally, B 𝑡 + Δ 𝑡 is calculated using e E 𝑡 + Δ 𝑡 / 2 and the particle positions and velocities are updated using e E 𝑡 + Δ 𝑡 / 2 and e B 𝑡 + Δ 𝑡 / 2 .", '2.8 Corrections for interpolated electromagnetic fields': 'In the absence of thermo-electric and resistive effects, the electric and magnetic fields are perpendicular to each other during their calculation on the computational grid. However, when these fields are interpolated to the particles, this may no longer be true due to interpolation errors. To address this issue, we perform the following correction to the interpolated electric field at the particle positions (Lehe et al. 2009; Kunz et al. 2014b; Achikanath Chirakkara et al. 2024),  \nE ∗ nores = E nores + [( E nores · B ) G → P -E nores · B ] B | B | 2 , (33)  \nwhere all the terms in the above equation are fields interpolated to the particle positions. ( E nores · B ) G → P is the dot product of E nores and B , calculated on the grid and interpolated to the particle positions. E ∗ nores is the corrected electric field, which the particles experience.  \nE  \n(  \nr  \n/u1D45D  \nr  \n/u1D454  \n)  \n)  \n/uni0394  \n/u1D466  \n/uni0394  \n/u1D466', '2.9 Hyper-resistivity': "Hyper-resistivity is introduced in Hybrid PIC codes to remove the propagation of grid-scale fluctuations by damping magnetic field noise on grid-cell scales (hence the use of a 4th-order spatial derivative). To do this, we add a hyper-diffusivity ( 𝜂 hyper ) term on the right-hand side of the Ohm's law, Eq. 11,  \nE = ( J -J i ) × B 𝜌 i -GLYPH<16> 𝑘 B 𝑇 e 𝑞 e GLYPH<17> ∇ 𝜌 i 𝜌 i + 𝜇 0 𝜂 J -𝜇 0 𝜂 hyper ∇ 2 J . (34)  \nThe choice for the resistivity and hyper-resistivity depends on the physical problem and the grid resolution (for more details, see Sec. 2.5.2 in Achikanath Chirakkara et al. 2024).", '2.10 Turbulent driving': 'To drive turbulence, we add a turbulent acceleration field f to Eq. 3,  \nd v 𝑝 d 𝑡 = 𝑞 i 𝑚 i GLYPH<0> E nores + v 𝑝 × B GLYPH<1> + f . (35)  \nThe turbulent acceleration field f is constructed in Fourier space and evolved by an Ornstein-Uhlenbeck process, using TurbGen (Federrath et al. 2010, 2022). We drive turbulence on large scales in our numerical simulations. This turbulent energy injection is restricted to large scales (usually half the computational domain). This energy then self-consistently cascades to smaller and smaller scales, where it is ultimately dissipated as heat (Frisch 1995; Federrath et al. 2021). Thus, the injection of turbulent energy via driving leads to an increase in the temperature of the ions.  \nThe turbulent velocity fluctuations excited via driving can be quantified in terms of the sonic Mach number,  \nM = 𝑣 turb 𝑣 th , (36)  \nwhere 𝑣 turb is the turbulent speed and 𝑣 th is the thermal speed. As a result of energy dissipation on small scales, the temperature increases in the absence of cooling, and therefore, the thermal speed increases and the Mach number of the plasma decreases. This means that statistically-steady turbulence cannot be maintained in the absence of cooling, which is problematic, because the properties of the turbulence are sensitive to the Mach number (Federrath et al. 2011; Seta & Federrath 2021; Achikanath Chirakkara et al. 2021; Seta & Federrath 2022). This problem is even more significant for turbulence in the supersonic regime ( M > 1), where 𝑣 th increases faster compared to the subsonic regime ( M < 1), as more turbulent energy is dissipated in the supersonic case, due to dissipation in shocks.', '2.11 Plasma cooling': "As explained in the previous section, maintaining isothermal plasma conditions is critical for certain physics applications and controlled numerical experiments. Therefore, to maintain an isothermal plasma, we introduce a cooling method, with a schematic diagram shown in Fig. 4. The goal of this cooling method is to achieve a constant (in space and time) target ion thermal speed,  \n𝜎 target = 𝑁 𝑑 𝑘 B 𝑇 target 𝑚 i ! 1 / 2 , (37)  \nwhere 𝑇 target is the target ion temperature and 𝑁 𝑑 = 1 , 2 or 3, depending on the dimensionality of the problem. We perform this cooling on a user-controlled timescale, ( Δ 𝑡 ) cool (in units of 𝑡 cool ), which is  \ntested and calibrated in Sec. 5. When 𝑇 target = 𝑇 e, cooling ensures that electrons and ions have the same temperature.  \nWe begin by noting that the total particle velocity can be decomposed as  \nv = ( v blk ) G → P + v th , (38)  \ni.e., the sum of the local bulk flow velocity, ( v blk ) G → P , and the particle thermal velocity, v th . To cool the plasma, we first perform P → G interpolations using Eq. 20 with 𝑄 ( r 𝑝 ) = 𝑚 i v 𝑝 and 𝑄 ( r 𝑝 ) = 𝑚 i v 2 𝑝 , where both quantities are vectors with each component representing the 3 spatial coordinate directions, i.e., v 𝑝 = ( 𝑣 𝑝,𝑥 , 𝑣 𝑝,𝑦 , 𝑣 𝑝,𝑧 ) and v 2 𝑝 = ( 𝑣 2 𝑝,𝑥 , 𝑣 2 𝑝,𝑦 , 𝑣 2 𝑝,𝑧 ) , respectively. Following this, the two vector quantities are interpolated back from the computational grid to the particle locations as ( 𝑚 i v ) G → P and ( 𝑚 i v 2 ) G → P , respectively. The ion mass density on the computational grid, 𝜌 𝑚 , is also interpolated from the grid to the particle positions (G → P) as ( 𝜌 𝑚 ) G → P using Eq. 20 with 𝑄 ( r 𝑔 ) = 𝜌 𝑚 . These back-and-forth interpolation operations are performed to obtain the bulk velocities and thermal speed at the particle positions. At each particle position, we then have  \n( v blk ) G → P = ( 𝑚 i v ) G → P ( 𝜌 𝑚 ) G → P , (39)  \nUsing the above, we can calculate  \n( v 2 ) G → P = ( 𝑚 i v 2 ) G → P ( 𝜌 𝑚 ) G → P . (40)  \nσ G → P = GLYPH<16> ( v 2 ) G → P - ( v blk ) 2 G → P GLYPH<17> 1 / 2 , (41)  \nwhere σ G → P is the thermal speed in the 𝑥 , 𝑦 and 𝑧 directions at the particle locations.  \nHaving gathered all this information, we are now ready to perform the cooling step. From Eq. (38), each particle's thermal velocity, v th = v - ( v blk ) G → P , which is re-scaled (corrected) by the proportional mismatch between the measured σ G → P from Eq. (41) and the target thermal speed, then the bulk speed is reattached, such that the total particle velocity after the cooling step is set to  \nv → v cooled = 𝜎 target | σ G → P | v th + ( v blk ) G → P = 𝜎 target | σ G → P | ( v - ( v blk ) G → P ) + ( v blk ) G → P . (42)  \nThe cooling is performed every ( Δ 𝑡 ) cool , in units of 𝑡 cool , which is defined as  \nThis re-scaling adjusts the thermal velocity component of each particle, while keeping the direction of the thermal velocity component at each particle position unchanged (i.e., using the magnitude, | σ G → P | of the three σ G → P components in Eq.(42), to maintain a constant target plasma temperature.  \n𝑡 cool = 𝑡 th M , (43)  \nusing the thermal crossing time, 𝑡 th = Δ 𝑙 / 𝑣 th , where Δ 𝑙 = min ( Δ 𝑥, Δ 𝑦, Δ 𝑧 ) is the minimum side length of a computational grid cell and 𝑣 th is the thermal speed. Cooling is done after the predictorpredictor-corrector integration scheme (see Sec. 2.7) is complete. This method ensures isothermal plasma conditions locally (on the scale of individual grid cells), and by extension also globally.", '2.12 Time-step constraints': 'To resolve the different physical timescales involved in weakly collisional plasmas, the following time-step constraints must be fulfilled:  \n<!-- image -->  \n<!-- image -->  \nCooling  \n-----→  \n(b)  \nFigure 4. Schematic diagram of the cooling method. A single grid cell with 3 particles is shown for simplicity. (a) Using interpolation operations, the particle velocity is decomposed into the bulk, ( v blk ) G → P, and the particle thermal velocity, v th, as described in Eq. 38. (b) Then v th = v - ( v blk ) G → P, is rescaled to the target value, 𝜎 target, and ( v blk ) G → P reattached (not shown in the schematic), as described in Eq. 42. This method ensures isothermal conditions locally and globally throughout the plasma.  \n<!-- image -->  \n(i) Larmor gyration: to resolve the Larmor motion of ions in the presence of a magnetic field, one needs to use a time step that is significantly smaller than the timescale required for a charged particle to perform one gyration about the magnetic field,  \n𝑡 Larmor = 2 𝜋𝑚 i 𝑞 i | B | , (44)  \nwhere 𝑚 i and 𝑞 i are the mass and charge of the ions, respectively. The required time step to resolve Larmor gyration is problem-dependent, and will be tested in Sec. 3.1.  \n(ii) Particle speed: the particle speeds set a time-step constraint determined by the time taken by the fastest moving particle in the computational domain to travel across the grid cell of minimum side length Δ 𝑙 = min ( Δ 𝑥, Δ 𝑦, Δ 𝑧 ) . This is given by  \n𝑡 part = Δ 𝑙 | v | max , (45)  \n(iii) Alfvén waves: the timescale for the propagation of Alfvén waves across a computational grid with cells of minimum side length Δ 𝑙 is  \nwhere | v | max is the speed of the fastest particle within the domain.  \n𝑡 Alfvén = Δ 𝑙 𝑣 A , (46)  \nwith the Alfvén speed 𝑣 A = | B |/ √ 𝜇 0 𝑚 i 𝑛 i , where 𝑛 i and 𝜇 0 are the number density and the vacuum magnetic permeability, respectively.  \n(iv) Whistler waves: the timescale for the propagation of whistler waves across a grid cell is  \n𝑡 whistler = Δ 𝑙 𝑣 whistler , (47)  \nwhere 𝑣 whistler = 𝜋 | B |/( 𝑞 i 𝜇 0 𝑛 i Δ 𝑙 ) is the speed of whistler waves with wavenumber, 𝑘 = 2 𝜋 / 2 Δ 𝑙 . Note that 𝑡 whistler ∝ Δ 𝑙 2 .  \n(v) Ohmic resistivity: the magnetic diffusivity of the plasma, 𝜂 , sets the resistive timescale  \n𝑡 𝜂 = ( Δ 𝑙 ) 2 𝜂 . (48)  \n(vi) Hyper-resistivity: the higher-order numerical diffusivity sets an additional dissipative time step, the hyper-resistive timescale  \n𝑡 𝜂 hyper = ( Δ 𝑙 ) 4 𝜂 hyper , (49)  \nwhere 𝜂 hyper is the hyper-diffusivity.  \n(vii) Thermal crossing: the thermal crossing timescale is the time taken by a sound wave to cross a grid cell, defined as  \n𝑡 th = Δ 𝑙 𝑣 th , (50)  \nwhere 𝑣 th is the thermal speed of the plasma.  \nThe overall time step is calculated as the minimum of the all above time-step constraints,  \nΔ 𝑡 = min { 𝑠 L 𝑡 Larmor , 𝑠 part 𝑡 part , 𝑠 A 𝑡 Alfvén , 𝑠 w 𝑡 whistler , 𝑠 𝜂 𝑡 𝜂 , 𝑠 𝜂 hyper 𝑡 𝜂 hyper , 𝑠 th 𝑡 th } , (51)  \nwhere 𝑠 L , 𝑠 part, 𝑠 A , 𝑠 w, 𝑠 𝜂 , 𝑠 𝜂 hyper , and 𝑠 th , are the safety factors for the Larmor, particle speed, Alfvén, whistler, resistive, hyper-resistive, and thermal time steps, respectively. We introduce these safety factors to ensure that each time step is sufficiently resolved. The default value for all the safety factors is 0.1, except for the particle-speed time step, which is 𝑠 part = 0 . 5. Depending on the physical problem, these safety factors can be adjusted.', '3 PARTICLE, WAVE, AND WAVE-PARTICLE INTERACTION TESTS': 'To test the accuracy of our code, we perform the following standard tests and compare our numerical results with the analytical solutions. First, we test the Larmor gyration of charged particles in Sec. 3.1. In Sec. 3.2 we test the propagation of Alfvén and whistler waves. Finally, we test whether the code can capture physical instabilities arising from particle-wave interactions in Sec. 3.3.', '3.1 Charged particle in a uniform magnetic field': "We test our implementation of the Boris integrator by studying the motion of a charged particle in a uniform magnetic field. We set up an external (non-evolving) magnetic field, B = 𝐵 0 ˆ 𝑧 , where 𝐵 0 = 10 -8 Tesla, and give the charged particles an initial velocity perpendicular to the direction of the magnetic field, v = 𝑣 0 ˆ 𝑥 , where 𝑣 0 = 1 . 0 m s -1 . Due to the effect of the Lorentz force, the charged particles will gyrate in a circle in the 𝑥 -𝑦 plane with the Larmor radius, 𝑟 Larmor = 𝑚 i 𝑣 0 / 𝑞 i 𝐵 0 , where 𝑚 i and 𝑞 i are the mass and charge of the particle, respectively. The analytical solution for the evolution of the 𝑥 -coordinate of the particle's position is 𝑥 ( 𝑡 ) = 𝑟 Larmor sin ( Ω i 𝑡 ) ,  \nFigure 5. The 𝑥 -coordinate of a charged particle undergoing gyrations in the presence of a uniform magnetic field, normalised to the Larmor radius ( 𝑟 Larmor), as a function of time in units of the Larmor time ( 𝑡 Larmor). We plot this for simulations with Δ 𝑡 = 0 . 1 𝑡 Larmor for the Boris integrator and a non-conserving particle integrator. The solid black line depicts the analytical solution. We find that the gyration radius of our numerical solutions does not change with the Boris integration method, while with the non-conserving particle integrator, the gyration radius increases over time.  \n<!-- image -->  \nwhere Ω i = 𝑞 i 𝐵 0 / 𝑚 i is the gyration frequency of the charged particle. The total kinetic energy of the particle does not change as the magnetic field does not do any work on the particle.  \nWe pick the time step, Δ 𝑡 = 0 . 1 𝑡 Larmor , where 𝑡 Larmor = 2 𝜋 / Ω i is the time taken by the charged particle to complete one gyration. Fig. 5 shows the 𝑥 -coordinate of the particle, normalised to 𝑟 Larmor , as a function of time, normalised to the Larmor time ( 𝑡 Larmor , defined in Eq. 44), along with the analytical solution. We also perform this test with a non-conserving cyclic leapfrog particle integrator (Holmstrom 2009).  \nWefindthat the amplitude (radius) of the Larmor gyration is stable in the Boris integration method, as the energy is conserved up to the machine precision level in this case. However, with a non-conserving particle integrator (dotted red line), the radius of gyration increases over time, due to the particle gaining energy as a result of integration inaccuracies. Thus, we conclude that using the Boris integrator provides the basic requirement for stably modelling Larmor gyration, and thus we use this integration scheme as the default for all PIC simulations.  \nWhile the Boris integrator provides stable Larmor (circular) orbits for any choice of time step Δ 𝑡 , for a uniform magnetic field it may not do so for a spatially varying magnetic field, and instead, the choice of time step may play a role. This is tested next, where we follow the motion of a charged particle in a non-uniform magnetic field,  \n𝐵 𝑧 ( 𝑥 ) = 𝐵 min GLYPH<26> 2 -exp GLYPH<20> -GLYPH<16> 𝑥 2 Δ 𝑥 GLYPH<17> 2 GLYPH<21> GLYPH<27> , (52)  \nwhere 𝑥 is the 𝑥 -coordinate of the grid cell centre, Δ 𝑥 is the cell length in the x-direction and 𝐵 min is the magnetic field strength at 𝑥 = 0. This profile is shown by the colour map in Fig. 6. The charged particle is given an initial velocity in the 𝑥 -direction, similar to the test in Fig. 5, and we plot the results up to 1 𝑡 Larmor , where 𝑡 Larmor is calculated based on the maximum magnetic field strength, 𝐵 max = 2 𝐵 min . We vary the time step, Δ 𝑡 / 𝑡 Larmor = 10 -3 , 10 -2 and 10 -1 , and show the results in Fig. 6. The near-perfect orbit is obtained for a very small time step-here Δ 𝑡 / 𝑡 Larmor = 10 -3 (dashed line). For a relatively  \nFigure 6. Orbit of a charged particle in a spatially varying magnetic field, 𝐵 𝑧 ( 𝑥 ) , described by Eq. 52, using the Boris integrator for varying time steps, Δ 𝑡 / 𝑡 Larmor = 10 -3 , 10 -2 , and 10 -1 . Using a relatively large Δ 𝑡 / 𝑡 Larmor = 10 -1 , we find that the trajectory of the particle deviates significantly from the cases that use Δ 𝑡 / 𝑡 Larmor = 10 -2 and 10 -3 . While the required time step depends on the physical problem, this test shows that it is important to pick a time step that can suitably resolve the trajectory of the particles.  \n<!-- image -->  \nlarge time step, Δ 𝑡 / 𝑡 Larmor = 10 -1 (blue crosses), we see that the orbit is significantly distorted, with a clear orbital precision. Using Δ 𝑡 / 𝑡 Larmor = 10 -2 (green pluses), the trajectory is fairly close to the optimal orbit. This test demonstrates that it is important to choose a suitable time step to achieve the desired accuracy in the trajectories of the charged particles, especially when the magnetic field is nonuniform, which is almost always the case for PIC applications.", '3.2 Propagation of waves in collisionless plasma': 'By using perturbation analysis on the Hybrid-kinetic equations, one can obtain the dispersion relation for the propagation of circularly polarized waves. Right-handed waves follow the dispersion relation 𝜔 2 R = 𝑘 2 ( 1 + 𝜔 R ) andleft-handed waves follow the dispersion relation 𝜔 2 L = 𝑘 2 ( 1 -𝜔 L ) (Kulsrud 2005). We note that in these equations, the frequency, 𝜔 , is normalised to the ion gyration frequency, Ω i , and the wavenumber, 𝑘 , is normalised to the ion inertial length, 𝑑 i . In the low-wavenumber limit, 𝑘 ≪ 1, both the right- and lefthanded waves follow the Alfvén wave dispersion relation, 𝜔 2 = 𝑘 2 . In the high-wavenumber limit, 𝑘 ≫ 1, the dispersion relation for the right-handed waves is that of whistler waves, 𝜔 R = 𝑘 2 , and for the left-handed waves, it is that of the ion-cyclotron waves, 𝜔 L → 1 asymptotically as 𝑘 →∞ .  \nTo test this with our numerical code, we initialise a mean magnetic field, 𝐵 0 , along the 𝑥 -axis and oscillatory magnetic fields along the 𝑦 and 𝑧 -axes, 𝐵 𝑦 = 𝛿𝐵 sin ( 𝑘𝑥 ) and 𝐵 𝑧 = 𝛿𝐵 cos ( 𝑘𝑥 ) , where 𝑘 = 2 𝜋 / 𝐿 𝑥 , 𝛿𝐵 = 10 -3 𝐵 0 and 𝐿 𝑥 is the length of the computational domain in the 𝑥 -direction. Our test simulations have cold ions and electrons, 𝑇 i = 𝑇 e = 0. The Ohmic dissipation and hyper-resistivity have been switched off in these tests. We choose 𝑁 ppc = 100 particles per cell and periodic boundary conditions. We also resolve the Alfvén and whistler wave time steps and use a safety factor of 0.1 for both these time steps (see Sec. 2.12 for details).  \nFig. 9 shows the squared amplitude of the density fluctuations as a function of time, normalised to the sound-crossing timescale, 𝑡 𝑐 s = 1 /( 𝑘𝑐 s ) , where 𝑐 s is the sound speed based on 𝑇 e. We set the ratio of the ion to electron temperature, 𝑇 i / 𝑇 e = 0 . 1. We find that the density fluctuations decay exponentially with time, by locating the local maxima of the density fluctuations and fitting an exponential decay  \n<!-- image -->  \ni  \nFigure7. Dispersion relationship of Alfvén, whistler and ion-cyclotron waves, i.e., frequency of the propagating wave, 𝜔 , normalised to the ion gyration frequency, Ω i, for different values of the wavenumber 𝑘 (in units of the ion inertial length, 𝑑 i), for a wide range of wavenumbers, 𝑘𝑑 i = [ 0 . 05 , 9 ] with 𝑘 Δ 𝑥 = 2 𝜋 / 64 (i.e. 64 grid cells per wavelength), where Δ 𝑥 is the grid-cell length in the direction of wave propagation. The red points show the numerical results with our new code, AHKASH . The black dashed curves show the analytical solutions for the left-handed and right-handed waves. In the low-wavenumber regime ( 𝑘𝑑 i ≪ 1), the Alfvén wave dispersion relation ( 𝜔 / Ω i = 𝑘𝑑 i) is followed and at higher wavenumbers ( 𝑘𝑑 i ≫ 1), the whistler waves ( 𝜔 / Ω i = ( 𝑘𝑑 i ) 2 ) and the ion-cyclotron waves ( 𝜔 = Ω i) appear. The measured errors in the frequencies are very small ( ≲ 4%). For the tests with 𝑘𝑑 i = 0 . 1 , 0 . 9 and 9, we vary the grid resolution, 𝑘 Δ 𝑥 = 2 𝜋 / 32 and 2 𝜋 / 128 and find that the numerical solutions are converged with 𝑘 Δ 𝑥 = 2 𝜋 / 32 in different 𝑘𝑑 i regimes. We find that our results are in excellent agreement with the analytical solutions and show stable and accurate propagation of waves.', '3.2.1 Wavelength dependence': 'The initial wave propagates across the computational domain and to accurately extract the frequencies of the waves, we perform a Fourier analysis on the time evolution data. We then use a lognormal fitting function to fit the frequency peaks in the data. The frequency we report is the mean of the lognormal distribution and the standard deviation of the distribution is a measure of its uncertainty. In Fig. 7, we show the frequency of the waves normalised to the ion gyration frequency, 𝜔 / Ω i , for different values of the wavenumber normalised to the ion inertial length, 𝑘𝑑 i = [ 0 . 05 , 9 ] . We also plot the analytical solutions for the dispersion relation and find excellent agreement with the analytical solutions across all 𝑘𝑑 i regimes. This shows that our code can stably and accurately capture the propagation of waves.  \nWe pick 𝑘 Δ 𝑥 = 2 𝜋 / 64 in our tests, where Δ 𝑥 is the grid-cell length in the direction of wave propagation and 𝑘 = 2 𝜋 / ℓ , where ℓ is the wavelength. This implies there are 64 cells per wavelength, ℓ = 64 Δ 𝑥 . The tests with 𝑘𝑑 i = 0 . 1 , 0 . 9 and 9 are repeated with varying grid resolution, Δ 𝑥 = ℓ / 32 or 𝑘 Δ 𝑥 = 2 𝜋 / 32 and Δ 𝑥 = ℓ / 128 or 𝑘 Δ 𝑥 = 2 𝜋 / 128. In Fig. 7, we show that the numerical results with 𝑘 Δ 𝑥 = 2 𝜋 / 32 are converged in different 𝑘𝑑 i regimes.', '3.2.2 Grid resolution criteria': 'In this section, we study in more detail the effect of the grid resolution on the wave test. We vary the grid-cell length in the direction of wave propagation, Δ 𝑥 = ℓ / 8 ℓ / 128, where ℓ = 2 𝜋 / 𝑘 is the wavelength, and measure the frequency of wave propagation, 𝜔 , normalised to the ion gyration frequency, Ω i , for the test with 𝑘𝑑 i = 0 . 9. The  \nFigure 8. The frequency of the propagating wave, 𝜔 , normalised to the theoretical frequency, 𝜔 theory, for the wave test with 𝑘𝑑 i = 0 . 9, as a function of grid resolution, ℓ / Δ 𝑥 , where ℓ = 2 𝜋 / 𝑘 is the wavelength and Δ 𝑥 is the grid-cell length. For tests with Δ 𝑥 = ℓ / 8 and Δ 𝑥 = ℓ / 16, the analytical solutions are not obtained for the right-handed waves. As the grid resolution is increased Δ 𝑥 ≳ ℓ / 32, the analytical solutions are recovered, and the numerical solutions are converged.  \n<!-- image -->  \nresults are shown in Fig. 8. We find that with ℓ = ( 32 - 128 ) Δ 𝑥 , the analytical solutions are recovered (to within one sigma) and the numerical solutions are converged with ℓ = 32 Δ 𝑥 . Using a lower grid resolution, the analytical solutions are not obtained, and we also find spurious frequencies start to dominate in the Fourier analysis.', '3.3 Wave-particle interaction': 'Here we test the decay of ion acoustic waves due to Landau damping using the 𝛿 𝑓 method (see Sec. 2.5). Landau damping of ion acoustic waves follows the dispersion relation 𝑑𝑍 ( 𝜁 )/ 𝑑𝜁 = 2 𝑇 i / 𝑇 e, where 𝑍 ( 𝜁 ) is the plasma dispersion function and 𝜁 is the phase velocity, which has a real (oscillatory) and an imaginary (decaying) part, and 𝑇 i and 𝑇 e are the ion and electron temperature, respectively. The decay rate of the wave depends on the ratio of the ion to electron temperature. We will compare our numerical solutions with the solutions of the analytical dispersion relation, 𝑑𝑍 / 𝑑𝜁 = 2 𝑇 i / 𝑇 e, which we solve numerically using the Newton-Raphson method for a given 𝑇 i / 𝑇 e.  \nThe length of our computational domain is 𝐿 𝑥 = 16 𝑑 i along the 𝑥 -axis, where 𝑑 i is the ion inertial length. We initialise the velocity of the ions (particles) using a Maxwell-Boltzmann distribution with temperature 𝑇 i and fixed electron temperature 𝑇 e. We initialise the density with the profile 1 + 𝑎 cos ( 𝑘𝑥 ) , where we choose a small perturbation amplitude, 𝑎 = 0 . 01 and the wavenumber of the perturbation, 𝑘 = 2 𝜋 / 𝐿 𝑥 . We used 𝑁 ppc = 10 6 particles per cell and periodic boundary conditions.', '3.3.1 Results': 'Figure 9. Landau damping of ion acoustic waves shown by the decay of the squared amplitude of the density fluctuation as a function of time, normalised to the sound-crossing time, 𝑡 𝑐 s , for a ratio of the ion temperature to the electron temperature of 𝑇 i / 𝑇 e = 0 . 1. The local maxima used for fitting are shown by the red dots. We fit for the Landau damping decay rate, shown by the dashed line and find a decay rate of 𝛾 fit = ( 4 . 5 ± 0 . 2 ) × 10 -2 𝑡 -1 𝑐 s , which shows good agreement with the theoretical value, 𝛾 theory.  \n<!-- image -->  \nFigure 10. SameasFig. 9, but for 𝑁 ppc = 10 2 , 10 3 , 10 4 , 10 5 and 10 6 particles per cell. For 𝑁 ppc ≳ 10 4 , the Landau damping of ion acoustic waves are resolved, while particle noise dominates below 𝑁 ppc ∼ 10 4 .  \n<!-- image -->  \nfunction. We find a decay rate of 𝛾 fit = ( 4 . 5 ± 0 . 2 ) × 10 -2 𝑡 -1 𝑐 s . We measure the decay rates between successive data points and the decay rate we report is the mean of all the intervals and the error reported is the standard deviation of the decay rate across all the intervals. We compare this with the solutions of the dispersion relation and find good agreement with the theoretical solution, 𝛾 theory = 0 . 042 𝑡 -1 𝑐 s .', '3.3.2 Particle resolution criterion': 'Here we repeat the Landau damping test with 𝑁 ppc = 10 2 , 10 3 , 10 4 and 10 5 particles per cell. The results are shown in Fig. 10. We find that with low particle count 𝑁 ppc ≲ 10 4 , particle noise dominates, while for higher 𝑁 ppc the Landau damping of ion acoustic waves are resolved.', '4 PARTICLE-GRID INTERPOLATION TESTS': 'After confirming that the code agrees with analytical solutions for standard particle, wave, and particle-wave problems, we perform collisionless turbulent dynamo tests to understand the effect of using different interpolation kernels and smoothing in our numerical simulations. The turbulent dynamo mechanism converts turbulent kinetic energy into magnetic energy, exponentially amplifying magnetic fields (Kazantsev 1968; Moffatt 1978; Ruzmaikin et al. 1988). The properties of the turbulent dynamo are sensitive to the Mach number, nature of turbulence driving and the magnetic and kinetic Reynolds number of the plasma (Schekochihin et al. 2004; Haugen et al. 2004a,b; Federrath et al. 2011; Schober et al. 2012; Seta et al. 2020; Seta & Federrath 2020; Achikanath Chirakkara et al. 2021; Kriel et al. 2022; Gent et al. 2023) and recently numerical studies have extended this to the collisionless regime using Hybrid-kinetic codes (Rincon et al. 2016; St-Onge & Kunz 2018; Achikanath Chirakkara et al. 2024).', '4.1 Turbulence driving setup and initial conditions': 'In all the numerical tests presented in this and the following sections, we drive turbulence with a random acceleration field, f , following an Ornstein-Uhlenbeck process (Eswaran & Pope 1988; Federrath et al. 2010, 2022). The acceleration field is constructed in Fourier space to contain only large-scale modes, 𝑘𝐿 / 2 𝜋 = [ 1 , 3 ] , where 𝐿 is the size of the cubic computational domain. The turbulence driving amplitude is modelled as a parabolic function that peaks at 𝑘𝐿 turb / 2 𝜋 = 2 and is zero at 𝑘𝐿 / 2 𝜋 = 1 and 3, where 𝐿 turb is the peak driving scale of the turbulence. The turbulent velocity dispersion, 𝑣 turb , is controlled by adjusting the amplitude of the driving. We inject purely solenoidal acceleration modes (∇· f = 0 ) (Federrath et al. 2010) into the plasma for these tests. For all the tests in this section, we initialize the velocity of the particles using a Maxwellian distribution with temperature 𝑇 i , which is equal to 𝑇 e.  \nAll the numerical tests in this section are performed on a periodic three-dimensional computational domain with 𝑁 3 grid = 128 3 grid cells. We drive the turbulence to produce a steady-state (fullydeveloped) turbulent Mach number of M = 0 . 2. Using the steadystate kinetic energy, we initialise the magnetic field such that the ratio of magnetic to turbulent kinetic energy is ( 𝐸 mag / 𝐸 kin ) init = 10 -8 , and the initial ratio of the Larmor radius to the box size is ( 𝑟 Larmor / 𝐿 ) init = 10 2 (see Sec. 2.5 of Achikanath Chirakkara et al. (2024) for further details). The convergence of the turbulent dynamo simulations with respect to the grid and particle resolution is studied in App. D of Achikanath Chirakkara et al. (2024). The magnetic Reynolds number of the plasma, defined as  \nRm = 𝑣 turb 𝐿 turb 𝜂 , (53)  \nis set to Rm = 500. The hyper-resistive magnetic Reynolds number is defined as  \nRm h = 𝑣 turb ( 𝐿 turb ) 3 𝜂 hyper , (54)  \nand is set to Rm h = 10 7 . For all the tests in this section, we use cooling (tested and discussed separately in detail in Sec. 5 below) with a cooling frequency of ( Δ 𝑡 ) cool = 0 . 1 𝑡 cool (see Eq. 43 and Sec. 2.11).  \nFigure 11. Particle → grid interpolation tests using collisionless turbulent dynamo simulations. The Mach number, M , is shown in the top panel, and the magnetic-to-kinetic energy ratio, 𝐸 mag / 𝐸 kin, is shown in the bottom panel, as a function of time normalised to the turbulent eddy turnover time, 𝑡 0. For the standard set of numerical parameters, we use the cloud-in-cell (CIC) interpolation kernel with 𝑁 ppc = 100 particles per cell and two smoothing passes, 𝑁 smooth = 2. We find that using nearest-grid-point (NGP) or triangularshaped-cloud (TSC) interpolation functions does not noticeably change the box-averaged solutions for the Mach number or the magnetic energy (provided 𝑁 smooth = 2). We also find that changing the number of particles per cell, 𝑁 ppc = 50 , 200, or 400, does not change the box-averaged solutions. Changing the number of smoothing operations when 𝑁 smooth ≥ 1 also has no significant impact on the box-averaged solutions. However, without any smoothing, 𝑁 smooth = 0 (red dashed line), the growth in magnetic energy is slightly underestimated.  \n<!-- image -->', '4.2 Comparison of global quantities in interpolation tests': "We test three interpolation kernels: the nearest-grid-point (NGP), the cloud-in-cell (CIC), and the triangular-shaped-cloud (TSC) interpolation functions, detailed in Sec. 2.4.1 and Appendix A. We also vary the number of smoothing passes, 𝑁 smooth = 0 , 1 , 2, and 4, as discussed in Sec. 2.4.2. Finally, we investigate the influence of different numbers of particles per cell, 𝑁 ppc = 50 , 100 , 200, and 400.  \nFig. 11 shows the time evolution (in units of the turbulent turnover time, 𝑡 0 = 𝐿 turb / 𝑣 turb ) of the sonic Mach number ( M ; top panel) and the magnetic-to-kinetic energy ratio ( 𝐸 mag / 𝐸 kin ; bottom panel), for all the interpolation tests considered here. We see that all test simulations yield very similar results, with the sonic Mach number reaching a steady state of M ∼ 0 . 2. We also see that after an initial transient phase of about 2-3 𝑡 0 , the magnetic energy grows exponentially in all tests. These are typical characteristics of small-scale turbulent dynamogrowth during the so-called 'kinematic phase' (Brandenburg & Subramanian 2005; Schekochihin et al. 2004; Federrath 2016; Seta", '& Federrath 2020; Kriel et al. 2022; Achikanath Chirakkara et al. 2024).': 'Comparing the different tests in Fig. 11, we find that all runs are virtually identical in terms of their box-averaged global evolution, except for the simulation that does not use any smoothing ( 𝑁 smooth = 0) after particle → grid interpolations (red dashed line).', '4.3 Spatial structure of interpolated fields': 'Fig. 12 shows the relative charge density at 𝑡 = 10 𝑡 0 for all the interpolation test simulations above. The top left panel shows the test with the standard set of numerical parameters (CIC, 𝑁 smooth = 2, 𝑁 ppc = 100), which serves as the main reference simulation. Comparing the tests in the first row (CIC vs. NGP vs. TSC), we find no significant difference in the spatial distribution of the charge density. While the NGP method on its own is generally inferior compared to the CIC and TSC interpolation methods, the fact that we perform 𝑁 smooth = 2 smoothing passes after these interpolations nullifies any specific differences between NGP, CIC, and TSC.  \nFinally, the last row of Fig. 12 shows the effect of varying the number of particles per cell, 𝑁 ppc = 50, 200, 400. While using a higher 𝑁 ppc improves the spatial smoothness of the charge density field, we found no significant impact on the global evolution of the system in Fig. 11. It is worth noting that the test with 𝑁 smooth = 2 and 𝑁 ppc = 400 (bottom right panel) is visually very similar to the case with 𝑁 smooth = 4 and 𝑁 ppc = 100 (middle row, righthand panel). This suggests that a balanced combination of 𝑁 ppc and 𝑁 smooth can lead to comparably good results for a desirably low 𝑁 ppc, as it saves computational resources. However, in general, we suggest that the number of smoothing passes should be kept sufficiently low, in particular in supersonic flows, which would result in excessive smearing of discontinuities (shocks). An optimal setting of 𝑁 smooth and 𝑁 ppc is therefore problem-dependent and requires experimentation and testing in any PIC simulation.  \nThe 2nd row of Fig. 12 shows a comparison of 𝑁 smooth = 0 (no post-interpolation smoothing), 𝑁 smooth = 1, and 4. We see that not smoothing after interpolation (middle row, left-hand panel) leaves a very noisy density distribution. This substantially improves already with a single smoothing pass (central panel). Using 4 smoothing passes ( 𝑁 smooth = 4; right panel) leads to a relatively smooth distribution. However, neither 𝑁 smooth = 1 nor 4 are substantially different from the standard 𝑁 smooth = 2 results. While the spatial details of the interpolated fields depend on the number of smoothing passes, we saw in Fig. 11 that the choice of smoothing passes (except for 𝑁 smooth = 0) has no significant impact on the global evolution of the system.', '5 COOLING TESTS': 'In this section, we demonstrate the working of our cooling method for ions described in Sec. 2.11.', '5.1 Cooling vs. no cooling': 'Here we investigate the requirement of plasma cooling to model subsonic and supersonic turbulent flows. We run the following numerical simulations on a triply periodic uniform computational domain with 𝑁 3 grid = 128 3 grid cells and 𝑁 ppc = 100 particles per cell for subsonic, M target = 0 . 2, and supersonic, M target = 2 target Mach numbers, without and with cooling. The turbulence driving used here is described in Sec. 4.1. We note that magnetic fields have not been  \nFigure 12. Slice plots of the relative charge density (colour map) at 10 𝑡 0 for the interpolation test simulations shown in Fig. 11. The superimposed vectors and streamlines show the velocity and magnetic field, respectively. For the standard test (top left panel), we use the cloud-in-cell (CIC) interpolation kernel with 2 smoothing passes ( 𝑁 smooth = 2), and 𝑁 ppc = 100 particles per cell. First row: comparison of different interpolation kernels, including the nearest-grid-point (NGP) and triangular-shaped-cloud (TSC) interpolation kernels, with 𝑁 smooth = 2 and 𝑁 ppc = 100. Second row: tests without smoothing ( 𝑁 smooth = 0), with 𝑁 smooth = 1, and 𝑁 smooth = 4, using the CIC interpolation kernel and 𝑁 ppc = 100. Third row: tests with different numbers of particles per cell, 𝑁 ppc = 50 , 200, and 400. These tests show that it is important to smooth the particle-sampled fields with at least 𝑁 smooth = 1 smoothing pass (somewhat better results still, are obtained here with 𝑁 smooth = 2), and to work with a sufficiently large number of particles per cell ( 𝑁 ppc ∼ 100) to limit the effects of particle interpolation noise.  \n<!-- image -->  \nincluded for all the tests described in this section as they play no role in the cooling functionality.', '5.1.1 Time evolution of global thermodynamic and turbulent quantities': 'The time evolution of the cooling experiments is shown in Fig. 13. The top panel shows the turbulent speed, normalised to the fixed target thermal speed, 𝑣 turb /( 𝑣 th ) target, in our numerical tests as a  \nfunction of time. We find that without cooling (red curves), it is difficult to maintain the target turbulent speeds, 𝑣 turb /( 𝑣 th ) target = 0 . 2 in the subsonic regime and 𝑣 turb /( 𝑣 th ) target = 2 in the supersonic regime. After an initial maximum, 𝑣 turb decreases in the tests without cooling (red curves). When cooling is included, we find that the target 𝑣 turb is reached and maintained throughout the simulations. The middle panel of Fig. 13 shows the thermal speed, normalised to the target thermal speed, 𝑣 th /( 𝑣 th ) target. For the tests without cooling, we find that 𝑣 th increases monotonically as turbulence energy dissipates,  \nFigure 13. Turbulent speed normalised to the target thermal speed, 𝑣 turb /( 𝑣 th ) target (top panel), thermal speed normalised to the target thermal speed, 𝑣 th /( 𝑣 th ) target (middle panel), and the sonic Mach number, M (bottom panel), as a function of time, normalised to the eddy turn-over time, 𝑡 0, for turbulence driving simulations on 𝑁 3 grid = 128 3 grid cells with 𝑁 ppc = 100 particles per cells with cooling (blue) and without cooling (red), in tests with a target Mach number of 0.2 (dashed lines) and 2 (solid lines), respectively. We find that without cooling it is impossible to maintain the target Mach number. This is primarily the result of heating due to turbulent dissipation and shocks (for the supersonic test), demonstrating the necessity for cooling to model steady-state turbulent flows.  \n<!-- image -->  \nt  \n/  \nt  \n0  \nwhile with cooling, 𝑣 th is maintained at its target level. We also find that the thermal speed of the plasma increases faster in the supersonic regime when cooling is not included. For both the tests with cooling in the subsonic and supersonic regime, ( Δ 𝑡 ) cool = 0 . 1 𝑡 cool , where 𝑡 cool is the cooling timescale (more details on the role of ( Δ 𝑡 ) cool below). The ratio of the data in the top and bottom panels is shown in the bottom panel of Fig. 13, i.e., the sonic Mach number, M , which reflects the issues discussed separately for 𝑣 turb and 𝑣 th without cooling. Over time, the plasma becomes increasingly subsonic, which makes studying maintained supersonic flows practically impossible without cooling. However, even for the subsonic regime, turbulence cannot be maintained at a steady Mach number without cooling. This demonstrates the necessity for cooling in applications such as steady-state turbulence.', '5.1.2 Statistical and spatial distribution of the temperature': 'Hereweinvestigate the distribution and spatial structure of the plasma temperature in more detail. Fig. 14 shows the probability density  \nFigure 14. Probability density functions (PDFs) of the thermal speed, 𝑣 th, normalised to the target thermal speed, ( 𝑣 th ) target, for simulations with target Mach number, M target = 0 . 2 (dashed lines) and M target = 2 (solid lines), without cooling (red) and with cooling at 𝑡 / 𝑡 0 = 0 . 5 (blue) and time-averaged in 𝑡 / 𝑡 0 = [ 2 , 10 ] (cyan). The tests without cooling show very broad distributions, i.e., not only is the target average temperature overestimated, but the plasma is also spatially at very different temperatures. For the tests with cooling, the target thermal speed is reached and maintained (see time-averaged curves shown in cyan) to within 0.6% accuracy in both the subsonic and supersonic regimes. The spatial fluctuations remain reasonably small with a standard deviation of 0.3% and 4.3% in the subsonic and supersonic regimes, respectively. The higher residual fluctuations in the supersonic regime are due to the presence of shocks.  \n<!-- image -->  \nfunction (PDF) of the thermal speed, 𝑣 th , normalised to the target thermal speed, ( 𝑣 th ) target. For the subsonic and supersonic test without cooling at 𝑡 = 0 . 5 𝑡 0 , we find 𝑣 th /( 𝑣 th ) target = 1 . 172 ± 0 . 010 and 1 . 189 ± 0 . 117, respectively. While the mean values for both these tests are similar with about 17-19% higher than the target, the spread of the thermal speed is much higher for the supersonic simulation. This can be attributed to more turbulent kinetic energy available to heat the plasma and the presence of shocks dissipating large amounts of energy in the supersonic regime. By contrast, when cooling is included, we find 𝑣 th /( 𝑣 th ) target = 1 . 001 ± 0 . 003 and 1 . 000 ± 0 . 002 for the same subsonic and supersonic test comparison, respectively. For the simulations with cooling, we also show the PDFs averaged between 2 -10 𝑡 0 , for which we find 𝑣 th /( 𝑣 th ) target = 1 . 000 ± 0 . 003 and 1 . 006 ± 0 . 043, respectively, i.e., the target thermal speeds are reasonably well maintained.  \nThe respective spatial distribution of the thermal speed is shown in Fig. 15. We see that without cooling, the distribution of the thermal speed is inhomogeneous. In the case of supersonic turbulence, this distribution is significantly more inhomogeneous as more energy is injected into the plasma and dissipation occurs primarily in shocks. Comparing the over- and under-densities in the bottom right-hand panels we see that they are associated with the under- and over-heated regions in the bottom left-hand panels, respectively, i.e., adiabatic heating and cooling occur in relatively over- and under-dense regions, respectively. Once cooling is introduced, the resulting thermal speed distribution is fairly homogeneous. This shows that it is important to cool the plasma locally (not just globally), to obtain a spatially homogeneous isothermal plasma.  \n(a) Thermal speed relative to the average thermal speed.  \n<!-- image -->  \n(b) Plasma density relative to the average density.  \n<!-- image -->  \nFigure 15. (a) Slice plots of the thermal speed, 𝑣 th, normalised to the mean thermal speed, ⟨ 𝑣 th ⟩ , at 𝑡 / 𝑡 0 = 0 . 5, for the same test simulations as in Fig. 13 (as indicated in the panel labels). We see that without cooling, the thermal speed is inhomogeneous, particularly in the supersonic regime. Cooling homogenises the thermal speed, showing that cooling is required locally to achieve reasonable isothermal conditions throughout the plasma. (b) Same as (a), but for the ion mass density, 𝜌 m, normalised to the mean ion mass density, ⟨ 𝜌 m ⟩ . The fluctuations in the mass density are much smaller in the subsonic tests when compared to the supersonic tests, as expected. Shocked structures with large density gradients are apparent in the supersonic regime, which correlates with regions of higher thermal speeds, showing that shocks dissipate kinetic energy efficiently.', '5.2 Choosing the cooling frequency': 'In order to keep the plasma isothermal, we need to apply the cooling operator at regular time intervals, ( Δ 𝑡 ) cool ∝ 𝑡 cool (see Sec. 2.11). Here we investigate the required cooling frequency, ( Δ 𝑡 ) -1 cool , to maintain approximate isothermal plasma conditions in the subsonic and supersonic regime of turbulence.', '5.2.1 Subsonic tests': 'In Fig. 16 we vary the cooling time for the subsonic turbulence simulations presented above, ( Δ 𝑡 ) cool = 0 . 01, 0 . 1, 1 and 10 𝑡 cool . We perform these simulations with a target M = 𝑣 turb /( 𝑣 th ) target = 0 . 2 and list the results in Tab. 1. The figure legend also lists the cooling timescale in units of the thermal crossing time, 𝑡 cool = 𝑡 th /M = 5 𝑡 th (see Eq. 43). The middle panel of Fig. 16 shows that for larger values of the cooling time, the thermal speeds start increasing, while for lower values of ( Δ 𝑡 ) cool , 𝑣 th is maintained close to the target thermal speed. We find that for the turbulent speed and the Mach number, shown by the top and the bottom panel of Fig. 16, respectively, the numerical solutions converge and the target 𝑣 turb and M are attained with ( Δ 𝑡 ) cool ∼ 0 . 1 𝑡 cool = 0 . 5 𝑡 th .', '5.2.2 Supersonic tests': 'We repeat the same experiment with a supersonic target Mach number, M = 𝑣 turb /( 𝑣 th ) target = 2, and show the results in Fig. 17. We list the results in Tab. 1. We find that for ( Δ 𝑡 ) cool = 1 and 10 𝑡 cool , the plasma heats up very quickly and the target thermal speed is not maintained. Performing the cooling with a higher frequency, by using ( Δ 𝑡 ) cool = 0 . 1 𝑡 cool = 0 . 05 𝑡 th , the target thermal speed is acquired and maintained and the target Mach number is achieved.  \nTable 1. List of simulations shown in Fig. 16 and Fig. 17 with the corresponding model name, turbulent speed normalised to the target thermal speed, 𝑣 turb /( 𝑣 th ) target, thermal speed normalised to the target thermal speed, 𝑣 th /( 𝑣 th ) target, and the sonic Mach number, M . These quantities are measured for 𝑡 / 𝑡 0 = 2 -10, where 𝑡 0 is the turbulent turnover time. This shows that the thermal speed and Mach number are converged (to within one sigma) as the cooling frequency is increased, in both the subsonic and the supersonic regimes of turbulence.', '5.2.3 Summary': 'The thermal crossing time, 𝑡 th , is fixed between the subsonic and supersonic tests, as both have been configured to reach the same target thermal speed. Since the Mach number is 10 times higher in the supersonic runs, we need a 10 times higher cooling frequency (in absolute terms) in the supersonic tests, compared to the subsonic tests. However, when expressed in terms of the cooling time, which takes the target Mach number into account (Eq. 43), the required cooling frequency is the same in the supersonic and subsonic regimes. That is, we find that for ( Δ 𝑡 ) cool ∼ 0 . 1 𝑡 cool , i.e., using a safety factor of ∼ 0 . 1, the target thermal speed and Mach number are reasonably achieved and maintained throughout.  \n〉  \nth  \n〈  \nth  \n〉  \n〈  \nFigure 16. Same as Fig. 13, but for subsonic turbulence with target M = 0 . 2 and different cooling times, ( Δ 𝑡 ) cool = 0 . 01, 0 . 1, 1 and 10 𝑡 cool, where 𝑡 cool is the cooling timescale. We see that the target thermal speed is overestimated with ( Δ 𝑡 ) cool = 1 and 10 𝑡 cool. By contrast, the thermal speed is converged and maintained close to the target value for ( Δ 𝑡 ) cool ≲ 0 . 1 𝑡 cool.  \n<!-- image -->', '6 CODE OPTIMISATION AND PERFORMANCE': "In this section, we discuss the performance and scalability of our code. We discuss the 'hybrid precision' method introduced in Federrath et al. (2021). We extend and test the method for Hybrid PIC in Sec. 6.1. Parallel scaling tests are provided in Sec. 6.2.", '6.1 Hybrid precision': "Hybrid PIC simulations can be extremely computationally demanding, as they require both high grid-cell counts and high particle counts. As we have seen in Sec. 4, 𝑁 ppc ∼ 100 particles per cell is a standard requirement for sufficient interpolation accuracy. With a grid-cell count of just 100 3 cells, this means that even at such moderate grid resolutions, we very quickly end up with an enormous number of resolution elements 4 . Therefore, it is important to investigate and develop numerical methods that will help reduce  \n4 For example, if we wanted to do a Hybrid PIC turbulence simulation with 1000 3 grid cells and 𝑁 ppc = 100, we would have 1000 3 × 100 = 10 11 resolution elements to compute on. This is not only significant in terms of memory requirements, but the compute time per time step is also significantly higher for particle operations compared to MHD grid-cell operations.  \nFigure 17. SameasFig.16,but for supersonic turbulence with a target M = 2. As for the subsonic case, a cooling time of ( Δ 𝑡 ) cool ∼ 0 . 1 𝑡 cool, is required to maintain 𝑣 th at its target level and therefore maintain the target Mach number. We also perform a test with ( Δ 𝑡 ) cool = 10 -3 𝑡 cool, to demonstrate the convergence of the thermal speed and the Mach number as the cooling frequency is increased in the supersonic regime and tabulate the results in Tab. 1.  \n<!-- image -->  \n0  \nthe computational resources required to perform these simulations, without compromising the accuracy of the numerical solutions.  \nOne such approach for optimisation is to store and compute solutions on single precision (4 byte) floating-point numbers instead of the standard double precision (8 byte per floating-point number). This reduces the memory footprint of an application by a factor of 2. The time for computations and parallel message passing interface (MPI) communications on single precision numbers is also approximately halved compared to double precision, as long as communication times are dominated by bandwidth and not latency 5 . However, using single precision without safeguards may lead to substantial inaccuracies in the numerical solution 6 . Thus, we will use the 'hybrid  \n5 This depends on the size of the problem and the number of cores used. However, for normal production applications, communication is often bandwidth limited, as communication is only required once per time step, and the data-package size scales with the number of particles, noting that practically all applications use 𝑁 ppc ≫ 1, and thus, the total number of particles determines the relevant problem size rather than the number of grid cells. Packing MPI messages and communicating many quantities per particle means that applications using this Hybrid PIC code are usually bandwidth-dominated and therefore benefit from using hybrid precision compared to double precision.  \nprecision' approach introduced in Federrath et al. (2021). The basic approach is to store and communicate all quantities that occupy large arrays (such as any quantity that is defined per grid cell or per particle) in single precision arrays. However, computations that require double precision for sufficient accuracy are carried out in double precision by explicit promotion to double precision during those operations. The operations that require double precision include summing up the mass, charge and momentum of all the particles. Additionally, the time step is promoted to double precision throughout the code. All global integral quantities for grid variables, such as electric and magnetic energy, as well as particle variables like the componentwise mean particle velocity and the particle Larmor radius, are also calculated in double precision. As shown in Federrath et al. (2021), this approach cuts down memory storage and compute time requirements by nearly a factor of ∼ 2 overall, while retaining sufficient (near double precision) accuracy for all relevant quantities.  \nHere we implement and test this hybrid precision approach for Hybrid PIC, and compare the results of this with pure single precision and pure double precision solutions. To do this, we use the same basic turbulence setup as in the previous two sections, as these represent complex and relevant use cases of the code. Fig. 18 shows the fluctuations in the average particle velocity over time, as a function of the total number of particles, 𝑁 , with single and hybrid precision. The fluctuations in the average particle velocity, 𝜎 (⟨ 𝑣 p ⟩) , are calculated as the standard deviation of the average particle speed, ⟨ 𝑣 p ⟩ , over time ( ≈ 0 -1 𝑡 0 ). The average particle speed and its fluctuations are expected to be zero to within machine precision, as there are no bulk flows in the plasma for the turbulent driving setup. Using single precision, we find that 𝜎 (⟨ 𝑣 p ⟩)/ 10 -6 increases strongly as the number of particles is increased, where 10 -6 is the expected precision level for single precision. This makes using single precision undesirable, as with a higher number of particles the numerical solutions become progressively worse and unreliable. On the other hand, using hybrid precision, we find 𝜎 (⟨ 𝑣 p ⟩)/ 10 -6 ≲ 1, and this value does not change as the number of particles is increased. Thus, the absolute error with hybrid precision remains bound to the level of single precision, while with pure single precision, the errors increase very quickly.  \nWenote that while Fig. 18 shows a specific metric for comparison, other metrics show a similar behaviour. Thus, hybrid precision provides a stable and accurate (to within ≲ 10 -6 ) solution, while single precision does not. We also note that pure double precision outperforms hybrid precision in this metric, as expected (the fluctuations in the average particle speed are zero to within double precision, as expected). However, the main goal here is to establish that in the hybrid-precision approach (while naturally only providing an absolute accuracy between that of single and double precision), the accuracy of the solution is stable with hybrid precision, while it quickly deteriorates with pure single precision.  \nTo further demonstrate that the relevant physical solution is practically as accurate with hybrid precision as it is with pure double precision, we show the Mach number and the magnetic-to-kinetic energy ratio in tests with double and hybrid precision in Fig. 19. We find excellent agreement between the double and hybrid precision runs. We also checked other quantities, including probability density functions and power spectra, which all show excellent agreement as well (see Appendix C). The following subsection compares the computational cost of hybrid vs. double precision, and quantifies the parallel scaling of the code.  \nprecision - the question is whether the quantity of interest is sufficiently accurate in the context of a specific application.  \nFigure 18. The time- and particle-averaged variations in the average particle velocity (averaged over all particles and across the 𝑥 , 𝑦 and 𝑧 components), normalised to the expected precision level with single precision ≈ 10 -6 , as a function of the total number of particles, 𝑁 , for numerical tests with single precision (red) and hybrid precision (blue). The error bars are calculated as the standard deviation of the 𝑥 , 𝑦 and 𝑧 components of 𝜎 (⟨ 𝑣 p ⟩) . Using single precision, the variation in the average particle speed increases, becoming worse with increasing number of particles. By contrast, with hybrid precision, the error is bounded, as desired.  \n<!-- image -->  \nFigure 19. Same as Fig. 11, but for simulations comparing double and hybrid precision. The numerical solutions for the Mach number, M , and the ratio of magnetic energy to kinetic energy, 𝐸 mag / 𝐸 kin, show excellent agreement.  \n<!-- image -->  \n0  \nFigure 20. Time taken by the code to perform one-time step (averaged over a total of 100 time steps), normalised by the total number of particles (top x-axis) as a function of the number of compute cores for the turbulence simulations described in Sec. 6.1, with double precision (magenta) and hybrid precision (blue). The respective dashed lines represent fits with a constant (ideal weak scaling). Both double and hybrid precision cases show excellent parallel scaling up to ∼ 10 4 compute cores. Moreover, hybrid precision simulations are a factor of ∼ 1 . 5 faster than pure double precision calculations. The error bars are obtained from the standard deviation of sub-intervals of 10 time steps over the total time-step interval.  \n<!-- image -->", '6.2 Parallel code scaling': 'Here we test how the performance of our code scales with an increasing number of computational cores for numerical tests with hybrid and double precision, discussed in Sec. 6.1. We show the time per particle per time step as a function of the number of compute cores in Fig. 20. We note that this is a weak scaling test, i.e., the number of particles grows proportionally with the number of compute cores (see top x-axis for the absolute number of particles). Therefore, ideal scaling is indicated by a constant time per particle per time step. We calculate the time by taking an average of over 100 time steps, where the error bars have been estimated as the standard deviation over intervals of 10 time steps in the total time step range. The dashed lines show fits with a constant, i.e., ideal scaling.  \nWe find excellent (near ideal) parallel scaling of the code with both double and hybrid precision up to ∼ 10 4 compute cores. Fig. 20 also shows that the hybrid precision simulations are a factor of ≈ 1 . 5 faster than the double precision calculations. The theoretical speed-up could be as high as factor 2, as discussed in Sec. 6.1 7 , but significant care must be taken in promoting relevant single-precision operations to double precision, in order to maintain sufficient overall accuracy of the code (c.f., Fig. 19) Therefore, a reduction of the required compute time by ∼ 33% with hybrid precision compared to pure double precision is highly beneficial for computationallyexpensive Hybrid PIC simulations.', '7 SUMMARY AND CONCLUSIONS': "We introduced AHKASH - a new Hybrid particle-in-cell (PIC) code within the FLASH astrophysics code framework (Fryxell et al. 2000; Dubey et al. 2008a,b) to study weakly collisional plasmas. The Hybrid PIC code includes state-of-the-art numerical methods for time integration of particle trajectories, i.e., the Boris integrator (Boris 1970; Zenitani & Umeda 2018, Sec. 2.3), and the predictor-predictorcorrector integration algorithm (Kunz et al. 2014b, Sec. 2.7). We use the constrained transport (CT) method to ensure that magnetic fields remain divergence-free (Yee 1966, Sec. 2.6). The code supports various grid → particle and particle → grid interpolation schemes, as well as post-interpolation smoothing to reduce finite particle noise (see Sec. 2.4). It further supports the ' 𝛿 𝑓 method' to study instabilities in weakly collisional plasmas (see Sec. 2.5). The code implements a series of time step constraints to ensure that physical plasma time scales are appropriately resolved (see Sec. 2.12).  \nWeperform several tests to demonstrate the ability and accuracy of the new code to model standard physical problems, such as the motion of a charged particle in the presence of magnetic fields (Sec. 3.1), the propagation of Alfvén and whistler waves (Sec. 3.2), and Landau damping of ion acoustic waves (Sec. 3.3).  \nIn Sec. 4 we compare the quality of various particle → grid interpolation schemes. We find that global quantities are largely independent of the details of the interpolation scheme, provided that at least 1 smoothing pass is performed after the particle deposition step. However, structural details can vary significantly with the interpolation scheme and with the number of particles per cell ( 𝑁 ppc). We suggest that a reasonable compromise between computational feasibility and accuracy is achieved with the cloud-in-cell interpolation and 2 post-deposition smoothing passes for 𝑁 ppc ∼ 100 particles per cell, noting that detailed requirements are case-specific and problemdependent.  \nThe new code further supports turbulence driving (Federrath et al. 2010, 2021), modelled with the Ornstein-Uhlenbeck process (Sec. 2.10). To study steady-state turbulence in weakly collisional plasma, we have to maintain isothermal conditions across the computational domain and to do so we introduced a novel cooling method in the code (Sec. 2.11). This allows the modelling of fully developed, steady-state turbulence, including turbulent dynamos in weakly collisional plasmas, whose properties are sensitive to the Mach number of the plasma. In particular, we show that the cooling method is essential for modelling steady-state supersonic turbulence in Hybrid PIC. We demonstrate and quantify the ability of the cooling method to keep the plasma both locally and globally isothermal (Sec. 5).  \nSince Hybrid PIC simulations can be computationally demanding, we extend the hybrid precision method introduced in Federrath et al. (2021) for grid-based hydrodynamics to our Hybrid PIC code (Sec. 6.1). This method uses a hybrid approach that promotes critical computations to double precision arithmetic (8 bytes per floatingpoint number), but otherwise stores and communicates all quantities in single precision (4 bytes per floating-point number). This provides a factor of ∼ 1 . 5 speed-up over a pure double precision calculation (c.f. Fig 20), but retains sufficient accuracy and precision throughout the calculations (c.f. Figs.18 and 19). Finally, Fig. 20 shows that the new code exhibits excellent parallel scalability up to 10,000 compute cores.  \nWeplan to use AHKASH to study astrophysical problems, especially the physics of the turbulent dynamo in supersonic plasmas and its application to the intracluster medium of galaxy clusters.", 'ACKNOWLEDGEMENTS': 'We thank Matthew W. Kunz for sharing his expertise on developing Hybrid PIC codes and for many useful discussions that have been very beneficial to this work. R. A. C. acknowledges that this work was supported by an NCI HPC-AI Talent Program 2023 Scholarship (project gp08), with computational resources provided by NCI Australia, an NCRIS-enabled capability supported by the Australian Government. C. F. acknowledges funding provided by the Australian Research Council (Discovery Project DP230102280), and the AustraliaGermany Joint Research Cooperation Scheme (UA-DAAD). We further acknowledge high-performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391), the Australian National Computational Infrastructure (grant ek9) and the Pawsey Supercomputing Centre (project pawsey0810) in the framework of the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme. The simulation software, FLASH , was in part developed by the Flash Centre for Computational Science at the University of Chicago and the Department of Physics and Astronomy at the University of Rochester.', 'DATA AVAILABILITY': 'The Hybrid particle-in-cell code AHKASH has been implemented in our private fork of the FLASH code. The code and results from our simulations will be shared on reasonable request to the corresponding author.', 'REFERENCES': 'Ruzmaikin A., Shukurov A., Sokoloff D., 1988, Magnetic Fields of Galaxies. Academic Press, Dordrecht  \nSchekochihin A. A., Cowley S. C., 2006, Physics of Plasmas, 13, 056501 Schekochihin A. A., Cowley S. C., Taylor S. F., Maron J. L., McWilliams J. C., 2004, ApJ, 612, 276  \nSchober J., Schleicher D., Federrath C., Klessen R., Banerjee R., 2012, Phys. Rev. E, 85, 026303  \nSeta A., Federrath C., 2020, MNRAS, 499, 2076  \nSeta A., Federrath C., 2021, Physical Review Fluids, 6, 103701  \nSeta A., Federrath C., 2022, MNRAS, 514, 957  \nSeta A., Bushby P. J., Shukurov A., Wood T. S., 2020, Physical Review Fluids, 5, 043702  \nShalaby M., Broderick A. E., Chang P., Pfrommer C., Lamberts A., Puchwein E., 2017, ApJ, 841, 52  \nSimionescu A., et al., 2019, Space Sci. Rev., 215, 24  \nSpitkovsky A., 2005, in Bulik T., Rudak B., Madejski G., eds, American Institute of Physics Conference Series Vol. 801, Astrophysical Sources of High Energy Particles and Radiation. pp 345-350 ( arXiv:astro-ph/0603211 ), doi:10.1063/1.2141897  \nSquire J., Kunz M. W., Arzamasskiy L., Johnston Z., Quataert E., Schekochihin A. A., 2023a, Journal of Plasma Physics, 89, 905890417  \nSquire J., Meyrand R., Kunz M. W., 2023b, ApJ, 957, L30  \nSt-Onge D. A., Kunz M. W., 2018, ApJ, 863, L25  \nStone J. M., Gardiner T. A., Teuben P., Hawley J. F., Simon J. B., 2008, ApJS, 178, 137  \nStone J. M., Tomida K., White C. J., Felker K. G., 2020, ApJS, 249, 4  \nVerscharen D., Klein K. G., Maruca B. A., 2019, Living Reviews in Solar  \nPhysics, 16, 5  \nWinske D., Karimabadi H., Le A., Omidi N., Roytershteyn V., Stanier A., 2022, arXiv e-prints, p. arXiv:2204.01676  \nYee K., 1966, IEEE Transactions on Antennas and Propagation, 14, 302  \nZenitani S., Umeda T., 2018, Physics of Plasmas, 25, 112110  \nZinger E., Dekel A., Birnboim Y., Nagai D., Lau E., Kravtsov A. V., 2018, MNRAS, 476, 56', 'APPENDIX A: WEIGHT FUNCTIONS': 'The Cartesian components of r and r 𝑙 are ( 𝑥, 𝑦, 𝑧 ) and ( 𝑥 𝑙 , 𝑦 𝑙 , 𝑧 𝑙 ) , respectively. Using these, we define h = ( ℎ 𝑥 , ℎ 𝑦 , ℎ 𝑧 ) , where ℎ 𝑥 = | 𝑥 -𝑥 𝑙 |/ Δ 𝑥 , ℎ 𝑦 = | 𝑦 -𝑦 𝑙 |/ Δ 𝑦 and ℎ 𝑧 = | 𝑧 -𝑧 𝑙 |/ Δ 𝑧 . Δ 𝑥 , Δ 𝑦 and Δ 𝑧 are the grid-cell size in the 𝑥 , 𝑦 and 𝑧 directions, respectively.  \nThe nearest-grid-point weight function, where the particle is assigned completely to its nearest grid cell, is defined as  \n𝑊 NGP ( r -r 𝑙 ) = ( 1 , if | r -r 𝑙 | = min (| h |) 0 , otherwise . (A1)  \nThe distance between the particle and each grid-cell centre is calculated and the minimum value, min (| h |) , gives us the grid-cell centre closest to the particle.  \nIn the cloud-in-cell (CIC) interpolation kernel, a \'cloud" is created by each particle in the shape of a grid cell and the weight function can be described as  \n𝑊 CIC ( r -r 𝑙 ) = ( ( 1 -ℎ 𝑥 )( 1 -ℎ 𝑦 )( 1 -ℎ 𝑧 ) , if ℎ 𝑥 , ℎ 𝑦 , ℎ 𝑧 < 1 0 , otherwise .  \n(A2)  \nThe triangular-shaped-cloud weight function is defined co-ordinate wise as  \n𝑊 TSC ( r 𝑖 - ( r 𝑙 ) 𝑖 ) =         1 -2 ℎ 2 𝑖 , if ℎ 𝑖 ≤ 0 . 5 2 ( 1 -ℎ 𝑖 ) 2 , if 0 . 5 < ℎ 𝑖 < 1 0 , otherwise . (A3)  \nThe total weight function, 𝑊 TSC ( r -r 𝑙 ) , is then given by  \nIn both CIC and TSC interpolation kernels, the size of the weight function is one grid cell ( ∼ Δ 𝑥 ). One can also construct larger stencils for the weight function, however, the computational cost of the interpolation operation increases as the extent of the weight function increases.  \n𝑊 TSC ( r -r 𝑙 ) = 𝑊 TSC ( x -x 𝑙 ) 𝑊 TSC ( y -y 𝑙 ) 𝑊 TSC ( z -z 𝑙 ) . (A4)', 'APPENDIX B: CELL-CENTERED FINITE DIFFERENCE METHOD': "Wealsohaveacell-centered approach to update magnetic fields using Faraday's law, Eq. 13. In this method, the electric and magnetic fields are defined and evolved on the cell-centers. This numerical stencil is depicted in Fig. B1.  \nIn this method, the magnetic field is updated in time in the following way  \n𝐵 𝑡 + Δ 𝑡 𝑥 ( 𝑖, 𝑗, 𝑘 ) = 𝐵 𝑡 𝑥 ( 𝑖, 𝑗, 𝑘 ) + Δ 𝑡 2 Δ 𝑧 GLYPH<2> e 𝐸 𝑡 + Δ 𝑡 / 2 𝑦 ( 𝑖, 𝑗, 𝑘 + 1 ) -e 𝐸 𝑡 + Δ 𝑡 / 2 𝑦 ( 𝑖, 𝑗, 𝑘 -1 ) GLYPH<3> -Δ 𝑡 2 Δ 𝑦 GLYPH<2> e 𝐸 𝑡 + Δ 𝑡 / 2 𝑧 ( 𝑖, 𝑗 + 1 , 𝑘 ) -e 𝐸 𝑡 + Δ 𝑡 / 2 𝑧 ( 𝑖, 𝑗 -1 , 𝑘 ) GLYPH<3> , (B1)  \n𝐵  \n+  \n𝑡  \n𝑦  \n𝑖, 𝑗, 𝑘  \n=  \n𝐵  \n𝑡  \n𝑦  \n𝑖, 𝑗, 𝑘  \n+ Δ 𝑡 2 Δ 𝑥 GLYPH<2> e 𝐸 𝑡 + Δ 𝑡 / 2 𝑧 ( 𝑖 + 1 , 𝑗 , 𝑘 ) -e 𝐸 𝑡 + Δ 𝑡 / 2 𝑧 ( 𝑖 -1 , 𝑗 , 𝑘 ) GLYPH<3> Δ 𝑡 GLYPH<2> e 𝐸 𝑡 + Δ 𝑡 / 2 𝑥 𝑖, 𝑗, 𝑘 1 e 𝐸 𝑡 + Δ 𝑡 / 2 𝑥 𝑖, 𝑗, 𝑘 1 GLYPH<3> , (B2)  \n(  \n)  \n(  \n-2 Δ 𝑧 ( + ) -( -)  \n𝐵 𝑡 + Δ 𝑡 𝑧 ( 𝑖, 𝑗, 𝑘 ) = 𝐵 𝑡 𝑧 ( 𝑖, 𝑗, 𝑘 ) + Δ 𝑡 2 Δ 𝑦 GLYPH<2> e 𝐸 𝑡 + Δ 𝑡 / 2 𝑥 ( 𝑖, 𝑗 + 1 , 𝑘 ) -e 𝐸 𝑡 + Δ 𝑡 / 2 𝑥 ( 𝑖, 𝑗 -1 , 𝑘 ) GLYPH<3> -Δ 𝑡 2 Δ 𝑥 GLYPH<2> e 𝐸 𝑡 + Δ 𝑡 / 2 𝑦 ( 𝑖 + 1 , 𝑗 , 𝑘 ) -e 𝐸 𝑡 + Δ 𝑡 / 2 𝑦 ( 𝑖 -1 , 𝑗 , 𝑘 ) GLYPH<3> . (B3)  \nThe divergence of the magnetic field in the cell-centered finite difference method can be written as  \n∇ · B 𝑡 + Δ 𝑡 ( 𝑖, 𝑗, 𝑘 ) = 1 2 Δ 𝑥 GLYPH<2> 𝐵 𝑡 + Δ 𝑡 𝑥 ( 𝑖 + 1 , 𝑗 , 𝑘 ) -𝐵 𝑡 + Δ 𝑡 𝑥 ( 𝑖 -1 , 𝑗 , 𝑘 ) GLYPH<3> + 1 2 Δ 𝑦 GLYPH<2> 𝐵 𝑡 + Δ 𝑡 𝑦 ( 𝑖, 𝑗 + 1 , 𝑘 ) -𝐵 𝑡 + Δ 𝑡 𝑦 ( 𝑖, 𝑗 -1 , 𝑘 ) GLYPH<3> + 1 2 Δ 𝑧 GLYPH<2> 𝐵 𝑡 + Δ 𝑡 𝑧 ( 𝑖, 𝑗, 𝑘 + 1 ) -𝐵 𝑡 + Δ 𝑡 𝑧 ( 𝑖, 𝑗, 𝑘 -1 ) GLYPH<3> . (B4)  \nUsing Eq. B1 - Eq. B3 in Eq. B4, we obtain that ∇ · B 𝑡 + Δ 𝑡 ( 𝑖, 𝑗, 𝑘 ) = 0. This is ensured by the construction of the cell-centered finite difference method similar to the constrained transport method (see Sec. 2.6), albeit on a stencil twice the size. The above calculations show that as long as the construction of a numerical method ensures that the analytical expression, ∇·∇× C = 0, is satisfied for any vector C , it can be used to ensure that magnetic fields are divergence-free.", 'APPENDIX C: HYBRID PRECISION - MAGNETIC SPECTRA AND PROBABILITY DENSITY FUNCTIONS': "Probability density functions (PDFs) and power spectra are key quantities in the study of turbulence and magnetic field amplification  \nΔ  \n𝑡  \n)  \ny  \nFigure B1. Cell-centered finite difference method to solve Faraday's law, Eq. 13. The charge density and ion current are deposited on the grid-cell center, depicted by the yellow circle. The electric and magnetic fields are also calculated on the cell-centers. The magnetic field of the central grid cell is updated using the electric field of the neighbouring grid cells as described by Eq. B1 - Eq. B3.  \n<!-- image -->", '(Federrath 2013; Federrath & Banerjee 2015; Seta & Federrath 2020,': '2021). In this section, we compare the PDFs and power spectra obtained in the hybrid vs. double precision tests described in Sec. 6.1, with Fig. C1 showing the results. We find excellent agreement between the double and hybrid precision tests for all the examined quantities, further demonstrating that the numerical solutions with hybrid precision are in practice as accurate as those with double precision.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \nFigure C1. Probability density functions (PDFs) of the natural logarithm of the mass density normalised to the mean mass density, ln ( 𝜌 m /⟨ 𝜌 m ⟩) (top left), the 𝑥 -component of the velocity normalised to the root mean square velocity, 𝑢 𝑥 / 𝑢 rms (top right), the 𝑥 -component of the magnetic field normalised to the root mean square magnetic field, 𝐵 𝑥 / 𝐵 rms (bottom left), and the power spectra of magnetic energy (bottom right), for the hybrid vs. double precision tests discussed in Sec. 6.1. The PDFs and power spectra are time-averaged in the kinematic regime of the turbulent dynamo ( ≈ 3 -15 𝑡 0). We find excellent agreement between the double and hybrid precision statistics.  \n<!-- image -->'}