We present a detailed guide to advanced collisionless fluid models that incorporate kinetic effects into the fluid framework, and that are much closer to the collisionless kinetic description than traditional magnetohydrodynamics. Such fluid models are directly applicable to modelling the turbulent evolution of a vast array of astrophysical plasmas, such as the solar corona and the solar wind, the interstellar medium, as well as accretion disks and galaxy clusters. The text can be viewed as a detailed guide to Landau fluid models and it is divided into two parts. Part 1 is dedicated to fluid models that are obtained by closing the fluid hierarchy with simple (non-Landau fluid) closures. Part 2 is dedicated to Landau fluid closures. Here in Part 1, we discuss the fluid model of Chew–Goldberger–Low (CGL) in great detail, together with fluid models that contain dispersive effects introduced by the Hall term and by the finite Larmor radius corrections to the pressure tensor. We consider dispersive effects introduced by the non-gyrotropic heat flux vectors. We investigate the parallel and oblique firehose instability, and show that the non-gyrotropic heat flux strongly influences the maximum growth rate of these instabilities. Furthermore, we discuss fluid models that contain evolution equations for the gyrotropic heat flux fluctuations and that are closed at the fourth-moment level by prescribing a specific form for the distribution function. For the bi-Maxwellian distribution, such a closure is known as the ‘normal’ closure. We also discuss a fluid closure for the bi-kappa distribution. Finally, by considering one-dimensional Maxwellian fluid closures at higher-order moments, we show that such fluid models are always unstable. The last possible non Landau fluid closure is therefore the ‘normal’ closure, and beyond the fourth-order moment, Landau fluid closures are required.

In Part 2 of our guide to collisionless fluid models, we concentrate on Landau fluid closures. These closures were pioneered by Hammett and Perkins and allow for the rigorous incorporation of collisionless Landau damping into a fluid framework. It is Landau damping that sharply separates traditional fluid models and collisionless kinetic theory, and is the main reason why the usual fluid models do not converge to the kinetic description, even in the long-wavelength low-frequency limit. We start with a brief introduction to kinetic theory, where we discuss in detail the plasma dispersion function $Z(\unicode[STIX]{x1D701})$, and the associated plasma response function $R(\unicode[STIX]{x1D701})=1+\unicode[STIX]{x1D701}Z(\unicode[STIX]{x1D701})=-Z^{\prime }(\unicode[STIX]{x1D701})/2$. We then consider a one-dimensional (1-D) (electrostatic) geometry and make a significant effort to map all possible Landau fluid closures that can be constructed at the fourth-order moment level. These closures for parallel moments have general validity from the largest astrophysical scales down to the Debye length, and we verify their validity by considering examples of the (proton and electron) Landau damping of the ion-acoustic mode, and the electron Landau damping of the Langmuir mode. We proceed by considering 1-D closures at higher-order moments than the fourth order, and as was concluded in Part 1, this is not possible without Landau fluid closures. We show that it is possible to reproduce linear Landau damping in the fluid framework to any desired precision, thus showing the convergence of the fluid and collisionless kinetic descriptions. We then consider a 3-D (electromagnetic) geometry in the gyrotropic (long-wavelength low-frequency) limit and map all closures that are available at the fourth-order moment level. In appendix A, we provide comprehensive tables with Padé approximants of $R(\unicode[STIX]{x1D701})$ up to the eighth-pole order, with many given in an analytic form.

Integrated modelling of electron runaway requires computationally expensive kinetic models that are self-consistently coupled to the evolution of the background plasma parameters. The computational expense can be reduced by using parameterized runaway generation rates rather than solving the full kinetic problem. However, currently available generation rates neglect several important effects; in particular, they are not valid in the presence of partially ionized impurities. In this work, we construct a multilayer neural network for the Dreicer runaway generation rate which is trained on data obtained from kinetic simulations performed for a wide range of plasma parameters and impurities. The neural network accurately reproduces the Dreicer runaway generation rate obtained by the kinetic solver. By implementing it in a fluid runaway-electron modelling tool, we show that the improved generation rates lead to significant differences in the self-consistent runaway dynamics as compared to the results using the previously available formulas for the runaway generation rate.

The condition of omnigenity is investigated, and applied to the near-axis expansion of Garren & Boozer (Phys. Fluids B, vol. 3 (10), 1991a, pp. 2805–2821). Due in part to the particular analyticity requirements of the near-axis expansion, we find that, excluding quasi-symmetric solutions, only one type of omnigenity, namely quasi-isodynamicity, can be satisfied at first order in the distance from the magnetic axis. Our construction provides a parameterization of the space of such solutions, and the cylindrical reformulation and numerical method of Landreman & Sengupta (J. Plasma Phys., vol. 84 (6), 2018, 905840616); Landreman et al. (J. Plasma Phys., vol. 85 (1), 2019, 905850103), enables their efficient numerical construction.

Dipole and stellarator geometries are capable of confining plasmas of arbitrary neutrality, ranging from pure electron plasmas through to quasineutral. The diocotron mode is known to be important in non-neutral plasmas and has been widely studied. However, drift mode dynamics, dominating quasineutral plasmas, has received very little by way of attention in the non-neutral context. Here, we show that non-neutral plasmas can be unstable respect to both density-gradient- and temperature-gradient-driven instabilities. A local shearless slab limit is considered for simplicity. A key feature of non-neutral plasmas is the development of strong electric fields, in this local limit of straight field line geometry, the effect of the corresponding $\boldsymbol{E}\times \boldsymbol{B}$ drift is limited to the Doppler shift of the complex frequency $\unicode[STIX]{x1D714}\rightarrow \unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}_{E}$. However, the breaking of the quasineutrality condition still leads to interesting dynamics in non-neutral plasmas. In this paper we address the behaviour of a number of gyrokinetic modes in electron–ion and electron–positron plasmas with arbitrary degree of neutrality.

We present three-dimensional direct numerical simulations and an analytic model of reflection-driven magnetohydrodynamic (MHD) turbulence in the solar wind. Our simulations describe transverse, non-compressive MHD fluctuations within a narrow magnetic flux tube that extends from the photosphere, through the chromosphere and corona and out to a heliocentric distance $r$ of 21 solar radii $(R_{\odot })$. We launch outward-propagating ‘$\boldsymbol{z}^{+}$ fluctuations’ into the simulation domain by imposing a randomly evolving photospheric velocity field. As these fluctuations propagate away from the Sun, they undergo partial reflection, producing inward-propagating ‘$\boldsymbol{z}^{-}$ fluctuations’. Counter-propagating fluctuations subsequently interact, causing fluctuation energy to cascade to small scales and dissipate. Our analytic model incorporates dynamic alignment, allows for strongly or weakly turbulent nonlinear interactions and divides the $\boldsymbol{z}^{+}$ fluctuations into two populations with different characteristic radial correlation lengths. The inertial-range power spectra of $\boldsymbol{z}^{+}$ and $\boldsymbol{z}^{-}$ fluctuations in our simulations evolve toward a $k_{\bot }^{-3/2}$ scaling at $r>10R_{\odot }$, where $k_{\bot }$ is the wave-vector component perpendicular to the background magnetic field. In two of our simulations, the $\boldsymbol{z}^{+}$ power spectra are much flatter between the coronal base and $r\simeq 4R_{\odot }$. We argue that these spectral scalings are caused by: (i) high-pass filtering in the upper chromosphere; (ii) the anomalous coherence of inertial-range $\boldsymbol{z}^{-}$ fluctuations in a reference frame propagating outwards with the $\boldsymbol{z}^{+}$ fluctuations; and (iii) the change in the sign of the radial derivative of the Alfvén speed at $r=r_{\text{m}}\simeq 1.7R_{\odot }$, which disrupts this anomalous coherence between $r=r_{\text{m}}$ and $r\simeq 2r_{\text{m}}$. At $r>1.3R_{\odot }$, the turbulent heating rate in our simulations is comparable to the turbulent heating rate in a previously developed solar-wind model that agreed with a number of observational constraints, consistent with the hypothesis that MHD turbulence accounts for much of the heating of the fast solar wind.

We show that when a solid plasma foil with a density gradient on the front surface is irradiated by an intense laser pulse at a grazing angle, ${\sim}80^{\circ }$, a relativistic electron vortex is excited in the near-critical-density layer after the laser pulse depletion. The vortex structure and dynamics are studied using particle-in-cell simulations. Due to the asymmetry introduced by non-uniform background density, the vortex drifts at a constant velocity, typically $0.2{-}0.3$ times the speed of light. The strong magnetic field inside the vortex leads to significant charge separation; in the corresponding electric field initially stationary protons can be captured and accelerated to twice the velocity of the vortex (100–200 MeV). A representative scenario – with laser intensity of $10^{21}~\text{W}~\text{cm}^{-2}$ – is discussed: two-dimensional simulations suggest that a quasi-monoenergetic proton beam can be obtained with a mean energy 140 MeV and an energy spread of ${\sim}10\,\%$. We derive an analytical estimate for the vortex velocity in terms of laser and plasma parameters, demonstrating that the maximum proton energy can be controlled by the incidence angle of the laser and the plasma density gradient.

I analyse and numerically evaluate the radiation field generated by an experimentally realized embodiment of an electric polarization current whose rotating distribution pattern moves with linear speeds exceeding the speed of light in vacuum. I find that the flux density of the resulting emission (i) has a dominant value and is linearly polarized within a sharply delineated radiation beam whose orientation and polar width are determined by the range of values of the linear speeds of the rotating source distribution, and (ii) decays with the distance $d$ from the source as $d^{-\unicode[STIX]{x1D6FC}}$ in which the value of $\unicode[STIX]{x1D6FC}$ lies between $1$ and $2$ (instead of being equal to $2$ as in a conventional radiation) across the beam. In that the rate at which boundaries of the retarded distribution of such a source change with time depends on its duration monotonically, this is an intrinsically transient emission process: temporal rate of change of the energy density of the radiation generated by it has a time-averaged value that is negative (instead of being zero as in a conventional radiation) at points where the envelopes of the wave fronts emanating from the constituent volume elements of the source distribution are cusped. The difference in the fluxes of power across any two spheres centred on the source is in this case balanced by the change with time of the energy contained inside the shell bounded by those spheres. These results are relevant not only to long-range transmitters in communications technology but also to astrophysical objects containing rapidly rotating neutron stars (such as pulsars) and to the interpretation of the energetics of the multi-wavelength emissions from sources that lie at cosmological distances (such as radio and gamma-ray bursts). The analysis presented in this paper is self-contained and supersedes my earlier works on this problem.

The $\unicode[STIX]{x1D6FF}f$-PIC method is widely used for electrostatic particle-in-cell (PIC) simulations. Its basic idea is the ansatz $f=f_{0}+\unicode[STIX]{x1D6FF}f$ ($\unicode[STIX]{x1D6FF}f$-ansatz) where the particle distribution function $f$ is split into a usually time-independent background $f_{0}$ and a time-dependent perturbation $\unicode[STIX]{x1D6FF}f$. In principle, it can also be used for electromagnetic gyrokinetic PIC simulations, but the required number of markers can be so large that PIC simulations become impractical. The reason is a decreasing efficiency of the $\unicode[STIX]{x1D6FF}f$-ansatz for the so-called ‘Hamiltonian formulation’ using $p_{\Vert }$ as a dynamic variable. As a result, the density and current moment of the distribution function develop large statistical errors. To overcome this obstacle we propose to solve the potential equations in the symplectic formulation using $v_{\Vert }$ as a dynamic variable. The distribution function itself is still evolved in the Hamiltonian formulation which is better suited for the numerical integration of the parallel dynamics. The contributions from the full Jacobian of phase space, a finite velocity sphere of the simulation domain and a shifted Maxwellian as a background are considered. Special care has been taken at the discretisation level to make damped magnetohydrodynamics (MHD) mode simulations within a realistic gyrokinetic model feasible. This includes devices like e.g. large tokamaks with a small aspect ratio.

In strictly collisionless electrostatic shocks, the ion distribution function can develop discontinuities along phase-space separatrices, due to partial reflection of the ion population. In this paper, we depart from the strictly collisionless regime and present a semi-analytical model for weakly collisional kinetic shocks. The model is used to study the effect of small but finite collisionalities on electrostatic shocks, and they are found to smooth out these discontinuities into growing boundary layers. More importantly, ions diffuse into and accumulate in the previously empty regions of phase space, and, by upsetting the charge balance, lead to growing downstream oscillations of the electrostatic potential. We find that the collisional age of the shock is the more relevant measure of the collisional effects than the collisionality, where the former can become significant during the lifetime of the shock, even for weak collisionalities.

We study the resistive evolution of a localized self-organizing magnetohydrodynamic equilibrium. In this configuration the magnetic forces are balanced by a pressure force caused by a toroidal depression in the pressure. Equilibrium is attained when this low-pressure region prevents further expansion into the higher-pressure external plasma. We find that, for the parameters investigated, the resistive evolution of the structures follows a universal pattern when rescaled to resistive time. The finite resistivity causes both a decrease in the magnetic field strength and a finite slip of the plasma fluid against the static equilibrium. This slip is caused by a Pfirsch–Schlüter-type diffusion, similar to what is seen in tokamak equilibria. The net effect is that the configuration remains in magnetostatic equilibrium whilst it slowly grows in size. The rotational transform of the structure becomes nearly constant throughout the entire structure, and decreases according to a power law. In simulations this equilibrium is observed when highly tangled field lines relax in a high-pressure (relative to the magnetic field strength) environment, a situation that occurs when the twisted field of a coronal loop is ejected into the interplanetary solar wind. In this paper we relate this localized magnetohydrodynamic equilibrium to magnetic clouds in the solar wind.

Stochastic heating refers to an increase in the average magnetic moment of charged particles interacting with electromagnetic fluctuations whose frequencies are smaller than the particles’ cyclotron frequencies. This type of heating arises when the amplitude of the gyroscale fluctuations exceeds a certain threshold, causing particle orbits in the plane perpendicular to the magnetic field to become stochastic rather than nearly periodic. We consider the stochastic heating of protons by Alfvén-wave (AW) and kinetic-Alfvén-wave (KAW) turbulence, which may make an important contribution to the heating of the solar wind. Using phenomenological arguments, we derive the stochastic-proton-heating rate in plasmas in which $\unicode[STIX]{x1D6FD}_{\text{p}}\sim 1$–30, where $\unicode[STIX]{x1D6FD}_{\text{p}}$ is the ratio of the proton pressure to the magnetic pressure. (We do not consider the $\unicode[STIX]{x1D6FD}_{\text{p}}\gtrsim 30$ regime, in which KAWs at the proton gyroscale become non-propagating.) We test our formula for the stochastic-heating rate by numerically tracking test-particle protons interacting with a spectrum of randomly phased AWs and KAWs. Previous studies have demonstrated that at $\unicode[STIX]{x1D6FD}_{\text{p}}\lesssim 1$, particles are energized primarily by time variations in the electrostatic potential and thermal-proton gyro-orbits are stochasticized primarily by gyroscale fluctuations in the electrostatic potential. In contrast, at $\unicode[STIX]{x1D6FD}_{\text{p}}\gtrsim 1$, particles are energized primarily by the solenoidal component of the electric field and thermal-proton gyro-orbits are stochasticized primarily by gyroscale fluctuations in the magnetic field.

Tokamak start-up is strongly dependent on the state of the initial plasma formed during plasma breakdown. We have investigated through numerical simulations the effects that the pre-filling pressure and induced electric field have on pure ohmic heating during the breakdown process. Three breakdown modes during the discharge are found, as a function of different initial parameters: no breakdown mode, successful breakdown mode and runaway mode. No breakdown mode often occurs with low electric field or high pre-filling pressure, while runaway electrons are usually easy to generate at high electric field or low pre-filling pressure (${<}1.33\times 10^{-4}$  Pa). The plasma behaviours and the physical mechanisms under the three breakdown modes are discussed. We have identified the electric field and pressure values at which the different modes occur. In particular, when the electric field is $0.3~\text{V}~\text{m}^{-1}$ (the value at which ITER operates), the pressure range for possible breakdown becomes narrow, which is consistent with Lloyd’s theoretical prediction. In addition, for $0.3~\text{V}~\text{m}^{-1}$, the optimal pre-filling pressure range obtained from our simulations is $1.33\times 10^{-3}\sim 2.66\times 10^{-3}$  Pa, in good agreement with ITER’s design. Besides, we also find that the Townsend discharge model does not appropriately describe the plasma behaviour during tokamak breakdown due to the presence of a toroidal field. Furthermore, we suggest three possible operation mechanisms for general start-up scenarios which could better control the breakdown phase.

The Arbitrary Linear Plasma Solver (ALPS) is a parallelised numerical code that solves the dispersion relation in a hot (even relativistic) magnetised plasma with an arbitrary number of particle species with arbitrary gyrotropic equilibrium distribution functions for any direction of wave propagation with respect to the background field. ALPS reads the background momentum distributions as tables of values on a $(p_{\bot },p_{\Vert })$ grid, where $p_{\bot }$ and $p_{\Vert }$ are the momentum coordinates in the directions perpendicular and parallel to the background magnetic field, respectively. We present the mathematical and numerical approach used by ALPS and introduce our algorithms for the handling of poles and the analytic continuation for the Landau contour integral. We then show test calculations of dispersion relations for a selection of stable and unstable configurations in Maxwellian, bi-Maxwellian, $\unicode[STIX]{x1D705}$-distributed and Jüttner-distributed plasmas. These tests demonstrate that ALPS derives reliable plasma dispersion relations. ALPS will make it possible to determine the properties of waves and instabilities in the non-equilibrium plasmas that are frequently found in space, laboratory experiments and numerical simulations.

In this paper, weak-turbulence theory is used to investigate the nonlinear evolution of the parametric instability in three-dimensional low-$\unicode[STIX]{x1D6FD}$ plasmas at wavelengths much greater than the ion inertial length under the assumption that slow magnetosonic waves are strongly damped. It is shown analytically that the parametric instability leads to an inverse cascade of Alfvén wave quanta, and several exact solutions to the wave kinetic equations are presented. The main results of the paper concern the parametric decay of Alfvén waves that initially satisfy $e^{+}\gg e^{-}$, where $e^{+}$ and $e^{-}$ are the frequency ($f$) spectra of Alfvén waves propagating in opposite directions along the magnetic field lines. If $e^{+}$ initially has a peak frequency $f_{0}$ (at which $fe^{+}$ is maximized) and an ‘infrared’ scaling $f^{p}$ at smaller $f$ with $-1<p<1$, then $e^{+}$ acquires an $f^{-1}$ scaling throughout a range of frequencies that spreads out in both directions from $f_{0}$. At the same time, $e^{-}$ acquires an $f^{-2}$ scaling within this same frequency range. If the plasma parameters and infrared $e^{+}$ spectrum are chosen to match conditions in the fast solar wind at a heliocentric distance of 0.3 astronomical units (AU), then the nonlinear evolution of the parametric instability leads to an $e^{+}$ spectrum that matches fast-wind measurements from the Helios spacecraft at 0.3 AU, including the observed $f^{-1}$ scaling at $f\gtrsim 3\times 10^{-4}~\text{Hz}$. The results of this paper suggest that the $f^{-1}$ spectrum seen by Helios in the fast solar wind at $f\gtrsim 3\times 10^{-4}~\text{Hz}$ is produced in situ by parametric decay and that the $f^{-1}$ range of $e^{+}$ extends over an increasingly narrow range of frequencies as $r$ decreases below 0.3 AU. This prediction will be tested by measurements from the Parker Solar Probe.

Two signatures of quantum effects on radiation reaction in the collision of a ${\sim}$GeV electron beam with a high intensity (${>}3\times 10^{20}~\text{W}~\text{cm}^{-2}$) laser pulse have been considered. We show that the decrease in the average energy of the electron beam may be used to measure the Gaunt factor $g$ for synchrotron emission. We derive an equation for the evolution of the variance in the energy of the electron beam in the quantum regime, i.e. quantum efficiency parameter $\unicode[STIX]{x1D702}\not \ll 1$. We show that the evolution of the variance may be used as a direct measure of the quantum stochasticity of the radiation reaction and determine the parameter regime where this is observable. For example, stochastic emission results in a 25 % increase in the standard deviation of the energy spectrum of a GeV electron beam, 1 fs after it collides with a laser pulse of intensity $10^{21}~\text{W}~\text{cm}^{-2}$. This effect should therefore be measurable using current high-intensity laser systems.

Although, for current laser pulse energies, the weakly nonlinear regime of laser wakefield acceleration is known to be the optimal for reaching the highest possible electron energies, the capabilities of upcoming large laser systems will provide the possibility of running highly nonlinear regimes of laser pulse propagation in underdense or near-critical plasmas. Using an extended particle-in-cell (PIC) model that takes into account all the relevant physics, we show that such regimes can be implemented with external guiding for a relatively long distance of propagation and allow for the stable transformation of laser energy into other types of energy, including the kinetic energy of a large number of high energy electrons and their incoherent emission of photons. This is despite the fact that the high intensity of the laser pulse triggers a number of new mechanisms of energy depletion, which we investigate systematically.

Systems of multiple current sheets arise in various situations in natural plasmas, such as at the heliospheric current sheet in the solar wind and in the outer heliosphere in the heliosheath. Previous three-dimensional simulations have shown that such systems can develop turbulent-like fluctuations resulting from a forward and inverse cascade in wave vector space. We present a study of the transition to turbulence of such multiple current sheet systems, including the effects of adding a magnetic guide field and velocity shears across the current sheets. Three-dimensional hybrid simulations are performed of systems with eight narrow current sheets in a triply periodic geometry. We carry out a number of different analyses of the evolution of the fluctuations as the initially highly ordered state relaxes to one which resembles turbulence. Despite the evidence of a forward and inverse cascade in the fluctuation power spectra, we find that none of the simulated cases have evidence of intermittency after the initial period of fast reconnection associated with the ion tearing instability at the current sheets. Cancellation analysis confirms that the simulations have not evolved to a state which can be identified as fully developed turbulence. The addition of velocity shears across the current sheets slows the evolution in the properties of the fluctuations, but by the end of the simulation they are broadly similar. However, if the simulation is constrained to be two-dimensional, differences are found, indicating that fully three-dimensional simulations are important when studying the evolution of an ordered equilibrium towards a turbulent-like state.

This paper summarises some of the recent progress that has been made in understanding astrophysical plasma turbulence in the solar wind, from in situ spacecraft observations. At large scales, where the turbulence is predominantly Alfvénic, measurements of critical balance, residual energy and three-dimensional structure are discussed, along with comparison to recent models of strong Alfvénic turbulence. At these scales, a few per cent of the energy is also in compressive fluctuations, and their nature, anisotropy and relation to the Alfvénic component is described. In the small-scale kinetic range, below the ion gyroscale, the turbulence becomes predominantly kinetic Alfvén in nature, and measurements of the spectra, anisotropy and intermittency of this turbulence are discussed with respect to recent cascade models. One of the major remaining questions is how the turbulent energy is dissipated, and some recent work on this question, in addition to future space missions which will help to answer it, are briefly discussed.

In recent papers, we show how accretion disks around massive black holes could act as dynamos producing magnetic jets similar to the jets that create spheromaks in the laboratory. In this paper, we discuss how these magnetic astrophysical jets might naturally produce runaway ion beams accelerated to $10^{20}$  eV or more, finally ejected as ultra high energy cosmic rays (UHECRs) long regarded as one of the mysteries of astrophysics. The acceleration is mainly due to the drift cyclotron loss cone kinetic instability known from plasma research. Experiments and simulations are suggested to verify the acceleration process.