2.1. Initial value problem

In the initial value problem, the spatial domain is infinite, and the wave envelopes evolve in time from their initial conditions

(2.8)\begin{equation} \boldsymbol{\alpha}(x,t=0) = \boldsymbol{A}(x), \end{equation}

where $\boldsymbol {A}=(A_2, A_3)^\mathrm {T}$. In the degenerate case $v_2=v_3=v$, after transforming to the comoving frame $x'=x-vt$ and $t'=t$ where $\partial _t+v\partial _x=\partial _{t'}$, the equations become ordinary differential equations (ODEs) in $t'$. The eigenvalues $\gamma _\pm =-\bar {\mu }\pm \varOmega$ of the linear ODEs are the $\kappa =0$ roots of (2.5), and the general solutions are of the form $\alpha (x',t')=A_+(x')\exp ({\gamma _+t'}) + A_-(x')\exp ({\gamma _-t'})$. The coefficients $A_\pm$ are determined by matching initial conditions, which give the solution

(2.9)\begin{equation} \boldsymbol{\alpha}=\exp(-\bar{\mu}t) \left(\begin{array}{@{}cc@{}} \cosh \varOmega t-\dfrac{\Delta\mu}{2\varOmega} \sinh \varOmega t & \dfrac{\gamma_0}{\varOmega} \sinh \varOmega t \\ \dfrac{\gamma_0}{\varOmega} \sinh \varOmega t & \cosh \varOmega t+ \dfrac{\Delta\mu}{2\varOmega} \sinh \varOmega t \end{array}\right)\boldsymbol{A}(x-vt). \end{equation}

The $\Delta v=0$ solution exhibits two features that are more general: a diagonal damping term and off-diagonal coupling terms that vanish when $\gamma _0=0$. The rest of this paper focuses on the non-degenerate case $\Delta v>0$, because in LPIs $v_2$ is close to the speed of light whereas $v_3$ is on the scale of thermal velocities, which are much slower.

In the non-degenerate case $\Delta v> 0$, because the equations are linear, they can be solved using Fourier transform $\mathcal {F}[p](k)=\int _{-\infty }^{+\infty } {{\rm d}x}\,p(x) \exp ({-{\rm i}kx})=\hat {p}(k)$. In momentum space, the equation becomes $\partial _t\hat {\boldsymbol {\alpha }}={\boldsymbol{\mathsf{K}}} \hat {\boldsymbol {\alpha }}$ where

(2.10)\begin{equation} {\it{\boldsymbol{\mathsf{K}}}} = \left(\begin{array}{@{}cc@{}} -{\rm i}kv_2-\mu_2 & \gamma_0 \\ \gamma_0 & -{\rm i}kv_3 -\mu_3 \end{array}\right). \end{equation}

As the matrix is time independent, the solution is $\hat {\boldsymbol {\alpha }}(t)=\exp (t {\boldsymbol{\mathsf{K}}} ) \hat {\boldsymbol {A}}$. The matrix exponential can be computed by diagonalising $\it{\boldsymbol{\mathsf{K}}}$ whose eigenvalues are $\lambda _\pm =-{\rm i}k\bar {v}-\bar {\mu }\pm \sqrt {\gamma _0^2-\omega ^2}$, where $\bar {v}=\frac {1}{2}(v_2+v_3)$ and $\omega =\frac {1}{2}(k\Delta v-{\rm i}\Delta \mu )$. Finding eigenvectors of $\lambda _\pm$ and diagonalising $\it{\boldsymbol{\mathsf{K}}}$, the solution map $\hat {\it{\boldsymbol{\varPhi }}}(k,t)=\exp (t \it{\boldsymbol{\mathsf{K}}} )$ can be written as

(2.11)\begin{equation} \hat {\it{\boldsymbol{\varPhi }}}(k,t)= \exp({-({\rm i}k\bar{v}+\bar{\mu})t}) \left(\begin{array}{@{}cc@{}} \partial_t-{\rm i}\omega & \gamma_0 \\ \gamma_0 & \partial_t+{\rm i}\omega \end{array}\right)\hat{G}(k,t), \end{equation}

where $\hat {G}(k,t) = \sinh (t\sqrt {\gamma _0^2-\omega ^2})/\sqrt {\gamma _0^2-\omega ^2}$. The solution map $\hat {\it{\boldsymbol{\varPhi }}}(k,t)$ satisfies matrix equation $\partial _t \hat {\it{\boldsymbol{\varPhi }}}(k,t)={\boldsymbol{\mathsf{K}}} \hat {{\boldsymbol{\varPhi }}}(k,t)$ and initial condition $\hat {\it{\boldsymbol{\varPhi }}}(k,t=0)=\it{\boldsymbol{\mathsf{I}}}$, where $\it{\boldsymbol{\mathsf{I}}}$ is the two-dimensional identity matrix. Note that the behaviour of $\hat {G}$ is regular when $\omega (k)\rightarrow \pm \gamma _0$. We can choose the branch cut for the square roots to lie between these two points. As $\sinh$ is an odd function, $\hat {G}$ is, in fact, analytic in the complex $k$ plane.

To find the solution in $x$ space, take inverse Fourier transform $\mathcal {F}^{-1}[\hat {p}](x) =\int _{-\infty }^{+\infty } ({{\rm d} k}/{2{\rm \pi} }) \hat {p}(k) \exp({\rm i}kx) =p(x)$. As $\omega$ depends on $k$, it is convenient to change the integration variable to $k'=2\omega /\Delta v$, which gives $k=k'+{\rm i}\Delta \mu /\Delta v$. Moreover, it is convenient to change the reference frame to

(2.12a)\begin{gather} \tau = \frac{\Delta v}{2} t, \end{gather}

(2.12b)\begin{gather}z = x-\bar{v}t, \end{gather}

which travels at the average velocity of the two waves. Note that $z-\tau =x-v_2t$ and $z+\tau =x-v_3t$. The phase factor simplifies as $\exp [{\rm i}kx-({\rm i}k\bar {v}+\bar {\mu })t]=\rho \exp ({{\rm i}k'z})$ where

(2.13)\begin{equation} \rho = \exp\left\{\frac{1}{\Delta v}[\mu_3(z-\tau)-\mu_2(z+\tau)]\right\}, \end{equation}

which is a damping factor independent of $k'$. When taking inverse Fourier transform of (2.11), all matrix elements can be expressed in terms of

(2.14a)\begin{align} g(z, \tau) & = \int_{-\infty-{\rm i}\epsilon}^{+\infty-{\rm i}\epsilon} \frac{{\rm d} k'}{2{\rm \pi}} \exp({{\rm i}k'z}) \frac{\sinh \tau\sqrt{m^2-k'^2}}{\sqrt{m^2-k'^2}} \end{align}

(2.14b)\begin{align} & = \frac{1}{2}\mathrm{sign}(\tau)I_0(m\sqrt{\tau^2-z^2})\varTheta(\tau^2-z^2), \end{align}

where $m=2\gamma _0/\Delta v$, $I_0$ is modified Bessel function and $\varTheta$ is Heaviside step function. The shift by $\epsilon = \Delta \mu /\Delta v$ is insignificant because the integrand is analytic. The above integral is calculated in Appendix A. The solution map $\it{\boldsymbol{\varPhi }}$ in configuration space is

(2.15)\begin{equation} \it{\boldsymbol{\varPhi}}(x,t)= \rho\left(\begin{array}{@{}cc@{}} \partial_\tau-\partial_z & m \\ m & \partial_\tau+\partial_z \end{array}\right) g(z, \tau). \end{equation}

The derivatives are evaluated using $I'_0(\xi )=I_1(\xi )$ and $\varTheta '(\xi )=\delta (\xi )$, and an explicit formula is given by (A7). The function $g$ satisfies the imaginary-mass Klein–Gordon equation $(\partial _\tau ^2-\partial _z^2-m^2)g(z,\tau )=0$, as shown in (A5). Consequently, the solution map satisfies $\it{\boldsymbol{\mathsf{L}}} \it{\boldsymbol{\varPhi}} (x,t)=\it{\boldsymbol{\mathsf{0}}}$ following (A6). The step function $\varTheta$ enforces the causality that information outside the light cone $\tau ^2-z^2=(v_2t-x)(x-v_3t)$ does not affect solutions.

Finally, to invert $\hat {\boldsymbol {\alpha }}=\hat {\it{\boldsymbol{\varPhi }}}\hat {\boldsymbol {A}}$, compute inverse Fourier transform of products, which are given by convolutions $\mathcal {F}^{-1}[\hat {p}\hat {q}](x)=\int _{-\infty }^{+\infty } {{\rm d}x}'\, p(x')q(x-x')$. As the phase factor is simpler in $k'=k-{\rm i}\Delta \mu /\Delta v$, in addition to the change of variables in (2.12), it is convenient to define $\hat {\boldsymbol {B}}(k')= \hat {\boldsymbol {A}}(k)$, which means that $\boldsymbol {B}(x) = \exp (x\Delta \mu /\Delta v)\boldsymbol {A}(x)$. With $\boldsymbol {\alpha }=\rho \boldsymbol {\beta }$, where $\rho$ is given by (2.13), the solution can be written as

(2.16a)\begin{gather} \beta_2(z,\tau) = B_2(z-\tau)+\frac{m}{2}\int_{-\tau}^{\tau}{\rm d} z' \left[B_2(z-z')\sqrt{\frac{\tau+z'}{\tau-z'}}I_1(\xi')+B_3(z+z')I_0(\xi')\right], \end{gather}

(2.16b)\begin{gather}\beta_3(z,\tau) = B_3(z+\tau)+\frac{m}{2}\int_{-\tau}^{\tau}{\rm d} z' \left[B_3(z+z')\sqrt{\frac{\tau+z'}{\tau-z'}}I_1(\xi')+B_2(z-z')I_0(\xi')\right], \end{gather}

where $\xi '=m\sqrt {\tau ^2-z'^2}$. An explicit expression of $\boldsymbol {\alpha }$ in $(x,t)$ coordinate is given later in (2.19). It is straightforward to check that the expression satisfies the differential equation and the initial conditions. When $\boldsymbol {B}$ only involves $\delta$ or step functions, the integrals can be readily evaluated. However, for general initial conditions, closed-form analytical expressions may not exist, so the integrals are evaluated numerically. Compared with numerical integration of the differential equations, which advances initial conditions step by step, the integral solution allows fast forwarding, which directly gives the solution at desired final time. Note that by rescaling $z'=\tau \zeta$, the numerical integration is always within the range $\zeta \in [-1,1]$, so the cost of evaluating the integral does not increase with time for smooth initial conditions.

2.2. Initial boundary value problem

Compared with the initial value problem, the boundary value problem is easier to set up in kinetic simulations. In the initial value problem, the distribution functions need to be specified meticulously in both the configuration space and the velocity space in order to ensure that only the desired eigenmodes are excited. In comparison, in a boundary value problem, only electromagnetic fields at the boundary need to be specified. Using wave frequencies and polarisations to select excited waves, the desired eigenmodes then propagate into the initially Maxwellian plasma where interactions occur. With a single boundary at $x=0$, the initial boundary value problem is specified by conditions

(2.17a)\begin{gather} \boldsymbol{\alpha}(x>0, t=0) = \it{\boldsymbol{A}}(x), \end{gather}

(2.17b)\begin{gather}\boldsymbol{\alpha}(x=0, t>0) = \it{\boldsymbol{h}}(t), \end{gather}

where $\boldsymbol {h}=(h_2,h_3)^\mathrm {T}$. In the forward-scattering case $v_2>v_3>0$, these conditions can be specified separately. However, in the backscattering case $v_2>0>v_3$, in order for the problem to be well-posed, only two conditions can be specified independently.

As the equation $\it{\boldsymbol{\mathsf{L}}} \boldsymbol {\alpha }=\boldsymbol {0}$ is linear, the initial boundary value problem can be solved using Laplace transform $\mathcal {L}[p](k)=\int _{0}^{+\infty } {{\rm d}x}\,p(x) \exp ({-{\rm i}kx}) =\tilde {p}(k)$. Using the property of Laplace transform that $\tilde {p}'(k)={\rm i}k\tilde {p}(k)-p(0)$, the equation becomes $\partial _t\tilde {\boldsymbol {\alpha }}={\boldsymbol{\mathsf{K}}} \tilde {\boldsymbol {\alpha }}+\it{\boldsymbol {H}}$, where $\it{\boldsymbol{\mathsf{K}}}$ is given by (2.10) and $\it{\boldsymbol {H}}=(v_2h_2,v_3h_3)^\mathrm {T}$. Using the Duhamel's principle, the inhomogeneous ODE is solved by $\tilde {\boldsymbol {\alpha }}(t)=\hat {\it{\boldsymbol{\varPhi }}}(t)\tilde {\boldsymbol {A}}+ \int _0^t {\rm d} t' \,\hat {\it{\boldsymbol{\varPhi }}} (t-t')\boldsymbol {H}(t')$, where $\hat {\boldsymbol{\varPhi }}$ is given by (2.11). To find the solution in configuration space, take inverse Laplace transform $\mathcal {L}^{-1}[\hat {p}\tilde {q}](x)=\int _\mathcal {C} ({{\rm d} k}/{2{\rm \pi} }) \mathcal {F}[p](k)\mathcal {L}[q](k) \exp({\rm i}kx)=\int _0^{+\infty } {{\rm d}x}' p(x-x')q(x')$, where the contour $\mathcal {C}$ runs below all poles. The solution is

(2.18)\begin{equation} \boldsymbol{\alpha}(x,t)=\int_0^{\infty} {{\rm d}x}'\,\it{\boldsymbol{\varPhi}}(x-x',t)\boldsymbol{A}(x')+ \int_0^t {\rm d} t'\,\it{\boldsymbol{\varPhi}}(x,t-t')\boldsymbol{H}(t'), \end{equation}

where $\it{\boldsymbol{\varPhi}}$ is given by (2.15). Using $\it{\boldsymbol{\mathsf{L}}} \it{\boldsymbol{\varPhi}} =\it{\boldsymbol{\mathsf{0}}}$, the expression clearly satisfies $\it{\boldsymbol{\mathsf{L}}} \boldsymbol {\alpha }=\boldsymbol {0}$. Moreover, using the explicit expression of $\it{\boldsymbol{\varPhi}}$ in (A7), it is easy to see that $\it{\boldsymbol{\varPhi}} (x,t=0)=\delta (x) \it{\boldsymbol{\mathsf{I}}}$, so the initial conditions are always satisfied. However, the situation for boundary conditions depends on the sign of $v_3$, as shown in the following.

Using the property that $\it{\boldsymbol{\varPhi}}$ is zero outside the light cone, the above integral solution, which is equivalent to $\boldsymbol {\alpha }(x,t)=\int _{-\infty }^x {{\rm d}x}'\,\it{\boldsymbol{\varPhi}} (x',t)\boldsymbol {A}(x-x')+\int _0^t {\rm d} t'\,\it{\boldsymbol{\varPhi}}(x,t')\boldsymbol {H}(t-t')$, is simplified in the three regions shown in figure 1. First, ahead of the light cone $x>v_2t$, effects of boundary conditions have not arrived, so only the initial conditions contribute. In terms of the rescaled variable $\boldsymbol {\beta }$, the solution is given by (2.16), and in terms of the original variables, the solution when $x>v_2t>0$ is

(2.19)\begin{equation} \boldsymbol{\alpha}(x,t)=\left( \begin{array}{@{}c@{}} \exp({-\mu_2t})A_2(x-v_2t) \\ \exp({-\mu_3t})A_3(x-v_3t) \end{array}\right)+\frac{\gamma_0}{\Delta v} \int_{v_3t}^{v_2t} {{\rm d}x}' \,\rho(x',t) \it{\boldsymbol{\varPhi}}_0 (x',t)\boldsymbol{A}(x-x'), \end{equation}

where $\rho$ is the damping factor given by (2.13) and $\it{\boldsymbol{\varPhi}} _0$ is the kernel within the light cone given in (A7). When $t=0$, the solution clearly satisfies the initial conditions. Second, behind the light cone $x< v_3t$, which is within the domain only when $v_3>0$, effects of the initial conditions have propagated away, so only the boundary conditions contribute. Therefore, when $v_3t>x>0$, the solution is

(2.20)\begin{align} \boldsymbol{\alpha}(x,t)=\left(\begin{array}{@{}c@{}} \exp({-\mu_2x/v_2})h_2(t-x/v_2) \\ \exp({-\mu_3x/v_3})h_3(t-x/v_3) \end{array}\right)+\frac{\gamma_0}{\Delta v}\int_{x/v_2}^{x/v_3} {\rm d} t' \,\rho(x,t') \it{\boldsymbol{\varPhi}} _0(x,t') \boldsymbol{H}(t-t'). \end{align}

The solution clearly satisfies the boundary conditions at $x=0$ for this forward-scattering case. Finally, inside the light cone, both the initial and boundary conditions contribute, and the solution when $v_3t< x< v_2t$ is given by

(2.21)\begin{align} \boldsymbol{\alpha}(x,t) & = \left( \begin{array}{@{}c@{}} \exp({-\mu_2x/v_2})h_2(t-x/v_2) \\ \exp({-\mu_3t})A_3(x-v_3t) \end{array}\right) \nonumber\\ & \quad + \frac{\gamma_0}{\Delta v}\int_{v_3t}^{x} {{\rm d}x}'\, \rho(x',t) \it{\boldsymbol{\varPhi}} _0(x',t) \boldsymbol{A}(x-x') \nonumber\\ & \quad + \frac{\gamma_0}{\Delta v}\int_{x/v_2}^{t} {\rm d} t' \,\rho(x,t') \it{\boldsymbol{\varPhi}} _0(x,t') \boldsymbol{H}(t-t'). \end{align}

In the forward-scattering case $v_3>0$, the future light cone is within the domain, so the initial and boundary conditions can be specified independently. However, in the backscattering case $v_3<0$, the future light cone is intercepted by the boundary. Intuitively, when information propagates towards left and arrives at the boundary, one cannot arbitrarily set values at the boundary.

Figure 1. Solutions are determined by initial and boundary conditions within the past light cone. For $(x_a, t_a)$ ahead of $x>v_2t$, the light cone (red) only intercepts the $x$ axis, so the solution is independent of boundary conditions. For $(x_b, t_b)$ behind $x< v_3t$, the light cone (blue) only intercepts the $t$ axis, so the solution is independent of initial conditions. For $(x_i, t_i)$ within $v_3t< x< v_2t$, the solution depends on both initial and boundary conditions. (a) In the forward-scattering case, initial and boundary conditions can be specified independently. (b) In the backscattering case, initial conditions arrive at the boundary so constraints must be satisfied.

2.3. Backscattering problem

When the initial conditions are zero, as is the case in kinetic simulations, the integral constraints $\boldsymbol {\alpha }(x=0,t>0)=\boldsymbol {h}(t)$ for the backscattering case $v_3<0$ can be solved explicitly. By setting $x=0$ and $\boldsymbol {A}=\boldsymbol {0}$ in (2.21), the constraints can be simplified as $(0, l_3(s))^\mathrm {T}=\int _0^s {\rm d} s'\,\it{\boldsymbol{\varPsi}} _0(s')\boldsymbol {l}(s-s')$, where $s=\gamma t$ is time normalised by $\gamma =2\gamma _0\sqrt {v_2|v_3|}/\Delta v$, $\boldsymbol {l}(s)={\rm e}^{\mu t}\boldsymbol{h}(t)$ is rescaled by damping $\mu =(\mu _3v_2+\mu _2|v_3|)/\Delta v$, and

(2.22)\begin{equation} {\boldsymbol{\varPsi}} _0(s)= \frac{1}{2}\left( \begin{array}{@{}cc@{}} I_1(s) & -\sqrt{\dfrac{|v_3|}{v_2}} I_0(s) \\ \sqrt{\dfrac{v_2}{|v_3|}} I_0(s) & -I_1(s) \end{array}\right). \end{equation}

Note that $\mu t=s\gamma _a/\gamma _0$, where $\gamma _a$ is the absolute instability threshold given by (2.7). As the constraint is a convolution, it becomes a product after taking Laplace transform, which gives $(0, \tilde {l}_3(\omega ))^\mathrm {T} =\tilde {\it{\boldsymbol{\varPsi}}}_0(\omega )\tilde {\boldsymbol {l}}(\omega )$. To compute $\tilde {\it{\boldsymbol{\varPsi}}}_0$, use integral representation of modified Bessel function (DLMF 2022, (10.32.2)) that $I_0(s)=({1}/{{\rm \pi} })\int _{-1}^1 {\rm d} t \exp(-st)/\sqrt {1-t^2}$ and perform the $s$ integral first, which gives $\tilde {I}_0(\sigma ) = 1/\sqrt {\sigma ^2-1}$ where $\sigma ={\rm i}\omega$. Then, $\tilde {I}_1(\sigma ) = \sigma /\sqrt {\sigma ^2-1}-1$ because $I_1(s)=I'_0(s)$. As $I_n(s)\simeq {\rm e}^s/\sqrt {2{\rm \pi} s}$ when $s\rightarrow \infty$, the Laplace transforms converge only when $\rm{Re}(\sigma )>1$. Solving the constraint in frequency domain gives a unique solution

(2.23)\begin{equation} \tilde{l}_3(\omega) = \sqrt{\frac{v_2}{|v_3|}}(\sigma-\sqrt{\sigma^2-1}) \tilde{l}_2(\omega), \end{equation}

which means that if we specify the boundary condition for $\alpha _2$, then the boundary condition for $\alpha _3$ is completely determined. Taking inverse Laplace transform, whose details are shown in Appendix B, the self-consistent boundary condition is

(2.24)\begin{equation} l_3(s)= \sqrt{\frac{v_2}{|v_3|}}\int_0^s {\rm d} s'\frac{I_1(s')}{s'} l_2(s-s'). \end{equation}

When the normalised time $s=0$, the boundary condition $l_3(0)=0$ is initially quiescent. At later time, $\alpha _3$ builds up due to wave growth and advection. As shown in Appendix B, we can also express $l_2$ in terms of $l_3$. However, in kinetic simulations, it is much easier to specify the boundary conditions for electromagnetic waves and let the plasma evolve to fulfill the above boundary condition.

Having expressed $h_3$ in terms of $h_2$, the solution of $\boldsymbol {\alpha }$ is a functional of $h_2$ only. The integrals simplify when $h_2$ is $\delta$ or $\varTheta$ functions. As $\delta$ functions cannot be resolved numerically, let us focus on the case $h_2(t)=h_0\varTheta (t)$, which is used later to set up kinetic simulations. Using (2.24) with $l_2(s)=h_0\exp ({s\gamma _a/\gamma _0}) \varTheta (s)$, $h_3(t)=h_0\sqrt {v_2/|v_3|}\int _0^{\gamma t} {\rm d} s' \exp ({-s'\gamma _a/\gamma _0}) I_1(s')/s'$. Substituting $\boldsymbol {h}$ into (2.21), changing integration variable to $\varphi '=\gamma (t'-x/v_2)$, and rotating the triangular double integral in a way that leads to (B5) gives

(2.25a)\begin{gather} \alpha_2(x< v_2t) = h_0\exp({-\mu_2x/v_2}) \left[1+\frac{1}{2}\int_0^{\gamma(t-x/v_2)} {\rm d}\varphi \exp({-\varphi\gamma_a/\gamma_0}) D_2(\varphi, \vartheta)\right], \end{gather}

(2.25b)\begin{gather}\alpha_3(x< v_2t) = h_0\exp({-\mu_2x/v_2})\sqrt{\frac{v_2}{|v_3|}}\frac{1}{2}\int_0^{\gamma(t-x/v_2)} {\rm d}\varphi\exp({-\varphi\gamma_a/\gamma_0}) D_3(\varphi, \vartheta), \end{gather}

whereas $\boldsymbol {\alpha }(x>v_2t)=\boldsymbol {0}$ ahead of the wave front. In these expressions, $\gamma =2\gamma _0\sqrt {v_2|v_3|}/\Delta v$, $\vartheta =\gamma _0x/\sqrt {v_2|v_3|}$, $D_2(\varphi, \vartheta ) = \sqrt {1+2\vartheta /\varphi }I_1(\sqrt {\varphi ^2+2\varphi \vartheta })-M_2(\varphi, \vartheta )$ and $D_3(\varphi, \vartheta ) = I_0(\sqrt {\varphi ^2+2\varphi \vartheta })-M_3(\varphi, \vartheta )$, where the kernel functions are

(2.26a)\begin{gather} M_2(\varphi, \vartheta) = \int_0^1 {\rm d} r\, I_0(\sqrt{r^2\varphi^2+2r\varphi\vartheta}) \frac{I_1(\varphi(1-r))}{1-r}, \end{gather}

(2.26b)\begin{gather}M_3(\varphi, \vartheta) = \int_0^1 {\rm d} r \frac{I_1(\sqrt{r^2\varphi^2+2r\varphi\vartheta})}{\sqrt{1+2\vartheta/r\varphi}} \frac{I_1(\varphi(1-r))}{1-r}. \end{gather}

The differential properties of the kernel functions (C3) and (C4) ensures that (2.25) satisfies $\it{\boldsymbol{\mathsf{L}}}\boldsymbol {\alpha }=\boldsymbol {0}$. Moreover, the special values $D_2(\varphi, 0)=0$ and $D_3(\varphi, 0)=2I_1(\varphi )/\varphi$ (C1) and (C2) ensure that the boundary conditions are satisfied. Using these special values, we see $\alpha _3\rightarrow +\infty$ when $t\rightarrow +\infty$ at $x=0$ if $\varsigma =\gamma _a/\gamma _0<1$, whereas $\alpha _3\rightarrow h_0\sqrt {v_2/|v_3|}(\varsigma - \sqrt {\varsigma ^2-1})$ if $\varsigma \ge 1$. More generally when $x\ge 0$, the solutions approach steady state if and only if $\gamma _0$ does not exceed the absolute instability threshold. The solutions are of the form $\alpha _2=h_0 \exp ({-\mu _2x/v_2})(1+\Delta _2)$ and $\alpha _3=h_0\exp ({-\mu _2x/v_2})\Delta _3$. Examples of the growth function $\Delta _j$ are plotted in figures 2 and 3.

Figure 2. The growth function $\Delta _j(x)$ in the step-function backscattering problem for (a) the seed laser and (b) the plasma wave at selected time slices when $v_2=3\times 10^8$ m s$^{-1}$, $\gamma _0=10^{13}$ rad s$^{-1}$ and $\mu _2=0$. As time increases, $\Delta _j$ propagates in space and grows in amplitude. The growth is always zero ahead of the wave front $x=v_2t$. Moreover, $\Delta _2(0)$ is always zero due to the boundary condition. In contrast, $\Delta _3(0)$ builds up from zero as time increases. Compared with the dampingless case (red lines), having an appreciable $\mu _3$ (blue lines) reduces the growth. Compared with the $v_3=0$ case (blue lines), having a small $v_3$ (cyan circles) slightly affects the solutions at early time. At later time, the discrepancies build up because $\gamma _0\approx 1.3\gamma _a$ exceeds the absolute instability threshold.

Figure 3. The growth function $\Delta _j(t)$ in the step-function backscattering problem for (a) the seed laser and (b) the plasma wave at $x=20\,\mathrm {\mu }\rm{m}$ when $v_2=3\times 10^8$ m s$^{-1}$, $\mu _3=5\times 10^{12}$ rad s$^{-1}$ and ${\mu _2=0}$. When $v_3=-v_2/10$ (symbols), the absolute instability threshold is $\gamma _a\approx 7.8\times 10^{12}$ rad s$^{-1}$. When $\gamma _0=7\times 10^{12}$ rad s$^{-1}$ (cyan) is below the threshold, the growth approaches steady state. In contrast, when $\gamma _0=10^{13}$ rad s$^{-1}$ (magenta) exceeds the threshold, $\Delta _j$ continues to increase. When $v_3=0$ (lines), $\gamma _a$ becomes infinite so the growth always saturates. The case with $\gamma _0=10^{13}$ rad s$^{-1}$ (red) has a larger steady-state value than the case with $\gamma _0=7\times 10^{12}$ rad s$^{-1}$ (blue). The effects of a small but finite $v_3$ only become significant at later time.

The integral solutions for the step-function problem are greatly simplified when $v_3\rightarrow 0$, which is a good approximation for LPIs where $v_2\gg |v_3|$. Inside the domain, $\vartheta \rightarrow \infty$ but $\gamma \rightarrow 0$, whereas $\gamma \vartheta =2\gamma _0^2x/v_2$ is finite. As the integrands in (2.26) approach zero when $\varphi \rightarrow 0$ at fixed $\varphi \vartheta$, both $M_2$ and $M_3$ become zero. Then, $D_2\rightarrow \xi I_1(\xi )/\varphi$ and $D_3\rightarrow I_0(\xi )$, where $\xi =\sqrt {2\varphi \vartheta }$. Writing integrals in (2.25) in terms of $\xi$ gives

(2.27a)\begin{gather} \Delta_2 = \int_0^\psi {\rm d}\xi\exp({-\nu\xi^2})I_1(\xi), \end{gather}

(2.27b)\begin{gather}\Delta_3 = \frac{v_2}{2\gamma_0 x}\int_0^\psi {\rm d}\xi\, \xi \exp({-\nu\xi^2})I_0(\xi), \end{gather}

where $\psi =2\gamma _0\sqrt {(t-x/v_2)x/v_2}$ and $\nu =\mu _3v_2/4x\gamma _0^2$. Observe that $\psi =\sqrt {\mu _3t_r/\nu }$ where $t_r=t-x/v_2$ is the retarded time since the wave front passes, and $1/4\nu$ is the spatial gain exponent in steady state when $\mu _2=0$. After integration by parts, $\Delta _3=(\gamma _0/\mu _3)[1-I_0(\psi ) \exp ({-\nu \psi ^2})+\Delta _2]$. It is straightforward to verify that (2.27) solves $\it{\boldsymbol{\mathsf{L}}} \boldsymbol {\alpha }=\boldsymbol {0}$ when $v_3=0$. When $t\rightarrow +\infty$, the solutions approach steady states, which are always finite because $\gamma _a\rightarrow \infty$. Using Gaussian integrals of modified Bessel function (DLMF 2022, (10.43.24)), $\Delta _2=\exp (1/4\nu )-1-\mathcal {R}$. As shown in Appendix C, the residual $\mathcal {R}=\int _\psi ^{+\infty }{\rm d}\xi \exp ({-\nu \xi ^2})I_1(\xi )$ decays as $\exp ({-\mu _3t_r})$ when $\mu _3t_r\gg \max (1,\gamma _0^2x/\mu _3v_2)$. The steady states $\alpha _2(x,+\infty )=(\mu _3/\gamma _0)\alpha _3(x,+\infty )=h_0 \,{\rm e}^{\kappa x}$, where $\kappa =(\gamma _0^2/\mu _3-\mu _2)/v_2$, are consistent with linear stability analysis of (2.5). The formulae in (2.27) are further simplified in two limiting cases. (1) When $\nu \rightarrow +\infty$, namely, when spatial gain is negligible, integrals are dominated by values near $\xi \simeq 0$. The growths $\Delta _2\simeq (\gamma _0^2x/\mu _3v_2)(1-\exp ({-\mu _3t_r}))\rightarrow 0$ and $\Delta _3\simeq (\gamma _0/\mu _3)(1-\exp ({-\mu _3t_r}))$. (2) When $\mu _3\rightarrow 0$, the Gaussian weight becomes unity. Using properties of modified Bessel function (DLMF 2022, (10.43.1)), the integrals are evaluated to $\Delta _2=I_0(\psi )-1$ and $\Delta _3=(v_2/2\gamma _0 x)\psi I_1(\psi )$. When $\psi \simeq 0$, which occurs near the boundary or the wave front, $\Delta _2\simeq \gamma _0^2t_rx/v_2$ and $\Delta _3\simeq \gamma _0 t_r$ grow linearly in time. At given $t$, the maximum of $1+\Delta _2$ is attained at $x=\frac {1}{2}v_2 t$, which propagates at half the wave group velocity. The maximum value attained at $\psi =\gamma _0t$ is $I_0(\gamma _0 t)\simeq {\rm e}^{\gamma _0t}/\sqrt {2{\rm \pi} \gamma _0t}$ when $t\rightarrow +\infty$. Although the exponential growth ${\rm e}^{\gamma _0t}$ is intuitive, the suppression by $1/\sqrt {2{\rm \pi} \gamma _0t}$ is perhaps not one would naïvely expect from linear instability analysis.

3.1. Launching linear eigenmodes

As three-wave equations describe amplitudes of linear eigenmodes, the lasers need to be launched with specified polarisations to excite targeted eigenmodes only. Although polarisation matching is trivial for unmagnetised plasmas, where the two electromagnetic eigenmodes are degenerate, special care needs to be taken in magnetised cases where the R and L elliptically polarised eigenmodes are non-degenerate. If the polarisation is not matched properly, both R and L waves will be excited, giving rise to four polarisation combinations with different couplings (Shi & Fisch Reference Shi and Fisch2019). As incident lasers are purely transverse in the vacuum region, only the perpendicular components need to be matched with plasma eigenmodes. Denoting $\psi$ the elliptical polarisation angle and $\theta$ is the wave phase,

(3.1)\begin{equation} \boldsymbol{E}_\perp{\propto}\, \hat{\boldsymbol{y}}\sin\psi\sin\theta + \hat{\boldsymbol{k}}\times\hat{\boldsymbol{y}}\cos\psi\cos\theta, \end{equation}

where $\hat {\boldsymbol {y}}$ and $\hat {\boldsymbol {k}}$ are unit vectors. The expression is symmetric about the $x$–$z$ plane where $\boldsymbol {B}_0$ lies. The special value $\psi =0$ means that the wave is linearly polarised along $\hat {\boldsymbol {k}}\times \hat {\boldsymbol {y}}$, whereas $\psi ={\rm \pi} /2$ means linear polarisation along $\hat {\boldsymbol {y}}$, $\psi ={\rm \pi} /4$ is R circular polarisation about $\hat {\boldsymbol {k}}$ and $\psi =-{\rm \pi} /4$ is L circular polarisation.

Due to numerical dispersion and PIC density fluctuations, the polarisation angles of numerical eigenmodes are slightly different from analytic results. Launching waves with analytical $\psi$ leads to a small but sometimes observable polarisation rotation of $\boldsymbol {E}_\perp$, revealing contaminations from the unintended eigenmode. To determine $\tilde {\psi }$ for numerical eigenmodes, a calibration step is performed where a linearly polarised wave at frequency $\omega$ is launched with a constant $\boldsymbol {E}_\perp$, which lies along the diagonal of $y$–$z$ plane. In most cases, this wave has large overlaps with both numerical eigenmodes, so

(3.2a)\begin{gather} E_y = E_R\sin\tilde{\psi}_R \sin(\tilde{k}_Rx +\theta_R) +E_L\sin\tilde{\psi}_L \sin(\tilde{k}_Lx +\theta_L), \end{gather}

(3.2b)\begin{gather}E_z = E_R\cos\tilde{\psi}_R \cos(\tilde{k}_Rx +\theta_R) +E_L\cos\tilde{\psi}_L \cos(\tilde{k}_Lx +\theta_L), \end{gather}

have unknown amplitudes $E_R$ and $E_L$ and numerical wave vectors $\tilde {k}_R$ and $\tilde {k}_L$ are close to their analytical values $k_R$ and $k_L$. The phases $\theta _R$ and $\theta _L$ vary deterministically with $\omega t$, but the initial phases depend on the unknown expansion coefficients when spanning the wave in terms of the eigenmodes. Note that each electric-field component is of the form $E=\mathcal {E}_R(\sin \tilde {k}_Rx \cos b_R + \cos \tilde {k}_Rx \sin b_R) + \mathcal {E}_L(\sin \tilde {k}_Lx \cos b_L + \cos \tilde {k}_Lx \sin b_L)$, where $b=\theta$ for $E_y$ and $b=\theta +{\rm \pi} /2$ for $E_z$. By fitting a region of the simulation data where the wave envelopes are constant, as shown by the example in figure 5, the unknown amplitudes $\mathcal {E}_R$ and $\mathcal {E}_L$, as well as the unknown phases $\theta _R$ and $\theta _L$, can be determined by linearly regressing $\boldsymbol {E}_\perp (x)$ against a four-dimensional basis $(\sin \tilde {k}_Rx, \cos \tilde {k}_Rx, \sin \tilde {k}_Lx, \cos \tilde {k}_Lx)$. In practice, the values of $\tilde {k}_R$ and $\tilde {k}_L$ are not easily known, so they are determined numerically by minimising the residual of linear regression using $k_R$ and $k_L$ as initial guesses. Finally, using $\mathcal {E}_{y}=E\sin \tilde {\psi }$ and $\mathcal {E}_{z}=E\cos \tilde {\psi }$ from the best fit, the polarisation angles $\tilde {\psi }$ of numerical eigenmode is computed. In most cases, $\tilde {\psi }$ is close to the analytical $\psi$. However, outliers can occur in limiting cases, such as when $B_0\rightarrow 0$ where the degeneracy is week, or when $\theta _B\rightarrow 90^\circ$ where the ellipticity is weak. In outlier cases, it is usually sufficient to launch the desired eigenmode using $\psi$, which is calculated using the code in Shi (Reference Shi2022b).

To verify that only the intended eigenmode is launched and to measure its transmitted amplitude, pump-only and seed-only eigen-runs are performed, where one laser is launched with numerically determined polarisation angle $\tilde {\psi }$ and the other laser is turned off completely. Example transverse fields $E_y(x)$ and $E_z(x)$ are plotted in figure 5. When the wave is not an eigenmode (figure 5c), the polarisation trajectory precesses. On the other hand, when the wave is close to an eigenmode, its trajectory is near a stationary polarisation ellipse. In figure 5( f), a slow precession can still be observed, but the contamination from the other eigenmode is small enough that it is not a major source of error. Using the pump-only eigen-run, the transmitted pump amplitude can be measured from data shown in figure 5(d). The transmitted amplitude is normalised by analytical wave energy coefficient to compute $a_1$, which enters $\gamma _0$ via (2.4). Combining with the analytical $\varGamma$, the expected growth rate is computed to compare with simulations.

In the following, we focus on eigenmodes that are right-handed with respect to the magnetic field, which usually have the largest pump–seed coupling among the four possible polarisation combinations. When $\cos \theta _B>0$, meaning that $\boldsymbol {B}_0$ points towards $+x$, the seed laser is right-handed with respect to its wave vector. On the other hand, because the pump laser propagates in $-x$ direction, the eigenmode is left-handed with respect to its wave vector. Another way to describe the mode selection is that at a given wave frequency $\omega$, the eigenmode with smaller $k$ will be launched. This mode has an electric-field vector that corotates with the gyrating electrons and is therefore more strongly coupled with the electron-dominant plasma waves, such as the Langmuir-like P wave and electron-cyclotron-like F wave.

3.2. Extracting pump and seed envelopes

In stimulated runs where both pump and seed eigenmodes are turned on, we need to separate their contributions from the total electromagnetic fields, which are what PIC simulations observe directly. A plausible way of separating the waves is to use Fourier-like filters. However, because the pump and seed envelopes evolve in space–time due to three-wave interactions, their line shapes have long tails. Attempting to filter waves in the Fourier space and then taking inverse transform often introduce spurious oscillations. As an alternative approach, the pump and seed are directly separated in the configuration space using both electric- and magnetic-field data. Due to Faraday's law in one spatial dimension, wave fields are related in simple ways. For a right-propagating wave, $B_y=-E_z/v$ and $B_z=E_y/v$, where $v=\omega /k>0$ is the phase velocity of the wave. On the other hand, for a left-propagating wave, $B_y=E_z/v$ and $B_z=-E_y/v$ have the opposite signs. When the signal contains a single right $(+)$ and left $(-)$ waves at constant amplitudes, the electric fields can be decomposed as $E_y = E^+_{y} +E^-_{y}$ and $E_z = E^+_{z} + E^-_{z}$, and the magnetic fields are $B_y = -E^+_{z}/v_+ + E^-_{z}/v_-$ and $B_z = E^+_{y}/v_+ - E^-_{y}/v_-$. Note that the right and left waves likely have different phase velocities $v_+$ and $v_-$. After removing static components of $\boldsymbol {B}_0$, the wave fields satisfy

(3.3a)\begin{gather} E^+_{y} = \frac{E_y+v_- B_z}{1+v_-{/}v_+},\quad E^+_{z} = \frac{E_z-v_- B_y}{1+v_-{/}v_+}, \end{gather}

(3.3b)\begin{gather}E^-_{y} = \frac{E_y-v_+ B_z}{1+v_+{/}v_-},\quad E^-_{z} = \frac{E_z+v_+ B_y}{1+v_+{/}v_-}. \end{gather}

Reconstructions using the above Faraday filter are exact for constant-amplitude waves if $v_+$ and $v_-$ are known. However, due to numerical dispersion and PIC noise, $v_\pm$ differ slightly from analytical values, causing leakage through the Faraday filter.

To reduce the leakage, a calibration step is performed to determine the numerical phase velocities. Summing data of a pump-only run with data of a seed-only run before the wave fronts reach the opposite plasma boundary, $\tilde {v}_+$ and $\tilde {v}_-$ are fitted using (3.3). The data are selected within the plasma region using similar criteria as in figure 5. The best-fit residuals are typically orders of magnitude smaller than the signals but a few times above the PIC noise level, as shown by the example in figure 6. The coherent errors are likely due to how electric- and magnetic-field data are combined. Note that when solving Maxwell's equations using a Yee-like mesh, transverse electric fields are discretised on cell boundaries whereas transverse magnetic fields are discretised on cell centres. Therefore, interpolations are needed before $B_{i\pm 1/2}$ can be summed with $E_i$. Here, a simple linear interpolation is used to estimate $B_{i}\approx (B_{i-1/2} + B_{i+1/2})/2$, which is likely the leading cause of coherent errors. Using interpolation schemes that better match the PIC algorithm may further reduce the coherent error.

Figure 6. (a) The optimal phase velocities for separating the right ($+$) and left ($-$) propagating waves are fitted using calibration data (black), which is the sum of a pump-only ($E_1$) and a seed-only ($E_2$) runs before the laser reaches the opposite plasma boundary. The amplitudes of the pump (dashed blue) and seed (dashed red) are known, and the fitting minimises the total difference between known and extracted pump (cyan) and seed (magenta) amplitudes. (b) The best-fit residuals are orders of magnitude smaller than the signals. However, the residuals have coherent leakages that are a few times larger than the PIC noise level. Fitting results for other field components are similar. The intensities are $I_1=10^{15}\,{\rm W}\,{\rm cm}^{-2}$ and $I_2=5\times 10^{12}\,{\rm W}\,{\rm cm}^{-2}$.

However, other leakage sources are more important when the calibrated Faraday filter is applied to data from stimulated runs. Note that the simple relation between $\boldsymbol{E}$ and $\boldsymbol{B}$ assumes plane waves with constant amplitudes. When the amplitude varies, additional terms arise, which depends on derivatives of the amplitude. In regimes where the amplitude varies slowly, these additional terms are small compared with the leading term but are large compared with the error from interpolation. Moreover, due to reflection and spontaneous scattering in PIC simulations, waves other than the pump and seed are present in the system, leading to even larger errors. For example, the pump laser reflects from the left plasma–vacuum boundary and propagates towards right together with the seed laser. Typically, the reflected pump amplitude is ${\sim }0.5\,\%$ of the incident amplitude, which is ${\sim }7\,\%$ of the seed amplitude when $I_1/I_2=200$. The Faraday filter is far from ideal for separating co-propagating waves, especially when they have comparable phase velocities. The leakage leads to spurious oscillations of the seed envelope, which are nevertheless usually less severe than caused by Fourier filters.

Having separated out pump and seed electric fields, the next step is to extract the wave envelopes that enter the three-wave equations. For a wave with slowly varying envelope, its electric field is $E(x)=\mathcal {E}(x)\sin (kx+b)+\epsilon (x)$, where $\epsilon$ is an error field. The envelope $\mathcal {E}(x)$ can be estimated by mixing $E(x)$ with a sinusoidal reference, which gives $E(x) \sin (kx+b') = \frac {1}{2}\mathcal {E}(x)[\cos (b-b')-\cos (2kx+b+b')]+\epsilon (x)\sin (kx+b')$. Averaging the mixed signal over $\lambda =2{\rm \pi} /k$ ideally eliminates the last two terms to give

(3.4)\begin{equation} \langle E(x) \sin(kx+b')\rangle_\lambda=\tfrac{1}{2}\mathcal{E}(x)\cos(b-b'). \end{equation}

The reference phase $b'$ is scanned to maximise the norm of the average, which is attained when the phases are matched $b'=b$ up to integer multiples of ${\rm \pi}$. Compared with a simple moving average, the above lock-in scheme reduces sensitivities to PIC noise and coherent errors. The numerical wave vector $\tilde {k}$, which is determined when fitting (3.2), is used to generate the sinusoidal reference. Moreover, to better remove the unwanted terms using the $\lambda$ average, a matching grid $\tilde {x}_i$ with spacing $\lambda /n$ is used, where the integer $n=\lceil \lambda /\Delta x\rceil$ and $\Delta x$ is the original spacing of the PIC grid $x_i$. The reference is generated on grid $\tilde {x}_i$ and mixed with $E(\tilde {x}_i)$, which is estimated from $E(x_i)$ using linearly interpolation. Averaging over $\lambda$ down samples the data, so the envelope $\mathcal {E}(x)$ lives on a grid that is more sparse than the $\tilde {x}_i$ grid by a factor of $n$. Note that due to their wavelengths difference, the pump and seed envelopes live on two separate grids. In addition to the transverse components, the $\mathcal {E}_x$ component is estimated using the lock-in scheme alone, without attempting to separate right and left contributions. The longitudinal component is small but non-zero when lasers propagate obliquely with $\boldsymbol {B}_0$. Combining all components, the total scalar wave envelope $\mathcal {E}=|\boldsymbol{\mathcal{E}}|$ is determined. Typical results of envelope extraction are shown in figure 7 for the $E_y$ component. Results are similar for the $E_z$ component, but are more noisy for the much smaller $E_x$ component. The procedure successfully separates out the pump and seed lasers, and identifies their envelopes that enclose fast-oscillating fields.

Figure 7. Raw data (a–c) is separated using a Faraday filter into pump (cyan) and seed (magenta) contributions (d–f), and envelopes of pump (blue) and seed (red) are extracted using a lock-in scheme. At $t\approx 0.25$ ps (a,d), the pump almost reaches the left plasma boundary and the seed has not entered. At later time $t\approx 0.48$ ps (b,e), the pump fills the plasma and the seed enters from the left. At final time $t\approx 0.67$ ps (c, f), the seed exits the right boundary and experiences noticeable growth, and the pump is slightly depleted. In this example, $T_0=10$ eV, $B_0=3$ kT, $\theta _B=30^\circ$, $I_1=10^{15}\,{\rm W}\,{\rm cm}^{-2}$ and $I_2=5\times 10^{12}\,{\rm W}\,{\rm cm}^{-2}$. The interaction is mediated by a P wave and the resonant seed wavelength in vacuum is $\lambda _2\approx 1.09\,\mathrm {\mu }$m.

3.3. Fitting data to analytical solutions

Before fitting extracted seed envelopes to analytical solutions, we need to ensure that the seed frequency is on resonance. Note that (2.1) is for resonant interactions that satisfy energy–momentum conservation $k^\mu _1=k^\mu _2+k^\mu _3$ where $k^\mu$ is the wave 4-momentum. When phase-matching conditions are not satisfied, the waves can still interact, but a detuning term $\exp ({\rm i}x^\mu \delta k_\mu )$ will appear in the three-wave equations. Analytical solutions in the detuned case are different from what is discussed in § 2, and usually exhibit less seed growth as a result of phase modulations. To determine the resonant $\omega _2^*$, the seed frequency is scanned near analytically expected resonances to maximise the growth of seed envelope $\mathcal {E}_2$ in kinetic simulations. To reduce the requisite number of runs, only a few frequencies are simulated and $\hat {\mathcal {E}}_2(\omega _2)\propto [(\omega _2-\omega _2^*)^2+\mu ^2]^{-1}$ is fitted using a Lorentzian profile. Here, $\hat {\mathcal {E}}_2$ is the spatial maximum of $\mathcal {E}_2$ at a given time slice and $\mu$ is a phenomenological linewidth. In most cases, the fitting behaves as expected, as shown in figure 8(a). Using $\hat {\mathcal {E}}_2$ at different time slices yields consistent estimates of $\omega _2^*$ (dashed vertical lines), as long as $\mathcal {E}_2$ has grown substantially above the noise level. In most cases, only nine $\omega _2$ values are scanned within a $3\times 10^{13}$-rad s$^{-1}$ window near the expected resonance. The fitting uses last four simulation outputs, when the seed has propagated across more than half the plasma length. The four estimates of $\omega _2^*$ are averaged using equal weight to give a final estimate (solid vertical line). However, in some cases, more complicated spectra emerge during the scan, as shown in figure 8(b). In this example, $B_0=200$ T and the electron gyro frequency is a few times smaller than the plasma frequency. In this regime, the kinetic wave dispersion relation involves multiple Bernstein-like waves near the expected warm-fluid resonance. In cases such as this, additional $\omega _2$ values are scanned to resolve individual peaks, and only data near the expected resonance (solid circles) is used to fit $\omega _2^*$. Using data at multiple time slices again yields consistent estimates. As an intriguing observation, spacings between the side peaks do not seem to have a simple dependence on $B_0$ and $\theta _B$, and often appears to change with time. These features remain to be investigated in the future.

Figure 8. The spatial maximum of seed envelope $\hat {\mathcal {E}}_2$ is fitted to seed frequency $\omega _2$ using a Lorentzian profile to determine the resonant seed frequency. (a) When $B_0=3$ kT and the interaction is mediated by a P wave, the data matches a simple Lorentzian, and fittings using $\hat {\mathcal {E}}_2$ at different time slices yield consistent results (dashed vertical lines), which is averaged to give a final estimate of $\omega _2^*$ (solid vertical line). This is representative of what happens in most cases. (b) When $B_0=200$ T and the interaction is mediated by an F wave, the data reveal multiple weak resonances. Data away from the expected warm-fluid peak are excluded from fitting (empty circles). The included data (solid circles) again yield consistent results across multiple time slices. In these examples, $T_0=10$ eV, $\theta _B=30^\circ$, $I_1=10^{15}\,{\rm W}\,{\rm cm}^{-2}$ and $I_2=5\times 10^{12}\,{\rm W}\,{\rm cm}^{-2}$.

Using a resonant seed to set up stimulated runs, the extracted seed envelopes are placed into analytic frame to prepare for fitting, as shown in figure 9. Note that in the step-function problem, the plasma–vacuum boundary is sharply at $x=0$ and a constant-amplitude seed arrives at the boundary when $t=0$. These idealised set-ups are necessarily softened in simulations, where the plasma expands into the vacuum forming a sheath region, and the seed is ramped up using a $\tanh$ profile whose temporal width equals to the seed period to smooth out numerical artefacts. For the spatial axis, to avoid boundary effects, data within $4\lambda _1\gg \lambda _D$ from the initial plasma–vacuum boundaries are excluded, where $\lambda _D$ is the Debye length. Hence, the data frame $x_d$ is shifted from the analytic frame $x_a=x_d+x_0$ by an offset $x_0$. For the time axis, because the seed is launched after the pump has propagated across the plasma, the data frame $t_d$ is delayed from the analytic frame $t_a=t_d-t_0$ for some offset $t_0$. To place data into the analytic frame, the seed wave front $x_a=\tilde {v}_2t_a$ is fitted to determine both offsets and the numerical group velocity $\tilde {v}_2$, which is close to but not equal to the analytical group velocity $v_2=\partial \omega /\partial k_x$ due to numerical dispersions. The seed wave front is identified as the location where the seed envelope $\mathcal {E}_2$ drops below half its boundary value, and the thickness of the wave front is identified as the spatial separation where $\mathcal {E}_2$ attains $10\,\%$ and $90\,\%$ of its boundary value. Only data behind the wave front by twice its thickness are used for fitting. Finally, we need to normalise the electric field envelope $\mathcal {E}_2$. Note that the analytical solution of the step-function problem can be written in terms of the ratio $\alpha _2/h_0$, where $h_0$ is the boundary value. In simulations, the value at the domain boundary is an input, but the value at the plasma–vacuum boundary is not controlled, because only a fraction of the incident seed laser is transmitted into the plasma. The boundary value inside the plasma is determined using the seed-only eigen-run. Using calibration data similar to what is shown in figure 5(d), the transmitted seed amplitude is measured and used to normalise envelopes in the stimulated run. Note that the calibration needs to be performed each time when plasma conditions are changed or the seed laser is varied.

Figure 9. Normalised seed amplitude $\alpha _2/h_0$ in the analytical frame $(x_a, t_a)$. The step-like feature near the wave front is a visualisation artefact when only nine time slices are plotted. The wave front $x_a=\tilde {v}_2 t_a$ (white line) is fitted to determine the numerical group velocity $\tilde {v}_2$. This example is the same run as in figure 7.

Having extracted $\alpha _2/h_0$ from the stimulated run, the data in the analytic frame are fitted to (2.27). Here, the simplification $v_3=0$ is made because $|v_3|\ll v_2$. From figure 2, we see that $v_3=0$ is already a good approximation even when $v_3=-v_2/10$. For a warm-fluid plasma with moderate temperature $T_0\sim 10$ eV, the ratio $v_3/v_2\sim O(10^{-2})$ is much smaller, so the limiting-form growth function provides an even better approximation. Evaluating the growth function requires a single numerical integral, which is much cheaper than computing the double integral in (2.25) that gives the exact solution at finite $v_3<0$. As another simplification, the physical collision module is turned off in PIC simulations, so lasers are undamped beyond spurious numerical collisions. Setting $\mu _2=0$ gives $\Delta _2=\alpha _2/h_0-1$. The right-hand side of this expression is determined from simulation data to give $\tilde {\Delta }_2$, and the left-hand side is the analytical formula, which contains three parameters $v_2$, $\mu _3$ and $\gamma _0$. As $\tilde {v}_2$ is already determined from fitting the wave front, only two parameters remain to be fitted. The fitting is treated as an optimisation problem, where parameters are searched near their expected values to minimise the residual $\|\tilde {\Delta }_2-\Delta _2\|$, where $\|f\| = ({1}/{n_xn_t}) [\sum _{i=1}^{n_x}\sum _{j=1}^{n_t}|f(x_i,t_j)|^2]^{1/2}$. The discrepancy $\tilde {\Delta }_2-\Delta _2$ is set to zero outside the region where data are valid. As discussed earlier, the data are valid when $x$ is sufficiently far away from both the plasma boundaries and the seed front, and when $t$ is after seed experiences noticeable growth but $a_1(x)$ and $f_e(v_x)$ remains largely unchanged. Typical best-fit results are shown in figure 10, where the simulation data (circles) are well matched by analytical solutions (lines), and the residuals are dominated by leakages during pump–seed separation.

Figure 10. (a) The growth function extracted from simulation data (circles) is well matched by the analytical solutions (lines) using best-fit parameters. Results are shown for three representative time slices as the seed propagates and grows. The fitting quality at other time slices are similar. (b) Residuals of best fit are an order of magnitude smaller than the signal, and are dominated by leakage during pump–seed separation. The reflected pump amplitude is ${\sim }7\,\%$ of the seed. In this example, $T_0=100$ eV, $B_0=3$ kT, $\theta _B=30^\circ$, $I_1=10^{15}\,{\rm W}\,{\rm cm}^{-2}$, $I_2=5\times 10^{12}\,{\rm W}\,{\rm cm}^{-2}$ and the interaction is mediated by a P wave. The best-fit parameters are $\tilde {v}_2\approx 2.97\times 10^8$ m s$^{-1}$, $\tilde {\mu }_3\approx 5.17\times 10^{12}$ rad s$^{-1}$ and $\tilde {\gamma }_0\approx 9.25\times 10^{12}$ rad s$^{-1}$.

Figure 11. Uncertainties of fitting parameters are quantified by scanning the fit residual. (a,c) When $T_0=10$ eV is cold, damping is weak and poorly constrained. The best fit is far below the vertical scale and is consistent with a negligible $\mu _3$. Nevertheless, $\gamma _0$ is well constrained and matches analytical growth rate (vertical solid line). (b,d) When $T_0=100$ eV is hotter, damping becomes stronger. In this case, both $\gamma _0$ and $\mu _3$ are well constrained. The best fit (blue dot) matches the analytical growth rate within error bar, which is defined as the range where fit residual doubles its minimum. (c,d) Using the best fit $\mu _3$, a one-dimensional scan along $\gamma _0$ gives marginal fit residuals and error bars (horizontal dashed line). In these examples, $B_0=3$ kT, $\theta _B=30^\circ$, $I_1=10^{15}\,{\rm W}\,{\rm cm}^{-2}$, $I_2=5\times 10^{12}\,{\rm W}\,{\rm cm}^{-2}$ and interactions are mediated by P waves. (b,d) are the same run as in figure 10.

To quantify uncertainty, the fit residual is scanned over a range of parameter values to estimate error bars. As parameters move away from their best-fit values, the fit residual increases. Error bars are defined as the parameter range beyond which the residual doubles its minimum. Using this definition, in the limit where data exactly match analytical solutions, the residual is zero and the error bar is also zero. In the opposite limit where simulation and theory poorly match, the residual is large and the error bar is wide. Two examples are shown in figure 11, where the resonant interactions are mediate by Langmuir-like P waves. When the temperature $T_0=10$ eV is cold (figure 11a), collisionless damping of the plasma wave is minuscule. In this case, the residual remains small for a wide range of $\mu _3$ values, and the best fit is $\tilde {\mu }_3\approx 1.36\times 10^{-2}$ rad s$^{-1}$, which is consistent with zero damping. On the other hand, the residual increases rapidly when $\gamma _0$ deviates from its best-fit value $\tilde {\gamma }_0\approx 9.22\times 10^{12}$ rad s$^{-1}$, which is close to the analytical value $\gamma _0^a\approx 9.53\times 10^{12}$ rad s$^{-1}$ (solid vertical line). The marginal error bar is shown in figure 11(c), where $\mu _3$ is fixed at its best-fit value and $\gamma _0$ is scanned. The residual is convex near $\tilde {\gamma }_0$ (dashed vertical line) and $\gamma _0^a$ is within the error bar (dashed horizontal line). When $T_0=100$ eV is hotter (figure 11b), collisionless damping is stronger, so both $\gamma _0$ and $\mu _3$ are constrained. As discussed earlier, magnetised collisionless damping is not easy to compute, so the procedure here provides a unique method for measuring the damping rate, which in this case is best fitted by $\tilde {\mu }_3\approx 5.17\times 10^{12}$ rad s$^{-1}$. The best-fit $\tilde {\gamma }_0\approx 9.25\times 10^{12}$ rad s$^{-1}$ matches the analytical value $\gamma _0^a\approx 9.05\times 10^{12}$ rad s$^{-1}$ within the marginal error bar shown in figure 11(d). The marginal residual is always convex and usually skewed towards smaller $\gamma _0$. The two-dimensional landscape of the residual is flat along a valley in the $\gamma _0$–$\mu _3$ space. This is intuitive because a larger damping cancels the effect of a larger growth. In the asymptotic limit $t\rightarrow +\infty$, the exponential solution ${\rm e}^{\kappa x}$ only depends on the ratio $\kappa =\gamma _0^2/\mu _3v_2$, so $\gamma _0$ and $\mu _3$ are not independently constrained. In other words, the separate constrains come from solutions at earlier time. As discussed after (2.27), the growth function $\Delta _2$ is symmetric about $x=\frac {1}{2}v_2t$ in the absence of damping. As $\mu _3$ increases, the growth behind $x=\frac {1}{2}v_2t$ is reduced more, leading to an asymmetric profile of $\Delta _2$. The asymmetry in the transient-time profile is what allows the effects of $\mu _3$ and $\gamma _0$ to be distinguished.

3.4. Scanning physical parameters

Having described in detail how simulations are set up and how data are processed, the protocol is applied to scan physical parameters and compare simulations with theory. First, it is useful to verify that the pump and seed intensities are selected properly. As discussed at the beginning of this section, the pump and seed intensities are constrained by ${10^{10}\,{\rm W}\,{\rm cm}^{-2}} \ll I_2\sim 10^{-2} I_1\ll {10^{14}\,{\rm W}\,{\rm cm}^{-2}}$. Within this range, the bare growth rate $\gamma _0\propto a_1\propto I_1^{1/2}$ scales with the pump intensity $I_1$ but is independent of $I_2$. This theory prediction is confirmed by results shown in figure 12, where black lines are from the theory and orange symbols with error bars are simulation results. In figure 12(a), the seed intensity is fixed at $I_2=5\times 10^{12}\,{\rm W}\,{\rm cm}^{-2}$ and the agreement is excellent within the range $I_1$ is scanned. Note that the black line is not a fit of data. A hint of deviation can be seen from the last data point at $I_1=8\times 10^{15}\,{\rm W}\,{\rm cm}^{-2}$, where data start to drop below the theory line. This behaviour is expected because a pump that is too intense causes large spontaneous scattering that depletes the pump. Signatures of spontaneous scattering are clearly visible in the inset, which shows the Fourier power spectrum of transverse electric fields within the plasma region. At $I_1=5\times 10^{14}\,{\rm W}\,{\rm cm}^{-2}$(blue), the spectrum is clean and dominated by the pump and seed near $ck=2\times 10^{15}$ rad s$^{-1}$ and the plasma wave near $ck=4\times 10^{15}$ rad s$^{-1}$. In comparison, at $I_1=8\times 10^{15}\,{\rm W}\,{\rm cm}^{-2}$(red), extra waves are excited due to spontaneous scattering. These extra waves remove energy that would otherwise be used to grow the seed, and they also constitute a messy background that jeopardises the data analysis. In figure 12(b), the agreement remains excellent within the range $I_2$ is scanned, where the pump intensity is fixed at $I_1=10^{15}\,{\rm W}\,{\rm cm}^{-2}$. The error bars remain largely constant due to two competing effects. At stronger seed intensities, pump depletion occurs earlier so more data are excluded from fitting, resulting in fewer constraints and larger error bars. On the other hand, at weaker seed intensities, the leakage during pump–seed separation more severely distorts the seed profile, causing poorer fits and larger error bars. The two competing effects cause error bars to remain roughly constant within this $I_2$ range. Effects of leakage can be clearly seen in the inset. The normalised envelope $\alpha _2/h_0$ is relatively smooth when $I_2=10^{13}\,{\rm W}\,{\rm cm}^{-2}$ (magenta), but the spurious oscillations become more severe when $I_2=6.25\times 10^{11}\,{\rm W}\,{\rm cm}^{-2}$ (cyan). The artefacts are primarily due to reflected pump that copropagates with the seed, which cannot be separated using (3.3) and (3.4), causing spurious oscillations. Limited by these effects, simulations below use $I_1=10^{15}\,{\rm W}\,{\rm cm}^{-2}$ and $I_2=5\times 10^{12}\,{\rm W}\,{\rm cm}^{-2}$.

Figure 12. The bare growth rate $\gamma _0\propto I_1^{1/2}$ scales with pump intensity (a) and is independent of the seed intensity (b) within the range where simulations are described by the linearised three-wave equations. Theory predictions (black lines) match simulation results (orange symbols) within error bars. At large pump intensity $I_1$, spontaneous scattering depletes the pump and excites additional waves, as shown in the inset of (a) where the Fourier power spectra are plotted for transverse electric fields within the plasma region. At larger seed intensity $I_2$, pump depletion occurs earlier, which reduces the amount of data that can be used for fitting analytical solutions. On the other hand, at smaller $I_2$, the leakage during pump–seed separation is more detrimental, as shown in the inset of (b) where the normalised envelopes $\alpha _2/h_0$ suffer from spurious oscillations. In this set of scans, $T_0=10$ eV, $B_0=3$ kT, $\theta _B=30^\circ$, the vacuum seed wavelength is $\lambda _2\approx 1.088\,\mathrm {\mu }$m and the interaction is mediated by a P wave.

Second, the protocol is used to scan plasma temperature $T_0$ to determine the range of validity of warm-fluid theory (figure 13). At low temperature, collisionless damping of plasma waves are negligible and the $\mu _3$ error bars are large. On the other hand, when $T_0>10^2$ eV, damping starts to dominate growth and data start to deviate from warm-fluid theory. At $T_0\sim 10^3$ eV, the seed barely grows unless the pump is stronger. However, a stronger pump suffers more from spontaneous scattering, leading to systematic errors and large uncertainties. Therefore, the applicability of the protocol does not far exceed the temperature range reported here. Within this range, the expected $\gamma _0$ from warm-fluid theory (line) agrees with simulation data (squares) up to $T_0=300$ eV, where hints of disagreements emerge. It is remarkable that even when collisionless damping has become significant, which is a purely kinetic phenomenon, the three-wave coupling remains well described by the warm-fluid theory. In subsequent scans, the plasma temperature is fixed at $T_0=10$ eV where damping is negligible.

Figure 13. The growth rate $\gamma _0$ (left axis) and damping rate $\mu _3$ (right axis) as functions of the plasma temperature $T_0$. Coloured symbols with error bars are simulation results and the orange line is the expected growth rate from warm-fluid theory. When $T_0$ increases by tenfold, $\mu _3$ increases by orders of magnitude. On the other hand, $\gamma _0$ decreases only slightly, and data are well matched by theory until the last two points. In this set of scans, $B_0=3$ kT, $\theta _B=30^\circ$, $I_1=10^{15}\,{\rm W}\,{\rm cm}^{-2}$, $I_2=5\times 10^{12}\,{\rm W}\,{\rm cm}^{-2}$ and the interactions are mediated by P waves.

Figure 14. (a) Growth rate $\gamma _0$ and (b) resonant seed frequency $\omega _2$ as functions of the background magnetic field $B_0$ when interactions are mediated by electron-dominant P and F waves. Analytical results (lines) match simulation data (symbols). Error bars are shown for $\gamma _0$ but are too small to show for $\omega _2$. The index $b_3$ labels the wave branch of $\alpha _3$. For example, $b_3=4$ means $\alpha _3$ is the fourth highest-frequency branch of the warm-fluid dispersion relation. In this set of scans, $T_0 =10\ {\rm eV}$, $\theta _B=30^\circ$, $I_1=10^{15}\,{\rm W}\,{\rm cm}^{-2}$ and $I_2=5\times 10^{12}\,{\rm W}\,{\rm cm}^{-2}$.

Figure 15. Same as figure 14, but for $\theta _B=\langle \boldsymbol {k},\boldsymbol {B}_0\rangle$ scan at fixed $B_0=3$ kT.

Third, the protocol is used to scan the background magnetic field strength $B_0$ for interactions that are mediated by electron-dominant resonances (figure 14). In weak magnetic fields, the $b_3=3$ resonance, namely, the $\alpha _3$ eigenmode whose frequency is the third highest at a given wave vector, is the Langmuir-like P wave. In particular, at $B_0=0$, this resonance is precisely mediated by the Langmuir wave that gives rise to Raman backscattering in unmagnetised plasmas. When $B_0\ne 0$ and $\theta _B\ne 0$, this eigenmode is modified due to cyclotron motion, but its frequency only weakly dependents on $B_0$, as shown in figure 14(b). Beyond $B_0\approx 10$ MG, the eigenmodes cross over, and the $b_3=3$ resonance becomes the electron-cyclotron-like F wave, whose frequency $\omega _3$ increases almost linearly with $B_0$. Complementarily, the $b_3=4$ resonance is the F wave in weak $B_0$ and crosses over to become the P wave in strong magnetic fields. Within the $B_0$ range scanned here, the interaction mediated by the F wave provides a smaller coupling than the P wave at the same $B_0$. In particular, at $B_0=0$, the F wave vanishes and provides zero coupling. The exception is near $B_0=10$ MG when the two eigenmodes strongly hybridise. Although they remain two distinct eigenmodes, their couplings during frequency crossover become equal. Immediately after the crossover, P-wave coupling recovers its unmagnetised strength and continues to increase, whereas F-wave coupling drops before starting to increase again. In stronger fields not shown here, $\gamma _0$ continues to increase and F-wave coupling surpasses P-wave coupling (Shi & Fisch Reference Shi and Fisch2019). In even stronger fields, the F-wave frequency approaches $\omega _1/2$ and two-magnon decays are encountered (Manzo et al. Reference Manzo, Edwards and Shi2022), for which kinetic effects are dominant. After the F-wave frequency increases beyond $\omega _1/2$, it can no longer mediate resonant interactions, so only the P resonance remains. As $B_0$ further increases, mediation by P wave continues to strengthen, and a strong resonance is encountered when electron gyro frequency approaches $\omega _1$ (Edwards et al. Reference Edwards, Shi, Mikhailova and Fisch2019). Beyond that point, the 1-$\mathrm {\mu }$m pump frequency is below the R-wave cutoff, so only L-wave pump can propagate, whose coupling is substantially weaker. In the range scanned in figure 14, kinetic effects are not overwhelming and the agreement is excellent. Note that the lines are not fits of data. Although error bars for $\gamma _0$ can perhaps be reduced, growth rates extracted from simulations already match detailed features of analytical predictions.

Finally, the protocol is used to scan $\theta _B$, the angle between the magnetic field and the wave vectors, when interactions are mediated by the two highest-frequency resonances (figure 15). The magnetic field is fixed at $B_0=3$ kT, where the electron cyclotron frequency is larger than the plasma frequency. In this case, the $b_3=3$ resonance is the electron-cyclotron wave when $\theta _B=0$ and becomes the UH wave when $\theta _B=90^\circ$. The other resonance with branch index $b_3=4$ is the Langmuir wave when $\theta _B=0$ and becomes the LH wave when $\theta _B=90^\circ$. The coupling provided by the electron-cyclotron wave is zero, because the interaction is both polarisation forbidden and energy forbidden (Shi et al. Reference Shi, Qin and Fisch2018). The coupling provided by the LH wave is finite, but smaller than UH mediation by a factor that is roughly proportional to $(\omega _{\mathrm {LH}}/\omega _{\mathrm {UH}})^{1/2}\ll 1$ (Shi et al. Reference Shi, Qin and Fisch2017b). On the other hand, the coupling provided by the electron-dominant UH wave is strong, and has a comparable strength as unmagnetised Raman scattering. As $\theta _B$ increases, the hybridisation between cyclotron motion and electrostatic oscillation changes. The growth mediated by the $b_3=3$ resonance increases, whereas the growth mediated by the $b_3=4$ resonance decreases. The simulation data (symbols) and analytical results (lines) agree within error bars. Note that the lines are not fits of data. The data point at $b_3=4$ and $\theta _B=90^\circ$ is missing because the LH coupling is too weak and the resonance cannot be unambiguously identified.