Microphysical Plasma Relations from Special-relativistic Turbulence

The microphysical, kinetic properties of astrophysical plasmas near accreting compact objects are still poorly understood. For instance, in modern general-relativistic magnetohydrodynamic simulations, the relation between the temperature of electrons $T_{e}$ and protons $T_{p}$ is prescribed in terms of simplified phenomenological models where the electron temperature is related to the proton temperature in terms of the ratio between the gas and magnetic pressures, or $\beta$ parameter. We here present a very comprehensive campaign of {two-dimensional} kinetic Particle-In-Cell (PIC) simulations of special-relativistic turbulence to investigate systematically the microphysical properties of the plasma in the trans-relativistic regime. Using a realistic mass ratio between electrons and protons, we analyze how the index of the electron energy distributions $\kappa$, the efficiency of nonthermal particle production $\mathcal{E}$, and the temperature ratio $\mathcal{T}:=T_{e}/T_{p}$, vary over a wide range of values of $\beta$ and $\sigma$. For each of these quantities, we provide two-dimensional fitting functions that describe their behaviour in the relevant space of parameters, thus connecting the microphysical properties of the plasma, $\kappa$, $\mathcal{E}$, and $\mathcal{T}$, with the macrophysical ones $\beta$ and $\sigma$. In this way, our results can find application in wide range of astrophysical scenarios, including the accretion and the jet emission onto supermassive black holes, such as M87* and Sgr A*.


INTRODUCTION
Considerable effort has been dedicated over the last few years to the modeling via general-relativistic simulations of plasma accreting onto supermassive black holes (Nathanail et al. 2020;Del Zanna et al. 2020;Ripperda et al. 2019;Younsi et al. 2020;Dihingia et al. 2022) and neutron stars (Parfrey & Tchekhovskoy 2017;Abarca et al. 2018;Das et al. 2022;Ç ıkıntoglu et al. 2022). Among the different approaches considered, surely general-relativistic magnetohydrodynamics (GRMHD) simulations have been the focus of many groups worldwide (see, e.g., Event Horizon Telescope Collaboration et al. 2019a;Porth et al. 2019;Event Horizon Telescope Collaboration et al. 2022a). While essential to make theoretical progress on these scenarios, these GRMHD simulations can only describe the dynamically important part of the fluid, the protons (or "ions" as they are sometimes referred to), leaving completely undetermined the physical Corresponding authors: Claudio Meringolo, Alejandro Cruz-Osorio claudiomeringolo@unical.it, osorio@itp.uni-frankfurt.de properties -such as the energy distribution, the number densities, and the temperatures -of the "lighter" part of the fluid, namely, the electrons. This represents a serious limitations for two different reasons. First, in hot, ionized plasma jets around black holes, the Coulomb coupling between electrons and protons is inefficient, so that protons and electrons are likely to have distinct temperatures, as it happens in the solar wind (Tu & Marsch 1997;van der Holst et al. 2010;Howes 2010;Dihingia et al. 2022). Second, a proper knowledge of the electron energy distribution is essential in order to obtain accurate imaging of supermassive black holes and hence compare with the observations (Davelaar et al. 2019;Mizuno et al. 2021;Cruz-Osorio et al. 2022).
To cope with this problem, a number of phenomenological prescriptions have been suggested in the literature to relate the electron temperature to the simulated proton temperature. In this context, a very commonly employed approach is the so-called R−β model (Mościbrodzka et al. 2016), where the electrons temperature is related to the protons temperature in terms of the plasma-β parameter, i.e., the ratio of the thermalto-magnetic pressure, and of two free parameters, R low and R high (see also Anantua et al. 2020, for a critical-β model, where two additional parameters are introduced). The R−β approach has been widely used by the Event Horizon Telescope (EHT) Collaboration to reconstruct theoretically the first images of supermassive black holes M87* (Event Horizon Telescope Collaboration et al. 2019b), and Sgr A* (Event Horizon Telescope Collaboration et al. 2022b). These investigations, in particular, have resorted to a simplified version of the R−β approach in which R low = 1 and spanning different values of R high (Event Horizon Telescope Collaboration et al. 2019a, 2022a. Taking into account a more realistic description of the plasma parameters using self-consistent kinetic models has shown that finer details of the image can appear, but also that the R − β approach is remarkably robust .
Clearly, it is essential to connect the microphysical properties of the plasma with the macrophysical ones β and σ, where hybrid-kinetic models might have some limitations (Arzamasskiy et al. 2019;Valentini et al. 2014;Cerri et al. 2017). To this scope, we have performed 38 large-scale fully kinetic (i.e., both protons and electrons are treated as particles) Particle-In-Cell (PIC) simulations of special-relativistic plasma in the so-called "trans-relativistic regime", that is, when the plasma magnetization σ -the ratio between the magnetic energy density to enthaply density (see below for a definition) -is of order unity (Ripperda et al. 2019;Mizuno et al. 2021;Bandyopadhyay 2022;Janssen et al. 2021), and covering four orders of magnitude in the plasma-β parameter (see Appendix for details on the various simulations). In all simulations, we employ a physical proton-to-electron mass ratio (see Rowan et al. 2017, for the importance of using a realistic mass ratio), and analyze the most important microphysical properties of the turbulent plasma, namely, the spectral index of the electron energy distributions κ, the efficiency in the production of nonthermal particles E, and the temperature ratio T := T e /T p . The parameter ranges explored here overlap and extend those considered in previous and influential works of astrophysical kinetic turbulence (Zhdankin et al. 2019;Zhdankin 2021;Kawazura et al. 2019Kawazura et al. , 2020. Exploiting the large coverage of the space of parameters, we are able to model via analytic fitting functions the behaviour of all of these quantities, thus providing a convenient tool to introduce kinetic effects in global GRMHD simulations of accretion onto compact objects, thus improving the modeling of radiatively inefficient accretion flows around black holes, such as M87* or Sgr A* (Tchekhovskoy & McKinney 2012;Qian et al. 2018;Porth et al. 2019;Cruz-Osorio et al. 2022;Chatterjee et al. 2021;Ripperda et al. 2022;Event Horizon Telescope Collaboration et al. 2022a).

SIMULATION SETUP
To carry out our investigation and to simulate the full development of relativistic turbulence in kinetic plasmas, we use the publicly available Zeltron code (Cerutti et al. 2015;Cerutti & Werner 2019). In particular, we employ a two-dimensional (2D) geometry in Cartesian coordinates retaining, as in any fully consistent plasma model, the threedimensional components of the magnetic and electric fields, of the current density, and of the pressure tensor. The temperature of each species, α = e, p for electrons and protons, is specified through the plasma-β parameter β α := 8πn α k B T α /B 2 0 , where n α and T α are the number densities and the temperatures, respectively. Here, B 0 = (0, 0, B 0 ) is the magnetic-field vector in the ambient plasma (B 0 = const.), and k B is the Boltzmann constant.
We initialize all of our simulations with the same number density for electrons and protons in a computational domain that is a Cartesian box of side L x = L y = L = 16384 dx, where dx = d e /3 is the cell resolution and d e := c/ω pe is the electron-skin depth. In the above, c is the speed of light and ω pe := (4πn e e 2 )/m e [1 + θ e /(Γ e − 1)] −1/2 is the electron plasma frequency, m e the electron mass, Γ e the electron adiabatic index, and θ e := k B T e /m e c 2 the dimensionless electron temperature. We have also carried out three more simulations with a smaller box, L = 2730 d e , and resolution of (8192) 2 , (16384) 2 , and (32768) 2 points (see Appendix for details). Furthermore, we set up our computational box so that it is periodic in the xand in the y-directions. Finally, in all our simulations, each computational cell is initialized with 10 particles per cell (i.e., 5 electrons and 5 protons). As a result, during our evolutions we follow the dynamics of ∼ 2.7 × 10 9 particles per simulation. We initialize our system as done in typical simulations of plasma turbulence (see, e.g., Servidio et al. 2012). The initial conditions consist of a relativistic plasma perturbed by a 2D spectrum of Fourier modes for the magnetic field. To avoid any compressive activity, neither density perturbations nor parallel variances (i.e., with components out of plane) are imposed at t = 0. In practice, we start from expressing the z-component of the vector potential in Fourier modes as A z (x, y) := kx,ky A k exp [i(k · x + φ k )], where k = (k x , k y ) is the wavevector with modulus k = |k| = (2π/L)m (m is the dimensionless wavenumber), and φ k are randomly chosen phases. The amplitude of the modes is set as A k = 1 + (k/k 0 ) 15/3 −1 , such that it is peaked at k 0 = (2π/L)m 0 with m 0 = 4. The spectrum is set to zero at m > 7 in order to construct initial conditions consistent with random large-scale structures. The magnetic-field components B x and B y are then computed by straightforward derivatives. Finally, to explore a regime of strongly perturbed field lines, we fix the amplitude of the fluctuations to be B ⊥ /B 0 ∼ 1, where B ⊥ is the root-mean-square value of the in-plane fluctuations. This choice leads to a broader particle energy distribution, while in Nättilä & Beloborodov (2022) the authors showed that when the amplitude is small, the particle energy distribution is quasi-thermal.
Other quantities that will be referred to in the rest of the paper are the Alfvén crossing time as t A := L/v A , where v A := c σ/(1 + σ) is the Alfvén speed. The plasma magnetization is instead defined as σ := B 2 0 /(4πw), where w is the enthalpy density of the plasma w := (ρ e +ρ p )c 2 +Γ e e + Γ p p , with ρ e,p and e,p being, respectively, the rest-mass densities and the internal energy densities of electrons and protons when following an ideal-fluid equation of state (Rezzolla & Zanotti 2013).
As the simulation proceeds, turbulent magnetic reconnection takes place, leading to a nonlinear change in magnetic topology and converting magnetic energy into kinetic and internal energy. This process strongly affects the dynamics of the plasma on all the scales we could reproduce with our simulations. This highly dynamical system evolves with magnetic flux ropes moving, colliding, and sometimes repelling each other depending on the magnetic-field polarity. This dynamics proceeds till a stationary state is achieved after about an Alfvén crossing time (see also the top panels of Fig. 2 and the Appendix for a detailed discussion).

RESULTS
At the initial time, after fixing β and σ for each simulation, we set the temperatures, which are uniform for both species. In particular, we first specify the proton-β p parameter and then obtain the electron temperature so as to have a specific initial temperature ratio T 0 := T e,0 /T p,0 . At any time during the simulation we measure the spatial distributions in the (x, y) plane of β tot := β e + β p and T , from which we compute the joint PDFs for the two quantities, namely, (β tot , T ). The temperature for each species is computed from T α := p α /n α k B , where p α is the isotropic pressure, i.e., p α := 1 3 (p xx α + p yy α + p zz α ) and p ij α is the pressure tensor. This is shown in Fig. 1, where we report the joint PDFs at the initial (t = 0) and final (t = 2 t A ) times for three representative simulations initialized respectively with T 0 = 0.1, 1.0, and 10. Clearly, the three initial setups have different joint PDFs narrowly distributed around the three initial values of the temperature ratio T 0 . Interestingly, however, at the final time they have all converged to the same equilibrium distribution, irrespective of the initial data. This can be best appreciated in the inset, which reports a zoom-in of the central region of the final distributions, with the color-coded contour reporting the 90%-value for each simulation, while the circles represent the maximum of each joint PDF. This convergence has been verified to take place for four different values of the initial temperature ratio (T 0 = 0.01, 0.1, 1.0 and 10.0), while keeping σ = 0.3 and β = 3 × 10 −4 . The behaviour in Fig. 1 induces us to conjecture that the choice of the initial temperature T 0 is effectively unimportant at least in the ranges explored here 1 as its memory is lost by the time the system has reached a steady state. In view of this, we set T 0 = 1.0 for the 35 simulations performed varying σ and β (note that with such initial temperature ratio, the plasma-β parameter is the same for electrons and protons, i.e., β e = β p =: β). The ranges of σ and β explored are compatible with previous kinetic studies, state-ofthe-art GRMHD simulations, and radiative-transfer calculations (Ball et al. 2018;Cruz-Osorio et al. 2022;Fromm et al. 2022). As also noted by Pecora et al. (2019), higher values of β would require a much higher number of particles to counter the statistical noise, making purely PIC calculations of this type computationally expensive with modern resources. Figure 2 provides a very compact but powerful overview of the fully developed turbulent state for a simulation with σ = 1.0 and β = 3 × 10 −3 , at time t = 1.5 t A . Each upper panel is split into two regions reporting different plasma properties. Panel (a) shows the electron number density n e normalized to the initial number of particles per cell n 0 (left), and the magnetization σ (right). Panel (b), instead, reports temperature ratio T (left) and the out-of-plane electric-current density J z (right). Note how, in analogy to nonrelativistic kinetic simulations, vortex-like and sheet-like structures corresponding to magnetic flux tubes are present at all the scales that are resolved in the simulation (Servidio et al. 2012;Comisso & Sironi 2018;Parashar et al. 2018;Pecora et al. 2019). High number-densities "magnetic islands" can be found in large-scale flux tubes, and in general, the density is larger in these coherent quasi-circular structures.
At the same time, the largest temperatures (ratios) are not achieved at the center of the islands, which are instead comparatively cooler. This is because the temperature is higher between flux tubes, where reconnection layers lead to the formation of plasmoids within narrow current sheets (Servidio et al. 2009;Comisso & Sironi 2018;Pezzi et al. 2021). Elongated unstable current sheets tend to fragment into chains of plasmoids and small-size current sheets are appear on a wide range of scales (Hellinger et al. 2015;Dong et al. 2018;Huang & Bhattacharjee 2016). Notice also that the out-ofplane electric-current density J z shows a variety of current sheets of different sizes. Some of these current layers break into smaller plasmoids and these regions are important for the heating of the plasma and the acceleration of the particles.
The various quantities shown in Fig. 2 are overlaid with the trajectories of some of the most energized particles that we tracked (protons in the left panels and electrons in the right ones). In particular, we track a sample of 500 electrons and 500 protons during the whole simulation, both randomly chosen. The starting-position of each particle is marked with a star. Note how, quite generically, and in addition to the basic gyrations at the corresponding Larmor radii, there are particles that have closed orbits as they are trapped in a flux rope, while others experience turnovers that suddenly bend the trajectory, similarly to what observed in nonrelativistic turbulence simulations (Pecora et al. 2018) 2 . Overall, when a particle experiences a reconnection process and is accelerated, it increases abruptly its Larmor radius, but also its Lorentz factor γ, and kinetic energy.
In the lower panels (c) and (d) of Fig. 2 we show the evolution of the Lorentz factor of the particles tracked in the upper panels (a) and (b), with protons being reported in panel (c) and electrons in panel (d). As expected, and shown by the different vertical scales of panels (c) and (d), electrons experience considerably larger accelerations when compared to protons. This is simply due to the different masses of the two species: electrons, which have smaller Larmor radius, are more efficiently accelerated by the thin current sheets where magnetic reconnection takes place. This stochastic acceler-ation mechanism of multi-reconnection events is very efficient and commonly observed in astrophysical plasma turbulence (Drake et al. 2009;Haynes et al. 2014).
The tracked particles start from γ 1, and most of them experience a sudden acceleration episode, and then a sequence of second-order Fermi-like processes of acceleration (Comisso & Sironi 2018. Particles trapped in magnetic islands show a Lorentz factor increasing in time (e.g., the red proton in the left panels). Other particles, instead gain energy only once and then reach a quasi-steady state as is typical for particles entering the magnetic island only for a short time and then being bounced in a stochastic manner between different structures.
Relativistic hydrodynamical turbulence naturally provides a landscape of intermittency and large spatial variance because the compressibility is enhanced by relativistic effects (Radice & Rezzolla 2013); in addition, relativistic magnetohydrodynamical turbulence provides the natural conditions to produce extreme-acceleration events and to generate a large population of particles -electrons in particular -with energy distributions that differ significantly from a thermal one (see, e.g., Zhdankin et al. 2017). This is summarized in Fig. 3, which reports the electron energy-distribution functions (spectra) (γ − 1)dN/dγ at t = 2 t A as a function of the Lorentz factor γ − 1, for some representative simulations. More specifically, the upper panel shows the electron spectra from simulations with σ = 0.3 and for a wide range of values of β; the black dashed line is a Maxwell-Jüttner distribution where the value of the dimensionless electron temperature θ e := k B T e /(m e c 2 ) = 45 is chosen to reproduce the lowenergy part of the spectrum for the case β = 0.11 and is obviously different for each simulation. Note that the highenergy part of the spectra is well approximated by a powerlaw dN/dγ ∝ γ −κ+1 (Davelaar et al. 2019;Fromm et al. 2021), whose index κ 3.2 is quite insensitive to the value of the plasma-β parameter in the range β . Also overplotted with different colors are representative particle trajectories, with protons on the left and electrons on the right of each panel (the initial position of each particle is marked with a star). The lower panels [(c) and (d)] report instead the evolution of the Lorentz factor for the same particles marked above.
dicate that turbulence promotes the particle acceleration, producing energy distributions that contain a considerable fraction of very energetic (suprathermal) particles. Given the kinetic behaviour of the plasmas described so far, it is essential to be able to express their properties via analytic fitting functions and in terms of the basic parameters of the plasma, namely, β, σ, so that the resulting expressions can then be employed directly in the GRMHD modelling of astrophysical plasma. A summary of this analytical modelling is presented in Fig. 4, where in the top row we show as a function of β and σ, respectively, the electron spectral index κ, the nonthermal energy efficiency E, and the temperature ratio T . Note that the data reported in the first two columns refers to simulations at t = 2 t A , while that in the right column is averaged over the time window 1.7 < t/t A < 2.3 to avoid the oscillations introduced by the stochastic behavior of turbulence. Similarly, the bottom row of Fig. 4 reports one-dimensional cuts of the same quantities, but at fixed val-ues of the magnetization (σ = 0.1 − 3.0), where each circle refers to a distinct simulation of our set. Note that for any fixed value of σ we explored plasma parameters up to the maximum one β max ∼ 1/(4σ) (Ball et al. 2018), where our estimates are inevitably less accurate.
Exploiting the large set of simulations performed, we can now construct analytical 2D fits to the various quantities, starting with the electron spectral index κ(β, σ), which can be expressed as where k 0 = 2.8, k 1 = 0.2, k 2 = 1.6 and k 3 = 2.25 (see top-left panel of Fig. 4). Note that Zhdankin et al. (2017) have proposed a similar but simpler fitting expression which depends σ only and thus does not account for variations in the plasma β. Overall, the spectral index shows two main features. First, at fixed σ, the spectral index is essentially independent of β, for β 10 −2 , but it increases at larger values of β, approaching a very steep tail. Second, at fixed β, the index becomes generally smaller for increasing values of σ.
Next, we quantify the efficiency in the production of particles with nonthermal energies in terms of the weighted average of the excess over a Maxwell-Jüttner distribution (Ball et al. 2018), namely where γ 0 denotes the peak of the spectrum, f MJ := γ 2 v/[c θ e K 2 (1/θ e )]e −γ/θe , with v the velocity and K 2 the modified Bessel function of the second kind. The corresponding 2D fit of the data can then be expressed as E(β, σ) = e 0 + e 1 √ σ + e 2 σ 1/10 tanh e 3 β σ 1/10 , where e 0 = 1.0, e 1 = −0.23, e 2 = 0.5 and e 3 = −10.18 (see top-middle panel of Fig. 4). Also in this case, the energy efficiency shows three main features. First, for β 10 −2 the efficiency saturates at a value that is independent of β, but systematically larger for higher values of σ. Second, for high values of β and low values of σ, it approaches E ∼ 0, because the electron spectrum becomes significantly softer. Third, for higher values of σ, the efficiency is the largest, since the spectra widen to larger electron energies. Interestingly, these results are similar to the ones found by Ball et al. (2018) when using different initial conditions. Finally, we consider what is arguably the most important quantity modelled in our simulations, namely, the dependence of the temperature ratio on the plasma properties. The corresponding 2D fit is given by where t 0 = 0.4, t 1 = 0.25, t 2 = 5.75, t 3 = 0.037, and τ 1 = −0.5, τ 2 = 0.95, τ 3 = −0.3, τ 4 = −0.05 (see topright panel of Fig. 4). Overall, it is easy to see that for low magnetizations, i.e., σ ∈ [0.1, 0.3], and small values of the β parameter, i.e., β 0.01, the temperature ratio is essentially constant and then starts to grow to values as large as T 1 for β 1.0. On the other hand, for high values of the magnetization, i.e., σ 3.0, the behavior is quite the opposite, the values of T are higher for lower β and decrease when increasing β. For intermediate values of the magnetization, i.e., σ = 1.0, the behavior is a combination of the two described above, showing a nonmonotonic dependence for β ∈ [0.01, 0.1]. Interestingly, in all cases, T ∼ 1.0 for β 1, independently of the value of σ, thus highlighting that, under these conditions, electrons and protons are fully coupled and have roughly the same temperature. Conversely, for β 10 −4 , the temperature ratio will depend on the plasma magnetization, being larger for larger magnetizations, as expected for regimes where electrons can be accelerated to suprathermal energies at reconnection sites. More importantly, expression (4) provides a compact and microphysically consistent description of the electron temperatures that can be employed in modern GRMHD codes of accretion flows onto black holes.
We conclude the discussion of our results by returning to the behaviour of the electron spectral index κ. As shown in the top-left panel of Fig. 4 and summarized in Eq. (1), electron acceleration is higher in low-β and high-σ turbulent plasmas. As suggested already by Drake et al. (2009), this behaviour may be due to the interaction of the electron orbits with small-sized current sheets; such a mechanism can then extract particles from the thermal population and bring them to very high energies via primary and secondary Fermi-like mechanisms (Pecora et al. 2018;Comisso & Sironi 2018).
In fully developed GRMHD turbulence, accelerating islands and current sheets are present on all scales and these could therefore provide the natural site for the accelerating mechanism.
In this simple picture, it is natural to expect that the larger the spectrum of fluctuations at small scales, the more efficient the accelerating mechanism (Haynes et al. 2014). To validate whether this applies also to trans-relativistic plasmas, we have computed the (not normalized) isotropic power spectrum of the magnetic field for three representative simulations and reported them in Fig. 5 as a function of the dimensionless kd e [the inset shows with colored squares the location in the (σ, β) plane of the three configurations, while the arrows mark the wavevectors associated to the proton-skin depth (kd p = 1) and to the proton Larmor radius (kρ p = 1)] and over a downsampled grid of (1024) 2 (see Appendix for a discussion). In essence, after assuming the turbulence to be isotropic and homogeneous, we integrate the 2D Fourier transforms B i over concentric shells (in this sense, the power spectrum is isotropic) to obtain one-dimensional spectra, whose sum we plot in Fig. 5 [note that the growth of the power spectrum at large wavenumbers is a typical noise effect of PIC simulations due to a finite number of particles, (see, e.g., Karimabadi et al. 2013)].
In general, Fig. 5 reveals a number of interesting features, moving in the parameter space from (low-β, high-σ) to (highβ, low-σ). First, the power spectrum is clearly higher in the case of the low-β, high-σ simulation, confirming a more efficient cascade process (Franci et al. 2016). Second, the spectrum is shallower in the sub-ion inertial range (Sahraoui et al. 2009) indicating a more developed turbulence. Finally, and more interestingly, the turbulent cascades terminate at much smaller scales for (low-β, high-σ) simulations, suggesting the existence of thinner current sheets at subproton scales that accelerate particles more efficiently (Pecora et al. 2018).

DISCUSSION AND CONCLUSIONS
With the goal of gaining a deeper understanding of the properties of plasmas near astrophysical compact objects, we have employed the PIC Zeltron code to carry out a large campaign of two-dimensional simulations of special-relativistic, decaying plasma turbulence in the transrelativistic regime. Particularly important in our analysis is the use of a physical mass ratio between electrons and pro-  Figure 5. Magnetic-field power spectra for three simulations sampling important locations in the (β, σ) space of parameters. Each simulation is marked with a different color and the corresponding location is shown in the inset, which reports also the electron spectral index. Black dashed lines indicate the turbulent power laws, while the circles delimit the boundaries of each turbulent range, which we define as the limits of the power-law scaling; the arrows mark the wavevectors associated to the proton-skin depth (kdp = 1) and to the proton Larmor radius (kρp = 1), which is outside the horizontal scale for the red line.
tons and the exploration of a wide range of values in the plasma-β parameter (β = 10 −4 − 1.5) and in the magnetization σ (σ = 0.1−3.0). Having simulated such a large portion of the space of parameters encountered in astrophysical plasmas has allowed us to derive analytical fitting functions for the behaviour of a number of important plasma quantities as a function of β and σ. More specifically, we have presented 2D fitting functions of the electron spectral index κ(β, σ), of the efficiency in generating nonthermal particles E(β, σ), and of the ratio between the electron and proton temperatures T (β, σ). These expressions provide compact and reasonably accurate descriptions of the behaviour of these microphysical plasma properties and can be employed in a number of scenarios involving compact objects and described by macro-physical plasma characteristics. Importantly, since they have been derived from first-principle calculations, they represent a considerable improvement over the rather crude and purely empirical expressions employed at the moment in GRMHD simulations. Finally, we have confirmed the suggestion that plasmas with low β and large σ naturally lead to broad turbulent scenarios and are the most efficient in extracting particles from the thermal population and accelerating them (Pecora et al. 2018;Comisso & Sironi 2018).
As these simulations represent one of the most systematic PIC explorations of trans-relativistic turbulence, can be employed in a wide range of astrophysical systems, such as jets and accretion disks around supermassive black holes, and, of course, their imaging (see, e.g., Event Horizon Telescope Collaboration et al. 2019aCollaboration et al. , 2022a. The formulas provided in this work can be improved by extending the present twodimensional treatment to three dimensions and thus assessing the role played by dimensionality in studies of this type.

ACKNOWLEDGMENTS
We thank the Referee for the useful comments that have improved our presentation. This research is supported by the ERC Advanced Grant "JETSET: Launching, propagation and emission of relativistic jets from binary mergers and across mass scales" (Grant No. 884631

APPENDIX
In what follows, we provide additional information on our analysis concentrating on three specific aspects: a detailed summary of the properties of the simulations carried out in the campaign, the evidence that stationarity is reached when extracting the spectral information, and a comparison of simulations with different resolutions.
As a first test, to show that our final configuration is independent of the initial electron-to-proton temperature, we have varied T 0 spanning in the range [0.001−10.0] (see Runs A1-A3 in Table 2, see Fig.1). Note that for these configurations, the plasma β is different for electrons and protons. Next, we checked that our results are insensitive to the choice of different (higher) resolutions in terms of d e /dx, increasing the resolution up to d e /dx = 12 (see Runs B1-B3 in Table 2). In the latter case, we have used a physical box of L/d e = 2730 in both directions and varied the number of mesh points from (8192) 2 up to (32768) 2 . In this last high-resolution configuration, we have followed the dynamics of ∼ 1.1 × 10 11 particles.
STATIONARITY OF SPECTRA Next, we provide evidence that the computed electron-energy spectra reach a steady state after t/t A 1.8 − 2.0, so that the extraction of the spectral index κ and of the efficiency E is both accurate and robust. Figure 6 shows four representative simulations having different values of σ (see Runs 7, 18, 27, and 31 in Table 1). In each case, we plot the electron-energy spectra at different times during the evolution as indicated by the colormap on the right of each of the four panels. Furthermore, marked with black vertical lines of various type are three different values of the Lorentz factor γ − 1 and the corresponding evolutions are shown in the bottom panels for each of the four simulations considered. Clearly, all cases show that by t/t A ∼ 2.0 the simulations have reached stationarity with relative time variations that are 1.5%, so that κ and E can be extracted reliably.
RESOLUTION TESTS Finally, we have verified that our results are insensitive to the choice of spatial resolution. In particular, we have performed three simulations using an increasing number of cells per electron skin depth, from d e /dx = 3 up to d e /dx = 12 (see Runs  Table  2). For each simulation, we report the spectra at different times during the evolution as indicated by the colormap on the right of each of the four panels. Marked with black vertical lines of various type are three different values of the Lorentz factor γ − 1 and the corresponding evolutions are shown in the bottom panels for each of the four simulations considered. Clearly, all cases show that by t/tA ∼ 2.0 the simulations have reached stationarity. Table 2). Figure 7 compares the electron-energy spectra for a case with σ = 0.3 and β = 3 × 10 −4 when varying the number of electron-skin depths per cell, i.e., d e /dx = 3 − 12. Clearly, the main features of the electron-energy spectra and in particular the slope are very similar for the three different resolutions. Indeed, the relative differences between the three spectra are 6.0% and thus even smaller than the variations due to the stochastic nature of turbulence, which can cause variation in κ up to ∼ 10.0% (Ball et al. 2018  . Electron-energy spectra with σ = 0.3 and β = 3 × 10 −4 for three different resolutions de/dx = 3, 6 and 12, using a physical box size of L/de = 2730. The spectra are computed at t/tA = 2.0 and clearly show to be nearly insensitive to the increased resolution.

B1-B3 in
In Figure 8 we show the joint PDFs for the ratio of temperatures T and the plasma β tot = β e + β p for the same runs. In the inset we report a zoom-in of the central region of the PDFs at the final time of t = 2 t A . The color-coded contour report the 90%-value for each distribution, while the circle represent the maximum of each joint PDF. One can see that for the three different resolutions we obtain similar final distributions, with a variation in T 5.0%.
As a concluding remark, we note that the power spectrum in Fig. 5 has been computed on a down-sampled grid of (1024) 2 points and not on the full-resolution data of (16348) 2 points. This coarse-graining operation is routinely done in such expensive simulations, and for two distinct reasons. First, the large particle noise due to the high temperatures reached essentially blurs out the smallest scales, so that using the full resolution does not really provide any additional information. Second, the downsampling allows us to reduce by a factor of 16 2 ∼ 250 the space needed for the ouput (we recall that we save data for 38 fields at very high cadence). As a result, while the simulation maximum wavenumber is k max d e = 9.4 and is not shown in the spectrum in Fig. 5, the maximum wavenumber in the downsampled spectrum is k max d e = 0.6 and is well-captured.