# Heat Flux Signatures of Finite Larmor Radius Ion Mixing Driven by Magnetopause Kelvin-Helmholtz Waves

Vertti Tarvus, Lucile Turc, Adriana Settino, Kevin Blasl, Giulia Cozzani, Yann Pfau-kempf, Jonas Suni, Minna Palmroth

► To cite this version:

Vertti Tarvus, Lucile Turc, Adriana Settino, Kevin Blasl, Giulia Cozzani, et al.. Heat Flux Signatures of Finite Larmor Radius Ion Mixing Driven by Magnetopause Kelvin-Helmholtz Waves. Journal of Geophysical Research Space Physics, 2025, 130 (11), 10.1029/2025JA034175 . insu-05352692

HAL Id: insu-05352692

https://insu.hal.science/insu-05352692v1

Submitted on 7 Nov 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
JGR Space Physics

RESEARCH ARTICLE

10.1029/2025JA034175

**Key Points:**

* Magnetopause Kelvin-Helmholtz waves under northward IMF conditions show enhanced ion heat flux transverse to the magnetic field

* The heat flux is consistent with being generated by ion finite Larmor radius mixing, which may contribute to magnetopause plasma transport

* Ion heat flux measurements could be useful for studying Kelvin-Helmholtz wave structure

**Correspondence to:**

V. Tarvus,
vertti.tarvus@helsinki.fi

**Citation:**

Tarvus, V., Turc, L., Settino, A., Blasl, K., Cozzani, G., Pfau-Kempf, Y., et al. (2025). Heat flux signatures of finite Larmor radius ion mixing driven by Magnetopause Kelvin-Helmholtz waves. *Journal of Geophysical Research: Space Physics*, 130, e2025JA034175. https://doi.org/10.1029/2025JA034175

Received 8 MAY 2025

Accepted 17 OCT 2025

<figure><img src="image_2.png" /></figure>

# Heat Flux Signatures of Finite Larmor Radius Ion Mixing Driven by Magnetopause Kelvin-Helmholtz Waves

Vertti Tarvus<fnref n="1" /> ID, Lucile Turc<fnref n="1" /> ID, Adriana Settino² ID, Kevin Blasl²,³ ID, Giulia Cozzani<fnref n="1" />,⁴
Yann Pfau-Kempf<fnref n="1" />,⁵ ID, Jonas Suni<fnref n="1" /> ID, and Minna Palmroth<fnref n="1" />,⁶ ID

**Abstract** The Earth's magnetopause flanks are commonly unstable to the Kelvin-Helmholtz (K-H) instability, which is excited by the velocity shear between the magnetosphere and the magnetosheath. The K-H instability generates magnetopause surface waves, which steepen into non-linear vortices that mix plasmas from the two regions, and transport magnetosheath plasma into the magnetosphere. One of the candidate mechanisms for the plasma mixing identified in numerical simulations are finite Larmor radius (FLR) effects, which become active along the magnetopause when it is deformed by K-H waves, and steepened to a thickness comparable to the thermal particle Larmor radius. Diffusive magnetosheath plasma transport into the magnetosphere occurs when the FLR mixing region is disrupted by turbulence. In this work, we use Magnetospheric Multiscale Mission data to show in situ evidence of an ion FLR ion mixing region during magnetopause K-H wave boundary crossings under northward interplanetary magnetic field conditions. We show that these crossings contain ion heat flux enhancements that are consistent with being generated by ion FLR mixing across a temperature gradient. In accordance with kinetic transport theory and our past numerical results, these enhancements are transverse to the local magnetic field and temperature gradient. Our present findings indicate that the K-H waves produce copious FLR ion mixing, which could have an important role in the K-H driven transport of magnetosheath plasma into the magnetosphere. We also demonstrate that heat flux measurements could be used for identifying K-H waves at different evolutionary stages.

## 1. Introduction

One of the fundamental processes mediating the interaction between the solar wind and the Earth's magnetosphere is the Kelvin-Helmholtz (K-H) instability (Hwang et al., 2023; Masson & Nykyri, 2018), which is excited by a velocity shear. Surface waves generated by the instability develop at the magnetopause, the boundary shielding the magnetosphere from the solar wind (Hasegawa et al., 2004). The growth rate is maximized for a magnetic field transverse to the shear (Chandrasekhar, 1961), and the most favorable conditions for growth occur at the low-latitude flanks of the magnetopause during northward IMF, when dayside magnetic reconnection is inefficient. Nonetheless, magnetopause K-H waves exist under a variety of IMF orientations (Blasl et al., 2022; Hwang et al., 2011; Kavosi & Raeder, 2015).

The waves are characterized by their non-linear steepening, which causes the surface waves to evolve into vortices. Due to this property, the instability can enhance the transport of solar wind momentum into the magnetosphere via a viscous interaction (Axford, 1964; Axford & Hines, 1961; Miura, 1984). Transport of solar wind plasma into the magnetosphere can occur if the frozen-in plasma condition is broken by magnetic reconnection. Reconnection may be triggered at the twisted magnetic field lines within the vortices (Nakamura et al., 2008; Nykyri & Otto, 2001; Pu et al., 1990), or at mid-latitudes, where the magnetic field lines become twisted due to the vortex motion at the lower latitudes (Faganello et al., 2012). Finite Larmor radius (FLR) effects enable additional plasma mixing and transport via particle diffusion across magnetic fields (Cowee et al., 2009; Ma et al., 2019; Matsumoto & Hoshino, 2006; Terasawa et al., 1992).

Magnetopause K-H waves contribute to the formation of the low-latitude boundary layer (LLBL), which consists of mixed magnetosheath and magnetosphere plasmas (Nakamura, 2021). The LLBL is particularly prominent under northward IMF conditions (Mitchell et al., 1987) and cold, dense plasma may be transported deep into the plasma sheet (Terasawa et al., 1997). High-latitude reconnection (Song & Russell, 1992) and kinetic Alfvén waves (KAW) (Chaston et al., 2008) also contribute to the formation of the LLBL, and they can work in parallel

<footnote n="1">¹University of Helsinki, Helsinki, Finland, ²Space Research Institute Graz, Plasma Physics, Graz, Austria, ³University of Graz, Institute of Physics, Graz, Austria, ⁴LPC2E, CNRS/CNES/University of Orléans, Orléans, France, ⁵CSC - IT Center for Science Ltd., Espoo, Finland, ⁶Finnish Meteorological Institute, Helsinki, Finland</footnote>

<footnote>© 2025. The Author(s).
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>
with the K-H instability. A pre-existing LLBL formed by high-latitude reconnection may promote the K-H instability (Nakamura et al., 2017), and KAWs may be generated by K-H waves through mode conversion (Chaston et al., 2007).

Identifying the contribution of the K-H instability to the formation of the LLBL and the associated solar wind entry into the magnetosphere is an ongoing effort (Hwang et al., 2023). In situ multi-spacecraft measurements have provided detailed evidence of the small-scale dynamics of magnetopause K-H waves. Magnetic reconnection driven by K-H waves was first observed by the Cluster spacecraft (Hasegawa et al., 2009; Nykyri et al., 2006), and later by the Magnetospheric Multiscale Mission (MMS) (Eriksson, Lavraud et al., 2016; Eriksson, Wilder et al., 2016; Li et al., 2016; Vernisse et al., 2016). MMS has also observed flux ropes (Hwang et al., 2020), flux transfer events (Kieokaew et al., 2020), turbulence (Hasegawa et al., 2020; Sorriso-Valvo et al., 2019; Stawarz et al., 2016), and large-amplitude electrostatic wave generation (Wilder et al., 2016) in association with magnetopause K-H waves. The high resolution of MMS measurements has also allowed identifying kinetic signatures characteristic of K-H waves. Settino et al. (2021) showed that rolled-up K-H vortices exhibit edges with an enhanced current density and agyrotropy (as defined by Swisdak, 2016), and an interior with a low current density and an enhanced temperature anisotropy. Settino et al. (2022) introduced a mixing parameter measuring the ratio of magnetosphere and magnetosheath particles, and showed that the behavior of the ion and electron mixing parameters during K-H wave encounters can be used to distinguish waves at different evolutionary stages.

While the picture of magnetopause K-H waves provided by in situ measurements is comprehensive, the role of FLR effects, and in particular their contribution to the K-H driven plasma transport, remain unknown in situ. In kinetic simulations (e.g., Ma et al., 2019; Terasawa et al., 1992), FLR mixing of magnetosheath and magnetosphere plasmas occurs at the magnetopause surface as it becomes steepened by K-H wave growth to a thickness of a Larmor radius, resulting in a continuous mixing layer along the deformed magnetopause. When the mixing layer is disrupted by turbulence, mixed plasma is diffusively transported into the magnetosphere (Cowee et al., 2009; Matsumoto & Hoshino, 2004, 2006; Matsumoto & Seki, 2010). K-H driven turbulence has been reported in situ (Blasl et al., 2023, 2024; Stawarz et al., 2016), but it is not known if it operates in tandem with FLR mixing, as there is no established method for identifying K-H driven FLR mixing from spacecraft data.

In our previous work (Tarvus et al., 2024, Paper I from now on), we looked for detectable signatures of ion FLR mixing in a hybrid-Vlasov simulation of the K-H instability in a northward IMF magnetopause setting. We used the kinetic transport theory (Braginskii, 1965) as a reference, and showed that a heat flux exists at the magnetopause FLR mixing layer due to a temperature gradient $\nabla_{\perp}T$ transverse to the magnetic field direction $\hat{b}$. The heat flux results from the kinetic energy difference between the hot and cold particles whose Larmor orbits cross at the temperature gradient. As stated by Braginskii (1965), the heat flux direction is given by $\hat{b} \times \nabla_{\perp}T$ (see Section 2.1 for further details). We showed in Paper I that this expression matches well with the heat flux direction at the magnetopause FLR mixing layer, and proposed that such transverse heat flux could be used as a signature of K-H driven ion FLR mixing. We also proposed that the variation of the heat flux direction at the magnetopause could be used to study K-H wave structure, complementing existing methods (Plaschke et al., 2016; Settino et al., 2021, 2022).

In the present work, our goal is to find whether the transverse ion heat flux signature identified in Paper I exists in situ, and if ion heat flux measurements provide information about K-H driven FLR mixing at the Earth's magnetopause. Addressing these questions will help determine if such mixing offers a viable path for the magnetosheath plasma entry into the magnetosphere, in conformity with kinetic simulations. We use MMS measurements from a well-known K-H event observed on 8 September 2015 at the dusk flank magnetopause under northward IMF conditions (Eriksson, Lavraud, et al., 2016). We utilize the four-spacecraft setup of MMS, and perform timing analysis (Plaschke et al., 2016; Shi et al., 2019) at magnetopause crossings to evaluate the K-H wave boundary orientation and the direction of $\nabla_{\perp}T$. This way, we calculate $\hat{b} \times \nabla_{\perp}T$ and compare it to the measured heat flux direction.

In the next section (Section 2), we provide an overview of the heat flux described by Braginskii (1965) and its properties in the context of K-H waves (Paper I). We then detail our method for evaluating the heat flux vector and $\hat{b} \times \nabla_{\perp}T$. In Section 3, we provide an overview of the ion heat flux data, and demonstrate the presence of heat flux peaks with a strong transverse component. We then perform an analysis of individual peaks, demonstrating
agreement between the measured heat flux and $\hat{b} \times \nabla_{\perp} T$ at K-H wave boundary crossings, and present evidence of active K-H driven FLR mixing. Lastly, we compare our results with our prior numerical work (Paper I). We demonstrate the usage of heat flux measurements for studying the K-H wave evolution, and present evidence of a temporally evolving FLR mixing layer. We conclude with a discussion of our findings and conclusions in Section 4.

## 2. Methods

### 2.1. Evaluating the Transverse Heat Flux

As shown by Braginskii (1965), particle Larmor motion at a temperature gradient $\nabla_{\perp} T$ transverse to a uniform background magnetic field $\vec{B}$ produces a heat flux transverse to both $\nabla_{\perp} T$ and $\vec{B}$. In the plasma rest frame, the unidirectional particle fluxes ($\sim nv_{th}$ with number density $n$ and thermal speed $v_{th}$) through the plane spanned by $\nabla_{\perp} T$ and $\vec{B}$ are equal, and there is no net particle flux. However, the unidirectional particle energy fluxes ($\sim nk_B(T \pm \delta T) v_{th}$ with Boltzmann constant $k_B$, and a small temperature change $\delta T$) are not equal, and their difference is a heat flux of the order of

$$ \sim nv_{th}k_B(T + \delta T) - nv_{th}k_B(T - \delta T) \quad (1) $$

$$ = 2nv_{th}k_B\delta T \quad (2) $$

$$ = 2 \frac{nk_B T}{m\omega_c} k_B \nabla_{\perp} T, \quad (3) $$

where the last expression is obtained by assuming $\delta T \approx n_L \nabla_{\perp} T$ ($n_L = v_{th}/\omega_c$ is the thermal Larmor radius, with Larmor frequency $\omega_c$), and rearranging using $v_{th} = \sqrt{k_B T/m}$ ($m$ is the particle species mass). The exact form of the transverse heat flux $\vec{q}_{\Lambda,i}$ due to ion FLR motion derived by Braginskii (1965) is

$$ \vec{q}_{\Lambda,i} = \kappa_{\Lambda,i} \hat{b} \times \nabla_{\perp} (k_B T_i), \quad (4) $$

where $\hat{b}$ is the magnetic field unit vector and $\kappa_{\Lambda,i}$ is the ion thermal conductivity

$$ \kappa_{\Lambda,i} = \frac{5nk_B T_i}{2m_i \omega_{ci}}. \quad (5) $$

In Paper I, we numerically demonstrated that such ion heat flux emerges along the K-H steepened magnetopause under northward IMF conditions, since $\hat{b}$ is relatively constant across the magnetopause and the steepening causes a large $\nabla_{\perp} T$. We identified factors that cause the measured heat flux to differ from $\vec{q}_{\Lambda,i}$, but concluded that the measured magnitude is typically within an order of magnitude of $|\vec{q}_{\Lambda,i}|$, and the angle relative to $\hat{b} \times \nabla_{\perp} T$ less than 15°.

The measured heat flux differs from $\vec{q}_{\Lambda,i}$ due to for example, $\nabla_{\perp} T$ being confined to a layer of a thermal Larmor radius thickness, which reduces the heat flux magnitude by limiting the maximum temperature of the particles contributing to the heat flux. The magnitude is also modulated by another FLR effect, a temperature anisotropy due to a velocity shear (see Appendix in Paper I). The heat flux vector is typically well-aligned with the magnetopause surface due to the low curvature of K-H waves (wavelength $\gg n_L$), but small-scale irregularities cause variations in the direction. A small magnetic field rotation ($\sim 5^\circ$) also causes small variations due to a non-constant $\hat{b}$, although the added magnetic tension can limit the effect of irregularities.

Since Equation 4 approximates well the heat flux direction and order of magnitude in a K-H active northward IMF magnetopause setting, we use it as a reference for our in situ heat flux measurements. $\kappa_{\Lambda,i}$ can be readily evaluated using ion and magnetic field measurements, and the heat flux direction can be found by estimating the temperature gradient direction via multi-spacecraft methods. We will further detail this procedure in the next section.
## 2.2. In-Situ Measurements

We use measurements from the MMS spacecraft during an outbound transit of the dusk flank magnetopause, which occurred under northward IMF conditions on 8 September 2015, between 9 and 12 UTC. During this event, the four MMS spacecraft flew in a tetrahedral formation with a separation of 150–185 km. Magnetopause K-H surface waves were encountered between 10 and 11:30 UTC (approximate period ~ 63 s). The event was first reported by Eriksson, Lavraud, et al. (2016), and has been further studied by for example, Li et al. (2016); Vernisse et al. (2016); Stawarz et al. (2016); Sorriso-Valvo et al. (2019); Settino et al. (2022).

For ion measurements, we make use of the full temporal resolution (150 ms) burst mode particle data provided by the Fast Plasma Investigation instrument (Pollock et al., 2016), available between 10:07:04–11:27:34 UTC. Unless stated otherwise, we use heat flux calculated directly from the ion VDF as

$$ \vec{q}_i = \int \frac{1}{2} m_i |\vec{v} - \vec{v}_i|^2 (\vec{v} - \vec{v}_i) f_i(\vec{v}) d\vec{v} \quad (6) $$

$$ = \vec{Q}_i - \vec{P}_i \cdot \vec{v}_i - \frac{3}{2} P_i \vec{v}_i - \frac{1}{2} m_i n v_i^2 \vec{v}_i, \quad (7) $$

where $\vec{v}$ is the velocity, and $\vec{v}_i$, $\vec{P}_i$, $\vec{Q}_i$, and $f_i$ are the ion bulk velocity, pressure tensor, total energy flux vector, and VDF, respectively. Equation 7 is the form used in spacecraft measurements as it is conveniently evaluated using the base velocity moments of the VDF. In this study, we exclude one-count VDF bins to reduce the uncertainty of the calculation. To check that the results are not affected by low particle counts, we have explicitly verified that consistent results are obtained with a 5 data point time-averaged VDF.

For the magnetic field measurements, we use data provided by the FluxGate Magnetometer instrument (Russell et al., 2016), with a 128 Hz temporal resolution. Whenever we do calculations involving both the particle and field data, the field data is downsampled to the particle data cadence. We will mainly focus on single spacecraft data from the MMS1 spacecraft, but use all four MMS spacecraft to calculate boundary normal vectors via the four-spacecraft timing method (Harvey, 1998; Shi et al., 2019).

We perform our analysis in an orthonormal coordinate system $(x', y', z')$, based on the magnetosheath velocity and magnetic field vectors reported by Eriksson, Lavraud, et al. (2016). These are given in the Geocentric Solar Ecliptic (GSE) coordinate system by $\vec{v}_{GSE} = 330(-0.76, 0.65, -0.04)$ km s⁻¹ and $\vec{B}_{GSE} = 72(0.35, -0.17, 0.92)$ nT, respectively. Using these, we define $(x', y', z')$ in the GSE system as

$$ \hat{x}' = \frac{\hat{y}' \times \vec{B}_{GSE}}{|\hat{y}' \times \vec{B}_{GSE}|} = (0.65, 0.75, -0.11); \quad (8) $$

$$ \hat{y}' = \frac{\vec{v}_{GSE}}{|\vec{v}_{GSE}|} = (-0.76, 0.65, -0.04); \quad (9) $$

$$ \hat{z}' = \frac{\hat{x}' \times \hat{y}'}{|\hat{x}' \times \hat{y}'|} = (0.04, 0.11, 0.99). \quad (10) $$

The coordinate system is analogous with the one used in our simulations in Paper I, allowing a straightforward comparison with our prior numerical results (Section 3.3). We also use an azimuthal angle $\phi$ to characterize directions in the $x' - y'$ plane, that is approximately transverse to the northward-oriented magnetic field. For a vector $\vec{u}$, $\phi$ is given by

$$ \phi(\vec{u}) = \tan2(u_y, u_x) \in [-180^\circ, +180^\circ], \quad (11) $$

with $\phi = 0^\circ$ fixed in the direction of the positive $x'$-axis.

Our coordinate system is illustrated in Figure 1, which depicts an idealized dusk flank magnetopause deformed by a K-H wave under northward IMF conditions. We also illustrate our method for finding the heat flux direction
<figure><img src="image_3.png" /><figcaption>**Figure 1.** An idealized cartoon of the dusk-flank magnetopause deformed by a K-H wave under northward IMF conditions. The magnetic field $\hat{b}$ is assumed to be in the $+z'$-direction. As the wave propagates in the $y'$-direction, a spacecraft along the trajectory shown by the dotted line crosses through leading/trailing wave boundaries. At each boundary crossing, the spacecraft measures a transverse temperature gradient $\nabla_{\perp} T$ along the boundary normal vector $\hat{n}$ (blue arrow), and a transverse heat flux (Equation 4) along $\hat{b} \times \hat{n}$ (red arrow).</figcaption></figure>

given by Equation 4 ($\hat{b} \times \nabla_{\perp} T$). We use an approximation $\hat{b} \times \hat{n}$ (red arrows), where $\hat{n}$ (blue arrows) is a boundary normal vector obtained via timing analysis at K-H wave boundary crossings. We use ion temperature data for the timing analysis, and assess the reliability of our result by repeating the analysis with other parameters including the ion density and magnetic field ($x'$, $y'$)-components.

Assuming the $\hat{n}$ obtained via the timing analysis is oriented in the wave propagation direction (+$y'$), we flip the sign of $\hat{n}$ when necessary to ensure $\hat{n} \parallel \nabla_{\perp} T$. We do this based on the sign of the time derivative $d\langle T \rangle/dt$ of the temporally averaged ion temperature $\langle T \rangle$, from which high-frequency fluctuations have been filtered out. The averaging window (7.5 s, 50 data points) is chosen to be large enough to cover a typical wave boundary crossing. A trailing (leading) wave boundary corresponds to a negative (positive) $d\langle T \rangle/dt$.

For the timing analysis, we adopt the methodology by Plaschke et al. (2016), who studied the steepening of magnetopause surface waves using MMS data. Accordingly, we obtain the time separation $\tau$ between a spacecraft pair ij by maximizing the cross-correlation $P_{ij}(\tau)$ between the equal length timeseries $u_i(t + \tau)$ and $u_j(t)$ obtained by each spacecraft, with the former shifted by $\tau$. $P_{ij}(\tau)$ is given by

$$P_{ij}(\tau) = \frac{\sum_t (u_i(t + \tau) - \mu_i)(u_j(t) - \mu_j)}{\sqrt{\sum_t (u_i(t + \tau) - \mu_i)^2 \sum_t (u_j(t) - \mu_j)^2}} \in [-1, +1], \quad (12)$$

where $\sum_t$ is the sum over the timeseries, and $\mu_i, \mu_j$ are the mean values of each timeseries. We use MMS1 as the reference spacecraft i, and calculate the boundary normal vector based on $\tau$ between the spacecraft pairs ij = 12, 13, 14. As done by Plaschke et al. (2016), we set a threshold $P_{ij}(\tau) > 0.9$ for a satisfactory result. We consider a $\sum_t$ with a length of 15 s, and a $\tau \in [-5, 5]$ s.

# 3. Results

## 3.1. Overview of the Ion Heat Flux Data

In Figure 2, we present fast mode ion data and survey mode magnetic field data from 8 September 2015, when the MMS spacecraft crossed through the dusk flank equatorial magnetopause into the magnetosheath. As reported by Eriksson, Lavraud, et al. (2016), MMS was located in the magnetosphere proper, characterized by steady measurements of high-energy ions (Figure 2a), until 09:21:24 UTC (first blue line). It then moved into the low-latitude boundary layer, characterized by mixed high-energy magnetosphere ions and low-energy magnetosheath ions. During the latter part of the boundary layer crossing (~10–11:30 UTC), MMS encountered magnetopause K-H surface waves, alternating between the boundary layer and the magnetosheath while observing periodic fluctuations in the ion density, ion velocity, ion temperature, and magnetic field (Figures 2b–2e). After 11:27:40 UTC (second blue line), MMS crossed into the magnetosheath, measuring a steady flow of low-energy ions.

The ion heat flux magnitude $|\vec{q}_i|$ (Figure 2f, unprocessed fast mode data shown) is notably enhanced during the crossing of the boundary layer, exhibiting recurring peaks of the order of $10^{-1}$ mW m$^{-2}$, well above the level in both the ambient magnetosphere and magnetosheath. Figure 2g shows a comparison of the parallel and transverse components of $\vec{q}_i$, revealing the presence of a significant transverse component. Figure 2h shows the transverse ion thermal conductivity $\kappa_{\Lambda,i}$ (Equation 5). We next consider a subinterval of ion burst mode data, in order to determine if the conductivity is sufficient for generating the measured heat flux magnitude.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** MMS1 fast mode ion data and survey mode magnetic field data throughout the magnetopause crossing of 8 September 2015. Panel (a): Ion differential energy flux, (b): ion number density, (c) ion velocity, (d): ion temperature, (e): magnetic field, (f): ion heat flux magnitude, (g): ion heat flux parallel and transverse to the magnetic field, (h): the transverse ion thermal conductivity $\kappa_{\Lambda,i}$ (Equation 5). For the calculation of $\kappa_{\Lambda,i}$, we have resampled the magnetic field data at the ion data cadence. The units are chosen such that the heat flux is in the units of $W m^{-2}$ when $\kappa_{\Lambda,i}$ is multiplied by a temperature gradient in the units of $eV m^{-1}$ ($\kappa_{\Lambda,i}$ has been multiplied by the square of the conversion factor $1.602176634 \times 10^{-19} J eV^{-1}$). The yellow highlighted interval shows where burst mode ion data is available (10:07:04–11:27:34 UTC), the gray highlighted interval is depicted in Figure 3.</figcaption></figure>

Following the format of Figure 2, Figure 3 shows a zoomed-in 10 min subinterval of burst mode measurements from the boundary layer (10:20–10:30 UTC, shaded gray in Figure 2). By comparing measurements of the ion temperature (Figure 3d) and the ion heat flux magnitude (Figure 3f), we clearly see that sharp peaks in the heat flux coincide with sudden transitions between hot boundary layer plasma and cold magnetosheath plasma at K-H wave boundary crossings. Figure 3f also illustrates the difference in the heat flux magnitude when calculated from the VDF with one-count bins included (in gray) and excluded (in black). The former, which is available as a ready data product provided by the MMS team, has a reported uncertainty of <20% at the peaks, low enough to confirm the peaks as robust features. The magnitude is lower in the latter case, since many high-energy, low count bins are
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Zoomed-in full-resolution MMS1 burst mode data from the gray highlighted interval in Figure 2. For the explanation of each panel (a–h), see the caption of Figure 2. Here, panel (f) includes a comparison of the ion heat flux magnitude calculated from the ion VDF with one-count bins included (gray line) and excluded (black line). The blue dashed (dotted) horizontal line shows the mean (mean + two standard deviations) value of the heat flux magnitude (with one-count bins excluded) over the full burst data interval (10:07:04–11:27:34 UTC). The black arrow in panel (f) denotes the heat flux peak analyzed in Figures 5 and 6. Panel (g) includes estimates of $|q_{\Lambda,i}|$ (Equation 4) at select wave boundary crossings (see text for details).</figcaption></figure>

excluded, highlighting the sensitivity of heat flux measurements to low particle count statistics at high energies. We will use the heat flux data calculated excluding the one-count bins for the rest of this work.

Figure 3g, comparing the parallel and transverse heat flux components, shows also the estimated transverse heat flux magnitude $|q_{\Lambda,i}| = \kappa_{\Lambda,i}|\nabla_{\perp}T|$ at select wave boundary crossings, which have been chosen using criteria detailed later in Section 3.2. The value of $\kappa_{\Lambda,i}$ is shown in Figure 3h, with an uncertainty of <5%. $|\nabla_{\perp}T|$ is approximated as the transverse part of ~$v_n^{-1}d\langle T \rangle/dt \hat{n}$, where $v_n$ and $\hat{n}$ are the boundary normal speed and normal vector obtained via the timing analysis of the ion temperature at the crossings. In the interval of Figure 3, we find $|\nabla_{\perp}T| \sim 0.3–1.7 eV km^{-1}$ at the leading wave boundary crossings, and $|\nabla_{\perp}T| \sim 1–2.6 eV km^{-1}$ at the trailing wave boundary crossings. Assuming temperature jumps of ~1000 eV (Figure 3d), the values yield boundary thicknesses of ~390–3330 km (the mean thermal proton Larmor radius over the burst data interval is 68 km), which are comparable to the current sheet thicknesses (~30–1330 km) reported by Eriksson, Lavraud, et al. (2016) during the event. We thus consider our estimates of $|\nabla_{\perp}T|$ to be representative of the magnetopause conditions.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** 2D histograms of MMS1 ion burst mode measurements during the interval of available ion burst data (10:07:04–11:27:34 UTC). The y-axis shows the heat flux magnitude. Panel (a): the x-axis shows the angle $\theta_{qb}$ between the heat flux vector and the magnetic field vector (resampled to the ion data cadence). Panel (b): the x-axis shows the time derivative $d\langle T\rangle/dt$ of the temporally averaged temperature (see Section 2.2 for details). Horizontal dashed (dotted) lines show the mean (mean + two standard deviations) value of the heat flux magnitude.</figcaption></figure>

The obtained values of $|q_{\Lambda,i}| \sim 0.03 - 0.13 \, \text{mW m}^{-2}$ are well below the measured heat flux peaks. The discrepancy is likely explained by the LLBL ion distribution not being a single $\sim 1000 \, \text{eV}$ distribution, but a mix of cold magnetosheath ($\sim 200 \, \text{eV}$) and hot magnetospheric ($\sim 10,000 \, \text{eV}$) ions (Figures 2a–2d). As Equation 4 only accounts for the combined temperature, the contribution of the hot ions is underestimated in $|q_{\Lambda,i}|$. Even a small fraction of hot ions in the LLBL could greatly enhance the heat flux, which is seen by comparing the thermal energy flux $nV_{th}k_B T$ in the ambient magnetosphere ($\sim 3.1 \, \text{mW m}^{-2}$, $n_i \sim 2 \, \text{cm}^{-3}$ from Figure 2b) and magnetosheath ($\sim 0.1 \, \text{mW m}^{-2}$, $n_i \sim 25 \, \text{cm}^{-3}$). For the earlier estimated boundary thicknesses, a layer with ambient magnetosheath ions at a pressure balance with magnetospheric ions ($n_i \sim 0.5 \, \text{mW m}^{-2}$) would have a $|q_{\Lambda,i}| \sim 0.09 - 0.75 \, \text{mW m}^{-2}$, a closer match with the measured peaks. With this consideration, we conclude that FLR effects are a plausible source for the transverse heat flux enhancements.

Figure 4 provides a statistical view of the burst mode ion heat flux data available during the event (10:07:04–11:27:34 UTC, yellow highlight in Figure 2) via two 2D histograms. Figure 4a illustrates the relationship between the heat flux angle relative to the magnetic field (x-axis) and the heat flux magnitude (y-axis). A strong transverse heat flux component is evident at large heat flux magnitudes. Angles between 45° and 135° account for the largest overall magnitudes, and for the majority of data points ($\sim 65\%$) at magnitudes exceeding two standard deviations above the mean ($>0.24 \, \text{mW m}^{-2}$, horizontal dotted line). As such magnitudes are chiefly found at distinct peaks (Figure 3f), we see that the peaks primarily consist of transverse enhancements.

Figure 4b illustrates the relationship between the time derivative $d\langle T\rangle/dt$ of the temporally averaged ion temperature (x-axis) and the heat flux magnitude (y-axis). Large heat flux magnitudes occur more commonly for negative $d\langle T\rangle/dt$ than for positive $d\langle T\rangle/dt$. At magnitudes exceeding two standard deviations above the mean (horizontal dotted line), $d\langle T\rangle/dt < 0$ accounts for $\sim 65\%$ of the data points. The difference indicates that the heat flux peaks occur preferentially at the trailing part of a K-H wave (sharp temperature drops in Figure 3d) as opposed to the leading part (gradual temperature increases). This is consistent with the expectation of larger temperature gradients generating a larger heat flux (Equation 4). We further compare the heat flux at trailing/leading wave parts in the next section, where we perform a statistical analysis of individual heat flux peaks.

## 3.2. Analysis of Wave Boundary Crossings

In this section, we perform a statistical analysis of heat flux peaks at individual K-H wave boundary crossings. By evaluating $\hat{b} \times \hat{n}$ at each boundary crossing (as outlined in Section 2.2), and comparing its direction with that of the measured heat flux, we evaluate whether the heat flux direction at the peaks is consistent with Equation 4 and the FLR generation mechanism described by Braginskii (1965). In Figure 5, we show 2D projections of an ion
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Projections of the ion velocity distribution function at a trailing K-H wave boundary crossing (black arrow in Figure 3f). A white cross shows the ion bulk velocity, while the arrows show the projections of the heat flux vector $\hat{q}$, magnetic field vector $\hat{b}$, and $\hat{b} \times \hat{n}_T$. The vectors are normalized to have equal total length.</figcaption></figure>

VDF from an example trailing wave boundary crossing at *t* = 10:27:21 UTC (black arrow in Figure 3f), together with projections of the measured heat flux vector $\hat{q}$ (white arrow), magnetic field vector $\hat{b}$ (pink arrow) and $\hat{b} \times \hat{n}_T$ (black arrow), where $\hat{n}_T$ is the boundary normal vector estimated via the timing analysis of ion temperature data.

As readily seen in Figure 5a ($v_{x'} - v_{y'}$ plane), the VDF exhibits a clear asymmetry around the bulk velocity $\vec{v}_{bulk} = 255 (0.13, 0.98, 0.13) \text{ km s}^{-1}$, (white cross), due to a net flux of hot ions in the $(-v_{x'}$, $-v_{y'}$) -direction. Consequently, a heat flux exists in the same direction. Comparing with Figures 5b and 5c ($v_{y'} - v_{z'}$ plane and

$v_{x'} - v_{z'}$ plane, respectively), the asymmetry and the heat flux appear strongest transverse to $\hat{b}$. The direction of $\hat{q}$ is in close agreement with $\hat{b} \times \hat{n}_T$, confirming that the heat flux is consistent with the mechanism described by Braginskii (1965). We have verified that our estimate of $\hat{b} \times \hat{n}_T$ is reliable by comparing it with different estimates of $\hat{b} \times \hat{n}$ obtained via the timing analysis of other quantities. The azimuthal angle $\phi$ for different estimates is visualized in Figure 6 (within ~9° of each other). We also note that the estimated direction of $\hat{b} \times \hat{n}$ is consistent with our depiction of a trailing K-H wave boundary in Figure 1.

We next assess the temporal variation of the heat flux direction across our example boundary crossing (not shown). To do this, we define the heat flux peak duration as the interval where the heat flux magnitude remains at ≥50% of the peak magnitude. At the example boundary crossing, the peak spans 10:27:19.3–10:27:22.7 UTC, with a total duration of 3.5 s. During this interval, the relative angle between the heat flux vector and the magnetic field vector (resampled to the ion data cadence) ranges between [75°, 104°] (mean 87°), revealing a consistent transverse heat flux across the peak.

The azimuthal angle $\phi$ of the measured heat flux vector varies between [-132°, -117°] (mean -127°), showing that the direction is more or less constant across the peak. The value of $\phi$ for $\hat{b} \times \hat{n}_T$ varies between [-138°, -131°] (mean -134°), while the difference $\Delta\phi$ between the measured heat flux vector and $\hat{b} \times \hat{n}_T$ varies between [0°, 15°] (mean 7°), confirming that the two vectors remain consistent across the peak.

Next, we perform a statistical analysis of the heat flux peaks found throughout the ion burst data interval (10:07:04–11:27:34 UTC). We use the following criteria to select our data set:

<figure><img src="image_8.png" /><figcaption>**Figure 6.** The azimuthal angle $\phi$ at the time shown in Figure 5 for the measured heat flux vector (white arrow) and different estimates of $\hat{b} \times \hat{n}$, where $\hat{n}$ is obtained via the timing analysis of ion temperature ($\hat{n}_T$, black arrow), ion density ($\hat{n}_i$, blue arrow) and the $x'$, $y'$-components of the magnetic field ($\hat{n}_{Bx'}$, $\hat{n}_{By'}$, green and red arrows, respectively).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Statistics from individual heat flux peaks during the ion burst data interval (10:27:19.3–10:27:22.7 UTC, yellow highlight in Figure 2). Panel (a): maximum heat flux magnitude. Panel (b): the angle $\theta_{qb}$ between the heat flux vector and the magnetic field vector. Panels (c) and (d): the azimuthal angle $\phi$ for the measured heat flux vector and $\hat{b} \times \hat{n}_T$, respectively. In panels (b)–(d), dots show the angle at the time of the maximum heat flux magnitude, while the vertical bars show the variation across the peak. A star marks the peak at the example boundary crossing analyzed in Figures 5 and 6.</figcaption></figure>

1. Only peaks at which the heat flux magnitude reaches at least two standard deviations above the mean magnitude of the burst data interval ($|\vec{q}_{peak}| \ge 0.24 \, \text{mW m}^{-2}$) are selected. Such peaks are unambiguously identified as distinct structures in the heat flux data (Figure 3f, horizontal dotted blue line).

2. The duration of a peak encompasses the interval where the heat flux magnitude is $\ge 50\%$ of the local maximum.

3. If the considered interval contains any larger local maxima, the peak is not selected to avoid overlapping peaks.

With these criteria, we obtain a data set of 100 heat flux peaks. We classify the heat flux peaks according to the sign of $d\langle T \rangle/dt$, in order to study the differences between the trailing wave boundaries (decreasing temperature) and leading wave boundaries (increasing temperature). We find 36 (29) peaks for which $d\langle T \rangle/dt < 0$ ($d\langle T \rangle/dt \ge 0$) over the peak duration, and 35 peaks across which $d\langle T \rangle/dt$ changes sign. We retain the first two categories for further analysis, identifying them with trailing and leading wave boundary crossings, respectively, and calculate $\hat{b} \times \hat{n}_T$ via the timing analysis of ion temperature data at each crossing. We select the peaks where the timing analysis yields a satisfactory result (as outlined in Section 2.2) for our final data set.

Our final data set contains 48 heat flux peaks, of which 33 (15) are associated with a boundary crossing with $d\langle T \rangle/dt < 0$ ($d\langle T \rangle/dt \ge 0$). We provide further analysis of these peaks in Figure 7. In Figure 7a, showing the peak magnitude, we readily see that such peaks are a common feature across the ion burst data interval (10:07:04–
<figure><img src="image_10.png" /></figure>

**Figure 8.** The azimuthal angle $\Delta\phi$ relative to $\hat{b} \times \hat{n}_T$ for the measured heat flux vector (thin colored vertical bars) and the range of different estimates of $\hat{b} \times \hat{n}$ (via the timing analysis of ion density and magnetic field ($x', y'$)-components, thick white vertical bars) over each heat flux peak shown in Figure 7. For the former, a dot shows $\Delta\phi$ at the time of maximum heat flux, and the latter is shown only for peaks where all estimates of $\hat{b} \times \hat{n}$ yield a satisfactory result. Panel (a) shows the peaks with $d\langle T \rangle/dt < 0$, and panel (b) those with $d\langle T \rangle/dt \ge 0$. Peaks at which the measured heat flux vector is within $|\Delta\phi| < 90^\circ$ for less than 50% of the peak duration are highlighted with a dotted line.

11:27:34 UTC). The mean separation between the peaks at the trailing and leading wave boundary crossings are 125 and 265 s, respectively. With a K-H wave period of 63 s (Eriksson, Lavraud, et al., 2016), the corresponding peak occurrence rates are thus 50% and 24%, respectively. The peak magnitude tends to be larger at boundary crossings with $d\langle T \rangle/dt < 0$ (mean and standard deviation $0.37 \pm 0.14 \text{ mW m}^{-2}$) than with $d\langle T \rangle/dt$ (mean and standard deviation $0.27 \pm 0.04 \text{ mW m}^{-2}$). This is consistent with the expectation of larger heat flux magnitudes occurring at the trailing wave boundaries, which display steep temperature gradients, as opposed to the leading wave boundaries with gentler temperature gradients (Figure 3d).

In Figure 7b, we show the angle of the heat flux vector relative to the magnetic field throughout each heat flux peak. The dots show the instantaneous angle at the time of the maximum heat flux, while the vertical bars show the variation throughout the full peak duration. For $d\langle T \rangle/dt < 0$, the angle is predominantly between $45^\circ$ and $135^\circ$ (mean and standard deviation $65^\circ \pm 27^\circ$), indicating that the peaks typically consist of transverse heat flux enhancements. For $d\langle T \rangle/dt \ge 0$, a strong parallel component is evident (angle mean and standard deviation $25^\circ \pm 25^\circ$). In 7 out of the 15 cases, the ion VDF contains a distinct field-aligned beam (not shown), suggesting that FLR mixing is not the dominant heat flux generating mechanism.

In Figures 7c and 7d, we show the azimuthal angle $\phi$ throughout each peak for the measured heat flux vector and $\hat{b} \times \hat{n}_T$, respectively. For $d\langle T \rangle/dt < 0$, both vectors commonly show a distinct, constant angle near $\phi \sim -135^\circ$, indicating recurring encounters with trailing wave boundaries inclined relative to the magnetosheath flow, similar to the example boundary analyzed in Figures 5 and 6. For $d\langle T \rangle/dt \ge 0$, $\phi$ is much more variable, but the value range ($\phi > -135^\circ$) clearly differs from that found for $d\langle T \rangle/dt < 0$. We will present an interpretation of the observed $\phi$-distribution through comparison with our prior K-H simulation (Paper I) in Section 3.3.

Lastly, in Figure 8 we visualize the azimuthal angle difference $\Delta\phi$ relative to $\hat{b} \times \hat{n}_T$ at each peak for the measured heat flux vector (colored vertical bars) and for the different estimates of $\hat{b} \times \hat{n}$ (white vertical bars, shown when all estimates yield a satisfactory result). To find if the measured heat flux vector is preferentially parallel (or anti-parallel) to $\hat{b} \times \hat{n}_T$, we have marked the peaks according to whether the measured heat flux vector is within $|\Delta\phi| < 90^\circ$ for more/less than 50% of the peak duration (solid/dotted vertical lines, respectively). Figure 8a, showing the peaks with $d\langle T \rangle/dt < 0$, indicates a clear parallel trend. Only at 3 out of 33 peaks, the heat flux vector is preferentially anti-parallel. We also calculate (not shown) that on average, the heat flux vector contains a component parallel to $\hat{b} \times \hat{n}_T$ 90% of the peak duration. At the peaks where all estimates of $\hat{b} \times \hat{n}$ yield a satisfactory result, generally close agreement is found with $\hat{b} \times \hat{n}_T$, indicating that the uncertainty in the direction $\hat{b} \times \hat{n}_T$ is minor. In Figure 8b, showing the peaks with $d\langle T \rangle/dt \ge 0$, the agreement with $\hat{b} \times \hat{n}_T$ is
weaker, but a parallel trend still exists (preferentially anti-parallel heat flux vector at 3 out of 15 peaks, parallel component present 72% of peak duration on average). In contrast to the peaks with $d\langle T\rangle/dt < 0$, there are no peaks where all estimates of $\hat{b} \times \hat{n}$ yield a satisfactory result. This suggests that the leading wave boundaries tend to have more irregular structure, which may explain the ambiguous heat flux signatures. In conclusion, the prevalence of a heat flux component parallel to $\hat{b} \times \hat{n}_T$ at both trailing and leading wave boundary crossings suggests that active ion FLR mixing is encountered during both types of crossings.

### 3.3. Kelvin-Helmholtz Wave Structure and Evolution

Lastly, we compare our results with a local 2-D hybrid-Vlasov simulation of the K-H instability (Paper I) in order to provide context for the in situ heat flux measurements. We investigate the link between the heat flux direction and K-H wave structure, and evaluate if the $\phi$-distribution found in Section 3.2 can be explained by the existence of a magnetopause FLR mixing layer created by K-H waves. We use the simulation to show how the $\phi$-distribution of the heat flux vector evolves in time at the FLR mixing layer, and look for corresponding patterns in the in situ data. As a reference for the K-H wave temporal evolution in the in situ data, we refer to the study of Settino et al. (2022), who provided a characterization of K-H wave structure during the event based on the spatial distribution of mixed electrons and ions within the waves.

The simulation (full description in Section 2 of Paper I), uses a 2D domain with two counter-streaming plasma flows (representing the magnetosphere and the magnetosheath, respectively) with a velocity shear at the center of the domain. A constant magnetic field is applied transverse to the flow. The coordinate system is analogous to the one in our present study (Section 2.2). The simulation parameters have been chosen to be representative of the present MMS event (see last row in Table 1 of Paper I).

In Figure 9, we present snapshots from the simulation that depict different stages of the K-H instability. The first column shows a colormap of the ion temperature, with an overlay identifying the area where $|\nabla_{\perp} T| > \max(|\nabla_{\perp} T|(t = 0))$. The area is colored based on the value of the y-component of $\nabla_{\perp} T$, which we use as a proxy for $d\langle T\rangle/dt$ in the MMS data. $dT/dy > 0$ is analogous to $d\langle T\rangle/dt < 0$, and vice versa. In Paper I, we showed that this area, corresponding to the FLR mixing layer formed by the K-H wave steepening, contains an enhanced heat flux (median magnitude ~0.3 – 0.5 mW m⁻², directed within $\Delta\phi < 15°$ of $\hat{b} \times \nabla_{\perp} T$). The second column shows the azimuthal angle $\phi$ of the heat flux vector in the selected area. The third and fourth columns show the distribution of $\phi$ in the areas with $dT/dy > 0$ and $dT/dy \le 0$, respectively.

The first three rows of Figure 9 depict the basic K-H wave types following the classification by Settino et al. (2022). Type I wave (first row) is a mostly linear surface wave, Type II (second row) a steepened non-linear wave, and Type III (third row) a rolled-up vortex. For completeness, the fourth row depicts the final stage of the K-H instability: a highly perturbed vortex of mixed plasma. Settino et al. (2022) showed that changing solar wind conditions during the present event resulted in different wave types being observed. The beginning of the event has the highest proportion of Type I waves, Type II waves are found steadily throughout the event, and the end of the event has the highest proportion of Type III waves. Thus, we might observe a corresponding temporal evolution of the $\phi$-distribution of the heat flux vector, reflecting these changes.

The Type I wave exhibits a narrow $\phi$-distribution for the heat flux vector. A trailing wave boundary inclined relative to the y-axis is evident from the narrow $\phi$-distribution for $dT/dy > 0$ (–135° < $\phi$ < –90°). A narrow $\phi$-distribution for $dT/dy \le 0$ (–90° < $\phi$ < –45°) indicates a weakly developed leading wave boundary. The Type II wave, in contrast, exhibits an extended $\phi$-distribution for $dT/dy \le 0$ (–90° < $\phi$ < 90°) due to a steepened leading wave boundary. The Type III wave $\phi$-distribution is similar to that of the Type II wave, but the $\phi$-distribution for $dT/dy \le 0$ is pronounced between 45° < $\phi$ < 90° due to the leading wave boundary folding upon itself. The final stage (fourth row) is characterized by an irregular wave structure following the growth of secondary K-H and Rayleigh-Taylor instabilities (Matsumoto & Hoshino, 2004). The $\phi$-distribution shares the main characteristics of the Type II and III waves, but contains anomalous features (e.g., the bifurcation between –135° < $\phi$ < –90° for $dT/dy > 0$) due to the irregularities.

In Figure 10, we show the $\phi$-distribution for the MMS heat flux vector over the full ion burst data interval (10:07:04–11:27:34 UTC), grouped according to the sign of $d\langle T\rangle/dt$. The distribution includes all heat flux data with a magnitude above two standard deviations from the mean magnitude (≥0.24 mW m⁻²), our threshold for the
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Snapshots of the 2-D hybrid-Vlasov simulation. First column: ion temperature *T*, and the area where $|∇T_⊥| > \max(|∇⊥(t=0)|)$, colored with the value of *dT/dy*. Second column: azimuthal angle $\phi$ for the heat flux vector in the colored area. Third and fourth columns: the distribution of $\phi$ in the parts of the colored area with $dT/dy > 0$ and $dT/dy \le 0$, respectively. The first three rows are labeled according to K-H wave type classification by Settino et al. (2022). The times are in the units of proton gyroperiod, and lengths relative to the initial shear layer thickness ($a = 220$ km).</figcaption></figure>
heat flux peaks studied in Section 3.2. The first row shows the $\phi$-distribution for the full interval, while the three rows below show the data divided into three sub-intervals of equal duration spanning 26 min and 50 s each.

The main features in the $\phi$-distribution compiled over the full interval are in agreement with the K-H wave features identified at all times in Figure 9. These include an enhancement between $-135^\circ < \phi < -90^\circ$ for $d\langle T \rangle/dt < 0$ (analogous to $dT/dy > 0$ in Figure 9), which is consistent with the trailing wave boundary, and an enhancement between $-90^\circ < \phi < -45^\circ$ for $d\langle T \rangle/dt \ge 0$ (analogous to $dT/dy \le 0$ in Figure 9), which is consistent with the leading wave boundary. By comparing the first two subintervals of the event (second and third rows of Figure 10, respectively), we find evidence suggesting K-H wave temporal evolution. The $\phi$-distribution for $d\langle T \rangle/dt < 0$ remains largely unchanged between the subintervals, while the $\phi$-distribution for $d\langle T \rangle/dt \ge 0$ evolves from a Type II-like distribution to a Type III-like distribution which includes an enhancement between $45^\circ < \phi < 90^\circ$. The evolution is consistent with the findings of Settino et al. (2022), who demonstrated that the number of Type III waves increases toward the end of the event.

The third subinterval (fourth row in Figure 10) contains only a few data points, likely due to MMS moving into the magnetosheath, where the temperature (and the heat flux) is lower. A possible reason for the irregular $\phi$-distribution (apart from a small sample size) could be the presence of highly evolved K-H waves (Figure 9m), which contain an irregular magnetosheath-side boundary, with a relatively small temperature gradient. Based on these findings, we conclude that the heat flux patterns observed during the event are consistent with the patterns exhibited by the FLR mixing layer in our simulation, and with the K-H wave temporal evolution described by Settino et al. (2022).

## 4. Discussion and Conclusions

In this work, we have demonstrated the presence of ion heat flux enhancements during magnetopause K-H wave boundary crossings by the MMS spacecraft, which took place at the dusk flank magnetopause under northward IMF conditions on 8 September 2015. The heat flux enhancements contain a significant component transverse to the magnetic field. This transverse heat flux is consistent with being generated by FLR ion mixing across a temperature gradient, in accordance with the mechanism described by Braginskii (1965). The measured heat flux enhancements are in the range ~0.2 – 0.9 mW m⁻², while the heat flux formula of Braginskii (1965) suggests only values comparable to the ambient heat flux fluctuations (0.14 ± 0.11mW m⁻²). The predicted heat flux magnitude is likely underestimated, due to the formula only accounting for the measured ion temperature, while the ion distribution at the low-latitude boundary layer is a mix of magnetosheath ions and magnetospheric ions with different temperatures. By separately considering the contribution of the hot magnetospheric ions, the predicted heat flux magnitude given by the formula becomes comparable to the measured enhancements.

By analyzing individual K-H wave boundary crossings over the event, we have shown that the transverse heat flux is preferentially along the direction given by $\hat{b} \times \nabla_{\perp} T$ as predicted by Braginskii (1965), where $\hat{b}$ is the magnetic field vector, and $\nabla_{\perp} T$ the temperature gradient transverse to $\hat{b}$, whose direction we have estimated via four-spacecraft timing analysis. We have identified differences between trailing and leading wave boundary crossings by classifying the crossings according to the sign of the time derivative of the temporally averaged ion temperature $d\langle T \rangle/t$. The trailing wave boundaries (negative $d\langle T \rangle/dt$) exhibit clearer transverse heat flux signatures than the leading wave boundaries (positive $d\langle T \rangle/dt$), which is consistent with the former exhibiting stronger temperature gradients. The heat flux direction is constant between trailing wave boundaries, while variation exists between leading wave boundaries.

By using our past hybrid-Vlasov simulation of the K-H instability (Tarvus et al., 2024) as a reference, we have shown that the observed behavior of the heat flux direction at the trailing and leading wave boundaries is consistent with the presence of an ion FLR mixing layer along the K-H wave boundaries. Such mixing layer is formed when the K-H wave growth steepens the magnetopause to a thickness comparable to an ion Larmor radius, and it has been previously identified in kinetic simulations of the K-H instability (Ma et al., 2019; Terasawa et al., 1992). Based on a classification of K-H wave structure and evolution during the present event (Settino et al., 2022), we have identified temporal changes in the variation of the heat flux direction, which are consistent with the FLR mixing layer growing in length as the K-H waves roll up into vortices.

Based on these findings, we conclude that the transverse ion heat flux enhancements found at the magnetopause K-H wave boundary crossings are likely caused by FLR ion mixing, which occurs at a mixing layer similar to that
Full Interval, 10:07:04-11:27:34 UTC

<figure><img src="image_12.png" /></figure>

10:07:04-10:33:54 UTC

<figure><img src="image_13.png" /></figure>

10:33:54-11:00:44 UTC

<figure><img src="image_14.png" /></figure>

11:00:44-11:27:34 UTC

<figure><img src="image_15.png" /></figure>

**Figure 10.** The distribution of the azimuthal angle $\phi$ for the heat flux vector in the MMS burst mode data, categorized according to $d\langle T \rangle/dt$. The data set includes all data with a heat flux magnitude above two standard deviations from the mean magnitude. The number of data points *N* is shown in each panel.
Acknowledgments

Data analysis was performed using the IRFU-Matlab and pyrfu analysis packages available at https://github.com/irfu/irfu-matlab, and https://github.com/louis-richard/irfu-python, respectively. The work of VT and LT is supported by the Research Council of Finland (Grant Number 322544, 328893 and 353197) and by the European Union (ERC Grant WAVESTORMS – 101124500). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. GC is supported by the Integration Fellowship of Le Studium Loire Valley Institute for Advanced Studies. YP acknowledges Research Council of Finland Grant Number 339756. MP acknowledges the Research Council of Finland Grant Number 361901. The work of JS was made possible by a doctoral researcher position at the Doctoral Programme in Particle Physics and Universe Sciences funded by the University of Helsinki. MP and VT acknowledge the Research Council of Finland for the Finnish Centre of Excellence in Research of Sustainable Space (Grant Number 352846). We acknowledge the European Research Council for Starting Grant 200141-QuESpace, with which Vlasiator was developed, and Consolidator Grant 682068-PRESTISSIMO, awarded to further develop Vlasiator and use it for scientific investigations. We acknowledge the CSC – IT Center for Science computational project 2004873 (WAVESTORMS), with which the simulation run used in this study was produced. The authors thank the Finnish Computing Competence Infrastructure (FCCI), the Finnish Grid and Cloud Infrastructure (FGCI) and the University of Helsinki IT4SCI team for supporting this project with computational and data storage resources. We express our gratitude to the research group at Space Research Institute Graz for hosting a research visit in March 2023, during which the initial work leading to the present results was conducted. We also express our gratitude to Dr. Takuma Nakamura for his role in arranging the research visit and in providing advice throughout the work. Open access publishing facilitated by Helsingin yliopisto, as part of the Wiley - FinELib agreement.

reported in kinetic simulations (Ma et al., 2019; Terasawa et al., 1992). We underline that such mixing layer can provide an important source of magnetosheath plasma transport into the magnetosphere when coupled with turbulence (Cowee et al., 2009; Matsumoto & Hoshino, 2004; Matsumoto & Seki, 2010). While K-H wave growth causes the FLR mixing layer to grow in length, turbulence causes the layer to also widen, leading to efficient diffusive transport of mixed plasma. Our present evidence of an FLR mixing layer, together with past reports of K-H driven turbulence (Stawarz et al., 2016), indicate that such transport process could play an important role in situ.

We also highlight the close relationship between K-H wave structure and the variation of the transverse heat flux direction. We propose that our method for evaluating the variation introduced in Section 3.3 (based on comparing the distributions of the heat flux direction between the different signs of $d\langle T \rangle/dt$) could be used for identifying K-H waves at different evolutionary stages, complementing existing methods (Settino et al., 2022). For reliable results, the measurements should cover the K-H wave interior, where the differences between the wave types are most apparent. Thus, the best option may be to consider a statistical sample over a full transit of the K-H active layer (i.e., the magnetopause). Care should be taken when interpreting partial transits by for example, comparing with alternative available methods, as done here. We also note that our method, which is based on the analytical results of Braginskii (1965), is only valid for an approximately transverse configuration (magnetic field transverse to the bulk flow) without a magnetic shear. This configuration however, corresponding to the low-latitude magnetopause under northward IMF conditions, is when magnetopause K-H waves are most common (Kavosi & Raeder, 2015). More complex settings warrant further numerical work to understand the resulting kinetic signatures. Future works could utilize an existing model for a velocity shear at a kinetic equilibrium state, such as the configurable hybrid-Vlasov equilibrium recently introduced by Malara et al. (2022).

Lastly, we note that the heat flux signature investigated here may be weaker at the dawn flank magnetopause due to another FLR effect, namely a temperature anisotropy due to the velocity shear (Braginskii, 1965). This anisotropy is dependent on the relative orientation between the shear and the magnetic field. Previously, in Paper I, we found the temperature to increase (decrease) along $\hat{b} \times \nabla_{\perp} T$ when the vorticity is aligned (opposite) to the magnetic field. This increases (decreases) the heat flux in the dusk (dawn) flank geometry, which should be taken into account when comparing the transverse heat flux and the associated FLR mixing between the flanks.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

The MMS data analyzed in this work is freely available at MMS science data center (https://lasp.colorado.edu/mms/sdc/public/about/browsewrapper). The Vlasiator simulation data used in this work is freely available at https://doi.org/10.5281/zenodo.17407342. The Vlasiator source code is freely available at https://github.com/fmihpc/vlasiator.

References

Axford, W. I. (1964). Viscous interaction between the solar wind and the Earth's magnetosphere. *Planetary and Space Science*, *12*(1), 45-53.
https://doi.org/10.1016/0032-0633(64)90067-4

Axford, W. I., & Hines, C. O. (1961). A unifying theory of high-latitude geophysical phenomena and geomagnetic storms. *Canadian Journal of Physics*, *39*(10), 1433-1464. https://doi.org/10.1139/p61-172

Blasl, K. A., Nakamura, T. K. M., Nakamura, R., Settino, A., Hasegawa, H., Vörös, Z., et al. (2023). Electron-scale reconnecting current sheet formed within the lower-hybrid wave-active region of kelvin-Helmholtz waves. *Geophysical Research Letters*, *50*(19), e2023GL104309.
https://doi.org/10.1029/2023GL104309

Blasl, K. A., Nakamura, T. K. M., Plaschke, F., Nakamura, R., Hasegawa, H., Stawarz, J. E., et al. (2022). Multi-scale observations of the magnetopause Kelvin-Helmholtz waves during southward IMF. *Physics of Plasmas*, *29*(1), 012105. https://doi.org/10.1063/5.0067370

Blasl, K. A., Settino, A., Nakamura, R., Hasegawa, H., Nakamura, T. K. M., & Hosner, M. (2024). Lower-Hybrid wave-induced plasma mixing related to kelvin-helmholtz vortices during Southward IMF. *Journal of Geophysical Research (Space Physics)*, *129*(11), 2024JA033152.
https://doi.org/10.1029/2024JA033152

Braginskii, S. I. (1965). Transport processes in a plasma. *Reviews of plasma physics*, *1*, 205.

Chandrasekhar, S. (1961). Hydrodynamic and hydromagnetic stability.

Chaston, C., Bonnell, J., McFadden, J. P., Carlson, C. W., Cully, C., Le Contel, O., et al. (2008). Turbulent heating and cross-field transport near the magnetopause from THEMIS. *Geophysical Research Letters*, *35*(17), L17S08. https://doi.org/10.1029/2008GL033601
Chaston, C., Wilber, M., Mozer, F. S., Fujimoto, M., Goldstein, M. L., Acuna, M., et al. (2007). Mode conversion and anomalous transport in kelvin-helmholtz vortices and kinetic Alfvén waves at the Earth’s magnetopause. *Physical Review Letters*, **99**(17), 175004. https://doi.org/10.1103/PhysRevLett.99.175004

Cowee, M. M., Winske, D., & Gary, S. P. (2009). Two-dimensional hybrid simulations of superdiffusion at the magnetopause driven by Kelvin-Helmholtz instability. *Journal of Geophysical Research (Space Physics)*, **114**(A10), A10209. https://doi.org/10.1029/2009JA014222

Eriksson, S., Lavraud, B., Wilder, F. D., Stawarz, J. E., Giles, B. L., Burch, J. L., et al. (2016). Magnetospheric Multiscale observations of magnetic reconnection associated with Kelvin-Helmholtz waves. *Geophysical Research Letters*, **43**(11), 5606–5615. https://doi.org/10.1002/2016GL068783

Eriksson, S., Wilder, F. D., Ergun, R. E., Schwartz, S. J., Cassak, P. A., Burch, J. L., et al. (2016). Magnetospheric multiscale observations of the electron diffusion Region of large guide field magnetic reconnection. *Physical Review Letters*, **117**(1), 015001. https://doi.org/10.1103/PhysRevLett.117.015001

Faganello, M., Califano, F., Pegoraro, F., & Andreussi, T. (2012). Double mid-latitude dynamical reconnection at the magnetopause: An efficient mechanism allowing solar wind to enter the Earth’s magnetosphere. *EPL*, **100**(6), 69001. https://doi.org/10.1209/0295-5075/100/69001

Harvey, C. C. (1998). Spatial gradients and the volumetric tensor. *ISSI Scientific Reports Series, 1*, 307–322.

Hasegawa, H., Fujimoto, M., Phan, T. D., Rème, H., Balogh, A., Dunlop, M. W., et al. (2004). Transport of solar wind into Earth’s magnetosphere through rolled-up Kelvin-Helmholtz vortices. *Nature*, **430**(7001), 755–758. https://doi.org/10.1038/nature02799

Hasegawa, H., Nakamura, T. K. M., Gershman, D. Å. J., Nariyuki, Y., Viñas, A. Å. F., Giles, B. Å. L., et al. (2020). Generation of turbulence in kelvin-helmholtz vortices at the Earth’s magnetopause: Magnetospheric multiscale observations. *Journal of Geophysical Research (Space Physics)*, **125**(3), e27595. https://doi.org/10.1029/2019JA027595

Hasegawa, H., Retinò, A., Vaivads, A., Khotyaintsev, Y., André, M., Nakamura, T. K. M., et al. (2009). Kelvin-Helmholtz waves at the Earth’s magnetopause: Multiscale development and associated reconnection. *Journal of Geophysical Research (Space Physics)*, **114**(A12), A12207. https://doi.org/10.1029/2009JA014042

Hwang, K. J., Dokgo, K., Choi, E., Burch, J. L., Sibeck, D. G., Giles, B. L., et al. (2020). Magnetic reconnection inside a flux rope induced by kelvin-helmholtz vortices. *Journal of Geophysical Research (Space Physics)*, **125**(4), e27665. https://doi.org/10.1029/2019JA027665

Hwang, K. J., Kuznetsova, M. M., Sahraoui, F., Goldstein, M. L., Lee, E., & Parks, G. K. (2011). Kelvin-Helmholtz waves under southward interplanetary magnetic field. *Journal of Geophysical Research (Space Physics)*, **116**(A8), A08210. https://doi.org/10.1029/2011JA016596

Hwang, K.-J., Wang, C.-P., Nykyri, K., Hasegawa, H., Tapley, M. B., Burch, J. L., et al. (2023). Kelvin-Helmholtz instability-driven magnetopause dynamics as turbulent pathway for the solar wind-magnetosphere coupling and the flank-central plasma sheet communication. *Frontiers in Astronomy and Space Sciences*, **10**, 1151869. https://doi.org/10.3389/fspas.2023.1151869

Kavosi, S., & Raeder, J. (2015). Ubiquity of Kelvin-Helmholtz waves at Earth’s magnetopause. *Nature Communications*, **6**(1), 7019. https://doi.org/10.1038/ncomms8019

Kieokaew, R., Lavraud, B., Foullon, C., Toledo-Redondo, S., Fargette, N., Hwang, K. J., et al. (2020). Magnetic reconnection inside a flux transfer event-like structure in magnetopause Kelvin-Helmholtz waves. *Journal of Geophysical Research (Space Physics)*, **125**(6), e27527. https://doi.org/10.1029/2019ja027527

Li, W., André, M., Khotyaintsev, Y. V., Vaivads, A., Graham, D. B., Toledo-Redondo, S., et al. (2016). Kinetic evidence of magnetic reconnection due to Kelvin-Helmholtz waves. *Geophysical Research Letters*, **43**(11), 5635–5643. https://doi.org/10.1002/2016GL069192

Ma, X., Delamere, P. A., Nykyri, K., Burkholder, B., Neupane, B., & Rice, R. C. (2019). Comparison between fluid simulation with Test particles and hybrid simulation for the kelvin-helmholtz instability. *Journal of Geophysical Research (Space Physics)*, **124**(8), 6654–6668. https://doi.org/10.1029/2019JA026890

Malara, F., Settino, A., Perrone, D., Pezzi, O., Guzzi, G., & Valentini, F. (2022). Exact shearing flow magnetized hybrid kinetic equilibria with inhomogeneous temperature. *The Astrophysical Journal*, **941**(2), 201. https://doi.org/10.3847/1538-4357/aca077

Masson, A., & Nykyri, K. (2018). Kelvin-Helmholtz instability: Lessons learned and ways forward. *Space Science Reviews*, **214**(4), 71. https://doi.org/10.1007/s11214-018-0505-6

Matsumoto, Y., & Hoshino, M. (2004). Onset of turbulence induced by a Kelvin-Helmholtz vortex. *Geophysical Research Letters*, **31**(2), L02807. https://doi.org/10.1029/2003GL018195

Matsumoto, Y., & Hoshino, M. (2006). Turbulent mixing and transport of collisionless plasmas across a stratified velocity shear layer. *Journal of Geophysical Research (Space Physics)*, **111**(A5), A05213. https://doi.org/10.1029/2004JA010988

Matsumoto, Y., & Seki, K. (2010). Formation of a broad plasma turbulent layer by forward and inverse energy cascades of the Kelvin-Helmholtz instability. *Journal of Geophysical Research (Space Physics)*, **115**(A10), A10231. https://doi.org/10.1029/2009JA014637

Mitchell, D. G., Kutchko, F., Williams, D. J., Eastman, T. E., Frank, L. A., & Russell, C. T. (1987). An extended study of the low-latitude boundary layer on the dawn and dusk flanks of the magnetosphere. *Journal of Geophysical Research*, **92**(A7), 7394–7404. https://doi.org/10.1029/JA092iA07p07394

Miura, A. (1984). Anomalous transport by magnetohydrodynamic Kelvin-Helmholtz instabilities in the solar wind-magnetosphere interaction. *Journal of Geophysical Research*, **89**(A2), 801–818. https://doi.org/10.1029/JA089iA02p00801

Nakamura, T. K. M. (2021). The earths low-latitude boundary layer. In R. Maggiolo, N. André, H. Hasegawa, & D. T. Welling (Eds.), *Magnetospheres in the solar system* (Vol. 2, pp. 177–191). https://doi.org/10.1002/9781119815624.ch12

Nakamura, T. K. M., Eriksson, S., Hasegawa, H., Zenitani, S., Li, W. Y., Genestreti, K. J., et al. (2017). Mass and energy transfer across the Earth’s magnetopause caused by vortex-induced reconnection. *Journal of Geophysical Research (Space Physics)*, **122**(11), 11505–11522. https://doi.org/10.1002/2017JA024346

Nakamura, T. K. M., Fujimoto, M., & Otto, A. (2008). Structure of an MHD-scale Kelvin-Helmholtz vortex: Two-dimensional two-fluid simulations including finite electron inertial effects. *Journal of Geophysical Research (Space Physics)*, **113**(A9), A09204. https://doi.org/10.1029/2007JA012803

Nykyri, K., & Otto, A. (2001). Plasma transport at the magnetospheric boundary due to reconnection in Kelvin-Helmholtz vortices. *Geophysical Research Letters*, **28**(18), 3565–3568. https://doi.org/10.1029/2001GL013239

Nykyri, K., Otto, A., Lavraud, B., Mouikis, C., Kistler, L. M., Balogh, A., & Rème, H. (2006). Cluster observations of reconnection due to the Kelvin-Helmholtz instability at the dawnside magnetospheric flank. *Annales Geophysicae*, **24**(10), 2619–2643. https://doi.org/10.5194/angeo-24-2619-2006

Plaschke, F., Kahr, N., Fischer, D., Nakamura, R., Baumjohann, W., Magnes, W., et al. (2016). Steepening of waves at the duskside magnetopause. *Geophysical Research Letters*, **43**(14), 7373–7380. https://doi.org/10.1002/2016GL070003

Pollock, C., Moore, T., Jacques, A., Burch, J., Gliese, U., Saito, Y., et al. (2016). Fast plasma investigation for Magnetospheric Multiscale. *Space Science Reviews*, **199**(1–4), 331–406. https://doi.org/10.1007/s11214-016-0245-4
Pu, Z. Y., Yei, M., & Liu, Z. X. (1990). Generation of vortex-induced tearing mode instability at the magnetopause. *Journal of Geophysical Research*, **95**(A7), 10559–10566. https://doi.org/10.1029/95JA095iA07p10559

Russell, C. T., Anderson, B. J., Baumjohann, W., Bromund, K. R., Dearborn, D., Fischer, D., et al. (2016). The magnetospheric Multiscale magnetometers. *Space Science Reviews*, **199**(1–4), 189–256. https://doi.org/10.1007/s11214-014-0057-3

Settino, A., Khotyaintsev, Y. V., Graham, D. B., Perrone, D., & Valentini, F. (2022). Characterizing satellite path through Kelvin-Helmholtz instability using a mixing parameter. *Journal of Geophysical Research (Space Physics)*, **127**(2), e2021JA029758. https://doi.org/10.1029/2021JA029758

Settino, A., Perrone, D., Khotyaintsev, Y. V., Graham, D. B., & Valentini, F. (2021). Kinetic features for the identification of kelvin-helmholtz vortices in in situ observations. *The Astrophysical Journal*, **912**(2), 154. https://doi.org/10.3847/1538-4357/abf1f5

Shi, Q. Q., Tian, A. M., Bai, S. C., Hasegawa, H., Degeling, A. W., Pu, Z. Y., et al. (2019). Dimensionality, coordinate System and reference frame for analysis of In-Situ space plasma and field data. *Space Science Reviews*, **215**(4), 35. https://doi.org/10.1007/s11214-019-0601-2

Song, P., & Russell, C. T. (1992). Model of the Formation of the low-latitude boundary layer for strongly Northward interplanetary magnetic field. *Journal of Geophysical Research*, **97**(A2), 1411–1420. https://doi.org/10.1029/97JA02377

Sorriso-Valvo, L., Catapano, F., Retinò, A., Le Contel, O., Perrone, D., Roberts, O. W., et al. (2019). Turbulence-Driven ion beams in the magnetospheric kelvin-helmholtz instability. *Physical Review Letters*, **122**(3), 035102. https://doi.org/10.1103/PhysRevLett.122.035102

Stawarz, J. E., Eriksson, S., Wilder, F. D., Ergun, R. E., Schwartz, S. J., Pouquet, A., et al. (2016). Observations of turbulence in a Kelvin-Helmholtz event on 8 September 2015 by the Magnetospheric Multiscale mission. *Journal of Geophysical Research (Space Physics)*, **121**(11), 11021–11034. https://doi.org/10.1002/2016JA023458

Swisdak, M. (2016). Quantifying gyrotropy in magnetic reconnection. *Geophysical Research Letters*, **43**(1), 43–49. https://doi.org/10.1002/2015GL066980

Tarvus, V., Turc, L., Zhou, H., Nakamura, T., Settino, A., Blasl, K., et al. (2024). Hybrid-Vlasov modelling of ion velocity distribution functions associated with the Kelvin–Helmholtz instability with a density and temperature asymmetry. *The Astrophysical Journal*, **974**(1), 62. https://doi.org/10.3847/1538-4357/ad697a

Terasawa, T., Fujimoto, M., Karimabadi, H., & Omidi, N. (1992). Anomalous ion mixing within a Kelvin-Helmholtz vortex in a collisionless plasma. *Physical Review Letters*, **68**(18), 2778–2781. https://doi.org/10.1103/PhysRevLett.68.2778

Terasawa, T., Fujimoto, M., Mukai, T., Shinohara, I., Saito, Y., Yamamoto, T., et al. (1997). Solar wind control of density and temperature in the near-Earth plasma sheet: WIND/GEOTAIL collaboration. *Geophysical Research Letters*, **24**(8), 935–938. https://doi.org/10.1029/96GL04018

Vernisse, Y., Lavraud, B., Eriksson, S., Gershman, D. J., Dorelli, J., Pollock, C., et al. (2016). Signatures of complex magnetic topologies from multiple reconnection sites induced by Kelvin-Helmholtz instability. *Journal of Geophysical Research (Space Physics)*, **121**(10), 9926–9939. https://doi.org/10.1002/2016JA023051

Wilder, F. D., Ergun, R. E., Schwartz, S. J., Newman, D. L., Eriksson, S., Stawarz, J. E., et al. (2016). Observations of large-amplitude, parallel, electrostatic waves associated with the Kelvin-Helmholtz instability by the magnetospheric multiscale mission. *Geophysical Review Letters*, **43**(17), 8859–8866. https://doi.org/10.1002/2016GL070404