Kinetic Turbulence in Astrophysical Plasmas: Waves and/or Structures?

The question of the relative importance of coherent structures and waves has for a long time attracted a great deal of interest in astrophysical plasma turbulence research, with a more recent focus on kinetic scale dynamics. Here we utilize high-resolution observational and simulation data to investigate the nature of waves and structures emerging in a weakly collisional, turbulent kinetic plasma. Observational results are based on in situ solar wind measurements from the Cluster and MMS spacecraft, and the simulation results are obtained from an externally driven, three-dimensional fully kinetic simulation. Using a set of novel diagnostic measures we show that both the large-amplitude structures and the lower-amplitude background fluctuations preserve linear features of kinetic Alfven waves to order unity. This quantitative evidence suggests that the kinetic turbulence cannot be described as a mixture of mutually exclusive waves and structures but may instead be pictured as an ensemble of localized, anisotropic wave packets or"eddies"of varying amplitudes, which preserve certain linear wave properties during their nonlinear evolution.

A series of compelling questions emerge in the above context such as, which of the two aspects, structures or waves, is more essential and/or fundamental to the turbulent dynamics? And also, are waves and structures mostly exclusive to each other and simply "coexist" in the turbulence or is there a deeper connection between the two? Here, we address these open questions in the context of kinetic-scale plasma turbulence, where the topic of waves and structures is hotly debated at present [15,23,25,26,[28][29][30][31][32][33][34][35][36]. In the following, we explain the astrophysical background and motivation for the present study. Our general approach and the qualitative conclusions drawn from it are, however, not exclusively limited to kinetic-scale plasma turbulence and could be relevant to a broad range of turbulent regimes and environments where waves and structures can be identified.

A. Astrophysical background and motivation
In magnetized astrophysical and space plasmas, the role of waves in turbulence has been widely debated since the early days of magnetohydrodynamic (MHD) turbulence research when the first models based on nonlinearly interacting Alfvén waves were proposed [37,38]. Soon after, in situ spacecraft measurements detected Alfvén wave signatures in the turbulent solar wind [39]. Observations and models were followed by pioneering numerical studies (see, e.g., Refs. [40][41][42]). These simulations revealed the development of intense electric current sheets, which are now believed to play a key role in the (Ohmic) dissipation of MHD turbulence [20,21,43]. It also soon became clear that magnetized plasma turbulence is inherently anisotropic, as evidenced by the preferential cascade of energy to small perpendicular scales with respect to the direction of the mean magnetic field [44][45][46].
Motivated by the observational and numerical evidence available at the time, Goldreich and Sridhar [47] devised their celebrated model of anisotropic, strong MHD turbulence based on the critical balance conjecture. In short, they assumed a balance between the linear (wave-crossing) and nonlinear (eddy turnover) time at each scale. This allows for a phenomenological prediction of the turbulence energy spectrum and its anisotropy. The latter is presently consistent with numerical simulations [48][49][50] and recent observations [9,[51][52][53]. Although the Goldreich and Sridhar model is by no means a rigorous description of MHD-scale turbulence, their phenomenological approach received a great deal of interest in the community and it later paved the way for many notable works on astrophysical plasma turbulence (see, e.g., Refs. [54][55][56][57]).
It is worth noting that the critical balance conjecture [47] is essentially a statement about the persistence of linear wave physics in a strongly turbulent wave system. The rather generic nature of the statement implies that critical balance may be relevant in a more general context outside the scope of MHD turbulence. Indeed, the concept has been more recently adapted in models of subion-scale plasma turbulence [56,57], and it even received attention in fields outside of plasma physics [4,6,12]. In particular, Ref. [12] proposed critical balance as a universal scaling conjecture for strong turbulence in wave systems, with application to MHD, rotating flows, and stratified flows.
The discussion on waves versus structures became more pronounced in the kinetic plasma turbulence context, since the two viewpoints lead to fundamentally different interpretations of how turbulent dissipation works. In a wavemediated type of turbulence, there is hope that at least a significant fraction of heating can be estimated from linear (collisionless) wave damping rates [60,62,[91][92][93][94], whereas a structure-mediated regime might first require an understanding of how structures form and dissipate energy during their evolution [25,28,32,33,80,81]. It is important to note that the two viewpoints are not a priori incompatible, if the kinetic-scale structures themselves were to preserve certain wavelike features. The latter possibility is explored in the present paper.

B. Our approach and relation to previous works
We employ the KAW turbulence phenomenology [56,57,95] as the main basis for making contact with linear wave predictions. It is not a priori obvious that the subionscale fluctuations are of the KAW type. Other possibilities, such as whistler waves [75,77,78], kinetic slow waves [96], or ion Bernstein modes [97], have been considered as well. Our choice is, however, favored by the presently available empirical evidence (see Podesta [98] for a detailed review). A significant portion of evidence is based on the measured statistical polarization properties of the turbulent fields [29,31,36,69,70,[99][100][101][102][103][104], which have been compared against linear KAW predictions. This includes, for instance, the measured spectral ratio between the parallel and perpendicular magnetic field component (see, e.g., Refs. [29,101,102]) and the spectral ratio between the electron density and the magnetic field [31]. Some works considered as well the magnetic helicity and related its observed spectral signatures to the right-hand elliptical polarization of the KAW [99,100,103]. Apart from the above, a number of computational and observational studies suggest that the subion-scale fluctuations are predominantly low frequency [91,105] and strongly anisotropic [36,[106][107][108], as qualitatively expected for KAW turbulence [56,57,88,95].
Previous works supporting wavelike interpretations of kinetic-scale plasma turbulence were often focused on turbulence spectral features, while paying relatively little attention to the characteristic structure of the fluctuations in real space and their intermittent statistics [21,23,28,33,80,89,109,110]. In order to account for turbulent structure formation, it is necessary to go beyond a standard spectral analysis. To this end, we introduce here appropriate generalizations of the frequently employed spectral field ratios. Unlike the standard field ratios, the generalized definitions are specifically designed to probe the statistical field polarizations within the large-amplitude turbulent structures. This is accomplished with the aid of a wavelet scale decomposition [23,34,102,111,112], which can be used to analyze a signal simultaneously in real space and spectral space.
Our newly defined diagnostic measures are applied jointly to high-resolution observational and kinetic simulation data. Observational results are based on in situ measurements from the Cluster and Magnetospheric Multiscale (MMS) spacecraft, whereas the numerical results are obtained from an externally driven, threedimensional fully kinetic particle simulation. The majority of previous kinetic simulations were carried out either in a two-dimensional geometry and/or with a reduced-kinetic model. The validity of both types of approximations has been questioned [72,89,95,[113][114][115]. Specifically for KAW turbulence, pioneering 3D kinetic simulations were carried out using the gyrokinetic formalism [30,60,92,116]. These works obtained results consistent with theoretical expectations. However, gyrokinetics assumes lowfrequency (compared to the ion cyclotron motion), strongly anisotropic (with respect to the mean field), and smallamplitude fluctuations [57,117], which naturally favors a KAW-type of turbulent cascade. Only recently, various aspects of KAW turbulence were studied in terms of 3D hybrid-kinetic (i.e., with fully kinetic ions and fluid electrons) and fully kinetic simulations [36,104,[118][119][120].
The rest of this paper is organized as follows. In Sec. II, we provide a short theoretical background. The methods employed are described in Sec. III. This includes a brief overview of our 3D fully kinetic simulation, a description of observational data, a summary of wavelet scale decompositions, and finally, the generalized field ratio definitions. Our main results are presented in Sec. IV. We conclude the paper with a discussion (Sec. V) and summary (Sec. VI) of our results.

II. THEORETICAL BACKGROUND
The essential plasma dynamics underlying kinetic Alfvén turbulence may be approximately described with a set of ideal (i.e., dissipation-free [121]) fluid equations [57,95,122] for the electron density n e and flux function ψ, defined via b ⊥ ¼ẑ × ∇ ⊥ ψ, where ∇ ⊥ ¼ ð∂ x ; ∂ y Þ and b ⊥ is the part of the perturbed magnetic field perpendicular to the mean field B 0 ¼ B 0ẑ . It is convenient to first normalize the density and magnetic field according to [31,95,122] where b k ¼ b kẑ is the parallel perturbed magnetic field, n 0 is the mean density, are the ion and electron betas, and T i and T e are the ion and electron temperatures (measured in energy units), respectively. From here on, we drop the prime signs, but it is to be understood that all fields are normalized according to Eqs. (1)-(3) unless noted otherwise. In appropriately chosen time and length units (see Ref. [95]), the fluidlike equations read Parallel magnetic fluctuations are not evolved explicitly but are instead determined by a linearized pressure balance condition: where n 0 is now the mean electron density in the normalized units. Note that Eqs. (4)-(6) are obtained for a β ∼ 1 plasma in the asymptotic limit: where k k is the parallel wave number, k ⊥ is the perpendicular wave number, ρ i is the ion thermal gyroradius, and ρ e is the electron thermal gyroradius. The above model is essentially nonlinear. By dropping the second term on the right-hand side of Eqs. (4) and (5), a linear wave system is obtained, the nonstationary solutions of which may be written as linear combinations of plane KAWs with polarizations [57,95]: for each mode with a wave vector k, where E KAW;k ¼ jb ⊥k j 2 þ jn ek j 2 is the KAW energy density. Adopting the convention k k ≥ 0, the frequency of the waves is ω k ¼ AEk k k ⊥ . Similarly to inertial waves in rotating fluids [10,12,123], KAWs are dispersive and their group velocity is preferentially oriented along the mean field (assuming k k ≪ k ⊥ ). Besides KAWs, solutions of the linear system are also static structures with k k ¼ 0.
Following previous works [29,31,36,101], we use relation (9) as the basis for identifying kinetic Alfvén wave features in a nonlinear, strongly turbulent regime [124]. Note that while Eqs. (8) and (9) relate the Fourier coefficients of a single KAW, they do not strictly guarantee an analogous relation in real space for a linear combination of such waves. An exception to the latter is the relation between n e and b k , which is given at any point in space by Eq. (6). The (linearized) pressure balance is, however, a quite general signature of low-frequency dynamics not exclusively limited to KAWs. For instance, it may include mirror structures [23,[126][127][128] or kinetic slow modes [96]. A less ambiguous method to identify KAW properties is to compare the relative spectral amplitudes between jn ek j 2 and jb ⊥k j 2 to the corresponding theoretical prediction, which is more exclusive to KAWs [31], and it is derived from the nonlinear equations by specifically assuming a (monochromatic) wave solution. If pressure balance is taken for granted, the latter is equivalent to measuring the relative amplitudes between jb kk j 2 and jb ⊥k j 2 . While pressure balance alone does not necessarily imply wave activity, it is still a fairly stringent condition at kinetic scales since it rules out high-frequency fluctuations such as whistler waves [75,77,78,108], which cannot be a priori neglected. Moreover, the KAW and whistler wave predictions of jb kk j 2 =jb ⊥k j 2 may become degenerate for β i ≳ 1 [101].
It is straightforward to show [57] that a single monochromatic KAW is an exact solution of the fully nonlinear system (4) and (5). But that is not all. In analogy with shear Alfvén waves in MHD [55] and inertial waves in rotating fluids [12,129], a few special linear combinations of plane KAWs are as well exact solutions. The exact wave solutions are generally found by requiring that the nonlinear terms vanish, which is satisfied whenever the contours of n e , ψ, and ∇ 2 ⊥ ψ are aligned in every perpendicular plane. Unless all waves share the same k ⊥ , the alignment between the contours of ψ and ∇ 2 ⊥ ψ restricts the geometry of the solutions. It then follows that the corresponding perpendicular profiles of n e and ψ are either (i) circularly symmetric or (ii) one-dimensional (i.e., with a spatial variation along a single perpendicular direction). In an unbounded domain, class (i) may be formally constructed as a Fourier integral of plane KAWs with cylindrically symmetric coefficients (with respect to the z axis), whereas (ii) corresponds to any combination of KAWs with a fixed direction of the perpendicular wave vector k ⊥ . The exact wave solutions of type (i) and (ii) differ from the special class of waves with a fixed magnitude of k ⊥ found in Ref. [57] in that they may be composed of counterpropagating plane KAWs. Moreover, the absence of a fixed k ⊥ constraint permits wave packets with a spatially localized envelope. Strictly speaking, even if the perpendicular envelope of a KAW is spatially localized at a given time, it will not remain well localized at later times. However, since KAWs disperse primarily in the parallel direction, the perpendicular spreading of localized wave packets is relatively slow. In conclusion, the above aspects suggest that KAWs are deeply rooted in the magnetized plasma dynamics, beyond the standard limits of a linear approximation.
The vanishing of nonlinear terms for solutions with either circularly symmetric or 1D perpendicular profiles is in fact a quite generic feature of turbulent wave systems that possess a strong mean field, together with perpendicular nonlinearities in the form of a Poisson bracket: ff; gg 1 z · ð∇ ⊥ f × ∇ ⊥ gÞ. This includes the strongly anisotropic (k k ≪ k ⊥ ) limits of the MHD and rotating Navier-Stokes equations (see, e,g., Ref. [12]). For these two systems, the same types of exact linear wave solutions [namely (i) and (ii)] may be found in the ideal (i.e., dissipation-free) regime. For rotating flows, the analogue of type (ii) solutions is derived explicitly in Ref. [129], even without assuming k k ≪ k ⊥ . Moreover, class (i) relates to the monopole (force-free) Alfvén vortex solution in MHD [130,131]. Because of these similarities, some of our results might be directly relevant to structure formation in strongly anisotropic MHD and rotating turbulence.

III. METHODS
Below, we describe the numerical simulation and observational data and the data analysis techniques and introduce the generalized field ratios. These aspects are essential for a complete understanding of this work. However, those interested only in the main outcome of the study may skip over to the results in Sec. IV.

A. Driven 3D fully kinetic simulation
We perform the simulation using the particle-in-cell code OSIRIS 3.0 [132,133]. The periodic domain size is ðL x ; L y ; L z Þ ¼ ð18.9; 18.9; 48.3Þd i , where d i ¼ ρ i = ffiffiffiffi β i p is the ion inertial length. The global mean field B 0 points in the z direction. The spatial resolution is ðN x ; N y ; N z Þ ¼ ð928; 928; 1920Þ with 150 particles per cell per species. Quadratic spline interpolation is used for the chargeconserving electric current deposit. A reduced ion-electron mass ratio of m i =m e ¼ 100 is used and the electron plasma to cyclotron frequency ratio is ω pe =Ω ce ¼ 2.86. To reduce particle noise, the output fields used for the analysis are short-time averaged over a window of duration where Ω ce is the electron cyclotron frequency. Following Ref. [134], the forcing is an adaptation of the Langevin antenna [135] and introduces a timedependent external electric current J ext . We apply J ext at wave numbers ð1; 0; AE1Þ, ð0; 1; AE1Þ, ð−1; 0; AE1Þ, and ð0; −1; AE1Þ in units of ð2π=L x ; 2π=L y ; 2π=L z Þ. The antenna current is divergence-free and drives lowfrequency Alfvénic fluctuations. To avoid a rapid transient response, we initialize the fluctuating field at t ¼ 0 to match the external current via ∇ is the Alfvén speed, and the decorrelation rate [135] is γ 0 ¼ 0.6ω 0 . The ion and electron velocity distributions at t ¼ 0 are isotropic Maxwellians with β i ≈ β e ≈ 0.5.
The approach towards a statistically steady state is depicted in Fig. 1. The simulation is run for about 2.24 Alfvén transit times, t A ¼ L z =v A , until the kinetic-scale spectra appear to be converged. Towards the end of the simulation, the mean ion and electron betas (based on their space-averaged local values) are β i ¼ 0.56 and β e ¼ 0.51, while the mean temperature anisotropy is T ⊥i =T ki ≈ 1.2 for the ions and T ⊥e =T ke ≈ 0.9 for electrons. The turbulence is strongly driven such that b rms ≈ L ⊥ =L z ≈ 0.4 towards the end of the simulation [ Fig. 1(a)], where b rms is the rootmean-square fluctuating field in units of B 0 . This corresponds to the critical balance (b rms ≈ k k =k ⊥ ) at the forced wave numbers.
The obtained 1D magnetic energy spectrum as a function of k ⊥ ≡ ðk 2 x þ k 2 y Þ 1=2 [ Fig. 1(b)] exhibits an approximate power law with a steepening of the spectral exponent close to electron scales. The spectrum is consistent with solar wind measurements, which typically show spectral exponents around −2.8 at subion scales, while steeper exponents are measured close to electron scales and beyond [71,105]. In the inset of Fig. 1(b) we show for reference the 1D k z spectrum. To further reduce contributions from particle noise in the subsequent analysis, we filter out the noisedominated modes with k z d i > 12 [119]. Employing the method of Cho and Lazarian [122], we consider in Fig. 2 the anisotropy relative to the local mean field. The subion range anisotropy k k ðk ⊥ Þ is scale dependent and has a slope of approximately 1=3 on the logarithmic graph [136], in agreement with a recent fully kinetic simulation of decaying turbulence [36]. According to the asymptotic KAW theory for a β ∼ 1 plasma [95], KAWs are expected to exist for wave numbers with k k =k ⊥ ≪ 1 and k k d i ≪ 1.
While these values are not asymptotically small in our simulation, the estimated anisotropy is within the range k k =k ⊥ < 1 and k k d i < 1 at subion scales. We also estimate the ratio of linear (KAW) to nonlinear timescales [ Fig. 2 [56,57,107]. To obtain the scale-dependent fluctuation δb ⊥;k ⊥ , we apply a bandpass filter between k ⊥ =2 and 2k ⊥ on b ⊥ and compute its root-mean-square value [36]. Note that τ L and τ NL are based here on Eqs. (4) and (5), but their ratio χ takes the same form as the nonlinearity parameter in MHD if dynamic alignment is neglected [47,50,55]. The ratio of linear to nonlinear timescales is close to unity for k ⊥ ρ i ≲ 1 and exhibits a slight decline throughout the subion range. This simple yet direct estimate of χ contrasts some previous works [14] and implies that the kinetic-scale nonlinear effects are not any more significant than linear wave physics, despite the fact that the turbulence fluctuation amplitudes are not small [ Fig. 1(a)]. These circumstances provide additional motivation for comparison with linear KAW predictions.
In order to facilitate a direct comparison with observational data (see Sec. III B), we construct a set of 1D traces from the simulation by mimicking a spacecraft fly-through with several passings over the periodic box [137]. A set of 100 virtual spacecraft trajectories are analyzed to improve statistics. In particular, we choose the directionn ¼ ð0.949; 0.292; 0.122Þ and extract the fluctuations along this given direction using cubic spline interpolation. Thus, the direction of extraction is quasiperpendicular to B 0 . The starting points of all linear trajectories are distributed uniformly in the z ¼ 0 plane. All trajectories end once they reach z ¼ L z . To avoid spurious edge effects, we skip the wavelet coefficients that are within a distance (along the 1D trace) of 19d i from the edges of the trace when calculating the generalized field ratios (see Sec. III C).

B. Spacecraft data
The solar wind data analysis is based on a 7-h interval from the Cluster spacecraft [138] and on a 159-s interval from the Magnetospheric Multiscale mission [139]. These intervals were previously studied by Chen et al. [140] and Gershman et al. [141], respectively. At the time of the measurement, the Cluster spacecraft were in the free solar wind far from Earth's foreshock, whereas the MMS spacacraft were in Earth's magnetosheath but well separated from the bow shock and the magnetopause. The mean plasma betas are β i ≈ 0.26 and β e ≈ 0.62 for Cluster and β i ≈ 0.27 and β e ≈ 0.03 for the MMS interval. MMS spacecraft can also measure the mean temperature anisotropy, which was T ⊥i =T ki ≈ 1.5 for ions and T ⊥e =T ke ≈ 0.8 for electrons during the interval. The analyzed data include magnetic measurements from Cluster [140] and the electron density, electron fluid velocity, and magnetic measurements from MMS [141]. We convert the spacecraft frame frequencies f sc to field-perpendicular wave numbers using Taylor's hypothesis [142]: k ⊥ ≈ 2πf sc =v 0 , where v 0 is a characteristic velocity. We take v 0 to be the magnitude of the mean solar wind speed V SW for the Cluster interval. Given the relatively small angle between B 0 and V SW during the MMS interval, we take v 0 ≈ jV SW j cosðθÞ for the MMS data, where θ ≈ 76°is the mean angle between V SW and the wave vector k, determined in Ref. [141] using the k-filtering technique. This angle was found to be relatively constant throughout the entire kinetic range. The inferred mean angle between k and B 0 , on the other hand, was very close to 90 deg [141]. To avoid edge effects, we skip a certain amount of points at the edges of the analyzed intervals when computing the wavelet-based diagnostics (see Sec. III C). In particular, we skip about 5 s on each side of the Cluster interval and about 7 s on each side of the MMS interval.
The MMS interval covers around 80 ion inertial lengths in the field-perpendicular direction. This is too short to allow for a statistically reliable analysis of intermittency, and the results are included here for reference. Nevertheless, we still use the MMS data in order to be able to analyze simultaneous magnetic field and density measurements, which is crucial for making direct contact with the KAW predictions, where density fluctuations play a major role [31,95]. Much longer suitable intervals from MMS are presently not available [141]. While the MMS interval is shorter than the typical large-scale turbulence correlation time (see, e.g., Ref. [32]), the Cluster interval features a relatively long continuous stream of measurements, covering several correlation times. It is thus more appropriate for studying intermittency, albeit with the limitation that only the magnetic measurements are available in this case.

C. Local scale extraction and generalized field ratios
The turbulent fields are decomposed locally among scales using the complex-valued, continuous Morlet wavelet transform [111,143]. The 1D Morlet wavelet basis functions ψ 1D s can be represented explicitly in spectral space as bandpass filters. Up to a normalization constant, the spectral representation is given bŷ where ΘðkÞ ¼ 1 for k > 0 and 0 otherwise, k ψ ¼ 6 is a dimensionless parameter [143], and k s is a characteristic wave number scale. The variable k s is related to the wavelet scale s via k s ¼ k ψ =s. For some field fðxÞ, the Morlet wavelet coefficientsf s ðxÞ at scale s are obtained from the inverse Fourier transform offψ s , wheref is the Fourier transform of f. The set of wave number scales fk s g is logarithmically spaced. We define a local, scale-dependent fluctuation as δf s ðxÞ ¼ c 1 Reff s g, where RefÁ Á Ág is the real part. Similarly, we define a local power spectral density as P f ðk s ; xÞ ¼ c 2 jf s j 2 =k s . The normalization constants c 1 and c 2 may be determined based on the exact parameters of the wavelet transform [111,143]. We set c 1 ¼ c 2 ¼ 1 since our results do not depend on such constant prefactors. In the following, we drop the scale subscript s, but it is to be understood that all quantities are scale dependent. When comparing the simulation results with spacecraft measurements, the 1D Morlet wavelet transform is applied to a set of virtual spacecraft trajectories extracted from the simulation (see Sec. III A). In Fig. 4, we also employ a 2D generalization of the Morlet transform [143,144]. The corresponding basis functions are represented in spectral space as where ϕ determines the direction of extraction, k ψ ¼ 6, and k s is defined the same as above. We also imposê ψ 2D s;ϕ ð0; 0Þ ¼ 0. Twelve directions for ϕ in the range from 0 to π are used and the final results are angle-averaged. The angular averaging of the local spectrum is performed after taking the squared modulus of the wavelet coefficient [144].
We use the wavelet decomposition to define a new set of statistical measures. In particular, we consider the spectral field ratios, frequently used to study wave properties (see, e.g., Refs. [29,31,36,77,104]), and generalize their definitions to investigate the impact of large-amplitude, turbulent structures on these ratios. Two sets of generalized ratios are defined. The first set is based on the scale-dependent moments [112] of the fluctuations: hjδn e j m i hjδb ⊥ j m i 2=m ; hjδn e j m i hjδb k j m i 2=m ; hjδb where m is the order of the moment and hÁ Á Ái represents a space average. In the simulation, we additionally average over a set of 100 virtual spacecraft trajectories (see Sec. III A). The parallel and perpendicular components are defined here relative to the local mean field B loc asb k ¼b ·b loc andb ⊥ ¼b −b kbloc , whereb loc ¼ B loc =jB loc j [101,102]. The local mean field is obtained from a Gaussian low-pass filter with a standard deviation σ s ¼ k s =k ψ in spectral space. Once the parallel and perpendicular components are determined, we normalize the fluctuations according to Eqs. (1)-(3). For m ¼ 2, the moment ratios yield the standard spectral ratios, defined in terms of the 1D k ⊥ spectra. On the other hand, for m > 2, the averaging becomes progressively more sensitive to the fluctuations at the tails of the probability distribution function, thus giving insight into the dependence of the ratios on high-amplitude events. We consider moments up to m ¼ 6.
Caution is needed when computing high-order statistics from finite datasets, since the tails of the probability distribution function may not be sufficiently sampled [145,146]. For example, an estimate of the maximal moment m max that can be determined accurately, based on the method presented in Ref. [145], yields typical values of m max between 3 and 4 for the Cluster interval with N ≈ 6.3 × 10 5 samples [147]. To obtain more reliable estimates, we employ the scheme of Kiyani et al. [146] (see also Ref. [110]) and remove a small fraction of the largest fluctuations at each scale until the moments appear reasonably converged. For Cluster, we find that removing 0.005% (i.e., about 30 samples) of largest fluctuations seems to be adequate. For consistency, we clip the same small fraction in the simulation when calculating the moments. Within the statistical uncertainties, the clipping does not significantly affect the results and only makes it easier to recognize true statistical trends. Because of the short duration of the MMS interval, the 0.005% clipping threshold has no effect and relatively large fractions would have to be removed to make the moments well behaved. Thus, no attempt is made to recover more reliable estimates via clipping for the MMS interval.
The ratios introduced in Eq. (12) are global measures in a sense that the average is taken over the entire ensemble. A more local measure can be obtained via conditional averaging of the local power spectral densities: hjñ e j 2 jLIM > ξi hjb ⊥ j 2 jLIM > ξi ; hjñ e j 2 jLIM > ξi hjb k j 2 jLIM > ξi ; hjb k j 2 jLIM > ξi where LIM is the local intermittency measure [143,148] and ξ is the threshold for the conditional average. We scan a range of different thresholds and study how the results depend on this choice. The LIM is defined as the local wavelet spectrum normalized to its mean at a given scale. Thus, LIM > 1 at the locations where the power spectral density exceeds its mean value. The LIM may be based on different quantities. We use everywhere the same type of LIM so that all conditional averages are constrained to the same spatial locations. In particular, we choose the LIM based on the KAW energy density (see Sec. II): LIM ≡ ðjb ⊥ j 2 þ jñ e j 2 Þ=hjb ⊥ j 2 þ jñ e j 2 i. This appears to be a reasonable choice given that the kinetic-scale structures carry both significant magnetic field and density fluctuations, as shown in what follows. For Cluster measurements, jb k j 2 is used as a proxy for jñ e j 2 to obtain the LIM under the assumption of pressure balance [Eq. (6)]. The conditional averages may eventually become energetically insignificant, since the averaging volume shrinks with the threshold ξ. To focus on the conditional ratios which are still of some energetic relevance, we require for any averaging subdomain to contain at least 1% of the total KAW energy at that scale. Estimates not satisfying the condition are omitted from the results. The typical volume fractions corresponding to the 1% energy fraction are naturally much smaller than 1%. The generalized ratios [Eqs. (12) and (13)] are evaluated from turbulence data and compared against the linear predictions (9). In the normalized, beta-dependent units [see Eqs. (1)-(3)], linear KAW theory predicts a numerical value of unity for all ratios considered above in the asymptotic limit [Eq. (7)]. Since the asymptotic range is vanishingly small in the simulation, we also compare our results against the linear wave predictions obtained from a fully kinetic numerical dispersion solver [149]. Note that the KAW predictions are formally obtained for a single plane wave. An ensemble of KAW packets with different propagation directions may exhibit slight statistical deviations from the expected values for those ratios in Eqs. (12) and (13) that depend on b ⊥ (i.e., for the ratios that are not based on the pair of fields n e and b k , which are directly related via pressure balance; see Sec. II).

IV. RESULTS
We now turn to the main results of this work. First, we characterize the statistical nature of fluctuations separately for each field in terms of the scale-dependent flatness [112]: where δf ∈ fδn e ; δb k ; δb ⊥ g represents a wavelet decomposed field (sec Sec. III C). For a field with a Gaussian probability distribution, we have F ¼ 3. Thus, high values of the flatness above 3 characterize the degree of non-Gaussianity, while a scale dependence of F points towards intermittency in the classical sense of the term as a departure from selfsimilarity [150]. The scale-dependent flatness results obtained from solar wind measurements and from the 3D fully kinetic simulation are compared in Fig. 3. To illustrate the statistical uncertainties, we add error bars to the flatness measurements. To obtain the error bars, we calculate the moments separately on a number of nonoverlapping subsets and then use these as input for a jackknife error estimate of the flatness [151].
As evident from the results presented in Fig. 3, the kinetic-scale fluctuations exhibit signatures of non-Gaussian statistics with flatness values above 3. The departure from Gaussian statistics is altogether largest for the Cluster interval. In agreement with previous works [102,152], the Cluster statistics are nearly scale independent at kinetic scales. This is in contrast with the simulation where the flatness is scale dependent and only gradually increases above the Gaussian value with decreasing scale, presumably due to the finite simulation domain size. All three analyzed fields in the simulation exhibit similar statistical properties, in agreement with the Cluster measurements of δb ⊥ and δb k below the ρ i scale. It is worth pointing out that a previous analysis of density intermittency in the free solar wind [110] found δn e flatness values comparable to our Cluster results, which do not include the electron density. As already mentioned (Sec. III B), the MMS interval is too short to allow for a statistically reliable analysis of intermittency. Nevertheless, within the relatively large uncertainties, the turbulence statistics appear to be mildly intermittent for the MMS interval.
Next, we inspect the spatial structure of fluctuations in the simulation. The fluctuations in a given x-y plane are visualized in Fig. 4. In Figs. 4(a)-4(c) we plot the fields in the range k ⊥ d i ¼ ½5; 10 by summing up the wavelet decomposed fluctuations in that range [23,87,100,143] using 8 logarithmically spaced scales. Note that k ⊥ d i ¼ 10 already corresponds to the electron inertial scale d e ¼ 0.1d i due to the reduced ion-electron mass ratio employed here. A distinct feature seen in Fig. 4 is the excellent matching of the δn e and δb k fluctuation profiles. Although not obvious considering the full range of kinetic effects retained in the simulation, the latter directly implies that the structures are pressure balanced according to Eq. (6). In Figs. 4(d)-4(i) we consider the local 2D wavelet spectra at k ⊥ d i ¼ 5 [Figs. 4(d)-4(f)] and at k ⊥ d i ¼ 10 [Figs. 4(g)-4(i)]. The spatial distribution of the spectral energy density is nonuniform and the peaks in the spectra at different scales tend to be concentrated around the same spatial locations. This local coupling across different scales is a characteristic feature of coherent structures [34,143,148]. It is also seen that the nonuniformity increases at smaller scales, consistent with the growth of the flatness with k ⊥ (Fig. 3). Finally, while the n e and b k wavelet spectra match very well, the local b ⊥ spectra match the former only in a rather loose sense. This is as well consistent with theoretical expectations, since no general linear relation exists between b ⊥ and n e in real space according to the KAW turbulence theory (see Sec. II). The fact that an apparent weak coupling is seen at all points towards the importance of nonlinear effects in shaping the local fluctuations.
Since aggressively scale-filtered results can be a somewhat misleading, we show in Fig. 5 the perpendicular magnetic and density fluctuations over the entire subion range and beyond, together with their corresponding (nonfiltered) perpendicular gradients. The perpendicular cross sections of the various subion-scale structures broadly resemble either sheets or circular shapes. A 3D inspection of the fluctuations (not shown) confirms that the structures are indeed elongated in the field-parallel direction, consistent with the anisotropy estimate presented in Fig. 2. It is also interesting to note that the small-scale perpendicular gradients of b ⊥ and n e tend to form structures that are aligned with respect to each other (see also Sec. V).
To further investigate the perpendicular pressure balance [Eq. (6)], we compute the wavelet cross-coherence [34,108,153] between n e and b k in the simulation and for the MMS interval, using the 1D Morlet wavelet transform [111]. High values of cross-coherence close to unity indicate a strong local phase synchronization between two signals. The results are compared in Fig. 6. Arrows are used to show the phase between the two fields. A strong phase synchronization is seen at subion scales of the simulation and of the MMS interval. With most arrows pointing to the left in Fig. 6, the results strongly suggest that the density and parallel magnetic field fluctuations tend to be anticorrelated and thus fulfill the pressure balance [Eq. (6)] to a good approximation [154]. This conclusion is in agreement with previous works based on MMS data [108,141]. Qualitatively similar results were also obtained by studies of pressure balance in the MHD-scale range (see, e.g., Refs. [155,156]). As noted in Sec. II, pressure balance is a rather generic feature of low-frequency dynamics and as such it rules out the high-frequency whistler waves. We mention that we also checked the wavelet coherence between n e and (different components of) b ⊥ . As expected (see Sec. II), n e and b ⊥ generally do not exhibit a strong cross-coherence (not shown). However, at rare times of high coherence we often observe a relative phase close to 90 deg, consistent with the elliptical polarization of the KAW [57,100].
Finally, we turn to the central subject of the present paper and present the generalized field ratios results. The generalized ratios obtained from the 3D fully kinetic simulation and from spacecraft measurements are plotted in Fig. 7. In Figs. 7(a)-7(f) we show the moment ratios [Eq. (12)], and in Figs. 7(g)-7(l) we display the conditional ratios [Eq. (13)], conditioned on the local KAW energy density (see Secs. II and III C). Dashed horizontal lines denote the linear asymptotic KAW predictions [Eq. (9)], while red lines show the more accurate linear predictions obtained from a fully kinetic plasma dispersion relation solver [149]. A wave propagation angle of 89.9 deg (with respect to B 0 ) was used for solving the numerical dispersion relations [157]. Dotted vertical lines indicate the characteristic ion and electron scales. The choice m ¼ 2 or ξ ¼ 0 for the moment ratios and conditional ratios, respectively, yields the standard spectral ratios. We mention that the maximal meaningful threshold ξ (see Sec. III C) grows with k ⊥ for the simulation and for MMS data since the fluctuations become more intermittent at smaller scales. The Cluster interval exhibits highly non-Gaussian statistics for all k ⊥ (Fig. 3) and thus allows for the use of high thresholds over the whole range. The higher the moment m or threshold ξ, the more sensitive the ratios are to the largeamplitude events. Note that the term "large amplitude" refers here to the locally enhanced fluctuations and power spectral densities, compared to their scale-dependent mean square values. Spatially, these intense bursts of local activity are not distributed incoherently but instead form distinct patterns (Figs. 4 and 5), typically associated with turbulent coherent structures. The locally enhanced, non-Gaussian fluctuations go hand in hand with the nonuniform power spectral densities, since it can be shown that the spatial variability of the energy spectrum and the scaledependent flatness are directly related [159]. Higher moments assign smaller statistical weights to the loweramplitude background fluctuations in favor of the large-amplitude events. Similarly, the conditional averages discard the locations with a spectral energy density below the threshold ξ and therefore measure the field ratios within the most energetic structures only.
As seen in Fig. 7, the generalized ratios either follow the linear wave predictions (red lines) throughout the transition region (k ⊥ ∼ 1=ρ i ), where shear Alfvén waves convert into KAWs, or converge toward the KAW predictions below the ρ i scale. This is consistent with KAW theory, where ρ i is the relevant transition scale [57,95,160]. A notable deviation from linear wave predictions is seen at subelectron scales of the simulation. A similar trend was previously seen in gyrokinetic simulations [101,161], which possibly implies that the wavelike activity is terminated below the ρ e scale via collisionless damping [60]. The moment ratios depend only weakly on the order m, which shows that the wave predictions are reasonably satisfied not only in terms of low-order statistics but also in the sense of higher-order moments related to intermittency. The conditional ratios between jñ e j 2 and jb k j 2 [Figs. 7(g)-7(h)] appear to be somewhat less scattered compared to the other two conditional ratios [Figs. 7(i)-7(l)], which exhibit a slight tendency to deviate further from linear predictions with increasing LIM threshold ξ. On the other hand, the subion-scale deviations attributed to large-amplitude events are altogether only moderate, such that the generalized ratios remain in reasonable agreement with the KAW polarization properties. Thus, the wavelike features are not exclusively limited to low-amplitude fluctuations but also carry over to the high-amplitude structures. Qualitatively similar results are obtained in the simulation as well as from spacecraft measurements, which suggests that the observed properties are a common feature of the kinetic turbulence. While some of the measured ratios might allow for alternative interpretations (see Sec. II), the fact that an agreement is seen for all generalized ratios and for different values of the ambient ion and electron plasma beta [recall that the normalizations (1)-(3) are beta dependent] implies that the fluctuations are indeed of the KAW type. Moreover, the hypothetical possibility that the field polarizations are not related to wavelike activity is contradicted by our estimate of the nonlinearity parameter (Fig. 2), and we are unaware of any other concrete theoretical prediction that could consistently explain all of our measurements. The observed order unity preservation of linear wavelike properties in kinetic-scale turbulent structures constitutes the main result of this work, together with the supplemental evidence presented in Figs. 2-6.

V. DISCUSSION
The general idea that intermittent turbulent structures may go hand in hand with wave physics was implied in a number of recent works. In a numerical study of MHD turbulence [50], a scale-independent probability distribution with an order unity mean was found for the local (in space) ratio between the linear wave timescale and nonlinear timescale, whereas all other quantities of interest exhibited scale-dependent, intermittent statistics. The former property was later successfully exploited to design a statistical model of intermittency in MHD turbulence [162]. Moreover, Howes [163] recently developed a dynamical model for current sheet formation via nonlinear Alfvén wave collisions. Subsequent numerical studies showed that current sheets emerge also in a more realistic configuration with two counterpropagating, initially separated wave packets [164], and that the modes produced during nonlinear wave interaction can be themselves characterized as Alfvén waves [165]. The above-described developments in astrophysical plasma physics are not a unicum in turbulence research. For instance, similar ideas have been put forward in the context of rotating, inertial wave turbulence [11].
In KAW turbulence, wave packet interactions are relatively complex due to the dispersive nature of the KAW, which allows for the nonlinear coupling between copropagating waves [57,166]. This prevents a direct application of the model by Howes [163], which is based on incompressible MHD equations, to subion-scale dynamics. On the other hand, our results shown in Fig. 7 imply that the subion-scale turbulent structures may as well be viewed as nonlinearly interacting wave packets. This quantitative evidence encourages the design of new models of turbulent heating in low-collisionality plasmas that incorporate intermittency by exploiting the wavelike character of coherent structures [167]. Additional evidence is also provided by a recent observational study [94] and by gyrokinetic turbulence simulations [30,35], which find the energy transfer between fields and particles to be dominated by particles satisfying the resonance condition for Landau damping of (kinetic) Alfvén waves, even though the transfer tends to be spatially concentrated in the vicinity of intense current sheets [28,32,33,93].
Since the nonlinear time of a turbulent eddy is inversely proportional to its fluctuation amplitude [14], one might be tempted to conclude that the large-amplitude nonlinear structures evolve on a timescale faster than the wave dynamics. A number of aspects may be pointed out here. Unlike the nonlinear time, the linear wave time is proportional to 1=k k [56,57]. For a realistic turbulent setup, the latter may be interpreted as a parallel correlation length l k ∼ 1=k k of an eddy [57,107]. In magnetized plasma turbulence, information is transmitted along the field lines by waves. It then follows from a causality principle that the maximal parallel correlation is set by the distance a wave packet can cover within the lifetime of an eddy [12,57,168]. This implies that the wave timescale is approximately limited from above by the nonlinear time, which is indeed consistent with our measurement of the anisotropy in Fig. 2. An analogous causality argument can be given for rotating turbulence in the strongly anisotropic (k k ≪ k ⊥ ) regime [12]. It was also recently argued, based on a set of 3D gyrokinetic simulations [165], that the nonpropagating k k ¼ 0 modes do not play a notable role in the mediation of energy transfer in MHD, provided that appropriate measures are taken to avoid periodic boundary artifacts. Secondly, it is worth considering the empirical fact that turbulence does not self-organize into structures of arbitrary geometry. Instead, a number of studies of kinetic range turbulence find that the structures tend to acquire relatively simple shapes, resembling either elongated sheets or tubes [23,30,33,83,88,91,169]. For k k ≠ 0, the idealized geometric versions of these two (i.e., sheets with a 1D perpendicular profile and circularly symmetric, fieldaligned tubes) correspond to exact wave solutions of the nonlinear Eqs. (4) and (5), as explained in Sec. II. The wandering of magnetic field lines in turbulent flows precludes such highly symmetric configurations [55,170,171]. However, even if a structure only resembles an idealized sheet or circular tube, the nonlinearity is locally weakened. Some of these circumstances were also appreciated by Terry and Smith [169,172], who argued that the magnetic shear at the edges of large-amplitude circular filaments and sheets prevents the structures from mixing with the background incoherent KAW fluctuations. Moreover, local depletions of the nonlinearity have been associated with the emergence of coherent structures in a variety of turbulent systems [173][174][175][176][177].
Is it possible to provide evidence for the above suggestion? To this end, we make use of the fact that in the first approximation, ∇ ⊥ n e ∝ẑ × u ⊥e and ∇ ⊥ ψ ∝ẑ × b ⊥ , at subion scales [178], where u ⊥e is the perpendicular electron fluid velocity. If u ⊥e × b ⊥ ¼ 0, the nonlinear coupling in Eq. (4) vanishes. We then consider the scale-dependent sine of the alignment angle, sin θ ≡ jδu ⊥e × δb ⊥ j=jδu ⊥e jjδb ⊥ j, conditionally averaged on a given LIM threshold for the MMS and simulation data. The same as in the rest of this paper, the scale decomposition is implemented using the Morlet wavelet transform and the perpendicular components are defined with respect to the local mean field. Our results are shown in Fig. 8. According to the above discussion, the alignment angle should be preferentially reduced within the energetic structures (corresponding to high LIM values). The results in Fig. 8 are indeed consistent with our expectation, although the preferential alignment within energetic structures is less pronounced for the MMS interval. The latter may be related to the fact that the fluctuations in the MMS interval are only mildly intermittent (see Fig. 3). We also add that an intermittent field alignment, analogous to the one studied here over the kinetic range of scales, was previously reported in the context of MHD turbulence (see, e.g., Refs. [179,180]). Here, we suggest that such intermittent, kinetic-scale alignment helps to preserve the linear KAW properties in large-amplitude structures via local nonlinearity depletions.
Finally, there exist numerous other examples of the importance of waves in complex nonlinear plasma phenomena. Classic examples include the nonlinear steepening of MHD waves into shocks [181][182][183][184][185] and solitary waves in magnetized plasmas [130,[186][187][188][189]. More recently, there has been a growing interest to relate soliton solutions with kinetic-scale structures in turbulent space plasmas [19,23,34,83]. Indeed, it cannot be ruled out that at least some turbulent structures are better described as solitons than dynamically evolving wave packets. Similar to the exact linear wave solutions discussed in Sec. II, soliton solutions typically require a rather idealized, highly symmetric and/or low-dimensional geometry. Thus, for the same reason as the turbulent structures cannot be exactly linear waves (namely, stochastic field line wandering), they are probably not precisely solitons either, unless perhaps if the system is far from a state of homogeneous, fully developed strong turbulence.

VI. CONCLUSIONS
In the present paper, we employ high-resolution observational and kinetic simulation data to study the relationship between wavelike physics and structure formation in astrophysical kinetic plasma turbulence. Observational data are based on in situ solar wind measurements from the Cluster and MMS spacecraft, and the simulation results are obtained from an externally driven, 3D fully kinetic simulation. We introduce appropriate generalizations of the frequently employed spectral field ratios to probe the statistical field polarizations within the large-amplitude energetic structures. We find that both the lower-amplitude background fluctuations as well as the intense structures satisfy linear KAW predictions to order unity. This implies that the turbulent structures themselves approximately preserve a linear wave footprint. In other words, our novel analysis highlights the possibility that wavelike features and turbulent structure formation are essentially inseparable from each other. As such, our results challenge one of the presently common views of the "coexistence" of waves and structures. Furthermore, it is suggested that a linear wave character in kinetic-scale structures is preserved with the aid of local nonlinearity depletions.
The only known framework capable of providing a reasonable theoretical basis for the interpretation of our results appears to be the KAW turbulence phenomenology [56,57,95]. Known alternatives, such as whistler turbulence [75,[77][78][79], cannot explain our results. In contrast to our results, the high-frequency whistlers are not pressure balanced and they carry only minor density fluctuations [31,108]. In a β ∼ 1 plasma, whistler waves are also expected to be rather strongly damped for k k d i ≲ 1 [95], which contradicts our anisotropy estimate in Fig. 2. In the KAW turbulence context, an interpretation similar to recent FIG. 8. Sine of the scale-dependent alignment angle between the perpendicular electron fluid velocity and magnetic field, conditionally averaged on a given LIM threshold ξ (see text for details). Dashed horizontal lines indicate the expected value for a randomly distributed angle. developments in MHD turbulence emerges (see, e.g., Refs. [162][163][164][165]). Namely, the kinetic-scale structures could perhaps be described as nonlinearly evolving, localized KAW packets. The possibility that at least a fraction of the subion-scale structures are sheetlike, as implied by our simulation results (see also Refs. [88,91]), lends credence to recent works emphasizing the role of reconnection in subion-scale turbulence [82,[190][191][192][193]. In our present understanding, the results of this work do not preclude the reconnection scenario, although the mere presence of sheetlike structures alone is not a sufficient condition for a reconnection-mediated type of turbulence.
Finally, we mention that the general approach employed here is not exclusively limited to kinetic range turbulence in astrophysical plasmas and we hope it might find exciting applications in a broad range of turbulent systems where waves and structures have been observed [2,7,11,17,18,20,22]. An immediate extension of the method lends itself in the context of MHD range plasma turbulence, where it could be used to study the interplay between structures and waves based on a generalized Alfvén ratio [9,27]. Another potential application includes rotating turbulence [7,11,12,22], where the generalized spectral field ratios could be employed to test the presence of inertial-wave polarization locally within the energetic columnar structures.