# Broadband Terahertz Generation in a Corrugated Waveguide with matched Phase
and Group Velocities

Sergey S. Siaber1,2  Jonathan Gratus1,2  Rebecca Seviour3  Steven P.
Jamison1,2 and Taylor Boyd1,2  1Dept Physics, Lancaster University, Lancaster,
UK
2Cockcroft Institute of accelerator science, Daresbury, Warrington, WA4 4AD,
UK
3Ion Beam Centre, University of Huddersfield, Huddersfield, UK
<EMAIL_ADDRESS>
###### Abstract

We show that it is possible to design corrugated waveguides where phase and
group velocities coincide at an inflection point of the dispersion relation,
allowing an extended regime of interaction with a charge particle beam. This
provides a basis for designing travelling slow-wave structures with a
broadband interaction between relativistic charged particle beams and
propagating terahertz waves allowing an energy exchange between beam and wave,
amplifying terahertz radiation.

We employ Fourier-Mathieu expansion, which gives approximate analytic
solutions to Maxwell equations in a corrugated waveguide with periodically
undulating cross-section. Being analytic, this enables quick design of
corrugated waveguides, determined from desirable dispersion relations.

We design a three dimensional waveguide with the desired dispersion and
confirm the analytical predictions of the wave profile, using numerical
simulations. Madey’s theorem is used to analyse the strength of the wave-beam
interaction, showing that there is a broad frequency interaction region.

††journal: oe††articletype: Research Article

##  1 Introduction

Driven by the growning demand of applications, from material science to
telecommunications, from biology to biomedicine, recent years have seen a
rapid rise in the development of coherent terahertz (THz) sources.
Technologies used to generate THz radiation include laser-driven emitters,
solid state oscillators, gas and quantum cascade lasers. Laser driven
emitters, the most widely used sources of pulsed THz radiation, are based on
frequency down-conversion from the optical region. THz pulses energies
exceeding 1 mJ and 10’s MW peak power have been demonstrated in period poled
nonlinear sources driven by high energy near-IR ultrafast lasers[1], while
more main stream laser systems are capable of 10’s of $\mu$J energies and
efficiencies in region of 0.1% [2]. In solid state oscillators the transit
time of carriers through semiconductor junctions limits the frequency and
power that can be generated, generating around 100 mW at 100 GHz, but the
power falls off as $f^{-2}$[3]. Optically pumped gas lasers are the oldest
technologies for THz generation, generating between 0.3 to 5 THz at around
100mW [4]. Quantum cascade lasers are a relatively new approach to generating
THz radiation generating between 1 to 4 THz at mW power levels [5].

We consider an electron beam driven approach to THz generation. Electron beam
approaches are either non-relativistic/moderately relativistic vacuum
electronics devices (VED) or ultra-relativistic accelerator based radiators,
where a modulated electron beam generates THz via transition, Cherenkov,
Smith-Purcell, or undulator radiation. We present a moderately relativistic
VED approach that uses a novel corrugated metallic waveguide as a slow-wave
structure to generate THz radiation via coherent spontaneous emission. More
broadly metallic waveguides with corrugations or modulations on the metallic
boundary are of wide interest for their ability to tune the EM propagation
characteristics through the corrugation structure. Corrugations in the form of
rectangular groves orthogonal to the propagation axis and sub-wavelength
separation form slow-wave structures that give rise to dispersion relations
similar to that of a dielectric lined waveguides with the groove geometry and
spacing determine the effective dielectric properties. In such an arrangement
the dispersion relation can be obtained through a continuing sequence of mode-
matching along the waveguide with step-wise changes in cross-section. The
mode-matching over the repeating structure leads to a global eigenvalue
problem from which the dispersion relation is obtained.

Figure 1: Rectangular Corrugated Waveguide, with periodically undulating
cross-section. It is designed to vary in such a way that phase and group
velocities are very close in an extended frequency range.

In this paper, we use a waveguide which has smooth sinusoidal undulations in
the waveguide cross-section, its form is shown in figure 1. This waveguide,
has symmetric undulations in cross-section height, given by the function
${\cal L}_{x}(z)$, unlike sine-waveguide, considered in [6, 7], where cross-
section shape and dimensions are constant, but its centre is oscillating. For
this structure we can find explicit approximate solutions for the EM fields.
The geometry of the waveguide structure is parametrised by four parameters.
Two of these are traverse dimensions, one, denoted as ${\cal L}_{x}(z)$, is
varied along the waveguide and oscillates around the average value, marked as
$L_{0}$, and the other is fixed (denoted as $L_{y}$). The other two are the
spatial period of corrugation, $L_{z}$, and depth of the corrugations $q$.
Theoretical technique that predicts EM propagating mode in an analytic form
enables us to quickly scan this parameter space for desired dispersion and EM
field properties. We applied this approach to find the parameters where group
and phase velocities coincide, and do so at a point of inflection of the
dispersion relation. Such a point corresponds to particle-wave phase matching
at a broad range of frequencies. We call this point the coincident inflection
point (CIP). Research presented in this manuscript originated from the ideas
of EM waves propagation in corrugated wire media[8, 9, 10].

The concept of engineering the dispersion of VEDs using a corrugated waveguide
is not new. For example, some gyrotron amplifier designs utilise corrugated
waveguides[11] to engineer dispersion, matching phase velocity with that of an
electron beam, using a longitudinal sinusoidal modulation in a cylindrical
waveguide wall together with a 3-fold helical twist in the modulation. This
dispersion relation can be obtained from a perturbative coupling between
uniform cross-section waveguide guide modes, with an anti-crossing between
mode dispersion arising from the coupling [12, 13]. While gyrotron structures
achieve the desired dispersion, matching wave-particle velocity’s, they also
require a complex electron beam structure and helical propagation in high
magnetic fields; the electron beam must have annular and helical propagation
matching to couple to the EM mode [14].

More recently researchers have considered uniform rectangular cross-section
waveguides to engineer the dispersion in VEDs, such as the piecewise sine
waveguide [6] for high-power terahertz (THz) travelling wave tubes (TWTs).
This topology results in multiple modes existing simultaneously point in the
dispersion curve (see figure 4 in Zhang [6]), which can result in mode
competition and noise on the generated RF. Other researchers have considered a
uniform cross-section dielectric lined waveguides have found application in
supporting electron acceleration driven by terahertz-frequency laser pulses
[15, 16, 17, 18]. The dielectric layers allow for phase-matching of the drive
laser (at THz frequencies) $LSM_{10}$ modes with sub-relativistic electron
beams. However, unlike in the helical gyrotron structures, the phase-velocity
matching intrinsically comes with a group-velocity mismatch, and the temporal
walk-off between the EM pulse and electron beam is a limiting factor in
application of these dielectric lined waveguides to particle acceleration.

In this article, we present the design of a rectangular corrugated, with
periodically undulating cross-section, waveguide that enables efficient THz
generation through a number of dispersion relation characteristics. Achieving
these characteristics in this waveguide is enabled by our analytical
technique. We seek and find, a structure that has a dispersion that provides
both group and phase-velocity matching like a gyrotron, yet has the potential
to support axial beams like uniform cross section waveguides, without mode
competition. Our approach provides a perturbative solution to the waveguide
modes, from which we show that the dispersion relation is a solution to
Mathieu’s Floquet equation with two free parameters: a cut-off frequency
${\hat{\omega}}_{{\textup{c}}}$, determined by the relative dimensions
$L_{0}/L_{z}$ and $L_{y}/L_{z}$ of the waveguide, and $q$ which determines the
relative height of the undulations. Solution of the associated eigenvalue
problem provides the dispersion relation, and the Floquet coefficients
describing the longitudinal spatial-harmonics of the propagating mode.
Solutions with group and phase velocity simultaneously being equal, and
subluminal, are shown to exist. This solution can be used for the coupling
system as well as the beam/wave waveguide interaction region.

Unlike previous approaches, in this paper we determine the waveguide profile
from the desired dispersion relation and not set a priori, i.e. we determine a
family of waveguide geometries that satisfy the dispersion relation. Solutions
are of the form of a slowly varying sinusoidal-like longitudinal modulation in
the waveguide cross-section, and an EM mode with a longitudinal electric field
on axis. In addition to matching the phase velocity and electron beam
velocity, to maximise energy transfer between the electron beam and EM wave we
require the spatial components of the electric field of the EM wave to be
parallel with the beam propagation.

The paper is structured as follows. In section 2 we give the explicit
approximate EM modes for the corrugated waveguide, based on Mathieu’s
functions. We then demonstrate, in section 3, a method of finding the
appropriate parameters ${\hat{\omega}}_{{\textup{c}}}$ and $q$, which we apply
this to our goal of finding the CIP. In section 4, we use numerical CST
Microwave studio simulations confirm the accuracy of the analytic
approximation and to find how the CIP moves slightly when going from this
approximation to numerical simulations. In section 5 we show the force that a
charged particle would experience, and apply Madey’s theorem [19, 20, 21, 22]
to estimate THz energy being emitted as a result of particle–wave interaction
in the waveguide. We conclude by discussing future directions. In addition, in
the Supplemental Document we show how to construct the other EM modes, analyse
further Mathieu’s equation and how it leads to interaction zones, perform
error analysis on the approximate solutions, show that there are no subluminal
modes in the first interaction zone, and expand on the subject of simulating
numerical field patterns and comparing it to analytical field patterns.

##  2 Explicit Solution for a corrugated waveguide profile

We consider a rectangular corrugated waveguide shown in figure 1. We assume it
has two flat walls set a $L_{y}$ distance apart. The other two walls have a
symmetric undulating profile a distance ${\cal L}_{x}(z)/2$ from the centre,
i.e. the two walls are a distance ${\cal L}_{x}(z)$ apart. The cross section
dimensions are given by the constant width $L_{y}$ and the variable height
${\cal L}_{x}(z)$ which depends on the positions $z$ along the waveguide. We
choose the $(x,y)$ coordinates to be $(0,0)$ in the centre of the waveguide so
that the faces of the waveguides are at $y=\pm L_{y}/2$ and $x=\pm{\cal
L}_{x}(z)/2$. The undulations in ${\cal L}_{x}(z)$ are periodic, with a period
$L_{z}$, oscillating around $L_{0}$ value. The exact form of undulating height
${\cal L}_{x}(z)$ of corrugations will be given below.

Solving Maxwell’s vacuum equations in the frequency domain (with factor
$e^{-i\omega t}$)

$\displaystyle\nabla\times{\tilde{\boldsymbol{B}}}+i\omega\,c^{-2}{{\tilde{\boldsymbol{E}}}}={\boldsymbol{0}},$
$\displaystyle\nabla\cdot{\tilde{\boldsymbol{B}}}=0,$ (1)
$\displaystyle\nabla\times{\tilde{\boldsymbol{E}}}-i\omega{{\tilde{\boldsymbol{B}}}}=\boldsymbol{\cal
E}_{\textup{Max}},$ $\displaystyle\nabla\cdot{\tilde{\boldsymbol{E}}}=0.$ (2)

together with the boundary conditions

$\displaystyle{\tilde{\boldsymbol{E}}}_{\parallel}|_{\textup{Bdd}}=\boldsymbol{\cal
E}_{\textup{bdd}},$
$\displaystyle{\tilde{\boldsymbol{B}}}_{\perp}|_{\textup{Bdd}}={\boldsymbol{0}},$
(3)

where the errors $\boldsymbol{\cal E}_{\textup{Max}}$ and $\boldsymbol{\cal
E}_{\textup{bdd}}$ for our waveguide, are small. Our solution is thus an
approximation so it does not solve Maxwell’s equations exactly.

We consider ${\textup{TM}}_{p_{x}p_{y}}$ modes, with positive odd integers
$p_{x}$ and $p_{y}$, as these are the modes with a longitudinal component of
${\tilde{\boldsymbol{E}}}$ in the centre of the waveguide that we require. We
start by specifying the ansatz for the magnetic field
${\tilde{\boldsymbol{B}}}$,

$\displaystyle{\tilde{\boldsymbol{B}}}$
$\displaystyle=B_{0}\,c^{-2}\,(-i\omega)\,\phi(\eta\,z)\big{(}\kappa_{y}\cos(\kappa_{x}\,x)\,\sin(\kappa_{y}\,y)\,{\boldsymbol{e}}_{x}-\kappa_{x}\,\sin(\kappa_{x}\,x)\,\cos(\kappa_{y}\,y)\,{\boldsymbol{e}}_{y}\big{)}$
(4)

where
${\left\\{{\boldsymbol{e}}_{x},{\boldsymbol{e}}_{y},{\boldsymbol{e}}_{z}\right\\}}$
are the unit vectors along the axes,

$\displaystyle\kappa_{x}(z)=\frac{\pi\,p_{x}}{{\cal
L}_{x}(z)},\quad\kappa_{y}$
$\displaystyle=\frac{\pi\,p_{y}}{L_{y}},\,\,\textrm{and}\quad\eta=\frac{\pi}{L_{z}}.$
(5)

Trivially, ${\tilde{\boldsymbol{B}}}$ satisfies divergence equation in (1).
Substituting ${\tilde{\boldsymbol{B}}}$ into curl equations and combining
these into one equation yields the 2nd order ODE for $\phi(\eta\,z)$

$\displaystyle\phi^{\prime\prime}(\eta\,z)+\eta^{-2}\,\big{(}c^{-2}\omega^{2}-\kappa_{x}^{2}-\kappa_{y}^{2}\big{)}\,\phi(\eta\,z)$
$\displaystyle=0.$ (6)

From Gauss’s equation (1) the electric field ${\tilde{\boldsymbol{E}}}$ is

$\displaystyle{\tilde{\boldsymbol{E}}}$
$\displaystyle=B_{0}\,\Big{(}\big{(}\kappa_{x}^{\prime}\,\sin(\kappa_{x}\,x)+\kappa_{x}^{\prime}\kappa_{x}\,x\,\cos(\kappa_{x}\,x)\big{)}\phi(\eta\,z)+\eta\kappa_{x}\,\sin(\kappa_{x}\,x)\,\phi^{\prime}(\eta\,z)\Big{)}\cos(\kappa_{y}\,y)\,{\boldsymbol{e}}_{x}$
(7)
$\displaystyle\qquad\qquad\qquad+B_{0}\,\kappa_{y}\Big{(}-\kappa_{x}^{\prime}\,x\,\sin(\kappa_{x}\,x)\,\phi(\eta\,z)+\eta\cos(\kappa_{x}\,x)\,\phi^{\prime}(\eta\,z)\Big{)}\sin(\kappa_{y}\,y)\,{\boldsymbol{e}}_{y}$
$\displaystyle\qquad\qquad\qquad-
B_{0}\,(\kappa_{x}^{2}+\kappa_{y}^{2})\,\cos(\kappa_{x}\,x)\,\cos(\kappa_{y}\,y)\phi(\eta\,z){\boldsymbol{e}}_{z}$

We see ${\tilde{\boldsymbol{B}}}$, ${\tilde{\boldsymbol{E}}}$, given by eqs.
(4) and (7), automatically satisfy divergence equation in (2) since
$\nabla\cdot{\tilde{\boldsymbol{E}}}=\nabla\cdot\big{(}c^{2}(i\omega^{-1})\nabla\times{\tilde{\boldsymbol{B}}}\big{)}=0$.
In Supplemental Document, in the section SD2, we give the formula for the
error $\boldsymbol{\cal E}_{\textup{Max}}$. We also calculate the error
$\boldsymbol{\cal E}_{\textup{bdd}}$ in the boundary, and we show that these
errors are small if waveguide profile have a shallow gradient. We note that in
the centre of the waveguide, where $x=0$, then $\boldsymbol{\cal
E}_{\textup{Max}}|_{x=0}=0$. The general solution to Maxwell’s equations for
our waveguide is given by the sum of TM and TE modes. These are briefly
discussed in the supplementary document, section SD4. We presume that most of
the electrons are travelling along the centre of the waveguide, where solution
is more accurate.

Equation (6) can be rewritten in the form of Mathieu equation (here we defined
scaled variable $\zeta:=\pi z/L_{z}=\eta\,z$ )

$\displaystyle\phi^{\prime\prime}(\zeta)+\big{(}a-2q\cos(2\zeta)\big{)}\phi(\zeta)=0.$
(8)

Indeed, by stipulating an identity between corresponding terms in (8) and (6),
and substituting explicit form of $\eta$, $\kappa_{y}$ and $\kappa_{x}$, we
get

$\displaystyle a-2q\\!\cdot\\!\cos(2\zeta)$
$\displaystyle\equiv\eta^{-2}\,\big{(}c^{-2}\omega^{2}-\kappa_{x}^{2}-\kappa_{y}^{2}\big{)}$
(9) $\displaystyle=L_{z}^{2}\,\big{(}c^{-2}\pi^{-2}\omega^{2}-{p_{x}^{2}}{\cal
L}_{x}(z)^{-2}-p_{y}^{2}L_{y}^{-2}\big{)}$
$\displaystyle=L_{z}^{2}\,\big{(}c^{-2}\pi^{-2}\omega^{2}-{p_{x}^{2}}L_{0}^{-2}-p_{y}^{2}L_{y}^{-2}\big{)}-L_{z}^{2}\,\big{(}{p_{x}^{2}}{\cal
L}_{x}(z)^{-2}-{p_{x}^{2}}L_{0}^{-2}\big{)},$

where we choose the separation into the dimensionless Mathieu parameters $a$
and $q$ as

$\displaystyle
a=L_{z}^{2}\,\big{(}c^{-2}\pi^{-2}\omega^{2}-{p_{x}^{2}}L_{0}^{-2}-p_{y}^{2}L_{y}^{-2}\big{)},$
(10) $\displaystyle 2q\cos(2\zeta)=L_{z}^{2}\,\big{(}{p_{x}^{2}}{\cal
L}_{x}(z)^{-2}-{p_{x}^{2}}L_{0}^{-2}\big{)},$ (11)

In the following we choose a constant $q$, which defines the wave-guide
profile

$\displaystyle{\cal
L}_{x}(\zeta)=\Big{(}L_{0}^{-2}+2L_{z}^{-2}\,p_{x}^{-2}\,q\,\cos(2\zeta)\Big{)}^{-1/2}$
(12)

Under the assumption that $q\,L_{0}^{2}\,L_{z}^{-2}$ is small, waveguide
profile ${\cal L}_{x}$ takes simply a form of sinusoidally undulating function

$\displaystyle{\cal L}_{x}(\zeta)\approx
L_{0}-L_{0}^{3}\,L_{z}^{-2}\,q\,\cos(2\zeta)$ (13)

It this form, it is evident that $L_{0}$ is the average height of the
waveguide and parameter $q$ determines the height of the corrugations. For the
example waveguides to be discussed here we will choose $q=0.1$ and
$L_{0}/L_{z}\approx 1$.

Solutions of the equation (8) are the Mathieu special functions. These
solutions, as stated in Floquet’s Theorem [23], can be presented in the form

$\displaystyle\phi_{a,q}(\zeta)=e^{i{\hat{k}}\zeta}\,{\Psi}_{a,q,{\hat{k}}}(\zeta),$
(14)

where ${\Psi}_{a,q,{\hat{k}}}(\zeta)$ is periodic,

$\displaystyle{\Psi}_{a,q,{\hat{k}}}(\zeta+\pi)={\Psi}_{a,q,{\hat{k}}}(\zeta),$
(15)

where we have introduced the subscripts on $\phi_{a,q}(\zeta)$ to make
explicit the dependencies. ${\hat{k}}$ is the Mathieu exponent, with a value
that is determined by $a$ and $q$. In this form we see that ${\hat{k}}$ may be
considered as a dimensionless wavevector describing the wave propagation,
related to the wavevector by $k=\eta{\hat{k}}=\pi{\hat{k}}/L_{z}$. The
relationship between $k$ and $\omega$, or equivalently between $k$ and $a$,
therefore defines the dispersion relationship for the traveling wave solutions
to equations (6) and (8), while ${\Psi}_{a,q,{\hat{k}}}(\zeta)$ provides a
structure function for the field variation within a single period or cell of
the undulating waveguide.

${\hat{\omega}}$1.51.2551.06${\hat{k}}$01234superluminal particle,
$1^{\text{st}}$ zonesuperluminal particle, $2^{\text{nd}}$
zonelightlinesubluminal particle, $3^{\text{rd}}$ zone Figure 2: Corrugated
waveguide dispersion relation (blue line). The light line is shown in pink and
the three black lines represent three particle beam velocities that would
interact with monochromatic EM wave in: the first zone (corresponding
hypothetical particle beam velocity would have been $1.8c$), second zone
(similarly, hypothetical $\beta_{0}$ would have been $1.29$), and third zone
where matching velocity is subluminal (particle beam $\beta_{0}=0.53$).

The shape of the corrugated waveguide is determined by four parameters: the
length of one period of the corrugation $L_{z}$, the width of the waveguide
$L_{y}$, the separation of the corrugated surfaces at the mid point of the
corrugation $L_{0}$ and $q$ a parameter related to the depth of the
corrugation. Once these four parameters are chosen this defines a dispersion
relation between the angular frequency $\omega$ and the wavevector $k$ for the
corresponding EM wave. In practice, rather than specifying a frequency
$\omega$, or equivalently $a$, and then evaluating the corresponding
${\hat{k}}$, we perform the calculation in reverse. Specifying ${\hat{k}}$ and
the structure geometry the corresponding value of $a$ can be determined from
an eigenvalue evaluation [24].

An example of the dispersion relation obtained for $q=0.1$ is shown in Figure
2. To highlight the generality we have displayed dispersion relation with the
following dimensionless parameters,

$\displaystyle{\hat{\omega}}$
$\displaystyle=\frac{L_{z}}{\pi\,c}\,\omega\qquad\textup{and}\qquad{\hat{\omega}}_{{\textup{c}}}=\bigg{(}\frac{L_{z}^{2}}{L_{0}^{2}}+\frac{L_{z}^{2}}{L_{y}^{2}}\bigg{)}^{1/2}.$
(16)

so that equation (10) the Mathieu parameter $a$ becomes

$\displaystyle
a_{{\hat{k}},q}={\hat{\omega}}^{2}-{\hat{\omega}}_{{\textup{c}}}^{2}$ (17)

Thus, the number of parameters for the dispersion relation in the waveguide
reduces to just two, ${\hat{\omega}}_{{\textup{c}}}$ and $q$. For small values
of $q$, such as we are considering here, ${\hat{\omega}}_{{\textup{c}}}$ may
be regarded as the normalised cut-off frequency. Given these two parameters we
can quickly generate the scaled dispersion relation for ${\hat{\omega}}$ and
${\hat{k}}$.

One of the main advantages of this approach is that we can quickly scan the
parameter space in order to find a desirable dispersion relationship. A wide
range of behavior can be found through the parameterisation of equations
(16-17) and a numerically straightforward eigenvalue evaluation of
$a_{{\hat{k}},q}$.

##  3 Particle-wave synchronism.

For extended wave-particle interaction it is necessary to find propagating EM
mode in which charged particle and electromagetic phase wave remain in phase,
with particle ($c\beta_{e}$) and phase-velocity
($c{{\hat{k}}}/{{\hat{\omega}}}$) being equal. We aim to find more stringent
solutions where, for ceratin
${\hat{k}},{\hat{\omega}}({\hat{k}}),{\hat{\omega}}_{{\textup{c}}},q$ values,
the the phase-velocity and group velocity both match the particle velocity

$\displaystyle\frac{{\hat{\omega}}}{{\hat{k}}}=\frac{d{\hat{\omega}}}{d{\hat{k}}}=\beta_{e}$
(18)

where ${\hat{\omega}}({\hat{k}})$ depends on the values of the parameters
${\hat{\omega}}_{{\textup{c}}}$ and $q$. In order for the interaction to occur
in a range of frequencies (as opposed to a single frequency), we seek to
simultaneously satisfy a third constraint: zero group velocity dispersion
(GVD), formulated as having an inflection point on the dispersion when phase
and group velocities are matched, i.e.

$\displaystyle\frac{{\hat{\omega}}}{{\hat{k}}}=\frac{d{\hat{\omega}}}{d{\hat{k}}}\quad\textrm{and}\quad\frac{d^{2}{\hat{\omega}}}{d{\hat{k}}^{2}}=0$
(19)

Such EM mode maintains in-phase interaction, is free of EM pulse walk-off, and
minimizes the pulse dispersion within the waveguide. We refer to frequencies
and wavenumbers that satisfy these conditions as the coincident inflection
point (CIP). We regard, the CIP as an important design and operating point
that enables enhanced broadband interaction with an electron beam, and
provided initial synchronicity, continuous energy transfer from beam to wave
or from wave to beam.

In the dispersion diagram of figure 2 we show the ${\hat{\omega}}$-${\hat{k}}$
relationship for particles that would be phase-synchronised at normalized
frequency of ${\hat{\omega}}=1.5$ (this corresponds to $f\simeq 474\,$GHz in
the example structure described in detail later), and for normalised cut-off
frequency ${\hat{\omega}}_{{\textup{c}}}=1.255$ and corrugation parameter
$q=0.1$.

The first zone corresponds to a particle traversing a single cell in one
oscillation period; the 2nd zone to two oscillations per cell traversal, and
the 3rd zone to three oscillations per traversal. Synchronism in the first
zone is generally as it requires a superluminal particle, as is also the case
for the 2nd zone in this example. For the 3rd zone however, synchronism can be
obtained for particle velocity of $\beta_{e}=0.53$ (corresponding
approximately 92 keV electrons). The intercept of the particle and wave
dispersion lines (matched phase-velocity), with parallel tangents (matched
group velocity), occurs at the point of inflection in the wave dispersion
(zero GVD). It is therefore represents a CIP for $\beta_{e}=0.53$, and
${\hat{\omega}}=1.5,{\hat{\omega}}_{{\textup{c}}}=1.255,q=0.1$.

More generally, for a given waveguide undulation scale parameter $q$, we find,
one unique CIP. That is, for a given $q$, we can find the corresponding unique
values of ${\hat{\omega}}$, ${\hat{k}}$ and ${\hat{\omega}}_{{\textup{c}}}$
for the CIP. For a $q$ range from 0.0 to 0.3, corresponding CIP synchronous
particle velocity extends from 0.56$c$ to 0.47$c$. Third zone of the
dispersion relation for the ends of this range are shown on Fig.3 (a)
toghether with synchronous particle lines. Correspondence between particle
velocities and parameters $q$ and ${\hat{\omega}}_{{\textup{c}}}$ for which
CIP occurs is shown on Fig.3 (b) and (c) accordingly. We also observe that for
$q\to 0$, there is an absolute maximum velocity for the CIP, $v=0.5754c$.

${\hat{\omega}}$${\hat{k}}$$1.65$$1.33$$1.08$$2$$3$$q\\!=0.02$${\hat{\omega}}_{{\textup{c}}}\\!=1.325$$q\\!=0.3$${\hat{\omega}}_{{\textup{c}}}\\!=1.1$$\beta_{e}\\!=0.56$$\beta_{e}\\!=0.47$$\beta_{e}$$q$$0.56$$0.47$$0.02$$0.3$${\hat{\omega}}_{{\textup{c}}}$$q$$1.36$$1.1$$0.02$$0.3$(a)(b)(c)
Figure 3: End of range dispersion relations (black lines) and CIP particle
lines (blue lines) are shown on subfigure (a). CIP velocity dependence on $q$
values is shown on subfigure (b). Relation between CIP parameters $q$ and
${\hat{\omega}}_{{\textup{c}}}$ is given on subfigure (c).

##  4 Numerical simulation

Using numerical simulations is an alternative way of calculating dispersion
relation of the waveguide, as well as checking the validity of the analytical
predictions. We employ commercial numerical simulation software package CST
that solves discretised Maxwell integral equations on a tetrahedral mesh. We
simulate corrugated waveguiding structure described above which is one or
several periods long, depending on the particular aim of the simulation. We
aim to numerically find the CIP theoretically predicted in the previous
section, and to numerically confirm behaviour of the EM modes suitable for the
particle interaction. We start by simulating ten period (10 $L_{z}$) long
waveguide using an eigenmode solver. We set boundary conditions as perfect
electric conductor (PEC) on the walls of the waveguide, and as quasiperiodic
conditions at the entrance and exit of the waveguide, i.e.

$\displaystyle{\tilde{\boldsymbol{E}}}_{\parallel}|_{\textup{Bdd}}=0,\qquad{\tilde{\boldsymbol{E}}}\\!\left(z=0\right)=e^{i\theta}{\tilde{\boldsymbol{E}}}\\!\left(z=10L_{z}\right).$
(20)

Firstly, in simulation results, we find field pattern that correpsonds to TM
mode, – E-field concentrated in the centre of the waveguide and collinear with
$z$ axis, and B-field confined in the cross-section of the waveguide. An
example of such field pattern, namely crossections of E-field in $x-z$ and
$y-z$ planes and B-field in $x-y$ plane, found in CST simulations is given on
Figure 4. Geometrical parameters of the waveguide used in the simulations
correspond to numerical CIP (parameter values are $L_{z}=0.475{\textup{mm}}$,
$L_{y}=1{\textup{mm}}$, $L_{0}=0.5{\textup{mm}}$ and $q=0.1$, and synchronous
particle velocity in this case is 0.46).

Figure 4: Example of the simulated electric and magnetic field for the
corrugated structure. We see that E-field is strongest at the centre of the
waveguide; it is also collinear with $z$-axis.

We proceed to compare analytical and numerical predictions for dispersion
relations and field profiles in the centre of the waveguide. Analytical
dispersion (shown as a black line on the figure 5(d)) is calculated as a
Mathiue exponent, while numerical dispersion (orange line) $\omega(k(\theta))$
is found by varying phase $\theta$ in quasiperiodic boundary conditions (20)
and using wavenumber-phase relation ${\hat{k}}=2+\theta/\left(10\pi\right)$.
The third zone of dispersion relation, as predicted by analytical calculation
and numerical simulations is given in the Figure 5(d). Line of synchronous
particle, travelling at $0.46$ speed of light, is shown in blue. We note good
agreement between numerical and analytical dispersion, with an exception of a
region close to the bandgap.

To compare field profiles in the centre of the waveguide obtained in numerical
calculations and those predicted by analytical theory we need set a definite
phase of the EM mode in the waveguide. To this end, we consider a 6-period
long corrugated structure with PEC conditions at the entrance and exit of the
waveguide, i.e. find standing modes of the resonator, rather than a waveguide.
This enables us to find shapes of the longitudinal electric field
$E_{z}\left(z\right)$ in the centre of the waveguide for several points in the
third zone of dispersion, with wavenumbers ${\hat{k}}=2.17$, ${\hat{k}}=2.5$
and ${\hat{k}}=2.83$. Results, calculated for analytical CIP parameters
($L_{z}=0.475{\textup{mm}}$, $L_{y}=1{\textup{mm}}$, $L_{0}=1{\textup{mm}}$,
and $q=0.1$), are shown in Figure 5 (a)-(c). We observe very good agreement
for ${\hat{k}}=2.17$ and ${\hat{k}}=2.5$, and acceptable agreement in the
field shape even for the point ${\hat{k}}=2.83$, where the discrepancy between
analytical and numerical dispersion curves becomes greater, as we are nearing
the bandgap.

${\tilde{\boldsymbol{E}}}$${\tilde{\boldsymbol{E}}}$${\tilde{\boldsymbol{E}}}$$\zeta$$\zeta$$\zeta$${\hat{k}}=2.83$${\hat{k}}=2.5$${\hat{k}}=2.17$${\hat{\omega}}$${\hat{k}}$2.02.172.52.833.01.071.151.31.34(a)(b)(c)(d)$6\pi$$6\pi$$6\pi$$\pi$$\pi$$\pi$$4\pi$$4\pi$$4\pi$
Figure 5: Longitudinal field profiles are shown along the centre of the
waveguide $E_{z}\left(x\\!=\\!0,y\\!=\\!0,z\right)$ on the subfigures (a)-(c).
On the left longitudinal component of the electric field along the centre of
the waveguide for the 6 periods long structure at 3 different wavenumber
values. Numerical (orange) versus predicted analytical (black dashed). Three
pairs of corresponding points on the numerical (orange) and analytical(black
solid) dispersion curves are shown on the right (d). We observe good agreement
between the field patterns for the most of the zone, apart from the region of
higher frequencies as we start approaching the bandgap.

Finally, to check how well the CIP conditions (19) are met, we calculate phase
and group velocities from analytical and numerical dispersions,– shown in
Figure 6(a). On Fig. 6(b) we show how numerical phase velocity (orange dashed
line), numerical group velocity (orange solid line), and particle velocity
(blue line) are matched in the CIP. Analytical results are calculated for the
parameters that correspond to numerical CIP($q=0.1$,
${\hat{\omega}}_{{\textup{c}}}=1.06$). It is evident that analytical phase
velocity (black dashed line) is close to the numerical results, whereas there
is a considerable discrepancy in the group velocity.

We observe that exact CIP is a strong condition on the shape of the dispersion
curve. Achieving exact CIP for numerical 3D dispersion we had to change value
of $L_{0}$ from $0.41{\textup{mm}}$ to $0.5{\textup{mm}}$ comapred to
analytical CIP (i.e. changing ${\hat{\omega}}_{{\textup{c}}}$ from $1.255$ to
$1.06$). This results in CIP synchronous particle velocities and frequencies
being different, $\beta_{e}=0.46$ and 394 GHz for numerical CIP, and
$\beta_{e}=0.53$ and 474 GHZ for analytical CIP. Both numerical and anlytical
estimations relying on different assumptions and are approximations of the
real system.

$f\\!$ (GHz)$k$(mm-1)$458$$394$$350$$332$$13.2$$18$$19.4$$f\\!$
(GHz)$V_{phase/group}$, $\beta_{e}$0.640.50.460$458$$394$$350$$332$(a)(b)
Figure 6: Dispersion relations, numerical (orange) and analytical (black), for
numerical CIP parameters are shown on subfigure (a). On subfigure (b)
numerical phase (dashed orange) and group (solid orange) velocities are shown,
and synchronous particle line (solid blue) are shown toghether with analytical
predictions ( paremeter values $q=0.1$, ${\hat{\omega}}_{{\textup{c}}}=1.06$)
of phase and group velocities (dashed and solid black lines correspondingly).
Numerical CIP itself is marked as a blue point.

##  5 Waveguide Mediated Particle Wave Interaction

In this section, we employ two approaches to examine the particle-wave
interaction as described in section 2 above. The first is based on a single
particle model that remains in phase with the EM wave. This demonstrates the
relationship between the dispersion relation and the velocity of the particle.
The second approach is based on a perturbation technique, commonly referred to
as Madey’s theory [21]. This approach combines a perturbation expansion of the
Lorentz equation and averaging over the initial phases to study energy
exchange between beam and wave.

$1^{\text{st}}$ zone, ${\hat{k}}=0.857$$2^{\text{nd}}$ zone,
${\hat{k}}=1.143$$3^{\text{rd}}$ zone, ${\hat{k}}=2.857$$4^{\text{th}}$ zone,
${\hat{k}}=3.143$-1.0-0.500.51.0force (eV m-1)$\tfrac{\pi}{2}$$\pi$$\zeta$
Figure 7: The force felt by an in-phase electron in the first four zones. The
strength of the electric field is 1Vm-1 at the peak part of the waveguide. The
zone is determined by the speed of the particle. The first and second zones
are inaccessible since they corresponds to superluminal particles. Here
$q=0.1$, ${\hat{\omega}}_{{\textup{c}}}=1.255$ and ${\hat{\omega}}=1.513$.
First zone (black), ${\hat{k}}=0.857$, $\beta_{e}=1.766$. Total Force
${f_{\textup{P}}}=4.81$. Observe that the $\text{force}>0$ for all $\zeta$.
Second zone (blue), ${\hat{k}}=1.143$, $\beta_{e}=1.324$,
${f_{\textup{P}}}=-0.516$. Third zone (red). ${\hat{k}}=2.857$,
$\beta_{e}=0.530$, ${f_{\textup{P}}}=0.244$. This is at the CIP. Fouth zone
(green), ${\hat{k}}=3.142$, $\beta_{e}=0.481$, ${f_{\textup{P}}}=-0.043$.
-1.2-1.001.0Energy change $10^{-5}$ eV050010001500distance ($\mu m$) Figure 8:
Change in kinetic energy of a single electron, from 92KeV(i.e. synchronous
beam, $\beta_{0}=0.53$, third zone in figure 7), in an RF field of 1Vm-1, over
four periods in the waveguide. Here the waveguide is at the analytical CIP for
$q=0.1$, and geometric parameters $L_{z}=0.475{\textup{mm}}$,
$L_{y}=1{\textup{mm}}$, $L_{0}=1{\textup{mm}}$. The blue line is the average
energy loss.

As suggested in the previous section $E_{z}$ field shape obtained from
analytical model closely resembles field shape produced in numerical
simulations, and in this section we will use analytical field expressions to
estimate wave particle interaction in the structure.

To examine the interaction between charged particles and EM waves in our
structure we start by considering a single particle (charge $Q$) moving along
a trajectory in the centre of the waveguide. In general the on-axis
($(x=0,\,y=0)$ electric field $E(t,z)$ within the waveguide at time $t$ and
position $z$ is given by inverse fourier transform

$E(t,z)=\frac{1}{2\pi}\int{\tilde{\boldsymbol{E}}}(\omega,z)e^{-i\omega
t}d\omega$

A particle that enters the waveguide at time $t_{o}$, the particle will be at
position $z_{p}$ at time $t=t_{o}+z/c\beta_{e}\equiv t_{p}$. Therefore the
field experienced by the particle as it travels through the waveguide is

$E(t_{p},z_{p})=\frac{1}{2\pi}\int{\tilde{\boldsymbol{E}}}(\omega,z_{p})e^{-i\omega(t_{0}+z/c\beta_{e})}d\omega$

Substituting the the on-axis field from equation 7,

${\tilde{\boldsymbol{E}}}(\omega,z_{p})=B_{0}(\kappa_{x}^{2}+\kappa_{y}^{2})\int
e^{ikz}{\Psi}_{a,q,{\hat{k}}}\big{(}\pi\,L_{z}^{-1}\,z\big{)},$

the field experienced becomes

$\displaystyle E(t_{p},z_{p})$
$\displaystyle=\frac{1}{2\pi}B_{0}(\kappa_{x}^{2}+\kappa_{y}^{2})\int
e^{ikz_{p}}{\Psi}_{a,q,{\hat{k}}}\big{(}\pi\,L_{z}^{-1}\,z_{p}\big{)}e^{-i\omega(t_{0}+z_{p}/c\beta_{e})}\rho(\omega)d\omega$
(21)

where we have introduced the spectral density of the waveguide field,
$\rho(\omega)$. For a monochromatic field with
$\rho(\omega)=\delta(\omega-\omega_{s})+\delta(\omega-\omega_{s})$ at the
phase-matched frequency $\omega_{s}$ we obtain the effective field to be

$\displaystyle E(t_{p},z_{p})$
$\displaystyle=\frac{1}{\pi}B_{0}(\kappa_{x}^{2}+\kappa_{y}^{2}){\textup{Re}}\left[{\Psi}_{a,q,{\hat{k}}}\Big{(}\frac{\pi\,z_{p}}{L_{z}}\Big{)}\right]\cos(\omega_{s}t_{0})$
(22) $\displaystyle=\pi
B_{0}L_{z}^{-2}\left({\hat{\omega}}_{{\textup{c}}}^{2}+2\,q\,\cos(2\pi\,L_{z}^{-1}\,z_{p})\right){\textup{Re}}\left[{\Psi}_{a,q,{\hat{k}}}\Big{(}\frac{\pi\,z_{p}}{L_{z}}\Big{)}\right]\cos(\omega_{s}t_{0})$

where we have used the velocity phase-matched condition
$k=\omega_{s}/c\beta_{e}$, and the expansion of $\kappa_{x}$ from equations
(5) and (12).

The energy gain of a particle traveling on-axis over one structure period is
therefore

$\displaystyle U$ $\displaystyle=Q\int_{0}^{L_{z}}E(t_{p},z_{p})\,dz_{p}$
$\displaystyle=\frac{Q\,B_{0}\pi\,\cos(\omega_{s}t_{0})}{L_{z}^{2}}\int_{0}^{L_{z}}\left({\hat{\omega}}_{{\textup{c}}}^{2}+2\,q\,\cos\Big{(}\frac{2\pi\,z_{p}}{L_{z}}\Big{)}\right){\textup{Re}}\left[{\Psi}_{a,q,{\hat{k}}}\Big{(}\frac{\pi\,z_{p}}{L_{z}}\Big{)}\right]\,dz_{p}$
(23)

From figure 7 we see the energy depends on which zone the interaction occurs
in, determined by the velocity of the particle. In all cases
${f_{\textup{P}}}$ is highest when we are in the first zone, since there is no
counter force. However, we show in Supplemental Document, in the section SD3
that this is impossible, as it requires superluminal particle velocities. The
total interaction of the backward wave (second zone) is comparable to the
forward wave (third zone). In this case the backward wave would also need
superluminal particle velocities. By comparison the forward wave is
subluminal, and has a very good interaction between wave and particles, equal
to $12\%$ of the first zone interaction. Over several periods the energy of
the particle shows a net decreases as seen in figure 8, and a particle
synchronised with the opposite phase would exhibit a net increase in energy.

We extend this examination by considering the energy transfer between a
charged particle beam and the EM wave using Madeys theorem[20, 21], which has
been used in systems, from free electron lasers[21, 25], to conventional
traveling wavetubes[21], to metamaterial vacuum electronic devices[22]. Madeys
theorem relates the phase averaged energy spread to the phase averaged energy
change experienced by a charged particle as it propagates through a system.
The theory uses the first two terms of a perturbation expansion of the Lorentz
force based on a single charged particle analysis where (ignoring space
charge) a uniformly distributed charged particle beam is injected into the
system. The theorem draws links between the spontaneous emission of photons by
a single electron passing through the structure to stimulated emission of
photons. The first perturbation energy term of the Lorentz equation
$\gamma_{1}$ is taken at entry and exit from the structure, the difference
$\Delta\gamma_{1}$ is averaged over the phase of the EM field to yield
equation (24). This change in energy $\Delta\gamma_{1}$ relates to the
classical spontaneous power spectrum from Maxwell’s equations from the
beam[21],

$\left<\Delta\gamma^{2}_{1}\right>=\left<\left(\int\\!\\!\\!dz\frac{-QE_{z}\beta_{0}}{m_{0}c}\right)^{2}\right>.$
(24)

The second order term relates the energy change in the beam due to a
stimulated emission response, a consequence of a generalized framework in
Hamiltonian mechanics [19], given by,

$\left<\Delta\gamma_{2}\right>=\frac{1}{2}\frac{d}{d\gamma}\left<\Delta\gamma^{2}_{1}\right>.$
(25)

In a simplistic model we can consider the electron beam as $N$ electrons
entering the system every second, this enables us to write the power change in
the EM wave as,

$\Delta
P=-\frac{1}{2}\frac{d}{d\gamma}\left<\Delta\gamma^{2}_{1}\right>m_{0}c^{2}N.$
(26)

Using (26) we calculate power change in the 10 period long waveguide (i.e.
structure is 10 $L_{z}$ long). We assume that electron beam is in the centre
of the waveguide and comprises of $N=10^{4}$ electrons; value EM field
amplitude is taken as 1 V/m. We perform calculation for two energies of the
beam: $\beta_{0}=0.53$ \- which corresponds to the CIP as on Fig.2, and
$\beta_{0}=0.55$ \- which correspond to a single intersection point in the 3rd
zone.

${\hat{k}}$$\Delta P$ (eV)${\hat{\omega}}$${\hat{\omega}}$$k$ (mm-1)$f$
(GHz)$f$ (GHz) Figure 9: On subfigures (a) and (b) beam lines (black) and
third zone waveguide dispersion Relations (blue) are shown for coincident
inflecton point and for a simple single intersection correspondingly. EM wave
Power Change, calculated using Madey Theorem, is shown for a CIP (c) and CIP
(d).

Third zone dispersion relation and beam line for $\beta_{0}=0.53$, shown in 9
(a), demonstrate considerable interval beam-wave interaction frequency
interval. This interaction corresponds to the beam generating EM wave, as can
be seen on 9 (c), where, in the ${\hat{\omega}}$ frequency interval from
approximately $1.32$ to $1.54$, positive sign of $\Delta P$ corresponds to the
power transferred from the electrons to EM wave. For the geometric dimensions
of the waveguide used in the CST simulation corresponding frequency generation
interval is 77GHz: from $419$ to $496$ GHz. Graph $\Delta
P\left({\hat{\omega}}\right)$ terminates as ${\hat{\omega}}$ frequency
approaches band gaps.

The beam line for $\beta_{0}=0.55$ beam, shown in 9 (b), crosses dispersion at
a single point, without considerable interaction frequency interval. The EM
power $\Delta P$ (shown on 9 (d)) is generated in a shorter frequency
interval.

##  6 Conclusion

We have shown it is possible to design corrugated waveguides where the phase
velocity and group velocity of a wave coincide at a point of inflection in the
phase velocity, which we refer to as the coincident inflection point (CIP).
Using Madey’s theorem we have shown the CIP creates an extended regime of
interaction between the EM wave and a charged particle beam. Even for the case
where the CIP is close to the band gap (as seen in figure 2) the benefit of
having a CIP is the frequency range of interaction extends over almost the
whole band, as seen in figure 9.

The CIP has been found using a novel approach, based on an explicit analytic
approximate solution (4) and (7) to Maxwell’s equations, which involves
Mathieu’s equation. This allows us to estimate the geometrical parameters for
waveguides with an engineered dispersion curve. We observe good agreement
between the analytical dispersion curve and the CST numerical simulations.

We show that waveguiding structure and synchronous particle beam at a CIP have
potential for efficient broadband THz generation. We note again that
waveguiding structures discussed in this manuscript have corrugations which
mirror each other, figure 1, and are distinct from sine waveguides in which
crossection remains constant in shape and size and its centre undulates along
the waveguide. In contrast, in our approach, the crossection of the waveguide
varies and its centre remains constant. We present a waveguide based traveling
slow-wave structure that allows for an interaction between a charged particle
beam with a velocity of 0.53c and a propagating EM wave with a longitudinal
electric field.

This work can be extended to more general waveguides since the EM field given
by equations (4) and (7) are still good approximate a solution to Maxwell’s
equations even if the corrugation height is not given by equation (13) or
equation (12). Replacing equation (13) with a more general periodic function
may offer the possibility of extending the CIP for higher velocities. One goal
being to extend the maximum velocity, figure 3 to $c$, in order to exchange
energy with ultra-relativistic particles. Furthermore, we can replace ${\cal
L}_{x}(z)$ with any other functions as long as ${\cal L}_{x}^{\prime}(z)$ and
${\cal L}_{x}^{\prime\prime}(z)$ are small, and hence we can apply these
solutions for the EM field in the system that couples the corrugated waveguide
to the external structure.

Corrugated waveguides could also be used for the acceleration of charged
particles. This work will be extended to investigate the use of our corrugated
structures for particle acceleration, which would only require a slight change
in the synchronisation between wave and beam to flip this interaction from EM
wave generation to particle acceleration [22].

##  7 Backmatter

Funding SSS, JG, SPJ and TB are grateful the support provided by STFC (the
Cockcroft Institute ST/P002056/1 and ST/V001612/1). JG is particularly
grateful to the Peter Ratoff, director of Cockcroft (2014-2023), for
supporting this research. RS gratefully acknowledges support from the AFRL
Directed Energy Chief Scientist Office and the EOARD, grant FA8655-20-1-7002.

Disclosures The authors declare no conflicts of interest.

Data availability Data underlying the results presented in this paper are not
publicly available at this time but may be obtained from the authors upon
reasonable request.

Supplemental document See Supplemental Document for supporting content.

Author contribution statement JG suggested the original idea and provided the
theoretical work. SJ and RS proposed applying the idea to THz generation. SJ
proposed the idea of searching for the CIP. RS proposed the idea of using
Madey’s theorem, with SSS implementing it. The numerical simulations was
principally undertaken by SSS based on initial work by TB. SSS lead the
writing of the article, with help from all the authors.

## References

  * [1] F. Lemery, T. Vinatier, F. Mayet, R. Aßmann, E. Baynard, J. Demailly, U. Dorda, B. Lucas, A.-K. Pandey, and M. Pittman, “Highly scalable multicycle thz production with a homemade periodically poled macrocrystal,” Communications Physics 3, 150 (2020).
  * [2] C. D. W. Mosley, D. S. Lake, D. M. Graham, S. P. Jamison, R. B. Appleby, G. Burt, and M. T. Hibberd, “Large-area periodically-poled lithium niobate wafer stacks optimized for high-energy narrowband terahertz generation,” Opt. Express 31, 4041–4054 (2023).
  * [3] H. Eisele, A. Rydberg, and G. I. Haddad, “Recent advances in the performance of inp gunn devices and gaas tunnett diodes for the 100-300-ghz frequency range and above,” IEEE Transactions on Microwave Theory and Techniques 48, 626–631 (2000).
  * [4] A. Crocker, H. A. Gebbie, M. F. Kimmitt, and L. E. S. Mathias, “Stimulated emission in the far infra-red,” Nature 201, 250–251 (1964).
  * [5] C. Fattinger and D. Grischkowsky, “Terahertz beams,” Applied Physics Letters 54, 490–492 (1989).
  * [6] L. Zhang, Y. Jiang, W. Lei, P. Hu, J. Guo, R. Song, X. Tang, G. Ma, H. Chen, and Y. Wei, “A piecewise sine waveguide for terahertz traveling wave tube,” Scientific Reports 12 (2022).
  * [7] X. Yi, H. Zeng, S. Gao, and C. Qiu, “Design of an ultra-compact low-crosstalk sinusoidal silicon waveguide array for optical phased array,” Optics Express 28, 37505–37513 (2020).
  * [8] J. Gratus and M. McCormack, “Spatially dispersive inhomogeneous electromagnetic media with periodic structure,” Journal of Optics 17 (2015).
  * [9] T. Boyd, J. Gratus, P. Kinsler, and R. Letizia, “Customizing longitudinal electric field profiles using spatial dispersion in dielectric wire arrays,” Optics Express 26, 2478–2494 (2018).
  * [10] T. Boyd, J. Gratus, P. Kinsler, R. Letizia, and R. Seviour, “Mode profile shaping in wire media: Towards an experimental verification,” Applied Sciences (Switzerland) 8 (2018).
  * [11] V. L. Bratman, A. W. Cross, G. G. Denisov, W. He, A. D. R. Phelps, K. Ronald, S. V. Samsonov, C. G. Whyte, and A. R. Young, “High-gain wide-band gyrotron traveling wave amplifier with a helically corrugated waveguide,” Physical Review Letters 84, 2746–2749 (2000).
  * [12] S. J. Cooke and G. G. Denisov, “Linear theory of a wide-band gyro-twt amplifier using spiral waveguide,” IEEE Transactions on Plasma Science 26, 519–530 (1998).
  * [13] G. G. Denisov, V. L. Bratman, A. W. Cross, W. He, A. D. R. Phelps, K. Ronald, S. V. Samsonov, and C. G. Whyte, “Gyrotron traveling wave amplifier with a helical interaction waveguide,” Physical Review Letters 81, 5680–5683 (1998).
  * [14] G. S. Nusinovich, _Introduction to the Physics of Gyrotrons_ (JHU Press, 2004).
  * [15] E. A. Nanni, W. R. Huang, K. . Hong, K. Ravi, A. Fallahi, G. Moriena, R. J. Dwayne Miller, and F. X. Kärtner, “Terahertz-driven linear electron acceleration,” Nature Communications 6 (2015).
  * [16] L. Zhao, Z. Wang, H. Tang, R. Wang, Y. Cheng, C. Lu, T. Jiang, P. Zhu, L. Hu, W. Song, H. Wang, J. Qiu, R. Kostin, C. Jing, S. Antipov, P. Wang, J. Qi, Y. Cheng, D. Xiang, and J. Zhang, “Terahertz oscilloscope for recording time information of ultrashort electron beams,” Physical Review Letters 122 (2019).
  * [17] L. J. Wong, A. Fallahi, and F. X. Kärtner, “Compact electron acceleration and bunch compression in thz waveguides,” Optics Express 21, 9792–9806 (2013).
  * [18] M. T. Hibberd, A. L. Healy, D. S. Lake, V. Georgiadis, E. J. H. Smith, O. J. Finlay, T. H. Pacey, J. K. Jones, Y. Saveliev, D. A. Walsh, E. W. Snedden, R. B. Appleby, G. Burt, D. M. Graham, and S. P. Jamison, “Acceleration of relativistic beams using laser-generated terahertz pulses,” Nature Photonics 14, 755–759 (2020). Cited By :47.
  * [19] P. Luchini and H. Motz, _Undulators and free electron lasers_ (Oxford University Press, 1990), pp. 106–112.
  * [20] J. M. J. Madey, “Relationship between mean radiated energy, mean squared radiated energy and spontaneous power spectrum in a power series expansion of the equations of motion in a free-electron laser,” Il Nuovo Cimento B Series 11 50, 64–88 (1979).
  * [21] S. Yan, “The gain calculation of media and electrostatic free-electron lasers by the madey theorem,” IEEE Journal of Quantum Electronics 23, 1642–1645 (1987).
  * [22] Y. S. Tan and R. Seviour, “Wave energy amplification in a metamaterial-based traveling-wave structure,” EPL 87 (2009).
  * [23] M. S. P. Eastham, _The spectral theory of periodic differential equations_ (Scottish Academic Press, 1973).
  * [24] F.M.Arscott, _Periodic differential equations_ (Pergamon, 1964), pp. 124–127.
  * [25] J. M. J. Madey, “Stimulated emission of bremsstrahlung in a periodic magnetic field,” Journal of Applied Physics 42, 1906–1913 (1971).