The Kelvin–Helmholtz Instability at the Boundary of Relativistic Magnetized Jets

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 (Blumen et al. 1975; Blandford & Pringle 1976; Turland & Scheuer 1976; Ferrari et al. 1978, 1980; Pu & Kivelson 1983; Kivelson & Zu-Yin1984; Sharma & Chhajlani 1998; Bodo et al. 2004, 2013,2016, 2019; Osmanov et al. 2008; Prajapati & Chhajlani 2010; Hamlin & Newman 2013; Sobacchi & Lyubarsky 2018; Berlok & Pfrommer 2019; Pimentel & Lora-Clavijo 2019; Rowan 2019), depending on whether the relative motion is nonrelativistic or ultrarelativistic, whether the two fluids have comparable or different properties (respectively, “symmetric” or “asymmetric” configuration), whether the flow is incompressible or compressible, and whether or not the fluids are magnetized.

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 flow-aligned 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 pre-energized by reconnection (e.g., Sironi & Spitkovsky 2014; Zhang et al. 2021; Sironi 2022) can then experience shear-driven acceleration (Rieger 2019; Wang et al. 2021, 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 limb-brightened relativistic jets (e.g., in Cygnus A—Boccardi et al. 2016; and in 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-pressure-dominated 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.

We consider a planar vortex sheet interface in the x–z plane at y = 0, as shown in Figure 1. The jet (y > 0) is cold and magnetized, with field B _0j = (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 . We adopt Gaussian units such that c = 4 π = 1 and define all velocities in units of c.

Zoom In Zoom Out Reset image size

We describe the flow with the equations of relativistic magnetohydrodynamics (RMHDs; 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 factor (), 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 Alfvén speed , where

and the jet enthalpy density is w_0j ≈ ρ_0j 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.

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 Equation (1), such that the perturbed variables take the form φ ≈ φ₀ + φ₁, where φ₀ and φ₁ are the background and the first-order perturbed variables respectively. The perturbed electric field in the jet is E ₁ = − v ₁ × B _0j in the ideal MHD limit, ⁴ where v ₁ is the perturbed velocity in the jet frame.

In the jet, we consider perturbed variables φ₁ of the form φ₁ ∝ 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 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 direction 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 φ₁ of the form , where 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, . By Lorentz transformations of and , we obtain

Since l_j and l_w are Lorentz invariant, by solving Equations (8) and (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 equation, Equation (1), together with the first-order pressure balance equation,

and the displacement matching condition at the interface,

yielding

Using Equation (5), we can eliminate w_0w /w_0j from Equation (13) and, finally, the dispersion relation for the surface wave at the interface can be obtained by equating Equation (10) and the square of Equation (13):

By introducing the following notations,

Equation (14) can be rewritten as (Sobacchi & Lyubarsky 2018; Rowan 2019)

The dispersion relation in Equation (16) holds for arbitrary values of c_sw, v_A, v, , and , subject only to the assumptions of a cold jet and an unmagnetized wind.

Since Equation (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 Equation (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 and . The expressions for l_w and l_j can be obtained from the derivation of Equation (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 Equation (16) can be obtained. We first note that the parameters in Equation (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₀ + c₁ + c₂ ², where c₀, c₁, and c₂ are constant with respect to and terms higher than ² are ignored. Substituting this into Equation (16) and comparing coefficients of various powers of on both sides, we can find an approximate solution for all six roots of Equation (16). If we define an effective Mach number,

, and recognize that , 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 Equation (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 Figures 2 and 3, we compare the numerical solution (left column) of Equation (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 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our analytical approximation allows us to determine the range of M_e where the system is unstable. If λ in Equation (22) is imaginary, then also Λ_± has a nonzero imaginary part. We then find the values of M_e that satisfy λ² = 0 and obtain the following unstable bounds: for M_e < 1,

whereas for M_e > 1,

where

Note that the condition is equivalent to the obvious requirement v < 1. The condition in Equation (23) can be equivalently cast as , 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 Figure 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.

Equations (23) and (24) fully characterize the instability boundaries in Figures 2 and 3. In particular, the vertical white dotted lines in the figures illustrate the upper bound in Equation (24) for the special case , 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.

4.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

In Appendix A, we demonstrate that this analytical approximation for the special case M_e = 1 is in good agreement with the numerical solution.

Equation (27) allows us to identify the range of M_e (near unity) where the diverging growth rate in Equation (19) should rather be replaced by Equation (27). By equating the imaginary parts of in Equation (19) and in Equation (27), and solving for M_e , we can obtain the upper bound for Equation (19) such that for . We expect to be close to unity, so we assume M_e = 1 in μ and λ for Λ₊ of Equation (19). The resulting expression for can then be written as

where we require < 3^3/22^1/5 ξ^−1/2 for real .

4.2. 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 Figure 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 Figure 4, where we show the peak growth rate as a function of M and , where we define

Zoom In Zoom Out Reset image size

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 (the white line), which comes from the instability condition in Equation (23). The range of unstable Mach numbers extends up to M < 1/v_A (the vertical white line), which simply corresponds to the requirement v < 1.

When the unstable mode propagates perpendicular to the magnetic field (), 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 Equation (16), similar to the one of Equation (2) in Blandford & Pringle (1976), i.e.,

where is the total enthalpy of the jet, namely the sum of the gas enthalpy w_0j and the magnetic enthalpy:

Then the dispersion relation Equation (16) can be equivalently written as

which, by setting , is exactly the same as Equation (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 , the instability behaves similarly to the case of a hydrodynamic jet. Here, the magnetic field provides pressure, but not tension.

We have studied the linear stability properties of the KHI for relativistic, asymmetric, magnetized flows, with a 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 (, as long as v < 1). For ψ = 0, the velocity shear needs to overcome the magnetic tension, and only super-Alfvénic jets are unstable (). 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 nucleus (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, gamma-ray 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, 2016, 2019; Sobacchi & Lyubarsky 2018). Second, 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, and 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 Figures 5 and 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 for Figure 5 and for Figure 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 (for Figure 5) and (for Figure 6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the main body of the paper, we have looked for an analytical approximation of the form ϕ ≈ c₀ + c₁ + c₂ ², where c₀, c₁, and c₂ are constant with respect to and terms higher than ² 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 Σ_± ² is required for identifying the physical solutions that satisfy the Sommerfeld condition. The explicit expression for Σ_± is