Evidence of entropy cascade in collisionless magnetized plasma turbulence

The turbulence of collisionless magnetized plasmas, as observed in space, astrophysical, and magnetically confined fusion plasmas, has attracted considerable interest for a long-time. The entropy cascade in collisionless magnetized plasmas is a theoretically proposed dynamics comparable to the Kolmogorov energy cascade in fluid turbulence. Here, we present evidence of an entropy cascade in laboratory plasmas by direct visualization of the entropy distribution in the phase space of turbulence in laboratory experiments. This measurement confirms the scaling laws predicted by the gyrokinetic theory with the dual self-similarity hypothesis, which reflects the interplay between the position and velocity of ions by perpendicular nonlinear phase mixing. This verification contributes to our understanding of turbulent heating in the solar corona, accretion disks, and magnetically confined fusion plasmas. Turbulence of collisionless magnetized plasmas is ubiquitous in space as well as laboratory plasmas, and as such is subject to intense study. The authors present experimental evidence of the existence of entropy cascade by direct visualization of entropy distribution in the phase-space of turbulence in laboratory experiments.

T he spectral universality in seemingly case-dependent fluid turbulence, which reflects self-similar dynamics in the inertial range of an energy cascade, is one of the most surprising findings in modern physics [1][2][3] . The presence or absence of universality in collisionless magnetized plasma turbulence exemplified by space, astrophysical, and magnetically confined fusion plasmas has been explored out of scientific interest. This question not only comes out of pure scientific curiosity but is also strongly associated with important practical problems, so-called turbulent heating in the solar wind 4,5 , solar corona 6,7 , accretion disks [8][9][10] , and energy transport in magnetically confined fusion plasmas [11][12][13] . These problems raise the question: at which scales are ions accelerated and thermalized? Some examples of turbulent heating include the coronal heating problem and the nonadiabatic temperature evolution of the solar wind 6 . One of the major models for explaining the heating mechanisms of the high-temperature corona/solar wind is the wave turbulence heating model (one is the microflare/nanoflare theory). The energy flux of the Alfvén waves/kinetic Alfvén waves propagating from the solar photosphere and chromosphere into the corona is considered to be sufficient to heat the corona temperature up to the observed temperature of~100 eV. However, this model encounters a complicated problem of what mechanisms (exemplified by ion cyclotron resonance, Alfvén resonance, linear phase mixing, and shock dissipation, among many others) are responsible for dissipating (partitioning) the wave energy of the Alfvén wave/kinetic Alfvén wave turbulence. The entropy cascade attributed to nonlinear phase mixing is a compelling hypothesis 4 for the mechanism. This hypothesis motivates us to experimentally study the validity of the entropy cascade via nonlinear phase mixing.
According to Boltzmann's H-theorem, entropy production, in other words, ion heating, is realized only by collisions in weakly collisional (collisionless) plasmas at the microscopic scale 14 much smaller than the ion gyro-scales (the inner scale) with a significant ∂f/∂v or ∂ 2 f/∂v 2 , where f and v are the ion velocity distribution function and the velocity vector of ions, respectively. However, macroscopic scales are typical energy inputs to the system (the outer scale), which are considerably separated from the dissipation scales (the inner scale). This fact strongly indicates that the injected energy is nonlinearly transferred as free energy (≈ negative signed entropy) from the outer scale to the inner scale in collisionless magnetized plasmas. This process is called an entropy cascade [14][15][16] . The gyrokinetic theory 15,17,18 predicts that the nonlinear phase-mixing 19 is responsible for the transfer of entropy in the inertial subrange 14,20 , and has been confirmed in a number of gyrokinetic numerical simulations regardless of twodimensional (2D) electrostatic 16 or three-dimensional (3D) electrostatic 21 or electromagnetic turbulence 4,22,23 . In addition, in a real plasma, an indication of a free energy cascade was observed in electromagnetic turbulence in the Earth's magnetosheath as fine-scale structures in the ion velocity distribution functions 24 . Regardless of the fact that the original concept of the entropy cascade by the nonlinear phase mixing considers only 2D electrostatic fluctuations 14,19 , signatures of the entropy cascade have been witnessed in magnetized plasma turbulence in various situations. This broad applicability of the entropy cascade picture to a wide range of collisionless magnetized plasma turbulence is attributed to the fact that turbulence of collisionless magnetized plasmas universally consists of highly anisotropic fluctuations with frequencies much lower than the ion cyclotron frequency in addition to the intrinsic electrostatic nature of perpendicular fluctuations.
Nonlinear phase mixing is a phase-randomizing mechanism of electrostatic fluctuations of quasi-2D magnetized plasmas by a finite-Larmor radius effect 19 . It couples ion dynamics in the position and the velocity spaces, resulting in entropy cascades in the phase space 14,25 .
Although the visualization of the entropy distribution in the phase space of ions has been implemented in some simulations 25,26 , no measurements have been conducted in real plasmas, including space and laboratory plasma experiments.
The gyrokinetic theory for treating gyroscale plasma phenomena considers kinetics of 'charged rings' of gyrating ions, whose statistics are represented by ring-averaged ion velocity distribution functions at a fixed guiding center. The gyrokinetic system is invariant under the scaling transformation; (g, ϕ, r, v ⊥ ) → (g, μ 2 ϕ, μr, μv ⊥ ), where g, ϕ, r, v ⊥ , and μ are fluctuation components of the ring-averaged ion velocity distribution function, the electrostatic potential, the space coordinate, the perpendicular velocity of ions, and a scaling factor, respectively 20 .
which is the perturbed distribution function, where F 0 is the background equilibrium distribution function that is steady on the timescale of turbulent fluctuations (a Maxwellian is assumed in the standard gyrokinetic theory), and R is the fixed guiding center position.
The dual self-similarity hypothesis for the position and velocity spaces based on the above scale invariance predicts the power laws of spectra of the free energy W g1 ¼ RR gðR;v ? Þ 2 2 dv ? dR and energy E g ¼ 1 2 Rφ ðrÞ 2 dr with the use of g(R, v ⊥ ) for the 2D electrostatic turbulence. The dual self-similar dynamics in the phase-space are expected to be accompanied by a dual cascade of two kinds of turbulence energy in 2D electrostatic turbulence: the forward cascade of free energy and an inverse cascade of electrostatic energy 25 . This dual self-similarity hypothesis extends the renowned Kolmogorov self-similarity hypothesis in fluid turbulence 27 to the phase space of ions in gyrokinetic turbulence. It argues that the turbulent behavior of ions at the sub-Larmor scales is independent of scale in both their position and velocity spaces. These are strongly coupled in the phase space by nonlinear phase mixing.
Therefore, in this experiment, we employed diagnostic instruments named Ring-averaged ion velocity distribution function probes (RIDFPs) 28 to measure g(R, v ⊥ ) for evaluation of the theoretically predicted power laws. Although our previous paper 29 showed agreement between measurements and the theoretical prediction on the power spectrum of E g as a function of the fluctuation wavenumber perpendicular to the background magnetic field k ⊥ , we still face a lack of information on the velocity space to conclude the validity of the dual-cascade theory of gyrokinetic turbulence.
Here, we show the direct visualization of the entropy distribution of ions in the phase space of laboratory-magnetized plasma turbulence. These results provide the evidence of the existence of dual self-similarity and the dual cascade of entropy and field energy in the phase space of turbulence.

Results and discussion
Experimental setup. In this experiment ("Methods": Magnetized plasma experiment (MPX) device, ring-averaged ion velocity distribution function probes, other measurement tools, and the plasma parameters), we prepared four states of electrostatic 2D turbulence driven by resistive drift-waves having varied driving perpendicular wavenumbers k drive ρ thi between~0 and 10 with control of the density of the background neutral particles, where ρ thi is the Larmor radius of ions at the thermal velocity v thi . While this control scheme of k drive ρ thi is not a well-established technique, we prepared the turbulent states with the control of the neutral particle density with reference to the numerical work 30 , which describes the numerical investigation results of an energy partition of turbulence in a cylindrical plasma. We use Fourier-Hankel components of the free energy: W g1 ðk ? ; pÞ p∑ jk ? j¼k ? jgðk ? ; pÞj 2 ¼ πk ? pjgðk ? ; pÞj 2 , where g(k ⊥ , p) is the Fourier-Hankel transform of g(R, v ⊥ ) concerning the guiding center position and the ion velocity 20 , and the parameter p indicates the reciprocal of the fluctuation scale of g(R, v ⊥ ) in the velocity direction corresponding to the wavenumber in the velocity space.
The gyrokinetic approximation requires ω ≪ Ω ci , where ω and Ω ci are the wave frequencies and the ion cyclotron frequency, respectively. The higher the axial magnetic field strength is, the more improved the gyrokinetic condition. However, our experiment does not completely satisfy ω ≪ Ω ci (ω/2π~0.1-4 kHz, Ω ci / 2π = 8.4 kHz), which is a consequence of the two following adverse impacts of the strong field on our plasma states: (i) the suppression of the excitation of the drift waves; in our plasmas, the high magnetic field stabilizes the drift waves because of a change in the radial profile of the density; and (ii) the reduction of the ion gyro-radius, which limits the performance of the diagnostic. In other words, the axial magnetic field strength is set at a weak value to fully use all the velocity channels of the RIDFPs by expanding the ion gyro radii. This improves the resolution of the RIDFP measurement in the velocity space. In the third paragraph from the end of this section, we explain the influence of this incomplete satisfaction of the gyrokinetic ordering on the interpretation of the experimental results.
Spectra of free energy in the k-v space. Figure 1 shows timeaveraged free energy spectra k ⊥ |g(k ⊥ , v ⊥ )| 2 as functions of k x ρ thi and v ⊥ up to~2v thi for ( Fig. 1a) k drive ρ thi~1 (W g1 injection at large scales) and (Fig. 1b) k drive ρ thi~1 0 (W g1 injection at small scales), respectively. We regard k x ≈ k ⊥ in evaluating k ⊥ |g(k ⊥ , v ⊥ )| 2 based on the assumption of isotropy and fast decay of |g(k ⊥ , v ⊥ )| 2 for k ⊥ . The factor k ⊥ indicates the circumferential length in the k ⊥ -space.
k ⊥ |g(k ⊥ , v ⊥ )| 2 for (A) k drive ρ thi~1 exhibits, in the entire range of v ⊥ , the large amplitude at k x ρ thi~1 −2 consistent with jφ f ð f ; k y Þj 2 measured by a fine-scale Langmuir probe (FSLP) for k drive ρ thi~1 (Fig. 2a), which shows excitation of k x ρ thi~1 components at f~1.0−1.3 kHz. Accordingly, the low k ⊥ components at k x ρ thi~1 -2 show variation in the v ⊥ -direction at the scale of~v thi in Fig. 1a. The amplitude of components for |k x ρ thi | > 2 decays in large k x with finer-scale variation in the v ⊥direction. On the other hand, k ⊥ |g(k ⊥ , v ⊥ )| 2 for (B) k drive ρ thi~1 0 exhibits significant power at |k x ρ thi | > 10. These |k x ρ thi | > 10 components exhibit structures in the v ⊥ -direction at sub-v thi scales. This indicates strong coupling between the position and velocity spaces. See the next subsection for a quantitative assessment of scales in v ⊥ space using spectra.
Comparison with theory. In the second-fourth rows of Fig. 2, we show variation in the spectra of the collisionless invariants of 2D electrostatic gyrokinetics, W g1 (k ⊥ , p) and E g (k ⊥ ), for varied k drive ρ thi approximately between 0 and 10, (the first column) k drive ρ thi~1 , (the second column) k drive ρ thi~3 −4, (the third column) k drive ρ thi~0 -4, and (the fourth column) k drive ρ thi~1 0. These values of k drive ρ thi are evaluated from jφ f ð f ; k y Þj 2 measured by FSLP as shown in the first row of Fig. 2a-d. In the respective cases of different k drive ρ thi , the following two distinct frequency bands can be identified in jφ f ð f ; k y Þj 2 : one showing significant power with nonstatistical/deterministic distribution corresponding to the drive, for example, k y ρ thi~1 component at f~1.3 kHz in Fig. 2a (k y ρ thi~3 -4 component at f~1.7 kHz in Fig. 2b, and k y ρ thi~1 0 component at f~3-4 kHz in Fig. 2d), and the other bands showing scattered, statistical distribution in the k y -direction corresponding to turbulent cascade. Nonlinear dynamics exhibited in these shots are local because no nonlocal resonant wave-wave interactions occurred except for the k drive ρ thi~0 -4 state (the vertical array of Fig. 2c, g, k and o), which exhibits a continuous driving spectrum in the f-k y ρ thi diagram and is a turbulent state originating from a nonlinear coherent state of drift waves called solitary drift waves 31 .
The plots shown in the second row ( Fig. 2e-h) are W g1 (k ⊥ , p) (negative signed entropy) in the k ⊥ -p space, namely, position-velocity space (see "Evaluation of 1D spectrum of free energy W g1 1D (k x , p) and relationship with W g1 (k ⊥ ) and W g1 (p)" in Methods). In the following results of varied k drive ρ thi , a position of the drive moves along the diagonal in the k ⊥ -p space, as shown in Fig. 2e-h.
When W g1 is injected at large scales (Fig. 2e), in which a significant power of W g1 can be seen at the left bottom of the diagram at k ⊥ ρ thi~p v thi~1 -4 corresponding to k drive ρ thi , W g1 broadly develops along the diagonal from the driving scale toward the upper-right direction, reaching the top right corner around k ⊥ ρ thi~p v thi~2 0. (A similar profile can be seen in Fig. 6 in the article by Tatsuno et al. 26 ). According to the gyrokinetic Poisson Eq., the potential fluctuationφðk ? Þ is generated solely by nonlinear phase mixing, that is, gðk ? ; p ¼ k ? Þ (Eq. (2.11) in the article by Plunk et al. 20 ). This provides a Fjørtoft-type relationship W g1 (k ⊥ , p = k ⊥ ) = k ⊥ E g (k ⊥ ) (Eq. (7.14) in the article by Plunk et al. 20 ). Therefore, diagonal transport of W g1 in the k ⊥ -p space is inevitable when a dual-cascade of W g1 and E g occurs. Fig. 1 Spectra of free energy in the wavenumber-velocity space. Time-averaged power spectra k ⊥ |g(k ⊥ , v ⊥ )| 2 measured by RIDFPs as functions of the wavenumber k x normalized by the ion Larmor radius ρ thi at the thermal velocity v thi and the perpendicular ion velocity v ⊥ up to~2v thi for a driving perpendicular wavenumbers k drive ρ thi~1 and b k drive ρ thi~1 0, respectively. This W g1 (k ⊥ , p) spectrum exhibits the generation of W g1 at finer scales than the drive. This indicates that the nonlinear phase mixing works, resulting in a forward cascade of W g1 . This cascading nature will be examined below by comparing their power laws obtained from the measurement and the theory.
As k drive ρ thi increased to ≈10, as shown in Fig. 2f-h, the significant W g1 part moved in the upper-right direction along the diagonal line toward higher k ⊥ ρ thi~p v thi locations corresponding to the increase in k drive ρ thi . In all the cases of k drive ρ thi (Fig. 2e-h), the spread of W g1 (k ⊥ , p) along with the diagonal directions (between the bottom left and the upper right) is obvious. In these cases, the Fjørtoft-type energetic transition 25 , namely, a transition involving only diagonal components, can be realized.
In the third row of Fig. 2i-l), (crosses) one-dimensional (1D) spectra W g1 1D of the free energy as functions of k x ρ thi and (open squares) p-spectra of W g1 (p), respectively, are shown. Note that when W g1 (k ⊥ ) ∝ k ⊥ −α , W g1 1D (k x ) ∝ k x −α (see "Evaluation of 1D spectrum of free energy W g1 1D (k x , p) and relationship with W g1 (k ⊥ ) and W g1 (p)" in Methods). Respective k drive ρ thi ranges are shown by purple bars. The blue and red bars represent the power laws predicted by the gyrokinetic theory for (blue) the inverse-cascade range k ⊥ < k drive and (red) the forward cascade range k ⊥ > k drive , respectively 20 . The respective power indices of W g1 1D (k x ) and W g1 (p) evaluated from the least-square fitting to the measured data points are indicated by numbers with the standard errors (SE). The noise levels of W g1 1D (k x ) and W g1 (p) estimated from a vacuum shot are included in Fig. 2i. W g1 (p) in Fig. 2i (injection at a large scale) shows power-law decay starting at k x ρ thi~kdrive ρ thi~1 up to the detection limit of k x ρ thi~p v thi~2 0, clearly indicating the existence of an inertial subrange of the entropy (free energy). The evaluated exponent for the measured W g1 (p), −1.0, approximately agrees with the theoretical value of −4/3 within 5 times SE. For smaller scales k x ρ thi > 3, W g1 1D (k x ) in Fig. 2i shows a similar decaying tendency to W g1 (p), which is a strong indication of the dual self-similarity dynamics, although the discrepancy of W g1 1D (k x ) from the theoretical exponent is larger than 25 times SE. Larger discrepancies between the theoretical exponents and the Fig. 2 Variation in spectra of the invariants of gyro-kinetics. The free energy W g1 (k ⊥ , p) and the electrostatic energy E g (k ⊥ ) for varied driving perpendicular wavenumbers k drive ρ thi between~0 and 10, (the 1st column) k drive ρ thi~1 , (the 2nd column) k drive ρ thi~3 −4, (the 3rd column) k drive ρ thi~0 −4, and (the 4th column) k drive ρ thi~1 0, respectively. The plots in the 1st row (a-d) are the power spectra of potential fluctuation jφ f ðf; k x Þj 2 measured by FSLP.
The plots shown in the 2nd row (e−h) are W g1 (k ⊥ , p) in the k ⊥ -p-space, namely, position-velocity-space. The 3rd row (i−l) are (cross) 1D spectra W g1 1D of the free energy as functions of k x ρ thi and (open square) p-spectra of W g1 (p), respectively. E g (k ⊥ ) spectra are shown in the 4th row (m−p). Respective k drive ρ thi ranges are shown by purple color bars. The blue and red bars represent power laws predicted by the gyrokinetic theory 20 for (blue) k ⊥ < k drive and (red) k ⊥ > k drive , respectively. measured ones in W g1 1D (k x ) than those in W g1 (p) are seen in all cases of k drive . As k drive ρ thi increases to 3-4 (Fig. 2j), W g1 (p) shows an inflection point at k x ρ thi~kdrive ρ thi~3 -4. Correspondingly, W g1 1D (k x ) follows this change in the spectrum. For an injection at small scales k drive ρ thi~1 0 (Fig. 2l), both W g1 (p) and W g1 1D (k x ) become independent of p and k x , indicating good agreement with the theory 20 . It should be noted that the decaying exponents of both W g1 1D (k x ) and W g1 (p) at k x ρ thi ≫ 1 exhibit similar values in the respective cases of k drive ρ thi and similar variation for varied k drive ρ thi in Fig. 2i-l. This indicates the formation of inertial subranges in both the position and velocity spaces and strong coupling between them. A similar result is obtained in the gyrokinetic simulation 16,26 .
In the case of k drive ρ thi~0 -4 (Fig. 2k), despite the driving source in the target k-range (i.e., noninertial range), strong k-p coupling occurs, indicating the universality of nonlinear phase mixing in magnetized plasma turbulence.
The abovementioned approximate agreement between the measurement and the theoretical prediction of W g1 (p) and W g1 (k ⊥ ) within~1 − 25 times SE is supported by the E g (k ⊥ ) spectra measured by FSLP shown in the fourth row of Fig. 2  (Fig. 2m-p). Note that when E g (k ⊥ ) ∝ k ⊥ −α , E g1 1D (k x ) ∝ k x −α ("Methods": Evaluation of the 1D spectrum of the electrostatic field energy E g 1D (k y ) and its relationship with E g (k ⊥ )). The blue and red bars represent the power laws predicted by kinetic theory for (blue) k ⊥ < k drive and (red) k ⊥ > k drive , respectively. The respective power indices of E g1 1D (k x ) evaluated from the leastsquare fitting to the measured data points are indicated aside by numbers with the standard errors. The noise levels of E g1 1D (k x ) estimated from a vacuum shot are included in Fig. 2m-p. The evaluated scaling for E g1 1D (k x ) in Fig. 2m is k x −3.0 , close to the theoretical power law k ⊥ −10/3 . The theoretical prediction k ⊥ −2 for the inertial range in the larger scale than the driving scale becomes dominant as k drive ρ thi increases. The exponent evaluated for the measured E g1 1D (k x ) in Fig. 2p is k x −1.8 . The observed behavior of these power laws in W g1 and E g is consistent with the theory based on the hypothesis of dual self-similarity and the dual cascade of the entropy and field energy of the gyrokinetic system 20 . Our result strongly supports the validity of the dual-self similar hypothesis of gyrokinetic turbulence.
Note that the nonlinear decorrelation time 16,26 for these states is evaluated as τ ρ~1 0 × 10 −5 s; hence, the dimensionless parameter 16,26 , which represents the scale separation analogous to Reynolds number, D~30-40. Accordingly, the cutoff wavenumber, above which the collisional dissipation dominates the cascade, k cut ρ thi~2 D 3/5~2 0. This guarantees that the observed power-law spectral range corresponds to the inertial subrange.
In what follows, we discuss the applicability of the gyrokinetic theory and its impact on the interpretation of the experimental results. The theoretical aspects of the applicability of gyro kinetics to space and astrophysical plasmas are discussed in previous works 15,18,23 . In brief, strong magnetization, anisotropy, small amplitude, low fluctuation frequencies, equilibrium Maxwellian distribution, and nonrelativistic effects are assumed.
As described in Experimental setup, our experiment does not completely satisfy the strong magnetization condition ω ≪ Ω ci in consideration of the two adverse impacts of the strong field on our plasma states. One possible influence of this incomplete satisfaction of the gyrokinetic ordering is the possibility of stochastic ion heating 32,33 . According to ref. 34 , the amplitude threshold for strong stochastic heating when β < 1 can be expressed in terms of the quantity ε i ∼ q i δΦ i /mv 2 , where β, q i , and δΦ i are the thermal pressure normalized by the magnetic field pressure, electric charge of the ions, and root mean square amplitude of the electrostatic potential fluctuation at k ⊥ ρ thi ∼ 1, respectively. In our case, ε i ∼ 0.2; therefore, the influence of stochastic ion heating is considered to be small. Ion cyclotron resonance heating is negligible because ω ⋡ Ω ci .
One of the weakest points of the gyrokinetic theory is the use of a Maxwellian as F 0 15 . Therefore, highly nonequilibrium turbulent states are not subjects of gyrokinetic theory. Nevertheless, the gyrokinetic theory is in effect in a wide range of space and turbulent astrophysical plasmas, as shown in numerous numerical simulations and theoretical considerations. Although our measurement has not had sufficient accuracy to be able to discuss a degree of deviation of F 0 from a Maxwellian distribution (Fig. 3d in "Methods" shows an example of the velocity distribution F 0 ), we believe that F 0 of our plasmas is not that different from a Maxwellian because of the following reasons: (i) no mechanisms that heat preferentially ions with specific velocity components exist; the collisional energy transfer caused by electron-ion collisions is the only possible heating source of ions; and (ii) ions are weak collisional, but not too weak. The meanfree path of ions for ion-ion collision is~0.5 m < 1.6 m~L, which is the length of the plasma column. For the same reasons, pressure anisotropy is not considered.
The influence of the electron behavior on the entropy cascade in the ion phase space is considered negligible according to gyrokinetic numerical simulations 4,16,23 . The work done by Tatsuno et al. showed that no essential difference was observed between the cases of applying Boltzmann-response (3D) and noresponse (2D) electrons in their simulations 4 . For more detailed electron energetics, please see, the work by Kawazura et al. 23 , which studied the partition of irreversible heating between ions and electrons in terms of a compressive drive using hybrid gyrokinetic simulations.

Conclusions.
In laboratory experiments, we present the measurement results of the Gibbs entropy distribution in the phase space of ions in electrostatic gyrokinetic turbulence. We confirm good agreement between the measurement and the theoretical prediction with a dual self-similar hypothesis in the scaling laws of the spectra of the entropy and the electric field energy. This result indicates an interplay between the position and velocity spaces of ions by perpendicular nonlinear phase mixing accompanied by a cascade of entropy in the phase space. Excitation of spectral components of the free energy W g1 with k ⊥ ∼ p in the phase space gives evidence of nonlinear phase mixing at sub-Larmor scales. This observation corroborates the evidence of the existence of dual self-similarity and the dual-cascade of entropy and field energy.
An entropy cascade via nonlinear phase mixing is considered universal in magnetized plasma turbulence because it has been witnessed in various numerical simulations in different magnetized turbulence setups, including 2D 16 and 3D electrostatic 21 , 3D electromagnetic turbulence, including the Alfvenic solar wind 4,22 , and a simulation targeting an Alfvenic hot accretion disk 23 . The essential physical components in entropy cascades via nonlinear phase mixing include the gyro motion of the particles at the sub-Larmor scales and the interaction between the particles and the electrostatic fluctuations in addition to the strong anisotropy (k ⊥ ≫ k || ). Our experiment prepares turbulent states with these essential components in a simplified setup. Therefore, although our turbulence setup is in a narrow range (electrostatic turbulence in a drift-kinetic regime), our result is crucial for verifying universality. We believe that verifying the entropy cascades via nonlinear phase mixing contributes to understanding turbulent heating in the solar corona, solar wind, accretion disks, and magnetically confined fusion plasmas.

Methods
Magnetized plasma experiment (MPX) device: experimental device. The experiment was conducted in the MPX device 35 (Fig. 3a), which can prepare qualified plasma configurations for our purposes, namely, quasi-2D configurations with the magnetic field (the strength is denoted by B 0 ), the density, and temperatures within the target ranges. The hot cathode generated plasmas with the introduction of argon gas as the working gas. g(R, v ⊥ ) was measured by two RIDFPs separated by 3 mm and 230 mm in the x-and the z-directions, respectively (Fig. 3a), to obtain perpendicular wavenumbers k ⊥ of g(R, v ⊥ ) with the use of a two-point correlation technique 36 . RIDFPs were positioned at the center axis of the target plasmas to minimize the influence of steady E × B rotation of the plasmas. Figure 3 shows the equilibrium ring-averaged ion velocity distribution functions F 0 for B 0 = 0.022 T and B 0 = 0.033 T, respectively. The red solid curves indicate Maxwellian with a temperature of 0.1 eV. It can be observed that the orbits of the ions (Larmor radii) shrink as B 0 change from 0.022 T to 0.033 T, indicating the validity of the RIDFP measurement.
Ring-averaged ion velocity distribution function probes 28 . The RIDFP is a set of ion collectors that detect different velocity components and is immersed into magnetized plasmas. The RIDFP achieves momentum selection of incoming ions by selection of the ion Larmor radii by the orbit filter. To nullify the influence of the sheath potential surrounding the RIDFP on the incoming ions' orbits, the RIDFP body's electrostatic potential is automatically adjusted to coincide with the space potential of the target plasma with the use of an emissive probe and a voltage follower 28 .
For precise measurement of the fluctuation g(R, v ⊥ ) in time, the ground potential of the current detection circuits is adjusted to that of the RIDFP body in order to eliminate capacitive coupling between the ion collectors and the RIDFP body. The length of RIDFP L = 79 mm satisfies the condition of v || f ci −1 > L, where v || is the parallel component of the thermal velocity of ions. This inequation imposes a condition for the plasma (ions) residing in the magnetic flux tube, in which an RIDFP is allocated, not to being severed by the RIDFP.
In the experiments, the resolution of the RIDFP measurement in the p-space Δ(pv thi ) can be evaluated as follows: Δpv upper_bound~Δ p(2v thi ) = j 0 ≒ 2.4, ∴Δp ≒ 1.2/ v thi ∴Δ(pv thi )~1.2, where v upper_bound and j 0 are the upper bounds of the measurable ion velocity by RIDFP and the first zero of the Bessel function of the first kind, respectively. The maximum (pv thi ) in the RIDFP measurement is obtained as: (pv thi ) max~j0 v thi /Δv~20.
The other measurement tools and the plasma parameters. In addition to RIDFPs, we employed a few sets of Langmuir probes (LPs) for different purposes. An LP measured the plasma density in addition to a microwave interferometer 37 .
The LP also measured the electron temperature and the space potential of the plasmas. A poloidal array of LPs named LPA diagnosed the macroscopic poloidal structures of fluctuations of ion saturation current and their floating potential. Another set of LPs named fine-scale LP (FSLP) was applied to measure the fluctuation of floating potentialφ f %φ as functions of frequencies f and k y (k x ), whose upper bound of measurable k ⊥ is~8 × 10 2 m −1 .
The hot cathode generated plasmas with the introduction of argon gas as the working gas. The typical plasma density and temperature are n e~0 .5-5 × 10 16 m −3 and T e~1 -5 eV, respectively, with a 15-20 cm diameter. Accordingly, the characteristic frequencies of the ion cyclotron f ci , drift-wave f DW , ion-ion Coulomb collisions f ii , and ion-neutral collisions f in have the following numbers: where Ω ci is the ion cyclotron angular frequency of singly ionized argon ion.

Iso
Evaluation of the 1D spectrum of the electrostatic field energy E g 1D (k y ) and its relationship with E g (k ⊥ ). From the FSLP measurement, the 1D spectrum of the electrostatic field energy E g 1D ðk y Þ R 1 À1 jφðk ? Þj 2 dk x is evaluated as follows: If turbulence is isotropic, E g 1D ðk y Þ ¼ R 1 À1 jφð ffiffiffiffiffiffiffiffiffi k 2 x þk 2 y p Þj 2 dk x : Applying a variable change tan θ ¼ k y k x and a cutoff θ cut for the upper bound of the integration, one can obtain When E g (k ⊥ ) ∝ k ⊥ −α , E g 1D (k y ) = ∫E g (k ⊥ )dk x ∝ k y −α .

Data availability
The data supporting this study's findings are available from the corresponding author upon reasonable request.
Received: 16 April 2022; Accepted: 9 December 2022; Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.