The Kelvin-Helmholtz instability at the boundary of relativistic magnetized jets

We study the linear stability of a planar interface separating two fluids in relative motion, focusing on conditions appropriate for the boundaries of relativistic jets. The jet is magnetically dominated, whereas the ambient wind is gas-pressure dominated. We derive the most general form of the dispersion relation and provide an analytical approximation of its solution for an ambient sound speed much smaller than the jet Alfv\'en speed $v_A$, as appropriate for realistic systems. The stability properties are chiefly determined by the angle $\psi$ between the wavevector and the jet magnetic field. For $\psi=\pi/2$, magnetic tension plays no role, and our solution resembles the one of a gas-pressure dominated jet. Here, only sub-Alfv\'enic jets are unstable ($0<M_e\equiv(v/v_A)\cos\theta<1$, where $v$ is the shear velocity and $\theta$ the angle between the velocity and the wavevector). For $\psi=0$, the free energy in the velocity shear needs to overcome the magnetic tension, and only super-Alfv\'enic jets are unstable ($1<M_e<\sqrt{(1+\Gamma_w^2)/[1+(v_A/c)^2\Gamma_w^2]}$, with $\Gamma_w$ the wind adiabatic index). Our results have important implications for the propagation and emission of relativistic magnetized jets.


INTRODUCTION
The Kelvin-Helmholtz instability (KHI) (Von Helmholtz & Monats 1868; Lord Kelvin 1871)-at the interface of two fluids in relative motion -is one of the most ubiquitous and well-studied instabilities in the Universe. Since the pioneering works of Chandrasekhar (1961), the linear theory of the KHI has been investigated under a variety of conditions (Turland & Scheuer 1976;Blandford & Pringle 1976;Ferrari et al. 1980;Pu & Kivelson 1983;Kivelson & Zu-Yin 1984;Bodo et al. 2004;Osmanov et al. 2008;Blumen et al. 1975;Ferrari et al. 1978;Sharma & Chhajlani 1998;Prajapati & Chhajlani 2010;Sobacchi & Lyubarsky 2018;Berlok & Pfrommer 2019;Rowan 2019;Hamlin & Newman 2013;Bodo et al. 2013Bodo et al. , 2016Bodo et al. , 2019Pimentel & Lora-Clavijo 2019), depending on whether the relative motion is non-relativistic or ultra-relativistic, whether the two fluids have comparable or different properties (respectively, "symmetric" or "asymmetric" configuration), whether the flow is The boundaries of relativistic astrophysical jets may be prone to the KHI, given the relative (generally, ultrarelativistic) shear velocity between the jet and the ambient medium (hereafter, the "wind"). In jet boundaries with flowaligned magnetic fields, KH vortices can wrap up the field lines onto themselves, leading to particle acceleration via reconnection (Rowan 2019;Sironi et al. 2021). Particles preenergized by reconnection (e.g., Sironi & Spitkovsky 2014;Zhang et al. 2021;Sironi 2022) can then experience sheardriven acceleration (Rieger 2019;Wang et al. 2021Wang et al. , 2023 -i.e., particles scatter in between regions that move toward each other because of the velocity shear, akin to the Fermi process in converging flows (Fermi 1949). The KHI may then constitute a fundamental building block for our understanding of the origin of radio-emitting electrons in limbbrightened relativistic jets (e.g., in Cygnus A (Boccardi et al. 2016) and M87 (Walker et al. 2018)), and for the prospects of shear-driven acceleration at jet boundaries in generating Ultra High Energy Cosmic Rays.
A study of the KHI in this context needs to account for the unique properties of the boundaries of relativistic jets. First, the relative motion between the jet and the wind can be ultrarelativistic; second, while the wind is likely gas-pressure dominated, relativistic jets are believed to be magnetically dominated (Blandford & Znajek 1977), i.e., an asymmetric configuration. The linear stability properties of the KHI in this regime (of relativistic, asymmetric, magnetized flows) are still unexplored. In this Letter, we derive the most general form of the dispersion relation for the KHI at the interface between a magnetized relativistic jet and a gas-pressuredominated wind. We also provide an analytical approximation of its solution for wind sound speeds much smaller than the jet Alfvén speed, as appropriate for realistic astrophysical systems.

SETUP
We consider a planar vortex-sheet interface in the x-z plane at y = 0, as shown in Fig. 1. The jet (y > 0) is cold and magnetized, with field B 0 j = (B 0x , 0, B 0z ) lying in the x-z plane, and Alfvén speed v A . The ambient wind (y < 0) is gas-pressure supported (with sound speed c sw ) and has a vanishing magnetic field. We use the subscript "j" for the jet and "w" for the wind. We solve the system in the jet rest frame, where the wind moves with velocity v = vx. We adopt Gaussian units such that c = 4π = 1 and define all velocities in unit of c. Figure 1. A 3D schematic diagram of the boundary of the relativistic jet. The boundary (grey color) is located in the x − z plane. Above and below the boundary are the magnetically-dominated cold jet and the unmagnetized gas-pressure-supported ambient wind, respectively. q ∥ is the projection of the wavevector q onto the boundary. The jet is at rest and the wind has a relative shear speed of v. The magnetic field in the jet is B. θ is the angle between q ∥ and v while ψ is the angle between B and q ∥ .
We describe the flow with the equations of relativistic magnetohydrodynamics (RMHD) (e.g., Mignone et al. 2018;Rowan 2019): supplemented with the divergence constraints Here, ρ, ρ e , J, v, γ, B, E, w and p are the rest-mass density, charge density, current density, fluid velocity, Lorentz fac- , magnetic field, electric field, gas enthalpy density and pressure, respectively. For an ideal gas with adiabatic index Γ, the enthalpy can be written as w = ρ + Γp/(Γ − 1).
We assume a cold and magnetically-dominated jet, with (3) and the jet enthalpy density is w 0 j ≈ ρ 0 j for a cold jet. The wind has negligible magnetic field and is gas-pressure supported, with sound speed (Mignone et al. 2018) where w 0w is the wind enthalpy density. From pressure balance across the interface, where Γ w is the wind adiabatic index.

DISPERSION RELATION
The dispersion relation of surface waves at the interface can be found from the dispersion relations of body waves in both the jet and the ambient wind, together with the displacement matching at the interface. The dispersion relations of body waves in each of the two fluids can be found by linearizing Eqs. (1), such that the perturbed variables take the form φ ≈ φ 0 +φ 1 , where φ 0 and φ 1 are the background and the firstorder perturbed variables respectively. The perturbed electric field in the jet is E 1 = −v 1 × B 0 j in the ideal MHD limit 1 , where v 1 is the perturbed velocity in the jet frame.
In the jet, we consider perturbed variables φ 1 of the form φ 1 ∝ e i(q·x−ωt) where q = (k, l j , m) is the complex wavevector and ω is the complex angular frequency, both defined in the jet rest frame. Note that Im(ω) > 0 implies that the amplitude of the wave grows exponentially, i.e., an instability. We define the angle θ between the projection of the wavevector onto the x-z plane and the directionx of the shear flow velocity such that Similarly, we define the angle ψ between the wavevector projection onto the x−z plane and the jet magnetic field such that For a magnetized cold jet, the dispersion relation of its body waves describes magnetosonic waves in the cold plasma limit: In the wind, we consider perturbed variables φ 1 of the form φ 1 ∝ e i(q·x−ωt) whereq = (k, l w , m) is the complex wavevector andω is the complex angular frequency, both defined in the wind rest frame. For an unmagnetized wind, the dispersion relation of its body waves reduces to the one of sound waves, Since l j and l w are Lorentz invariant, by solving Eq. (8) and Eq. (9) for l j and l w respectively, we can construct a Lorentz invariant ratio: An independent way of obtaining l w /l j is to simultaneously solve the linearized RMHD equations, Eqs.
(1), together with the first order pressure balance equation and the displacement matching condition at the interface 1 Resistive effects are likely important in the non-linear stages ), but not for the linear analysis presented here. yielding Using Eq. (5), we can eliminate w 0w /w 0 j from Eq. (13) and finally, the dispersion relation for the surface wave at the interface can be obtained by equating Eq. (10) and the square of Eq. (13): By introducing the following notations, Eq. (14) can be rewritten as (Sobacchi & Lyubarsky 2018; Rowan 2019) The dispersion relation in Eq. (16) holds for arbitrary values of c sw , v A , v, cos θ and cos ψ, subject only to the assumptions of a cold jet and an unmagnetized wind.
Since Eq. (16) is a sextic equation in ϕ, it has a total of six (generally, complex) roots. However, not all of them may be acceptable. First, not all of the solutions will satisfy Eq. (13), since we have introduced spurious roots when squaring it. Also, by the Sommerfeld radiation condition (Sommerfeld 1912), only outgoing waves should be retained. This requires Im(l w ) < 0 and Im(l j ) > 0. The expressions for l w and l y can be obtained from the derivation of Eq. (13), so the Sommerfeld condition can be expressed as Since in general a sextic equation has no algebraic roots (Abel 1826), only approximate solutions of ϕ in Eq. (16) can be obtained. We first note that the parameters in Eq. (15) are chosen such that for a realistic wind with c sw ≪ v A , we have ϵ ≪ 1, whereas the other parameters do not depend on c sw . We then expand ϕ as a power series of ϵ of the form ϕ ≈ c 0 + c 1 ϵ + c 2 ϵ 2 , where c 0 , c 1 and c 2 are constant with respect to ϵ and terms higher than ϵ 2 are ignored. Substituting this into Eq. (16) and comparing coefficients of various powers of ϵ on both sides, we can find an approximate solution for all six roots of Eq. (16). If we define an effective Mach number µ ≡ cos 2 ψ − M 2 e , and recognize that γ −2 = 1 − M 2 v 2 A , then the approximate roots that correspond to the unstable modes can be written as where We find that the first order term Λ ± ϵ generally provides a good approximation of the numerical solution for ϕ. However, the second order term (which we write explicitly in Appendix B) is required for identifying the physical solutions that satisfy Eq. (13) and the Sommerfeld condition. At zeroth order in ϵ, the real part of the solution (i.e., the phase speed of unstable modes) is ϕ = M e , or equivalently ω/k = v, i.e., unstable modes are purely growing in the wind frame. In Fig. 2 and Fig. 3, we compare the numerical solution (left column) of Eq. (16) with our analytical approximation (right column). We fix c sw = 0.005 and consider v A = 0.2 and 0.8, so the assumption c sw /v A ≪ 1 of our analytical approximation is well satisfied. The analytical solution for Im(ϕ) displayed in the figures only employs the first order terms (as discussed above, we also use the second order terms to check the Sommerfeld constraint), yet it provides an excellent approximation of the numerical results, apart from M e = 1. For M e = 1, the first-order term Λ ± of our analytical approximation diverges. We discuss below this special case.
Our analytical approximation allows to determine the range of M e where the system is unstable. If λ in Eq. (22) is imaginary, then also Λ ± has nonzero imaginary part. We then find the values of M e that satisfy λ 2 = 0 and obtain the following unstable bounds: for M e < 1, whereas for M e > 1 where ν 1 = 2 cos 2 ψ + (1 + v 2 A )Γ 2 w (25) Figure 2. Dependence of the instability growth rate Im(ϕ) on θ and M e , for two choices of v A and two choices of cos ψ, as indicated in the plots. The left and right columns represent the numerical and analytical solutions, respectively. For cos ψ = 0, the maximum growth rate of the analytical solution is capped at its numerical counterpart to avoid the divergence at M e = 1. In all the panels, Im(ϕ) is then normalized to its maximum value, which is quoted in the panels themselves. The vertical dotted lines show the analytical upper bound on M e when cos ψ = 1, see Eq. (26). The vertical solid white lines indicate M e = 1.
Note that the condition M e < cos θ/v A is equivalent to the obvious requirement v < 1. The condition M e > cos ψ in Eq. (23) can be equivalently cast as v cos θ > v A cos ψ, which has a simple interpretation. The system is unstable if the projection of the shear velocity onto the direction of q ∥ (which we defined as the projection of the wavevector q on the x − z plane, see Fig. 1) is larger than the projection of the Alfvén speed onto the same direction. In other words, the shear is able to overcome magnetic tension. Eq. (23) and Eq. (24) fully characterize the instability boundaries in Fig. 2 and Fig. 3. In particular, the vertical white dotted lines in the figures illustrate the upper bound in Eq. (24) for the special case cos ψ = 1, which yields It follows that the unstable range at M e > 1 shrinks for v A → 1, but never disappears as long as v A < 1.

The special case M e =1
In the case M e = 1, our analytical approximation diverges. The singular case M e = 1 can be solved by expanding ϕ with a Puiseux series (Wall 2004;Wolfram Research 2020). Among the six approximate solutions of ϕ at M e =1, the only unstable one is where ξ = (cos 2 ψ − 1) 2 (cos 2 θ − v 2 A ) Γ 2 w cos 2 θ .
In Appendix A we demonstrate that this analytical approximation for the special case M e = 1 is in good agreement with the numerical solution.
Eq. (27) We expect M * e to be close to unity, so we assume M e = 1 in µ and λ for Λ + of Eq. (19). The resulting expression for M * e can then be written as where we require ϵ < 3 3/2 2 1/5 ξ −1/2 for real M * e .

Maximum growth rate
The results presented so far retain the explicit dependence on the angle θ between the projected wavevector q ∥ and the flow velocity v, and on the angle ψ between q ∥ and the magnetic field B (see Fig. 1). In practice, for a given Mach number M = v/v A and a fixed magnetic field orientation (e.g., with respect to the shear direction), one can determine the maximum growth rate, irrespective of the specific value of θ at which it is attained. This is presented in Fig. 4, where we show the peak growth rate as a function of M and cos Ω, where we define The plot shows that, for most magnetic field orientations, the peak growth rate is achieved at M ∼ 1. The exception is the case of fields nearly aligned with the shear velocity, where magnetic tension pushes the unstable region to higher M. The region of stability in the upper left corner is delimited by M = cos Ω (white line), which comes from the instability condition M e > cos ψ in Eq. (23). The range of unstable Mach numbers extends up to M < 1/v A (vertical white line), which simply corresponds to the requirement v < 1.

COMPARISON TO THE HYDRODYNAMIC CASE
When the unstable mode propagates perpendicularly to the magnetic field (cos ψ = 0), we expect magnetic tension to have no effect, and the solution should resemble the hydrodynamic asymmetric case discussed by Blandford & Pringle (1976). We demonstrate this by choosing a different parameterization in Eq. (16), similar to the one of Eq. (2) in Bland- where w * 0 j is the total enthalpy of the jet, namely the sum of the gas enthalpy w 0 j and the magnetic enthalpy: Then the dispersion relation Eq. (16) can be equivalently written as which, by setting cos ψ = 0, is exactly the same as Eq.
(1) in Blandford & Pringle (1976), where both the jet and the wind were assumed to be unmagnetized. We conclude that, even though our jet is magnetized, in the case cos ψ = 0 the instability behaves similarly to the case of a hydrodynamic jet.
Here, the magnetic field provides pressure, but not tension.

DISCUSSION AND CONCLUSIONS
We have studied the linear stability properties of the KHI for relativistic, asymmetric, magnetized flows, with focus on conditions appropriate for the interface between a magnetized relativistic jet and a gas-pressure-dominated wind. We derive the most general form of the dispersion relation and provide an analytical approximation of its solution for ϵ = c sw /v A ≪ 1. The stability properties are chiefly determined by the angle ψ between the jet magnetic field and the wavevector projection onto the jet/wind interface. For ψ = π/2, magnetic tension plays no role, and our solution resembles the one of a gas-pressure dominated jet. Here, only sub-Alfvénic jets are unstable (0 < M e ≡ (v/v A ) cos θ < 1, as long as v < 1). For ψ = 0, the velocity shear needs to overcome the magnetic tension, and only super-Alfvénic jets ). At zeroth order in ϵ, the phase speed of unstable modes is ω/k = v in the jet frame, i.e., they are purely growing in the wind frame.
Our analytical results are valuable for both theoretical and observational studies. They can be easily incorporated into global MHD simulations of jet launching and propagation, to identify KH-unstable surfaces (Chatterjee et al. 2020;Sironi et al. 2021;Wong et al. 2021). On the observational side, claims have been made that the KHI is observed along Active Galactic Nuclei (AGN) jets, based on the geometry of the outflow (Lobanov & Zensus 2001;Issaoun et al. 2022). Our formulae can place this claim on solid grounds, if estimates of the field strength and orientation and of the flow velocities are available. Besides AGNs, our results have implications for other jetted sources such as, but not limited to, gammaray bursts, tidal disruption events, X-ray binaries, and pulsar wind nebulae.
We conclude with a few caveats. First, the plane-parallel approach we employed is applicable only if the jet/wind interface is much narrower than the jet radius (for studies of surface instabilities in force-free cylindrical jets see, e.g., Bodo et al. 2013;Sobacchi & Lyubarsky 2018;Bodo et al. 2016Bodo et al. , 2019. Secondly, our local description implicitly assumes that the flow properties do not change as the KHI grows. Third, we have assumed the jet plasma to be cold, and the surrounding medium to be unmagnetized. These assumptions will be relaxed in a future work. We are grateful to R. Narayan for many inspiring discussions and collaboration on this topic. We would like to thank G. Bodo for many useful discussions and suggestions. L.S. acknowledges support from the Cottrell Scholars Award and the DoE Early Career Award DE-SC0023015. L.S and J.D. acknowledge support from NSF AST-2108201, NSF PHY-1903412 and NSF PHY-2206609. J.D. is supported by a Joint Columbia University/Flatiron Research Fellowship, research at the Flatiron Institute is supported by the Simons Foundation. For the singular case M e = 1, our analytical solutions take the form of the first order Puiseux series. Here we compare the analytical and numerical solutions. In Fig. 5 and Fig. 6, we plot the instability growth rate for M e = 1, comparing analytical and numerical solutions. We choose the same parameters as in the figures of the main paper, namely c sw = 0.005, and v A = 0.2 or 0.8. We fix cos ψ = 0 for Fig. 5 and cos θ = 1 for Fig. 6. We use solid and dashed lines to represent numerical and analytical solutions, respectively. The figures show that our analytical solutions in Puiseux series provide a good approximation to the numerical ones across the entire range of cos θ (for Fig. 5) and cos ψ (for Fig. 6).

B. THE SECOND ORDER TERMS
In the main body of the paper, we have looked for an analytical approximation of the form ϕ ≈ c 0 + c 1 ϵ + c 2 ϵ 2 , where c 0 , c 1 and c 2 are constant with respect to ϵ and terms higher than ϵ 2 are ignored. For the unstable solutions, we find that the first order term Λ ± ϵ generally provides a good approximation of the numerical solution. However, the second order term Σ ± ϵ 2 is required for identifying the physical solutions that satisfy the Sommerfeld condition. The explicit expression for Σ ± is Σ ± = M e µ[(1 − M 2 e )(cos 2 ψ(1 − v 2 A ) + M 2 e (1 + 3v 2 A ) − 2 − 2M 4 e v 2 A )Γ 2 w + 2(cos 2 ψ + M 2 e − 2)(µ 2 ± λ)] γ 2 (1 − M 2 e ) 2 Γ 2 w λ (B1) Figure 6. Dependence of the instability growth rate Im(ϕ) on cos ψ for two choices of v A and a fixed value of cos θ = 1 in the singular case M e = 1. Solid lines represent the numerical solutions while dashed lines represent the analytical solutions obtained by Puiseux series expansion in the main text.