Linear Stability Analysis of Relativistic Magnetized Jets: The Minimalist Approach

A minimalist approach to the linear stability problem in fluid dynamics is developed that ensures efficiency by utilizing only the essential elements required to find the eigenvalues for given boundary conditions. It is shown that the problem is equivalent to a single first-order ordinary differential equation, and that studying the argument of the unknown complex function in the eigenvalue space is sufficient to find the dispersion relation. The method is applied to a model for relativistic magnetized astrophysical jets.


Introduction
Understanding the stability properties of magnetized plasma flows is important for unraveling the basic characteristics of many phenomena we observe in astrophysics, solar, and space physics.Unstable perturbations may significantly alter the dynamics of these flows and be the reason behind changes in the shape; bulk velocity; magnetization; heating; particle acceleration; and consequently, the emitted radiation.
Linear stability analysis is a tool that enables simplifying the problem as much as possible, focuses on the main ingredients of each mechanism, and corroborates more complicated and expensive work, such as numerical simulations or laboratory experiments.It has its own difficulty, mainly due to the complex mathematics involved, especially in the relativistic regime, and using curvilinear coordinates.This remains true even if one neglects dissipative effects related to finite viscosity and resistivity, and assumes a simplified unperturbed state of a cylindrical flow with helical magnetic field, which approximates astrophysical jets in the propagation phase well.
All existing works on ideal fluids (neglecting dissipative effects) share a common procedure: they arrive in a system of two first-order differential equations in the complex domain, or equivalently, one second-order differential equation.The requirement that the solution of this system needs to satisfy certain boundary conditions gives the dispersion relation.In this paper, we show that this procedure may be simplified by reducing the number of equations by half.Following this novel approach, which we dub minimalist, one needs to solve only one differential equation in order to find the eigenvalues of the problem.In Section 2, we present the equations of the approach, and in Section 3, we explain how to solve the boundary condition, again following a minimal path by using only the argument of the eigenfunction.In Section 5, we apply the method to a model of a magnetized jet, extending a result of [6] to relativity.We conclude with a discussion in Section 6.

The Minimalist Approach
The linearization of the relativistic ideal magnetohydrodynamic equations in cylindrical coordinates (ϖ, ϕ, z), assuming a cylindrical unperturbed state and perturbations of the form Q 1 (ϖ) exp[i(mϕ + kz − ωt)] for all physical quantities, yields the following system of two complex first-order differential equations: Here y 1 is related to the Lagrangian displacement in the radial direction, and y 2 is related to the total pressure in the perturbed locations of fluid elements.There are known algebraic relations giving all the other quantities in terms of y 1 and y 2 , and thus, these two functions fully determine the solution.The pair of functions y 1 and y 2 is the most convenient choice because they are continuous everywhere, even at cylindrical surfaces, where the undisturbed state is discontinuous, e.g., in the interface between the jet and its environment.Details on the form of the unperturbed state, the derivation of the system of Equation ( 1), the expressions for the various F ij /D, and the proof that y 1 and y 2 should be continuous everywhere can be found in [16].
The above system is linear, and if our purpose is to find the eigenvalues and eigenfunctions, the proportionality constants are free.The only constraints that need to be satisfied are the boundary conditions.All kinds of boundary conditions are discussed in [16] and they refer to either on the axis, an interface where the unperturbed state is discontinuous, infinity, or a solid boundary (wall).In all these cases, the boundary conditions give the ratio and not the functions y 1 and y 2 separately.Thus, a minimalist approach is to work with this ratio, which has the advantage that it is uniquely defined for each eigenvalue, but most importantly, that the number of differential equations is reduced by half.The new equation is non-linear, but this does not burden the procedure if one satisfies the boundary conditions using a shooting method.

Differential Equation
The differential equation for the new unknown function Y can be derived by direct differentiation and using the system (1).It is the following single complex equation: Any solution of this equation that satisfies certain boundary conditions, which will be discussed below, and is continuous everywhere is an eigenfunction corresponding to a particular eigenvalue.The original pair of functions, if needed, can be found a posteriori using the equations (that are equivalent to the original system) (Note that only one variable, i.e., y 2 , needs to be found by solving the related differential equation, the other is given by y 1 = Yy 2 .)An instructive approximate solution of Equation (3), which shows its typical expected behavior, is presented in Appendix A.

Boundary Conditions on the Axis
On the axis, there are regularity conditions on y 1 and y 2 that are translated in a condition on their ratio Y.The detailed derivation of the former can be found in [16].In addition, the behavior of Y near the axis is discussed in Appendix B.1.
The resulting boundary condition at the symmetry axis is

Boundary Conditions at Infinity
If at large distances from the axis, the medium is static and homogeneous with zero B 0ϕ , as it is usually assumed, the solution for Y can be found analytically; see Appendix B.2.In any case, whatever the unperturbed state of the jet environment, care must be taken that the perturbations vanish as ϖ → ∞, and if they are oscillating, they correspond to outgoing waves in the radial direction.

Boundary Conditions at Interfaces
The function Y is everywhere continuous, including possible interfaces where the unperturbed state is discontinuous.
In the subcase of a solid boundary (wall), we simply require Y = 0 at the boundary.

Finding the Eigenvalues
The dispersion relation is found by solving Equation (3) subject to boundary conditions.In the following, we use the temporal approach (given a real k and integer m, search for complex eigenvalues ω), although the procedure is the same in the spatial approach (given a real ω and integer m, search for complex eigenvalues k).
Also, to be specific, suppose we have a jet with radius ϖ j and we know Y in its environment as a function of ϖ, k, m, and ω.Its value Y| ϖ=ϖ + j on the boundary is denoted as Y BC and is an analytic function of the complex variable ω.For given m, k, and all possible ω, we can integrate in the interior of the jet, starting from the axis using the appropriate boundary condition there, and when we reach the interface ϖ = ϖ j from the left, we find the value Y| ϖ=ϖ − j , which, for brevity in this section, we simply call Y.This is also an analytical function of the complex variable ω.The accepted eigenvalues ω are the ones for which Y = Y BC .
This procedure with obvious modifications can be applied to any other case.For example, if we have discontinuities in the jet interior, we simply cross them, keeping Y continuous and again work with Y = Y BC at the jet surface.If we do not know the solution in the environment, we simply continue the integration in the ϖ > ϖ j regime, and the boundary condition that will determine the dispersion relation is at infinity.Alternatively, we begin the integration from a large distance, reach ϖ j from the right, and match with the solution from the left.In any case, a complex equation of the form Y = Y BC at some surface will determine the eigenvalues.

Roots and Poles as Positive and Negative Line Charges
It is essential to work with the difference Y − Y BC in the complex plane ω.Since we are interested in unstable modes, we consider only the half plane ℑω > 0. The properties of this analytic function help us to find the eigenvalues, which are its roots.It is convenient to write it as Y − Y BC = e −Φ+iΨ , where −Φ and Ψ are the real and imaginary parts of ln(Y − Y BC ).
There is a direct analogy between the complex plane ω and a Cartesian xy plane by writing ω = x + iy.The Cauchy-Riemann conditions give ∇Φ⊥∇Ψ and that Φ and Ψ satisfy the Laplace equation at all points y > 0, except the positions where Φ becomes infinity, i.e., the roots and poles of the function Y − Y BC .
Suppose this function has roots ω rn , n = 1, 2, . . ., and poles ω pm , m = 1, 2, . . ., with the corresponding positions in the xy plane r rn and r pm .The function Φ(x, y) can be thought of as an electric potential associated with line charges, i.e., positive at the positions of the roots and negative at the positions of the poles (sources of the potential are also line charges located at positions with ℑω ≤ 0, but we are interested in finding the line charges only at ℑω > 0).The function Ψ is the stream function of the corresponding electric field −∇Φ = ∇Ψ × ẑ (the field lines are isocontours of Ψ and are normal to the isopotentials Φ = constant). 1 These can be proved by looking at the form of the potential/stream function near roots and poles.Writing the function as Y − Y BC = C rn (ω)(ω − ω rn ) q rn (with q rn as the multiplicity of the root), we indeed see that the dominant contributions near a root are Φ ≈ −q rn ln |r − r rn | + C rn and Ψ ≈ q rn arctan y − y rn x − x rn + D rn (where C rn and D rn are constants), corresponding to an electric field of a positive line charge −∇Φ = ∇Ψ × ẑ ≈ q rn r − r rn |r − r rn | 2 .Similarly, near a pole, by writing the function as Y − Y BC = C pm (ω) (ω − ω pm ) q pm , we find that the dominant contributions near the pole are Φ ≈ +q pm ln |r − r pm | + C pm and Ψ ≈ −q pm arctan y − y pm x − x pm + D pm (where C pm and D pm are constants), corresponding to an electric field of a negative line charge Near each positive line charge (a root), the relations Φ ≈ −q rn ln |r − r rn | + C rn and Ψ ≈ q rn arctan y − y rn x − x rn + D rn mean that the polar coordinates in a system with the axis at the line charge are e −Φ+C rn and Ψ − D rn .It is important to note, first, that Φ is +∞ at the location of the line charge and decreases as we move away, and second, that Ψ increases as we move counterclockwise around the line charge.Working similarly, near a negative line charge (a pole), we conclude that we have the opposite behavior: Φ is −∞ at the location of the line charge and increases as we move away, and Ψ decreases as we move counterclockwise around the line charge.Evidently, the potential and the stream function are not independent.We can find one from the other using the Cauchy-Riemann conditions and their isocontours are normal to each other.In the spirit of the minimalist approach, we can use only one.The stream function is a better choice since it does not involve infinities.
Concluding, and returning to the ω plane, in order to find the eigenvalues, we only need the isocontours of Ψ = Arg[Y − Y BC ].An example is shown in the upper panel of Figure 1.The points where Ψ experiences jumps equal to π are either roots or poles.Moving counterclockwise around such a point, if Ψ increases, it is a root, and if Ψ decreases, it is a pole.Moving from right to left, we see six line charges at locations ω ≈ 8.1 + 0.1i, ω ≈ 7.74 + 0.29i, ω ≈ 6.82 + 0.47i, ω ≈ 6.71 + 0.57i, ω ≈ 6.46 + 0.28i, and ω ≈ 6.35 + 0.25i.Checking the change in Ψ, we conclude that the first is a pole, the second a root, etc.
In fact, even the full map of isocontours of Ψ = Arg[Y − Y BC ] is not necessary.If we are in a position at the ω plane, find its Ψ, and then find the ∇Ψ at this point (by looking the Ψ of its neighbors), we know that a positive charge can be found in the direction opposite to the field ∇Ψ × ẑ.We can move from point to point on the curve Ψ = constant and we will reach a positive charge.In another variant that is probably even more economic, we can plot the isocontours of just two neighboring values −π < Ψ 1 < Ψ 2 ≤ π.The roots/poles are the endpoints of the contours, while the direction of ∇Ψ (it is normal to the contours pointing from the Ψ 1 to the Ψ 2 > Ψ 1 contour) is sufficient to understand which end is the root (it is in the direction opposite to the field ∇Ψ × ẑ).An example is shown in the upper panel of Figure 2.   The "Spectral Web" method was developed and presented in [2,21] for nonrelativistic magnetohydrodynamic flows.According to this method, in cases where y 1BC = 0, the eigenvalues can be found in the ω complex plane as intersections of the "solution paths", where {ℜ This system gives the roots but also "spurious roots" corresponding to y 2BC = 0.They also discuss the "complex oscillation theorem", according to which, along the solution path and in between the spurious roots, y 1 ϖy 2 , which they call the alternator, is real and monotonic function of the arc length, and similarly along the conjugate path, the alternator is purely imaginary and monotonic function of the arc length.
The analogy with the minimalist approach that uses the analytical function Y is obvious: the "solution paths these paths are shown in the middle panel of Figure 2), the roots correspond to positive line charges, and the spurious roots to negative line charges (poles).As we move along a field line Ψ = Arg[Y − Y BC ] = constant (any constant, not only 0, ±π/2, and π) approaching a positive line charge, the potential Φ = − ln |Y − Y BC | monotonically increases.At the position of the positive line charge, the potential becomes +∞ and the argument Ψ experiences a jump equal to π (since it corresponds to the angular polar coordinate in the local system with the line charge at the axis).As we move away from the positive line charge, the potential decreases, reaching −∞ if we meet a negative line charge.
As suggested by [2], this property helps with counting the eigenmodes that correspond to roots we meet as we move along a single field line.However, there is not always a single line connecting all the roots, making the counting of solutions impossible in general.In fact, in the particular case, there is a single line connecting all roots and poles, but it is not the "solution path" nor the "conjugate path".The "Spectral Web" corresponds to the field lines Arg[Y − Y BC ] = 0, π/2, as well as their continuations after jumps Arg[Y − Y BC ] = π, −π/2.Although these values do not seem to have any special significance relative to others, it is convenient to include these "paths" in the ω plane to show the location of line charges through the crossings of the paths (as already discussed, any other choice of isocontours will also show these locations).We include them in the following and continue to name the plot a "Spectral Web", although it has important additional information.The paths themselves are not enough; we need to know the values of Ψ and separate roots from poles knowing whether Ψ increases or decreases when moving counterclockwise around a line charge (or equivalently, understand the direction of ∇Ψ and the electric field ∇Ψ × ẑ).
The middle panel of Figure 1 shows the potential for illustrative purposes (since the eigenvalues were already found from the contours of Ψ alone).We verify that Φ = +∞ at the roots and Φ = −∞ at the poles.We also see that the isopotentials are normal to the field lines shown in the upper panel as contours of Ψ (to see the angles correctly, it is necessary to have the same scaling in the ℜω and ℑω axes).The bottom panels of Figure 1 show the eigenfunctions for the three eigenvalues found.All the results of that figure correspond to the case examined in Section 6 of [16].
A last point to discuss is related to the possibility to encounter infinities of Y when we integrate Equation (3).Since Y = y 1 /y 2 is defined as a ratio, one could expect infinities at points where y 2 = 0.However, working in the complex domain, an infinite value of Y requires both the real and imaginary parts of y 2 to vanish simultaneously with perfect accuracy, which is something that never happens during a numerical integration.Even if we encounter such a point, we automatically pass through it without a problem.An example is shown in Appendix C.

Energy Consideration
In the minimalist approach, where we do not solve an equation reminiscent of Newton's second law, one may think that the connection with the energy principle that is often used to discuss instabilities (through the increase/decrease of kinetic energy due to a decrease/increase of potential when moving from an equilibrium, in analogy to a ball on a hill/valley; see, e.g., [22]) has been lost.However, there is still some connection, as one could expect.
The energy of a particular mass of the plasma in general evolves in time according to the relation • da, which can be easily proved using the equations of motion 2 .Considering the volume of the jet and using the boundary condition on its perturbed surface, according to which the magnetic field remains always normal to the boundary, we find This has the simple meaning that the change in energy of the jet is due to the work done by the total pressure of its environment.On the perturbed boundary, we also have (see Using the functions y 1 and y 2 , we can express Since we need to keep quadratic (nonlinear) terms in Equation ( 7), we should carefully replace the real parts of the functions and take the real part of the product.Taking the mean value of the result, the zeroth-and first-order terms disappear. 3The mean value has a double meaning; over a length of the jet equal to multiples of the wavelength 2π/k, or over a time period equal to multiples of 2π/ℜω, provided that the growth time is much larger.The resulting expression for the mean value of the energy that is transferred from a length ∆z of the jet to its environment is evaluated at ϖ j .This is related to the "complementary energy" derived in [2] in connection with the force operator and its non-self-adjointness.The continuity of Y and y 2 ensures that the opposite energy per time is found when one integrates in the volume of the environment (because the area vector is opposite), and thus, the total energy remains constant.
Finding the function Y in the minimalist approach is sufficient to understand whether energy flows from the jet toward the environment, which is the case if ℜYℑω + ℑYℜω > 0, or the opposite.If one needs the exact value, first, y 2 needs to be found from Equation (4).Obviously, Equation ( 8) also shows that the system is unstable if there exist at least one mode with positive ℑω.
The result can be directly generalized to any cylindrical distance, and the mean power transferred from the interior to the exterior (algebraically) at a radius ϖ over a length ∆z is We can also integrate this equation and find the mean energy contained between the axis and the cylindrical distance ϖ over a length ∆z: where E 0 (ϖ) is the energy of the unperturbed state and we assume ⟨E (ϖ)⟩| t=0 = E 0 (ϖ).

An Example Case
In this section, the results of applying the minimalist approach to a model for relativistic magnetized jets are presented.In particular, we chose to explore how relativity changes the results of the constant pitch magnetic field model considered in [6].This is a simple model containing the basic characteristics of a jet and its environment.The goal is not only to discuss the physics involved-without, of course, being able to exhaust this interesting topic in this connection-but also to investigate the form of the eigenfunctions corresponding to different instability mechanisms using the new formalism.
Relativity is included in the dynamics in three ways: by allowing the bulk, Alfvén, and sound velocities to be relativistic in general.Defining the Alfvén velocity on the axis U Aa = B a √ ρ 0a and the corresponding Mach number M A = γ 0 V 0 U Aa , the dimensionless parameters that fully define the unperturbed state are U Aa , M A , ϖ 0 , and η.
We use units defined in [16], i.e., lengths in ϖ j , wavelengths in 1/ϖ j , velocities in c, frequencies in c/ϖ j , and the factor √ 4π absorbed into the magnetic field.We additionally set B a = 1, and thus, we measure densities, energy densities, and pressures in units of B 2 a .We assume in the following that U Aa = 1, M A = 1, ϖ 0 = 0.33, and η = 1; therefore, the bulk four velocity is and the Alfvén three velocity on the axis is Figure 3 shows the Spectral Web for k = 2 and m = −1, 0, 1 from top to bottom.The top panel of Figure 3 shows two eigenvalues, one at ω ≈ 1.44 + 0.14i (shown more clearly in the upper right panel) and a second at ω ≈ 0.84 + 0.31i.(It also shows two poles; we recall that Arg[Y − Y BC ] increases as we move counterclockwise around the eigenvalues, while it decreases as we move counterclockwise around poles.) The m = −1 has a sign opposite to the azimuthal magnetic field, making it possible to fulfill the resonant surface relations k co • B co = 0 and ω co = 0 (for which D = 0).This is related to the current-driven instability (CDI).Thus, the eigenvalue ω ≈ 1.44 + 0.14i corresponds to CDI since it is absent for m = 0 and m = +1.
The second eigenvalue corresponds to the Kelvin-Helmholtz instability (KHI).It is present for all m (at a slightly shifted position) and is the so-called surface mode (SM), ordinary mode, or fundamental mode.
Figure 4 shows the Spectral Web for k = 10 and m = −1, 0, 1 from top to bottom.We observe that the CDI is no longer present, which is something that is expected for sufficiently large k beyond the value that satisfies the resonant surface relation.
We also observe that new modes appear on the left of the SM (two modes in this case).These are connected to the Kelvin-Helmholtz instability and are called body modes (BMs) or reflection modes.
The three panels are very similar, which is something that is expected for sufficiently large k (compared with m).Repeating the process for all k, we find the dispersion relation shown in Figure 5 for the three values of m.Regarding the CDI, which is the main focus in [6], for the chosen parameters, relativity only slightly modify the results.The maximum growth rate is ℑω ≈ 0.24 at k ≈ 3.2, where ℜω = 2.28.In the frame comoving with the jet, these translate to ℜω co ≈ 0, ℑω co = γℑω ≈ 0.34, and k coz ≈ 2.26 − 0.24i.The ℑω co ϖ j v Aa ≈ 0.48 is compared with the maximum growth rate shown in Figure 1 of [6], i.e., ℑω co ϖ j v Aa ≈ 0.4 for k coz ≈ 2.3.
As expected, the CDI is advected with the flow (ℜω co ≈ 0) and the comoving wavevector along the motion is and ω co = 0 are shifted to k = −γ 0 m ϖ 0 ≈ 4.24, in agreement with the dispersion relation for CDI shown in the upper panel of Figure 5 in blue (for k larger than the one satisfying the resonant surface relation the mode is stable-this limiting k depended on m).
The eigenfunctions for k = 2 are shown in Figure 6.The first row corresponds to the CDI (for m = −1).The following three rows correspond to the KHI modes for the three values of m.
The corresponding y 1 and y 2 are shown in Figure 7.The absolute values of the amplitudes of the Lagrangian displacement |y 1 | ϖ and the total pressure |y 2 | are also shown.These plots definitely provide additional information for the perturbation compared with what is shown in Figure 6 through Y.For m = −1, the CDI is mostly concentrated near the axis (and displaces the axis, as expected), while the KHI is near the interface.The modes for m = 0 and m = 1 show mixed behavior.For m = 0, the displacement is maximum at the interface, but the perturbation of the pressure is maximum on the axis.For m = 1, we find displacement on the axis similarly to the CDI, but also on the interface.Figure 8 shows the eigenfunctions for k = 10.There are three eigenvalues for each value of m, which correspond, from larger to smaller ℜω, to SM, BM1, and BM2.For large k, all m give approximately the same result (eigenvalues and eigenfunctions).Counting seems to work well since more oscillations appear as we move to the next body mode.The general characteristics of the eigenfunctions, which are verified for even larger k as well, are as follows: a region near the axis that is controlled by the boundary condition; then a region with oscillations, the number of which depends on the number of the mode; and then the region near the interface controlled by the boundary condition.
The corresponding eigenfunctions y 1 and y 2 are shown in Figure 9 (the case m = −1 is shown, the others do not differ significantly due to the relatively large value of k).The numbers of oscillations are not as clear as in Y.They show, however, that the perturbations are more important near the interface.
-0.6 -0.Clearly, the Y on one hand and the y 1 and y 2 on the other give complementary information for the physics of the various unstable modes.However, the fact that the eigenvalues can be found in relation to Y alone means that this function needs to be more carefully analyzed in the various cases.Its relation to the number of oscillations was already shown.Another way to see it is through the trajectory that the eigenfunction Y follows in the complex Y plane as ϖ increases.For sufficiently large k, the number of windings equals the number of the mode.An example is shown in Figure 10.

Discussion
The minimalist approach presented in this paper offers a more economic way to find the eigenvalues of a linear stability problem by solving a single first-order differential equation for the complex function Y. Needless to say that although the presented examples are related to relativistic magnetohydrodynamic jets because this is perhaps the most challenging and unexplored area of stability, the minimalist approach can be applied in all kinds of fluid dynamic settings, simply by choosing the appropriate forms of the functions F ij /D and boundary conditions.For problems in cylindrical geometry, the formulas are given in [16] (these cover the nonrelativistic and the hydrodynamic limits).In other geometries, one can find these functions by linearizing the full equations.The method can also be applied to any other system that concerns the growth of perturbations in the linear regime.
A method to solve a complex boundary condition equation was developed by taking advantage of the fact that the function Y is an analytic function of the complex eigenvalue; using an analogy with electrostatics; and, in essence, generalizing the "Spectral Web" method presented in [2] by using only Arg[Y] and finding a way to distinguish roots and poles.It is interesting that the minimalist process of finding the eigenvalues depends only on the phase difference between y 1 and y 2 -which is the Arg[Y]-and not the phases themselves nor the energies, although a connection is discussed in Section 4.
The problem always leads to a complex equation Y − Y BC = 0, but there are many variations of this equation.For example, in the case presented in [16], the mathematical expression of the equation is different if we apply it to the wall interface ϖ = 1 (having integrated the equation from the axis to this interface), or in the interface ϖ = 0.5 (having integrated the equation from the axis to this interface and from the wall to this interface separately).Although the poles are different, the roots are always the same and equal to the eigenvalues of the problem.
To obtain the full picture of an instability, we of course also need y 1 or y 2 (in addition to their ratio Y).For example, for the energy consideration, we saw that by knowing Y, we find the direction of the energy flux but not the magnitude and its radial dependence.Also, the eigenfunction Y cannot say whether, e.g., there is a common dropping factor in both y 1 and y 2 since it is their ratio and the factor cancels out.We provide such an example in Appendix A. In addition, some connection with y 1 and y 2 is hidden inside the boundary conditions on the axis and at large distances.When we derived these conditions, we used information for y 1 and y 2 (to not diverge on the axis and to vanish at infinity, representing outgoing waves).However, this was done once, and since we know the boundary conditions, we are set to use the minimalist approach and the rest of the process uses only Y.This is advantageous because it more efficiently solves the most difficult part of the problem, which is to find the eigenvalues.Once we know an eigenvalue, it is trivial to return to the full system and calculate the perturbations of all quantities to obtain an overall picture of the eigenstate.
The fact that Y does not contain common factors of y 1 and y 2 also has its advantages.By dropping the common factors that affect the amplitudes, we concentrate on the phase of the oscillations, and this is why Y is apparently the function more closely related to the characterization of the body modes through the number of their oscillations.The oscillations are shown more clearly through Y, not only because the dropping factors are missing but also because the wavelength is smaller compared with the wavelength of y 1 and y 2 , as shown in the example in Appendix A.
Note also that the oscillatory behavior is expected in general for large k and for Kelvin-Helmholtz instability modes, but not necessarily in current-driven instability modes.There are interesting topics to be explored further, e.g., the connection between the solutions for Y and the function κ, and the phase difference between y 1 and y 2 in relation to the phase of each function separately, but we leave these for future studies.Although we already commented about the physics of the perturbations in various aspects, the goal of this paper was mainly to present the minimalist approach and the way to find the eigenvalues using it, which can be seen as the starting point for the rest.
We recall that the functions y 1 and y 2 are solutions of the linear problem, and thus, they can be freely multiplied with a complex constant (the same for both).Thus, we can freely multiply their amplitudes with a constant number and shift their phases with a constant angle.Their ratio, namely, the eigenfunction Y, is uniquely defined.In case we need to recover the units, y 1 has units of length squared, and thus, we should multiply it by ϖ 2 j , and y 2 has units of pressure, and thus, we should multiply it by B 2 a (which includes a factor of 4π, and thus, it is twice the magnetic pressure on the axis).
Funding: This research received no external funding.
Data Availability Statement: This research is analytical; no new data were generated or analyzed.If needed, more details on the study and the numerical results will be shared upon reasonable request to the author.

Conflicts of Interest:
The author declares no conflicts of interest.

Appendix A. Typical Behavior of Y
Noting the dominant dependences in F ij /D, as given by Equations ( 54)-( 57) in [16], and assuming , and A, B, C, and κ, we obtain the analytical solution Y = −ϖF 12 /D using the Bessel asymptotics, essentially corresponding to a harmonic oscillator with a complex frequency.For y 1 and y 2 , we obtain The function Y is proportional to cot( κϖ + ϕ 0 ) κ , which can be written as cot( κϖ We can see that the ℜ κϖ dependence creates oscillations in Y with wavelength π/ℜ κ and the ℑ κϖ dependence affects their amplitude.The imaginary part of cot( κϖ + ϕ 0 ) κ has a sign that is controlled by the sign of ℑ κ2 , and its absolute value peaks whenever 2 κϖ + 2ϕ 0 approaches even multiples of π.
The functions y 1 and y 2 have an additional common factor 1/ √ ϖ in their amplitudes and they oscillate with wavelength 2π/ℜ κ and some difference in their phases.
An example can be seen in Figure A1.This behavior turns out to be a very good approximation of exact solutions with | κ|ϖ larger than unity.Actually, the parameters for the solution in Figure A1 were chosen to fit an exact solution of the problem from [16], as shown in the bottom-right panel of Figure 1 (the bottom-left panel of Figure A1 and the bottom-right panel of Figure 1 are practically indistinguishable in region ϖ < 0.5).
Of course in the general case, κ is a function of ϖ and the wavelength of the oscillations is variable.In this case, the basic characteristics of the solution can be understood by replacing the argument of the tan with κ dϖ.
If | κ|ϖ is smaller than unity, the above approximations are not valid.This is the case at least in a small region near the axis, in which we know the behavior of the solution from the analysis of the boundary condition in Appendix B.1.).

Appendix B. Details on the
Actually, the general solution of ϖ |m| , and thus, the acceptable solution F = −1 corresponds to ϖ 0 = 0.
Concluding, for m ̸ = 0, the boundary condition at the symmetry axis is the one given in Equation ( 5) of the main text.
In more detail, for m = 0, the constant limits are b 11 = lim The boundary condition can be seen as a regularity to choose the acceptable solution Actually, we can find the general solution by noting that the term with ϖY is always negligible compared with the largest between the Y 2 /ϖ and ϖ, and thus, the differential equation can be approximated as b 21 ϖ . The acceptable solution corresponds to C = 0 and the unacceptable solution to C = ∞.
Concluding, for m = 0, the boundary condition at the symmetry axis is the one given in Equation ( 6) of the main text.

Appendix B.2. Boundary Conditions at Infinity
At large distances from the axis, assuming zero velocity and homogeneous medium with zero B 0ϕ , we have the case of Section 5.1 of [16].The solution is More generally, whatever is the unperturbed state of the jet environment, the perturbation should vanish at infinity, and if it is oscillating, it should correspond to waves propagating toward larger ϖ.

Appendix C. Integration through Infinities of Y
Infinities of Y are rare since they correspond to the vanishing of ℜy 2 and ℑy 2 simultaneously, but in general, it is possible at any distance ϖ for particular values of ω.They correspond to poles in the ω plane, and require fine tuning to find them, similar to the process we follow at the distance ϖ j for the function Y − Y BC .The needed perfect accuracy means that Y never becomes infinity and the numerical integration passes through such points without a problem.Figure A2 shows an example.It corresponds to the case of Figure 1.As seen in the top-left panel of that figure, by integrating from the axis, we find a pole at ϖ = 1 for ω ≈ 8.1 + 0.1i.To explore how the numerical integration behaved around the pole, we continued the integration for ϖ > 1.The resulting Y is shown in Figure A2.Even when we fine tuned the value of ω, approaching as much as possible to the value corresponding to the pole, the integration continued without a problem.The ln |Y| became large, but not infinity, and the Arg[Y] approached a step function, but its variation was smooth, as can be seen by zooming close to the point ϖ = 1.Comparing the result of the integration with the integration of the linear system (1)-or with the analytical expressions that exist for this particular model, see Ref. [16]-we found indistinguishable results (shown in the Figure A2 as dotted gray lines).
This example shows that even if the integration encounters a pole at some distance, it is capable of resolving the differential equation for Y around this point and gives the correct solution at larger distances.3) and the dotted gray lines (that are practically on top of the solid lines) to the integration of the linear system (1).

Notes 1
There are other ways to make connections with other physical We can think of the roots as line sources of incompressible fluid and the poles as line sinks.In another analogy, we can think of the roots/poles as line vortices with positive/negative circulation, respectively.Another possibility is to treat the real/imaginary parts of Y − Y BC as a potential/stream function.In this picture, there is an electric cylindrical dipole at the location of each pole and the potential/stream function vanishes at the positions of the roots.
The proof can be performed by writing the equations of motion as T µν ;ν = 0 and elaborating the energy momentum tensor, whose components in Cartesian coordinates are T 00 = γ 2 ξρ 0

Figure 1 . 2 πFigure 2 .
Figure 1.The two upper panels show the parts of the function Y − Y BC in the complex ω plane, corresponding to the case examined in[16].Y is the value at the jet radius ϖ = 1 as the result of the integration from the axis, and Y BC = 0 in this particular case.The bottom row shows the eigenfunctions for the three eigenvalues that satisfy Y − Y BC = 0.
It is the line in which (Y − Y BC )e iπ/4 is real, consisting of the parts Arg[Y − Y BC ] = −π/4, 3π/4 shown in the lower panel of Figure 2 with blue (the normal field lines corresponding to Arg[Y − Y BC ] = π/4, −3π/4 are also shown with orange lines).

2 πFigure 3 .
Figure 3.The Spectral Web.The parameters for each panel are shown at the top of the panel.The upper-right panel is a zoomed region of the upper-left panel.

Figure 4 .
Figure 4. Similarly to Figure 3 but for k = 10.Compared with that figure, here we see the new body modes emerging at the left.

Figure 6 .Figure 7 .
Figure 6.The eigenfunctions for k = 2.The parameters for each panel are shown at the top of the panel.

Figure 10 .
Figure 10.The trajectory in the Y plane for the BM2 for k = 10 and m = −1.

11
For large | κ|ϖ, we can further simplify the expressions F and obtain the solution Y = − F 12 κD

Figure A2 .
Figure A2.Integration through a pole at the distance ϖ = 1.The two parts of the function Y are shown: the ln|Y|, which becomes infinity at the pole, and the Arg[Y], whose π jump corresponds to the change in sign of y 2 .The solid (orange and blue) lines correspond to the integration of Equation (3) and the dotted gray lines (that are practically on top of the solid lines) to the integration of the linear system (1).

Boundary Conditions
•For m ̸ = 0, the boundary condition at the symmetry axis is Y(ϖ = 0) = 1 ∝ ϖ −|m| and y 2 ∝ ϖ −|m| .The boundary can be seen as a regularity condition to choose the acceptable solution Y = − It is acceptable if ℑλ ≥ 0 (such that the amplitude of y 2 does not diverge for ϖ → ∞) and ℜλ has the sign of ℜω, corresponding to outgoing waves.Note that asymptotically, the ratio of the Hankel functions in the above equations approaches −i; thus, for |λ|ϖ ≫ 1, we have j − E i E j + B i B j + P + E 2 + B 2 2 δ ij , with i, j = 1, 2, 3.The equation for the energy is ∂T 00 ∂t + ∇ • T 0i xi = 0. Its integral form in a volume whose boundary is moving with velocity V s (and thus, each part of the boundary creates a volume V s dt • da in the time interval dt) is d dt T 00 dτ + T j0 xj − T 00 V s • da = 0. Following the volume of a given mass, each point of the boundary moves with V s = V , and substituting the components of the tensor, we obtain dE dt+ E × B + PV − E 2 + B 2 2 V • da = 0. Substituting E = −V × B, we arrive at dE dt − [ΠV − (V • B)B] • da, with Π = P + B 2 − E 2 2 .For two complex functions A and B that are proportional to exp(iψ), the mean value of the product ⟨ℜA ℜB⟩ in an interval 3 * B].