2 Plasma mode analysis

Consider a cylindrical column of radius $a$. This cylindrical column is filled with a plasma of density $n_e$. The plasma in this cylindrical plasma column is supported by a longitudinal magnetic field $B_s$. A rotating electron beam is launched axially to this column for plasma interaction. Consider the plasma equilibrium is perturbed by an electrostatic perturbation with electrostatic potential given by

(2.1)\begin{equation} \phi=\phi(r)\exp[-\iota(\omega t-\beta\theta-k_z z)], \end{equation}

where $\phi (r)=-({n_{0b}(r)e}/{\epsilon _0 k_z^2})$, $\beta$ is the azimuthal mode number, $\epsilon _0$ is electric permittivity of free space, $k=\widehat{z}k_z$ is the propagation vector, $e$ is the electron charge, $\omega$ is the angular frequency of the perturbation and $\theta$ is the angular co-ordinate. The electron beam density profile is taken as $n_{0b}(r)=({N_0^0}/{2{\rm \pi} a})\delta (r-a)$, where $N_0^0$ is the number of beam electrons per unit axial length, $\delta$ is the Dirac delta function. The static magnetic field is applied along the $z$-direction as $B_s\parallel \hat {z}$.

The response of plasma electrons to the fields can be governed by the equation of motion

(2.2)\begin{equation} m\frac{{\rm d}\boldsymbol{V}}{{\rm d} t}=e\boldsymbol{\nabla}\phi - e(\boldsymbol{V}\times \boldsymbol{B}_s) \end{equation}

and the continuity equation

(2.3)\begin{equation} \frac{\partial n}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}{n\boldsymbol{V}}=0. \end{equation}

The linearization of these equations gives the velocity and density perturbations as (Liu and Tripathi Reference Liu and Tripathi1994)

(2.4)\begin{gather} \boldsymbol{V}_{1\perp e}=\frac{e}{m}\frac{\boldsymbol{\nabla}_\perp\phi\times \boldsymbol{\omega}_c +\iota\omega\boldsymbol{\nabla}_\perp\phi}{\omega^2-\omega_c^2}, \end{gather}

(2.5)\begin{gather}\boldsymbol{V}_{1z}=\frac{-ek_z\phi}{m\omega}, \end{gather}

(2.6)\begin{gather}n_{1e}={-}\frac{n_0^0e}{m}\left(-\frac{\nabla_\perp^2\phi}{\omega^2-\omega_c^2}+ \frac{k_z^2\phi}{\omega^2}\right), \end{gather}

where $-e$ and $m$ are the electron charge and mass, respectively. The sign $\perp$ indicates the direction perpendicular to the magnetic field. Here, the electron cyclotron frequency is defined as ${\omega _c={e{B_s}}/{m}}$.

The electrostatic potential excited in the plasma can be obtained using Poisson's equation. We write Poisson's equation in cylindrical coordinates (assuming $\partial /\partial z=0$) as

(2.7)\begin{equation} \frac{\partial^2\phi(r)}{\partial r^2}+\frac{1}{r}\frac{\partial\phi(r)}{\partial r}+ \frac{1}{r^2}\frac{\partial^2\phi(r)}{\partial\theta^2}={-}\frac{\rho_1}{\epsilon_0}, \end{equation}

where $\rho _1$ is the charge density associated with the electron density given by (2.6). Using (2.6) in (2.7), we get

(2.8)\begin{equation} \frac{\partial^2\phi(r)}{\partial r^2}+\frac{1}{r}\frac{\partial\phi(r)}{\partial r}+ \left(\frac{\left(\dfrac{\omega_p^2}{\omega^2}-1\right)k_z^2}{\left(1-\dfrac{\omega_p^2}{\omega^2-\omega_c^2}\right)}- \frac{\beta^2}{r^2}\right)\phi(r)=0, \end{equation}

where $\omega _p^2={n_e e^2}/{\epsilon _0}{m}$.

The solution of above equation can be considered as

(2.9)\begin{equation} \phi(r)=\alpha \varGamma_{\beta,n} J_\beta(pr), \end{equation}

where $\varGamma _{\beta,n}=[({2}/{a^2})J_{\beta +1}^2(p_n a)]^{1/2}$ is the normalization constant such that $\phi$ must vanish at $r=a$; thus, $(J_\beta (pa)=0)$. Here, $p=p_n={\chi _n}/{a}$ (with $n=1,2,3$), where $\chi _n$ is the zeroth of the Bessel function $J_\beta (x)$.

The electrostatic potential excited in a plasma can be given by (2.9). The expression of different modes of the excited potentials can be obtained using different orders of the Bessel function. In a cylindrical geometry, plasma can support various modes of plasma potential due to geometry effects. This can be seen from the results shown in figure 1. Figure 1 shows the electrostatic potential ($\phi (r)$ in a.u.) of the plasma modes supported in a magnetized plasma column. The solution shown by (2.9) shows that the potential of plasma mode characterizes the extent of localization of the plasma mode and increases with $\omega /\omega _p$. This means that the potential of the plasma modes is more strongly localized near $\omega _p$, which can support the excitation of electromagnetic fields by coupling with the spiralling electron beam in the presence of a magnetic field. We study the coupling of plasma modes with the electron beam for THz field excitation in the next section.

Figure 1.

Potential of plasma mode for different azimuthal quantum numbers $\beta =0$, $\beta =1$, $\beta =2$, $\beta =3$ with normalization constant $\alpha =2\times 10^6\,{\rm V}\,{\rm cm}^{-1}$, plasma column radius $a=3.6\times 10^{-4}\,{\rm cm}$, $n_b= 16\times 10^{16}\,{\rm cm}^{-3}$ and $n_e= 8\times 10^{16}\,{\rm cm}^{-3}$.

3 The THz field generation

We consider the propagation of an electron beam of density $n_{0b}$ and velocity $V_b$ with spreading in the azimuthal direction. The beam acquires an oscillatory velocity due to the excited plasma modes. The governing equation for electron motion in the presence of an electron beam can be written as

(3.1)\begin{equation} m\frac{{\rm d}}{{\rm d} t}(\gamma {\boldsymbol{V}})=e{\boldsymbol{\nabla}\phi}-\frac{e}{c}(\boldsymbol{V}\times\boldsymbol{B_s}), \end{equation}

where $\gamma =(1-{{V}^2}/{c^2})^{-({1}/{2})}$, c is the velocity of light in vacuum.

We solve the above equation by using the perturbation approach in the relativistic regime. The physical quantities are expanded in the form of power series expansions in term of fields. We expressed these expansions as $X=X^{(0)}+\sum _{n=1}^{\infty } \alpha ^{(l+1)} X^{(n)}$, where $X=(v)$ and $X^{(0)}=(1)$. Here, $\alpha$ is a normalized parameter which defines the scale of nonlinearity. The source of nonlinearity is the beam electric field for this calculation, hence, $\alpha$ can be expressed in terms of the beam electric field. We apply a first-order perturbation to obtain the perturbed electron velocity as $\boldsymbol {V}={\boldsymbol {V}}_{0b}+{\boldsymbol {V}}_{1b}$. Here, subscript $b$ is for the beam and $\boldsymbol {V}_{1b}$ has components $\boldsymbol {V}_r$ and $\boldsymbol {V}_\theta$.

To obtain the perturbed quantities, we use $\gamma =\gamma _0+\gamma _0^3{\boldsymbol {V}}_{0b}\boldsymbol {\cdot }{\boldsymbol {V}}_{1z}/c^2$ and $n=n_{0b}+n_{1b}$, where ${\boldsymbol {V}}_{1br}=-\iota (\bar {\omega }-{q\omega _c}/{\gamma _0}){\boldsymbol {r}}_1,$ ${\boldsymbol {V}}_{1b\theta }={\boldsymbol {r}}_1({\omega _c}/{\gamma _0}) + {\boldsymbol {r}}_0\dot {\boldsymbol {\theta }_1},$ ${\boldsymbol {V}}_{1bz}=-\iota (\bar {\omega }-{\beta \omega _c}/{\gamma _0}){\boldsymbol {z}}_1$ with

(3.2)\begin{gather} \boldsymbol{r}_1={-}\frac{e}{m_b}\frac{\left[\gamma_0 \left(\bar{\omega}-\dfrac{\beta\omega_c}{\gamma_0}\right)\boldsymbol{E}_{1r}-\iota \left(\dfrac{\omega_c}{\gamma_0}\right)\boldsymbol{E}_{1\theta}\right]}{\left(\bar{\omega}- \dfrac{\beta\omega_c}{\gamma_0}\right)\left[\dfrac{\omega_c^2}{\gamma_0^2}-\gamma_0^2 \left(\bar{\omega}-\dfrac{\beta\omega_c}{\gamma_0}\right)^2\right]}, \end{gather}

(3.3)\begin{gather}\boldsymbol{\theta}_1={-}\frac{e}{m_b{\boldsymbol{r}}_0} \frac{\left[\gamma_0\left(\bar{\omega}-\dfrac{\beta\omega_c}{\gamma_0}\right) \boldsymbol{E}_{1\theta}+\iota\left(\dfrac{\omega_c}{\gamma_0}\right) \boldsymbol{E}_{1r}\right]}{\left(\bar{\omega}-\dfrac{\beta\omega_c}{\gamma_0}\right) \left[\dfrac{\omega_c^2}{\gamma_0^2}-\gamma_0^2\left(\bar{\omega}- \dfrac{\beta\omega_c}{\gamma_0}\right)^2\right]}, \end{gather}

(3.4)\begin{gather}{\boldsymbol{z}}_1={-}\frac{e\boldsymbol{E}_{1z}}{m_b\gamma_0^3 \left(\bar{\omega}-\dfrac{q\omega_c}{\gamma_0}\right)^2}, \end{gather}

where $\bar \omega =\omega -k_z V_{0\parallel }$.

Solving the equation of continuity ${\partial n}/{\partial t}+ \boldsymbol {\nabla }\boldsymbol {\cdot }{(n\boldsymbol {V}})=0$, one obtains

(3.5)\begin{equation} n_{1b}=\frac{n_{0b}\boldsymbol{k_z}\boldsymbol{\cdot}\boldsymbol{V}_{1b}}{\omega}. \end{equation}

Now the perturbed current density can be written as

(3.6)\begin{equation} {\boldsymbol{J}}_{1b}={-}e(n_{1b}{\boldsymbol{V}}_{0b}+n_{0b}{\boldsymbol{V}}_{1b}). \end{equation}

We define $\boldsymbol {E_T}$ and $\boldsymbol {H_T}$ as THz field structures in the absence of a beam. In the presence of a beam current, let the fields be ${\boldsymbol {E}}=A_1(t){\boldsymbol {E}}_T$ and magnetic field ${\boldsymbol {B}}=A_2(t){\boldsymbol {B}}_T$. Here, the subscript $T$ is used for THz. The above fields satisfy the following Maxwell's equations:

(3.7)\begin{gather} \boldsymbol{\nabla}\times{\boldsymbol{E}}={-}\frac{1}{c}\frac{\partial {\boldsymbol{B}}}{\partial t}, \end{gather}

(3.8)\begin{gather}\boldsymbol{\nabla}\times{\boldsymbol{B}}=\frac{4{\rm \pi}}{c}{\boldsymbol{J}}_{1b}+\frac{\epsilon}{c}\frac{\partial {\boldsymbol{E}}}{\partial t}. \end{gather}

Here, $\epsilon =1-\omega _p^2/\omega ^2$.

Using the above solutions in (3.7) and (3.8), we obtain

(3.9)\begin{gather} \frac{\partial A_2}{\partial t}={-}\iota\omega(A_1-A_2), \end{gather}

(3.10)\begin{gather}\left[\epsilon\frac{\partial A_1}{\partial t}-\iota\omega\epsilon(A_1-A_2)\right]\boldsymbol{E}_T={-}4{\rm \pi}{\boldsymbol{J}}_{1b}. \end{gather}

Solving (3.9) and (3.10), assuming ${\partial A_1}/{\partial t}\approx {\partial A_2}/{\partial t}$ (for high-frequency electromagnetic radiation), we obtain

(3.11)\begin{equation} {-}2\iota\epsilon\omega(A_1-A_2){\boldsymbol{E}}_T={-}4{\rm \pi}{\boldsymbol{J}}_{1b}. \end{equation}

Taking the dot product of this equation by ${\boldsymbol {E}}_{T}r\,{\rm d} r$ and integrating over $r$ from 0 to $\infty$, we get

(3.12)\begin{equation} \frac{\partial A_1}{\partial t}=\frac{2{\rm \pi}}{\epsilon} \frac{\displaystyle\int_{0}^{\infty}{\boldsymbol{J}}_{1bz}{\boldsymbol{E}}_{Tz}^*r\,{\rm d} r}{ \displaystyle\int_{0}^{\infty} {\boldsymbol{E}}_T\boldsymbol{\cdot}{\boldsymbol{E}}_T^* r\,{\rm d} r}. \end{equation}

By simplifying, we get all components of the generated THz electric field as follows:

(3.13)\begin{gather} {\boldsymbol{E}}_{Tr}=\frac{2{\rm \pi} e^2n_{0b}l}{\omega\epsilon m_b(A_1-A_2)} \left(\frac{\boldsymbol{k}{\boldsymbol{V}}_{0b}}{\omega}+1\right) \frac{\left[\gamma_0\left(\bar{\omega}-\dfrac{\beta\omega_c}{\gamma_0}\right) \boldsymbol{\nabla}\phi_{r}+\iota\left(\dfrac{\omega_c}{\gamma_0}\right) \boldsymbol{\nabla}\phi_{\theta}\right]}{\left[\dfrac{\omega_c^2}{\gamma_0^2}-\gamma_0^2 \left(\bar{\omega}-\dfrac{\beta\omega_c}{\gamma_0}\right)^2\right]}, \end{gather}

(3.14)\begin{gather} {\boldsymbol{E}}_{T\theta}= \frac{2{\rm \pi}\iota e^2n_{0b}}{\omega\epsilon m_b (A_1-A_2)\left(\bar{\omega}- \dfrac{\beta\omega_c}{\gamma_0}\right)\left[\dfrac{\omega_c^2}{\gamma_0^2}- \gamma_0^2\left(\bar{\omega}-\dfrac{l\omega_c}{\gamma_0}\right)^2\right]} \left(\frac{{\boldsymbol{k}}{\boldsymbol{V}}_{0b}}{\omega}+1\right) \nonumber\\ \times \left[\omega_c\left(\bar{\omega}-\frac{\beta\omega_c}{\gamma_0}\right)\boldsymbol{\nabla}\phi_{r}+ \iota\frac{\omega_c^2}{\gamma_0^2}\boldsymbol{\nabla}\phi_\theta+\iota\omega(-\gamma_0) \left(\bar{\omega}-\frac{\beta\omega_c}{\gamma_0}\right)\boldsymbol{\nabla}\phi_{\theta}- \iota\frac{\omega_c}{\gamma_0}\boldsymbol{\nabla}\phi_{r}\right], \end{gather}

(3.15)\begin{gather} {\boldsymbol{E}}_{Tz}={-}\frac{2{\rm \pi} e^2n_{0b}}{\omega\epsilon m_b\gamma_0^3(A_1-A_2)} \frac{\left(\dfrac{\boldsymbol{k}{\boldsymbol{V}}_{0b}}{\omega}+1\right)}{ \left(\bar\omega-\dfrac{\beta\omega_c}{\gamma_0}\right)}\boldsymbol{\nabla}\phi_z. \end{gather}

Equations (3.13)–(3.15) show the excited electromagnetic fields generated by the spiral electron beam interaction with a plasma. We estimate the transverse component of the electromagnetic field (THz field) for various electron beam parameters. The dependency of the excited THz fields on the electron beam velocity and electron beam density is illustrated in figure 2. Figure 2(a) shows the spectral variation of the THz field (in GV m$^{-1}$) for different electron beam velocities. The results depicted in figure 2(a) show that the beam velocity plays an important role in THz field excitation. The THz field increases with the electron beam velocity due to large transverse current generation during the electron beam interaction with plasmas. The large beam velocity obviously contributes profoundly to the transverse current density because the current density is directly proportional to the charge carrier velocity. Thus, the THz field strength is enhanced by the beam velocity, as confirmed by these results. We know that the THz field peaks near $\omega _p$ due to resonance. Our results shows that the THz field is maximized near $\omega _p$, as shown in the results of figure 2. The deviation in plasma frequency diverts the pump energy, hence, the transverse current may be varied accordingly. On the other hand, the applied magnetic field supports the transverse oscillations, which contribute via electron cyclotron resonance. Figure 2(b) shows the spectral distribution of the THz field for different beam velocities. It is clear from (3.13) that the radial electric field increases with the beam velocity. Consequently, the THz field increases with the beam velocity. The higher electron beam density enhances the plasma oscillations via strong excitation of the plasma density perturbation. In a result, the transverse current associated with the plasma oscillations will be stronger for the large beam density. A stronger peak THz field for the large beam density can be confirmed from the results shown in figure 2(b).

Figure 2.

Spectral distribution of THz field (in GV m$^{-1}$) for (a) different electron beam velocities $V_{0b}=2\times 10^{10}\,{\rm cm}\,{\rm s}^{-1}$, $V_{0b}=2.3\times 10^{10}$ cm s$^{-1}$, $V_{0b}=2.7\times 10^{10}$ cm s$^{-1}$, $V_{0b}=3\times 10^{10}$ cm s$^{-1}$ with $n_e=8\times 10^{16}\,{\rm cm}^{-3}$, $n_{0b}=16\times 10^{16}\,{\rm cm}^{-3}$ and (b) different electron beam densities $n_{0b}=10\times 10^{16}\,{\rm cm}^{-3}$, $n_{0b}=16\times 10^{16}\,{\rm cm}^{-3}$, $n_{0b}=20\times 10^{16}\,{\rm cm}^{-3},\ n_{0b}=24\times 10^{16}\,{\rm cm}^{-3}$ with $n_e=8\times 10^{16}\,{\rm cm}^{-3}$. The plasma is filled in a plasma column of radius $a=3.6\times 10^{-4}\,{\rm cm}^{-1}$ with a magnetic field of 0.1 T.

4 The PIC simulation results

To validate the proposed theoretical model, we carried out quasi-three-dimensional (quasi-3-D) simulations. In this study, we investigate THz field generation using the interaction of an electron beam with a plasma in a cylindrical column. For this study, these simulations are performed using the spectral, quasi-3-D PIC code FBPIC (Lehe et al. Reference Lehe, Kirchen, Andriyash, Godfrey and Vay2016). This code uses a spectral cylindrical representation, which decreases the computational load in the laser interaction process. FBPIC contains a set of 2-D radial grids which represents an azimuthal mode. The fields are decomposed into two azimuthal modes ($N_m=2$). The simulation grid size is set to be $0.14\,\mathrm {\mu }{\rm m}\times 0.4\,\mathrm {\mu } {\rm m}$ in the longitudinal and transverse directions, respectively (with 8 particles per cell). The moving window with size $170\,\mathrm {\mu }{\rm m}\times 60\,\mathrm {\mu }{\rm m}$ is large enough to simulate plasma oscillations in the first cycle. The parameters for an electron beam used in the simulations are as follows: beam size ($\sigma _r=5\,\mathrm {\mu }{\rm m}$, $\sigma _z=10\,\mathrm {\mu }{\rm m}$), beam energy 28 GeV and beam density range is $n_b= 16\unicode{x2013}24\times 10^{16}\,{\rm cm}^{-3}$. The plasma density taken for these simulations is $n_e= 8\times 10^{16}\,{\rm cm}^{-3}$ with a magnetic field of 0.1 T.

Figure 3 shows the PIC simulation results of the THz field excitation corresponding to the theoretical model given § 3. The excited THz field and the corresponding transverse current density both are given in this figure. The average low-frequency ponderomotive force associated with the electron beam drives plasma oscillations, generating a dynamic plasma current. The finite rotation of the electron beam generates a transverse component of this current, which excites THz radiation. It is evident from the simulation results that the transverse plasma current is generated by a rotating electron beam while interacting with the plasma. Thus the corresponding THz field is finite. To clearly demonstrate the role of the electron beam parameters, we have depicted the THz field evolution in figure 4 for different beam densities and plasma densities. These results are obtained from the simulation data.

Figure 3.

Snapshots of (a) THz electric field ${E_{Tr}}$ (GV m$^{-1}$) and (b) transverse current density ${J_{x}}$ (normalized by $n_c ec$, where $n_c=10^{21} \lambda _0^{-2}\,{\rm cm}^{-3}$ is the critical plasma density) for $n_b= 24\times 10^{16}\,{\rm cm}^{-3}$.

Figure 4.

Simulation results of THz electric field ${E_{Tr}}$ (GV m $^{-1}$) for different (a) electron beam densities with $n_e=8\times 10^{16}\,{\rm cm}^{-3}$ and (b) plasma densities $n_e$ with $n_{0b}=24\times 10^{16}\,{\rm cm}^{-3}$. The other numerical parameters are same as those of figure 3.

THz field increases with the electron beam density, as predicted previously by the theoretical model. The PIC simulation results show that the maximum THz field is approximately 10 GV m$^{-1}$ for the electron beam density of $n_{0b}=16\times 10^{16}\,{\rm cm}^{-3}$. This value of the peak THz field is also very close to the data obtained from the theoretical model (see figure 2b). This peak THz field corresponds to the resonance condition ($\omega \thickapprox \omega _p$). The peak THz field also depends on the plasma density. The maximum THz field is enhanced by the plasma density. The peak THz field position is also shifted due to the variation of $\omega _p$, as seen in figure 4(b). The radiation field is peaked at $\omega _p$ for a given plasma density. Thus, the peak radiation field is tuneable with the beam density and velocity. These simulation results are consistent with the theoretical findings. However, a deviation in the THz field profile can be observed as the simulation results are inconsistent with the envelope model.