Dispersion relations

Dispersion relations are a valuable tool for code verification.

Base classes

Base classes for dispersion relations.

class struphy.dispersion_relations.base.DispersionRelations1D(*branch_names, velocity_scale='alfvén', **params)

Bases: object

The base class for analytic 1d dispersion relations of the form \(\omega_i(k) \in \mathbb C\), where \(i=1,\ldots,n\) denote the differnt branches of the spectrum.

Parameters:

• branch_names (str) – Branche names (str) of the spectrum.

• velocity_scale (str) – Determines the unit of \(\omega\). Must be one of alfvén, cyclotron or light.

Notes

Analytic dispersion relations should be added to analytic.

property branches

Dictionary of branch names holding the numpy arrays of \(\omega_i(k) \in \mathbb C\).

property params

Dictionary of parameters necessary to compute the dispersion relation.

property velocity_scale

Determines the unit of \(\omega\). Must be one of alfvén, cyclotron or light.

property k_crit

Dictionary of critical k-values (plotted as vertical lines with self.plot()).

plot(k)

class struphy.dispersion_relations.base.ContinuousSpectra1D(*branch_names, **params)

Bases: object

The base class for analytical continuous spectra in one spatial dimension.

Parameters:

• *branch_names – Names (str) of the continuous spectra.

• **params – Physical parameters necessary to compute the continuous spectra (e.g. profiles for magnetic field components).

property branches

List of branch names in the spectrum.

property nbranches

number of branches.

Type:

Integer

property params

Dictionary of parameters necessary to compute the continuous spectra.

Available dispersion relations

Analytic dispersion relations.

class struphy.dispersion_relations.analytic.Maxwell1D(c=1.0)

Bases: DispersionRelations1D

Dispersion relation for Maxwell’s equation in vacuum:

\[\omega = c k \,.\]

Parameters:

c (float) – Speed of light. Remark that \(c=1.0\) in Struphy units, see Maxwell.

class struphy.dispersion_relations.analytic.MHDhomogenSlab(B0x=0.0, B0y=1.0, B0z=1.0, p0=0.1, n0=1.0, gamma=1.6666666666666667)

Bases: DispersionRelations1D

Dispersion relation for linear MHD equations in homogeneous background \((n_0,p_0,\mathbf B_0)\) and wave propagation along \(z\)-axis in Struphy units (see LinearMHD):

\[ \begin{align}\begin{aligned}\textnormal{shear Alfvén}:\quad &\omega = v_\textnormal{A} k\frac{B_{0z}}{|\mathbf B_0|}\,,\\\textnormal{fast (+) and slow (-) magnetosonic}:\quad &\omega = k\sqrt{\frac{1}{2}(c_\textnormal{S}^2+v_\textnormal{A}^2)(1\pm\sqrt{1-\delta})}\,,\quad\delta=\frac{4B_{0z}^2c_\textnormal{S}^2v_\textnormal{A}^2}{(c_\textnormal{S}^2+v_\textnormal{A}^2)^2|\mathbf B_0|^2}\,,\end{aligned}\end{align} \]

where \(v_\textnormal{A}^2=|\mathbf B_0|^2/n_0\) is the Alfvén velocity and \(c_\textnormal{S}^2=\gamma\,p_0/n_0\) is the speed of sound.

Parameters:

• B0x (float) – x-component of magnetic field (default: 0.).

• B0y (float) – y-component of magnetic field (default: 1.).

• B0z (float) – z-component of magnetic field (default: 1.).

• p0 (float) – Plasma pressure (default: 0.1).

• n0 (float) – Ion number density (default: 1.).

• gamma (float) – Adiabatic index (default: 5/3).

class struphy.dispersion_relations.analytic.ExtendedMHDhomogenSlab(B0x=0.0, B0y=0.0, B0z=1.0, p0=0.1, n0=1.0, gamma=1.6666666666666667, eps=0.1)

Bases: DispersionRelations1D

Dispersion relation for linear extended MHD equations in homogeneous background \((n_0,p_0,\mathbf B_0)\) and wave propagation along \(z\)-axis in Struphy units (see LinearExtendedMHD). The linear mode analysis is performed on the following system:

\[\begin{split}\begin{align} -i \omega m n_0 \tilde{\mathbf U} + i \mathbf k \tilde p &= i\frac{B_0}{\mu_0}(\mathbf k \times \tilde{\mathbf B}) \times \mathbf e_z\,, \\ -i\omega \tilde{\mathbf B} - B_0 i \mathbf k \times ( \tilde{\mathbf U} \times \mathbf e_z ) \color{red} + i \mathbf k \times \left[\frac{B_0}{q n_0\mu_0} (i\mathbf k \times \tilde{\mathbf B} ) \times \mathbf e_z \right] \color{default} &= 0\,, \\[3mm] -i\omega \tilde p + i \gamma p_0 \mathbf k \cdot \tilde{\mathbf U} &= 0\,. \end{align}\end{split}\]

Computed are the four roots of

\[(x - \omega_0^2) \Big\{ x^3 - x^2\Big[(c_s^2 + v_A^2)k^2 + v_A^2 k_z^2 + \omega_0^2\Big] + x\Big[c_s^2 k^2 (2 v_A^2 k_z^2 + \omega_0^2) + v_A^4 k_z^2 k^2 \Big] - c_s^2 v_A^4 k_z^4 k^2\Big\} = 0\]

where \(x = \omega^2\), \(v_\textnormal{A}^2=|\mathbf B_0|^2/n_0\) denotes the Alfvén velocity, \(c_s^2=\gamma\,p_0/n_0\) is the speed of sound, and \(\omega_0 = v_A^2 k_z k/\Omega_i\) stands for the first root where \(\Omega_i = |\mathbf B_0|/\epsilon\) is the cyclotron frequency with \(\epsilon = 1/(\hat \Omega_i \hat t)\).

class struphy.dispersion_relations.analytic.FluidSlabITG(vstar=10.0, vi=1.0, Z=1.0, kz=1.0, gamma=1.6666666666666667)

Bases: DispersionRelations1D

Dispersion relation for ion fluid equations with adiabatic electrons in homogeneous background \((n_0,p_0,\mathbf B_0)\), with a constant density and a temperature gradient in \(x\)-direction. The Braginskii closure for drift cancellations is taken into account (strong magentic field). The dispersion relation is

\[\omega^3 - (Z + \gamma)v_i^2 k_z^2 \omega + Z \frac{v_i^4}{v^*} k_z^2 k_y = 0\,,\]

where \(v_i^2=k_B T_0/m\) denotes the ion thermal velocity and

\[v^* = \frac{\Omega_i}{\partial_x T_0/T_0}\,,\qquad \Omega_i = \frac{Ze B_0}{m}\,.\]

The dispersion relation is calculated as a function of \(k_y\) for fixed \(k_z\). Instabilites occur when

\[\frac{k_y}{k_z} > \frac{2}{3^{3/2}}\frac{\sqrt{Z + \gamma}^3}{Z^2} \frac{v^*}{v_i} \,, \qquad \frac{v^*}{v_i} = \frac{\partial_x T_0/T_0}{\rho_i}\,.\]

class struphy.dispersion_relations.analytic.ColdPlasma1D(**params)

Bases: DispersionRelations1D

Dispersion relation for cold plasma model for homogeneous background \((n_0,\mathbf B_0)\) and wave propagation along z-axis \((\mathbf k = k \mathbf e_z)\) in Struphy units (see ColdPlasma in Models):

\[\left[ \left( \omega^2 - |k|^2 \right) \mathbb I + \mathbf k \otimes \mathbf k + i \frac{\alpha^2}{\varepsilon_c} \omega \sigma_c \right] \mathbf E = 0\,,\]

where \(\left( \omega^2 - |k|^2 \right) \mathbb I + \mathbf k \otimes \mathbf k + i \frac{\alpha^2}{\varepsilon_c} \omega \sigma_c = \epsilon\) is the dielectric tensor, \(\alpha\) is the plasma frequency in units of the electron cyclotron frequency, \(1/\varepsilon_c\) is the electron cyclotron frequency in struphy units, and \(\sigma_c = \left( \mathbb I - i Q / \varepsilon_c \omega \right)^{-1} i n_0 / \varepsilon_c \omega\), with \(Q\) being an operator which, if applied to vector \(\mathbf v\), returns \(\mathbf v \times \mathbf B_0\).

class struphy.dispersion_relations.analytic.CurrentCoupling6DParallel(**params)

Bases: DispersionRelations1D

Dispersion relation for linearized hybrid MHD-Vlasov model (current coupling scheme) in Struphy units for homogeneous background \((n_0=1,p_0,\mathbf B_0=B_0\mathbf e_z)\), wave propagation along z-axis and EP distribution function

\[f_0=\frac{1}{\pi^{3/2}v_\textnormal{th}^3}\exp\left[-\frac{(v_\parallel-v_0)^2+v_\perp^2}{v_\textnormal{th}^2}\right]\,.\]

The two branches of the dispersion relation are circularly polarized shear Alfvén waves (R-wave and L-wave)

\[D_\textnormal{R/L}(\omega,k)=\omega^2-B_0^2k^2\pm\nu_\textnormal{h}\omega\frac{Z_\textnormal{h}B_0}{A_\textnormal{b}}\kappa\mp\nu_\textnormal{h}\frac{Z_\textnormal{h}B_0}{A_\textnormal{b}}\kappa v_0 k+\nu_\textnormal{h}\frac{Z_\textnormal{h}^2B_0^2}{A_\textnormal{h}A_\textnormal{b}}\kappa^2\frac{\omega-kv_0}{kv_\textnormal{th}}Z(\xi^\pm)\,,\]

and standard sound waves

\[\omega^2=\gamma p_0 k^2\,,\]

where \(\xi^\pm=(\omega-kv_0\pm Z_\textnormal{h}B_0\kappa/A_\textnormal{h})/kv_\textnormal{th}\) and \(Z(\xi)=\sqrt{\pi}\exp(-\xi^2)[i-\textnormal{erfi}(\xi)]\) is the plasma dispersion function.

Parameters:

**params –

Keyword arguments that characterize the dispersion relation.

• B0float

  Magnetic field strength (default: 1.).

• p0float

  Plasma pressure (default: 0.5).

• gammafloat

  Adiabatic index (default: 5/3).

• Abint

  Bulk species mass number (default: 1).

• Ahint

  Energetic species mass number (default: 1).

• Zhint

  Energetic species charge number (default: 1).

• vthfloat

  Energetic species thermal velocity (default: 1.).

• v0float

  Energetic species shift of Maxwellian (default: 2.).

• nuhfloat

  Ratio of energetic/bulk number densities (default: 0.5).

• nbfloat

  Bulk species number density in units of 1e20 / m^3 (default: 0.0005185219355).

D_RL(w, k, pol, der=0)

Dispersion relation \(D_\mathrm{R/L}(\omega,k)=0\) (or its first derivative with respect to \(\omega\)) for R- and L- shear Alfvén waves.

Parameters:

• w (list) – The complex frequencies at which to evaluate the dispersion relation. w[0] is the real part, w[1] the imaginary part.

• k (array_like) – The real wave numbers at which to evaluate the dispersion relation.

• pol (int) – The polarization of the wave (+1 : R-wave, -1 : L-wave).

• der (int, optional) – Whether to evaluate the dispersion relation (der = 0) or its first derivative with respect to w (der = 1).

Returns:

• d_real (ndarray) – The real part of the evaluated dispersion relation.

• d_imag (ndarray) – The imaginary part of the evaluated dispersion relation.

class struphy.dispersion_relations.analytic.PressureCouplingFull6DParallel(params)

Bases: DispersionRelations1D

Dispersion relation for linear MHD equations coupled to the Vlasov equation with Full Pressure Coupling scheme for homogeneous background \((n_0,p_0,\mathbf B_0)\), wave propagation along z-axis in Struphy units and space-homogeneous shifted Maxwellian energetic particles distribution \(f_h = f_{h,0} + \tilde{f_h}\) where \(f_{h,0}(v_{\parallel}, v_{\perp}) = n_0 \frac{1}{\sqrt{\pi}} \frac{1}{\hat{v_{\parallel}}} e^{- (v_{\parallel} - u_0)^2 / \hat{v}^2_{\parallel} } \frac{1}{\pi} \frac{1}{\hat{v^2_{\perp}}} e^{- v^2_{\perp} / \hat{v}^2_{\perp}}\) here, \(u_0\) is a velocity shift in the parallel direction (see PC_LinMHD_6d_full in Models):

\(\textnormal{shear Alfvén (R) and (L) wave}\) :

\[ \begin{align}\begin{aligned}\omega^2 = v_\textnormal{A}^2 k^2\frac{B_{0z}^2}{|\mathbf B_0|^2} + \omega k \nu_h &\left[ \frac{\omega_c}{\omega} \left\{ \left( 1 - \frac{\hat{v}^2_\perp}{\hat{v}^2_\parallel}\right) \hat{v}_\parallel \left( \pm Y_3 \mp \frac{\omega - \omega_c}{\hat{v}_\parallel k_\parallel} Y_2 \right) + u_0 \frac{\hat{v}^2_\perp}{\hat{v}^2_\parallel} \left( \pm Y_2 \mp \frac{\omega - \omega_c}{\hat{v}_\parallel k_\parallel} Y_2 \right) \right\} \right.\\&- \left. \frac{\hat{v}^2_\perp}{\hat{v}^2_\parallel} \left( Y_3 - \frac{\omega \mp \omega_c}{\hat{v}_\parallel k_\parallel} Y_2 - \frac{u_0}{\hat{v}_\parallel} Y_2 + \frac{\omega \mp \omega_c}{\hat{v}^2_\parallel k_\parallel} u_0 Y_1 \right)\right]\,,\end{aligned}\end{align} \]

\(\textnormal{sonic wave}\) :

\[\omega^2 =v_\textnormal{A}^2 k^2 - 2 \omega k_\parallel \nu_h \hat{v}_\parallel X_4 \,\]

where \(v_\textnormal{A}^2=|\mathbf B_0|^2/n_0\) is the Alfvén velocity and \(c_\textnormal{S}^2=\gamma\,p_0/n_0\) is the speed of sound.

Variaous integrals are defined as follows

\[ \begin{align}\begin{aligned}X_4(\xi_0, a_0):= \frac{1}{\sqrt{\pi}} \int^\infty_\infty \frac{(t+a)^3 t }{t - \xi_0} e^{- t^2} dt \, \quad \qquad &= \frac{5}{4} \xi_0 + \frac{3}{2} a_0 + (\xi_0 + a_0)^3 [1 + \xi_0 Z(\xi_0)] \,,\\Y_1(\xi_-, \xi_+, a_0) := \frac{1}{\sqrt{\pi}} \int^\infty_\infty \frac{t+a_0}{(t-\xi_-)(t-\xi_+)} e^{-t^2} dt &= Z(\xi_-) + (\xi_+ + a_0) \frac{Z(\xi_-) - Z(\xi_+)}{\xi_- - \xi_+} \,,\\Y_2(\xi_-, \xi_+, a_0) := \frac{1}{\sqrt{\pi}} \int^\infty_\infty \frac{(t+a)^2}{(t-\xi_-)(t-\xi_+)} e^{-t^2} dt &= 1 + (\xi_- + \xi_+ + 2a_0) Z(\xi_-)\\&+ (\xi_+ + a_0)^2 \frac{Z(\xi_-) - Z(\xi_+)}{\xi_- - \xi_+} \,,\\Y_3(\xi_-, \xi_+, a_0) := \frac{1}{\sqrt{\pi}} \int^\infty_\infty \frac{(t+a)^3}{(t-\xi_-)(t-\xi_+)} e^{-t^2} dt &= \xi_- + \xi_+ + 3a_0\\&+ [\xi_-^2 + \xi_- \xi_+ + \xi_+^2 + 3a_0(\xi_- + \xi_+) + 3a_0^2] Z(\xi_-)\\&+ (\xi_+ + a_0)^3 \frac{Z(\xi_-) - Z(\xi_+)}{\xi_- - \xi_+} \,,\end{aligned}\end{align} \]

where \(\xi_0 = \frac{\omega / k_\parallel - u_0}{\hat{v}_\parallel}, \quad \xi_\pm = \frac{(\omega \pm \omega_c) / k_\parallel - u_0}{\hat{v}_\parallel}, \quad a_0 = \frac{u_0}{\hat{v}_\parallel}\) and \(Z(\xi) = \frac{1}{\sqrt{\pi}} \int^\infty_\infty \frac{e^{- t^2}}{t - \xi} dt = i \sqrt{\pi} e^{- \xi^2} ( 1 + \text{erf}(i\xi))\) is the plasma dispersion function.

D_RL(w, k, pol)

Dispersion relation \(D_\mathrm{R/L}(\omega,k)=0\) for R- and L- shear Alfvén waves.

Parameters:

• w (list) – The complex frequencies at which to evaluate the dispersion relation. w[0] is the real part, w[1] the imaginary part.

• k (array_like) – The real wave numbers at which to evaluate the dispersion relation.

• pol (int) – The polarization of the wave (+1 : R-wave, -1 : L-wave).

Returns:

• d_real (ndarray) – The real part of the evaluated dispersion relation.

• d_imag (ndarray) – The imaginary part of the evaluated dispersion relation.

D_sonic(w, k)

Dispersion relation \(D_\mathrm{sonic}(\omega,k)=0\) for sonic waves.

Parameters:

• w (list) – The complex frequencies at which to evaluate the dispersion relation. w[0] is the real part, w[1] the imaginary part.

• k (array_like) – The real wave numbers at which to evaluate the dispersion relation.

Returns:

• d_real (ndarray) – The real part of the evaluated dispersion relation.

• d_imag (ndarray) – The imaginary part of the evaluated dispersion relation.

class struphy.dispersion_relations.analytic.MhdContinousSpectraShearedSlab(**params)

Bases: ContinuousSpectra1D

Continuous shear Alfvén and slow sound spectra along x-direction in slab geometry with side lengths \(L_x=a,L_y=2\pi a, L_z=2\pi\,R_0\) in Struphy units.

The profiles in Cartesian coordinates \((x, y, z)\) are

\[ \begin{align}\begin{aligned}\mathbf B_0 &= \mathbf B_0(x) = B_{0z}(x)\left( \mathbf e_z + \frac{a}{q(x) R_0}\mathbf e_y \right)\,,\\p_0 &= p_0(x)\\n_0 &= n_0(x)\,.\end{aligned}\end{align} \]

The continuous spectra are then given by

\[ \begin{align}\begin{aligned}\textnormal{shear Alfvén}:\quad & \omega^2(x)=\frac{B_{0z}(x)^2}{n_0(x)}\frac{1}{R_0^2}\left(n+\frac{m}{q(x)}\right)^2\,\\\textnormal{slow sound}:\quad & \omega^2(x)=\frac{\gamma p_0(x)B_{0z}(x)^2}{n_0(x)\,[\gamma p_0(x) + B_{0y}(x)^2 + B_{0z}(x)^2]}\frac{1}{R_0^2}\left(n+\frac{m}{q(x)}\right)^2\,.\end{aligned}\end{align} \]

Parameters:

**params –

Keyword arguments that characterize the dispersion relation.

• afloat

  ”Minor” radius (must be compatible with \(L_x=a\) and \(L_y=2\pi a\), default: 1.).

• R0float

  ”Major” radius (must be compatible with \(L_z=2\pi R_0\), default: 3.).

• gammafloat

  Adiabatic index (default: 5/3).

• Bzcallable

  Profile of axial magnetic field Bz=Bz(x) (default: 1. - 0*x).

• pcallable

  Pressure profile p=p(x) (default: 0.5 - 0*x).

• rhocallable

  Profile of mass density rho=rho(x) (default: 1. - 0*x).

• qcallable

  Safety factor profile q=q(x) (default: 1.1 + 0.7*x**2).

class struphy.dispersion_relations.analytic.MhdContinousSpectraCylinder(**params)

Bases: ContinuousSpectra1D

Continuous shear Alfvén and slow sound spectra along radial direction in cylindrical geometry with radius \(a\) and length \(2\pi\,R_0\) in Struphy units.

The profiles in cylindrical coordinates \((r, \theta, z)\) are

\[ \begin{align}\begin{aligned}\mathbf B_0 &= \mathbf B_0(r) = B_{0z}(r)\left( \mathbf e_z + \frac{r}{q(r) R_0}\mathbf e_\theta \right)\,,\\p_0 &= p_0(r)\\n_0 &= n_0(r)\,.\end{aligned}\end{align} \]

The continuous spectra are then given by

\[ \begin{align}\begin{aligned}\textnormal{shear Alfvén}:\quad & \omega^2(r)=\frac{B_{0z}(r)^2}{n_0(r)}\frac{1}{R_0^2}\left(n+\frac{m}{q(r)}\right)^2\,\\\textnormal{slow sound}:\quad & \omega^2(r)=\frac{\gamma p_0(r)B_{0z}(r)^2}{n_0(r)\,[\gamma p_0(r) + B_{0\theta}(r)^2 + B_{0z}(r)^2]}\frac{1}{R_0^2}\left(n+\frac{m}{q(r)}\right)^2\,.\end{aligned}\end{align} \]

Parameters:

**params –

Keyword arguments that characterize the dispersion relation.

• R0float

  ”Major” radius (must be compatible with \(L_z=2\pi R_0\), default: 3.).

• gammafloat

  Adiabatic index (default: 5/3).

• Bzcallable

  Profile of axial magnetic field Bz=Bz(x) (default: 1. - 0*r).

• pcallable

  Pressure profile p=p(x) (default: 0.5 - 0*r).

• rhocallable

  Profile of mass density rho=rho(x) (default: 1. - 0*r).

• qcallable

  Safety factor profile q=q(x) (default: 1.1 + 0.7*r**2).