Parametric decay instability of circularly polarized Alfvén wave in magnetically dominated plasma
Wataru Ishizaki Kunihito Ioka


APS/123-QED





 wataru.ishizaki@yukawa.kyoto-u.ac.jp
Astronomical Institute, Graduate School of Science, Tohoku University, Sendai 980-8578, Japan
Frontier Research Institute for Interdisciplinary Sciences, Tohoku University, Sendai 980-8578, Japan

 Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kitashirakawa Oiwake-cho, Sakyo, Kyoto 606-8502, Japan 

 
 kunihito.ioka@yukawa.kyoto-u.ac.jp

 Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kitashirakawa Oiwake-cho, Sakyo, Kyoto 606-8502, Japan 






















             








We investigate parametric decay instability (PDI) of circularly polarized Alfvén wave into daughter acoustic wave and backward Alfvén wave in magnetically-dominated plasma, in which the magnetization parameter σ (energy density ratio of background magnetic field to matter) exceeds unity.
We analyze relativistic magnetohydrodynamics (MHD), focusing on wave frequencies sufficiently lower than the plasma and cyclotron frequencies.
We derive analytical formulae for the dispersion relation and growth rate of the instability as a function of the magnetization σ, wave amplitude η, and plasma temperature θ.
We find that PDI persists even in high magnetization σ, albeit with a decreased growth rate up to σ→∞.
Our formulae are useful for estimating the decay of Alfvén wave into acoustic wave and heat in high magnetization σ plasma, which is a ubiquitous phenomenon such as in pulsars, magnetars, and fast radio bursts.











Parametric decay instability of circularly polarized Alfvén wave in magnetically dominated plasma
    Kunihito Ioka0000-0002-3517-1956
    April 29, 2024
=================================================================================================








§ INTRODUCTION
 

Fast radio bursts (FRBs) are the most luminous radio transients in the universe with millisecond duration <cit.>.
The identification of the host galaxies <cit.> and the excess of the dispersion measure over the free electron density in the Galaxy <cit.> have shown that almost all of the sources are at cosmological distances.
FRBs have very high brightness temperatures, suggesting that such emission is due to a coherent mechanism, i.e., a nonlinear process in which many particles emit photons of the same phase simultaneously, rather than a superposition of individual particle emissions <cit.>.
The mechanism of such high-intensity emission is still unresolved and remains an important problem in astrophysics <cit.>.

The mechanism and the source of FRBs are not known <cit.>.
The discovery of an FRB-like burst, FRB200428, coinciding with an X-ray burst from the Galactic magnetar SGR 1935+2154 has made highly magnetized neutron stars promising candidates for the source objects of FRBs <cit.>.
Two models are intensively debated regarding the emission region.
One is that the origin of FRBs is the maser emission in strongly magnetized relativistic shock waves generated outside the magnetosphere <cit.>.
It has been pointed out that this model has the difficulty that the model conflicts with the observation of the polarization <cit.> and the short time variability in some bursts <cit.>.
The other model is that the energy of the magnetic field of the neutron star is released and converted into emission inside the magnetosphere <cit.>.
Since the energy density of magnetic fields is strong near the surface of stars, magnetic reconnection <cit.> or starquake <cit.> at the surface of stars are discussed as possible energy sources.
Even though we consider phenomena occurring within the magnetosphere of a neutron star, the emission site remains uncertain.
If FRBs are generated in the magnetosphere near the stellar surface, the difficulty has been pointed out that their highly intense photons cannot escape from the magnetosphere <cit.>.
Therefore, the energy released at the stellar surface needs to be transported to an outer emission region and converted into radio waves <cit.>.



On the surface of the neutron star, the energy should be released in the form of magnetic pulse (including Aflvén wave) <cit.> or thermalized fireball <cit.>.
This paper specifically examines the cases where the release occurs in the form of Alfvén wave. 
Fluctuations in the magnetic field propagate through the magnetosphere as Alfvén waves (although fast magnetosonic waves are also conceivable, they can traverse the magnetic field and, thus, should disperse). 
In the neutron star magnetosphere, the energy density of the magnetic field significantly surpasses that of the plasma, so the Alfvén waves are expected to heat the plasma to relativistic energy as they dissipate, and then, to generate a bunch of charged particles, which might lead to coherent emission via some plasma instabilities <cit.>.
This expectation has motivated several FRB models using Alfvén waves.
Since the amplitude of Alfvén waves in the magnetosphere increases as they propagate outward, nonlinear effects during the propagation are important <cit.>. 
It remains a non-trivial question whether these large amplitude waves can traverse the magnetosphere without significant dissipation before reaching the FRB emission region. 


Finite amplitude Alfvén waves decay through the excitation of daughter acoustic waves and daughter Alfvén waves <cit.>.

This process is one of the parametric instability, particularly one of three-wave resonant interactions.
In the context of plasma physics, it is known as induced Brillouin scattering.
In this paper, we call this process parametric decay instability (PDI).
In the non-relativistic regime, this phenomenon has been extensively explored in the context of the solar coronal heating problem and the mechanism for launching the solar wind <cit.>.

When Alfvén waves propagate from the chromosphere into the corona, the PDI generates daughter acoustic waves, which immediately heat the coronal plasma through some dissipative process (e.g., shock heating) <cit.>.
In the magnetosphere of a neutron star, in contrast to solar physics, a magnetic field energy density is much larger than the matter energy density <cit.>. 
Although such systems are often analyzed using the force-free approximation <cit.>, this approach is inadequate for discussing Alfvén wave stability because this formulation lacks acoustic waves (i.e., slow waves). 
Additionally, within the force-free approximation, an Alfvén wave alone does not decay into lower-frequency Alfvén or fast waves through three-wave interaction (except for the interaction with ω=0 waves <cit.>
). 
Plasma heating due to the decay of the Alfvén wave is also considered important in the context of magnetar X-ray bursts (e.g., <cit.>, and in particular, <cit.> for calculation without force-free approximation[
Their calculation <cit.> is mainly focused on the macroscopic view of how much the plasma in the magnetosphere can be heated, rather than on the fundamental processes of the decay of the Alfven wave, which is our focus.
See Section <ref> for a comparison of their results with ours.
]).
Consequently, an investigation of the stability of Alfvén waves in magnetically-dominated plasmas, without imposing the force-free approximation, is essential.


The PDI in magnetically-dominated plasmas remains largely unexplored.  
Matsukiyo & Hada (2003) <cit.> investigated the stability of Alfvén waves in two-fluid pair plasmas, providing analytical expressions for the growth rate of the instability. 

However, as we will see later, the expression they obtained is inaccurate when the sound speed of plasma is sufficiently slow or fast (see equation (<ref>) and the following text for details).
Lopez et al. (2012) <cit.> incorporated full relativistic effects, but their work lacks an analytical discussion of the instability growth rate, impeding a systematic understanding. 
Deeply in the neutron star magnetosphere, the plasma frequency ω_p and cyclotron frequency ω_c significantly exceed the frequency of interest, even in comparison to the typical frequency of FRB emission <cit.>.
Thus, the magnetohydrodynamic (MHD) model, rather than the two-fluid model, is easy to handle for describing Alfvén wave propagation in the neutron star magnetosphere. 



In this paper, we study the stability of Alfvén waves in magnetically-dominated plasmas based on the MHD model, aiming to achieve a systematic understanding of the timescale of the instability in the neutron star magnetosphere.
We investigate the stability of Alfvén waves propagating parallel to the background magnetic field in the fluid's rest frame.
For simplicity, we analyze the 1-D case (i.e., both parent and daughter waves are propagating parallel to the background magnetic field) and focus on circularly polarized plane waves as parent Alfvén waves, which are easy to handle in analytical treatment due to the existence of exact analytical solutions.
We determine the growth rate of the PDI and elucidate its dependence on system parameters, particularly the magnetization parameter σ (the ratio of magnetic to matter energy density).
The structure of this paper is as follows:
In Section 2, we formulate the governing equations for perturbations added to the exact solution of Alfvén waves with finite amplitude.
Section 3 derives the dispersion relation from the perturbation equations, leading to the determination of the growth rate of the PDI.
Section 4 summarizes the obtained results, discusses the interpretation and application of our findings, and outlines future prospects.




§ FORMULATION





 §.§ Basic equations

In this paper, we analyze the stability of the Alfvén wave in the initially rest frame of plasma (we call it the reference frame).
We assume that the frequencies of waves considered in this paper are sufficiently smaller than both the cyclotron frequency and the plasma frequency, and that the plasma is quasi-neutral.
Therefore, we adopt the ideal magnetohydrodynamic (MHD) equation as the governing equation.
Furthermore, for simplicity, in this paper, we only consider the case where the waves propagate in a direction parallel to a constant magnetic field oriented in the z direction.
When there is no change of physical quantities in the x and y directions (i.e., we consider a 1D problem), the special relativistic MHD equations can be written as follows:


    1/c∂/∂ t[(ϵ+p) γ^2-p+1/8 π(E^2+B^2)]+∂/∂ z[(ϵ+p) γ^2 β_z+1/4 π(E×B)_z]=0,




    1/c∂/∂ t[(ϵ+p) γ^2 β_z+1/4 π(E×B)_z]+∂/∂ z[(ϵ+p) γ^2 β_z^2-1/4 π(E_z^2+B_z^2)]+ ∂/∂ z[p+E^2+B^2/8 π]=0   ,




    1/c∂/∂ t[(ϵ+p) γ^2 β_x y+1/4 π(E×B)_x y]+∂/∂ z[(ϵ+p) γ^2 β_xyβ_z-1/4 π(E_x y E_z+B_x y B_z)]=0,




    1/c∂B_x y/∂ t=[∇×(β×B)]_x y,


where ϵ is the energy density (including the rest mass energy density) of matter in the fluid rest frame, p is the pressure, β is the 3-velocity normalized by the speed of light c, γ=(1-β^2)^-1/2 is the Lorentz factor, E and B are the electric and magnetic fields in the reference frame, respectively.
The subscript z means the component in the z direction and xy means the projection onto the xy plane.
The electric field satisfies the ideal MHD condition:


    E=-β×B.


Note that the magnetic field in the z direction should be constant independent of t and z, as we can see from the z component of the induction equation and the solenoidal condition ∇·B=0.
Moreover, note that a stationary and uniform plasma (B=B_0e_z, ϵ= const, p= const, and β=0) is a solution to equations (<ref>)–(<ref>).
We also assume that the fluid is adiabatic.
As the thermodynamic equation relates between the energy density ϵ and the pressure p, we adopt the following:


    β_s^2≡(pϵ)_ ad,


where β_s is the sound speed normalised by c and the subscript ad means the adiabatic change.



 §.§ Circularly polarized finite amplitude Alfvén wave

There is an eigenmode in equations (<ref>)–(<ref>), which is a circularly polarized Alfvén wave.
A circularly polarised Alfvén wave is a solution to the MHD equation not only at small amplitudes but also at finite amplitudes.
The circularly polarized Alfvén waves propagating in the direction of the background magnetic field B_0 are expressed as


    B=B_0+δB,




    {δ B_x   =δ B cos (k z-ω t)
    δ B_y   =δ B sin (k z-ω t) 
    δ B_z   =0
    .,




    β=δβ=-χ_0β_A,ηδB/B_0,




    ϵ=ϵ_0,   p=p_0,


where β_A,η is the phase velocity of the finite amplitude Alfvén wave, δB and δβ are the magnetic field and velocity field fluctuations associated with the Alfvén wave respectively, and χ_0 is the direction of the wave propagation:


    χ_0= sgn(k).


Furthermore, a dimensionless number representing the amplitude of the wave is introduced as follows:


    δ B=η B_0.


Substituting these ansatz into equations (<ref>)–(<ref>), we obtain the following from the x component of the induction equation (<ref>):


    ω=χ_0β_A,η ck.


From the x component of the equation of momentum conservation (<ref>), we obtain


    [(ϵ_0+p_0) δγ^2+B_0^2/4 π] χ_0β_A,ηω-B_0^2/4 π c k=0,


where δγ=(1-δβ^2)^-1/2=(1-η^2 β_A,η^2)^-1/2.
Note that the same equations can be obtained from the y component of the equations (<ref>) and (<ref>).
By eliminating ω and k from equations (11) and (12), we obtain the phase velocity β_A,η as follows:


    β_A,η^2=β_A^2[1+β_A^2 η^2/2(1+√(1-(2 β_A^2 η/1+β_A^2 η^2)^2))]^-1,


where β_A is a phase velocity of the linear (i.e. infinitesimal amplitude) Alfvén wave normalized by the light speed, which can be written as 


    β_A^2 = B_0^2/4π(ϵ_0+p_0)+B_0^2.


From equation (<ref>), it can be seen that a circularly polarised Alfvén wave with finite amplitude expressed in the form of equations (<ref>)–(<ref>) is a wave propagating with a phase velocity expressed in equation (<ref>).
The phase speed of the wave decreases as η increases.
For η≪1, the phase velocity can be written as:


    β_A,η = β_A[1-1/2β_A^2(1-β_A^2) η^2+𝒪(η^4)].


As can be seen from this, the phase speed reduction appears on the order of 𝒪(η^2).





 §.§ linearization for daughter waves






As seen in the previous section, a circularly polarised Alfvén wave of finite amplitude is an exact solution of the MHD equation.
In this section, we investigate the stability of Alfvén waves by adding perturbations b_⊥, β_⊥, ϵ_ and β_ of magnitude 𝒪(ε), expressed in Fig. <ref> or in equation (<ref>) below, to the solution of the MHD equations (<ref>)–(<ref>):


    {[    B =B_0  +δB +b_⊥     ;    β    =   δβ +β_⊥  +β_;    ϵ =ϵ_0            +ϵ_ ].,


where b_⊥ and β_⊥ are the magnetic field and velocity field fluctuations of the component perpendicular to the z direction (i.e., b_⊥⊥B_0 and β_⊥⊥B_0), respectively, meaning Alfvén wave-like perturbations.
On the other hand, ϵ_ and β_ are fluctuations of the internal energy density and velocity parallel to the z direction (i.e., β_B_0), respectively, meaning acoustic wave-like perturbations.

Substitute equations (<ref>) into equations (<ref>)–(<ref>) and equate the coefficients of each order of ε with 0.
At the 0-th order of ε, the non-trivial relations is obtained from equations (<ref>) and (<ref>) as follows:


    [δγ^2(ϵ_0+p_0)+B_0^2/4π]
    	1/c∂/∂ t( δβ)
    	-B_0/4 πz(δB)
    	=0,




    1/ct(δB)-z(B_0δβ)=0.


These equations constitute the wave equation of the parent wave of magnitude 𝒪(η).
Note that equations (<ref>) and (<ref>) are trivial at the 0-th order.
The 1-st order equations for ε are as follows (see Appendix <ref> for derivation):


    1+β_s^2 δβ^2/1-δβ^21/c∂ϵ_/∂ t+(δγ^2(ϵ_0+p_0)) ∂β_/∂ z=
    
        -1/c∂/∂ t[(2 δγ^4(ϵ_0+p_0)+B_0^2/4 π)(δβ·β_⊥)+χ_0 δβδ B B_0 β_/4 π]
        +B_0/4 π[β_⊥·∂/∂ z(δB)+δβ·∂b_⊥/∂ z],





    [δγ^2(ϵ_0+p_0)+δ B^2/4 π] 1/c∂β_/∂ t+β_s^2 ∂ϵ_/∂ z=-∂/∂ z(δB·b_⊥/4 π)+B_0/4 π(δB·1/c∂β⊥/∂ t+b_⊥·1/c∂/∂ t(δβ)),





    [δγ^2(ϵ_0+p_0)+B_0^2+δ B^2/4 π] 1/c∂β_⊥/∂ t-B_0/4 π∂b_⊥/∂ z=-δγ^2(ϵ_0+p_0) β_∂/∂ z(δβ) 
    
    -δγ^2(1+β_s^2) ϵ_1/c∂/∂ t(δβ)-δββ_s^2 1/c∂ϵ_/∂ t+B_0/4 π1/c∂/∂ t(β_δB) 
    
    -2 δγ^4(ϵ_0+p_0)(δβ·β_⊥) 1/c∂/∂ t(δβ)+1/4 π∂/∂ z[ẑ·{(δβ×b_⊥)+(β_⊥×δB)}](δβ×B_0) 
      -B_0/4 πδβ[β_⊥·∂/∂ z(δB)+δβ·∂b_⊥/∂ z-B_0 1/c∂/∂ t(δβ·β_⊥)-χ_0 δβδ B 1/c∂β_/∂ t] 
    
    -δB·b_⊥/4 π1/c∂/∂ t(δβ)-(δB/4 π·1/c∂b_⊥/∂ t) δβ+χ_0δ B δβ/4 π1/c∂b_⊥/∂ t+δB/4 π(1/c∂β_⊥/∂ t·δB)+δB·β_⊥/4 π1/c∂/∂ t(δB),




    1/c∂b_⊥/∂ t-B_0 ∂β_⊥/∂ z=-∂/∂ z(β_δB).






 §.§ linearlization for parent wave


Hereafter, for simplicity, we assume a situation where the amplitude of the parent wave is not very large (i.e., η≲ 1), and assume that terms of magnitude 𝒪(εη^2) are ignored.
In this case, from equation (<ref>), β_A,η∼β_A for 1-st order equations.
Furthermore, ignoring the term 𝒪(εη^2) in equations (<ref>)-(<ref>), we obtain the following equations:


    1/c∂ϵ_/∂ t+(ϵ_0+p_0) ∂β_/∂ z=-1/c∂/∂ t[(2(ϵ_0+p_0)+B_0^2/4 π)(δβ·β_⊥)]+B_0/4 π[β_⊥·∂/∂ z(δB)+δβ·∂b_⊥/∂ z],




    (ϵ_0+p_0) 1/c∂β_/∂ t+β_s^2 ∂ϵ_/∂ z=-∂/∂ z(δB·b_⊥/4 π)+B_0/4 π(δB·1/c∂β_⊥/∂ t+b_⊥·1/c∂/∂ t(δβ)),




    ℰ1/c∂β_⊥/∂ t-B_0/4 π∂b_⊥/∂ z=-(ϵ_0+p_0) β_∂/∂ z(δβ)-(1+β_s^2) B_0 ϵ_/4 πℰ∂/∂ z(δB)-δββ_s^2 1/c∂ϵ_/∂ t+B_0/4 π1/c∂/∂ t(β_δB),




    1/c∂b_⊥/∂ t-B_0 ∂β_⊥/∂ z=-∂/∂ z(β_δB),


where 


    ℰ≡ϵ_0+p_0+B_0^2/4 π.


Note that all the terms in the right-hand sides of these equations are 𝒪(εη), and this implies that ε≪η is assumed because 𝒪(ε^2) terms are neglected.
On the right-hand side of these equations, there are some terms that did not appear in the analysis of parametric decay instabilities based on the non-relativistic MHD equations <cit.>.

These terms stem from (1) the velocity squared (i.e., 𝒪(β^2)), (2) the internal energy density that can be comparable to the rest mass energy density, and (3) the displacement current, which were neglected for the non-relativistic regime.




 §.§ dimensionless equations


To clarify the dominated terms in the equation, we normalize the equation and obtain a non-dimensionalized equation.
The quantity representing the amplitude of the parent wave is normalized as follows:


    δu≡δβ/β_A,   δe=δB/B_0.


Note that this non-dimensionalisation yields |δu|^2=|δe|^2=η^2.
Also, since 𝒪(η^2) is ignored here, the normalization scale of δβ is β_A rather than β_A,η.
The amplitude of the daughter Alfvén wave of order 𝒪(ε) is normalized as follows:


    u_⊥≡β_⊥/β_A,   e_⊥=b_⊥/B_0.


For the daughter acoustic wave of order 𝒪(ε):


    u_≡β_/β_s,    e_≡ϵ_/ϵ_0+p_0.




Here, we define two dimensionless parameters that characterize this system.
One is the magnetization parameter σ, defined as follows:


    σ≡B_0^2/4 π(ϵ_0+p_0).


By using σ, the Alfvén speed can be written as β_A^2=σ/(1+σ).
The other is the ratio of the sonic velocity to the Alfvén velocity:


    θ≡β_s/β_A.


Note that this corresponds to approximately the 1/2 power of a quantity called "plasma beta" in the non-relativistic case [
The plasma beta (β_ pl) is typically defined as the ratio of the plasma pressure to the magnetic pressure (i.e., β_ pl≡ p_0/(B_0^2/8π)).
The value derived from this definition and the square of the ratio of the sound speed to the Alfvén velocity is nearly equal, in fact for a polytropic gas with a specific heat ratio Γ (i.e., satisfying Γ-law), θ^2=(2/Γ)β_ pl.
].

Using equations (<ref>)–(<ref>), we obtain following normalized equations:


    1/c∂ e_/∂ t+β_s ∂ u_/∂ z=-β_A^2 1/c∂/∂ t(δu·u_⊥)-σδu·(1/c∂u_⊥/∂ t-β_A ∂e_⊥/∂ z),




    1/c∂ u_/∂ t+β_s ∂ e_/∂ z=-θ^-1β_A ∂/∂ z(δe·e_⊥)+σθ^-1δe·(1/c∂u_⊥/∂ t-β_A ∂e_⊥/∂ z),




    1/c∂u_⊥/∂ t-β_A ∂e_⊥/∂ z=θβ_A^2 1/c∂/∂ t(u_δe)-1/1+σ[β_s u_∂/∂ z(δu)+β_A e_∂/∂ z(δe)+β_s^2 1/c∂/∂ t(e_δu)],




    1/c∂e_⊥/∂ t-β_A ∂u_⊥/∂ z=-θβ_A ∂/∂ z(u_δe).


Here we used the following normalized equation of the order 𝒪(ε^0):


    1/c∂/∂ t(δu)
    	=
    	β_A∂/∂ z(δe),




    1/c∂/∂ t(δe)
    	=
    	β_A∂/∂ z(δu).


The second term on the right side of equation (<ref>) is found to be a small quantity of 𝒪(εη^2) using equation (<ref>).
Therefore, the equation (<ref>) becomes


    1/c∂ e_/∂ t+β_s ∂ u_/∂ z=-σ/1+σ1/c∂/∂ t(δu·u_⊥).


Similarly for equation (<ref>), we have


    1/c∂ u_/∂ t+β_s ∂ e_/∂ z=-θ^-1β_A ∂/∂ z(δe·e_⊥).






§ DISPERSION RELATION

As seen in the previous section, the 𝒪(ε) perturbation quantities u_⊥, e_⊥, u_, e_ follow the linear differential equations (<ref>), (<ref>), (<ref>), and (<ref>).
In this section, we examine the stability of the wave solutions in these equations by finding the dispersion relation.



 §.§ mode expansion

In order to find the dispersion relation, each physical quantity of the wave is divided into Fourier components.
Since the parent Alfvén wave is a circularly polarized wave with a given wavenumber k_0 and frequency ω_0, it can be written as


    δe=1/√(2)(δ e_0 exp(i ϕ_0) e_ R+),




    δu=-δe,


where ϕ_0=k_0z-ω_0t is the phase of the parent Alfvén wave, e_ R=(e_x+i e_y) / √(2) is the unit polarization vector for right-handed polarization, e_j (j=x,y,z) is the unit vector, and c.c. denotes the complex conjugate.
Note that |δ e_0|^2=η^2.
Here, although equation (<ref>) seems to represent only right-handed polarization waves, equation (<ref>) simultaneously represents not only right-handed polarized waves with positive ω_0 and positive k_0 but also negative ω_0 and negative k_0 at the same time <cit.>.
Furthermore, in the MHD regime (i.e., where the wave frequency is much smaller than the cyclotron frequency and plasma frequency), right- and left-handed polarizations are symmetric, so that what happens with one polarization also occurs with the other polarization.
Therefore, it is sufficient to discuss wave stability for the parent wave expressed by the equation (<ref>).


For the acoustic wave perturbation, we write the frequency as ω and the wavenumber as k as follows:


    e_=1/2(e_k exp (i ϕ)+),




    u_=1/2(u_k exp (i ϕ)+),


where ϕ=kz-ω t is the phase of the acoustic wave.
Since there are two daughter Alfvén waves produced by the interaction of the acoustic wave and the parent Alfvén wave, corresponding to Stokes and anti-Stokes waves, they are written as follows:


    e_⊥=1/√(2)(e_+exp(i ϕ_+)e_ R+)+1/√(2)(e_-exp(i ϕ_-)e_ R+),




    u_⊥=1/√(2)(u_+exp(i ϕ_+)e_ R+)+1/√(2)(u_-exp(i ϕ_-)e_ R+),


where ϕ_±=ϕ_0±ϕ=k_± z-ω_± t is the phase of the daughter Alfvén waves (- for the Stokes wave, and + for the anti-Stokes wave), and k_±=k_0± k and ω_±=ω_0±ω are wavenumbers and frequencies of daughter Alfvén waves, respectively.


Substituting equations (<ref>)–(<ref>) into equation (<ref>), the following equation for the Fourier component is obtained from the term proportional to exp(iϕ):


    ω e_k-cβ_s k u_k=β_A^2 ω(u_+δ e_0^*+u_-^* δ e_0).


From equation (<ref>), the following equation is obtained from the exp(iϕ) component:


    ω u_k-cβ_s k e_k=θ^-1 cβ_A k(e_+δ e_0^*+e_-^* δ e_0).


Similarly from the exp(iϕ_+)e_ R component of equation (<ref>), we obtain the following:


    ω_+ u_++cβ_A k_+ e_+=1/2δ e_0[(θβ_A^2 ω_+-1/1+σθ cβ_A k_0) u_k
    +1/1+σ(cβ_A k_0+θ^2 β_A^2 ω_+) e_k],


and from the exp(-iϕ_-)e^*_ R component:


    ω_- u_-+cβ_A k_- e_-=1/2δ e_0^*[(θβ_A^2 ω_--1/1+σθ cβ_A k_0) u_k
    +1/1+σ(cβ_A k_0+θ^2 β_A^2 ω_-) e_k].


Similarly, from the exp(iϕ_+)e_ R component of equation (<ref>), we obtain


    ω_+ e_++ cβ_A k_+ u_+=1/2θ cβ_A k_+ u_k δ e_0,


and from the exp(-iϕ_-)e^*_ R component, 


    ω_- e_-^*+ cβ_A k_- u_-^*=1/2θ cβ_A k_- u_k δ e_0^*.


Equations (<ref>)–(<ref>) are six equations for six Fourier components e_k, u_k, e_+, u_+, e_-^*, u_-^*.



 §.§ dispersion relation

Hereafter, we adopt the following unit system:


    ω_0=1,   k_0=1.


From the conditions for the existence of nontrivial solutions to e_k, u_k, e_+, u_+, e_-^*, u_-^* in equations (<ref>)–(<ref>), the dispersion relation is obtained.
By equating the determinant of the coefficient matrix of equations (<ref>)–(<ref>) with 0, the dispersion relation, which is the equation relating ω and k, is obtained as follows:












    (ω-k) (ω^2-θ^2 k^2) {(ω+k)^2-4}   =   η^2/(1+σ)^3(S_0 + S_1 σ + S_2 σ^2 + S_3 σ^3)
      
       +   η^4 k ωσ/(1+σ)^4(T_0 + T_1 σ + T_2 σ^2),


where


    S_0   =   k^2(ω^3+kω^2-3ω+k),
    
      S_1   =   2S_0-ω^3+kω^2-θ^2 k^2
      [2kω^2+(2k^2-3)ω-k],
    
      S_2   =   S_0-ω^3+kω^2-θ^2
      [ω^5+2kω^4+(k^2-3)ω^3-3k(k^2+1)ω^2-k^2(3k^2-7)ω+k^3],
    
      S_3   =   θ^2 k ω (ω-k) [(ω+k)^2-4],
    
      T_0   =   k,
    
      T_1   =   -ω + 2k - θ^2 (kω^2+k^2ω-k),
    
      T_2   =   -ω + k - θ^2 (kω-1)(ω+k).


Note that since we have only calculated correctly up to the lowest order of η in this paper, expressions of T_j (𝒪(η^2) times higher order terms than S_j (j=0,1,2,3)) may change when the calculation correctly incorporates higher order terms of η.
Note also that although we used the dispersion relation ω_0=cβ_A k_0 for the parent Alfvén wave propagating in the positive z axis in deriving equation (<ref>), the above expression is identical for a wave propagating in the opposite direction (i.e. the negative z direction).
Equation (<ref>) is an algebraic equation for ω and k, and its solution represents the wave that can propagate through the plasma in the presence of parent Alfvén wave.





For η=0, equation (<ref>) recovers the dispersion relation for waves propagating parallel to the magnetic field in a uniform plasma.
Figure <ref> shows the solution of equation (<ref>) for η=0 and θ=0.2.
Since the parent Alfvén wave is propagating in the z-axis positive direction, ω-k=0 can be written as ω_±=cβ_A k_± (in the physical unit), which represents the daughter Alfvén waves (Stokes and anti-Stokes waves) propagating in the same direction as the parent wave.
From ω+k-2=0, this can be written in terms of the frequency and wavenumber of the daughter Alfvén wave as ω_-=-cβ_Ak_-, indicating that it is the daughter Alfvén wave (Stokes wave) propagating in the opposite direction of the parent wave.
Similarly, ω+k+2=0 becomes ω_+=-cβ_Ak_+, which represents the daughter Alfvén wave (anti-Stokes wave) propagating in the opposite direction to the parent wave.
Furthermore, ω^2-θ^2k^2=0 represents acoustic waves propagating backward or forward in the z direction, as shown by the fact that ω^2=β_s^2k^2 in the physical unit.
















Figure <ref> shows the dispersion relation obtained by solving the equation (<ref>) under the parameters σ=10, η=0.2, and θ=0.2.
The curves shown in blue in the figure represent solutions with stable waves, while the curves shown in red dotted lines represent solutions that have non-zero imaginary parts (i.e., solutions with unstable growth).

In the figure, the instability occurs at the intersection point of the dispersion relations for the forward acoustic wave and the backward daughter Alfvén wave (Stokes waves).
Figure <ref> shows the real (upper panel) and imaginary (lower panel) parts of the frequency when the dispersion relation (<ref>) is solved with the same parameters as in Figure <ref>.
As can be seen from these, even for plasma with relativistic magnetization (in this case σ=10), the three-wave interaction which satisfies the resonance conditions ω_0=ω+ω_- and k_0=k+k_- is unstable.



 §.§ Non-relativistic limit

In order to check the validity of the obtained dispersion relation, we take the non-relativistic limit and compare it with the known results.
The non-relativistic limit here means the limit where the velocity scale of the system is much slower than the speed of light and the rest mass energy density is much larger than the internal energy density.
This can be achieved by taking the limit of Alfvén velocity β^2_A≪1, or since σ=β_A^2/(1-β_A^2), the limit of σ→0.
By taking the limit in equation (<ref>), we obtain


    (ω-k)(ω^2-θ^2 k^2){(ω+k)^2-4}=η^2 k^2(ω^3+k ω^2-3 ω+k).


This is the same dispersion relation as previous results obtained for non-relativistic MHD equation system <cit.>.
Therefore, in this regime (σ≪1), based on this dispersion relation, the daughter Alfvén and the acoustic waves grow exponentially.



 §.§ σ→∞ limit

As another extreme case, let us consider the case where σ→∞.
The dispersion relation (<ref>) leaving only the highest order of σ is as follows:


    (ω-k){(ω+k)^2-4}(ω^2-η^2 θ^2 k ω-θ^2 k^2)=0.


Since there are only real solutions of ω to this equation for any k, θ, and η, it can be seen that instability does not occur in the limit of σ→∞.



 §.§ 1/σ expansion

As seen in the previous two sections, in the limit of small σ (i.e., non-relativistic limit), the PDI occurs, but in the limit of σ→∞, no instability occurs.
This fact implies that the parametric decay instability (PDI) is suppressed as σ increases.
In order to investigate this quantitatively, we expand the dispersion relation with σ^-1 as a small quantity and investigate how the instability is suppressed up to the first order of σ^-1.
Hereafter, the term η^4 in the dispersion relation is ignored as small, and θ≪1 is assumed for simplicity.
The dispersion relation (<ref>) expanded to the first order of σ^-1 is as follows:


    (ω-k)(ω^2-θ^2 k^2)[(ω+k)^2-4]
    =
    η^2[S_3+1/σ(S_2-3 S_3)].


Since this equation is a 6-th-order algebraic equation, it cannot be solved analytically in general.
Therefore, in order to derive the conditions for the existence of unstable solutions, we set η to be small and expand the dispersion relation around the following resonance frequency ω_1 and wavenumber k_1, which satisfies both the dispersion relation ω=θ k for forward acoustic waves and ω+k-2=0 for backward Alfvén waves (Stokes waves).
The resonance frequency ω_1 and wavenumber k_1 can be written as


    ω_1=2θ/1+θ,   k_1=2/1+θ.


We take the frequency ω and the wavenumber k as follows:


    ω=ω_1+δω,  k=k_1+δ k.


We assume that δω and δ k are small quantities of magnitude about 𝒪(η).
This will be confirmed later to be correct in the case where the resonance condition is satisfied.
Expanding the left-hand side of equation (<ref>) to the lowest order of η, we obtain


    (ω-k)(ω^2-θ^2 k^2){(ω+k)^2-4}∼ -32 θδω(δω+δ k).


Note that here we used θ≪1 and keep only the lowest order of θ, but assuming that θ is greater than 𝒪(η).
Similarly, the following equation can be obtained for the right-hand side by expanding to the lowest order of η:


    η^2[S_3+1/σ(S_2-3 S_3)]
    ∼8η^2/σ.


Equating the above results, we thus obtain:


    δω^2+δ kδω + η^2/4σθ=0.


If equation (<ref>) has a solution with an imaginary part, equation (<ref>) has an unstable growing solution.
Regarding equation (<ref>) as a quadratic equation of δω, the condition that its discriminant is negative gives the following condition that δω has an imaginary part, namely:


    -ησ^-1/2θ^-1/2<δ k<ησ^-1/2θ^-1/2.


The imaginary part Γ≡ Im(δω) takes its maximum value Γ_ max when the resonance condition is satisfied (i.e. δ k=0) and is as follows:


    Γ_ max=1/2ησ^-1/2θ^-1/2.


Note that the validity of the assumption that δω and δ k are 𝒪(η) can be confirmed from equations (<ref>) and (<ref>).
As can be seen from equations (<ref>) and (<ref>), the instability ceases in the limit of σ→∞.



 §.§ Decay rate for general σ

By using the same method used in section <ref> to equation (<ref>), we can obtain the growth rate and the conditions under which instability occurs for general σ.
In this section, we continue to use 𝒪(η) as a small quantity and ignore the term of T_j on the right-hand side of the equation (<ref>).
Substituting equation (<ref>) into equation (<ref>) and deriving the lowest order equation for η, we obtain:


    δω^2
    +(1-θ)δ kδω
    -θδ k^2
    +η^2(1-θ)/4θ(1+θ)^2[1+σ(1+θ^2)]^2/(1+σ)^3
    =0.


The condition for δω to have a solution with an imaginary part is as follows:


    -ηθ^-1/2√(1-θ)/(1+θ)^21+σ(1+θ^2)/(1+σ)^3/2
    <
    δ k
    <
    ηθ^-1/2√(1-θ)/(1+θ)^21+σ(1+θ^2)/(1+σ)^3/2.


The growth rate is at maximum Γ_ max when δ k=0, and it can be written as:


    Γ_ max=1/2ηθ^-1/2√(1-θ)/1+θ1+σ(1+θ^2)/(1+σ)^3/2.


In the non-relativistic limit (i.e. σ≪1) and the low-beta limit (i.e. θ≪1), this expression becomes Γ_ max=ηθ^-1/2/2 and reproduces the well-known expression of the growth rate of the PDI for non-relativistic plasma.

This formula (<ref>) for the growth rate agrees with Matsukiyo & Hada (2003) <cit.> within the range of 𝒪(θ) [Although Matsukiyo & Hada (2003) assumed σ≲1 when deriving the growth rate, the σ-dependence of their formula is consistent with our result for general σ.
They justified neglecting some terms when expanding the dispersion relation around the resonance point by assuming that σ is not too large.
However, considering only points well close to the resonance point is sufficient to neglect these terms, making the assumption about σ unnecessary.
].
The discrepancy beyond 𝒪(θ) arises because their formulation implicitly assumes that the internal energy is significantly smaller than the rest mass energy.
Furthermore, the above result is also approximately derived from the growth rate of the induced Brillouin scattering for unmagnetized pair plasma, equation (A3) in Iwamoto et al. (2023) <cit.>, by multiplying (ω/ω_ c)^2, where ω_ c is the electron cyclotron frequency.
This is because the oscillation amplitude of electrons perpendicular to B_0 is suppressed by ω/ω_ c in the presence of a background magnetic field compared to the scattering of electromagnetic waves without a magnetic field.

In the bottom panel of Figure <ref>, the maximum growth rate obtained from equation (<ref>) is shown by a horizontal thick black dotted line.
The upper and lower limits of the range of wavenumber k in which instability occurs, represented by the equation (<ref>), are also indicated by the vertical dashed-dotted lines in the upper and lower panels.
As is clear from the figure, the obtained analytical expressions reproduce the numerical results well.
Figure <ref> shows the dispersion relation when σ=10^2 and η and θ are the same as in figure <ref>.
As σ increases, not only the maximum growth rate becomes smaller, but also the range of instability becomes narrower.
Takamoto et al. (2014) <cit.> performed simulations of Alfvén wave propagation based on the relativistic MHD equation and obtained a qualitative trend that Alfvén waves tend to be more stable as σ increases.
Our result, equation (<ref>), is consistent with their results and gives a quantitative interpretation of their nature.







Figure <ref> shows the dependence of the maximum growth rate Γ_ max on θ when σ=10^3 and η=0.2.
The red curve is the result of numerically and directly solving the dispersion relation (<ref>), and the blue dotted curve is the growth rate calculated from equation (<ref>).
It can be seen that the analytical formula (<ref>) is in good agreement with the numerical solution in θ≳10^-2.
However, for θ≲10^-2, equation (<ref>) fails to reproduce the numerical result.
This is because the approximation of |δω|≪ω_1 is broken.
The correct growth rate in such a small limit of θ will be obtained in the next section.



 §.§ Strong coupling limit: decay rate for small θ

As seen in Figure <ref>, when θ is small, equation (<ref>) fails to reproduce the numerical solution.
This is due to the fact that the condition |δω|≪ω_1, which was assumed in obtaining equation (<ref>), is violated.
In the context of plasma physics, such a situation is called strong coupling limit <cit.>.
The typical θ=θ_c where the approximation begins to break down is estimated from ω_1∼|δω| using equations (<ref>) and (<ref>) as follows:


    θ_c ∼ 2^-4 / 3η^2 / 3(1+σ)^-1 / 3.


Since Γ_ max∼Γ_ max(θ=θ_c) at θ≪θ_c (see Figure <ref>), we can presume that the growth rate at θ≪θ_c is proportional to η^2/3 from equation (<ref>).
Therefore, we assume that δω and δ k are 𝒪(η^2/3) and expand the dispersion relation again.
Note that θ in this regime is smaller than 𝒪(η^2/3) because θ_c=𝒪(η^2/3).
Expanding the equation (<ref>), we obtain the following equation from the lowest order 𝒪(η^2) terms:


    δω^2(δω+δ k)+η^2/1+σ=0.


The condition that the solution δω of equation (<ref>) has an imaginary part is that the wavenumber k is greater than the following value k_ min:


    k>k_ min=2-√(27/4)η^2/3(1+σ)^-1/3.


When δ k=0, the growth rate takes the maximum value Γ_ max as follows:


    Γ_ max=√(3)/2η^2 / 3(1+σ)^-1 / 3.



The black dotted line in Figure <ref> represents the growth rate calculated by equation (<ref>).
In the parameters of Figure <ref>, θ_c∼0.01, it can be seen that equation (<ref>) reproduces well the growth rate obtained by the numerical calculation for θ≪θ_c.
Note that, in the non-relativistic limit (i.e. σ→0), it reproduces the growth rate of PDI in non-relativistic plasmas at θ∼0 <cit.>.



§ DISCUSSION


The obtained expression for the growth rate of the instability reveals that the maximum growth rate Γ_ max is a decreasing function of σ, and in the limit of σ→∞, the instability ceases.
Since σ is the ratio of the energy densities of the matter and electromagnetic fields, the limit of large σ corresponds to a situation close to the force-free limit.
Our findings are consistent with the fact that an Alfvén wave is stable to resonant three-wave interactions in the force-free limit.
While it is not evident whether the methods using force-free equations and the MHD equations employed in our study are mathematically equivalent, we can demonstrate that they are equivalent within the scope of this study.
By eliminating u_⊥ from equations (<ref>) and (<ref>), and taking the limit of σ→∞, we obtain


    1/c^2∂^2 b_⊥/∂ t^2-∂^2 b_⊥/∂ z^2=0.


This is the equation satisfied by the magnetic field of the free electromagnetic wave.
Similarly, by eliminating e_⊥ from equations (<ref>) and (<ref>), and taking the limit of σ→∞, we obtain


    1/c^2∂^2/∂ t^2(u_⊥-θ u_δe)
    -∂^2/∂ z^2(u_⊥-θ u_δe)
    =0.


The quantity u_⊥-θ u_δe that satisfies the wave equation in vacuum can be interpreted as the electric field E_⊥,xy associated with the daughter Alfvén wave perpendicular to the background magnetic field B_0.
By multiplying both sides of equation (<ref>) by β_AB_0 as an outer product from the left side, we obtain


    1/c^2∂^2 E_⊥,xy/∂ t^2-∂^2 E_⊥,xy/∂ z^2=0,


where 


    E_⊥,xy=
    -(β_⊥×B_0+β_×δB).


This corresponds to the electric field in the xy-plane, except for the one associated with the parent Alfvén wave (see equation (<ref>)).
In this way, equations (<ref>) and (<ref>) describing daughter Alfvén waves reduce to the wave equation of free electromagnetic waves in the limit σ→∞.
On the other hand, the wave equation for acoustic waves still has a non-zero source term in the limit of σ→∞:


    1/c^2∂^2 ϵ_/∂ t^2-β_s^2∂^2 ϵ_/∂ z^2
    =
    -w_0[
    1/c^2∂^2/∂ t^2(δβ·β_⊥)
    -∂^2/∂ z^2(δB·b_⊥/B_0^2)
    ],




    1/c^2∂^2 β_/∂ t^2-β_s^2∂^2 β_/∂ z^2
    =
    -1/c∂^2/∂ t ∂ z(
    δB·b_⊥/B_0^2
    -β_s^2δβ·β_⊥).


Although the acoustic waves are generated by the interaction between parent and daughter Alfvén waves even for σ→∞, the PDI is considered to be suppressed because the acoustic wave is not possible to excite daughter Alfvén waves.

In this work, we have treated the amplitude η of the parent Alfvén wave as constant and calculated the growth rate of infinitesimal amplitude daughter waves.
Consequently, the decay rate of the parent wave cannot be determined from this analysis.
A formulation that considers the case ϵ∼η is necessary for a more precise derivation.
However, according to the general consideration of three-wave interactions, there is an implication that the decay rate of the parent wave differs only by a logarithmic factor from the growth rate of the daughter wave, namely the decay time scale of the parent wave can be written as


    t_ decay∼Γ_ max^-1ln(η/ϵ_ d,0),


where ϵ_ d,0 is the normalized initial amplitude of the daughter waves (e.g., for the Alfvén wave, ϵ_ d,0=b_⊥/B_0 at the initial state) <cit.>.
Therefore, the growth rate obtained in this study is effective for estimating the decay rate of the parent wave.


We can apply our results to relativistically magnetized plasmas, such as the plasma of the magnetosphere of a neutron star.
Let us estimate the timescale for the decay of Alfvén waves generated on the surface of a neutron star, assuming they serve as the energy source for the X-ray flare of SGR1935+2154 and FRB 20200428A <cit.>.
The estimated amplitude of the parent wave, based on the luminosity L_ X∼10^41  erg s^-1 of the X-ray burst (i.e., L_ X∼4π R_ X^2 c η^2 B_ NS^2/8π where R_ X∼10^4  cm is the size of the energy release region[
The energy release might occur in one part of the stellar surface rather than the entire surface.
Here we adopt a radius of R_ X∼10^4  cm which is estimated by using the luminosity L_ X∼10^41  erg s^-1 and effective temperature T_ eff∼80  keV of the X-ray burst <cit.>.
] <cit.>, and B_ NS∼ 10^15  G is the surface magnetic field), is η∼3×10^-4 <cit.>.
Although the magnetization parameter σ in the magnetosphere is not precisely known and depends on the efficiency of particle production in the magnetosphere, an estimation of σ∼ 10^4 has been given by previous studies <cit.> [
If the e^± fireball is formed in the X-ray burst, σ should be much smaller <cit.>
].

In the magnetosphere, where there are a large number of photons from the X-ray burst, we can expect the temperature of the magnetosphere to be T∼𝒪( keV) due to the Compton heating <cit.>, so that we evaluate the sound speed of magnetospheric plasma as β_ s∼(1  keV/511  keV)^1/2∼0.05.
According to equation (<ref>), since θ_c∼8×10^-5, we find that θ≫θ_c and the growth rate is appropriate to be estimated by equation (<ref>).
This suggests that the Alfvén waves propagate approximately 2× 10^4 wavelengths before decaying.
Assuming a wave period comparable to the rotation period of the central neutron star, ω_0∼1 Hz, λ∼10^10 cm, it implies that the waves do not decay significantly within the light cylinder.
However, if we consider waves excited at a frequency similar to the observed FRB frequency (around 10^9 Hz), the waves would decay at a distance of approximately 7×10^5 cm, indicating limited escape from the magnetosphere.
While these estimates are rough approximations, our results may be crucial for discussing the radiation mechanisms of FRBs.

In this study, we analyze the idealized setting of an infinitely extended uniform plane wave in a dissipationless system.
However, when applying these results to real systems, it is essential to evaluate the impacts of the following effects.
The first is the finite bandwidth effect of the parent wave.
Although stability analysis is conducted for monochromatic waves in this study, it is more natural to consider waves excited in general systems as wave packets with a finite bandwidth Δω.
In the context of the parametric instability, it is known that when the bandwidth of the parent wave is larger than the growth rate of daughter waves in the monochromatic case (Δω≫Γ_ max), the effective growth rate of the daughter waves becomes Γ_ eff∼Γ_ max^2/Δω, suppressed by Γ_ max/Δω <cit.>.
On the other hand, if Δω≪Γ_ max, the phase of the parent wave and the daughter waves can be considered synchronized during the time scale of daughter wave growth, allowing us to use the analysis of the monochromatic parent wave.
The second is the effect of dissipation in the plasma.
The dominant dissipation process in the plasma can vary depending on the system under consideration.
Here, we assume an appropriate timescale for the dissipation process as Γ_ diss^-1.
For instability to occur in the wave, the timescale of instability without considering dissipation should be shorter than the dissipation timescale (i.e., Γ_ max>Γ_ diss) <cit.>.
In this case, the growth rate of instability can be estimated as Γ_ eff∼Γ_ max-Γ_ diss.
When the effect of dissipation is weak, so that the dissipation proceeds over a much longer timescale compared to the growth timescale of the instability, it is justified to use the growth rate obtained from the analysis of the monochromatic parent wave.

In this study, we have only investigated the PDI within the linear stage (i.e., 𝒪(ε)≪1), so that further research is necessary to understand the development in the nonlinear stage.
In systems where σ≫ 1, meaning the electromagnetic energy density greatly exceeds the matter energy density, non-thermal high-energy phenomena such as non-thermal particle acceleration may occur due to the transfer of electromagnetic energy to low-density regions.
Additionally, for the generated acoustic waves, damping effects including the Landau damping are expected to occur, thus requiring an analysis incorporating kinetic effects to further quantitatively understand the PDI.
Magnetohydrodynamic simulation encounters difficulty in dealing with large σ and cannot analyze kinetic phenomena such as particle acceleration, making Particle-in-Cell (PIC) simulations more suitable for such studies.
While this study also applied linearization to η, further extensions are possible.
Equations (19)-(22) are formulated for general η, and we are currently conducting ongoing research to directly analyze these equations.
When η reaches unity, higher-order wave-wave interactions (e.g., 4-wave interaction) may become significant, although this phenomenon lies beyond the scope of this study.
Exploring higher-order wave-wave interactions through simulations, particularly PIC simulations, will provide valuable insights into understanding the stability of Alfvén waves.



§ SUMMARY

We have investigated the parametric decay instability (PDI) of circularly polarized Alfvén waves in a plasma with relativistic magnetization, i.e., σ≫1, based on the special relativistic magnetohydrodynamic (MHD) equations.
Our results indicate that even in plasmas with relativistic magnetization, the Alfvén wave is unstable to three-wave interactions where the parent wave generates two daughter waves - a forward propagating acoustic wave and a backward propagating Alfvén wave.
We have obtained the analytical expression of the maximum growth rate (see equation (<ref>)) as follows:


    Γ_ max/ω_0=1/2ηθ^-1/2√(1-θ)/1+θ1+σ(1+θ^2)/(1+σ)^3/2,


where η=δ B/B_0 is the amplitude of the parent Alfvén wave, σ=B_0^2/(4π(ϵ_0+p_0)) is the magnetization, θ=β_ A/β_ s is the ratio of the sound speed to the Alfvén velocity, and ω_0 is the frequency of the parent Alfvén wave.
This is only valid when the ratio θ of the sound speed to Alfvén velocity is greater than θ_c, which is expressed by equation (<ref>).
For θ≪θ_c, the growth rate can be written as follows:


    Γ_ max/ω_0=√(3)/2η^2 / 3(1+σ)^-1 / 3.


As can be seen from these equations, we have observed that the instability is suppressed with increasing magnetization σ.
Our results for the instability growth rates are consistent with the known non-relativistic limits.
Such a process may provide a clue to the energy conversion mechanism from electromagnetic field to plasma in environments around compact objects, where strongly magnetized plasmas are expected to exist, e.g., in the magnetospheres of neutron stars or black holes.


We are grateful to 
Toshio Terasawa, Kenji Toma, Takayoshi Sano, Yohei Kawazura, Shuichi Matsukiyo, Masanori Iwamoto, Shoma Kamijima, Rei Nshiura, Koutarou Kyutoku, Tomoki Wada, Shota Kisaka, Kazumi Kashiyama, and Shuta Tanaka
for fruitful discussion and valuable comments.
This work is supported by Grants-in-Aid for Scientific
Research No. 23H01172, 23H05430, 23H04900, 22H00130, 20H01901, 20H01904, 20H00158 (KI), 21J01450 (WI) from the Ministry
of Education, Culture, Sports, Science and Technology (MEXT) of Japan.









§ DERIVATION OF LINEARLIZED EQUATIONS




 §.§ Fluid equations with electromagnetic field as external field

In order to derive equations (<ref>)-(<ref>), we first simplify equations (<ref>)-(<ref>).
Using the Maxwell equations, the energy-momentum conservation law can be rewritten into the fluid equations in the external electromagnetic field as follows:


    1/c∂/∂ t[(ϵ+p) γ^2-p]
    	+z[(ϵ+p) γ^2β_z]=1/cj·E,




    1/c∂/∂ t[(ϵ+p).γ^2β]
    	+z[(ϵ+p) γ^2β_z β]
        +pz
        =ρ_eE+1/cj×B,


where ρ_e is the charge density and j is the current density.
The electromagnetic field is calculated from the Maxwell equations and the ideal MHD conditions.
Since the electric field satisfies the ideal MHD condition E=-β×B, the evolution of the magnetic field from the Maxwell-Faraday equation is described by the following induction equation:


    1/cBt=∇×(β×B).


The charge density ρ_e and the current density j in the right-hand side of equations (<ref>) and (<ref>), respectively, can be written as follows by using Maxwell's equations:


    ρ_e
    	=1/4π∇·E,




    j
        =c/4π(∇×B-1/cEt).


Note that since the electric field E satisfies the ideal MHD conditions, the charge density ρ_e and the current density j can be obtained from the magnetic field B and velocity β.



 §.§ Derivation of linearized equations

The quantities required in the subsequent calculations up to the order of 𝒪(ε) are calculated as follows.
Square of the Lorentz factor of the fluid:


    γ^2=1/1-β·β
    	
    	∼δγ^2(1+2δγ^2δβ·β_⊥),


where δγ=(1-δβ^2)^-1/2=(1-η^2 β_A,η^2)^-1/2.
Enthalpy density:


    ϵ+p∼(ϵ_0+p_0)+ϵ_(1+β_s^2).


z component of the electric field:


    E_ze_z∼-(δβ×b_⊥+ β_⊥×δB).


Note that E_z is a quantity of 𝒪(ε).
In-plane xy component of the electric field:


    E_xy∼-(δβ×B_0+β_⊥×B_0+β_×δB).


Square of the absolute value of the electric field:


    E^2∼δβ^2 B_0^2
    	+2(δβ·β_⊥)B_0^2
    	+2χ_0δβδ BB_0β_.


Square of the absolute value of the magnetic field:


    B^2
     
    	∼ B_0^2+δ B^2+2δB·b_⊥.



Using equations (<ref>)–(<ref>), calculate the right-hand side of equations (<ref>) and (<ref>).
First, the work rate of the electromagnetic field j·E is as follows:


    j·E∼c/4π[
    	B_0β_⊥·z(δB)
    	+B_0δβ·b_⊥z
    	-B_0^21/ct(δβ·β_⊥)
    	-χ_0δβδ BB_01/cβ_t.
    	]


Here we used the following relations, which hold for circularly polarised Alfvén waves:


    δB·z(δB)=0,  δβ·z(δB)=0.



Next, let us calculate the electric force ρ_e E.
From the equation (<ref>), the charge densityρ_e is as follows.


    ρ_e
        
    	=1/4πE_zze_z
    	∼
    	-1/4πz[
    	B_0/B_0·(δβ×b_⊥+ β_⊥×δB)
    	].


As can be seen from this equation, ρ_e=𝒪(ε).
Since E_z is also a quantity of 𝒪(ε), the electric force in the z direction ρ_e E_z is a second order of ε:


    ρ_e E_z∼ 0.


The xy in-plane component of the electric force is as follows:


    ρ_e E_xy∼1/4π(δβ×B_0)
    	z[
    		B_0/B_0·(δβ×b_⊥+ β_⊥×δB)
    	].



From equation (<ref>), the z component of the magnetic force j×B can be calculated as follows:


    [1/cj×B]_z=
    -∂/∂ z(δB·b_⊥/4π)
    +B_0/4π[
    δ B·1/c∂β_⊥/∂ t
    +b_⊥·1/c∂/∂ t(δβ)
    ]
    -δ B^2/4π1/c∂β_/∂ t.


The xy in-plane component of the magnetic force is as follows:


    [1/cj×B]_xy=
    B_0/4π∂/∂ z(δB+b_⊥)
    -B_0^2/4π1/c∂/∂ t(δβ)
    -B_0^2/4π1/c∂β_⊥/∂ t
    +B_0/4π1/c∂/∂ t(β_δB)
    -δB·b_⊥/4π1/c∂/∂ t(δβ)
    
    -δβ(δB/4π·1/c∂b_⊥/∂ t)
    -χ_0δ B/4πδβ1/c∂b_⊥/∂ t
    -δ B^2/4π1/c∂β_⊥/∂ t
    +δB/4π(1/c∂β_⊥/∂ t·δB)
    +δB·β_⊥/4π1/c∂/∂ t(δB).


By substituting the above equations (<ref>)–(<ref>) into equations (<ref>), (<ref>) and (<ref>), equations (<ref>)–(<ref>) can be obtained by straightforward calculation.










