# Adiabatic theory of slow extraction and its comparison with non-linear
resonance crossing experiment at VEPP-4M

S.A. Nikitin, BINP SB RAS, Novosibirsk, Russia

# Adiabatic theory of crossing third order resonance and its comparison with
experiment at VEPP-4M

S.A. Nikitin, BINP SB RAS, Novosibirsk, Russia

###### Abstract

The universal analytical approach was developed to describe the process of
particle slow extraction from a synchrotron using adiabatic crossing the
betatron resonance of the third order. The obtained formulas make it possible
to calculate the time diagrams for the current of the extracted beam (”spill
profile”), taking into account the chromaticity, particle momentum spread and
synchrotron oscillations. On this basis, a method for monochromatization of
extracted beam is indicated. The theory is compared with some results obtained
in the experiment on crossing the third-order resonance at the VEPP-4M storage
ring.

## Introduction

Since the 1960s, one of the widely known and used methods for extracting a
beam from a synchrotron is the slow crossing of a third-order nonlinear
betatron resonance [1]. In general terms, the theory of the method is based on
consideration of the region of stability of particle motion on the phase
plane, which decreases as the frequency of betatron oscillations approaches
resonance. Different interpretations of this approach are possible (see, for
example, [2, 3, 4]). The approach described in this paper with an in-depth
analysis of the adiabatic nature of the resonance crossing was developed in
the mid-90s, during the author’s participation in the research program on
obtaining polarized deuterons at the Nuclotron facility (JINR, Dubna) [5].
This approach was previously reported in abbreviated form in [6]. Its
distinguishing feature is clarification of the magnitude of the adiabatic
interval of motion, which lasts until the rapid increase in the amplitude of
particle oscillations. This is achieved through a rigorous analysis of the
phase integral (action integral or adiabatic invariant), which makes it
possible to define more precisely in what time a particle with given initial
conditions flies out from a beam. We find the adiabatic interval distribution
functions for the cases, which differ with respect to accounting the
chromaticity, the particle momentum spread and synchrotron oscillations. Using
these distribution functions, formulas are obtained to calculate the time
diagrams for the current of the extracted beam (”spill profile”) and the
momentum spread of particles in it. The theory indicates the possibility of
obtaining a monochromatic extracted beam due to a small, definitely dosed
acceleration of the particles in the process of crossing the resonance. In
some numerical examples, we use the Nuclotron parameters from [5].

On the storage ring VEPP-4M, a study was made of particle losses when crossing
a third-order nonlinear resonance in vertical motion [7]. In this paper, we
compare the form of the measured time diagrams of particle losses with those
calculated according to the presented theory.

##  1 Parameters of motion and non-linearity

Let list the known relations describing orbital motion of particles in a
synchrotron with quadratic nonlinearity of the guide field
$\frac{\partial^{2}B_{y}}{\partial x^{2}}$. The influence of the octupole
component of the field, which can stabilize the beam when crossing resonance,
is assumed to be negligible. Under these conditions, the amplitude of radial
betatron oscillations can increase indefinitely near the nonlinear resonance

$\nu_{x}=\frac{k}{3}+\delta,$ $None$

where $\nu_{x}$ is the betatron frequency of $x$-oscillations; $\delta<<1$ is
the detuning; $k$ is the number of the resonance harmonic of perturbations
$h(\theta)=\frac{1}{<B_{y}>}\frac{\partial^{2}B_{y}}{\partial x^{2}}$ induced
by a system of sextupoles ($<B_{y}>$ is the average guide field). The function
of azimuth $x(\theta)$, which presents radial oscillations of a particle with
a zero momentum deviation ($\delta p=0$), is found by the method of averaging
and has the Floquet form:

$x=a_{x}f_{x}e^{\imath\nu_{x}\theta}+c.c..$

Here, $|f_{x}|=\sqrt{\beta_{x}};$ $c.c.$ is the complex conjugate summand; the
amplitude $a_{x}=|a_{x}|e^{\imath\eta_{x}}$ is changing ’slowly’ by the
following law

$\frac{da_{x}}{d\theta}=a^{\prime}_{x}={\cal P}\bar{a}_{x}^{2}e^{\imath
3\delta},$ ${\cal P}=|{\cal P}|e^{\imath\,arg{\cal P}}=\left<\frac{\imath
h}{4}\bar{f}_{x}^{3}e^{-\imath k\theta}\right>.$

The angle brackets mean averaging over the synchrotron storage ring and the
bar denotes complex conjugation. In the variables

$I_{x}=\frac{|a_{x}|^{2}}{2},$ $w=\frac{1}{3}\left(\mbox{arg}{\cal
P}-3\eta_{x}-3\int\limits_{0}^{\theta}\delta d\theta\right)$

such motion corresponds to the perturbation Hamiltonian

$H=-I_{x}\delta-\frac{(2I_{x})^{3/2}}{3}|{\cal P}|\sin{3w}.$ $None$

Under conditions of slow extraction, Hamiltonian (1.2) depends on the azimuth
explicitly since the detuning varies, generally, by the law

$\delta=\delta_{0}-\int\limits_{0}^{t}\delta^{\prime}dt.$

##  2  Perturbation Theory method

Applying the expression for the x-motion Hamiltonian near the resonance
$3\nu_{x}=k$ (1.1), we find the phase integral (the action integral)

$\Gamma=\oint Idw$

for the case of the slow changing detuning $\delta=3\nu_{x}-k$:

$\Bigg{|}\frac{dH}{d\theta}\Bigg{|}=|I\delta^{\prime}|<<|H\delta|,$

where $\delta^{\prime}=\frac{d\delta}{d\theta}$. At a distance from the
resonance, $I=-H/\delta$ and thus this condition takes the following form:

$|\delta^{\prime}|<<\delta^{2}.$ $None$

Here $\Gamma$ can be considered an adiabatic invariant ($\Gamma=const$) and
one can determine the dependence of the radial oscillation amplitude
$a_{x}^{2}/2=I$ on $\delta(\theta)$ during slow extraction of particles from a
synchrotron.

At a distance from the resonance, at the initial amplitude of oscillations
$I=I_{i}$ and detuning

$|\delta_{i}|>>\sqrt{2I_{i}}|{\cal P}|,$

a zero approximation is true (${\cal P}\rightarrow 0$)

$I^{(0)}=<I>=-\frac{H}{\delta}=I_{i},$
$w^{\prime}_{0}=\frac{dw_{0}}{d\theta}=-\delta,$
$w_{0}=-\theta\delta-\frac{1}{3}\mbox{ arg}{\cal P}.$

$<I>$, a function of $\delta$, denotes the value $<|a_{x}|^{2}/2>$ averaged
over beatings of the closed phase trajectory $H=const$ in the plane $I,w$. Let
call this function the ’action’. Evidently, in a zero approximation
$\Gamma=2\pi<I>.$ Beatings taken into account, the Hamiltonian $H^{(1)}$ in a
first approximation has the following form

$H^{(1)}=-I^{(1)}\delta(1+\frac{\chi}{\tilde{\delta}}\sin{3w_{0}}),$ $None$

where

$\chi=\frac{2}{3}\frac{\sqrt{2I_{i}}}{\delta_{i}}|{\cal P}|$ $None$

is the perturbation parameter (the relative width of the resonance);
$\tilde{\delta}=\delta/\delta_{i}$, $\delta_{i}$ is the initial value of
detuning. The corresponding amplitude and phase are found from

$I^{(1)}=<I>(1-\frac{\chi}{\tilde{\delta}}\sin 3w_{0}),$
$w_{1}=w_{0}-\frac{\chi}{2\tilde{\delta}}\cos 3w_{0}.$

The phase integral in a first approximation has the following form:

$\Gamma=\int\limits_{0}^{2\pi}I^{(1)}dw_{1}=\int\limits_{0}^{2\pi}I^{(1)}\frac{dw_{1}}{dw_{0}}dw_{0}\simeq
2\pi<I>.$

Moreover, one can directly verify that

$<I>=-\frac{H_{i}+\int\limits_{0}^{\theta}<I>\delta^{\prime}d\theta}{\delta_{i}+\int\limits_{0}^{\theta}\delta^{\prime}d\theta}=const$

(we have used the relation $\sin 3w_{1}\simeq\sin
3w_{0}-\frac{3\chi}{2\tilde{\delta}}\cos^{2}3w_{0})$. Thus, the first
approximation allows one to find the oscillation amplitude

$<I>\frac{\chi}{\tilde{\delta}}=-\frac{H}{\delta}\frac{\chi}{\tilde{\delta}}=\frac{2}{3}<I>\frac{\sqrt{2I_{i}}}{\delta}|{\cal
P}|.$

The ’action’ remains unchanged: $<I>^{(1)}=I_{i}$.

In the second approximation,

$H^{(2)}=-I^{(2)}\delta(1+\frac{2}{3}\frac{\sqrt{2I^{(1)}}}{\delta}|{\cal
P}|\sin 3w_{1}),$ $w_{1}=w_{0}-\frac{\chi}{2\tilde{\delta}}\cos 3w_{0},$
$I^{(1)}\simeq I_{i}(1-\frac{\chi}{\tilde{\delta}}\sin 3w_{0}).$ $None$

Whence

$I^{(2)}=<I>^{(2)}\left[1-\frac{\chi}{\tilde{\delta}}\sin
3w_{0}+\frac{\chi^{2}}{\tilde{\delta}^{2}}(\frac{1}{2}+\cos^{2}3w_{0})\right],$
$w_{2}^{\prime}=-\delta-\frac{3}{2}\chi\frac{\delta}{\tilde{\delta}}(1-\frac{\chi}{2\tilde{\delta}}\sin
3w_{0})(\sin 3w_{0}-\frac{3}{2}\frac{\chi}{\tilde{\delta}}\cos^{2}3w_{0}),$
$\Gamma=\int\limits_{0}^{2\pi}I^{(2)}dw_{2}=2\pi<I>^{(2)}(1-\frac{10}{8}\frac{\chi^{2}}{\tilde{\delta}^{2}})=2\pi
I_{i}.$

In the second approximation, the ’action’ varies as

$<I>^{(2)}\simeq\frac{I_{i}}{1-\frac{10}{8}\frac{\chi^{2}}{\tilde{\delta}^{2}}}.$
$None$

This expression approximately describes how the beat-averaged squared
amplitude of x-oscillations grows with decrease of $\delta$ outside the
resonance band $\delta\sim\sqrt{2I_{i}}|{\cal P}|$. Below we consider the same
effect, using a more accurate method.

##  3 Phase Integral method

In order to analyze in more details the conditions of amplitude growth at the
band boundary, we will find an exact expression of the phase integral. Taking
the type of symmetry of the trajectories (1.2) at $H=const$ into account, we
can write

$\Gamma=6\int\limits\limits_{c}^{b}\frac{I}{dI/dw}dI,$ $None$

where $b$ and $c$ are, correspondingly, the maximum and minimum of $I$ at
periodic motion. Within the interval $0\leq w\leq 2\pi$ the closed trajectory
undergoes 3 full oscillations in accordance with the condition of extremum

$\frac{dI}{dw}=\frac{(2I)^{3/2}|{\cal P}|\cos 3w}{-\delta-\sqrt{2I}|{\cal
P}|\sin 3w}=0,$

or $\cos 3w=0$. Expressing the dependence of the integrand on $w$ via $I$ and
$H$, we obtain

$\Gamma=\frac{3}{2\sqrt{2}|{\cal
P}|}\int\limits_{c}^{b}\frac{(3H+I\delta)}{\sqrt{I^{3}-\frac{9}{8|{\cal
P}|^{2}}(H+I\delta)^{2}}}dI.$ $None$

The radicand in the integrand has a cubic polynomial form and can be presented
as

$I^{3}-\frac{9}{8|{\cal P}|^{2}}(H+I\delta)^{2}=(a-I)(b-I)(I-c),$ $None$

where $a>b>c$ are the real roots of the polynomial. The roots $b$ and $c$ are
chosen as the integration limits because of the type of symmetry of closed
phase trajectories mentioned above ($c<I<b$), and also because the root $a$
corresponds to the minimum of non-closed trajectories. The latter will be
shown below.

In the mentioned expression, integral (3.2) is a table one and at $a>b\geq
I>c$ equals

$\Gamma=\frac{3}{\sqrt{2}|{\cal
P}|}\left\\{\frac{(3H+a\delta)}{\sqrt{a-c}}K(k)-\delta\sqrt{a-c}E(k)\right\\},$
$None$

where $K$ and $E$ are complete elliptic integrals of the first and second kind
respectively; $k=\sqrt{\frac{b-c}{a-c}}$. Let write (3.4) using the
dimensionless variables $\tilde{H}=H/(I_{i}\delta_{i})$, $\tilde{I}=I/I_{i}$,
$\tilde{\delta}=\delta/\delta_{i}$,
$\tilde{a},\tilde{b},\tilde{c}=a/I_{i},b/I_{i},c/I_{i}$, as well as the
perturbation parameter

$\chi^{2}=\frac{8}{9}\frac{I_{i}|{\cal P}|^{2}}{\delta_{i}^{2}}.$

We obtain

$\Gamma=\frac{2I_{i}}{\chi}\left\\{\frac{3\tilde{H}+\tilde{a}\tilde{\delta})}{\sqrt{\tilde{a}-\tilde{c}}}K(k)-\tilde{\delta}\sqrt{\tilde{a}-\tilde{c}}E(k)\right\\},$
$None$

besides $\Gamma=2\pi I_{i}$. We will also use (3.3) in a form of equation of
extremums:

$\tilde{I}^{3}-\frac{\tilde{\delta}^{2}}{\chi^{2}}\tilde{I}^{2}-2\frac{\tilde{H}\tilde{\delta}}{\chi^{2}}\tilde{I}-\frac{\tilde{H}^{2}}{\chi^{2}}=0.$
$None$

The latter expression can be replaced by equivalent combined equations that
are derived from the Vieta theorem:

$\tilde{a}+\tilde{b}+\tilde{c}=\frac{\tilde{\delta}^{2}}{\chi^{2}},\ \
\tilde{a}\tilde{b}\tilde{c}=\frac{\tilde{H}^{2}}{\chi^{2}},\ \
\tilde{a}\tilde{b}+\tilde{b}\tilde{c}+\tilde{a}\tilde{c}=-\frac{2\tilde{H}\tilde{\delta}}{\chi^{2}}.$
$None$

We will add forth equation (3.5) to these combined equations

$-\tilde{H}=\frac{1}{3}\left[\frac{\pi\chi\sqrt{\tilde{a}-\tilde{c}}}{K(k)}+\tilde{a}\tilde{\delta}-\tilde{\delta}(\tilde{a}-\tilde{c})\frac{E(k)}{K(k)}\right].$
$None$

The algebraic system (3.7-3.8) is solved numerically for the decreasing series
of detuning values $\tilde{\delta}$. For example, by the iteration method with
the following initial approximations for the unknowns:

$\tilde{\delta}=1;\ \ \ \tilde{H}=-1;\ \ \ \tilde{a}=1/\chi^{2}-2;\ \ \
\tilde{b}=1+\chi;\ \ \ \tilde{c}=1-\chi.$

Normalized evolution of resonance detuning looks like
$\tilde{\delta}=1-2\tilde{\theta},$ where
$\tilde{\theta}=\theta/(\omega_{0}T_{ex})$, $T_{ex}$ is the time interval when
detuning changes by the value $2\delta_{i}$ with the rate
$\delta^{\prime}=const$.

Fig.1 shows an example of behavior of the solutions
$\tilde{b}(\tilde{\theta})$ and $\tilde{c}(\tilde{\theta})$ at $\chi=0.04$.
For the sake of comparison, the values $I_{max}/I_{i}=1+\chi/\tilde{\delta}$
and $I_{min}/I_{i}=1-\chi/\tilde{\delta}$ derived by the perturbation theory
(see Eq. (2.2) at $\sin{3w}_{0}=\pm 1$) method have also been plotted. The
third root $\tilde{a}$ (far from the resonance
$\tilde{a}\sim\chi^{-2}>>1$)decreases with detuning, approaching to the
growing root $\tilde{b}$ (Fig.2). When they finally approach, the discriminant
of equation (3.6) changes sign from positive to negative, which means that
there is only one real root (and two more complex conjugate roots). On the
phase plane, the trajectory, which before this critical moment looked like a
figure close to a triangle, opens at the trianagle vertices.

Figure 1: Variation of the maximal, minimal and average values of the action
for closed phase trajectories under adiabatic decreasing the resonance
detuning at $\chi=0.04$. The dotted line designates computation in terms of
the perturbation theory: $\tilde{I}_{max}=1+\chi/\tilde{\delta}$;
$\tilde{I}_{min}=1-\chi/\tilde{\delta}$;
$<\tilde{I}>=1+10\chi^{2}/(8{\tilde{\delta}}^{2})$. Solid lines indicate the
use of the exact expression for the phase integral. Figure 2:  Convergence of
the roots $\tilde{a}$ and $\tilde{b}$ with decreasing detuning ($\chi=0.04$).
Figure 3: Variation of the Hamiltonian in units of its initial value with
adiabatic decreasing rresonance detuning ($\chi=0.04$).

It follows from the determination of $dI/dw$ (see above) that, with that, it
is near the extremum $I_{e}=I_{max}=b$ ($\sin{3w}=-1$) where a singularity
arises because of the denominator going to zero:

$-\delta-\sqrt{2I}|{\cal P}|\sin{3w}\rightarrow 0.$

Actually, that means that the amplitude $I$ of this interval of phase
trajectory increases rapidly while the phase $w$ changes relatively slowly.
Otherwise, the condition for such regime of motion to arise can be expressed
as

$1-\frac{3\chi}{2\tilde{\delta}}\sqrt{1+\frac{\chi}{\tilde{\delta}}}=0$ $None$

or

$\chi/\tilde{\delta}=-1+\sqrt{\frac{7}{3}}\approx 0.5.$

For the plot in Fig.1, in particular, that takes place near the resonance
band, whose halfwidth is $\Delta\tilde{\delta}\approx 2\chi=0.08$. The value
$\chi=0.04$ corresponds to the initial amplitudes $I_{i}=\Sigma_{x}/8\pi$
where $\Sigma_{x}=2\pi\epsilon_{x}$ is the circulating beam emittance
multiplied by $2\pi$ (in the nuclotron, $\Sigma_{x}=13\pi\cdot 10^{-4}$
cm$\cdot$rad at 0.2 GeV/nucleon).

Figure 3 shows how the Hamiltonian function calculated in this way in units of
its unperturbed value ($I_{i}\delta_{i}$) changes in the process of
adiabatically slow decrease in the resonant detuning. It would be interesting
to note linearity of this dependence: $\tilde{H}\propto\tilde{\delta}$. It
follows from the given results that particles can be extracted slowly before
violation of adiabaticity condition (2.1), which as applied to the beam as a
whole takes the following form (here,$<..>$ mean averaging over an ensemble of
particles):

$\delta^{\prime}<<4\left<\chi^{2}\delta_{i}^{2}\right>\simeq\frac{1}{9}|{\cal
P}|^{2}\Sigma_{x}$ $None$

This conclusion can be easily verified by the given numerical example if the
adiabaticity condition is written in normalized variables as

$|\tilde{I}|<<|\tilde{H}\Delta\psi_{b}|,$

where $|\Delta\psi_{b}|=|\delta_{i}\tilde{\theta}|>>1$ is the betatron phase
incursion in the time $T_{ex}$.

Such is an evolution of the picture of motion on the ”action-phase” plane. Far
from resonance, it is represented by a circle, which, as the detuning
decreases, turns into a three-beam star which is close in shape to a triangle
with a closed trajectory. When crossing the border of the resonance band, the
trajectory opens at the vertices of the ”triangle”. In the well-known
”triangular separatrix” approach, the beam extraction begins when the phase
area occupied by it becomes less than the phase area occupied by the beam at a
certain critical value of the resonant detuning. For example, in [4] the
critical detuning is given by the equation

$\delta_{\Delta}\leq\frac{1}{4\sqrt{3\sqrt{3}}}|S_{k}|\frac{|a_{x}|}{\beta_{x}},$

where the ”sextupole harmonic” has a form

$S_{k}=\frac{1}{2B\rho}\int_{0}^{L}\frac{\partial^{2}B_{y}}{\partial
x^{2}}\beta_{x}^{3/2}\exp{[-ik\phi(s)]}ds.$

Here, $B\rho$ is the magnetic rigity; $L$ is a machine circumference;
$\phi=\int_{0}^{s}ds^{\prime}/(\nu_{x}\beta_{x})$. In the smooth approximation
$|S_{k}|/|{\cal P}|\approx 4\pi$ and $|a_{x}|^{2}\approx{\cal E}_{x}/4$. In
this case, it is easy to compare the value of the critical resonant detuning
from [4] with that given in (2.3):

$\frac{\delta_{\Delta}}{2\chi\delta_{i}}=\frac{3}{8\sqrt{3\pi\sqrt{3}}}\frac{|S_{k}|}{|{\cal
P}|}\approx 1.17.$

Thus, with respect to the specified parameter, the phase integral method is
quite consistent with the theory of the ”triangular separatrix”.

##  4 Adiabatic Interval method

Adiabatic motion starts at $\theta=0$, when the excitation is ’switched on’
and the resonance detuning begins lessening. This mode ends approximately at
$\theta=\Theta$, when the amplitude begins a rapid growth. From (3.9) we find
that the adiabatic motion interval $\Theta$ for an arbitrary particle in a
beam meets the equation

$2\chi=\tilde{\delta}=\frac{1}{\delta_{i}}\left[\delta_{i}-\int\limits_{0}^{\Theta}\left(\delta^{\prime}-\xi_{x}\frac{d}{d\theta}\frac{\Delta
p_{0}}{p_{0}}\right)d\theta+\xi_{x}\frac{\delta p}{p_{0}}\right],$ $None$

Here $\chi=2\rho|{\cal P}|/3\delta_{i}$, $\rho=|a_{x}|$;
$\delta_{i}=\delta_{0}+\xi_{x}(u_{0}+u)|_{\theta=0}$ , the initial detuning,
and $\delta_{0}$ is its average value in a circulating beam;
$\xi_{x}=\partial{\nu_{x}}/\partial{u}$ is the chromaticity coefficient;
$\delta^{\prime}=\Delta\nu_{x}/\omega_{0}T_{ex}$ is the rate of detuning
change in the range $\Delta\nu_{x}=2\delta_{0}$. The term
$u_{0}=\int\frac{d}{d\theta}\frac{\Delta p_{0}}{p_{0}}d\theta$ is a small
increment of the average beam momentum relative to the equilibrium one due to
a special beam acceleration, which, we assume, occurs at a constant speed
$(\Delta p_{0}/p_{0})/\omega_{0}T_{ex}$ during the time $T_{ex}$ of the
resonant detuning change. In the case of heavy particles, for this, in
principle, one can try to switch on an additional accelerating RF resonator
during the retuning of the betatron frequency, while the field in the bending
magnets does not change. In the electron synchrotron, for this purpose, it is
possible to apply a gradual change in the frequency of the master oscillator
of the RF system. The term $\xi_{x}u=\xi_{x}\frac{\delta p}{p_{0}}$ describes
the evolution of the detuning of an arbitrary particle, determined by the
deviation of its momentum without taking into account the contribution of the
total additional acceleration. It is owing to either synchrotron oscillations
or, as in the case of a debunched beam of heavy particles, a ”frozen” momentum
spread.

###  4.1 Zero momentum spread approximation

In (4.1), the parameters $\delta_{i}$ and $\chi$ depend on the magnitude of
the action $I_{i}$. Therefore, the particles in the beam will have different
values of the adiabaticity interval $\Theta$ in which they reach the edge of
the resonance band. First of all, we find the distribution function
$f_{\Theta}$ of the value $\Theta$ in the beam as well as its first two
moments - the average $<\Theta>$ and the variance $<\Theta^{2}>$ \- in the
monochromatic beam approximation ($\sigma_{p}=\delta p/p_{0}=0$). Using (4.1)
let express $\Theta$ in relative units for the case of $p_{0}$ linearity
($\Delta p_{0}^{\prime}=\Delta p_{0}/\omega_{0}T_{ex}$)

$\tilde{\Theta}=\frac{\Theta}{\omega_{0}T_{ex}}=\frac{1}{2q}\left(1-\frac{4}{3}\rho\frac{|{\cal
P}|}{\delta_{0}}\right),$ $None$ $q=1-\frac{\xi}{2\delta_{0}}\frac{\Delta
p_{0}}{p_{0}},$ $None$

where $\Delta p_{0}$ is the magnitude of $p_{0}$ increment in the time
$T_{ex}$. From (4.2) it follows, in particular, that at $\Delta
p_{0}\rightarrow 2p_{0}\delta_{0}/\xi$ and $\xi>0$ the real particle
extraction is stretched in time compared to the case of $\Delta p_{0}=0$. When
determining $f_{\Theta}$ ($\int\limits f_{\Theta}d\Theta=1$), it would be
convenient to express the value $\rho=|a_{x}|$ in the Twiss parameters [8]

$\rho=\frac{1}{2}\left[\frac{x^{2}}{\beta_{x}}+\beta_{x}(x^{\prime}+\frac{\alpha_{x}x}{\beta_{x}})^{2}\right]^{1/2},\
\ \alpha_{x}=-\frac{1}{2}\frac{d\beta_{x}}{ds},\ \ \
x^{\prime}=\frac{dx}{ds};$

$s=R\theta$, $R$ is the average radius of a synchrotron. Let introduce the
characteristic $w=4\rho^{2}\beta_{x},$ proportional to the radial oscillation
energy, whose average value is linked with the variances of a trajectory
deviation ($\sigma_{x}^{2}=<x^{2}>$) and a tilt
($\sigma_{x^{\prime}}^{2}=<x^{\prime 2}>$):

$<w>=4<\rho^{2}>\beta_{x}=\sigma_{x}^{2}(1-\alpha_{x}^{2})+\beta_{x}^{2}\sigma_{x^{\prime}}^{2}.$

In the smooth approximation
($\beta_{x}=R/\nu_{x},\alpha_{x}=0,\beta_{x}^{2}\sigma_{x^{\prime}}^{2}=\sigma_{x}^{2}$)
the probability density function of a random variable $w$ for the gaussian
beam is [8]:

$f_{w}=\frac{\exp{(-w/<w>)}}{<w>},$ $None$

$<w>=2\sigma_{x}^{2}$. The probability density function of $\rho$ can be
derived from the last equation and has the Rayleigh’s function form:

$f_{\rho}=\frac{\rho}{d}\exp{\left(-\frac{\rho^{2}}{2d}\right)},\ \ \
d=\frac{\sigma_{x}^{2}}{4\beta_{x}}.$ $None$

For completeness we present also the expressions for the mean and variance of
$\rho$

$<\rho>=\frac{1}{4}\sqrt{\frac{2\pi}{\beta_{x}}}\sigma_{x},\ \ \
D\rho=<(\rho-<\rho>)^{2}>=\left(2-\frac{\pi}{2}\right)d.$ $None$

Using the function $f_{\rho}$ and the definition (4.2) we obtain the
ptrobability density function of the interval $\tilde{\Theta}$ and then its
mean value and variance:

$f_{\Theta}=\frac{9\pi\delta_{0}^{2}q}{|{\cal
P}|^{2}\Sigma_{x}}|1-2\tilde{\Theta}q|\exp{\left\\{-\frac{9\pi\delta_{0}^{2}}{4|{\cal
P}|^{2}\Sigma_{x}}(1-2\tilde{\Theta}q)^{2}\right\\}},$ $None$
$<\tilde{\Theta}>=\frac{1}{2q}(1-\frac{1}{3}\sqrt{\frac{2\pi}{\beta_{x}}}\frac{|{\cal
P}|}{\delta_{0}}\sigma_{x})=\frac{1}{2q}(1-\frac{\sqrt{\Sigma_{x}}}{3}\frac{|{\cal
P}|}{\delta_{0}}),$
$<(\tilde{\Theta}-<\tilde{\Theta}>)^{2})>=<\delta\tilde{\Theta}^{2}>=\frac{(4-\pi)}{36\pi}\frac{|{\cal
P}|^{2}}{\delta_{0}^{2}q^{2}}\Sigma_{x},$ $None$

Here $\Sigma_{x}=2\pi\sigma_{x}^{2}/\beta_{x}\approx
2\pi\sigma_{x}^{2}\nu_{x}/R$ is characteristic area on the phase plane of
radial betatron oscillations (differs from the emittance by the factor $2\pi$)
.

Current of particles extracted is determined by the function
$f_{\Theta}=f_{\Theta}[\tilde{\Theta}(t)]$ and depends on the time ($t>0$) as

${\cal
J}(t)=eNf_{\Theta}\left|\frac{d\tilde{\Theta}}{dt}\right|=\frac{eN}{T_{ex}}f_{\Theta},$
$None$

$e$ is the elementary charge, $N$ is the number of deuterons accelerated. In
particular, in the ’smooth’ approximation ($\delta p=0$)

${\cal J}(t)=\frac{eN}{T_{ex}}\frac{9\pi\delta_{0}^{2}}{|{\cal
P}|^{2}\Sigma_{x}}\left|1-\frac{2tq}{T_{ex}}\right|\exp{\left\\{-\frac{9\pi\delta_{0}^{2}}{4|{\cal
P}|^{2}\Sigma_{x}}(1-\frac{2tq}{T_{ex}})^{2}\right\\}}.$ $None$

The average time of extraction and its variance equal to
$<t_{*}>=<\tilde{\Theta}>T_{ex}$ and $\sqrt{<\delta
t_{*}^{2}>}=\sqrt{<\delta\tilde{\Theta}^{2}>}T_{ex}$. The extracted beam gets
the spread in momentum

$<\delta p_{*}^{2}>^{1/2}=\Delta p_{0}<\delta\tilde{\Theta}^{2}>^{1/2},$
$None$

though the circulating beam was monochromatic in the approximation considered.
In reality such a case could approximately correspond to the condition
$\sigma_{p}<<\Delta p_{0}$.

Figure 4 presents an example of computation from (4.10) of the time dependence
of relative-unit current of extracted particles
($\tilde{J}=JT_{ex}/eN;\tilde{\Theta}=t/T_{ex}$). The spill time structure is
plotted for the following Nuclotron parameters: $E=0.2$ GeV/nucleon;
$\Sigma_{x}=13\pi\cdot 10^{-4}$ rad$\cdot$cm; $\delta_{0}=0.035;|{\cal
P}|=0.07$ cm-1/2. It is given that $\Delta p_{0}=0\hskip 5.69054pt(q=1)$.

Figure 4:  Extracted particle current $\tilde{J}={\cal J}T_{ex}/eN$ (in units
of average current $eN/T_{ex}$) as a function of time in units of $T_{ex}$ for
zero momentum spread in a circulating beam (or zero chromaticity).

###  4.2 Zero chromaticity

Practically, the above-considered approximation $\sigma_{p}=0$ can correspond
to the case of zero chromaticity of $x$ -oscillations at a special choice of
the system to correct quadratic non-linearity of the guide field
($\xi_{x}=0$). Then $\Theta$ does not depend on $\delta p$ and variance of
momenta in the extracted beam is given by the following expression:

$<\delta p_{*}^{2}>=\Delta p_{0}^{2}<\delta\tilde{\Theta}^{2}>+<\delta
p^{2}>,$

where $<\delta\tilde{\Theta}^{2}>$ is computed with the help of (4.6) and
(4.8) ($q=1$); $<\delta p^{2}>^{1/2}=\sigma_{p}$ is the momentum spread in the
circulating beam.

###  4.3 Consideration of momentum spread

Below, we consider taking into account the momentum spread and synchrotron
oscillations for the case of heavy particles, when the effect of radiative
diffusion can be neglected. If the beam is extracted with the RF field shut
down (no synchrotron oscillations), the condition for a particle with the
momentum deviation $p_{0}u$ to ’outfly’ can be written in the following form:

$2\chi\approx 1-\frac{\delta^{\prime}}{\delta_{i}}\Theta,$ $None$

where $\delta_{i}=\delta_{0}+\xi_{x}u$, or ($q\equiv 1$)

$\tilde{\Theta}=\tilde{\Theta}_{1}+\tilde{\Theta}_{2},$ $None$
$\tilde{\Theta}_{1}=\frac{1}{2q}\left(1-\frac{4\rho}{3}\frac{|{\cal
P}|}{\delta_{0}}\right),$
$\tilde{\Theta}_{2}=\frac{\xi_{x}}{2q\delta_{0}}\frac{\delta p}{p_{0}}.$

The function of distribution $f_{1}(\tilde{\Theta}_{1})$ coincides with (4.7).
The spread in momentum is gaussian, with the variance $\sigma_{p}^{2}$:

$f_{2}(\tilde{\Theta}_{2})=\frac{|\Lambda|}{\sqrt{2\pi}\sigma_{p}}\exp{\left(\frac{-\Lambda^{2}\tilde{\Theta}_{2}^{2}}{2\sigma_{p}^{2}}\right)}.$
$None$

Spread of the interval value $\tilde{\Theta}$ can be found via composition of
the functions $f_{1}$ and $f_{2}$:

$f(\tilde{\Theta})=\int\limits_{\tilde{\Theta}-1/2q}^{\tilde{\Theta}}f_{1}(\tilde{\Theta}_{1}-\tilde{\Theta}_{2})\cdot
f_{2}(\tilde{\Theta}_{2})d\tilde{\Theta}_{2}.$ $None$

Omitting bulky calculation, we will write the equation for current of
extracted particles:

${\cal
J}(\tilde{\Theta})=\frac{eN}{2T_{ex}}\sqrt{\frac{2}{\pi}}\frac{\mu|\Lambda|}{\sigma_{p}(\mu+\tau)}\left\\{\sqrt{\frac{\pi}{\mu+\tau}}\tau\ae\cdot
e^{-\frac{\mu\tau\ae^{2}}{\mu+\tau}}\left[\mbox{erfc}\left(-\frac{\tau\ae}{\sqrt{\mu+\tau}}\right)\right.\right.+$
$\mbox{}+\left.\left.\mbox{erfc}\left(\frac{\mu+\tau+\tau\ae}{\sqrt{\mu+\tau}}\right)-2\right]+e^{-\tau\ae^{2}}-e^{-\tau(\ae+1)^{2}-\mu}\right\\},$
$None$

Here and next we use the following designations:

$\Lambda=\frac{2q\delta_{0}}{\xi_{x}},\ \
\tau=\frac{\Lambda^{2}p_{0}^{2}}{8q^{2}\sigma_{p}^{2}},\ \
\mu=\frac{9\pi\delta_{0}^{2}}{4|{\cal P}^{2}\Sigma_{x}},\ \
\ae(\tilde{\Theta})=2q\tilde{\Theta}-1.$

In this case, spread of momenta in the extracted beam is defined completely by
that in the circulating one.

Figure 5 presents the time diagram of the current extracted from the Nucletron
for conditions in the example in section 4.1. A non-zero spread of momenta
$\delta p/p=10^{-3}$ is taken into account at the chromaticity $\xi_{x}=10$.
Other parameters are similar to those from the example in Figure 4.

Figure 5: Timing ($\tilde{\Theta}=t/T_{ex}$) of the extracted particle current
at the momentum spread of $10^{-3}$ in circulating beam (RF is shut down).

###  4.4 Non-coherent synchrotron oscillations

Generally ($\xi_{x}\neq 0$), the betatron frequency $\nu_{x}$ is modulated by
the law

$\nu_{x}=<\nu_{x}>+\xi_{x}\frac{\delta p}{p}\sin{(\nu_{s}\theta+\phi)},$

which leads to the splitting of the resonance (1.1) into a series of
resonances of the form

$3<\nu_{x}>+m\nu_{s}=k_{m}+\delta_{m}.$

Here $<\nu_{x}>$ is the average over the period of synchrotron oscillations
$2\pi/\nu_{s}$; $m$ is an order of the modulation resonance; $k_{m}$ is a
harmonic number; $\delta_{m}$ is a detuning. The power of side-band resonances
($m\neq 0$)

${\cal P}_{m}=(-i)^{m}J_{m}\left(\frac{3\xi_{x}}{\nu_{s}}\frac{\delta
p}{p}\right)\langle\frac{ih}{4}f_{x}^{*3}e^{-ik_{m}\theta}\rangle,$

$J_{m}(y)$ is the $m$-order Bessel function. At a small synchrotron frequency

$\nu_{s}<<6\langle\ \chi\delta_{i}\rangle\simeq|{\cal P}|\Sigma_{x}^{1/2},$
$None$

the nearest modulation resonances are fully inside the band of the principal
resonance (the modulation index $\sim 3\xi_{x}\sigma_{p}/\nu_{s}$ can be of
the order of 1 and larger). At an energy of 0.2 GeV/nucleon in the Nuclotron
this requires

$\nu_{s}<<4.5\cdot 10^{-3},\hskip
8.53581ptf_{s}=\nu_{s}\omega_{0}<<19\mbox{kHz}$

($\omega_{0}=4.2\cdot 10^{6}$ rad/s). For comparison, in the Nuclotron working
energy range, $f_{s}<10^{2}\mbox{Hz}$ ($\nu_{s}<2\cdot 10^{-5}$). So, for slow
extraction with RF voltage switched on, consideration of momentum spread is
limited to assessment of the corresponding initial spread of detuning from
resonance $\xi_{x}\delta p/p_{0}$. Since $\nu_{s}$ is small, the modified
condition of adiabaticity is met111In electron synchrotrons, the frequency
$\nu_{s}$ is much higher as compared to proton or ion machines, due to which
the conditions (4.17) and (4.18) may not be satisfied.:

$|\delta^{\prime}|+|\xi_{x}\nu_{s}\sigma_{p}|<<4\langle\chi\delta_{i}\rangle^{2}.$
$None$

At $\sigma_{p}\neq 0$ particles ’escape’ when ($\Delta p_{0}^{\prime}=$const,
$\delta^{\prime}=$const)

$2\chi\approx
1-\frac{\delta^{\prime}}{\delta_{i}}\Theta+\frac{\xi_{x}}{\delta_{i}}\frac{d}{d\theta}\frac{\Delta
p_{0}}{p_{0}}\Theta,$ $None$

where $\delta_{i}=\delta_{0}-\xi_{x}|\delta p/p_{0}|$, or

$\tilde{\Theta}=\frac{1}{2q}\left(1-\frac{4\rho}{3}\frac{|{\cal
P}|}{\delta_{0}}\right)-\frac{\xi_{x}}{2q\delta_{0}}\biggl{|}\frac{\delta
p}{p_{0}}\biggr{|}.$ $None$

In this case, the modulus of momentum deviation (i.g. an amplitudue value of
synchrotron oscillations) implies a positive initial detuning ($\delta_{0}>0$)
and, simultaneously, a positive chromaticity ($\xi_{x}>0$). Besides, the
frequency of synchrotron oscillations should be significantly larger than the
inverse time of passing the resonance band with a width of the order of
$\chi$. These features constitute an approximation in which we propose to take
into account synchrotron oscillations. With the chosen sign of the initial
detuning and chromaticity, it is the amplitude of these oscillations, together
with the changing values of the detuning and $\Delta p_{0}$, that determines
the moment when the critical resonance band is reached.

We write (4.4.3) as a sum
$\tilde{\Theta}=\tilde{\Theta}_{1}+\tilde{\Theta}_{2}$. The distribution
function $f_{1}(\tilde{\Theta}_{1})$ coincides with the similar one found for
the case of the spread of momenta of the debunched beam. The distribution
function $f_{2}(\tilde{\Theta}_{2})$ taking into account an amplitude
$|dp/p_{0}|$ differs from (4.3.3) only by a factor of 2:

$f_{2}(\tilde{\Theta}_{2})=\frac{2|\Lambda|}{\sqrt{2\pi}\sigma_{p}}\exp{\left(\frac{-\Lambda^{2}\tilde{\Theta}_{2}^{2}}{2\sigma_{p}^{2}}\right)}.$
$None$

Applying the composition of distribution functions and the definition (4.9),
we now find the expression for the spill structure in consideration of
synchrotron oscillations:

${\cal
J}(\tilde{\Theta})=\frac{eN}{T_{ex}}\sqrt{\frac{2}{\pi}}\frac{\mu|\Lambda|}{\sigma_{p}(\mu+\tau)}\left\\{\sqrt{\frac{\pi}{\mu+\tau}}\tau\ae\cdot
e^{-\frac{\mu\tau\ae^{2}}{\mu+\tau}}\left[\mbox{erfc}\left(-\frac{\tau\ae}{\sqrt{\mu+\tau}}\right)\right.\right.+$
$\mbox{}-\left.\left.\mbox{erfc}\left(\frac{\mu\ae}{\sqrt{\mu+\tau}}\right)\right]+e^{-\tau\ae^{2}}-e^{-\mu\ae^{2}}\right\\}.$
$None$

Unlike the case of the RF field shut down ($\nu_{s}=0$, Fig.5), the obtained
spill profile is asymmetric at similar other parameters (see the left curve in
Fig.6).

Figure 6:  The spill profile with taking into account the synchrotron
oscillations in two cases: without ($\Delta p_{0}=0$) and with a definitely
dozed acceleration ($\Delta p_{0}=p_{0}\delta_{0}/\xi_{x}$). The initial
momentum spread is $10^{-3}$; the chromaticity $\xi_{x}=10$.

Under these conditions, a dispersion of momenta in the extracted beam is found
from the expression

$<\delta p_{*}^{2}>={\cal D}p_{*}={\cal D}|\delta p|+(\Delta p_{0})^{2}{\cal
D}\Theta_{*}+2\Delta p_{0}\mathrm{Cov}\\{\Theta_{*},|\delta p|\\}.$

The symbols ${\cal D}...$ and $\mathrm{Cov}\\{...,...\\}$ mean a variance and
a covariance, respectively. The value $\Theta_{*}$ is defined by the condition
of ’outfly’ (4.19). As a result, we obtain the relation for the momentum
spread of extracted particles:

$<\delta p_{*}^{2}>\approx\frac{\pi-2}{\pi}\sigma_{p}^{2}\left(1-\frac{\Delta
p_{0}}{p_{0}\Lambda}\right)^{2},$ $None$

where $\Lambda=2q\delta_{0}/\xi_{x}$. As seen, the presence of slight
acceleration ($0<\Delta p_{0}/p_{0}<\Lambda$) at non-zero chromaticity
($\xi_{x}>0$) is advantageous for minimization of the momentum spread of the
extracted beam222In this case, apparently, the extraction should be carried
out in the area with zero dispersion. . The amplitude of the synchrotron
oscillations and the acquired additional momentum averaged over the beam are
related by the condition that the edge of the resonance band is reached. With
positive values of the initial detuning, chromaticity and additional
acceleration, the greater the amplitude, the lower the current increase in the
average momentum, and the earlier this event occurs. But the lower the
amplitude, the greater the current increase in the average momentum, the
critical band is reached later. In both cases, at the moment of crossing the
border of the band, the particle energy is the sum of two mentioned
components. In principle, a monochromatic extraction could be possible at
exact equality $\Delta p_{0}=p_{0}\Lambda$. The latter can be alternatively
written using the definition for $q$ (4.12) like

$\frac{\Delta p_{0}}{p_{0}}=\frac{\delta_{0}}{\xi_{x}}.$ $None$

In Fig.7, a particle with an amplitude of synchrotron oscillations $u_{1}$
reaches the boundary of the critical resonance band, ahead of a particle with
a smaller amplitude $u_{2}$. If (4.4.8 ) is satisfied, the second particle at
the band boundary will have a momentum increment equal to the amplitude
difference

$\Delta u=u_{1}-u_{2}=\frac{\Delta
p_{0}}{p_{0}}\frac{\Delta\delta}{\delta_{0}}.$

As a result, the momentum of each particle is $p_{0}+\Delta p_{0}$. It should
be kept in mind that we neglect the contribution of the additional
acceleration during the motion of particles inside the resonant band, since it
occurs in a relatively short time.

Figure 7: The positions of an equilibrium particle and two particles with
nonzero amplitudes of synchrotron oscillations are shown on the resonance
detuning scale.

In Fig.6, the spill diagram at right describes the case of the monochromatic
extraction. In contrast to the case of no dozed acceleration (the left curve),
the extraction with acceleration ends not at the moment when the resonant
detuning achieves 0 ($\tilde{\Theta}=0.5$), but at $\tilde{\Theta}=1$ (the
finish of the detuning change process).

Figure 8: The ratio of the momentum spread in extracted beam to that in
circulating one depending on the momentum increment due to the additional
acceleration at $\xi_{x}=10$ and $\delta_{0}=0.035$.

A numerical example of the dependence of the ratio of the momentum spreads of
the extracted and circulating beams on the magnitude of the additional
acceleration is shown in Fig.8. At $\Delta p_{0}\rightarrow
2p_{0}\delta_{0}/\xi_{x}$ (i.g. $q\rightarrow 0$) a time to reach the exact
resonance becomes infinitely large.

##  5 Experiment at VEPP-4M

At the storage ring collider VEPP-4M the experiment [7] was performed at 1.85
GeV to study the behavior of the electron beam near the nolinear resonance
$3\nu_{y}=23$, where $\nu_{y}$ is the vertical betatron tune. It was a
crossing of the resonance in various modes. In one of the modes the effective
nonlinearity related to the octupole lens was tuned to be close to zero. In
this case, no additional stability islands in the phase space emerge. Besides,
the vertical chromaticity was rather low: $\xi_{y}\approx 0.5$. In such
conditions one can estimate a law of the particle loss change in a time ($t$)
using the results of the adiabatic theory as applied for the case of vertical
oscillations and zero energy spread. In terms of slow extraction, the time-
dependent electron current is described by the equation (4.10) in which the
index $x$ is replaced by $y$:

${\cal J}(t)=\frac{eN}{T_{ex}}\frac{9\pi\delta_{0}^{2}}{|{\cal
P}|^{2}\Sigma_{y}}\left|1-\frac{2t}{T_{ex}}\right|\exp{\left\\{-\frac{9\pi\delta_{0}^{2}(1-\frac{2t}{T_{ex}})^{2}}{4|{\cal
P}|^{2}\Sigma_{y}}\right\\}}.$ $None$

Here,

${\cal P}=\left<\frac{\imath h}{4}\tilde{f}_{y}^{3}e^{-\imath
23\theta}\right>,$

the strength of the resonance; $h=(\partial^{2}B_{x}/\partial y^{2})$, the
skew sextupole perturbation in units of an average guide field;
$\tilde{f}_{y}$, the conjugated vertical Floquet function; $\Sigma_{y}$, the
vertical emittance multiplied by a factor of $2\pi$. Regading our experiment,
we replace the current of extracted beam with the rate of beam losses due to
crossing the resonance in arbitrary units. Beam losses were measured using a
synchrotron radiation sensor and a scintillation counter moved in the vertical
plane.

To adjust the required quadratic nonlinearity, it was planned to use the
existing skew sextupole lens with an estimated maximum resonance strength of
about $0.06$ mm-1/2 at the energy of experiment. However, as it was
established in the course of the experiment, in general, the strength of the
resonance was determined by the interference of the contributions from that
skew sextupole and from the chromaticity correction system in the arcs [7]. As
numerical simulation shows, the contribution of chromatic correction is
related to inclination of the axes of the eigenmodes of the transverse
oscillations with respect to the ideal directions along the vertical and
horizontal due to the finite betatron coupling333Recent experiments at VEPP-4M
on resonant excitation of transverse oscillations with their registration at
turn-by-turn pickup stations revealed the form of eigenmodes in the form of
ellipses with different inclinations of the axes, reaching a few tens of
degrees depending on the coupling strength.. The coupling effect was enhanced
using an skew-quadrupole corrector. A significant part of the experiments was
performed without the use of the skew-sextupole lens. These are the results
discussed below.

Figure 9:  Measured (in black) and calculated (in red) normalized losses of
particles in the process of crossing the resonance from bottom to top when
changing the frequency for 10 sec in the interval $\nu_{y}=0.6652\rightarrow
0.6687$ ((top figure), and from top to bottom for the same time in the
interval $\nu_{y}=0.6685\rightarrow 0.6653$ (bottom figure).

The resonant harmonic was estimated in the isolated resonance approximation
from the distortion of phase trajectories, as well as from the value of the
dynamic aperture near the resonance at a certain detuning [7, 9]. To do this,
the impact of the inflector excited vertical oscillations of the beam, and
according to the data of the pickup station, either a change in their
amplitude between the maximum and minimum values was observed, or the critical
value of the amplitude was fixed at which the beam was lost during a certain
number of turns. It was estimated that the amplitude of the resonant harmonic
could lie in the range ($0.02\div 0.07$) mm-1/2 or ($0.06\div 0.22$) cm-1/2 in
other units [7]. In the described experiment, an old electronic system with a
single-turn pickup resolution of 100 $\mu$m was used444In 2013–2016, new
pickup signal processing modules with a higher accuracy were developed and
installed at the VEPP-4M storage ring [10]., while the amplitude of the
oscillations excited by the impact was fractions of a millimeter. Therefore,
the accuracy of perturbation harmonic estimation was quite low. For this
reason, when comparing the experiment with the calculation by formula (4.25),
we choose the parameter ${\cal P}$ from the conditions of its closeness to the
specified range and the best fit to the measured data .

Fig. 9 shows the examples of comparing theory with experiment with a
significant coupling and relatively slow tuning of the frequency (by 0.003 per
10 s), in one case, towards its increase, in the other, towards decrease555A
similar figure is shown in [7], but the directions of detuning change are
confused in the caption.. The vertical emittance is estimated at the level of
5$\%$ of the calculated horizontal one, which is $2.5\cdot 10^{-6}$ cm $\cdot$
rad at 1850 MeV. When tuning from the bottom to up, the harmonic ${\cal
P}=0.3$ cm-1/2 was chosen, when tuning down ${\cal P}=0.5$cm-1/2 . The jagged
structure of the measured particle loss graphs is apparently associated with
the instability (albeit small - on the order of $10^{-4}$) of the power
sources of the quadrupole magnets. The given estimates of the harmonic [7]
refer to the definition of the Hamiltonian in the form $H=\delta\cdot
I_{y}+I_{y}^{3/2}A_{3}\cos{3w}$, where $A_{3}=2^{3/2}|{\cal P}|/3\approx|{\cal
P}|$. At $\delta_{0}=0.003$, $T=10$ seconds, $|{\cal P}|=0.3$ cm-1/2,
$\Sigma_{y}=2\pi\cdot 0.05\cdot 2.5\cdot 10^{-6}$ cm$\cdot$ rad the adiabatic
condition is fulfilled with a large margine:

$\frac{d\delta}{dt}<<\frac{\omega_{0}}{9}|{\cal
P}|^{2}\Sigma_{y}~{}~{}\rightarrow~{}~{}0.003~{}s^{-1}<<0.041~{}s^{-1}.$

Other loss diagrams obtained in a series of adiabatic resonance crossing, as a
rule, had a similar asymmetric shape - with a sharper leading edge at the
resonance point, as in Fig. 9, in accordance with the prediction of the theory
for the case of low chromaticity. At the same time, the loss of particles was
small - about 10% (the insertion depth of the scintillation counter was 7.3 mm
into a chamber with a diameter of about 8 cm). This contradicts the theory
which implies the total loss of particles. Radiative damping (disregarded by
the theory) and residual cubic nonlinearity for particles with large
oscillation amplitudes can be considered as an explanation of this fact

##  6 Discussion

For comparison, let us briefly note other cases in the experiment on crossing
a third-order resonance at VEPP-4M, which differ from those considered above.
If the condition of adiabaticity is violated, for example, by reducing the
rate of change in detuning by an order of magnitude (to thousandths of 1 s or
less) or by reducing the controlled coupling of oscillations, the shape of the
loss diagrams becomes close to symmetrical. In addition, the fraction of lost
particles dropped noticeably. At a high chromaticity value ($\xi_{y}=9.5$) ,
in addition to the main bell-shaped distribution (without sharper leading
edge), the similar, but lower amplitude distributions of losses on modulation
resonances with the frequency of synchrotron oscillations appeared. Since,
according to the definition (4.18 ), the adiabaticity condition is deeply
violated, this case is also not described by the theory. In the case of low
chromaticity, it can be concluded that the shape of the loss distribution
predicted in the theory is in good agreement with the experiment. At the same
time, the quantitative verification of the theory is not yet sufficiently
precise. In this regard, it would be interesting to repeat the experiment with
the use of a new existing system of turn-to-turn measurements, which will
provide a higher accuracy in determination of the harmonic amplitude compared
to that achieved earlier.A complete verification of the theory is apparently
possible at some ion synchrotron, in which the adiabaticity condition can also
be satisfied with allowance for synchrotron oscillations.

##  7 Acknowledgements

The author thanks A.M. Kondratenko for his interest and some remarks made at
the initial stage of the theoretical part of the work. The author is grateful
to all colleagues, participants of the experiment at VEPP-4M [7], especially,
E.B. Levichev and P.A. Piminov - for useful discussions. Special thanks to
N.A. Vinokurov and V.A. Vostrikov, who made valuable comments on the
manuscript.

## References

  * [1] A.W.Chao, M.Tigner, Handbook on Accelerator Physics and Engineering, 1998.
  * [2] V.P. Kartashev, V.I. Kotov. Metody formirovaniia puchkov chastits na uskoritliakh vysokikh energii. Moskva, Energoatomizdat, 1989.
  * [3] M. Pullia. Transverse aspects of the slowly extracted beam, CERN/PS 2000-06 (DR).
  * [4] S.V. Ivanov, O.P. Lebedev. Stokhasticheski medlennyi vyvod iz U-70. Novosti i problemy fundamental’noi fiziki, N 2 (6) 2009.
  * [5] B.V.Vasilishin, et al., Communication of the JINR, 9-86-512, Dubna 1986.
  * [6] S. Nikitin, Adiabatic theory of slow extraction of particles from synchrotron, in Proceedings of the European Particle Accelerator Conference (EPAC-2006), Edinburgh, United Kingdom, 2006 (Edinburgh, 2006), p.2119.
  * [7] O.V. Anchugov et al., Experiments on the physics of charged particle beams at the VEPP-4M electron-positron collider, Journal of Experimental and Theoretical Physics, 2009, vol. 109, No, 4 pp. 590-601.
  * [8] M. Sands. The physics of electron storage rings, SLAC-121, UC-28 (ACC), 1970.
  * [9] V. Smaluk, Nucl. Instr. Meth. A 578, 306314 (2007).
  * [10] E.A. Bekhtenev, G.V. Karpov, P.A. Piminov. Instruments and Experimental Techniques, 2017, No. 5, pp 74-80.


# Transmission in graphene through a double laser barrier

Rachid El Aitouni  Laboratory of Theoretical Physics, Faculty of Sciences,
Chouaïb Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Miloud
Mekkaoui  Laboratory of Theoretical Physics, Faculty of Sciences, Chouaïb
Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Ahmed Jellal
<EMAIL_ADDRESS>Laboratory of Theoretical Physics, Faculty of Sciences,
Chouaïb Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Canadian
Quantum Research Center, 204-3002 32 Ave Vernon, BC V1T 2L7, Canada

###### Abstract

In this work, we will study the transmission probability of Dirac fermions
through a double laser barrier. As part of the Floquet approximation, we will
determine the spinors in the five regions. Due to the continuity of the wave
function at the barrier edges, we find eight equations, each with infinity
modes. To simplify, we use the matrix formalism and limit our study to the
first three bands, the central band, and the first two side bands. From the
continuity equation and the spinors in the five regions, we will determine the
current density in each region, which makes it possible to determine the
expression of the transmission probability corresponding to each energy band.
The time-dependent laser fields generate several transmission modes, which
give two transmission processes: transmission with zero photon exchange
corresponds to the central band $\varepsilon$, and transmission with emission
or absorption of photons corresponds to the first two sidebands
$\varepsilon\pm\varpi$. One of the two modes can be suppressed by varying the
distance between the two barriers or the barrier width. The transmission is
not permitted if the incoming energy is below an energy threshold
$\varepsilon>k_{y}+2\varpi$. Increasing the intensity of the laser fields
makes it possible to modify the sharpness and amplitude of the transmission.

###### pacs:

78.67.Wj, 05.40.-a, 05.60.-k, 72.80.Vp
Keywords: Graphene, laser fields, Dirac equation, transmission channels, Klein
tunneling.

##  I Introduction

Graphene is a two-dimensional carbon-based material with a thickness of one
atom [1], its atoms are arranged in the form of a hexagonal like honeycomb
[2], it is isolated for the first time by the two researchers Giem and
Novoselov in 2004 [3], it has incredible electronic properties [4], exhibits
the quantum Hall effect [5], its has high mobility [7, 6] its electrons move
with a speed 300 times smaller than the speed of light, they are considered as
massless Dirac fermions [8], in addition to these electronic properties
graphene also has optical and photonic properties, the most attractive
property of graphene is the ability to absorb 2.3% of light on all ultraviolet
and infrared rays [9]. It is studied in the framework of tight-binding
Hamiltonian [10]. The dispersion relation in the vicinity of the Dirac points
is linear [11], that is to say the conduction and valance bands are in
contact, which implies that the fermions pass from one band to the other
easily with no effect, so they are uncontrollable. That’s the problem in
graphene, which delayed its use in the technological domain in addition to the
Klein tunnel effect [12, 13]. This effect is achieved experimentally by [14],
it shows that fermions of normal incidence cross the barrier even if their
energy is lower than the height of the barrier, whatever the width of the
barrier. At present, most research is focused on creating a band gap between
the valence and conduction bands, which makes it possible to control the
passage of fermions for technological use.

Several methods are created, for example, the deposition of the graphene sheet
on a substrate [15], the doping by another type of atom [16], the deformation
of the sheet [17, 18], the application of an external electric, magnetic, or
laser field. The effect of a potential barrier on fermions was already studied
in [19] but the Klein paradox is still present. The magnetic barrier [20, 21]
also shows the quantization of the energy spectrum with the appearance of
Landau levels, but the paradox is still present. The potential barrier
oscillates in time [22, 23] also allows to create energy sub-bands, each band
corresponds to a mode of transmission. The inclination of the barrier [24]
also allows to create a forbidding band but is not sufficient to confine all
the fermions. The irradiation of the barrier by laser field [26, 27, 25] makes
it possible to generate several energy bands, which give two types of
transmission, the transmission with zero photon exchange between the barrier
and the fermions, and the transmission process with photon exchange. The
increase in laser field intensity is capable of suppressing the transmission
inside the barrier for all transmission modes [28], so-called anti-Klein. In
the presence of a magnetic field, the intensity of the laser field makes it
possible to suppress the transmission process with zero photon exchange, but
it activates the process with photons exchange [29].

We investigated the behavior of Dirac fermions in graphene using double laser
barriers with varying amplitudes shifted by $\beta$. We begin by determining
the wave functions in each region by resolving the eigenvalue equation in each
region using the Floquet approximation, and then we employ the conditions at
the barrier edges, yielding eight equations, each for an infinite mode. We use
the matrix formalism to construct an infinite-order transfer matrix, and we
limit our analysis to the first three bands, where the central band
corresponds to $l=0$ and the first two sidebands correspond to $l=\pm 1$. The
vibration of the barrier over time generates many modes of transmission, the
transmission with zero photon exchange corresponds to the central band
$\varepsilon$ and the transmission of the side bands corresponds to the energy
side bands $\varepsilon+l\varpi$. To have a transmission, the following
conditions must be verified: $\varepsilon>k_{y}+2\varpi$, so we have an energy
threshold that must be moved for the fermions to cross the barrier. The
barrier parameters allow you to change the transmission mode and its
amplitude. Changing the space between the two barriers and the barrier width
allows you to change the amplitude of the two transmission modes or suppress
one of them. The variation of the laser field amplitude has a direct effect on
the amplitude of each transmission mode, and the two modes of transmission
vary sinusoidally.

The paper is organized as follows: After the introduction, we present in Sec.
II our theoretical model, which describes the movement of an electron in the
graphene sheet. To determine the spinors in the five regions, we solve the
eigenvalue equations. In Sec. III, we use the boundary conditions and also the
current density to find the expression of the transmission corresponding to
each energy band. We numerically present our results and discuss the basic
features of the system in Sec. IV. We close by concluding our results in Sec.
V.

##  II THEORETICAL MODEl

We consider a graphene sheet divided into five regions indexed by $j=1,2\cdots
5$. In the three regions $1$, $3$ and $5$, there is pristine graphene, but in
the two regions 2 and 4, we apply two different laser fields of amplitudes
$A_{2}$ and $A_{4}$, but phase shifted by the angle $\beta$, as shown in Fig.
1.

Figure 1: (Color online) Schematic of a graphene sheet in the presence of a
double laser barriers

Due to time independence of the Dirac Hamiltonian, the wave equation that
describes the movement of an electron through this double laser barrier is
given by $i\hbar\partial_{t}\Psi_{j}=H_{j}\Psi_{j}$ such that

$H=v_{F}\vec{\sigma}\cdot\left(\vec{p}+\frac{e}{c}\vec{A}_{j}(t)\right)$ (1)

where $\vec{\sigma}$ are the Pauli matrices $(\sigma_{x},\sigma_{y})$, $v_{F}$
is the Fermi velocity, $\vec{p}$ is the momentum operator, $\vec{A}_{j}(t)$
are the applied laser fields

$\vec{A}_{j}(t)=(0,A_{j}\cos\Phi_{j},0)$ (2)

with $\Phi_{j}$ and $A_{j}$ being the phase and amplitude of the laser field
in each region, which are defined by

$\displaystyle\Phi_{j}=\left\\{\begin{array}[]{l}\Phi_{1}\\\ \Phi_{2}\\\
\Phi_{3}\\\ \Phi_{4}\\\ \Phi_{5}\end{array}=\left\\{\begin{array}[]{c}0\\\
\omega t\\\ 0\\\ \omega t+\beta\\\ 0\end{array}\right.,\quad
A_{j}=\left\\{\begin{array}[]{l}A_{1}\\\ A_{2}\\\ A_{3}\\\ A_{4}\\\
A_{5}\end{array}=\begin{array}[]{c}0\\\ A_{2}\\\ 0\\\ A_{4}\\\
0\end{array}\right.\right.$ (23)

and $\beta$ is a phase shift. To determine the wave function corresponding to
each region, we solve the wave equation. The Hamiltonian contains two parts,
one spatial and the other temporal such that we can write it in the following
form:

$\displaystyle H_{j}=H_{0}+\widetilde{H}_{j}$ (24)

where we set

$\displaystyle
H_{0}=v_{F}\left(\sigma_{x}p_{x}+\sigma_{y}p_{y}\right),\quad\widetilde{H}_{j}=v_{F}\sigma_{y}A_{j}(t).$
(25)

Since $H_{0}$ is time independent and $\widetilde{H}_{j}$ is coordinate
independent, then the total wave function
$\Psi_{j}(x,y,t)=\psi_{j}(x,y)\phi_{j}(t)$ is a tensor product of the two
eigenvectors $\psi_{j}(x,y)$ and $\phi_{j}(t)$ associated with $H_{0}$ and
$\widetilde{H}_{j}$, respectively. In the framework of the Floquet
approximation [30], the temporal part is written as (in system unit
$\hbar=e=c=1$) $\phi_{j}(t)=\chi_{j}(t)e^{-i\varepsilon t}$ with
$\varepsilon=\frac{E}{v_{F}}$ is the Floquet energy and $\chi_{j}(t)$ is a
periodic function over time. To determine the expression of $\phi_{0}$, we
solve the eigenvalue equation $H_{j}\Psi_{j}(x,y,t)=E\Psi_{j}(x,y,t)$ and set
$\psi_{j}(x,y)=\dbinom{\varphi_{j}^{+}(x,y)}{\varphi_{j}^{-}(x,y)}$. This
yields

$\displaystyle\left[\partial_{x}+k_{y}-A_{j}\cos\Phi_{j}\right]\varphi_{j}^{-}(x,y)\chi_{j}(t)=\varphi_{j}^{+}(x,y)\frac{\partial}{\partial
t}\chi_{j}(t)$ (26)
$\displaystyle\left[\partial_{x}-k_{y}+A_{j}\cos\Phi_{j}\right]\varphi_{j}^{+}(x,y)\chi_{j}(t)=\varphi_{j}^{-}(x,y)\frac{\partial}{\partial
t}\chi_{j}(t).$ (27)

Unfortunately, we cannot solve the above system because there are three
unknown functions ($\varphi_{j}^{+},\varphi_{j}^{-},\chi_{j}$). To overcome
this situation, we may proceed with an approximation by assuming that inside
the barrier the laser-free coupled differential equations are satisfied by
$\varphi_{j}^{+}$ and $\varphi_{j}^{-}$. As a result, (26) and (27) reduce to
the following

$\displaystyle-A_{j}\cos\Phi_{j}\
\varphi_{j}^{-}(x,y)\chi_{j}(t)=\varphi_{j}^{+}(x,y)\frac{\partial}{\partial
t}\chi_{j}(t)$ (28) $\displaystyle A_{j}\cos\Phi_{j}\
\varphi_{j}^{+}(x,y)\chi_{j}(t)=\varphi_{j}^{-}(x,y)\frac{\partial}{\partial
t}\chi_{j}(t)$ (29)

which gives rise to the following second order differential equation

$\left({\partial_{t}^{2}}+\omega\tan\Phi_{j}\
\partial_{t}+A_{j}^{2}\cos^{2}\Phi_{j}\right)\chi_{j}(t)=0$ (30)

having the solution

$\chi_{j}(t)=e^{-i\alpha\sin\Phi_{j}}=\sum_{m=-\infty}^{+\infty}J_{m}(\alpha_{j})e^{-m\Phi_{j}}$
(31)

where $J_{m}$ is the Bessel functions. Combing all to write the spinor of $H$
as

$\Psi_{j}(x,y,t)=\psi_{j}(x,y)\sum_{m=\infty}^{\infty}J_{m}(\alpha_{j})e^{-i(\varepsilon
t+m\Phi_{j})}.$ (32)

To get a complete derivation of (32), we have to determine $\psi_{j}(x,y)$.
Indeed, in the regions 1, 3 and 5 there is only pristine graphene, then the
corresponding spinors can be written as [29, 25]

$\displaystyle\Psi_{1}(x,y,t)$ $\displaystyle=$
$\displaystyle\sum_{l,m=\infty}^{\infty}\left[\begin{pmatrix}1\\\
\gamma_{l}\end{pmatrix}\delta_{m,0}e^{ik^{0}_{x}x}+r_{l}\begin{pmatrix}1\\\
-\gamma^{*}_{l}\end{pmatrix}e^{-ik^{l}_{x}x}\right]e^{ik_{y}y}\delta_{m,l}e^{-iv_{F}(\varepsilon+m\varpi)t}$
(33) $\displaystyle\Psi_{3}(x,y,t)$ $\displaystyle=$
$\displaystyle\sum_{l,m=\infty}^{\infty}\left[c_{1l}\begin{pmatrix}1\\\
\gamma_{l}\end{pmatrix}e^{ik^{l}_{x}x}+c_{2l}\begin{pmatrix}1\\\
-\gamma^{*}_{l}\end{pmatrix}e^{-ik^{l}_{x}x}\right]e^{ik_{y}y}\delta_{m,l}e^{-iv_{F}(\varepsilon+m\varpi)t}$
(34) $\displaystyle\Psi_{5}(x,y,t)$ $\displaystyle=$
$\displaystyle\sum_{l,m=\infty}^{\infty}\left[t_{l}\begin{pmatrix}1\\\
\gamma_{l}\end{pmatrix}e^{ik^{l}_{x}x}+\mathbb{0}_{l}\begin{pmatrix}1\\\
-\gamma^{*}_{l}\end{pmatrix}e^{-ik^{l}_{x}x}\right]e^{ik_{y}y}\delta_{m,l}e^{-iv_{F}(\varepsilon+m\varpi)t}$
(35)

and the corresponding energy is

$\displaystyle\varepsilon+l\varpi=s_{l}\sqrt{(k^{l}_{x})^{2}+k_{y}^{2}}$ (36)

where
$\gamma_{l}=s_{l}\frac{k^{l}_{x}+ik_{y}}{\sqrt{(k^{l}_{x})^{2}+k_{y}^{2}}}=s_{l}e^{i\theta_{l}}$,
$\theta_{l}=\arctan\frac{k_{y}}{k_{x}^{l}}$, $\delta_{m,l}=J_{m-l}(0)$,
$s_{l}=\text{sgn}(\varepsilon+l\varpi)$, $\mathbb{0}_{l}$ is the null vector,
and $c_{il}$ $(i=1,2)$ are two constants. The coefficients $r_{l}$ and $t_{l}$
are, respectively, the reflection and transmission amplitudes, which can be
obtained from the boundary conditions. Here we have set with
$A_{j}=\frac{F_{j}}{\omega}$, $\alpha_{j}=\frac{F_{j}}{\omega^{2}}$ and
$\varpi=\frac{\omega}{v_{F}}$.

For both regions 2 and 4 there are the applied laser fields, and from the
eigenvalue equation we get

$\displaystyle\left(-i\partial_{x}-i(k_{y}-m\varpi)\right)\varphi_{2j}(x,y)$
$\displaystyle=$ $\displaystyle(\varepsilon+m\varpi)\varphi_{1j}(x,y)$ (37)
$\displaystyle\left(-i\partial_{x}+i(k_{y}-m\varpi)\right)\varphi_{1j}(x,y)$
$\displaystyle=$ $\displaystyle(\varepsilon+m\varpi)\varphi_{2j}(x,y).$ (38)

These can be worked out to end up with the solutions

$\displaystyle\Psi_{2}(x,y,t)$ $\displaystyle=$
$\displaystyle\sum_{l,m=\infty}^{\infty}\left[a_{1l}\begin{pmatrix}1\\\
\Gamma_{l}\end{pmatrix}e^{iq^{l}_{x}x}+a_{2l}\begin{pmatrix}1\\\
-\Gamma^{*}_{l}\end{pmatrix}e^{-iq^{l}_{x}x}\right]e^{ik_{y}y}J_{m-l}(\alpha_{2})e^{-iv_{F}(\varepsilon+m\varpi)t}$
(39) $\displaystyle\Psi_{4}(x,y,t)$ $\displaystyle=$
$\displaystyle\sum_{l,m=\infty}^{\infty}\left[b_{1l}\begin{pmatrix}1\\\
\Gamma_{l}\end{pmatrix}e^{iq^{l}_{x}x}+b_{2l}\begin{pmatrix}1\\\
-\Gamma^{*}_{l}\end{pmatrix}e^{-iq^{l}_{x}x}\right]e^{ik_{y}y}J_{m-l}(\alpha_{4})e^{-iv_{F}(\varepsilon+m\varpi)t}e^{-i(m-l)\beta}$
(40)

associated with the energy

$\displaystyle\varepsilon+l\varpi=s_{l}\sqrt{(q^{l}_{x})^{2}+(k_{y}-l\varpi)^{2}}$
(41)

where we have defined
$\Gamma_{l}=s_{l}\frac{q^{l}_{x}+i(k_{y}-l\omega)}{\sqrt{(q^{l}_{x})^{2}+(k_{y}-l\varpi)^{2}}}=s_{l}e^{i\theta_{l}^{\prime}}$.
The coefficients $a_{il}$ and $b_{il}$ $(i=1,2)$ are four constants.

##  III Analyzing transmission channels

In Appendix  A, we have determined all transmission channels and the
associated total transmission. Then, according to (73) and (74), we have

$\displaystyle T_{m}=\frac{\cos\theta_{m}}{\cos\theta_{0}}|t_{m}|^{2}$ (42)
$\displaystyle T=\sum_{m=-N}^{N}T_{m}$ (43)

where $\cos\theta_{m}=\frac{k_{x}^{m}}{\sqrt{(k_{x}^{m})^{2}+k_{y}^{2}}}$ and
$k_{x}^{m}=\sqrt{(\varepsilon+m\varpi)^{2}-k_{y}^{2}}$. To fully understand
the effect of the two laser barriers on the behavior of Dirac fermions through
the graphene sheet, we numerically represent our results. The oscillation of
the barrier over time generates several energy bands, which implies infinity
transmission modes: transmission with zero photon exchange corresponding to
the energy band $\varepsilon$, and transmission with photon exchange
corresponding to the sub-energy bands $\varepsilon+l\varpi$. To make the
graphic representation simpler, we focus only on the first three bands: the
central band and the first two side bands.

(a)

(b)

(c)

(d)

Figure 2: (Color online) Transmission probabilities as a function of the
incident energy $\varepsilon$ for $\omega=1.5$, $k_{y}=1$,
$\beta=\frac{\pi}{8}$, $F_{2}=F_{4}=1.5$, and different distances (a):
$d_{1}=3$, $d_{2}=5$, (b): $d_{1}=3$, $d_{2}=10$, (c): $d_{1}=1$, $d_{2}=5$,
(d): $d_{1}=1$, $d_{2}=10$. With $T_{1}$ (red line), $T_{-1}$ (green line),
$T_{0}$ (blue line), and $T$ (magenta line).

Fig. 2 depicts the transmission probability as a function of the incoming
energy $\varepsilon$ for various values of $d_{1}$, the spacing between the
two barriers, and $d_{2}$, the width of the two barriers. To have a
transmission, it is necessary that the following condition:
$\varepsilon>k_{y}+2\varpi$ should be fulfilled. As a result, we can say that
the quantity $k_{y}+2\varpi$ plays the role of an effective mass [31]. We
observe that the transmission varies in an oscillatory way, and the total
transmission oscillates in the vicinity of unity. Fig. 2a plotted for
$d_{1}=3$ and $d_{2}=5$, we see that the transmission with photon emission
$T_{1}$ (red line) decreases exponentially, the transmission with absorption
$T_{-1}$ (green line) increases along the $\varepsilon$-axis, and the
transmission $T_{0}$ with zero photon exchange (blue line) decreases rapidly
towards zero. When the width of the barrier is increased, $T_{-1}$ becomes
null from $\varepsilon=7$ but $T_{1}$ increases, and $T_{0}$ shows different
behaviors varying between decreasing then increasing along the
$\varepsilon-$axis as depicted in Fig. 2b. For the particular values $d_{1}=1$
and $d_{2}=10$, $T_{1}$ is almost zero for all incident energies, $T_{-1}$
increases then decreases exponentially, and $T_{0}$ shows the opposite
behavior (decreases then increases exponentially). From the energy $E=10$, we
notice that the transmission is carried out only with zero photon exchange
($T_{0}$) and then all fermions cross the barrier, showing the Klein tunnel.
When the width of the barrier is increased to $d_{2}=10$ in Fig. 2d, we
observe that $T_{0}$ is dominant between $\varepsilon=6$ and $\varepsilon=9$,
and after that, $T_{1}$ becomes more dominant but $T_{-1}$ is almost null. As
a conclusion, varying the two distances makes it possible to control the mode
of transmission, meaning that we can vary the mode of transmission by changing
the two distances, but the Klein paradox is still present.

(a)

(b)

Figure 3: (Color online) Transmission probability with zero photon exchange
$T_{0}$ as a function of $\varepsilon$ for $\beta=0$, $k_{y}=0.1$,
$\varpi=1.5$, $F_{2}=F_{4}=1.5$, and different distances ($d_{1},d_{2}$) such
that (a): $d_{2}=10$, $d_{1}=1$ (blue line), $d_{1}=5$ (red line), $d_{1}=9$
(green line), and (b): $d_{1}=3$, $d_{2}=4$ (blue line), $d_{1}=6$ (red line),
$d_{1}=9$ (green line).

In Fig. 3 we have plotted the transmission of the central band corresponding
to $l=0$ as a function of the incident energy of the fermions for different
values of $d_{1}$ and $d_{2}$ for well-determined values of other variables.
Fig. 3a depicts for various spacing values of $d_{1}$ between the two
barriers. We can see that as the distance increases, the transmission process
with zero photon exchange $T_{0}$ becomes more dominant because the two
barriers become like two peaks of width $d=d_{2}-d_{1}$. In this case, most of
the incident fermions cross the barrier with zero photon exchange, as clearly
seen in the green curve, which corresponds to $d=9$. Fig. 3b is plotted for
different values of the barrier width $d_{2}$, and we notice that for $d_{2}$,
which is almost equal to $d_{1}$ (red line), $T_{0}$ is very weak and
decreases in an oscillatory way. For $d_{2}=2d_{1}$ (green line), $T_{0}$
increases for low energies, then decreases exponentially towards zero in the
vicinity of $\varepsilon=11$. For $d_{2}=3d_{1}$ (blue line), $T_{0}$ becomes
more important than the other transmission process because the sum of the
three transmission modes is close to unity, as we have seen in the previous
figures. We can conclude that increasing the barrier width suppresses the
transmission of the side bands and increases the transmission with zero photon
exchange.

(a)

(b)

(c)

(d)

Figure 4: (Color online) Transmission probabilities as a function of $d_{1}$
for $d_{2}=10$, $k_{y}=1$, $\varpi=2$, $\varepsilon=20$, $\beta=0$, and
different amplitudes $F$ ($F_{1}=F_{4}=F$) such that (a): $F=0.5$, (b):
$F=1.9$, (c): $F=2.9$, (d): $F=3.9$. With $T$ (magenta line), $T_{0}$ (blue
line), $T_{1}$ (red line), $T_{-1}$ (green line).

In Fig. 4, we present the transmission probabilities as a function of the
distance between the two barriers $d_{1}$ for different values of the
amplitude $F$ of the laser field. Fig. 4a is plotted for $F=0.5$, and we
observe that the total transmission $T$ (magenta line) almost equals unit
whatever the distance $d_{1}$ because the laser fields are very weak and they
have almost a negligible effect. The transmission with zero photon exchange
$T_{0}$ oscillates in the vicinity of unit, and the transmissions ($T_{1}$,
$T_{-1}$) with photon exchange oscillate in the vicinity of zero. The laser
fields are very weak but allow for quantifying the energy, even though the
majority of the fermions cross the barrier with zero photon exchange. Fig. 4b
is plotted for $F=1.9$, and hence that the laser effect is very clear because
we observe that the transmissions with photon exchange vary periodically,
showing that $T_{1}$ (red line) decreases and $T_{-1}$ increases along the
$d_{1}$-axis. However, $T_{0}$ varies in phase opposition with the two other
transmission modes, with an increase in the amplitude of the oscillations
along the $d_{1}$ axis. The total transmission always oscillates in the
vicinity of the unit, meaning the Klein effect is present. Now, by increasing
to the value $F=2.9$ in Fig. 4c, $T_{0}$ oscillates along the $d_{1}$-axis,
and with an increase in the amplitude of the oscillations, it vanishes for
very precise values. The intervals at which this transmission is obtained are
maximal, and the transmissions with photon exchanges are zero. Fig. 4d is
plotted for $F=3.9$, that is to say $\alpha=\frac{F}{\omega^{2}}$ is closed to
$1$, we observe a decrease in the interval where $T_{1}$ and $T_{-1}$ get
canceled with an increase in the number of peaks. We see that $T_{1}$
oscillates with a decrease in amplitude, and in contrast, $T_{-1}$ increases
along the $d_{1}$-axis. As a result, it appears that increasing the amplitude
$F$ of the laser fields makes it possible to decrease the interval of
cancellation of the transmissions, but it also increases the number of
oscillations.

(a)

(b)

(c)

(d)

Figure 5: (Color online) Transmission probabilities as a function of $d_{2}$
for $d_{2}=10$, $k_{y}=1$, $\varpi=2$, $\varepsilon=20$, $\beta=0$, and
different amplitudes $F$ ($F_{2}=F_{4}=F$) such that (a): $F=0.5$, (b):
$F=1.9$, (c): $F=2.9$, (d): $F=3.9$. With $T$ (magenta line), $T_{0}$ (blue
line), $T_{1}$ (red line), $T_{-1}$ (green line).

Fig. 5 is similar to the previous figure with a variation of the barrier
widths and the same parameter values. Fig. 5a is plotted for $F=0.5$, and then
the generation of transmission modes is observed, but the effect of the laser
fields is very weak. The transmissions $T_{1}$ and $T_{-1}$ oscillate around
zero, which can be neglected in comparison to $T_{0}$, implying a very clear
Klein effect. For the value $F=1.9$ in Fig. 5b, $T_{1}$ and $T_{-1}$ vary
periodically with the increase in amplitude, but $T_{0}$ varies in an
oscillatory way with the decrease in the amplitude along the $d_{2}$-axis. For
$F=2.9$, Fig. 5c shows that $T_{0}$ varies regularly with the appearance of
peaks in the minimum part. For $F=3.9$ in Fig. 5d, we observe that the total
transmission $T$ oscillates around the unit, while $T_{0}$ is more dominant,
which is oscillating between zero and one, canceling out at several points,
$T_{1}$ and $T_{-1}$ also oscillate with the increase in amplitude along the
$d_{2}$-axis. For example, in the vicinity of the value $d_{2}=4$, $T_{1}$ and
$T_{-1}$ almost null, but $T_{0}$ is equal to the unit, which implies that all
the fermions cross the barrier without exchanging photons. There are several
points where transmissions get canceled. From the width of the barriers, it is
possible to control the passage by which transmission mode.

Fig. 6 presents the transmission probabilities as a function of the phase
shift $\beta$ for $k_{y}=1$, $\varpi=2$, $\varepsilon=12$, $F_{2}=2$,
$d_{1}=1.3$, $d_{2}=3$, and different values of the amplitude $F_{4}$ of the
second laser barrier. For $F_{4}=0.9$ in Fig. 6a, we notice that the three
transmission modes vary in an oscillating way with the same frequency for
different amplitudes. The transmission with zero photon exchange $T_{0}$ is
more dominant, while the transmission with photon emission $T_{1}$ and
absorption $T_{-1}$ vary in the same way. A similar result is obtained in [32]
for a double oscillating potential. The total transmission $T$ oscillates
around the unit, which implies the presence of the Klein effect. Fig. 6b is
plotted for $F_{4}=2$ that is to say $\alpha_{1}=\alpha_{2}$, in this case,
the amplitude of the oscillations increases, but $T_{0}$ is still more
dominant and $T$ is still fluctuating around the unit. In Fig. 6c with
$F_{4}=2.9$, we see that the transmissions become non-sinusoidal but vary
periodically. $T_{0}$ is always more dominant, and its amplitude is almost
equal to $0.9$, whereas $T_{1}$ and $T_{-1}$ vary symmetrically. In addition,
we observe that the Klein tunnel presents periodically for very precise values
of $\beta$. When $\alpha_{4}$ tends towards $1$, Fig. 6d shows that the
amplitude of $T_{0}$ decreases with the increase in the number of
oscillations, and the Klein tunnel is exhibited periodically as before.

(a)

(b)

(c)

(d)

Figure 6: (Color online) Transmission probabilities as a function of the phase
shift $\beta$ for $k_{y}=1$, $\varpi=2$, $\varepsilon=12$, $F_{2}=2$,
$d_{1}=1.3$ $d_{2}=3$, and different value of amplitude $F_{4}$ such that (a):
$F_{4}=0.9$, (b): $F_{4}=2$, (c): $F_{4}=2.9$, (d): $F_{4}=3.9$. With $T$
(magenta line), $T_{0}$ (blue line), $T_{1}$ (red line), $T_{-1}$ (green
line).

##  IV Conclusion

We have studied the behavior of Dirac fermions through double laser barriers
generated by two electric fields, the first of amplitude $F_{2}$ and the
second of amplitude $F_{4}$, and of frequency $\omega$ shifted by a phase
shift $\beta$. The two barriers divide the graphene sheet into five regions.
In regions 1, 3, and 5, there is only pristine graphene, and the other two
regions are irradiated by laser fields. Within the framework of Floquet’s
approximation, we solved the eigenvalue equation to determine the wave
functions corresponding to each region. The oscillation of the barrier over
time generates several energy bands. Then, using the boundary conditions of
each barrier, we obtained eight equations, each of which has several modes. To
simplify the calculation, we used the matrix formalism to arrive at a transfer
matrix of infinite order. The latter is difficult to solve, and for this
reason we have limited our study to the first three energy bands: the central
band corresponds to the energy $\varepsilon$ and the first two side bands
correspond to the $\varepsilon+\varpi$. The current densities are also used to
determine the transmission coefficient corresponding to each energy band.

The numerical analysis of our theoretical results shows that the transmission
exists if the incidence energy of the Dirac fermions satisfy the condition
$\varepsilon>k_{y}+2\varpi$, such that this threshold plays the role of an
effective mass. The vibration of the barrier over time produces two
transmission processes: the transmission process with zero photon exchange and
the process with photon exchange. The variation of the distance between the
two barriers makes it possible to cancel one of the two processes, and the
variation of the width of the barriers makes it possible to change the
transmission process. When the distance between the two barriers $d_{1}$ is
increased, the number of oscillations increases, and the transmission with
zero photon exchange becomes more dominant. The decrease in the width of the
barrier $d_{2}$ makes it possible to reduce the transmission $T_{0}$ even if
the incident energy increases. The increase in laser field amplitude increases
the number of oscillations and their amplitude. The Klein tunnel is observed
periodically, and it is obtained for very precise values of the phase shift
$\beta$.

## Acknowledgment

We thank Prof. A. H. Alhaidari for valuable discussions.

## References

  * [1] A. K. Geim and K. S. Novoselov, Nat. Mater. 6, 183 (2007).
  * [2] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
  * [3] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * [4] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * [5] Y. Zhang, Y. W. Tan, H. L. Stormer, and P. Kim, Nature 438, 201 (2005).
  * [6] S. Morozov, K. Novoselov, M. Katsnelson, F. Schedin, D. Elias, J. Jaszczak, and A. Geim, Phys. Rev. Lett. 100, 016602 (2008).
  * [7] K. I. Bolotin, K. J. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim, and H. L. Stormer, Solid State Commun. 351. 146 (2008).
  * [8] A. K. Geim, Science 324, 1530 (2009).
  * [9] Qiaoliang Bao, Han Zhang, Bing Wang, Zhenhua Ni, Candy Haley Yi Xuan Lim, Yu Wang, Ding Yuan Tang, and Kian Ping Loh, Nature Photon 5, 411 (2011).
  * [10] S. Reich, J. Maultzsch, C. Thomsen, and P. Ordejon, Phys. Rev. B 66, 035412 (2002).
  * [11] N. M. R. Peres, J. Phys.: Condens. Matter 21, 323201 (2009).
  * [12] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Nat. Phys. 2, 620 (2006).
  * [13] C. W. J. Beenakker, Rev. Mod. Phys. 80, 1337 (2008).
  * [14] N. Stander, B. Huard, and D. Goldhaber-Gordon, Phys. Rev. Lett. 102, 026807 (2009).
  * [15] A. N. Sidorov, M. M. Yazdanpanah, R. Jalilian, P. J. Ouseph, R. W. Cohn, and G. U. Sumanasekera, Nanotechnology 18, 135301 (2007).
  * [16] G. A. K. P. A. Giovannetti, P. A. Khomyakov, G. Brocks, V. V. Karpan, J. van den Brink, and P. J. Kelly, Phys. Rev. Lett. 101, 026803 (2008).
  * [17] F. Guinea, M. I. Katsnelson, and A. K. Geim, Nat. Phys. 6, 30 (2010).
  * [18] G.-X. Ni, Y. Zheng, S. Bae, H. R. Kim, A. Pachoud, Y. S. Kim, C.-L. Tan, D. Im, J.-H. Ahn, B. H. Hong, and B. Ozyilmaz, ACS Nano 6, 1158 (2012).
  * [19] H. Bahlouli, E. B. Choubabi, A. El Mouhafid, and A. Jellal, Solid State Communi. 151, 1309 (2011).
  * [20] L. Dell’Anna and A. De Martino, Phys. Rev. B 79, 045420 (2009).
  * [21] A. Jellal and A. El Mouhafid, J. Phys. A: Math. Theo. 44, 015302 (2011).
  * [22] A. Jellal, M. Mekkaoui, E. B. Choubabi, and H. Bahlouli, Eur. Phys. J. B 87, 123 (2014).
  * [23] M. Ahsan Zeb, K. Sabeeh, and M. Tahir, Phys. Rev. B 78, 165420 (2008).
  * [24] A. El Mouhafid and A. Jellal, J. Low Temp. Phys. 173, 264 (2013).
  * [25] R. El Aitouni, M. Mekkaoui, A. Jellal, and M. Schreiber, arXiv:2307.03999 (2023).
  * [26] R. Biswas, S. Maitty, and C. Sinha, Physica E 84, 235 (2016).
  * [27] R. Biswas and C. Sinha, Appl. Phys. 114, 183706 (2013).
  * [28] R. El Aitouni, M. Mekkaoui, and A. Jellal, Ann. Phys. (Berlin) 535, 2200630 (2023).
  * [29] R. El Aitouni and A. Jellal, Phys. Lett. A. 447, 128288 (2022).
  * [30] C. W. J. Beenakker, Rev. Mod. Phys, 80, 1337 (2008).
  * [31] M. Mekkaoui, A. Jellal, and H. Bahlouli, Solid State Communi. 358, 114981 (2022).
  * [32] M. Mekkaoui, E. B. Choubabi, A. Jellal, and H. Bahlouli, Mater. Res. Express 4, 035002 (2017).

##  Appendix A Determining transmission channels

The continuity of the wave function at the barrier edges makes it possible to
find eight equations, each of which has an infinity mode. The analytical
resolution of this system is difficult. For simplification, we use the matrix
formalism, which allows us to write this system in matrix form with an
infinite order transfer matrix. To simplify, we limit our study to the three
first bands: the central band corresponds to $l=0$ and the two first side
bands correspond to $l=\pm 1$. For the first barrier, we have
$\Psi_{1}(-d_{2},y,t)=\Psi_{2}(-d_{2},y,t)$ and
$\Psi_{2}(-d_{1},y,t)=\Psi_{3}(-d_{1},y,t)$, which give

$\displaystyle\delta_{m,0}e^{-ik^{0}_{x}d_{2}}+r_{m}e^{+ik^{m}_{x}d_{2}}=\sum_{l=-\infty}^{+\infty}\left(a_{1l}e^{-iq^{l}_{x}d_{2}}+a_{2l}e^{+iq^{l}_{x}d_{2}}\right)J_{m-l}(\alpha_{2})$
(44)
$\displaystyle\delta_{m,0}\gamma_{m}e^{-ik^{0}_{x}d_{2}}-r_{m}\gamma^{*}_{m}e^{+ik^{m}_{x}d_{2}}=\sum_{l=-\infty}^{+\infty}\left(a_{1l}\Gamma_{l}e^{-iq^{l}_{x}d_{2}}-a_{2l}\Gamma^{*}_{l}e^{+iq^{l}_{x}d_{2}}\right)J_{m-l}(\alpha_{2})$
(45) $\displaystyle
c_{1m}e^{-ik^{m}_{x}d_{1}}+c_{2m}e^{+ik^{m}_{x}d_{1}}=\sum_{l=-\infty}^{+\infty}\left(a_{1l}e^{-iq^{l}_{x}d_{1}}+a_{2l}e^{+iq^{l}_{x}d_{1}}\right)J_{m-l}(\alpha_{2})$
(46) $\displaystyle
c_{1m}\gamma_{m}e^{-ik^{m}_{x}d_{1}}-c_{2m}\gamma^{*}_{m}e^{+ik^{m}_{x}d_{1}}=\sum_{l=-\infty}^{+\infty}\left(a_{1l}\Gamma_{l}e^{-iq^{l}_{x}d_{1}}-a_{2l}\Gamma^{*}_{l}e^{+iq^{l}_{x}d_{1}}\right)J_{m-l}(\alpha_{2}).$
(47)

For the second barrier we have $\Psi_{3}(d_{1},y,t)=\Psi_{4}(d_{1},y,t)$ and
$\Psi_{4}(d_{2},y,t)=\Psi_{5}(d_{2},y,t)$, which give

$\displaystyle
c_{1m}e^{ik^{m}_{x}d_{1}}+c_{2m}e^{-ik^{m}_{x}d_{1}}=\sum_{l=-\infty}^{+\infty}\left(b_{1l}e^{iq^{l}_{x}d_{1}}+b_{2l}e^{-iq^{l}_{x}d_{1}}\right)e^{-i(m-l)\beta}J_{m-l}(\alpha_{4})$
(48) $\displaystyle
c_{1m}\gamma_{m}e^{ik^{m}_{x}d_{1}}-c_{2m}\gamma^{*}_{m}e^{-ik^{m}_{x}d_{1}}=\sum_{l=-\infty}^{+\infty}\left(b_{1l}\Gamma_{l}e^{iq^{l}_{x}d_{1}}-b_{2l}\Gamma^{*}_{l}e^{-iq^{l}_{x}d_{1}}\right)e^{-i(m-l)\beta}J_{m-l}(\alpha_{4})$
(49) $\displaystyle
t_{m}e^{ik^{m}_{x}d_{2}}+\mathbb{0}_{m}e^{-ik^{m}_{x}d_{2}}=\sum_{l=-\infty}^{+\infty}\left(b_{1l}e^{iq^{l}_{x}d_{2}}+b_{2l}e^{-iq^{l}_{x}d_{2}}\right)e^{-i(m-l)\beta}J_{m-l}(\alpha_{4})$
(50) $\displaystyle
t_{m}\gamma_{m}e^{ik^{m}_{x}d_{2}}-\mathbb{0}_{m}\gamma^{*}_{m}e^{-ik^{m}_{x}d_{2}}=\sum_{l=-\infty}^{+\infty}\left(b_{1l}\Gamma_{l}e^{iq^{l}_{x}d_{2}}-b_{2l}\Gamma^{*}_{l}e^{-iq^{l}_{x}d_{2}}\right)e^{-i(m-l)\beta}J_{m-l}(\alpha_{4}).$
(51)

We can write all sets of equations in matrix form as follows

$\begin{pmatrix}\delta_{m,0}\\\
r_{m}\end{pmatrix}=\begin{pmatrix}\mathbb{S}_{1,1}&\mathbb{S}_{1,2}\\\
\mathbb{S}_{2,1}&\mathbb{S}_{2,2}\end{pmatrix}\begin{pmatrix}t_{m}\\\
\mathbb{0}_{m}\end{pmatrix}=\mathbb{S}\begin{pmatrix}t_{m}\\\
\mathbb{0}_{m}\end{pmatrix}$ (52)

where $\mathbb{S}$ is the total transfer matrix. The transfer matrix
$\mathbb{S}$ has infinite order to simplify we go to finite order, so we take
$m$ included between -$N$ and $N$, knowing that $N>\alpha_{j}$ [22], whit
$\mathbb{S}$ is the matrix product of the transfer matrix for each region
(53), and $\mathbb{S}(j,j+1)$ the transfer matrix of the region $j$ to region
$j+1$

$\mathbb{S}=\mathbb{S}(1,2).\mathbb{S}(2,3).\mathbb{S}(3,4).\mathbb{S}(4,5)$
(53)

with

$\displaystyle\mathbb{S}(1,2)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\mathbb{G}^{0}_{-d_{2}}&\mathbb{G}_{d_{2}}\\\
\mathbb{R}_{+,-d_{2}}^{+1}&\mathbb{R}_{-,d_{2}}^{-1}\\\
\end{pmatrix}^{-1}\begin{pmatrix}\mathbb{H}_{+,-d_{2}}^{2}&\mathbb{H}_{+,d_{2}}^{2}\\\
\mathbb{L}_{+,-d_{2}}^{2,+1}&\mathbb{L}_{-,d_{2}}^{2,-1}\\\ \end{pmatrix}$
(54) $\displaystyle\mathbb{S}(2,3)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\mathbb{H}_{+,-d_{1}}^{2}&\mathbb{H}_{+,d_{1}}^{2}\\\
\mathbb{L}_{+,-d_{1}}^{2,+1}&\mathbb{L}_{-,d_{1}}^{2,-1}\end{pmatrix}^{-1}\begin{pmatrix}\mathbb{G}_{-d_{1}}&\mathbb{G}_{d_{1}}\\\
\mathbb{R}_{+,-d_{1}}^{l,+1}&\mathbb{R}_{-,d_{1}}^{l,-1}\\\ \end{pmatrix}$
(55) $\displaystyle\mathbb{S}(3,4)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\mathbb{G}_{d_{1}}&\mathbb{G}_{-d_{1}}\\\
\mathbb{R}_{+,d_{1}}^{+1}&\mathbb{R}_{-,-d_{1}}^{-1}\\\
\end{pmatrix}^{-1}\begin{pmatrix}\mathbb{K}^{4}_{+,-d_{1}}&\mathbb{K}^{4}_{+,d_{1}}\\\
\mathbb{F}^{4,+}_{+,d_{1}}&\mathbb{F}^{4,+}_{-,-d_{1}}\\\ \end{pmatrix}$ (56)
$\displaystyle\mathbb{S}(4,5)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\mathbb{K}^{4}_{+,d_{2}}&\mathbb{K}^{4}_{+,-d_{2}}\\\
\mathbb{F}^{4,+}_{+,d_{2}}&\mathbb{F}^{4,+}_{-,-d_{2}}\\\
\end{pmatrix}^{-1}\begin{pmatrix}\mathbb{G}_{d_{2}}&\mathbb{G}_{-d_{2}}\\\
\mathbb{R}_{+,d_{2}}^{+1}&\mathbb{R}_{-,-d_{2}}^{-1}\\\ \end{pmatrix}$ (57)

where $z=\pm d_{1,2}$ and we have set

$\displaystyle(\mathbb{G}^{0}_{z})=e^{ik^{0}_{x}z}$ (58)
$\displaystyle(\mathbb{G}_{z})_{m,l}=e^{ik^{l}_{x}z}\delta_{m,l}$ (59)
$\displaystyle(\mathbb{R}_{\pm,z}^{\pm 1})_{m,l}=\pm\gamma_{l}^{\pm
1}e^{ik^{l}_{x}z}\delta_{m,l}$ (60)
$\displaystyle(\mathbb{H}^{j}_{\pm,z})_{m,l}=\pm
e^{iq^{l}_{x}z}J_{m-l}(\alpha_{j})$ (61)
$\displaystyle(\mathbb{L}_{\pm,z}^{j,\pm 1})_{m,l}=\pm\Gamma_{l}^{\pm
1}e^{iq^{l}_{x}z}J_{m-l}(\alpha_{j})$ (62)
$\displaystyle(\mathbb{K}^{j}_{\pm,z})_{m,l}=\pm
e^{iq^{l}_{x}z}J_{m-l}(\alpha_{j})e^{-i(m-l)\beta}$ (63)
$\displaystyle(\mathbb{F}_{\pm,z}^{j,\pm 1})_{m,l}=\pm\Gamma_{l}^{\pm
1}e^{iq^{l}_{z}}J_{m-l}(\alpha_{j})e^{-i(m-l)\beta}.$ (64)

We assume the electron propagation direction is from left to right with energy
$\varepsilon$, therefore, the transmission coefficient of the m-$th$ band is
written as follows:

$t_{m}=\mathbb{S^{\prime}}\delta_{m,0}$ (65)

with $\mathbb{S^{\prime}}=\mathbb{S}^{-1}_{1,1}$. Recall that $m$ varies from
-$N$ to $N$, that is to say

$\begin{pmatrix}t_{-N}\\\ .\\\ .\\\ t_{-1}\\\ t_{0}\\\ t_{1}\\\ .\\\ .\\\
t_{N}\end{pmatrix}=\mathbb{S^{\prime}}\begin{pmatrix}0\\\ .\\\ .\\\ 0\\\ 1\\\
0\\\ .\\\ .\\\ 0\end{pmatrix}.$ (66)

To simplify, we have limited our study to the three first bands: the central
band corresponding to the energy $E$ and the two first side bands
corresponding to the energy $\varepsilon\pm\varpi$. In this case, we obtain

$\displaystyle t_{-1}=\mathbb{S^{\prime}}[1,2]$ (67) $\displaystyle
t_{0}=\mathbb{S^{\prime}}[2,2]$ (68) $\displaystyle
t_{1}=\mathbb{S^{\prime}}[3,2].$ (69)

On the other hand, the transmission probability can be obtained from the
electric current density, which is also found using the continuous equation
$J=\Psi^{\dagger}\sigma_{x}\Psi$. Thus, the different current densities are
given as follows:

$\displaystyle J^{0}_{\text{inc}}=\gamma_{0}+\gamma^{*}_{0}$ (70)
$\displaystyle J^{m}_{\text{ref}}=r_{m}r_{m}^{*}(\Gamma_{m}+\Gamma^{*}_{m})$
(71) $\displaystyle
J^{m}_{\text{tra}}=t_{m}t_{m}^{*}(\gamma_{m}+\gamma^{*}_{m})$ (72)

where
$\gamma_{m}=s_{m}\frac{k^{m}_{x}+i(k_{y}-m\omega)}{\sqrt{(k^{m}_{x})^{2}+(k_{y}-\varpi)^{2}}}$
and
$\Gamma_{m}=s_{m}\frac{q^{m}_{x}+i(k_{y}-m\omega)}{\sqrt{(q^{m}_{x})^{2}+(k_{y}-m\varpi)^{2}}}$.
The transmission corresponds to the $m$-th side-band, is given by the
following relation

$T_{m}=\frac{\cos\theta_{m}}{\cos\theta_{0}}|t_{m}|^{2}$ (73)

where $\cos\theta_{m}=\frac{k_{x}^{m}}{\sqrt{(k_{x}^{m})^{2}+k_{y}^{2}}}$ and
$k_{x}^{m}=\sqrt{(\varepsilon+m\varpi)^{2}-k_{y}^{2}}$. Finally, the total
transmission probability is given by the sum over all modes

$T=\sum_{m=-N}^{N}T_{m}.$ (74)