2 Results of hybrid simulations

We conducted simulations using a 2.5-dimensional hybrid code in which the energetic electrons are treated as particles, while the background electrons are treated as a fluid (Taguchi, Antonsen Jr & Mima Reference Taguchi, Antonsen and Mima2004). The ions are assumed to be immobile. The basic equations of our hybrid code are the same as described in Taguchi et al. (Reference Taguchi, Antonsen and Mima2004) with the exception that the current code uses the full set of Maxwell’s equations including the displacement current, whereas the Darwin approximation was made in the reference.

Figure 1. (a) Line plots of field quantities averaged over the transverse direction in the case of no external magnetic field. (b) The spatial distribution of the beam electron density at $\unicode[STIX]{x1D714}_{pe}t=1600$ . (c) Contour plots of electron distribution function in the longitudinal phase space ( $z$ – $p_{z}$ ).

Figure 2. (a) The line plots of field quantities averaged over the transverse direction in the case that an external magnetic field ( $\unicode[STIX]{x1D714}_{c}/\unicode[STIX]{x1D714}_{pe}=0.3$ ) is applied. (b) The spatial distribution of the beam electron density at $\unicode[STIX]{x1D714}_{pe}t=1600$ , (c) Contour plot of the electron distribution function in the longitudinal phase space ( $z$ – $p_{z}$ ).

Results of sample simulations are shown in figures 1(a–c) and 2(a–c). The simulation domain is $0<z/d_{e}<750$ in the direction of beam propagation ( $z$ ) and $0<x/d_{e}<125$ in one direction transverse to the direction of beam propagation ( $x$ ), where $d_{e}$ is the collisionless skin depth $d_{e}=c/\unicode[STIX]{x1D714}_{pe}$ and the plasma frequency $\unicode[STIX]{x1D714}_{pe}=\sqrt{e^{2}n_{0}/\unicode[STIX]{x1D700}_{0}m_{e}}$ is based on the background electron density $n_{0}$ and the rest mass of an electron $m_{e}$ . The boundary conditions in $x$ are that all quantities are periodic.

The boundary conditions in $z$ are more complicated. Two damping layers are added: one to each end of the simulation box. In the layers, all field quantities, such as the electromagnetic field components and the fluid variables describing the background plasma, are gradually relaxed to ‘equilibrium values’ by adding a damping term in the time evolution equation for each variable of the form:

Here, $f$ is one of the quantities and $f_{0}$ is represents its equilibrium value. The quantity $D_{damp}$ is the width of the damping layer and the maximum damping rate, $\unicode[STIX]{x1D708}_{max}$ is empirically determined for each value so that each quantity smoothly approaches its equilibrium value at both boundaries. The coupling strength between the beam electrons and the electromagnetic fields is also gradually reduced in the damping layer.

An energetic beam of electrons is introduced in the layer near $z=0$ , the boundary on the left, with a distribution in momentum given by, $f(\boldsymbol{p})=A\exp [(m_{e}\unicode[STIX]{x1D6FE}c^{2}-\boldsymbol{p}\boldsymbol{\cdot }\boldsymbol{u}_{0})/T_{b}]$ . Here, $\boldsymbol{p}$ is the momentum of each particle, $\unicode[STIX]{x1D6FE}=\sqrt{1+\boldsymbol{ p}^{2}/m_{e}^{2}c^{2}}$ is its Lorentz factor, $\boldsymbol{u}_{0}$ is the drift velocity and $T_{b}$ is the temperature of the beam electrons. During each time step of the simulation particles sampled from this distribution are placed in the first cell.

Parameters in the simulation are set as follows: beam velocity, $u_{0}=0.95c$ , injected beam density, $n_{b}=n_{0}/10$ , beam temperature, $T_{b}=100$  keV, while the background electron temperature $T_{e}$ is 10 keV. In a plasma realized in a fast ignition scheme, the density of the background plasma varies in the range from 1 to 100 times the critical density for the incident laser pulse. Thus, the density ratio we chose corresponds to a background electron density of $n_{0}=10^{22}~\text{cm}^{-3}$ for a $1\unicode[STIX]{x03BC}$ -laser. For these parameters the background electron collision time is estimated to be approximately $3\times 10^{-11}/Z$  (s), where $Z$ is an average charge state of ions, and the plasma frequency is $5.6\times 10^{15}~\text{s}^{-1}$ . According to these estimates the collision time is longer than the duration of the simulation ( $1600/\unicode[STIX]{x1D714}_{pe}=3\times 10^{-13}$  s) even for a gold cone and as a result the background is treated as being collisionless. Of course this estimate is influenced by the assumed value of the background temperature. Collisions will become important if the temperature is significantly lowered. The processes we observe and will describe occur on the time scale of a hundred plasma periods. The temperature for which this equals the collision time is approximately 10 eV. Studying the effects of collisions is important, but will be left for the future.

Figure 1(b) displays the beam density (colour scale indicates density relative to the background density) at $\unicode[STIX]{x1D714}_{pe}t=1600$ with no applied magnetic field. This is well after the front of the beam has passed through the simulation domain and a statistically steady state has formed. The injected beam distribution is Weibel unstable and the beam forms filaments and heats in the transverse direction. This is also illustrated by the line plots attached in figure 1. These are line plots of the following $x$ -averaged quantities: the normalized beam density ( $\bar{n}_{h}$ ), the normalized in-plane ( $\bar{B}_{x}$ ) and out-of-plane ( $\bar{B}_{y}$ ) magnetic field fluctuations and the longitudinal electric field fluctuations ( $\bar{E}_{z}$ ). Here the magnetic fields are normalized by $m_{e}\unicode[STIX]{x1D714}_{pe}/e$ and the electric field is normalized by $m_{e}\unicode[STIX]{x1D714}_{pe}c/e$ . We also display the $z{-}p_{z}$ phase space distribution of the beam electrons as a contour plot in figure 1(c). As shown in the figure, after the beam passes through the primary unstable region of the Weibel instability ( $50<z/d_{e}<120$ ), it heats and decelerates. Since the Weibel instability becomes very weak for high transverse beam temperature, the magnetic field fluctuations in the downstream region are relatively small compared with those in the upstream region. While there is the potential for a residual Weibel instability in the downstream region involving the cold electrons (Fiuza et al. Reference Fiuza, Fonseca, Tonge, Mori and Silva2012), it does not lead to large magnetic fields in this region. As a result, the heated beam crosses the entire simulation domain as evidenced by the line plot of the beam density. At the time of this image the ratio of beam density leaving the right boundary to that injected is 1.25. This small density increase is caused by the decrease of the flow velocity due to beam electron scattering by the Weibel instability. This is confirmed by the fact that the ratio of the fluxes of positively propagating beam electrons at the two ends of the simulation is close to unity, $R_{flux}(=nu_{h+}(z/d_{e}=400)/nu_{h+}(z/d_{e}=0))$ , is 0.988.

The situation is different when a strong magnetic field is applied along the direction of beam propagation. This is illustrated in figure 2 for the case of the normalized external $B$ -field, $\unicode[STIX]{x1D714}_{c}/\unicode[STIX]{x1D714}_{pe}=0.3$ , where $\unicode[STIX]{x1D714}_{c}=eB_{0}/m_{e}$ for the external magnetic field $B_{0}$ . This value corresponds to 9.6 kT in the case of $n_{0}=10^{22}~\text{cm}^{-3}$ . In the simulation, the Weibel instability grows and heats the beam, as evidenced by the out-of-plane magnetic field fluctuations. However, after the beam propagates to a distance $z/d_{e}\sim 200$ , both in-plane and out-of-plane magnetic fields grow and the beam electrons are reflected and their density accumulates as shown in figure 2(c). This is further illustrated by the line plot of the beam density, where the downstream density is half that of figure 1. For this case, the ratio of the electron fluxes entering on the left and leaving on the right ( $R_{flux}$ ) at the time of figure 2 is 0.292. This value continuously decreases in time as the simulation progresses.

The in-plane and out-of-plane components of the magnetic field for $180<z/d_{e}<300$ in figure 2 have the characteristics of a circularly polarized wave: they vary sinusoidally in $z$ with wavelength approximately $\unicode[STIX]{x1D706}\approx 16d_{e}$ , and they are $90^{\circ }$ out of phase. We thus identify the disturbance as a whistler wave. The dispersion relation for long-wavelength ( $\unicode[STIX]{x1D706}\gg d_{e}$ ) low-frequency ( $\unicode[STIX]{x1D714}\ll \unicode[STIX]{x1D714}_{c}$ ) whistler waves propagating at an angle to the magnetic field in cold plasma is (Kennel Reference Kennel1966)

where $k=\sqrt{k_{x}^{2}+k_{z}^{2}}$ . Here $\unicode[STIX]{x1D714}$ is the wave frequency and $\boldsymbol{k}=k_{x}\hat{\boldsymbol{x}}+k_{z}\hat{\boldsymbol{z}}$ is the wave vector. For parallel propagation ( $k_{x}=0$ ), as applies to the $x$ -averaged fields, the dispersion relation implies a positive phase velocity and a positive group velocity along the magnetic field. The group velocity increases with magnetic field strength, the implications of which will be discussed subsequently.

3 Mechanism of an excitation of a large amplitude whistler wave

Direct excitation of the purely parallel-propagating whistler by wave particle energy exchange involving the energetic stream of electrons is not expected. Wave–particle interaction is possible when the resonance condition

is satisfied, where the integer $n=1,0,-1$ denotes cyclotron, Cherenkov and anomalous cyclotron resonances respectively (Kennel Reference Kennel1966). These resonances are mediated by different components of the electromagnetic field. The Cherenkov resonance is mediated by the axial electric field (which is small for whistler waves), and the cyclotron and anomalous cyclotron resonances are mediated by one or the other circular polarizations of the transverse electric field. The whistler wave is circularly polarized in the same sense of rotation as is the gyro-motion of electrons, which means that only the $n=1$ resonance is active for purely parallel propagation. Thus, for the low-frequency whistler wave ( $\unicode[STIX]{x1D714}\ll \unicode[STIX]{x1D714}_{c}$ ) only counter-propagating electrons can be resonant, and since the wave phase velocity is positive, energy will flow from particles to fields only if the perpendicular temperature of these counter-propagating electrons is greater than parallel temperature. This is the mechanism of the heat flux instability studied extensively by Gary (Gary et al. Reference Gary, Feldman, Forslund and Montgormery1975; Gary Reference Gary1985). To tap the energy of the forward-propagating hot electrons it is necessary to consider off-angle propagation, for which the transverse wave electric field becomes elliptically polarized and non-zero Larmor radius effects enter. Both of these effects activate the $n=-1$ resonance and energy can be transferred from electrons to fields as the electrons lower their energy while increasing their perpendicular momentum in the presence of a wave with a positive phase velocity (Kennel Reference Kennel1966).

In our simulations the mechanism of excitation of the parallel-propagating whistler wave, once the Weibel has stabilized, is by the growth of hot electrons driven by obliquely propagating whistlers excited through the anomalous Doppler resonance. The oblique whistlers then nonlinearly couple to the weakly damped parallel-propagating whistler wave, which grows with a wavelength determined by the condition that its group velocity is sufficiently small so that the wave remains close to the point of injection of the electron beam. This situation arises specifically in our simulations due to the fact that the beam is continually injected from one end. In support of this point we note that if we increase the magnetic field in our simulation to a value $\unicode[STIX]{x1D714}_{c}/\unicode[STIX]{x1D714}_{pe}=1.0$ the whistler wave does not form. We attribute this to the increased longitudinal group velocity of the whistler so that the stagnation structure is no longer stationary in the simulation frame.

In order to investigate the beam dynamics theoretically we have performed a linear stability analysis based on solution of the linearized, relativistic Vlasov equation coupled with Maxwell’s equations (Ichimaru Reference Ichimaru1973). This system leads to a $3\times 3$ tensor equation for the components of the electric field. The vanishing of the determinant of the tensor gives rise to a dispersion relation,

In forming the tensor $\unicode[STIX]{x1D63F}$ we assume that there are two populations of electrons (hot and cold) having different drift velocities and temperatures. To help in the evaluation of the elements of the tensor $\unicode[STIX]{x1D63F}$ we use a simplified relativistic Maxwell’s distribution (Watson, Bludman & Rosenbluth Reference Watson, Bludman and Rosenbluth1960; Hao et al. Reference Hao, Ding, Sheng, Ren and Zhang2009),

Here, the quantities $\unicode[STIX]{x1D6FE}_{s}T_{s}$ and $\unicode[STIX]{x1D6FE}_{s}^{3}T_{s}$ represent perpendicular and parallel temperatures of species $s$ , and $n_{s}$ and $p_{s}$ are the species density and mean parallel momentum. The advantage of using this distribution is that the elements of $\unicode[STIX]{x1D63F}$ can be approximated using the familiar non-relativistic plasma dispersion function, and using modified Bessel functions to account for non-zero Larmor radius effects. Details of the evaluations may be found in the references Watanabe (Reference Watanabe1991) and Hao et al. (Reference Hao, Ding, Sheng, Ren and Zhang2009).

Figure 3. Growth rates of theoretically predicted unstable modes in two-dimensional wavenumber space for different external magnetic fields, $\unicode[STIX]{x1D714}_{c}/\unicode[STIX]{x1D714}_{pe}=$  (a) 0.01, (b) 0.30 and (c) 0.70, and (d) growth rate and frequency of oblique whistler versus $k_{z}$ for selected values of $k_{x}$ for parameters of case (b).

Solving (3.2) numerically, we obtain the growth rate, $\text{Im}[\unicode[STIX]{x1D714}]$ , of the most unstable mode as a function of wave vector ( $k_{z},k_{x}$ ). These results are plotted as false colour images in figure 3 for the case of (a)  $\unicode[STIX]{x1D714}_{c}/\unicode[STIX]{x1D714}_{pe}=0.01$ , (b)  $\unicode[STIX]{x1D714}_{c}/\unicode[STIX]{x1D714}_{pe}=0.3$ and (c)  $\unicode[STIX]{x1D714}_{c}/\unicode[STIX]{x1D714}_{pe}=0.7$ . In the three cases the parameters of the electron beam, such as the drift velocity, temperature and density are the same as those in figures 1 and 2. The drift velocity and the number density of the background electrons are determined by the charge and current neutrality conditions. We then pick the parameters for the simplified distribution, equation (3.3), to populate a range of parallel momentum similar to what is observed in the simulation, specifically we use $\unicode[STIX]{x1D6FE}_{b}=2$ , which corresponds to $p_{s}=1.9m_{e}c$ .

The beam temperature and the background electron temperature are set to be 100 keV and 10 keV, respectively. These values are the same as in the simulation runs.

As shown in figure 3 there are unstable modes occupying different regions of wavenumber space. One mode occupies a region of wavenumber space with growth rate maximum for $ck_{z}/\unicode[STIX]{x1D714}_{pe}\approx 1.3$ and $k_{x}=0$ , and is identified as a two-stream instability. The phase velocity of the mode with peak growth rate is found to be $\unicode[STIX]{x1D714}/k_{z}=0.7c$ . Thus, the mode is in Cherenkov resonance with the hot electrons. This mode is suppressed but not eliminated by the applied magnetic field. Instability also appears in the region where $ck_{z}<\unicode[STIX]{x1D714}_{pe}$ and $ck_{x}\approx \unicode[STIX]{x1D714}_{pe}$ . Two different mode types are present, depending on the strength of the magnetic field. For low magnetic field values, as in figure 3(a), growth is peaked at $k_{z}=0$ and the real part of the mode frequency is small. We identify this mode as the Weibel instability. For larger magnetic field values, as indicated in figure 3(b), the Weibel instability is suppressed. However, growth occurs for values of wavenumber $ck_{z}/\unicode[STIX]{x1D714}_{pe}\approx 0.1$ and $ck_{x}/\unicode[STIX]{x1D714}_{pe}=0.5$ . A plot of the real frequency and growth rate of this mode is shown in figure 3(d) and indicates that the mode is an off-angle whistler wave as described by (2.2). It has a maximum growth rate in the range $0.01\unicode[STIX]{x1D714}_{pe}$ and a real frequency in approximate agreement with the dispersion curves in figure 3(d). One thing to note from figure 3(d) is that the group velocity in the transverse direction for the off-angle whistler wave is small.

We note that for the large beam speed, $u_{0}=0.95c$ , considered here, the two-stream and Weibel/whistler mode branches are completely separated in wavenumber space, and the two-stream growth is confined to small transverse wavenumbers. For other parameters, specifically lower beam speeds, a different class of modes, ‘oblique modes’, has been found (Bret, Firpo & Deutsch Reference Bret, Firpo and Deutsch2004, Reference Bret, Firpo and Deutsch2005; Taguchi et al. Reference Taguchi, Antonsen and Mima2004), for which two-stream growth extends to large transverse wavenumber. The method we use here, in particular (3.2), returns these modes when the appropriate beam and plasma parameters are imposed.

Figure 4. Temporal evolution of the spatial spectrum of the out-of-plane magnetic field in the case of $\unicode[STIX]{x1D714}_{c}/\unicode[STIX]{x1D714}_{pe}=0.3$ , which are obtained by Fourier transform in the domain $50<z/d_{e}<175$ .

In order to further investigate the interaction of unstable modes in the simulation, we have performed spatial, two-dimensional Fourier transforms of the out-of-plane magnetic field $\overline{B}_{y}$ in the region near the beam entrance ( $50<z/d_{e}<175$ ). Figure 4(a–d) shows the results for the case of $\unicode[STIX]{x1D714}_{c}/\unicode[STIX]{x1D714}_{pe}=0.3$ at different times in the simulation. As shown in figure 4(a,b), longitudinal unstable modes ( $ck_{z}/\unicode[STIX]{x1D714}_{pe}\sim 1$ ) grow during the first stage. (While the spectrum of longitudinal modes appears to be peaked off axis, $k_{x}\neq 0$ , it must be recalled that we are plotting the out of plane magnetic field spectral density, which carries an additional factor $k_{x}^{2}$ relative to the longitudinal electric field spectral density.) These modes generate a large modulation in the beam density, and the nonlinear coupling generates higher harmonic modes as shown in figure 4(b).

Figure 4(b) shows the emergence at small longitudinal wavenumber ( $ck_{z}/\unicode[STIX]{x1D714}_{pe}\sim 0.1$ ) of disturbances that we interpret to be the off-angle whistler waves. These modes grow slowly and begin to dominate the longitudinal modes as shown in figure 4(c). By the time of figure 4(d), the original longitudinal modes are no longer discernible and the spectrum of fluctuations has condensed into a peaked distribution at small wavenumbers corresponding to the fields in figure 2.

Figure 5. Transverse system size dependence of the flux of positively propagating beam electrons at $z/d_{e}=400$ in the case of $\unicode[STIX]{x1D714}_{c}/\unicode[STIX]{x1D714}_{pe}=0.3$ .

The role of fields with transverse structure and associated mode couplings can be examined by varying the transverse size of the simulation domain. We note that the whistler mode that grows to stop the beam has essentially zero transverse wavenumber. However, whether or not it appears depends on the transverse size of the simulation domain. Figure 5 shows the dependence of the flux ratio, $R_{flux}$ , (flux of hot electrons leaving on the right to flux entering on the left) on the transverse size of the simulation box, $L_{x}/d_{e}$ for the case of $\unicode[STIX]{x1D714}_{c}/\unicode[STIX]{x1D714}_{pe}=0.3$ . The flux ratio is unity in the absence of the large amplitude whistler wave. The figure shows that as the transverse dimension is increased, the flux ratio drops, indicating the effect of fields with transverse structure on the growth of the longitudinal whistler wave. The ratio stops changing when $L_{x}/d_{e}>50$ , which corresponds to about $2.5~\unicode[STIX]{x03BC}\text{m}$ for $n_{0}=10^{22}~\text{cm}^{-3}$ , and indicating that the fields have important spectral content for $k_{x}d_{e}>0.02$ . This result further supports contention that the generation of the obliquely propagating modes plays a crucial role for the amplification of the large whistler wave.