The excitation of Terahertz (THz) radiation by the interaction of an ultrashort laser pulse with the modes of a miniature corrugated plasma waveguide is considered. The axially corrugated waveguide supports the electromagnetic modes with appropriate polarization and subluminal phase velocities that can be phase matched to the ponderomotive potential associated with the laser pulse, making significant THz generation possible. This process is studied via full format Particle-in-Cell simulations that, for the first time, model the nonlinear dynamics of the plasma and the self-consistent evolution of the laser pulse in the case where the laser pulse energy is entirely depleted. It is found that the generated THz is characterized by lateral emission from the channel, with a spectrum that may be narrow or broad depending on the laser intensity. A range of realistic laser pulse and plasma parameters is considered with the goal of maximizing the conversion efficiency of optical energy to THz radiation. As an example, a fixed drive pulse (0.55 J) with a spot size of 15 *μ*m and a duration of 15 *fs* produces a THz radiation of 37.8 mJ of in a 1.5 cm corrugated plasma waveguide with an on axis average density of 1.4 × 10^{18 }cm^{−3}.

## I. INTRODUCTION

Terahertz radiation (THz) lies between microwave and infrared wavelengths in the electromagnetic (EM) spectrum and spans frequencies from 300 GHz to 20 THz. A wide variety of applications for THz radiation^{1} can be found including time domain spectroscopy (TDS),^{2} medical and biological imaging,^{3,4} and remote detection^{5–7} among others. Most airports, for instance, use millimeter wave/THz scanners for security checking. Existing small scale THz sources based on laser-solid interactions are limited to *μJ*/pulse levels due to material damage^{8} although a recent discovery using optical rectification (OR)^{9–11} in organic crystals can exceed this limit. The need for higher power sources has led to the consideration of THz generation via laser-plasma interactions,^{12–14} in which THz peak energies of tens of *μJ* can be achieved. Higher energy THz pulses can be generated at accelerator facilities via synchrotron^{15} or transition radiation.^{16} Such facilities are relatively large and expensive to operate. This motivates research interest in small-scale, high efficiency terahertz sources.

Research has been actively conducted to investigate THz radiation generation by laser pulses propagating in plasmas since it was first demonstrated by Hamster *et al.*^{17,18} In this case, the source of the radiation is the current driven by the ponderomotive force of a laser pulse. However, generation of radiation by a laser pulse propagating through a uniform plasma is generally minimal. In order for electromagnetic modes to efficiently couple to the driving source, which travels at the group velocity of the laser pulse, the plasma must be inhomogeneous or immersed in a strong background magnetic field. A THz generation scheme involving laser pulses (or possibly electron beams) propagating through axially corrugated plasma channels has been proposed by Antonsen *et al.*^{19,20} In this scheme, the corrugated plasma channel acts as a slow wave structure that supports electromagnetic modes that can be phase matched with the driver. The scheme offers the possibility of much higher efficiency of THz generation under conditions of full laser pulse depletion. This prospect is investigated here.

In this paper, we report theoretical, numerical, and simulation results of ponderomotively driven THz generation by laser pulses propagating through corrugated plasma waveguides. Such waveguides have already been realized in the laboratory.^{21,22} The experimental set up is shown in Fig. 1(a). A Nd:YAG laser pulse is line-focused onto a cluster jet creating a plasma that hydrodynamically expands, leading to the formation of a channel with a transversely parabolic density profile. The periodic structure is created by spatially modulating the laser intensity on the back side of the axicon or by periodically obstructing the cluster flow. A second ultra-short Ti:Sapphire laser pulse is then injected into the channel following the channel formation pulse. The second pulse drives the terahertz generation. Shown in Fig. 1(b) is a snap shot of the experimentally generated axially modulated plasma density profile.

The organization of this paper is as follows: In Sec. II, we introduce the mechanisms underlying THz generation in corrugated plasma waveguides. A mathematical model of the channel is considered, and the dispersion relation is analyzed. In Section III, we conduct full format Particle-in-Cell (PIC) simulations to investigate the self-consistent evolution of the laser pulse in the strongly nonlinear regime and optimize the conversion efficiency from the laser pulse to THz energy. Tunability of the THz spectrum is also discussed. A numerical method to approximately find the frequencies of modes and coupling to the vacuum is also provided. The dependence of THz radiation on plasma density, driver pulse intensity, laser pulse duration, channel length, and other channel parameters is investigated and discussed in detail. Two different types of plasma channels are considered and compared. In Section IV, we present our conclusions and discuss future directions.

## II. EXCITATION OF THZ MODES IN CORRUGATED PLASMA CHANNELS

### A. Ponderomotive driver

We consider the interaction of a drive pulse with a tenuous plasma. The central frequency of the drive pulse is significantly higher than the generated electromagnetic waves. The pulse produces a cycle-averaged low frequency ponderomotive force on the electrons, which induces an electron current that can produce radiation. The force on electrons is the gradient of the ponderomotive potential of the pulse, $ F p = \u2212 \u2207 V p $, where *V _{p}* =

*mc*

^{2}

*a*

^{2}/4, the normalized laser vector potential is

*a*=

*eE*/

_{a}*mcω*

_{0},

*m*is the electron mass,

*e*is the electron charge,

*c*is the speed of light in vacuum,

*ω*

_{0}is the laser carrier frequency, and

*E*is the electric field amplitude with both spatial and temporal dependences. The plasma ions are relatively heavy and can be treated as a stationary background during the interaction. In order for the excited current to generate electromagnetic radiation, the plasma must be inhomogeneous so that in no reference frame, the driven currents are static.

_{a}Further to have significant THz emission, the phase velocity of the excited channel modes must match the group velocity of the laser pulse in the plasma. In a uniform plasma, the modes typically have superluminal phase velocities according to the dispersion relation $ \omega 2 = k 2 c 2 + \omega p 2 $, where $ \omega p = 4 \pi e 2 N / m $ is the plasma frequency. However, the phase velocity in a corrugated plasma channel can be subluminal, allowing phase matching. In particular, the periodic spatial modulations of the channel parameters cause the electromagnetic modes of the channel to have a Floquet type representation, which implies subluminal partial waves.^{23} Because of the superposition of spatial harmonics, these EM modes can have components with subluminal phase velocities. Thus, significant THz emission can be achieved.

### B. Corrugated plasma channels

We consider the corrugated plasma waveguides to be cylindrically symmetric, with electron densities *N*(*r*, *z*) described by the following:

where *N*_{0} is a normalization density. The quantities *n*_{0}, *n*_{1}, *r*_{0}, and *r _{c}* are all potentially periodic functions of

*z*with period

*λ*and the modulation wavenumber

_{m}*k*= 2

_{m}*π*/

*λ*. The quantity

_{m}*n*

_{0}(

*z*) is the normalized on-axis density, and $ n 1 ( z ) $ is the normalized density at

*r*=

*r*. For our studies, we take $ n 0 = 1 + \u2009 \delta \u2009 sin ( k m z ) $ and $ n 1 = n \xaf 1 + \delta 1 \u2009 sin ( k m z ) $, respectively,

_{c}*δ*is the density modulation amplitude of the on-axis density

*n*

_{0}, $ n \xaf 1 $ is the average density at

*r*=

*r*, and

_{c}*δ*

_{1}is the density modulation amplitude of

*n*

_{1}. The density has a parabolic transverse profile to guide the laser pulse during propagation. The quantity

*r*is the radius at which the density reaches

_{c}*n*

_{1}, and then the density decreases linearly to zero from

*r*to

_{c}*r*

_{0}. The quantities

*r*and

_{c}*r*

_{0}may also be axially modulated, $ r c , 0 = r \xaf c , 0 + \Delta c , 0 \u2009 cos ( k m z + \theta c , 0 ) $. Figures 2(a) and 2(b) are false color images of two different density profiles, based on parameters given in the caption. Both profiles are similar to those realized experimentally.

^{21,22}The channel in Fig. 2(a) has its peak density off axis, while the channel in Fig. 2(b) has its peak density on axis. We refer to these two types as the “on-axis peak” (Fig. 2(b)) and the “off-axis peak” (Fig. 2(a)) channels, respectively.

### C. Dispersion relation and mode excitation

We review the linear theory of channel modes presented in Ref. 23. The plasma is taken to be a cold fluid with a linear response. Radially polarized, azimuthally symmetric TM modes (*E _{r}*,

*E*, and

_{z}*B*) of the channel are considered. In this case, an approximate wave equation assuming $ | \u2207 \xb7 E | \u226a | E | / r c $ (Refs. 19 and 23) for the radial electric field

_{θ}*E*can be derived

_{r}where *ω _{p}*

_{0}is the plasma frequency evaluated for the normalization density

*N*

_{0}.

In the case in which the electron density profile is parabolic as *r* → *∞*, i.e., it is characterized by the first line of Eq. (1), and we take *δ*_{1} = *δ*, we can find analytic expressions for the eigenmodes. Specifically, the *γ*th order radial eigenmode of the channel is described by

where *w _{ch}* is the mode width given by $ 4 / w c h 4 = \omega p 0 2 ( n 1 \u2212 n 0 ) / r c 2 c 2 $. The function

*H*is the

_{γ}*γ*th polynomial defined as $ \u2211 n = 1 n = 2 \gamma \u2212 1 \alpha n ( r / w c h ) n $ with coefficients determined by $ \alpha n / \alpha n \u2212 2 = 4 ( n \u2212 1 \u2212 2 \gamma ) / ( n 2 \u2212 1 ) $ and

*α*

_{1}= 1. Given the dependence of the fields on the radial coordinates, the function

*f*(

*z*) satisfies the following Mathieu equation:

where the value of *k*_{0} is related to the frequency by

The dispersion relation is found by solving Eq. (4) with Floquet boundary conditions, $ f ( z + \lambda m ) = exp ( i k z \lambda m ) f ( z ) $. This then determines the dependence *k*_{0}(*k _{z}*), which is inserted in Eq. (5) and gives

*ω*(

*k*) the dispersion relation.

_{z}Shown in Fig. 3 is the dispersion relation for the lowest (*γ* = 1) radial mode of the model channel with parameters given in the caption. The dependence of *ω* on phase advance demonstrates the characteristic periodicity of frequencies in *k _{z}* for periodic structures. The laser pulse, represented by a straight line in the plot, moves at its group velocity (

*v*≃

_{g}*c*) in the plasma channel. At places where the dashed pulse line and the dispersion relation curve intersect, phase matching occurs, and THz excitation can be expected at these frequencies.

Equations (4) and (5) apply as long as *w _{ch}* <

*r*, that is, the THz or laser spot size is smaller than the channel width. Equation (5) is a good approximation to the dispersion relation regarding the cycle-averaged density profile and assuming that the transverse parabolic shape extends to infinity. However, to determine the exact frequencies of the excited modes, one can numerically evaluate the wave equation, Eq. (2), using an exact electron density profile. A more accurate calculation of the dispersion relation, including the comparison with Eq. (5), is provided in the Appendix.

_{c}## III. SIMULATION RESULTS

### A. THz mode excitation

THz generation in corrugated plasma waveguides is simulated using the full format PIC code TurboWAVE.^{24} The simulations, performed in 2D planar geometry, feature a finite sized plasma channel illuminated by an ultrashort, intense laser pulse incident from the simulation boundary. Figure 4(a) shows an example of an off-axis density peak plasma channel of 10 periods with a modulation wavelength of *λ _{m}* = 50

*μ*m. To quantify the radiation emitted from the plasma channel, we calculate the Poynting flux and its spectral density through each of the surfaces outside the plasma region indicated by the dashed lines. This captures emission in the forward, backward, and lateral directions. The initial 2D simulations are performed in the lab frame in a domain of dimensions 205.9 × 753.8

*μ*m with 1024 × 20 480 cells in the

*x*and

*z*directions, respectively. A laser pulse, with parameters detailed in Table I, traverses the plasma channel from left to right. The initial normalized vector potential is

*a*

_{0}= 0.4 (pulse energy

*U*= 66 mJ) with

_{L}*a*

_{0}defined as

*a*

_{0}=

*eE*

_{0}/

*mcω*

_{0}, where

*E*

_{0}is the peak field amplitude. Figure 4(b) shows a false color image of the transverse component of Poynting flux

*P*∝

_{x}*E*after the laser pulse traverses the channel. The excited plasma wave can be observed as the rapid oscillations of

_{z}B_{y}*P*inside the channel; the alternate positive and negative values indicate a small average flux. However, at both lateral boundaries (top and bottom of the image), one can observe the THz emission in the form of the red and blue streaks in the lateral Poynting flux.

_{x}Laser pulse . | Central wavenumber λ 800 nm
. | Spot size r 15 _{L}μm
. | Pulse duration τ, FWHM 50 fs
. | Normalized vector potential a_{0} 0.4
. | ||||||
---|---|---|---|---|---|---|---|---|---|---|

Channel type | λ[_{m}μm] | $ r \xaf c [ \mu m ] $ | $ r \xaf 0 [ \mu m ] $ | δ | $ n \xaf 1 $ | δ_{1} | Δ_{c} | θ _{c} | Δ_{0} | θ_{0} |

Off-axis peak (Fig. 2(a)) | 15 | 30 | 40 | 0.9 | 3 | 0.9 | 0 | 0 | 0 | 0 |

On-axis peak (Fig. 2(b)) | 15 | 22.5 | 37.5 | 0.7 | 1.3 | 0.1 | 7.5 | π/2 | 7.5 | π/2 |

Laser pulse . | Central wavenumber λ 800 nm
. | Spot size r 15 _{L}μm
. | Pulse duration τ, FWHM 50 fs
. | Normalized vector potential a_{0} 0.4
. | ||||||
---|---|---|---|---|---|---|---|---|---|---|

Channel type | λ[_{m}μm] | $ r \xaf c [ \mu m ] $ | $ r \xaf 0 [ \mu m ] $ | δ | $ n \xaf 1 $ | δ_{1} | Δ_{c} | θ _{c} | Δ_{0} | θ_{0} |

Off-axis peak (Fig. 2(a)) | 15 | 30 | 40 | 0.9 | 3 | 0.9 | 0 | 0 | 0 | 0 |

On-axis peak (Fig. 2(b)) | 15 | 22.5 | 37.5 | 0.7 | 1.3 | 0.1 | 7.5 | π/2 | 7.5 | π/2 |

^{a}

Note: Unless otherwise stated, all parameters are identical to those in Table I.

To investigate the spectrum of the lateral THz emission, the fields are Fourier transformed in time. The radiated energy per unit length $ U \u2032 $ through each diagnostic surface in Fig. 4(a) can be obtained in the form of a spectral density in $ ( \omega , r ) $, where ** r** denotes the 2D spatial coordinates (

*x*,

*z*)

where $ n \u0302 $ is the unit vector normal to the surface *A*. The spectral density is given by

To quantify the THz radiation emitted across each diagnostic surface, we calculate the *z* component of the spectral density, *S _{z}*, for the left and right diagnostic boundaries and the

*x*component,

*S*, for the lateral boundary. Figure 5(b) is the laterally radiated spectral density,

_{x}*S*, from a PIC simulation of the off-axis type density profile. The low frequency, broad band THz radiation observed at the entrance of the channel, shown in Fig. 5(a), is due to resonant transition radiation as discussed in Refs. 25–27. Lateral THz radiation

_{x}^{28}is also observed in the corrugated plasma channel and characterized by a coherent, narrow band spectrum as shown in Fig. 5(c). The frequencies of the first three excited THz modes based on the phase matching condition predicted by the simplified model of Sec. II C are 13.18 THz, 15.3 THz, and 16.5 THz, respectively. In this case, the simulation results in Fig. 5(b) show that the THz radiation leaves the channel, creating an intensity pattern with maxima at 9 separate locations along the longitudinal distance where the fundamental mode emission dominates.

### B. Dependence on plasma density

We now explore the dependence of THz generation on plasma density in the case of the off-axis peak channel of Fig. 2(a). Figure 6 displays the radiated spectral density $ d U \u2032 / d \omega $, for several average, on-axis plasma densities. The central frequency of the emission peak increases with increasing density, in accord with the dependence of the channel mode frequencies on density. For all the plasma densities considered in Fig. 6, the fundamental mode provides the dominant contribution to the radiation energy. As shown in Fig. 6, the radiation generated for this channel peaks around a plasma density of 1.75 × 10^{17 }cm^{−3}. As discussed in Ref. 19, the energy extracted from the drive pulse is converted into both electromagnetic radiation (EM) and plasma waves (PW). For example, the simulation predicts that the THz energy radiated in a 1.5 cm channel for the 1.4 × 10^{17 }cm^{−3} density is 0.048 mJ, while the energy extracted from the laser pulse over the same distance is 4.1 mJ (∼1.17% of depletion). To efficiently deplete the laser pulse energy within a shorter channel, a higher plasma density is preferred. Therefore, to efficiently generate THz radiation, there must be a trade-off between the increasing density to increase the plasma current and the decreasing density to increase the lateral output coupling as shown in Fig. 6.

### C. Dependence on laser intensity

The spatial and spectral dependences of the laterally radiated THz are illustrated for different laser intensities in Fig. 7. Ponderomotively driven radiation is expected to scale with the laser amplitude as $ U L 2 \u223c a 0 4 $ (quadratic in ponderomotive potential) for *a*_{0} ≪ 1. One can see that the THz energy follows this scaling for small *a*_{0} but is enhanced above this scaling for *a*_{0} > 1. Because the relativistic ponderomotive potential scales as $ a 0 2 / \gamma $, the observed enhancement is even stronger than might be expected. The enhancement phenomenon is accompanied by a broadening of the THz spectrum to the point that individual modes can no longer be identified. For large *a*_{0}, we find that the higher frequency channel modes are excited by nonlinear currents. This, along with the broadening of the spectrum, provides the enhancement over the linear scaling. For this short channel (10 periods), the modification of the temporal profile of the laser pulse is small. Thus, the broad spectrum can be attributed to nonlinearity in the excited plasma response.

Our goal is to optimize the conversion efficiency of optical laser pulse energy to THz. Channel lengths in our simulations are limited by computation time, and so, we are not able to simulate every channel for a distance long enough to substantially deplete the laser pulse. We thus first consider short channels (10 periods) and define an efficiency that is the fraction of the depleted laser energy ( $ | \Delta E L a s e r | $) transferred to THz, $ \eta = E T H z / | \Delta E L a s e r | $. By maximizing this efficiency, less power is expended driving the plasma oscillations, thus freeing it to drive THz over longer distances. We note that this efficiency does not depend on the laser intensity for *a*_{0} < 1, based on the linear theory for which both *E _{THz}* and Δ

*E*scale as $ a 0 4 $. However, for large

_{Laser}*a*

_{0}, higher frequency THz modes are excited by nonlinear currents, enhancing the efficiency scaling. We have achieved an energy conversion efficiency of approximately 3% for the THz generation as displayed in Fig. 7(b), and this can be further optimized by varying the corrugated plasma density profiles. For example, the conversion efficiency for the weakly relativistic case can be efficiently enhanced by finding an optimum plasma density as discussed in regard to Fig. 6. Further, as the laser pulse propagates through the channel, its energy depletes. The accompanying spectral red-shifting results in compression and increases the pulse amplitude,

^{29,30}which contributes to the enhancement of conversion efficiency. This phenomenon will be discussed in detail in Sec. III E.

### D. Dependence on laser pulse duration

The frequencies of the linear channel modes are determined by the plasma profile. Thus, the pulse duration does not affect the THz mode frequency, but it does determine the amplitude of the driving current at each frequency. THz emission can be expected for frequencies given by the intersections of the channel dispersion curves and the laser pulse “light line.” An additional requirement for the generation of THz is that the temporal spectrum of the laser pulse envelope includes the mode frequency. If the laser pulse has a Gaussian temporal profile $ exp [ \u2212 ( t \u2212 z / c ) 2 / \tau p 2 ] $ with $ 2 l n ( 2 ) \tau p $ as the pulse duration (FWHM), the amplitude of the ponderomotive driver (for fixed *a*_{0}) at a mode frequency *ω* is given by $ \tau p \u2009 exp ( \u2212 \omega 2 \tau p 2 / 4 ) $, and thus, minimal radiation is expected for frequencies *ωτ _{p}* ≫ 1. Therefore, the value of

*τ*can be adjusted to excite a specific range of channel mode frequencies. For excitation of the fundamental mode of 13.18 THz, the desired value is

_{p}*τ*= 16.9 fs corresponding to a FWHM pulse duration of 28.4 fs. The simulation results further verify our estimation. Shown in Fig. 8 is a comparison of the radiated spectral density $ d U \u2032 / d \omega $ for pulse durations of 100 fs, 30 fs, and 15 fs, respectively. The initial normalized vector potential

_{p}*a*

_{0}= 0.4 is kept fixed for all 4 cases. For the case of a 100 fs laser pulse, the amplitude of the ponderomotive driver for any channel mode is small, such that minimal THz generation is observed. For the fundamental mode (13.18 THz) of the same plasma channel, the desired pulse duration is 30 fs, and the simulation result shows that the THz radiation maximizes at this frequency. In addition, higher order radiation is observed as the pulse duration is shortened to 15 fs as shown in Fig. 8, where there is an enhancement of the channel mode near 20 THz.

### E. Scaling with channel length

We now investigate for the first time THz generation in corrugated plasma waveguides of sufficient length to deplete the laser pulse. As the laser pulse propagates through the plasma channel, the pulse envelope is modified and its energy is depleted through conversion into both electromagnetic radiation and plasma waves.^{19,31–34} For example, the energy of a laser pulse with *a*_{0} = 2.0 and a pulse duration of 15 fs propagating through the 10 period channel shown in Fig. 4(a) is depleted by 1.12%. At the same time, the central frequency red-shifts from 375 THz to 371 THz. These changes in energy and central frequency are consistent with action conservation:^{29–32} as the pulse depletes, the spectrum red-shifts and the normalized vector potential *a*_{0} associated with the laser pulse increases. As shown in Sec. III C, the pulse depletion rate increases dramatically with intensity.

We thus simulate an 800 nm laser pulse of a duration of 15 fs, a transverse spot size of 15 *μ*m, and *a*_{0} = 2.0 (pulse energy 0.55 J). The channel parameters are the same as in Fig. 4(a), except the channel length extended to 1.5 cm. The simulation is conducted with a moving window of length 750 *μ*m to collect all the THz emission following the pulse. The 2D moving frame has a size of 205.9 × 753.8 *μ*m with 1024 × 20 480 cells in the *x* and *z* directions, respectively. The energy stored in the laser pulse is displayed in Fig. 9 as a function of propagation distance. Within the propagation distance of 1.5 cm, 80% of the pulse energy is depleted. Fig. 9 also displays the radiated THz energy versus propagation distance. The rate ( $ d U \u2032 / d z $) of THz energy generation increases with distance as the normalized vector potential *a*_{0} increases during propagation due to action conservation. As a result, after propagation of 1.5 cm, more than 8% of the total pulse energy is converted into THz radiation.

For comparison, THz generation using a relatively low intensity laser pulse with a pulse duration of 30 fs, a transverse spot size of 15 *μ*m, and *a*_{0} = 0.4 (pulse energy, 44.5 mJ) is also simulated. In this case, we found the optimum plasma density for short propagation distance (10 periods) to be 1.4 × 10^{17 }cm^{−3} as shown in Sec. III B. The lower plasma density and laser intensity lead to a much longer pulse depletion length than for that of the previous case displayed in Fig. 9. Consequently, due to computational restrictions, we were not able to simulate a channel long enough to substantially deplete the pulse energy. Instead, we simulate pulse propagation through a corrugated plasma channel for a distance of 1.5 cm. The simulation results displayed in Fig. 10 indicate that only 10% of the energy stored in the laser pulse is depleted within that distance. The variation with *z* of the depletion rate is due to the mismatch between the transverse pulse width and the matched width for the guiding channel. This leads to variations in spot size and intensity, with the depletion rate depending strongly on intensity. Shown in Fig. 10 is the radiated THz energy versus plasma channel length. About 50 *μ*J of THz energy is generated within a 1.5 cm interaction distance. As a result, it can be concluded that the conversion efficiency for *a*_{0} = 0.4 is much lower than that of a higher intensity pulse. However, a much longer plasma channel is needed to deplete the pulse energy for *a*_{0} = 0.4.

### F. On axis peak density channel

In this section, we consider the corrugated plasma channel of the type shown in Fig. 2(b). The channel length is 10 periods with a modulation wavelength of 50 *μ*m. The other channel parameters are $ \delta = 0.7 , \u2009 n \xaf 1 = 1.3 , \u2009 \delta 1 = 0.1 $ (see Eq. (1)). In order to match the density profile shown in Fig. 1(b), both cut-off radius and channel radius are also modulated according to $ r c [ \mu m ] = 22.5 \u2212 7.5 \u2009 cos ( k m z ) $ and *r*_{0}[*μ*m] = *r _{c}* + 15. Simulation results for laser pulses, with normalized vector potential

*a*

_{0}= 0.4 and

*a*

_{0}= 2.0, are shown in Fig. 11. Both laser pulses have the same transverse spot size of 15

*μ*m and a pulse duration (FWHM) of 50 fs with a central wavenumber of 800 nm. The case with

*a*

_{0}= 0.4 shows that narrow band THz radiation is excited around the fundamental frequency of 13.6 THz. As displayed in Table II, the amount of generated THz energy and conversion efficiency for this channel are significantly higher than the case shown in Fig. 5 for the same laser pulse. This could be explained by the excitation of a higher electron current due to the relatively greater density inhomogeneity experienced by the driver. In addition, the axially averaged density profile has a lower radial barrier that allows the generated THz waves to escape the channel. For a higher intensity laser pulse with

*a*

_{0}= 2.0, the generated THz shown in Fig. 11 is characterized by a different spectrum compared with Fig. 7(a). Although the amount of THz energy is still enhanced relative to

*a*

_{0}= 0.4, the spectrum is confined in a relatively narrow band near the fundamental frequency, while in the case of Fig. 7(a), higher order THz modes are significantly generated and consequently modify the spectrum.

Channel type . | a_{0}
. | Pulse duration [fs] . | On-axis density [×10^{18 }cm^{−3}]
. | Channel length [mm] . | Energy depletion into the channel (%) . | THz energy [mJ] . | efficiency η
. |
---|---|---|---|---|---|---|---|

Fig. 2(a) | 0.4 | 30 | 0.14 | 15 | 10.11 | 0.048 | 1.11% |

0.4 | 50 | 1.4 | 0.5 | 0.0175 | 3.36 e-5 | 0.28% | |

2.0 | 15 | 1.4 | 15 | 81.8 | 38 | 8.44% | |

2.0 | 50 | 1.4 | 0.5 | 0.39 | 0.2 | 3.25% | |

Fig. 2(b) | 0.4 | 30 | 1.4 | 0.5 | 0.033 | 0.0016 | 12.03% |

0.4 | 50 | 1.4 | 0.5 | 0.024 | 0.0013 | 8.17% | |

0.4 | 50 | 0.14 | 0.5 | 0.006 | 4.68 e-4 | 11.7% | |

2.0 | 30 | 1.4 | 15 | 48.2 | 16 | 3% | |

2.0 | 50 | 1.4 | 0.5 | 0.4 | 0.6493 | 9.84% |

Channel type . | a_{0}
. | Pulse duration [fs] . | On-axis density [×10^{18 }cm^{−3}]
. | Channel length [mm] . | Energy depletion into the channel (%) . | THz energy [mJ] . | efficiency η
. |
---|---|---|---|---|---|---|---|

Fig. 2(a) | 0.4 | 30 | 0.14 | 15 | 10.11 | 0.048 | 1.11% |

0.4 | 50 | 1.4 | 0.5 | 0.0175 | 3.36 e-5 | 0.28% | |

2.0 | 15 | 1.4 | 15 | 81.8 | 38 | 8.44% | |

2.0 | 50 | 1.4 | 0.5 | 0.39 | 0.2 | 3.25% | |

Fig. 2(b) | 0.4 | 30 | 1.4 | 0.5 | 0.033 | 0.0016 | 12.03% |

0.4 | 50 | 1.4 | 0.5 | 0.024 | 0.0013 | 8.17% | |

0.4 | 50 | 0.14 | 0.5 | 0.006 | 4.68 e-4 | 11.7% | |

2.0 | 30 | 1.4 | 15 | 48.2 | 16 | 3% | |

2.0 | 50 | 1.4 | 0.5 | 0.4 | 0.6493 | 9.84% |

## IV. CONCLUSIONS AND DISCUSSION

We have investigated THz generation in corrugated plasma channels accounting for nonlinear excitation of plasma waves and laser pulse depletion. Theoretical analysis and Full format PIC simulations were conducted. A range of laser pulse parameters and plasma channel structures were considered with the goal of maximizing the conversion efficiency of optical pulse energy to THz energy.

Table II displays the simulation results for different pulse and plasma parameters for the two types of corrugated channel displayed in Fig. 2. Most of the simulations were conducted using a 10-period channel to investigate the conversion efficiency of the depleted optical pulse energy to generated THz. For these simulations, the pulse energy only depleted a small percentage. Three examples of longer channels were conducted to examine the consequences of significant pulse depletion, and the results are included in Table II. Our general conclusions are as follows: Generally, THz generation increases with laser amplitude *a*_{0}. For fixed *a*_{0}, THz generation at frequency *ω* maximizes, where *ωτ _{p}* ∼ 1, where

*τ*is the pulse duration. This applies for both channel types. For a short channel with only 10 periods (channel length 0.5 mm), simulation results indicate that more THz energy is generated for the on-axis peak density channel type shown in Fig. 2(b). Efficient THz generation involves a trade-off between the increasing density to increase the plasma current and the decreasing density to increase the lateral output coupling. However, since the excited THz mode depends on the channel structure according to the dispersion relation discussed in Sec. II C, lower plasma density also modifies the generated THz spectrum. In addition, to efficiently deplete the laser pulse energy within a shorter channel, a higher plasma density is preferred.

_{p}As an example, we choose a laser pulse with *a*_{0} = 2.0 and a pulse duration of 30 fs. The channel type shown in Fig. 2(b) is used with an averaged on axis electron density of 1.4 × 10^{18 }cm^{−3} and the channel length is extended to 1.5 cm. The simulation results as displayed in Table II show after the laser pulse propagates through the channel for 1.5 cm, around 48% of the laser pulse energy is depleted and 16 mJ of this energy is converted to THz energy with a narrow spectrum (Fig. 11). The conversion efficiency is around 3% and less than the case shown in Sec. III E. This is probably due to the fact that the channel displayed in Fig. 2(b) no longer remains a transverse parabolic structure capable of guiding. In fact, laser energy leaks laterally, and the wave action is not conserved as the laser pulse propagates through the channel. Therefore, the normalized vector potential *a*_{0} decreases with the propagation distance and less THz energy is generated.

In order for the present mechanism to be a useful high power source of THz radiation, the spectrum should be tunable. Since the excited THz consists of a superposition of the channel modes of the corrugated plasma structure, all the parameters used in Eq. (1) can be tuned to modify the THz spectrum. More specifically, varying the averaged on axis density *N*_{0} is the most straightforward way to tune the frequency in the experiment. For example, Figure 6 shows that varying the average density *N*_{0} from 1.75 × 10^{17 }cm^{−3} to 1.4 × 10^{18 }cm^{−3} results in a shift in the central frequency from 5 THz to 14 THz. For the channel type of Fig. 2(a), the THz spectrum changes dramatically with laser intensity going from a narrow spectrum at *a*_{0} = 0.4 to a broad spectrum at *a*_{0} = 2.0 with enhanced energy as well. For the channel type of Fig. 2(b), the spectrum remains strongly peaked over the same range of laser amplitudes. Since the generated THz waves are emitted laterally, one needs an optical system to collect the radiation. This might be realized with a conical mirror to focus the THz radiation to one direction for practical uses. Overall, the mechanism, using realistic corrugated plasma structures, presented in this paper provides a potential high power source of THz with a tunable spectrum and a conversion efficiency of over 8%.

## ACKNOWLEDGMENTS

The authors would like to acknowledge Dr. Daniel Gordon for the use of TurboWAVE and thank Luke Johnson and Thomas Rensink for fruitful discussions. Part of the simulations was performed on NERSC Edison and Deepthought2 clusters at UMD. This work was supported by the Office of Naval Research and the U.S. Department of Energy.

### APPENDIX: CALCULATING FREQUENCIES OF THE RADIAL EIGENMODES

In this appendix, we will discuss the numerical method to calculate the exact frequencies of the axially averaged excited modes. Equation (2) can be further written as

where the cut off wavenumber *k _{c}* of EM modes and plasma wavenumber

*k*are defined as $ k c 2 = \omega 2 / c 2 \u2212 k 0 2 $ and $ k p 2 ( r ) = \omega p 2 ( r ) / c 2 $, respectively. Since we consider radially polarized TM modes,

_{p}*E*must vanish on axis, i.e.,

_{r}*E*(0) = 0. We know that outside the channel,

_{r}*N*(

*r*>

*r*

_{0}) = 0, and thus Eq. (A1) yields Bessel's differential equation. The far field

*E*outside the channel must match the properties of an outgoing wave, which has the form of the first kind Hankel function $ H 1 ( 1 ) ( k c r ) $ and asymptotically behaves as $ E r \u223c 1 r exp ( i k c r ) $. The boundary condition allows us to know the ratio of

_{r}*E*to its derivative outside the channel. As a result, we can numerically integrate Eq. (A1) using the shooting method to determine the

_{r}*k*values that satisfy the on axis boundary condition

_{c}*E*(0) = 0.

_{r}To calculate the radial eigenmodes numerically for an arbitrary, but given the transverse density profile, one can numerically evaluate *E _{r}* by Eq. (A1) using the shooting method for a set of

*k*

_{⊥}and determine what

*k*

_{⊥}satisfies the on axis boundary condition

*E*(0) = 0. For mathematical simplicity, we set $ \Phi = r E r , \u2009 \beta ( r ) = k p 2 ( r ) $ and Eq. (A1) yields

_{r}One can also find that as *r* → 0, Φ ∼ *a* ⋅ *r*^{2} + *b*; the on axis boundary condition is satisfied only if *b* → 0. The finite difference (FD) shooting method can be implemented by

where *h* is the step size.

To find out the number of different radial modes, i.e., *k _{c}*, which can be supported by a channel with finite transverse size, one can scan the parameter

*k*in Eq. (A3) and apply the Nyquist Theory illustrated in Fig. 12.

_{c}*F*(

*s*) is an analytic function defined in a closed region of the complex

*s*-plane shown on the left. As

*s*travels a clockwise path in the

*s*-plane,

*F*(

*s*) encircles the origin on the complex

*F*(

*s*)-plane

*N*times

where *Z* and *P* denote the number of zeros and poles of the function *F*(*s*) in the closed region, respectively. For our shooting method, as shown in Eq. (A3), Φ(0) has no poles and the result yields *N* = *Z*.

This method predicts the frequency of radiation for any given transverse density profile. For example, Fig. 13 shows that for a density profile shown in Fig. 4(a), Φ(0) encircles the origin twice as *k _{c}* scans from 0.301

*μ*m

^{−1}to 0.402

*μ*m

^{−1}, which implies that the channel can support two radial eigenmodes according to Cauchy's principle. Further, one can apply linear interpolation to narrow the range of

*k*for each mode to find the exact value of

_{c}*k*for the field to satisfy the boundary condition. Figures 14(a) and 14(b) are two figures indicating the range of

_{c}*k*during the interpolation to find the exact value

_{c}*k*for first and second radial eigenmodes, respectively.

_{c}The calculation of *k _{c}* for the radial eigenmodes is displayed in Table III and matches closely with the estimate $ k c = \omega p 0 2 / c 2 + 8 \gamma / w c h 2 $ from Eq. (5).