How the turbulent energy cascade develops below the magnetohydrodynamic scales in space and astrophysical plasmas is a major open question. Here, we measure the power spectrum of magnetic fluctuations in Parker Solar Probe's observations close to the Sun and in state-of-the-art numerical simulations of plasma turbulence. Both reveal a power-law behavior with a slope compatible with −11/3 at scales smaller than the ion characteristic scales, steeper than what is typically observed in the solar wind and in the Earth's magnetosheath. We explain such behavior by developing a simple two-fluid model which does not require any kinetic processes nor electron-inertia effects. This is characterized by a significant contribution of the ion kinetic energy to the total turbulent energy cascade at sub-ion scales, although the dynamics is driven by the magnetic field through the current density. We expect that this regime may be relevant for a broad class of low-beta plasmas, e.g. the solar corona, non-relativistic magnetized jets and disks, and laboratory plasmas.
The solar wind, the collisionless plasma stream of charged particles that continuously flows outward from the Sun, is heated from thousands to millions of degrees over a distance of just a solar radius and accelerated up to speeds of 800 km/s. The mechanisms responsible for such a behavior are not yet well understood. In-situ spacecraft measurements close to the Earth suggest that plasma turbulence is a key ingredient for the solar wind dynamics and therefore a possible candidate for heating and accelerating its particles [e.g., [1] [2] [3] [4] [5] . The observed power spectra of the plasma and electromagnetic fluctuations exhibit a power-law behavior spanning many decades in frequency [e.g., 6, 7] , thus providing evidence that an extended turbulent cascade is at play. This is responsible for transferring energy via nonlinear coupling from large down to smaller and smaller scales, through those characteristic of ions and electrons, where it is eventually dissipated by kinetic heating mechanisms other than particle collisions.
Particle heating and energy dissipation in the absence of collisions are also important phenomena in a number of different plasma environments over an extremely large range of scales, both inside the heliosphere (e.g., in the solar corona [8] ) and outside (e.g., in accretion disks [9] , in the interstellar medium [10] , and in galaxy clusters [11] ). The solar wind is, however, unique among such systems, since it can be directly probed by means of spacecraft missions.
The NASA Parker Solar Probe (PSP) spacecraft mission [12] has been approaching the Sun in an attempt to shed light on the origin of the solar wind. During its first * l.franci@qmul.ac.uk five orbits, PSP has measured the plasma and electromagnetic fluctuations closer than any previous mission. PSP observations have already led to important discoveries, including dramatic magnetic field reversals on time scales of seconds or minutes [13, 14] , in-situ generated electromagnetic waves at ion resonant scales [15, 16] , a significant increase in the energy transfer rate at magnetohydrodynamic scales with respect to that observed at 1 au [17] , an enhanced proton heating associated with large-amplitude turbulent fluctuations [18] , and a steepening of the magnetic field spectrum below the magnetohydrodynamic scales [19] .
Here we present the first evidence, both observational and numerical, for a plasma regime that occurs at scales just below the ion inertial length and is characterized by a power spectrum of the magnetic fluctuations with a slope (spectral index) compatible with −11/3. We also provide a theoretical interpretation based on a simple non-relativistic two-fluid model, which does not require any ion and electron kinetic processes (e.g., wave-particle resonances) or any electron-inertia effects. The model predictions are in agreement with recent spacecraft observations and are supported by results of numerical simulations. The predicted regime occurs when the ion kinetic energy is not negligible with respect to the magnetic energy at sub-ion scales and provides an additional contribution to the energy cascade. The dynamics, however, is still guided by the evolution of the magnetic field via the current density. In this sense, we define this as a "current-mediated" turbulent plasma regime. Based on our simulations, we expect that this is more likely to occur when the turbulent fluctuations are large and the plasma beta, i.e., the ratio of the plasma to the magnetic pressure, is low. As such conditions are more typically met in the near-Sun environment rather than further away, this might explain why it has not been identified before, while it seems to be frequent in PSP observations. The current-mediated plasma regime can account for an enhanced sub-ion-scale energy transfer rate, due to the additional contribution of the ion kinetic energy, and possibly for a stronger turbulent heating close to the Sun than that predicted by alternative turbulent cascade models. This could shed light on the solar wind origin and, more in general, on particle heating and acceleration in collisionless space and astrophysical plasmas.
The observational dataset consists of high-resolution magnetic field measurements by the FIELDS fluxgate magnetometer [22] on-board PSP [12] at 36 solar radii from the Sun during its first perihelion on November 5, 2018 (00:00-24:00 UT) (see Sec. III B for further details). The time series of the magnetic field magnitude and of its tangential, normal, and radial components for the full day are shown in Fig. 1a , from top to bottom. We choose to divide the dataset into 32 consecutive 45-minute intervals, as we intend to focus our attention on the spectral behavior at frequencies larger than ∼ 10 −2 Hz. The spectral properties at larger frequencies for the first two PSP orbits, including the first perihelion, have been recently analysed and discussed in Chen et al. [23] . The first four of our intervals correspond to the ones discussed in Fig. 2 of [13] . Most, although not all, intervals contain reversals of the radial magnetic field component, corresponding to jets (sometimes also called "switchbacks") [13, 14] . Some intervals show instead a more quiet radial-field wind, with very small fluctuations of the radial magnetic field component and a quite constant magnitude. Such behavior is evident in Fig. 1b , which show the same quantities as in Fig. 1a , limited to the time interval #10 (07:30-8:15 UT).
The trace power spectrum of the magnetic field, P B , for interval #10 is shown in the top panel of Fig. 1d (see Sec. III C for the definition of power spectrum). Its slope (spectral index) is well approximated by −5/3 above the ion characteristic scales ("supra-ion" scales), over at least one full decade in frequency, and by −11/3 below (sub ion-scales). The latter behavior holds up to ∼ 10 Hz, where the measurements from the fluxgate magnetometer get close to its noise floor. In the bottom panel of the same figure, the local spectral index, α, obtained by performing power-law fits over sliding windows of half a decade, is shown. The averaged spectrum, i.e., the result of averaging all the 32 spectra, each normalized to the mean field in the corresponding interval, further confirms this behavior (Fig. 1c) : the local spectral index between ∼ 3 Hz and ∼ 8 Hz is −3.56.
The numerical dataset consists of a 2D and a 3D highresolution hybrid-kinetic simulations of decaying plasma turbulence performed with the code CAMELIA [26] . Both simulations are initialized with a uniform plasma, embedded in an ambient magnetic field and perturbed with large isotropic Alfvénic-like fluctuations, i.e., ion bulk velocity and magnetic fluctuations with the same amplitude in Alfvénic units (cf. Sec. III A for all normalization units). The ion and electron plasma betas, i.e., the ratio of particle thermal pressure to the magnetic pressure, were set to be comparable to their average observed values on the day of PSP first perihelion (β i = 0.2, β e = 0.5), in order to mimic similar plasma conditions.
The same fitting procedure for the observational spectra has been applied to the numerical spectra obtained from the 2D (Fig. 2c ) and the 3D case (Fig. 2c) . Both numerical spectra show a power-law range at sub-ion scales with a spectral index compatible with −11/3, although this is less extended in the latter case due to the larger computational cost which limits the spatial resolution. Fig. 3 also shows the spectra of the ion and electron bulk velocities and of the current density for the 2D simulation. The corresponding spectra from the 3D simulation, not shown here, are very similar, especially at kd i 1. The power spectrum of the ion bulk velocity fluctuations, P ui , starts decreasing around kd i 1, consistent with solar wind observations and previous simulations [e.g., 24, 27] . It does not, however, drop off steeper than the magnetic fluctuations, as previously observed [25, 28, 29] , but reduces steadily until the noise level is reached. On the contrary, corresponding to the scale of the transition in the magnetic field spectra, it starts behaving as a power law, with a spectral index comparable to that of P B . An almost constant ratio P B /P ui ∼ 8 is observed, i.e., the ion velocity fluctuations are of the same order as the magnetic fluctuations in Alfvén units. The ion motion is here still contributing to the plasma dynamics at these scales, providing a non-negligible contribution to the total current density, J. This is also evident from the fact that the spectra of the current density and of the electron bulk velocity, u e = u i − J/n exhibit a slightly different level, although with a similar power-law behavior. The slope is compatible with −11/3 + 2 = −5/3, consistent with J = ∇ × B, which can be assumed here at all observed scales due to the non-relativistic regime of fluctuations.
The power spectra of the magnetic fluctuations in both observational and numerical data exhibit a power-law behavior at scales just below the ion-scale transition, with slopes compatible with −11/3. spectra, each normalized to its mean magnetic field. As a reference, power laws are superimposed with slope −5/3, classically predicted for the inertial-range turbulent cascade at supra-ion scales, and with −11/3, predicted by our model at sub-ion scales. Bottom panel: Local spectral index α. Horizontal black dashed lines mark the classical asymptotic predictions for the turbulent cascade in the MHD inertial range, −5/3, and at sub-ion scales for a regular cascade, −7/3 [20] , for an intermittent one, −8/3 [21] , and our prediction, −11/3. The horizontal red line represents the fit of the local slopes over the range of frequencies limited by the two dot-dashed vertical red lines. (d) The same as in (c), but for time interval #10, (07:30-8:15 UT). The spectrum has been normalized to the ambient magnetic field B0 in the time interval.
short-interval spectra, of the 24h-averaged PSP spectrum, and of the spectra of the 2D and the 3D hybridkinetic simulations. From the distribution of their values, we see that: i) there are no spectral indices around the theoretically predicted −7/3 or −8/3 and moreover no interval exhibits a power law flatter than ∼ −3; ii) most of the intervals have a slope which is very close and compatible with −11/3. We can conclude that the −11/3 spectral index of the PSP averaged spectrum is not just the result of a few individual spectra dominating over the others. On the contrary, it seems to be the manifestation of a very frequent behavior over the whole day of PSP first perihelion. Fig. 4b directly compares the observational and the numerical magnetic field spectra, overlapping Fig. 1d and Fig. 2c . The former has been arbitrarily shifted along the x-axis so that the scale of the transition would coincide, with no pretence to use the Taylor hypothesis to accurately convert the frequency in wavenumber.
FIG. 2. Numerical dataset: 2D and 3D hybrid-kinetic simulations of plasma turbulence initialized with values of ion and electron plasma beta similar to those encountered on average at PSP first perihelion, here shown at the time of the maximum turbulent activity. (a) Contour plot of the squared magnitude of the magnetic fluctuations, |B| 2 , for the 2D simulation. The plasma is evolved in a 256 di × 256 di domain and the ambient magnetic field B0 = B0ẑ is out of the plane. The system looks qualitatively similar to what have been previously obtained with simulations reproducing plasma conditions typical of the solar wind [24] and of the Earth's magnetopause [25] . The main difference is the larger amplitude of the magnetic field gradients, i.e., stronger current sheets, and their larger filling factor. (b) Contour plot of |B| 2 for the 3D simulation on a 64 di × 64 di × 64 di domain, with the same physical parameters as the 2D one. Despite the initial fluctuations being isotropic in wavevetor domain, small-scale gradients mainly develop in the x, y plane perpendicular to B0 and they look elongated along it, as in previous similar 3D simulations [24] . (c) Isotropized 1D power spectrum of the total magnetic field fluctuations for the 2D simulation. (d) The same as in Fig. 1c , but for the 3D simulation. .
Here we provide a theoretical interpretation for the −11/3 spectral index in terms of a simple non-relativistic two-fluid model which does not require any ion and electron kinetic processes (e.g., wave-particle resonances) or any electron-inertia effects. In the following equations, we will normalize the magnetic field to the magnitude of the ambient magnetic field, B 0 , the ion and electron bulk velocities to the Alfvén velocity, V A , the ion and electron densities to the ambient density, n 0 , and lengths to the ion inertial length, d i . We will also use this constraint between the magnetic field and the current density
which follows from neglecting the displacement current in Ampère's law in the non-relativistic regime. Here and in the following equations, the subscript k denotes the Fourier transform of a quantity. In the limit of small density fluctuations, n i,e ∼ n 0 , the current density can be written as
(2)
In the same limit, we can express the total energy density at sub-ion scales in normalized units as
The scale-invariant transfer rate is then E k /τ nl ∼ ε, where τ nl is a characteristic nonlinear time expressing the scaleto-scale energy coupling. At supra-ion scales, both observations and simulations indicate that the fluctuations are Alfvénic-like, i.e., u 2 i,k ∼ B 2 k . The nonlinear time associated with the scale-to-scale coupling is therefore τ sup ∼ 1/(ku i,k ) ∼ 1/(kB k ). Demanding a scale-independent energy transfer rate yields the well-known −5/3 power-law magnetic field spectrum. Indeed, we observe this power law in both the PSP and the simulation spectra at the largest scales here investigated.
At sub-ion scales, Fig. 3 shows that the ion bulk velocity spectrum follows a power-law behavior which is similar to the one of the magnetic field, although with a lower level. Using Eq. 1, in the sub-ion range we have
As a consequence, the current density is mainly due to the electrons although the ions still provide a nonnegligible contribution to the plasma dynamics. Based on this and on the similar spectral behavior of P ue and P J observed at sub-ion scales (cf. Fig. 3) , we define the ratio Fig. 1a and Fig. 1c) , for the averaged PSP spectrum (ID 33, blu square, see Fig. 1d ), for the 2D simulation (ID 34, red diamond, see Fig. 2c ), and for the 3D simulation (ID 35, pink triangle, see Fig. 2d ). The grey points mark the intervals where the automatic fit interval has less than a factor of 2 between its lower and upper frequencies. In this case, a power-law behavior is either not present or it's extension is very small and so not significant enough. The horizontal dashed black lines mark the theoretically predicted spectral indices −7/3 and −8/3, while the blue dashed line marks the prediction of our model, −11/3. The two light blue shaded areas mark 5% and 10% around this value. (b) Direct comparison between the observational and the numerical magnetic field spectra as the result of overlapping Fig. 1d and Fig. 2c . The former has been arbitrarily shifted along the x-axis for the purpose of comparing the slope, with no pretence to use the Taylor hypothesis to accurately convert the frequency in wavenumber.
and assume it to be weakly dependent on k in that range, such that it tends to a negative constant of the order of 1 for sufficiently large k, i.e.,
Note that it has been implicitly assumed that w k is the FIG. 5. Schematic summary of the theoretical model for the plasma regime characterized by power-law magnetic field spectrum with a spectral index of −11/3. The main requirement for the regime is that ion bulk velocity fluctuations are not negligible with respect to the magnetic fluctuations, so their spectra are relatively close to each other. As a consequence, the ion motion provides a non-negligible contribution to the current density, which is however still dominated by the electron motion. At sub-ion scales, the spectrum of J exhibits a spectral index of −5/3 and the one of ue is close to it although not coincident. As a consequence of the constraint expressed by Eq. 1 (black dash-dotted line), this results in a −11/3 spectral index for the magnetic field spectrum.
same for all vector components. This has been verified in the 2D simulation, where we find that P u i,⊥ /P J ⊥ = P u i, /P J = P ui /P J at scales below the break, which is the range we are interested in. Using Eq. 2, at sub-ion scales we can write the ion bulk velocity fluctuations in terms of the current density as
By substituting Eq. 7 and Eq. 1 into Eq. 3, one can express the total energy density in the sub-ion range as a function only of the scale k and of the magnetic fluctuations at such scale through
The two terms on the RHS of Eq. 8 can now be compared to each other to determine the scaling of the magnetic field spectrum.
In the limit ǫ 2 k 2 ≪ 1, we can approximate the energy at sub-ion scales as E k ≈ B 2 k . We obtain the nonlinear time associated to the evolution of the magnetic field from Faraday's law as τ nl ∼ (ku e,k ) −1 ∼ (k 2 B k ) −1 , where we have used that at sub-ion scales u e,k ∼ J k = kB k . By requesting the energy flux E k /τ nl ∼ k 2 B 3 k to be scaleindependent, we obtain P (B) = B 2 k /k ∼ k −7/3 . This includes the particular case ǫ → 0, when u i,k → 0 and the current density is all due to the electrons, i.e., the electron-MHD limit [e.g., 30, 31] .
If, on the contrary,
we can approximate the energy at sub-ion scales as E k ≈ ǫ 2 k 2 B 2 k . This term, although coming from the ion kinetic energy, is expressed in terms of the magnetic field via the current density. Therefore, the relevant nonlinear time is again the one associated to the evolution of the magnetic field, i.e., τ nl ∼ (k 2 B k ) −1 . Repeating the same argument as above, we now obtain P (B) ∼ k −11/3 . The condition expressed by Eq. 9 includes the case ǫ ≡ −1, when u e,k = 0 and the current is all due to the ions, i.e., ion-MHD [32] . This, however, represents an extreme case which does not correspond to what we observe in the numerical spectra (see Fig. 3 ).
More in general, we can expect the −11/3 power-law regime to occur below a threshold k thr ∼ 1/ǫ. The threshold will move towards smaller scales the closer ǫ gets to 0. In other words, the smaller the ion contribution to the current, the smaller the chances that such regime will occur at scales close to the ion-scale transition. As a consequence, we should expect the currentmediated regime to be observed at "large sub-ion scales" only when the current density is much larger than the ion bulk velocity which, in turn, is not negligible with respect to the magnetic field.
It is important to note that, although Eq. 6 and the last term of Eq. 7 strictly hold only for k ≫ 1, we use them in the range 3 k 20, where the dependence of w k on k might actually not be negligible. If we assume that such dependence is the form w k = −(1 + ǫ) + a/k + O(1/k 2 ) (where a is a constant), we obtain u e,k ∼ (−1 + a/k − ǫ)J k and u i,k ∼ (a/k − ǫ)J k . At the leading order, this yields u e,k ∼ −J k and u i,k ∼ aJ k /k ∼ aB k . This would explain why the spectrum of u i,k has a similar slope to that of B k rather than J k (cf. Fig. 3 ). Although we do not investigate here the physical meaning of w k and of its dependence on k, we speculate that this might be somehow related to the alignment between the ion and the electron bulk velocity fluctuations, which could be scale-dependent. Fig. 5 summarizes the main features of the theoretical model. Further details about it can be found in [33] .
All the above arguments has been made in k-space. Whenever the Taylor hypothesis holds, however, they can be also applied to the frequency power spectra measured by spacecraft, by using f = ck to link the frequency, f , to the wavenumber, k, through a constant that depends on the wind speed.
This work provides the first combined observational and numerical evidence of a new plasma turbulent regime at scales just below the ion inertial length. This is characterized by: (i) a power-law spectrum of the magnetic field with a spectral index of −11/3, and (ii) non-negligible ion bulk velocity fluctuations with respect to the magnetic fluctuations.
We are able to predict and explain such steep spectra with a simple two-fluid model for the turbulent cascade in low-beta plasma that does not include corrections due to intermittency of structures, instabilities, and/or waveparticle resonance, and which is in agreement with both in-situ spacecraft data and hybrid-kinetic simulations.
The dynamics in this regime is driven by the current and therefore by the magnetic field through the constraint expressed by Eq. 1. This is consistent with the fact that in almost all the time intervals exhibiting the −11/3 power-law strong gradients in the magnetic field components are observed, including reversals of the radial component (also called jets or switchbacks). The only two exceptions are intervals #29 and #30, which are characterized by an almost constant radial component. The main requirement, and the key difference with respect to alternative models, is that the ion kinetic energy cannot be neglected in the total energy density, although the ion motion only provides a very small contribution to the total current. This regime is therefore very different from previous theoretical models predicting a spectral index of −11/3, which typically require one of the following conditions: (i) the kinetic energy overtakes the magnetic energy [34] , (ii) the current is completely carried by the ions [32] , or (iii) the ion temperature is much larger that the electron temperature ("inertial kinetic Alfvèn waves") [29] . In our theoretical model none of these represents a necessary condition nor is observed to be satisfied in the numerical simulations presented here.
Additional simulations with varying physical parameters, currently under investigation, suggest that the combination of the ion plasma beta and the level of magnetic fluctuations with respect to the ambient field play a fundamental role in determining favorable conditions for the current-mediated regime to occur. Both a small beta and a large turbulence amplitude can determine an extended contribution of the "MHD term" u i × B to the electric field towards smaller scales below d i and, as a consequence, to the magnetic field evolution through Faraday's law. Firstly, in a low-beta plasma the ion bulk velocity fluctuations are observed to remain coupled to the magnetic field up to larger values of k (e.g., see Fig.  1 of Hellinger et al. [35] ). Secondly, if |δB| ∼ |B 0 |, then the term δu i × δB, which is often negligible with respect to δu i × B 0 , will also contribute in Faraday's law. Moreover, it is interesting to note that the contribution from the MHD turbulent cascade flux to the total cascade rate, computed via the third-order structure function, is observed to extend to smaller scales for smaller values of the ion plasma beta (see Fig. 3 of Hellinger et al. [35] ).
The importance of the small plasma beta in determining a steeper magnetic field spectrum at sub-ion scales has already been hinted in [36] . There, simulations with very small values (β i = 1/32 and 1/100), exhibit a spectral index compatible with −11/3 (see the top panel of Fig. 2 ), despite the much lower level of turbulent fluctuations with respect to the cases presented here (B rms /B 0 ∼ 1/8 and 1/16, respectively). This suggests that the large-amplitude condition may be necessary only when the plasma beta approaches unity and not for low-beta plasmas.
Our speculation on the role of the plasma beta and the level of fluctuations is also compatible with previous spacecraft observations by Bruno et al. [e.g., 37 ]. There, a steepening of the magnetic field spectrum just below d i , with a spectral index compatible with −11/3, has been observed in the solar wind, in intervals with lower values of the plasma beta and with a higher power of the fluctuations. This might suggest that the current-mediated regime can be important in other environments rather than just close to the Sun, as long as such conditions are met. The fact that these are more typical of the near-Sun environment, however, might explain why the current-mediated regime can be observed more clearly by PSP and it has not been specifically investigated before in observations. Both in PSP observations at its first perihelion and in [37] , the steep part of the magnetic field spectrum seems to be limited to the range of scales close to d i , i.e., slightly less than a decade in frequency. This is then followed by a flatter part, whose spectral index is more compatible with standard models which predict −7/3 [20, 31] or −8/3 [21] . Our model provides a simple explanation for this. If the ion velocity fluctuations are not negligible with respect to the magnetic fluctuations at the ion-scale transition, the −11/3 regime occurs. If, however, at a smaller scale their contribution to the plasma dynamics become negligible, we recover E k ∼ ✓ ✓ u 2 i + B 2 k and therefore the −7/3 spectral index. A last but not least important aspect to stress is that the current-mediated regime we propose is not a "kinetic" regime, in the sense that wave-particle resonance effects, such as ion cyclotron or ion Landau damping, are not necessarily required in order to explain the steepening of the spectrum. Since the dynamic is driven by the current, dissipation is likely rather associated with the disruption of intense current sheets through magnetic reconnection. It is worth mentioning that other features of the simulations presented here (not shown) include enhanced reconnection, enhanced proton heating, and enhanced intermittency. At present, we have not investigated the presence of ion cyclotron waves, which have been detected in PSP observations at its first perihelion [15, 16] . These are preferentially observed during time intervals characterized by a radial mean magnetic field and cause a bump in the magnetic field spectrum at frequencies corresponding to the ion scales. This is consistent with the results of our analysis: in correspondence with the intervals with a quasi-radial field (grey shaded areas in Fig. 1b) we do not detect clear power laws in the magnetic field spectrum and this is typically very steep, (see grey dots in Fig. 4a ) compatible with the drop that follows a bump around the ion scales. This further supports the idea that ion cyclotron resonances are unlikely to be responsible for the spectrum with a spectral index of −11/3. Future investigation of the above mentioned aspects will provide a fundamental insight for understanding the role of the energy cascade and current structures on energy dissipation and particle heating in turbulent plasmas. In this regard, it is important to note that we can expect the current-mediated regime to be relevant not only for space plasmas, such as the solar wind and the solar corona, but also for a broad class of other low-beta plasma environments. These can include astrophysical plasmas, e.g., in non-relativistic magnetized jets [38] and disks [39] , as well as laboratory plasmas, e.g. in magnetic confinement devices if the Hall term in Ohm's law is not negligible [40] , and possibly in turbulent dynamo experiments induced by laser-plasma interactions [41, 42] .
The numerical dataset was produced using the highperformance 3D hybrid particle-in-cell code CAMELIA [24] on the Peta4 machine at the Cambridge Service for Data Driven Discovery (CSD3), within the UK DiRAC infrastructure. Additional preliminary and complementary simulations (not shown here but important as a further confirmation of the results) were also performed on the Theta HPC system at the Argonne Leadership Computing Facility and on the Marconi HPC system at the Italian SuperComputing centre CINECA.
CAMELIA is based on the Current Advance Method and Cyclic Leapfrog (CAM-CL) code of Matthews [43] and solves the Vlasov-Maxwell equations, comprising the equations of motions for individual ions, and the electron fluid equations. The ions are modelled as macroparticles, i.e., fluid elements representing portions of the particles velocity distribution functions in phase space, and are tracked in a Lagrangian frame, whereas moments of the distribution such as densities and currents are computed simultaneously on Eulerian (stationary) mesh points. Electrons act as a massless, charge-neutralizing fluid. The magnetic field is expressed in units of the magnitude of the ambient magnetic field, B 0 , the ion and electron bulk velocities are in units of the Alfvén velocity, V A , time is in units of inverse ion gyrofrequency, Ω −1 i , and lengths are in units of the ion inertial length,
The initial condition consists of a homogeneous plasma embedded in a uniform ambient magnetic field, B 0 = B 0ẑ , which in the 2D case is out of the (x, y) simulation plane. We perturb the plasma with Alfvénic-like magnetic and ion bulk velocity fluctuations which are perpendicular to B 0 . These are the superposition of modes with equal amplitude, random phases, energy equipartition between magnetic and kinetic energy (i.e., with negligible initial residual energy) and they are balanced (i.e., with zero initial cross-helicity). Such initialization does not provide a realistic and comprehensive representation of the environment met by PSP at its first perihelion. We try, however, to reproduce similar plasma conditions by setting the ion and electron plasma betas to values compatible with the their average observed value, i.e., β i = 0.2 and β e = 0.5. The reader can find the detailed time variation of the electron core beta close to the first PSP perihelion in [44] . The other parameters characterizing the physical initial condition for the 2D and 3D simulations are, respectively: spectrum of initial magnetic fluctuations proportional to k and k 2 at k inj ⊥ d i 0.4 and k inj d i 0.5, amplitude of the initial magnetic fluctuations B rms /B 0 ∼ 0.44 and B rms /B 0 ∼ 0.33. The numerical setting for the 2D and 3D simulations consists of, respectively: 4096 2 and 512 3 grid points, box size L x,y = 256 d i and L x,y,z = 64 d i , spatial resolution δ x,y = d i /16 and δ x,y,z = d i /8, 1024 and 512 particle-percell, and resistivity η = 1.5 × 10 −3 and η = 5 × 10 −3 in units of 4πV A c −1 Ω −1 i . All parameters were chosen on the basis of our previous studies, which recovered a qualitative agreement with observations in the solar wind [e.g., 36, 45] and a quantitative agreement with observation at the Earth's magnetopause [25] . Futher details about our numerical setup and its implementation can be found in Franci et al. [46] and [26] for the 2D and the 3D simulations, respectively. Simulation data are here analysed and shown at t = 50 Ω −1 i , i.e., the time when the RMS value of the current density, J, reached its maximum. This is universally regarded as a proxy for the development of a fully-developed quasi-stationary turbulent state. We have verified that the spectral properties do not change significantly around that time. The time evolution of global quantities (e.g., the rms values of J and of the magnetic fluctuations, B, not shown here) is qualitatively the same as in previous 2D [46] and 3D [24] hybridkinetic simulations, the quantitative difference being the time needed for the maximum turbulent activity to develop. This is inversely proportional to the initial rms of fluctuations, B rms , and to the injection-scale wavenumber, k inj , through the non-linear time τ nl ∼ (k inj B rms ) −1 . This work represents a further evidence that a 2D simulation may be representative enough of the spectral behavior in a more realistic 3D geometry, due to the quasi-2D nature of the turbulent cascade [e.g., 24, 25] . Based on this, we only show and analyse here the spectra from the 2D simulation, which employs a much higher spatial resolutions and therefore exhibits more extended power laws.
The in-situ spacecraft data shown in Fig. 1 have been obtained from the FIELDS Fluxgate Magnetometer (MAG) onboard Parker-Solar-Probe at its first perihelion (from 00:00 UT to 23:59 UT on 6 November 2018). The data are in the RTN (radial-tangential-normal) reference frame, which is centered at the spacecraft and oriented with respect to the line connecting it to the Sun: the R axis is directed radially away from the Sun, the T axis is the cross product of the Sun's spin vector and the R axis, and the N axis completes the right handed set. The time series for the three components of the magnetic field have been divided in 32 subintervals of 45mins. For each time interval, the time series have been normalized to the amplitude of their respective mean vector magnetic field. The corresponding multitaper power spectral densities have been computed using 16 tapers.
For the simulations, we define the 1D omnidirectional power spectrum of a quantity Ψ(x) as
where δΨ(k) is the amplitude of the fluctuations of Ψ at the scale corresponding to k. In the main text, we will omit both the δ and the hat symbol and express the fluctuations δΨ 2 (k) just as Ψ k .
In the 2D case, this is obtained by integrating over circles in the 2D plane (k x , k y ) delimited by k = |k ⊥ | = k 2 x + k 2 y and k + dk. Due to the numerical discretization of the grid, the integral is replaced with the sum
whereΨ are the (complex) Fourier coefficients.
In the 3D case, this is obtained by integrating over spherical shells in the 3D space (k x , k y , k z ) delimited by k = |k| = k 2 x + k 2 y + k 2 z and k + dk, which corresponds to P Ψ (k) = √ k 2 x +k 2 y +k 2 z =k |Ψ| 2 (k x , k y , k z ).
For observations, the power spectrum of a quantity Ψ(t) is
where δΨ(f ) is the amplitude of the fluctuations of Ψ at the frequency f . The local spectral index, α, is computing by fitting magnetic field spectrum over a sliding window of half a decade in wavenumber for the simulations and in frequency for observations. Once this is computed in the whole range of scale, its value at sub-ion scales is obtained through an automatic fitting procedure, where the fit interval is chosen automatically including the minimum value and all the surrounding points within 10% from that value. This automatic fitting procedure is applied to the PSP spectra of all 32 time intervals, to the 24h-averaged PSP spectrum, and to the 2D and the 3D numerical spectra.
The observational data used in this study are available at the NASA Space Physics Data Facility (SPDF), https://spdf.gsfc.nasa.gov/index.html. The numerical data will be uploaded to EUDAT (https://b2share.eudat.eu/) and made publicly available before the manuscript is published, as soon as the title and the journal are definitive. This will allow us to include those details in the description of the dataset. This section will be then updated accordingly, including the corresponding doi of the dataset. 