Excitation of wakefields in carbon nanotubes: a hydrodynamic model approach

Carbon nanotubes (CNTs) were discovered by Iijima in 1991 [1] and they can be thought of as a sheet of graphene (a hexagonal lattice of carbon) rolled into a cylinder. CNTs can exhibit metallic or semiconductor properties depending on their rolling pattern (i.e. on their radius and geometrical angle). Thus, as a consequence of their unique thermo-mechanical and electronic properties and dimensional flexibility, CNTs have been widely studied in both theoretical and experimental aspects. As a hollow structure, CNTs may be used for transporting and focusing charged particles similar to crystal channeling. In particular, experimental results on 2 MeV $\textrm{He}^{+}$ ions [2] and 300 keV electrons [3] channeling in CNTs have been reported.

On the other hand, solid-state wakefield acceleration using crystals was proposed in the 1980s and 1990s by Tajima et al [4–6] as an alternative particle acceleration technique to sustain TV m⁻¹ acceleration gradients. In the original Tajima's conceptual scheme [4], a longitudinal electric wakefield is excited by a laser (laser driven) in the crystal so that a properly injected witness beam may be accelerated. Similarly, the ultrashort charged particle bunches (beam driven) can excite electric wakefields so that the energy loss of the driving bunch can be transformed into an increment of energy of a witness bunch. However, the angstrom-size channels of natural crystals pose a limitation for the beam intensity acceptance and the dechanneling rate. In this context, CNTs can obtain wider channels in two dimensions and longer dechanneling lengths [7, 8], which together with their remarkable electronic properties, larger degree of dimensional flexibility and thermo-mechanical strength, make them a robust candidate for TeV m⁻¹ acceleration. Consequently, carbon nanostructures (CNTs or even graphene layers) are currently being widely studied for wakefield acceleration [9, 10].

Wakefields in CNTs are excited through the collective oscillation of electrons on the nanotube's surface, often referred to as plasmons. This excitation arises from the interaction between the driving bunch and the CNT, leading to the generation of these wakefields. The electronic excitations on the single-wall CNT surface produced by the interaction with charged particles have been theoretically studied using a dielectric theory [11–13], a hydrodynamic model [14–16], a two-fluid model [17, 18], a quantum hydrodynamic model [19] and a kinetic model [20–23]. However, while these articles mostly explore properties such as the energy loss and stopping power or, at most, evaluate the induced surface electron density and/or induced potential [16, 22, 23], they do not address the induced longitudinal and transverse wakefields, which could provide acceleration and focusing, respectively, for a witness charge.

Thus, this article is motivated by the need to study the wakefields excited by charged particles moving paraxially inside CNTs. Furthermore, this article presents analytical expressions that enable the rapid optimization of CNT parameters. This optimization aims to achieve the maximum longitudinal wakefield for applications in particle acceleration.

This work is organized as follows. In section 2 the general expressions are derived for the longitudinal and transverse wakefields excited by the interaction of a charged particle with a CNT in the realm of the hydrodynamic model. This model has been chosen because of its simplicity and its good agreement with the dielectric formalism in random-phase approximation [15]. In section 3, after investigating the influence of different model parameters in terms of dispersion relations, the CNT parameters are optimized to achieve the highest longitudinal wakefield for a given driving velocity. Finally, the main conclusions of this study are presented in section 4.

In this work, a linearized hydrodynamic theory [14, 15, 18] is adopted, in particular a model which includes single-electron based excitations on nanotube surfaces and was described by Wang and Mišković in [15], although we will use the SI units instead of atomic units. In this theory, a single-wall CNT is modelled as an infinitesimally thin and infinitely long cylindrical shell with a radius a. The delocalized electrons of the carbon ions are considered as a two-dimensional free-electron gas that is confined over the cylindrical surface of the CNT with a uniform surface density n₀. It is considered a driving point-like charge Q travelling parallel to the z-axis inside the tube with a constant velocity v (see figure 1). Consequently, its position as a function of time t is $\mathbf{r}_0(t) = \left(r_0, \varphi_0, vt\right)$ in cylindrical coordinates. As a consequence of the presence of the driving charge Q, the homogeneous electron gas will be perturbed and can be modelled as a charged fluid with a velocity field $\mathbf{u}(\mathbf{r}_a ,t)$ and surface density $n(\mathbf{r}_a ,t) = n_0+n_1(\mathbf{r}_a ,t)$, where $\mathbf{r}_a = (a, \varphi_a,z_a)$ are the coordinates of a point at the cylindrical surface of the tube and $n_1(\mathbf{r}_a ,t)$ is the perturbed density per unit area. In the linearized hydrodynamic model, it is assumed that the perturbed density n₁ and the fluid velocity u are small perturbations (the validity of this linear approximation is demonstrated in appendix A). As the electron gas is confined to the cylindrical surface, the normal component to the surface of the tube of the velocity field u is zero. Since the carbon ions are much heavier than the electrons, the time scale associated with the ionic motion is orders of magnitude slower than that of the electronic motion. Hence, for the purpose of investigating the wakefield dynamics, the ionic motion can be neglected [24, 25].

Zoom In Zoom Out Reset image size

In the linearized hydrodynamic model, the electronic excitations on the tube wall can be described by three differential equations: (i) the continuity equation

$$\begin{equation} \frac{\partial n_1\left(\mathbf{r}_a, t\right)}{\partial t}+n_0 \nabla_{\|} \cdot \mathbf{u}\left(\mathbf{r}_a, t\right) = 0, \end{equation} \tag{ 1 }$$

(ii) the Poisson's equation

$$\begin{equation} \nabla^2 \Phi\left(\mathbf{r}, t\right) = \frac{1}{\varepsilon_0}\left[e n_1\left(\mathbf{r}_a, t\right) \delta\left(r-a\right)-Q \delta\left(\mathbf{r}-\mathbf{r}_0\right)\right], \end{equation} \tag{ 2 }$$

and (iii) the momentum-balance equation

$$\begin{equation} \begin{split} \frac{\partial \mathbf{u}\left(\mathbf{r}_a, t\right)}{\partial t} = \frac{e}{m_e} \nabla_{\|} \Phi\left(\mathbf{r}_a, t\right)-\frac{\alpha}{n_0} \nabla_{\|} n_1\left(\mathbf{r}_a, t\right)+\frac{\beta}{n_0} \nabla_{\|}\left[\nabla_{\|}^2 n_1\left(\mathbf{r}_a,t\right)\right]-\gamma \mathbf{u}\left(\mathbf{r}_a, t\right). \end{split} \end{equation} \tag{ 3 }$$

where we have retained only the first-order terms in n₁ and u. In these equations, $\mathbf{r} = (r,\varphi,z)$ is the position vector, $\nabla = \hat{\mathbf{r}}\frac{\partial}{\partial r} +\hat{\boldsymbol{\varphi}}\frac{1}{r} \frac{\partial}{\partial \varphi} +\hat{\mathbf{z}}\frac{\partial}{\partial z}$, $\nabla_{\|} = \hat{\boldsymbol{\varphi}}\frac{1}{a} \frac{\partial}{\partial \varphi} +\hat{\mathbf{z}}\frac{\partial}{\partial z}$ differentiates only tangentially to the tube surface, Φ the electric scalar potential, e the elementary charge, m_e the rest mass of the electron, ε₀ the vacuum electric permittivity and δ the Dirac delta. Equation (3) shows the sum of four different contributions. The first term in the right-hand side is the force on electrons on the nanotube surface due to the tangential component of the electric field generated by the driving charge Q and the consequent perturbed density n₁. The second and third terms are related to the parts of the internal interaction force in the electron gas. In particular, the second term takes into account the possible coupling with acoustic modes defining the parameter $\alpha = v^2_F/2$ (in which $v_F = \hbar(2\pi n_0)^{1/2}/m_e$ is the Fermi velocity of the two-dimensional electron gas; $\hbar$ is the reduced Planck constant), and the third term is a quantum correction that arises from the functional derivative of the Von Weizsacker gradient correction in the equilibrium kinetic energy of the electron fluid [18] and describes single-electron excitations in the electron gas, where the parameter $\beta = \frac{1}{4}(\frac{\hbar}{m_e})^2$ has been defined. The last term is introduced to satisfy the non-conservation of the system and represents a frictional force on electrons due to scattering with the ionic-lattice charges, where γ is the friction parameter. The friction parameter may be also used as a phenomenological parameter to take into account the broadening of the plasmon resonance in the excitation spectra of different materials [11].

Taking into account that the electric potential vanishes at $r \rightarrow \infty$ and is finite at the origin r = 0, the potential can be expanded in terms of the modified Bessel functions $I_m(x)$ and $K_m(x)$ of integer order m (i.e. a Fourier-Bessel expansion). The total potential inside the nanotube (r < a) can be calculated as $\Phi_{\textrm{in}} = \Phi_0+\Phi_{\textrm{ind}}$ with the Coulomb potential Φ₀ due to the driving charge and the induced potential $\Phi_{\textrm{ind}}$ due to the perturbation of the electron fluid on the CNT surface. The total potential outside the nanotube (r > a) will be denoted as $\Phi_{\textrm{out}}$. Thus, the Fourier-Bessel expansion of these three components is

$$\begin{equation}\Phi_0\left(r, \varphi, \zeta\right) = \frac{1}{4 \pi \varepsilon_0} \frac{Q}{\left\|\mathbf{r}-\mathbf{r}_{\mathbf{0}}\right\|} = \frac{Q}{4 \pi^2 \varepsilon_0} \sum_{m = -\infty}^{+\infty} \int_{-\infty}^{+\infty} \mathrm{d} k\, e^{i k \zeta+i m\left(\varphi-\varphi_0\right)} I_m\left(|k| r_{\textrm{min}}\right) K_m\left(|k| r_{\textrm{max}}\right),\end{equation} \tag{ 4 }$$

$$\begin{equation} \begin{split} \Phi_{\textrm{ind}}\left(r, \varphi, \zeta\right) = \frac{Q}{4 \pi^2 \varepsilon_0} \sum_{m = -\infty}^{+\infty} \int_{-\infty}^{+\infty} \mathrm{d} k\,e^{i k \zeta+i m\left(\varphi-\varphi_0\right)} I_m\left(|k| r_0\right) I_m\left(|k| r\right) A_m\left(k\right), \end{split} \end{equation} \tag{ 5 }$$

$$\begin{equation} \begin{split} \Phi_{\textrm{out}}\left(r, \varphi, \zeta\right) = \frac{Q}{4 \pi^2 \varepsilon_0} \sum_{m = -\infty}^{+\infty} \int_{-\infty}^{+\infty} \mathrm{d} k \,e^{i k \zeta+i m\left(\varphi-\varphi_0\right)} I_m\left(|k| r_0\right) K_m\left(|k| r\right) B_m\left(k\right), \end{split} \end{equation} \tag{ 6 }$$

where k is the wavenumber, a comoving coordinate $\zeta = z-vt$ has been defined and $r_{\textrm{min}} = \textrm{min}(r,r_0)$, $r_{\textrm{max}} = \textrm{max}(r,r_0)$. The unknown coefficients $A_m(k)$ and $B_m(k)$ can be calculated if the following boundary conditions are imposed: (i) the continuity of the electric potential at the nanotube surface

$$\begin{equation} \left.\Phi_{\textrm{in}}\left(r, \varphi, \zeta\right)\right|_{r = a} = \left.\Phi_{\textrm{out}}\left(r, \varphi, \zeta\right)\right|_{r = a}, \end{equation} \tag{ 7 }$$

and (ii) the discontinuity of the radial component of the electric field due to perturbed density n₁ of the electron fluid

$$\begin{equation} \left.\frac{\partial \Phi_{\textrm{out}}\left(r, \varphi,\zeta\right)}{\partial r}\right|_{r = a}-\left.\frac{\partial \Phi_{\textrm{in}}\left(r, \varphi, \zeta\right)}{\partial r}\right|_{r = a} = \frac{en_1\left(a,\varphi, \zeta\right)}{\varepsilon_0}.\end{equation} \tag{ 8 }$$

Thus, the coefficients $A_m(k)$ and $B_m(k)$ are given by the following non-dimensional functions:

$$\begin{equation} A_m\left(k\right) = \frac{\Omega_p^2 a^2\left(k^2+m^2 / a^2\right) K_m^2\left(|k| a\right)}{D_m\left(k\right)}, \end{equation} \tag{ 9 }$$

$$\begin{equation} B_m\left(k\right) = \frac{D_m\left(k\right)+\Omega_p^2 a^2\left(k^2+m^2 / a^2\right) K_m\left(|k| a\right) I_m\left(|k| a\right)}{D_m\left(k\right)}, \end{equation} \tag{ 10 }$$

where $\Omega_p = \sqrt{\frac{e^2 n_0}{\varepsilon_0 m_e a}}$ is the plasma frequency,

$$\begin{equation} D_m\left(k\right) = k v\left(k v+i \gamma\right)-\omega_m^2\left(k\right), \end{equation} \tag{ 11 }$$

and

$$\begin{equation} \begin{split} \omega_m^2\left(k\right) = \alpha\left(k^2+\frac{m^2}{a^2}\right)+\beta\left(k^2+\frac{m^2}{a^2}\right)^2+ \Omega_p^2 a^2\left(k^2+\frac{m^2}{a^2}\right) K_m\left(|k| a\right) I_m\left(|k| a\right). \end{split} \end{equation} \tag{ 12 }$$

Consequently, the resonant excitations occur when $kv = \omega_m(k)$ if the damping or friction factor γ vanishes.

On the other hand, the longitudinal and transverse electric wakefields inside the tube are, respectively,

$$\begin{equation} \begin{aligned} W_z\left(r, \varphi, \zeta\right) & = -\frac{\partial \Phi_{\mathrm{in}}}{\partial z} = \frac{Q}{4 \pi^2 \varepsilon_0} \sum_{m = -\infty}^{+\infty} \int_{-\infty}^{+\infty} \mathrm{d} k\, k e^{i m\left(\varphi-\varphi_0\right)} \left(I_m\left(|k| r_{\textrm{min}}\right) K_m\left(|k| r_{\textrm{max}}\right) \sin \left(k \zeta\right)\right. \\ &\quad +I_m\left(|k| r_0\right) I_m\left(|k| r\right)\left[\textrm{Re}\left[A_m\left(k\right)\right] \sin \left(k \zeta\right)\right. \left.\left.+\textrm{Im}\left[A_m\left(k\right)\right] \cos \left(k \zeta\right)\right]\right) \\ & = W_{z0}+W_{z1}+W_{z2}, \end{aligned} \end{equation} \tag{ 13 }$$

$$\begin{equation} \begin{aligned} W_r\left(r, \varphi, \zeta\right) & = -\frac{\partial \Phi_{\mathrm{in}}}{\partial r} = -\frac{\partial \Phi_{0}}{\partial r}-\frac{Q}{4 \pi^2 \varepsilon_0} \sum_{m = -\infty}^{+\infty} \int_{-\infty}^{+\infty} \mathrm{d} k |k| e^{i m\left(\varphi-\varphi_0\right)} I_m\left(|k| r_0\right) I_m^{^{\prime}}\left(|k| r\right)\\ & \quad \times \left[\textrm{Re}\left[A_m\left(k\right)\right] \cos \left(k \zeta\right) -\textrm{Im}\left[A_m\left(k\right)\right] \sin \left(k \zeta\right)\right] \\ & = W_{r0}+W_{r1}+W_{r2}, \end{aligned} \end{equation} \tag{ 14 }$$

where $I^{^{\prime}}_m(x) = \mathrm dI_m(x)/\mathrm dx$ and the properties $\textrm{Re}\left[A_m(k)\right] = \textrm{Re}\left[A_m(-k)\right]$ and $\textrm{Im}\left[A_m(k)\right] = -\textrm{Im}\left[A_m(-k)\right]$ were considered; $\textrm{Re}$ and $\textrm{Im}$ denote the real and imaginary part, respectively. The previous integrals have been separated in three different terms: the first terms $W_{z0}$ and $W_{r0}$ come from the Coulomb potential and the other terms $W_{z1}, W_{z2}$ and $W_{r1}, W_{r2}$ from the induced potential. To reduce the computational time and prevent artificial numerical errors, a cutoff for large wavenumbers k has been introduced in the numerical integration. It is important to note that the third terms $W_{z2}$ and $W_{r2}$ can be analytically integrated (cf appendix B) if the damping factor vanishes ($\gamma \rightarrow 0^+$):

$$\begin{equation} \begin{aligned} W_{z2}\left(r, \varphi, \zeta\right)& = \frac{Q}{4 \pi^2 \varepsilon_0} \sum_{m = -\infty}^{+\infty} \int_{-\infty}^{+\infty} \mathrm{d} k\, k e^{i m\left(\varphi-\varphi_0\right)} I_m\left(|k| r_0\right) I_m\left(|k| r\right)\textrm{Im}\left[A_m\left(k\right)\right] \cos \left(k \zeta\right) \\ & = \frac{-Q}{2 \pi \varepsilon_0} \sum_{m = -\infty}^{+\infty} e^{i m\left(\varphi-\varphi_0\right)} k_m I_m\left(k_m r_0\right) I_m\left(k_m r\right) \Omega_p^2 a^2\left(k_m^2+\frac{m^2}{a^2}\right) \\ & \quad \times K_m^2\left(k_m a\right)\left|\frac{\partial Z_m}{\partial k}\right|_{k = k_m}^{-1} \cos \left(k_m \zeta\right), \end{aligned} \end{equation} \tag{ 15 }$$

$$\begin{equation} \begin{aligned} W_{r2}\left(r, \varphi, \zeta\right)& = \frac{Q}{4 \pi^2 \varepsilon_0} \sum_{m = -\infty}^{+\infty} \int_{-\infty}^{+\infty} \mathrm{d} k\, |k| e^{i m\left(\varphi-\varphi_0\right)} I_m\left(|k| r_0\right) I_m^{^{\prime}}\left(|k| r\right)\textrm{Im}\left[A_m\left(k\right)\right] \sin \left(k \zeta\right) \\ & = \frac{-Q}{2 \pi \varepsilon_0} \sum_{m = -\infty}^{+\infty} e^{i m\left(\varphi-\varphi_0\right)} k_m I_m\left(k_m r_0\right)I_m^{^{\prime}}\left(k_m r\right) \Omega_p^2 a^2\left(k_m^2+\frac{m^2}{a^2}\right) \\ &\quad \times K_m^2\left(k_m a\right)\left|\frac{\partial Z_m}{\partial k}\right|_{k = k_m}^{-1} \sin \left(k_m \zeta\right), \end{aligned} \end{equation} \tag{ 16 }$$

where it has been defined the quantity $Z_m(k) = \textrm{Re}[D_m(k)] = (k v)^2-\omega_m^2(k)$ and k_m are the (positive) roots of $Z_m(k)$, i.e. the condition of the plasma resonance $k_mv = \omega_m(k_m)$.

As it can be deduced from equations (15) and (16), the roots k_m given by the plasma resonance are essential to describe the behaviour of the wakefields. For this reason, this section begins with a detailed analysis of the dispersion relation.

3.1. Dispersion relation

In the following calculations, unless otherwise indicated, it is assumed that the surface electron density of a single-wall CNT can be approximated by the electron-gas density of a graphite sheet: $n_0 = n_g = 1.53\times10^{20}$ m⁻² [15, 26]. Figure 2 shows the dispersion curves $\omega_m(k)$ for the first modes at different CNT radii. If the radius is too small, the first mode m = 0 does not satisfy the resonance condition for high velocities, while modes with m > 0 have a solution k_m for those velocities, as seen in figure 2(a) at a = 1 nm. If the CNT radius increases, the resonance condition can be satisfied for a wider range of high velocities, as it is depicted in figure 2(b) for a = 100 nm. Furthermore, it can be seen that the modes $\omega_m(k)$ converge for a sufficiently large value of the wavenumber k. Therefore, if the resonance condition k_m is sufficiently large, then all modes will have a similar value of k_m . Note that a large value of k_m indicates that the associated wavelength of the wakefield ($\lambda_m = 2\pi/k_m$) will be smaller. Moreover, if the surface density n₀ increases, then the dispersion curves increase and the resonance conditions are obtained for higher k_m .

Zoom In Zoom Out Reset image size

Furthermore, it is analysed the fundamental mode $\omega_0(k)$, since it is the only mode which contributes for a particle travelling on axis (or if the wakefield is calculated on axis). Figure 3 depicts the dependence of the fundamental mode $\omega_0(k)$ of the resonant frequency on the wavenumber k for a CNT with a radius a = 1 nm as well as the contribution of the three addends in equation (12). It can be observed that the contribution from the $\Omega_p$ term dominates for low values of k, whereas the contribution from the β term is dominant for larger values of k. Nevertheless, there is an intermediate region, $10^{10}\,\textrm{m}^{-1}\lesssim k\lesssim 10^{11}\,\textrm{m}^{-1}$, where the three contributions exhibit considerable similarity. The resonance condition k_m for a given velocity v is given by the intercept of (kv)² (a parallel line to the α contribution) with $\omega_m^2(k)$. Hence, the resonance condition cannot be satisfied if $v\lt\alpha$, where $\alpha = v_F/\sqrt{2}$. In general, two resonances k_m can exist (as seen in figure 3 for $v = 0.1c$) and then the contribution from both resonances must be summed in equations (15) and (16) (cf appendix B). However, the contribution from the resonance with a larger value of k_m is, in general, totally negligible taking into account the exponential behaviour of the Bessel function $K_m(x) = \sqrt{\frac{\pi}{2x}}e^{-x}$ for $x\rightarrow \infty$ and the factor $\left|\frac{\partial Z_m}{\partial k}\right|_{k = k_m}^{-1}$ that decreases rapidly for high values of k_m . Therefore, by resonance condition we practically mean the first root k_m (the resonances shown in figure 2). Consequently, for high velocities v, the resonance condition k_m is obtained (if it is satisfied) in the region where the $\Omega_p$ term dominates and contribution from the parameters α and β is negligible. Thus, the parameters α and β only become important if the resonance condition is obtained in the intermediate region.

Zoom In Zoom Out Reset image size

3.2. Electric wakefields

In the following calculations, it is considered that the point-like charged particle is a proton, i.e. Q = e. Figure 4 shows the three different contributions to the longitudinal wakefield (cf equation (13)) for different driving velocities. It can be seen that the Coulomb term is only important near the driving particle, whereas $W_{z1}$ and $W_{z2}$ are responsible for the plasmonic excitations and are practically identical, except in the proximity of the driving particle. Moreover, the wavelength of the wakefield increases with the velocity v as it was deduced from the dispersion relation. Besides, the wakefield amplitude decreases as the proton speed increases since we are approaching to the velocities that do not satisfy the resonance condition. The friction parameter γ produces an exponential decay of the wakefield ($W_{z1}+W_{z2}$) with the distance behind the driving charged particle. For this reason, the value of γ has been diminished in figures 4(c) and (d) in order to see the plasmonic excitations. Thus, when the friction parameter γ is very small, $W_{z2}$ follows a cosine pattern and the longitudinal wakefield can be approximated by $W_z\approx W_{z0}+2W_{z2}$ ($W_{z2}$ calculated using equation (15)), as shown in figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 depicts an example of the longitudinal and transverse wakefield inside a CNT. Here, to better appreciate the wakefield details, only the contribution from the plasmonic excitation is shown, taking into account that the Coulomb contribution is negligible except near the driving particle. In figure 6 one remarkable thing is that plasmonic wakefields increase with the radial distance r because of the dependence on $I_m(|k|r)$ and $I_m^{\prime}(|k|r)$ (both increasing functions with the argument) of W_z and W_r , respectively. A similar reasoning can be used to show that the wakefields are higher if the particle travels off axis ($r_0\neq0$), although in this case higher order modes ($|m|\gt0$) should be computed. These effects become more important for lower velocities, since the value of the resonance wavenumber increases.

Zoom In Zoom Out Reset image size

Furthermore, it is interesting to note that there is a phase offset of $\pi/2$ between the longitudinal and transverse wakefields. These results agree with the Panofsky-Wenzel theorem [27] and are similar to what is observed for wakefields excited in homogeneous plasmas in the linear regime [28]. As a consequence, there are periodical regions where the witness charged particles can simultaneously experience both acceleration and focusing (if they travel off-axis), as it is depicted in figure 7.

Zoom In Zoom Out Reset image size

3.3. Optimization of the CNT parameters

As pointed out in the previous section, the plasmonic excitations can be approximated by $2W_{z2}$ when the friction parameter γ converges to zero in the considered system. Therefore, equation (15) can be used to efficiently optimize key parameters, such as n₀, a and v, to enhance the longitudinal wakefield amplitude. Concretely, this section is focused on the plasmonic excitation created on axis by a proton travelling on axis, i.e. $r = r_0 = 0$. Figure 8(a) depicts the amplitude of $W_{z2}$ as a function of the radius a. The maximum of the wakefield increases with the surface density n₀ and moves to smaller radii a. It is also worth noting that these plots do not depend on the surface density if both the radius and the wakefield amplitude are normalized to the plasma wavelength $\lambda_p = 2\pi c/\Omega_p$ and the peak maximum $W_{z2}^{\textrm{max}}$, respectively (see figure 8(b)). Thus, there is an optimum radius (in units of λ_p ) for a given driving velocity regardless of the surface density n₀. This optimum radius $(a/\lambda_p)^{\textrm{max}}$ is proportional to the driving velocity v, as it can be seen in figure 9(a). On the other hand, figure 9(b) shows the peak maximum $W_{z2}^{\textrm{max}}$ as a function of v for different surface densities. The maximum wakefield is obtained for lower velocities (as long as they satisfy the resonance condition) and increases with the surface density. Thus, driving particles with low velocities can excite more efficiently plasmonic modes in CNTs. However, at very low velocities the resonance condition cannot be satisfied. As a result, the curves in figure 9(b) do not originate at v = 0.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In summary, for a given velocity, the maximum longitudinal wakefield is obtained for the optimum CNT radius given by figure 9(a) and a surface density as high as possible. Finally, it is important to remark that for high velocities we should take into account the relativistic effects which are not considered in this manuscript. Nevertheless, it is worth mentioning that a behavior qualitatively similar to that shown in figure 8(b) was also observed in figure 4 of [9], where particle-in-cell (PIC) simulations were used with a driving bunch composed of 1 GeV electrons (with an energy spread of 1$\%$). In that case, the maximum wakefield was obtained for $a/\lambda_p\approx0.10$, a value quite similar to what was obtained in figure 9(a) as $v\rightarrow c$. Therefore, the CNT radius optimization carried out in this manuscript might provide a good approximation of its value.

The linearized hydrodynamic model in conjunction with the Poisson's equation has been used to study the electric wakefields generated by a point-like charge travelling parallel to the axis in a CNT. General expressions have been derived for the longitudinal and transverse wakefields and their dependencies on the surface density, the CNT radius and the velocity of the driving charged particle have been numerically studied and related to the dispersion relation. It has been shown that the friction parameter produces an exponential decay of the excited plasmonic modes. If the friction is negligible, the plasmonic excitations can be approximated by twice the equations (15) and (16). This approximation for the longitudinal wakefield was used to perform an optimization of the CNT radius (in units of the plasma wavelength) for a given driving velocity. Interestingly, at least at a qualitative level, the results agree with those obtained in [9] through PIC models. Hence, the linearized hydrodynamic model might be used to obtain an approximation of the optimum radius without requiring time-consuming PIC simulations. While in the hydrodynamic model the CNT wall is a cylindrical shell, and the driver is a point particle, PIC simulations in a 3D (or quasi-3D) geometry require a volumetric region populated with particles, and the beam with a finite size. Hence, as one tries to set the PIC simulation aiming to match the scenario described by the hydrodynamic model, a thinner wall thickness and a smaller beam (if compared to the CNT internal radius r_in) are required. For a constant volumetric charge density, reducing the CNT wall thickness ($w = r_{\textrm{out}} - r_{\textrm{in}}$) implies reducing the number of available electrons as well. Then, for a constant internal radius r_in, the amplitude of the excited wakefield is damped as $r_{\textrm{out}} \rightarrow r_{\textrm{in}}$. Because there is no clear correspondence between the CNT radius and the surface density in the hydrodynamic model and the tube aperture and volumetric density in the PIC code, approximations and assumptions might be required to properly benchmark the hydrodynamic model with PIC simulation results. In future works, a systematic quantitative comparison will be performed to analyze similarities, differences and limitations of both methods.

It is worth mentioning that the wakefields generated by a point-like charge could be used as a Green's function to compute the wakefields excited by a driving bunch with an arbitrary charge distribution. On the one hand, it has been shown that a single proton may excite GV m⁻¹ wakefields in CNTs with nanometric radii. However, if we consider a beam with a certain number of protons with low velocities (provided that the beam size is small enough to simply consider the sum of the contribution of each single particle), the linear approximation may not be valid (see appendix A). Furthermore, the obtention of beams with such small sizes can prove to be extremely challenging. On the other hand, an ultra-relativistic proton may excite wakefields on the order of 1 MV m⁻¹ in CNTs with radii ∼100 nm. Consequently, an ultra-relativistic, high-charge proton beam may be able to excite ∼1 TV m⁻¹ wakefields (i.e. more than three orders of magnitude higher than those obtained with conventional RF cavities), as long as its energy is enough to prevent Coulombian explosion from space-charge effects. As the optimal radius and wavelength associated with these wakefields are estimated to be $\gtrsim$100 nm, it would be also easier to obtain beams which can be propagated inside the CNT and witness beams which fit in the periodic regions where they can experience acceleration and focusing for ultra-relativistic driving velocities. Moreover, ultra-relativistic driver particles would allow to mitigate space-charge effects and avoid the dephasing between the witness beam and the wakefield excited by the driver during the acceleration. Hence, the use of nano-structures may open new possibilities to obtain ultra-high particle acceleration gradients.

This work has been supported by Ministerio de Universidades (Gobierno de España) under Grant Agreement FPU20/04958, and the Generalitat Valenciana under Grant Agreement CIDEGENT/2019/058.

All data that support the findings of this study are included within the article (and any supplementary files).

The linearized hydrodynamic model is based on the assumption that the perturbed quantities n₁ and u are very small (cf equations (1) and (3)). In this appendix we are going to study if the linear condition $n_1 \ll n_0$ is fulfilled throughout the results shown in the manuscript. The perturbed density n₁ can be obtained as a Fourier-Bessel expansion as [15]:

$$\begin{equation} n_{1}\left(r, \varphi, \zeta\right) = \frac{Q n_0 e}{4 \pi^2 \varepsilon_0 m_e} \sum_{m = -\infty}^{+\infty} \int_{-\infty}^{+\infty} \mathrm{d} k \,e^{i k \zeta+i m\left(\varphi-\varphi_0\right)} I_m\left(|k| r_0\right) C_m\left(k\right), \end{equation} \tag{ 17 }$$

where

$$\begin{equation} C_m\left(k\right) = -\frac{\left(k^2+m^2 / a^2\right) K_m\left(|k| a\right)}{D_m\left(k\right)}. \end{equation} \tag{ 18 }$$

For instance, the perturbed densities associated with the wakefields depicted in figure 4 are shown in figure 10. It is shown that the perturbed density satisfies the condition $n_1\lt2\times10^{-3} n_0$ in all the cases (i.e. the linear approximation is fulfilled). Note that, as the driver is travelling on-axis, the perturbed density just depends on the comoving coordinate ζ.

Zoom In Zoom Out Reset image size

Evidently, if the perturbed density increases, the electric potential and the excited wakefields are higher. Thus, figure 11 depicts the amplitude of the perturbed density oscillation for the optimum radii shown in figure 9(a). It can be seen that the linear approximation is fulfilled except for the lower velocities when $n_1 \lesssim n_0$. Nevertheless, it is known that the analytical expressions that describe the beam-driven wakefields in homogeneous plasmas in the linear regime hold reasonably well even in cases of nonlinear regime [29]. Therefore, the linear assumption's range of validity for deriving approximations of the excited wakefields might be more extensive, with the linear regime serving, at the very least, as a predictive tool. Consequently, we can conclude that the linear approximation may be a good prediction even for the results obtained with the lower velocities. As the perturbed density n₁ is proportional to the driving charge Q, the main limitation of the linear model will appear if we consider a driving charge representing a macroparticle with the charge of a certain number of protons (or electrons), since in that case the perturbed density n₁ can become much larger than n₀. However, it is worth noting that the linear approximation will be valid even for an ultra-relativistic driver with the charge of 10⁶ protons if $n_0 = n_g$.

Zoom In Zoom Out Reset image size

We are going to demonstrate the expression for $W_{z2}$, since the demonstration for $W_{r2}$ is totally analogous but exchanging $I_m(k_mr) \rightarrow I_m^{^{\prime}}(k_mr)$ and $\cos(k_m\zeta) \rightarrow \sin(k_m\zeta)$. The resolution of the integral is as follows:

$$\begin{equation} \begin{aligned} W_{z2}& = \lim _{\gamma \rightarrow 0^+}\frac{Q}{4 \pi^2 \varepsilon_0} \sum_{m = -\infty}^{+\infty} \int_{-\infty}^{+\infty} \mathrm{d} k\, k e^{i m\left(\varphi-\varphi_0\right)} I_m\left(|k| r_0\right) I_m\left(|k| r\right)\textrm{Im}\left[A_m\left(k\right)\right] \cos \left(k \zeta\right) \\ & = \lim _{\gamma \rightarrow 0^+}\frac{Q}{2 \pi^2 \varepsilon_0} \sum_{m = -\infty}^{+\infty} \int_{0}^{+\infty} \mathrm{d} k\, k e^{i m\left(\varphi-\varphi_0\right)} I_m\left(k r_0\right) I_m\left(k r\right)\textrm{Im}\left[A_m\left(k\right)\right] \cos \left(k \zeta\right) \\ & = \lim _{\gamma \rightarrow 0^+}\frac{-Q}{2 \pi^2 \varepsilon_0} \sum_{m = -\infty}^{+\infty} \int_{0}^{+\infty} \mathrm{d} k\, k e^{i m\left(\varphi-\varphi_0\right)} I_m\left(k r_0\right) I_m\left(k r\right) \Omega_p^2 a^2\left(k^2+m^2 / a^2\right) K_m^2\left(k a\right) \\ &\quad \times\frac{\textrm{Im}\left[D_m\left(k\right)\right]}{\left(\textrm{Re}\left[D_m\left(k\right)\right]\right)^2+\left(\textrm{Im}\left[D_m\left(k\right)\right]\right)^2} \cos \left(k \zeta\right) \\ & = \frac{-Q}{2 \pi \varepsilon_0} \sum_{m = -\infty}^{+\infty} \int_{0}^{+\infty} \mathrm{d} k\, k e^{i m\left(\varphi-\varphi_0\right)} I_m\left(k r_0\right) I_m\left(k r\right) \Omega_p^2 a^2\left(k^2+m^2 / a^2\right) K_m^2(k a) \\ &\quad \times\lim _{\textrm{Im}\left[D_m(k)\right] \rightarrow 0^+}\frac{1}{\pi} \frac{\textrm{Im}\left[D_m(k)\right]}{(\textrm{Re}\left[D_m(k)\right])^2+(\textrm{Im}\left[D_m(k)\right])^2} \cos (k \zeta) \\ & = \frac{-Q}{2 \pi \varepsilon_0} \sum_{m = -\infty}^{+\infty} \int_{0}^{+\infty} \mathrm{d} k\, k e^{i m\left(\varphi-\varphi_0\right)} I_m\left(k r_0\right) I_m(k r) \Omega_p^2 a^2\left(k^2+m^2 / a^2\right) K_m^2(k a) \\ &\quad \times \delta((\textrm{Re}\left[D_m(k)\right])) \cos (k \zeta) \\ & = \frac{-Q}{2 \pi \varepsilon_0} \sum_{m = -\infty}^{+\infty} e^{i m\left(\varphi-\varphi_0\right)} \sum_{n} k_{m,n} I_m\left(k_{m,n} r_0\right) I_m\left(k_{m,n} r\right) \Omega_p^2 a^2\left(k_{m,n}^2+\frac{m^2}{a^2}\right) \\ &\quad \times K_m^2\left(k_{m,n} a\right)\left|\frac{\partial Z_m}{\partial k}\right|_{k = k_{m,n}}^{-1} \cos \left(k_{m,n} \zeta\right), \end{aligned} \end{equation} \tag{ 19 }$$

where we have used that in the limit $\gamma \rightarrow 0^+$ the function to integrate is an even function, the relation $\textrm{Im}\left[D_m(k)\right] = kv\gamma\propto\gamma$ and the following properties of the Dirac delta:

$$\begin{equation} \delta\left(x\right) = \lim _{y \rightarrow 0^+} \frac{1}{\pi}\left(\frac{y}{x^2+y^2}\right), \end{equation} \tag{ 20 }$$

$$\begin{equation} \delta\left(f\left(x\right)\right) = \sum_n\left|f^{^{\prime}}\left(x_n\right)\right|^{-1} \delta\left(x-x_n\right), \quad \operatorname{with} \quad f\left(x_n\right) = 0, \quad f^{^{\prime}}\left(x_n\right) \neq 0. \end{equation} \tag{ 21 }$$

In equation (19), $k_{m,n}$ are the positive roots of $Z_m(k) = \textrm{Re}[D_m(k)] = (k v)^2-\omega_m^2(k)$. As it is explained in section 3.1 (cf figure 3), the function $Z_m(k)$ has a maximum of two roots, but the contribution from the root with a large value is totally negligible. Thus, equation (19) becomes equation (15) if only the first root is considered.