# Probability of high intensities of the light wave
propagating in turbulent atmosphere

I.V. Kolokolov, V.V. Lebedev  Landau Institute for Theoretical Physics,
Russian Academy of Sciences,
142432, Chernogolovka, Semenova 1A, Moscow region, Russia;
and National Research University Higher School of Economics,
101000, Myasnitskaya 20, Moscow, Russia

###### Abstract

We examine statistics of fluctuations of the laser beam intensity at its
propagating in turbulent atmosphere. We are interested in relatively large
propagating distances and the remote tail of the probability density function.
The tail is determined by the stretched exponent, we find its index.

The subject of our work is the theoretical study of the physical properties of
a laser beam propagating in a turbulent atmosphere or, more generally, an
electromagnetic wave in a turbulent medium. The main effect that is being
investigated in the framework of this problem is the diffraction of a laser
beam by fluctuations of the refractive index, which are induced by pressure
fluctuations of turbulent pulsations. The fluctuations of the refractive index
are a random field whose properties are described statistically. Therefore,
theoretical predictions of the behavior of the laser beam concern the average
values (or the corresponding probability density functions), which are
obtained by averaging over the statistics of fluctuations.

The problem of laser beam propagation in the atmosphere has a long history.
The basic theoretical results concerning the propagation of a laser beam in a
turbulent atmosphere were obtained back in the sixties and seventies of the
last century, they are summarized in a number of monographs Tat71 ; Goodman ;
Tat75 ; StrohbehnBook1978 ; AndrewsPhillipsBook1998 . Recently, there has been
a revival of interest in this topic Vorontsov2010 ; Vorontsov2011 ;
Lachinova2016 ; Lushnikov2018 ; Fahey2021 , connected mainly with numerical
modeling of the propagation of a laser beam (electromagnetic wave) in a
turbulent medium, which allows obtaining detailed information about the
process.

The known theoretical results describe mainly the typical behavior of a laser
beam in the atmosphere. At the same time, in the equation for the
electromagnetic field, the random refractive index is a multiplicative noise,
which leads to very non-trivial statistical properties of the electromagnetic
field. It can be expected that the probability of rare events when one or
another value (say, the intensity of the electromagnetic field) has an
abnormally large value will be significantly higher than naive estimates based
on the analysis of typical processes. As an example of such behavior, we can
cite the statistics of quantum particles in a random potential (see, for
example, quenched ), which is the multiplicative noise in the Schrödinger
equation for a quantum particle.

We theoretically study the propagation of a monochromatic electromagnetic wave
in an unlimited turbulent medium. All the characteristic scales of the problem
(the size of the wave packet, the length of its propagation) are assumed to
significantly exceed the wavelength, so that the description of the wave in
terms of the complex envelope $\Psi$ is applicable. Due to the large value of
the speed of light, we can assume that the state of the medium does not change
during the propagation of the wave and use a stationary approximation to
describe its envelope, that is, consider $\Psi$ as a function of the
coordinate $z$ in the direction of wave propagation, and a two-dimensional
radius vector ${\bm{r}}$ in the plane transverse to the direction of
propagation. Of course, the state of the turbulent medium varies with time. In
the accepted approximation, the envelope $\Psi(\bm{r},z)$ adiabatically
adjusts to the current state of the medium.

In this paper, we consider the intensity of the electromagnetic wave to be
quite small, so that nonlinear effects (such as self-focusing) are
insignificant. Then the equation for the envelope $\Psi(\bm{r},z)$ in suitable
units of measurement has the form of a two-dimensional Schrödinger equation

$i\partial_{z}\Psi+\nabla^{2}\Psi+\xi\Psi=0,$ (1)

where $\nabla$ is the two-dimensional gradient in the plane transverse to the
direction of propagation, and $\xi({\bm{r}},z)$ is the fluctuating component
of the refractive index. By virtue of the equation (1), the envelope $\Psi$ is
a field that depends on the implementation of $\xi$. The latter varies over
time, leading to variations in the envelope. Therefore, $\Psi$ can be
considered as a random field whose statistical properties can be extracted by
averaging over long times.

The distance traveled by the wave is assumed to be greater than the integral
turbulence scale, whereas the transverse size of the wave packet is assumed to
be smaller than the integral scale. In this case, the refractive index $\xi$
fluctuates rapidly along the direction of wave propagation. We are interested
in the integral characteristics associated with $\xi$, therefore, by virtue of
the central limit theorem, the random field $\xi$ can be assumed to have
Gaussian probability distribution. It is determined by the pair correlation
function of fluctuations in the refractive index:

$\langle\xi(\bm{r}_{1},z_{1})\xi(\bm{r}_{2},z_{2})\rangle=\delta(z_{1}-z_{2})\left[\mathrm{const}-r_{12}^{c}\right],$
(2)

where angle brackets mean averaging over implementations of the state of the
medium, $r_{12}=|\bm{r}_{1}-\bm{r}_{2}|$ and $c$ is a number. The expression
(2) is valid for distances $r_{12}$ lying in the inertial interval of
turbulence. The constant in the ratio (2) is determined by turbulent
fluctuations on the integral scale, and the power correction is determined by
fluctuations on the scales $\sim r_{12}$. Fluctuations of the refractive index
$\xi$ in a turbulent medium are determined mainly by pressure fluctuations.
For the Kolmogorov spectrum, the exponent of $c$ in the expression (2) is
$c=5/3$. Further, the exponent $c$ is treated as an arbitrary number lying in
the interval $1<c<2$.

We focus mainly on the statistical properties of the intensity $I$ of the
laser beam (of the electromagnetic wave) at some observation point, which we
choose as the origin: $I=|\Psi({\bm{r}}=0,z)|^{2}$. The distance $z$ traveled
by the beam is assumed to be large, $z\gg 1$ in our units. This means that the
effects associated with random diffraction are strong. As the initial state,
we choose a plane wave, $\Psi(\bm{r},0)=1$. A generalization of our
calculation scheme to other cases, say, to the initial Gaussian beam, will be
published elsewhere. We will be interested in the moments of $\langle
I^{n}\rangle$ for large values of $n$. These values determine the probability
of events with an abnormally high intensity value.

The average value $\langle I^{n}\rangle$ can be expressed as $\langle
I^{n}\rangle=F_{2n}(\bm{0},\dots,\bm{0},z)$, where $F_{2n}$ is the $2n$-point
correlation function of the envelop:

$\displaystyle F_{2n}(\bm{r}_{1},\dots,\bm{r}_{2n},z)=\qquad$
$\displaystyle\langle\Psi(\bm{r}_{1},z)\dots\Psi(\bm{r}_{n},z)\Psi^{\star}(\bm{r}_{n+1},z)\dots\Psi^{\star}(\bm{r}_{2n},z)\rangle.$
(3)

This correlation function is represented as a convolution with the Green
function ${\mathcal{G}}_{2n}$:

$\displaystyle F_{2n}=\int d^{2}x_{1}\dots
d^{2}x_{2n}{\mathcal{G}}_{2n}\Psi_{in}(\bm{x}_{1})\dots\Psi_{in}^{\star}(\bm{x}_{2n}).$
(4)

Here $\Psi_{in}(\bm{r})=\Psi({\bm{r}},z=0)$ is the initial value of the
envelop, and the Green function

${\mathcal{G}}_{2n}={\mathcal{G}}_{2n}(\bm{r}_{1},\dots,\bm{r}_{2n},\bm{x}_{1},\dots,\bm{x}_{2n},z),$
(5)

depends on $z$ and $4n$ radius-vectors.

The Green function can be represented as a path integral, which is derived
from the Schrödinger equation (1). After averaging over the fluctuations of
the refractive index in accordance with the expression (2), we arrive at an
integral over the variables $\bm{y}_{j}$, which are functions of the
coordinate $\zeta$, $0<\zeta<z$:

$\displaystyle{\mathcal{G}}_{2n}=\int\prod_{j=1}^{2n}D\bm{y}_{j}\exp\left\\{\int_{0}^{z}d\zeta\,\left(\frac{i}{4}K+W\right)\right\\},$
(6) $\displaystyle
K=\left(\frac{d\bm{y}_{1}}{d\zeta}\right)^{2}+\dots+\left(\frac{d\bm{y}_{n}}{d\zeta}\right)^{2}\qquad$
$\displaystyle-\left(\frac{d\bm{y}_{n+1}}{d\zeta}\right)^{2}-\dots-\left(\frac{d\bm{y}_{2n}}{d\zeta}\right)^{2},\qquad$
(7) $\displaystyle
W={\sum_{i=1}^{n}}\sum_{j=i+1}^{n}y_{ij}^{c}+\sum_{i=n+1}^{2n}\sum_{j=i+1}^{2n}y_{ij}^{c}-\sum_{i=1}^{n}\sum_{j=n+1}^{2n}y_{ij}^{c}.$
(8)

One takes arguments $\bm{x}_{j}$ and $\bm{r}_{j}$ of the Green function as
boundary conditions for the trajectories $\bm{y}_{j}(\zeta)$ at $\zeta=0$ and
$\zeta=z$, see Eq. (5). The “potential” $W$ (8) is the function of the
variables $y_{ij}=|\bm{y}_{i}-\bm{y}_{j}|$. The constant appearing in the
expression (2) falls out of consideration, as it should be.

It is possible to find an explicit expression for the pair Green function
${\mathcal{G}}(\bm{r}_{1},\bm{r}_{2},\bm{x}_{1},\bm{x}_{2},z)$:

$\displaystyle{\mathcal{G}}=\frac{\theta(z)}{16\pi^{2}z^{2}}\exp\left[\frac{i}{2z}(\bm{r}-\bm{x})(\bm{R}-\bm{X})\right.$
$\displaystyle\left.-{z}\int_{0}^{1}d\chi\,|\chi\bm{x}+(1-\chi)\bm{r}|^{c}\right],$
(9)

which determines the behavior of the pair correlation function. Here the
designations $\bm{R}=(\bm{r}_{1}+\bm{r}_{2})/2$,
$\bm{r}=\bm{r}_{1}-\bm{r}_{2}$, $\bm{X}=(\bm{x}_{1}+\bm{x}_{2})/2$,
$\bm{x}=\bm{x}_{1}-\bm{x}_{2}$ are introduced.

For the initial state in the form of a plane wave, when $\Psi_{in}=1$, the
pair correlation function has the simple form TatarskiiJETP1969 ; ZKT77 :

$F(\bm{r}_{1},\bm{r}_{2},z)=\int
d^{2}x\,d^{2}X{\mathcal{G}}=\exp(-z|\bm{r}_{1}-\bm{r}_{2}|^{c}).$ (10)

Thus, the characteristic distance between the points $1$ and $2$ at a distance
$z$ is $r_{ph}=z^{-1/c}\ll 1$. This value makes sense of the phase corruption
length due to fluctuations of the refractive index. The characteristic value
of the variable $X$ is estimated as $z^{1/c+1}$.

The analysis of the $2n$-point Green function (5) shows that for $z\gg 1$ it
has sharp maxima on configurations where there are $n$ pairs of close points
$\bm{x}_{j}$ located at distances of the order of $r_{ph}$ separated by much
larger distances ZKT77 ; KLL20 . In this case, the first point of the pair is
taken from the set $\bm{x}_{1},\dots,\bm{x}_{n}$, and the second is taken from
the set $\bm{x}_{n+1},\dots,\bm{x}_{2n}$. The trajectories $\bm{y}_{j}(\zeta)$
in the path integral (6) starting from close points $\bm{x}_{j}$ also remain
close. In this case, the “potential” $W$ (8) can be approximated by the sum of
$y_{ij}^{c}$ for the pairs, and the remaining summands in the sum (8) mutually
reduce each other. Then the Green function is factorized. For example, for
close couples $(\bm{x}_{1},\bm{x}_{n+1})$, …$(\bm{x}_{n},\bm{x}_{2n})$

$\displaystyle{\mathcal{G}}_{2n}={\mathcal{G}}(\bm{r}_{1},\bm{r}_{n+1},\bm{x}_{1},\bm{x}_{n+1},z)$
$\displaystyle\dots\mathcal{G}(\bm{r}_{n},\bm{r}_{2n},\bm{x}_{n},\bm{x}_{2n},z).$
(11)

Obviously, there is $n!$ of similar contributions, by the number of splits of
$\bm{x}_{j}$ into pairs.

Thus, the space of integration by the initial coordinates in the expression
for the moments of intensity

$\displaystyle\langle I^{n}\rangle=\int d^{2}x_{1}\dots
d^{2}x_{2n}{\mathcal{G}}_{2n}(\bm{0},\dots,\bm{0},\bm{x}_{1}\dots\bm{x}_{2n}).$
(12)

for the case of an initial plane wave is divided into $n!$ areas corresponding
to different point splits into pairs. The contributions of all such regions
are the same and in the main approximation for intensity moments we get
$\langle I^{n}\rangle=n!$, which means the exponential probability density
function $P(I)=\exp(-I)$. Note that $P(I)=\exp(-I)$ corresponds to the
Gaussian statistics of the $\Psi$ field, natural for a complex field with a
random phase.

It is demonstrated in the works ZKT77 ; KLL20 ; Char94  that corrections to
the approximation (11) related to the discarded terms in the expression (8)
for $W$ are proportional to the parameter $\alpha$:

$\alpha=z^{-a},\quad a=\frac{4}{c}-c,$ (13)

it is small at $z\gg 1$. For the Kolmogorov spectrum of turbulence, the
exponent $a$ is equal to $a=11/15$. Corrections to the value $\langle
I^{n}\rangle=n!$ grow as $n$ increases and become essential when $\alpha n\sim
1$. In this paper, we determine the value of $\langle I^{n}\rangle$ at $\alpha
n\gg 1$ and demonstrate that they significantly exceed $n!$. In other words,
at $I\gg\alpha^{-1}$, the tail of the probability density function $P(I)$
arises, significantly exceeding $\exp(-I)$.

If $\alpha n\gg 1$ then the integration space in the integral (12) is still
divided into $n!$ regions, corresponding to well separated pairs of the
trajectories $\bm{y}_{j}(\zeta)$, the regions give identical contributions to
$\langle I^{n}\rangle$. Let us take for definiteness the following subsetting
into the pairs: ($\bm{y}_{1},\bm{y}_{n+1}$), …($\bm{y}_{n},\bm{y}_{2n}$).
Introducing the variables $\bm{Y}_{j}=(\bm{y}_{j}+\bm{y}_{j+n})/2$ and
$\bm{\rho}_{j}=\bm{y}_{j}-\bm{y}_{j+n}$, where $j$ runs through values from
$1$ to $n$, we get a functional representation for the Green function with
arguments in the selected area in the form:

$\displaystyle{\mathcal{G}}=\int\prod_{j=1}^{n}D\bm{Y}_{j}\,D\bm{\rho}_{j}\,\exp\left(-{\mathcal{S}}\right).$
(14)

In the path integral (14), some non-zero initial values are implied
$\bm{Y}_{j}(0)=\bm{X}_{j}$ and zero final values (at $\zeta=z$) are implied
for the trajectories $\bm{Y}_{j}$ and $\bm{\rho}_{j}$, the values are dictated
by examining $\langle I^{n}\rangle$.

As it follows from the representation (6), the action $\mathcal{S}$ in Eq.
(14) is equal to

$\displaystyle{\mathcal{S}}=-\frac{i}{2}\sum_{j}\int_{0}^{z}d\zeta\,\frac{d\bm{Y}_{j}}{d\zeta}\frac{d\bm{\rho}_{j}}{d\zeta}+\sum_{j}\int_{0}^{z}d\zeta\,\rho_{j}^{c}$
$\displaystyle+\sum_{j>k}\int_{0}^{z}d\zeta\,U(\bm{Y}_{j}-\bm{Y}_{k},\bm{\rho}_{j},\bm{\rho}_{k}),$
(15)

In this expression, the indices $j,k$ run through the values from $1$ to $n$.
The quantity $U$ at $\bm{\rho}_{j}$, small in comparison with the
characteristic values of $\bm{Y}_{j}-\bm{Y}_{k}$, is determined by the formula
following from the expression (8):

$\displaystyle
U(\bm{R},\bm{\rho}_{1},\bm{\rho}_{2})\approx\rho_{1,\alpha}\rho_{2,\beta}V_{\alpha\beta}(\bm{R}),$
(16) $\displaystyle
V_{\alpha\beta}(\bm{R})=cR^{c-2}\left[\delta_{\alpha\beta}+(c-2)\frac{R_{\alpha}R_{\beta}}{R^{2}}\right].$
(17)

Here the Greek indices $\alpha,\beta$ designate the components of the vectors
in a plane transverse to the direction of wave propagation.

As we will see below, the values of $\bm{Y}_{j}$ in the integral (14) are
parametrically large by $n$. This allows us to calculate the integral over
$\bm{Y}_{j}$ in the expression (14) in the saddle-point approximation. The
integral over $\bm{\rho}_{j}$ cannot be calculated in the approximation.

We first investigate the saddle-point value of $\bm{Y}_{j}$, which is
determined by the extremum condition $\delta{\mathcal{G}}/\delta\bm{Y}_{j}=0$.
Substituting here the expression (14), we find the equation

$\frac{i}{2}\frac{d^{2}\bar{\bm{\rho}}_{j}}{d\zeta^{2}}+\sum_{k\neq
j}\frac{\partial
V_{\alpha\beta}}{\partial\bm{Y}_{j}}(\bm{Y}_{j}-\bm{Y}_{k})\bar{\rho}_{j\alpha}\bar{\rho}_{k\beta}=0.$
(18)

Here $\bar{\bm{\rho}}_{j}$ is the average value of ${\bm{\rho}}_{j}$, it is
non-zero because of non-zero initial values $\bm{X}_{j}$ of the trajectories
$\bm{Y}_{j}$. When deriving the equation (18), we replaced the average of the
product of ${\rho}_{j\alpha}{\rho}_{k\beta}$ with the product of averages. The
reason for this is that the fluctuations of $\bm{\rho}_{j}$ are determined
mainly by the second term in the action (15), which is diagonal by $j$.
Therefore, the contribution of fluctuations to the mean
${\rho}_{j\alpha}{\rho}_{k\beta}$, where $j\neq k$, turns out to be
negligible.

In field theory, the average value of the fluctuating field is found from the
condition of the extremum of the so-called quantum effective action (see, for
example, Weinberg ). In the problem under study, this condition is equivalent
to the relation $\delta{\mathcal{G}}/\delta\bar{\bm{\rho}}_{j}=0$, which
determines the average $\bar{\bm{\rho}}_{j}$. Utilizing the expression (14),
we find the equation

$\displaystyle\frac{i}{2}\frac{d^{2}Y_{j\alpha}}{d\zeta^{2}}+\sum_{k\neq
j}V_{\alpha\beta}(\bm{Y}_{j}-\bm{Y}_{k})\bar{\rho}_{k\beta}=0.$ (19)

At deriving the equation we have neglected the contribution of the second term
in the action (15). The reason is that the fluctuations of $\bm{\rho}_{j}$
(which are determined mainly by this term) turn out to be much larger than
$\bar{\bm{\rho}}_{j}$, which is why the dependence of the second term on
$\bar{\bm{\rho}}_{j}$ is weak. Further, we will justify this neglect. Note
that the extremum condition for the initial value of $\bar{\bm{\rho}}_{j}(0)$
leads to the condition $d\bm{Y}_{j}/d\zeta(0)=0$, due to the structure of the
action (15).

The equations (18,19) show that the saddle-point value of the expression (14)
is written as

$\displaystyle\ln{\mathcal{G}}_{sp}=\frac{i}{2}\sum_{j}\int_{0}^{z}d\zeta\,\frac{d\bm{Y}_{j}}{d\zeta}\frac{d\bar{\bm{\rho}}_{j}}{d\zeta}$
$\displaystyle-\sum_{j>k}\int_{0}^{z}d\zeta\,\bar{\rho}_{1,\alpha}\bar{\rho}_{2,\beta}V_{\alpha\beta}(\bm{Y}_{j}-\bm{Y}_{k})$
$\displaystyle=\sum_{j>k}\int_{0}^{z}d\zeta\,V_{\alpha\beta}(\bm{Y}_{j}-\bm{Y}_{k})\bar{\rho}_{j,\alpha}\bar{\rho}_{k,\beta}.$
(20)

The last equality in (20) is obtained after integration in parts, taking into
account the equation (19) and the boundary conditions. The value (20) is
negative because $\bm{Y}_{j}$ are real quantities, whereas the averages
$\bar{\bm{\rho}}_{j}$ are purely imaginary quantities, and the function
$V_{\alpha\beta}(\bm{Y}_{j}-\bm{Y}_{k})$ is positive for the values of the
exponent $c$ considered here: $1<c<2$.

The system of equations (18,19) enables one to evaluate the saddle-point
quantity (20). Assuming that all values of $Y_{j}$ together with their
boundary values have the same order of $Y_{j}\sim X$, as well as the average
$\bar{\rho}_{j}\sim i\bar{\rho}$, we arrive at the estimates

$\displaystyle\bar{\rho}\sim(nz^{2})^{-1}X^{3-c},\quad\ln{\mathcal{G}}_{sp}\sim-\frac{1}{z^{3}}X^{4-c}.$
(21)

Note the absence of dependence on $n$ in the estimation for
$\ln{\mathcal{G}}_{sp}$.

Now we turn to accounting for the fluctuation contribution to
${\mathcal{G}}_{2n}$. Fluctuations of $\bm{\rho}_{j}$ can be estimated as
$z^{-1/c}$, the value is much larger than the average $\bar{\rho}$ (21).
Therefore after shifting $\bm{\rho}_{j}$ by its mean value and shifting
$\bm{Y}_{j}$ by its saddle-point value, the effects associated with the
“interaction potential” (16) can be neglected (estimation of the accuracy see
below). As a result, we arrive at a factorized approximation of the type (11).
Thus, the expression for the Green function is factorized
${\mathcal{G}}_{2n}={\mathcal{G}}_{sp}{\mathcal{G}}_{fl}$, where
${\mathcal{G}}_{sp}$ is determined by the expression (20) and the fluctuation
factor is given by the product of the pair Green functions

$\displaystyle{\mathcal{G}}_{fl}=\prod_{j=1}^{n}\frac{1}{16\pi^{2}z^{2}}\exp\left[-\frac{z}{c+1}|\bm{x}_{j}-\bm{x}_{j+n}|^{c}\right].$
(22)

The expression standing here for the pair Green functions is derived from the
general formula (9) after substituting $\bm{r}=\bm{R}=\bm{0}$ and taking into
account that the initial values of the shifted variables
$\bm{Y}_{j}-\bm{X}_{j}$ are zero.

Returning to the expression (12) for the moments of intensity, we find

$\displaystyle\frac{\langle I^{n}\rangle}{n!}=\int d^{2}x_{1}\dots
d^{2}x_{2n}{\mathcal{G}}_{sp}{\mathcal{G}}_{fl}.$ (23)

The factor ${\mathcal{G}}_{sp}$ depends solely on the coordinates of the
centers of the pairs $\bm{X}_{j}=(\bm{x}_{j}+\bm{x}_{j+n})/2$, therefore,
after integration over differences $\bm{x}_{j}-\bm{x}_{j+n}$ we find

$\frac{\langle I^{n}\rangle}{n!}\sim z^{-2(1+1/c)}\int d^{2}X_{1}\dots
d^{2}X_{n}\,{\mathcal{G}}_{sp}.$ (24)

For a large $n$, the value of the integral (24) is estimated as:

$\displaystyle\ln\frac{\langle I^{n}\rangle}{n!}\sim 2n\ln
X-2n\left(1+\frac{1}{c}\right)\ln z-\frac{C}{z^{3}}X^{4-c},$ (25)

where $C$ is a constant of the order of one. Here $X$ is the characteristic
value of variables $\bm{X}_{j}$ introduced above.

Optimizing the expression (25) by $X$, we find:

$X^{4-c}\sim nz^{3}.$ (26)

The substitution of Eq. (26) into Eq. (25) gives the desired asymptotics of
the high intensity moments:

$\displaystyle\ln\frac{\langle I^{n}\rangle}{n!}\sim\frac{2}{4-c}n\ln
n-\frac{2(4-c^{2})}{(4-c)c}n\ln z.$ (27)

The corresponding asymptotics of the probability density function $P(I)$ has
the form

$\displaystyle\ln P\sim-\frac{1}{\alpha}(\alpha
I)^{\beta},\quad\beta=\frac{4-c}{6-c}.$ (28)

For the Kolmogorov spectrum $\beta=7/13$. Since $\beta<1$, we conclude that
there is a higher probability of large values of $I$ in comparison with the
exponential probability density $\exp(-I)$. Formula (28) is in accordance with
the results of recent numerical simulations Lushnikov . Note that the
expression (28) means a significant non-Gaussianity of the statistics of the
wave field $\Psi$, formed by the chaotic scattering on the turbulent
fluctuations.

The asymptotic law (27) turns to $\langle I^{n}\rangle=n!$ at $\alpha n\sim
1$, or at $z\sim n^{c/(4-c^{2})}$, as it should be. A similar statement is
true for the probability density function, which is matched with the exponent
$\exp(-I)$ at $I\sim 1/\alpha$. Eq. (26) enables one to find the
characteristic size of the region at the front of the initial wave defining
$I^{n}$ at $n\gg 1/\alpha$:

$X\sim(\alpha n)^{1/(4-c)}z^{1/c+1}.$ (29)

At $\alpha n\gg 1$ the quantity (29) is much larger than the estimate
$z^{1/c+1}$ for $X$ obtained in the Gaussian regime. Physically, this means
that to create high intensity values, it is required to collect energy from a
large area of the original wave. The expression (29) implies an estimate for
the average values $\bar{\rho}$ (21)

$\bar{\rho}\sim\frac{1}{nz^{2}}X^{3-c}\sim r_{ph}(n\alpha)^{-1/(4-c)}.$ (30)

These averages are small in comparison with the characteristic amplitude of
fluctuations: $\bar{\rho}\ll r_{ph}$, which justifies the approximation (22).
Note, however, that the saddle-point action (20) is determined precisely by
these averages, and in action (20) this smallness is compensated by a large
number of terms.

The found asymptotic law (27) is bounded from above by $n$. This limitation is
determined by the mean square of the fluctuation component of the interaction
of pairs between themselves that we have discarded. The assessment is fair for
him

$z^{2}n^{2}r_{ph}^{4}X^{2(c-2)}\sim n\left(n^{c/4}\alpha\right)^{4/(4-c)}.$

The smallness of this value in comparison with the main value of the
fluctuation action, estimated as $\sim n$, leads to the condition of
applicability of our approach $n\ll\alpha^{-4/c}$. The inverse limit case
requires a special analysis.

We have theoretically established that when a laser beam propagates in a
turbulent medium at distances where diffraction on random fluctuations of the
refractive index plays a significant role, the probability of abnormally large
intensity fluctuations is significantly higher than estimates made on the
basis of Gaussian statistics of the electromagnetic wave envelope $\Psi$. This
is due to the multiplicativity of the refractive index in the equation for the
envelop and it is a universal property of such stochastic systems. We have
found the form of the probability distribution function of abnormally large
intensity values. Our conclusions are made for the simplest case of the
initial plane wave, although the analysis scheme itself is applicable to other
initial forms of the laser beam. The results of their research will be
published elsewhere. The case of extremely high intensity values requires a
separate analysis. A separate analysis also requires to take into account the
nonlinearity (for example, the effect of self-focusing), which may be
significant for the large intensities of the electromagnetic wave we are
considering.

The work is supported by the scientific program of the National Center for
Physics and Mathematics, the project “Physics of high energy densities”, stage
2023-2025. The authors thank P.M. Lushnikov for useful discussions.

## References

  * (1) V. I. Tatarskii, _The Effects of the Turbulence Atmosphere on Wave Propagation_ , (Jerusalem: Israel Program for Scientific Translations, Jerusalem, 1971),
Translated from: V. I. Tatarskii, Rasprostranenie voln v turbulentnoi
atmosfere, Moskva, Nauka, 1967.

  * (2) J.W.Goodman, Statistical Properties of Laser Speckle Pattern, pp. 9-75, in Laser Speckle and Related Phenomena, ed. J. C. Dainty, Springer-Verlag, Berlin, 1975.
  * (3) M. S. Rytov, Yu. A. Kravtsov, and V. I. Tatarskii, _Principles of Statistical Radiophysics 4: Wave Propagation through Random Media_ , Springer-Verlag, Berlin, 1989.
Translated from: M. S. Rytov, Yu. A. Kravtsov, V. I. Tatarskii, Vvedenie v
statisticheskuyu radiofiziku, chast’ II, Sluchainye polya, Moskva, Fizmatlit,
1978

  * (4) J. W. Strohbehn, ed., _Laser Beam Propagation in the Atmosphere_ Springer, New York, 1978.
  * (5) L. C. Andrews and R. L. Phillips, _Laser Beam Propagation Through Random Media_ , SPIE Press Monograph Series, PM152 (2005).
  * (6) M. A. Vorontsov, G. W. Carhart, V. S. R. Gudimetla, T. Weyrauch, E. Stevenson, S. L. Lachinova, L. A. Beresnev, J. Liu, K. Rehder, and J. F. Riker, Proc. 2010 AMOS Conf., (2010).
  * (7) M. A. Vorontsov, V. R. Gudimetla, G. W. Carhart, T. Weyrauch, S. L. Lachinova, E. Polnau, J. R. L. A. Beresnev, J. Liu, and J. F. Riker, Proc. 2011 AMOS Conf., (2011).
  * (8) S. L. Lachinova and M. A. Vorontsov, J. Opt. 18, 025608 (2016).
  * (9) P.M. Lushnikov, N. Vladimirova, Pis’ma v ZhETF, 108, 611-613 (2018) [JETP Lett., 108, 571-576 (2018)].
  * (10) T.Fahey, M.Islam, A.Gardi, R.Sabatini, Atmosphere 12, 918 (2021).
  * (11) F. Evers, A.D.Mirlin, Rev. Mod. Phys. 80, 1355 (2008).
  * (12) V. I. Tatarskii, ZhETF, 56, 2106 (1969) [Sov. Phys JETP, 29, 1133 (1969)].
  * (13) V.U.Zavorotnyi, V.I.Klyatskin, V.I.Tatarskii, ZhETF 73, 481 (1977) [Sov. Phys. JETP, 46, 252 (1977)].
  * (14) I. Kolokolov, V. Lebedev, and P. Lushnikov, Phys. Rev. E 101, 042137 (2020).
  * (15) M. I. Charnovskii, Waves in Random Media, 4, 243 (1994).
  * (16) S. Weinberg, The quantum theory of fields, v. 2, Cambridge University Press, Cambridge, 2001.
(Russian translation: S. Weinberg, Kvantovye polya, Moskva, FIZMATLIT, 2015)

  * (17) P.M.Lushnikov, private communication. Home page: https://math.unm.edu/ plushnik/


# Terahertz plasmonic resonances in coplanar graphene nanoribbon structures

V. Ryzhii1, C. Tang1,2, T. Otsuji1, M. Ryzhii3, and M. S. Shur4 1Research
Institute of Electrical Communication, Tohoku University, Sendai 980-8577,
Japan
2Frontier Research Institute for Interdisciplinary Studies, Tohoku University,
Sendai 980-8577, Japan
3Department of Computer Science and Engineering, University of Aizu, Aizu-
Wakamatsu 965-8580, Japan
4Department of Electrical, Computer, and Systems Engineering, Rensselaer
Polytechnic Institute, Troy, New York 12180, USA

###### Abstract

We analyze plasmonic oscillations in the coplanar graphene nanoribbon (GNR)
structures induced by the applied terahertz (THz) signals and calculate the
GNR impedance. The plasmonic oscillations in the CNR structures are associated
with the electron and hole inductances and the lateral inter-CNR capacitance.
A relatively low inter-GNR capacitance enables the resonant excitation of the
THz plasmonic oscillations in the CNR structures with long GNRs. The GNR
structures under consideration can be used in different THz devices as the
resonant structures incorporated in THz detectors, THz sources using resonant-
tunneling diodes, photomixers, and surface acoustic wave sensors.

##  I Introduction

Due to the remarkable properties of graphene layers (GLs) (high intrinsic
mobility enabling an elevated electrical conductivity, high thermal
conductivity, transparency to light, and mechanical strength), GLs are very
attractive for different device applications. In particular, graphene micro-
and nanoribbons can be used in the detectors and sources of electromagnetic
radiation, surface acoustic wave sensors, power and flexible devices [1, 2, 3,
4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. The possibility of the terahertz (THz)
plasmonic wave excitation in the graphene structures [2, 3, 4, 7, 8, 9, 10,
15, 16, 17, 18, 19, 20, 21, 22, 23] markedly enhances the functionality of the
graphene-based devices. The structures with GLs partitioned into arrays of
graphene nanoribbons (GNRs) have also been already studied (for example, [24,
25, 26, 27, 28, 29]). In this work, we focused mainly on the plasmonic waves
propagating perpendicular to the GNRs aiming the achievement of the plasmonic
resonances at relatively high frequencies, particularly, in the frequencies
corresponding to the mid-infrared range.

In the present paper, we explore the structures comprising two co-planar GNRs
on a dielectric substrate focusing on their terahertz applications. In such
structures the plasmonic wave number is determined by the GNR length, which is
much longer than their widths. As a result, the perpendicular plasmonic modes
are not excited in the frequency range under consideration (the THz range).

It is assumed that the bias dc voltage $V_{G}$ and the THz signal voltage
$\delta V_{\omega}$ with the frequency $\omega$ are applied between the GNR
edges. The bias voltage results in the formation of the two-dimensional
electron and hole systems, 2DES and 2DHS. The signal voltage can be produced
by the impinging electromagnetic radiation and received by an antenna
integrated with the GNR structure. The THz signals excite the plasmonic
standing waves (plasmonic oscillations) in the GNRs. Examining the plasmonic
properties of the co-planar GNRs, we show that with a relatively low
geometrical inter-GNR capacitance, the plasmonic resonances can be in the THz
range even for fairly long GNRs. This property of the coplanar GNRs, as well
as their lateral periodic arrays, can find applications in different devices
using the THz plasmonic resonant structures.

##  II Equations of the model

Figure 1: Schematic top view of the co-planar GNR structures with different
types of voltage bias. Figure 2: Plasmonic frequency $\Omega/2\pi$ versus bias
voltage $V_{G}$ for GNR structures width different length $2H$ , and widths
$2L_{G}$ ($a=L_{G}/L=2$). Figure 3: Frequency dependences of (a, b)) the real
part Re $Y_{\omega}^{P}$, and the imaginary part Im $Z_{\omega}^{P}$ of the
P-GNR admittance and (c, d) the real part Re $Z_{\omega}^{P}$, and the
imaginary part Im $Z_{\omega}^{P}$ of the P-GNR impedance for different
plasmonic frequencies $\Omega/2\pi$ ( $\nu 1.0=$ ps-1, $a=L_{G}/L=2$). Figure
4:  The same as in Fig. 3, but for the S-GNRs.

Considering the co-planar GNR structures, we assume two types of signal
voltage bias, P-GNR and S-GNR, as shown in Fig. 1(a) and Fig. 1(b). The
lengths and widths of the GNRs are $2H$ and $2L_{G}$, respectively, whereas
the spacing between the GNRs is equal to $2L$ being much smaller than other
geometrical parameters ($L_{G},L\ll H$ ).

The GNR potentials are
$\varphi(z,t)^{\pm}=\delta\varphi_{\omega}^{\pm}(z)e^{-i\omega t}$, where the
signs ”+” and ”-” correspond to the GNRs to which the voltages of different
polarities are applied. The displacement current density between the GNRs and
the continuity equations for the electrons and holes governing their transport
along the GNRs can be presented as

$\displaystyle\delta
j_{\omega}=-i\omega\,C(\delta\varphi_{\omega}^{+}-\delta\varphi_{\omega}^{-}).$
(1) $\displaystyle
2L_{G}\sigma_{G,\omega}\frac{d^{2}\delta\varphi_{\omega}^{\pm}}{dz^{2}}=\mp
i\omega\,C(\delta\varphi_{\omega}^{+}-\delta\varphi_{\omega}^{-}).$ (2)

Here

$\displaystyle\sigma_{G,\omega}=\sigma_{G}\frac{i\nu}{(\omega+i\nu)}$ (3)

and

$\displaystyle C=\frac{\kappa}{2\pi^{2}}c$ (4)

are the GNR ac conductivity and the geometrical inter-GNR geometrical
capacitance, respectively, where $\sigma_{G}=e^{2}\mu/\pi\hbar^{2}\nu$ is the
dc Drude conductivity of a GNR with the carrier Fermi energy $\mu$ and the
carrier collision frequency $\nu$ (scattering by impurities and acoustic
phonons), $\kappa=(\kappa_{S}+1)/2$ and $\kappa_{S}$ are the effective
dielectric constant and the substrate dielectric constant (the dielectric
constant of the media above the structure is assumed to be unity), the factor

$c=a\tan^{-1}\biggl{(}\frac{1}{\sqrt{a^{2}-1}}\biggr{)}+\ln(a+\sqrt{a^{2}-1})$

reflects the GNR blade-like shape [30] (see also Refs. [31, 32, 33, 34, 35])
with parameter $a=L_{G}/L$.

The electron and hole Fermi energy in the pertinent GNR is determined by the
bias voltage $V_{G}$:
$\mu\simeq\hbar\,v_{W}\sqrt{\pi\Sigma_{G}}\simeq\hbar\,v_{W}\sqrt{\kappa\,c\,V_{G}/2\pi\,eL_{G}}$,
where $v_{W}\simeq 10^{8}$ cm/s is the characteristic electron and hole
velocity in graphene, $\Sigma_{G}$ the steady-state electron and hole density,
and $\hbar$ is the Planck constant.

Using Eqs. (1) - (3), we arrive at the following equations governing the
potentials $\delta\varphi_{\omega}^{\pm}$:

$\displaystyle\frac{\partial^{2}\delta\varphi_{\omega}^{+}}{\partial
z^{2}}\pm\frac{(\omega+i\nu)\omega}{s^{2}}(\delta\varphi_{\omega}^{+}-\delta\varphi_{\omega}^{-})=0,$
(5)

$\displaystyle\frac{\partial^{2}\delta\varphi_{\omega}^{-}}{\partial
z^{2}}-\frac{(\omega+i\nu)\omega}{s^{2}}(\delta\varphi_{\omega}^{+}-\delta\varphi_{\omega}^{-})=0,$
(6)

with $s=\sqrt{2e^{2}\mu\,L_{G}/\pi\,C\hbar^{2}}$ being the characteristic
velocity of the plasmonic wave along the GNRs. Depending on the signal voltage
application, the boundary conditions for Eqs. (5) and (6) are:

$\displaystyle\delta\varphi_{\omega}^{\pm}|_{z=-H}=\pm\delta
V_{\omega}/2,\quad\partial\delta\varphi_{\omega}^{\pm}/\partial z|_{z=H}=0.$
(7)

for the situation corresponding to the device corresponding to structure Fig
1(a) P-GNR, and

$\displaystyle\delta\varphi_{\omega}^{\pm}|_{z=\mp H}=\pm\delta
V_{\omega}/2,\qquad\partial\delta\varphi_{\omega}^{\pm}/\partial z|_{z=\pm
H}=0.$ (8)

for the device shown in Fig. 1(b) - interdigital GNR connection S-GNR.

##  III Plasmonic oscillations

Equations (5) and (6) with the boundary conditions given by Eqs.(7) and (8)
for the local signal voltage swing between the GNRs
$\delta\varphi_{\omega}=(\delta\varphi_{\omega}^{+}-\delta\varphi_{\omega}^{-})$
result in

$\displaystyle\delta\varphi_{\omega}^{P}=\delta
V_{\omega}\frac{\cos(\gamma_{\omega}z/H)}{\cos\gamma_{\omega}}$ (9)

and

$\displaystyle\delta\varphi_{\omega}^{S}=\delta
V_{\omega}\frac{\cos(\gamma_{\omega}z/H)}{[\cos\gamma_{\omega}-\gamma_{\omega}\sin\gamma_{\omega}]}.$
(10)

for the P- and S GNR structures, respectively.

Here

$\displaystyle\gamma_{\omega}=\frac{H\sqrt{(\omega+\nu)\omega}}{S}=\frac{\pi\sqrt{(\omega+i\nu)\omega}}{2\Omega},$
(11)

where

$\displaystyle\Omega=\frac{e}{H\hbar}\sqrt{\frac{\pi\mu\,L_{G}}{4C}}\simeq\frac{e}{H\hbar}\sqrt{\frac{\pi^{3}\mu\,L_{G}}{2\kappa\,c}}$
(12)

is the characteristic frequency of the plasmonic oscillation fundamental mode
in the GNRs. Considering the dependence of $\mu$ on $V_{G}$, we obtain

$\displaystyle\frac{\Omega}{2\pi}=\frac{1}{4H}\biggl{(}\frac{2\pi\,e^{3}v_{W}^{2}}{\hbar^{2}}\frac{L_{G}V_{G}}{\kappa\,c}\biggr{)}^{1/4}\propto\frac{1}{H}\biggl{(}\frac{L_{G}V_{G}}{\kappa\,c}\biggr{)}^{1/4}.$
(13)

In particular, for $\kappa_{S}=4$ ($\kappa=2.5$), $\mu=140$ meV,
$L_{G}=0.5~{}\mu$m, $L=0.25~{}\mu$m ($c\simeq 2.5$), $H=(2.0-4.0)~{}\mu$m, Eq.
(15) yields $\Omega/2\pi\simeq(1.0-2.0)$ THz.

Equations (9) and (10) describe the spatial distributions of the self-
consistent potential and carrier densities in the GNRs (along them), i.e., the
plasmonic oscillations (the standing plasmonic waves with the frequency
$\omega$ and the normalized wave number $\gamma_{\omega}$. AS follows from
Eqs. (9) and (10), the relative amplitudes of the ac potential
$|\delta\varphi_{\omega}^{P}|/\delta V_{\omega}$ and
$|\delta\varphi_{\omega}^{S}|/\delta V_{\omega}$ can be markedly larger than
unity when the denominators in right-hand sides of these equations
$\cos\gamma_{\omega}|\ll 1$ and
$\cos\gamma_{\omega}\gamma_{\omega}\sin\gamma_{\omega}|\ll 1$ at certain
signal frequencies (the plasmonic resonant frequencies).

Deriving Eqs. (5),(6), (9), and (10), we disregarded the ac potential
nonuniformity across the GNRs (in the $x$-direction), i.e., neglected the
plasmonic oscillation longitudinal modes with the wave vectors $k_{x}\propto
1/L_{G}$ directed perpendicular to the GNR. Such modes correspond to higher
harmonics. This restricts our consideration by the condition $L_{G}\ll H$
(assumed above).

Figure 2 shows the dependences of the plasmonic frequency $\Omega/2\pi$ on the
bias voltage $V_{G}$ for different GNR geometries calculated using Eq. (13)
with Eq. (4). It is assumed that $\kappa_{S}=4$, $H=10~{}\mu$m, and
$L_{G}/L=2$. In line with Eq. (12), an increase in the bias voltage $V_{G}$
leads to an increase in the characteristic plasmonic frequency $\Omega$ due to
the carrier density $\Sigma_{G}$ and the carrier Fermi energy $\mu$ rise when
$V_{G}$ increases. The plasmonic frequency $\Omega$ is sensitive to the
geometrical parameters $L_{G}$ and $L$, which determine the carrier density
and the GNR ac conductivity. The specifics of the GNR structures under
consideration lead to relatively small inter-GNR capacitance and elevated
values of the plasmonic wave velocity $s$. Consequently, rather high values of
plasmonic frequency $\Omega$ (in the THz range) can be realized at fairly
large lengths of the GNRs.

##  IV GNR structure admittance and impedance

The net inter-GNR displacement current is given by

$\displaystyle\delta J_{\omega}=\int_{-H}^{H}dz\delta j_{\omega},$ (14)

with the current density given by Eq. (9) or Eq. (10). For the admittance of
GNR structure under consideration $Y_{\omega}=\delta J_{\omega}\delta
V_{\omega}$, we obtain (for the P- and S-devices, respectively)

$\displaystyle
Y_{\omega}^{P}=-i\frac{2HC{\omega}}{\gamma_{\omega}\,\cot\gamma_{\omega}},\qquad
Y_{\omega}^{S}=-i\frac{2HC\omega}{\gamma_{\omega}(\cot\gamma_{\omega}-\gamma_{\omega})}.$
(15)

The impedances of the P-GNRs and S-GNRs are given by

$\displaystyle Z_{\omega}^{P}=\frac{1}{Y_{\omega}^{P}},\qquad
Z_{\omega}^{S}=\frac{1}{Y_{\omega}^{S}}.$ (16)

At low signal frequencies ($\omega\ll\nu,\Omega$) the admittances and
impedances of both GNR structures tend to zero and infinity, respectively,
provided the inter-GNR leakage current is insignificant. In the range of
elevated frequencies comparable with the characteristic plasmonic frequency
$\Omega$, the admittances and impedances exhibit the oscillations associated
with the plasmonic response of the carriers in the GNRs.

Figures 3 and 4 show the P-GNRs and S-GNRs spectral characteristics.

As follows from the analysis of the frequency characteristics of the P-GNRs
described by the first Eq. (15), Re $Y_{\omega}^{P}$ exhibits maxima at
$\omega=\omega_{n}^{P}\simeq(2n-1)\pi\Omega/2$, where $n=1,2,3,...)$ is the
plasmonic resonance index, i.e., when $\cot\gamma_{\omega_{n}^{P}}\ll 1$. At
these frequencies, the peak height is proportional to the plasmonic
oscillations quality factor $Q=4\Omega/\pi\nu\gg 1$ being weakly dependent on
the resonance index. The first maxima of the curves in Fig. 3(a) correspond to
$\omega_{1}/2\pi\simeq=0.72$ and 1.09 THz (dashed and solid lines). The second
maximum for $\Omega/2\pi=0.72$ THz corresponds to $\omega_{2}^{P}/2\pi\simeq
3\Omega/2\pi=2.16$ THz. Comparison of Figs. 3(a) and 3(b) shows that the
maxima of Re $Y_{\omega}^{P}$ at $\omega=\omega_{n}^{P}$ correspond to a sharp
change in the sign of Im $Y_{\omega}^{P}$ (so that Im $Y_{\omega_{n}}^{=}0$).
In contrast, Re $Z_{\omega}^{P}$ versus $\omega/2\pi$ dependences, given by
the first Eq. (16), exhibit maxima at $\omega=2\omega_{n}^{P}$ [see Fig.3(c)],
i.e., at $\omega/2\pi\simeq(2n-1)\Omega/\pi$ ($\tan\gamma_{2\omega_{n}^{P}}\ll
1$) with the peak height proportional to $Q$. Simultaneously, Im
$Z_{\omega}^{P}$ turns to zero when $\omega=2\omega_{n}^{P}$. This is in line
with the plots in Fig. 3(d).

The frequency dependences of the S-GNR admittance shown in Fig. 4 differ from
those for the P-GNRs. This is due to different conditions of the plasmonic
resonances. Indeed, considering the second Eq. (15), one can see that the
maxima of Re $Y_{\omega}^{S}$ correspond to resonant frequencies
$\omega=\omega_{n}^{S}$, which satisfy the condition
$\cot\gamma_{\omega_{n}^{S}}=\gamma_{\omega_{n}^{S}}$. This condition yields:
$\omega_{0}^{S}\simeq 1.72\Omega/\pi$ (”soft” mode) and $\omega_{n}^{S}\simeq
2\Omega(n+1/n\,\pi^{2})$ with $n=1,2,...$. As seen from Figs. 4(a) and 4(b),
the positions of the maximum of Re $Y_{\omega}^{S}$ and the point, where Im
$Y_{\omega}^{S}=0$ are close to $\omega_{0}^{S}/2\pi\simeq 0.4$ (for
$\Omega/2\pi=0.72$ THz) and to $\omega_{0}^{S}/2\pi\simeq 0.4$THz (for
$\Omega/2\pi=1.09$ THz). As can be derived using the second Eq. (15), the
height of the Re $Y_{\omega}^{S}$ maxima is proportional to
$(4\Omega/\pi\nu)/[2+(\pi\omega_{n}^{S}/2\Omega)^{2}]=Q/(2+\pi^{2}n^{2})$,
demonstrating a substantial decrease with increasing plasmonic index $n$ and
the fact that the peak height is generally lower than the pertinent peaks of
the P-GNR admittance. Such a dependence of Re $Y_{\omega}^{S}$ on $n$ is
clearly seen in Fig. 4(a). In particular, the ratio
$\biggl{(}{\mathrm{R}e}~{}Y_{\omega_{1}^{S}}^{S}\biggr{/}{\mathrm{R}e}Y_{\omega_{0}^{S}}^{S}\biggr{)}\simeq
0.23$.

It is interesting that in both P- and S-GNRs the Re $Z_{\omega}^{P}$ and
$Z_{\omega}^{S}$ exhibit maxima at approximately the same frequencies having
approximately the same height (compare Figs. 3(c) and (4c).

One notes that the resonant frequencies at which Im $Y_{\omega}^{P}$, Im
$Y_{\omega}^{S}$, Im $Z_{\omega^{P}}$, and Im $Z_{\omega}^{S}$ equal to zero
correspond to very low values of Re $Z_{\omega}^{P}$ and Re $Z_{\omega}^{S}$.

The sharpness of the resonant peaks substantially depends on the carrier
collision frequency $\nu$. An increase in the ratio $\nu/\Omega$ (i.e., the
inverse plasmonic oscillations quality factor $Q^{-1}$) leads to the pertinent
lowering of the resonant peak heights. This implies that the realization of
the GNR structures exhibiting pronounced plasmonic resonances requires using
sufficiently perfect GNRs or lowering the operation temperature down from room
temperature.

##  V Discussion and comments

The structures based on the vertically stacked separately contacted GLs and
GNRs with the hBN, WS2, Al2O3, and other insulating layers were studied and
used in different devices (optical and infrared modulators and detectors, THz
detectors, THz detectors and photomixers, and others[36, 37, 38, 39, 40, 41,
42, 43, 44, 45]. The resonant plasmonic effects in such structures can lead to
a substantial enhancement of device characteristics. The resonant plasmonic
response in the lateral co-planar structures with separately contacted GNRs,
considered above, is akin to that in the vertically stacked multiple-GL
structures [40]. However, as predicted above, the lateral coplanar GNR
structures can exhibit a pronounced plasmonic response controlled by the bias
voltage even when the GNRs are relatively long (about several micrometers).
This is because the plasmonic frequency in such structures is determined by
the inter-GNR geometrical capacitance. This capacitance can be much smaller
than the capacitance of the vertically stacked GLs. Indeed, the geometrical
capacitance (per unit of the GNR length) of the lateral structure with two
coplanar GNRs considered above is given by $C=(\kappa/2\pi^{2})c$, whereas the
geometrical capacitance of the vertically stacked GNRs of the same lateral
sizes separated by the spacing $W$ is equal $C_{V}=(\kappa_{S}L_{G}/2\pi\,W)$,
hence

$\displaystyle\frac{C}{C_{V}}=\frac{2cW}{\pi\,L_{G}}\frac{\kappa}{\kappa_{S}}=\frac{2cW}{\pi\,L_{G}}\frac{(\kappa_{S}+1)}{2\kappa_{S}}.$
(17)

Comparing the plasmonic frequencies of co-planar and vertically stacked GNRs,
$\Omega$ and $\Omega_{V}$, and assuming that GNR lengths are equal to $2H$ and
$2H_{V}$, respectively, we obtain

$\displaystyle\frac{\Omega}{\Omega_{V}}\biggr{|}_{\mu}=\frac{H_{V}}{H}\biggl{(}\frac{\pi
L_{G}}{2cW}\frac{2\kappa_{S}}{(\kappa_{S}+1)}\biggr{)}^{1/2}$ (18)

if the carrier Fermi energy $mu$ (carrier density) in the GNR structures is
the same, and

$\displaystyle\frac{\Omega}{\Omega_{V}}\biggr{|}_{V_{G}}=\frac{H_{V}}{H}\biggl{(}\frac{\pi
L_{G}}{2cW}\frac{2\kappa_{S}}{(\kappa_{S}+1)}\biggr{)}^{1/4}$ (19)

for equal bias voltages, $V_{G}$, applied between the GNRs. Considering that
in the devices with the vertically stacked GNRs the thickness of the inter-GNR
dielectric layer $W$ is relatively small and setting $2L_{G}=0.75~{}\mu$m,
$L=0.375~{}\mu$m ($c\simeq 2.36$), $W=(0.01-0.03)~{}\mu$m, and $\kappa_{S}=4$,
we find $(\Omega/\Omega_{V})|_{\mu}\simeq(4.1-7.1)H_{V}/H$. This implies that
in the coplanar GNR structures the same plasmonic frequency can be realized at
much longer GNRs ($H\gg H_{V}$).

As has been shown previously, the plasmonic resonances in the structures with
the arrays of several vertically stacked GNRs (GNR stacks) can be more
pronounced [17, 18, 27]. The results obtained above can be applied for the
lateral arrays vertically stacked GNRs the proper renormalization of the
characteristic plasmonic frequency $\Omega$ given by Eq. (12). In the case
when the GNR stack comprises $K$ non-Bernal stacked GNRs, the electron and
hole density in each GNR is equal to $\Sigma_{G,K}=C\kappa\,V_{G}/4\pi\,e\,K$.
Hence, for the carrier Fermi energy $\mu_{K}$, which determines the GNR Drude
conductivity, we obtain $\mu_{K}\propto\sqrt{\Sigma_{G,K}}\propto 1/\sqrt{K}$.
Accounting for that the net GNR stack conductivity is proportional to $K$, for
the renormalized plasmonic frequency $\Omega_{K}$ we find
$\Omega_{K}\propto\sqrt{K}$.

In principle, the plasmonic waves along the lateral depleted channel formed in
the graphene p-n- junction predicted previously [46] resemble those considered
above. However, there is a substantial difference between the plasmonic waves
associated with the oscillating charges localized in the isolated GNRs,
considered by us, and the plasmonic waves associated with the residual charges
in the depleted area in graphene lateral p-n junctions.

The voltage-controlled coplanar GNR plasmonic resonators can be used in the
THz detectors exploiting the signal rectification in p-i-n diodes and the THz
radiation sources based on the resonant-tunneling diodes. In the first case,
an enhanced real part of the GNR structures admittance at the plasmonic
resonances promotes the rectified current amplification, which leads to an
elevated detector responsivity. In the second case, low values of the real
part of the GNR resonant impedance, can support the THz oscillation self-
excitation and radiation emission.

Using periodic coplanar S-GNR structures with multiple GNR pairs can simplify
the problem matching of these structures with the external circuits (i.r.,
with a load resistance or a THz antenna). These structures might bealso used
for the THz photomixers exploiting the plasmonic resonances, GNR optical
transparency, and strong light absorption in underlying layers, as well as for
the acoustic and photoacoustic transducers.

## Conclusions

We explored the plasmonic response of the lateral coplanar GNR structures to
the impinging radiation. The pertinent resonant frequencies are determined by
the plasmonic frequencies, which can fall into the THz range in structures
with fairly long GNRs (about several micrometers). Due to the sensitivity of
the plasmonic frequency to the bias voltage, the resonant response can be
effectively voltage-controlled. The coplanar GNR structures and their version
with multiple interdigital GNRs can serve as the resonant cavities for
different passive and active THz devices, including THz detectors, photomixers
using the interband transitions in the structure substrate (photoconductive
antenna), oscillators with the diodes exhibiting negative dynamic conductivity
(for example, resonant-tunneling diodes), and the acoustic transducers.

## Acknowledgments

The work at RIEC, FRIIS, and UoA was supported by the Japan Society for
Promotion of Science (KAKENHI Nos. 21H04546, 20K20349), Japan. The work at RPI
was supported by AFOSR (contract number FA9550-19-1-0355).

## Data availability

All data that support the findings of this study are available within the
article.

## References

  * [1] P. Avouris, Z. Chen, and V. Perebeinos, “Carbon-based electronics,” Nat. Nanotech. 2, 605 (2007).
  * [2] V. Ryzhii, T. Otsuji, and M. Shur, “Graphene based plasma-wave devices for terahertz applications,” Appl. Phys. Lett., 116, 140501 (2020).
  * [3] M. Shur, G. Aizin, T. Otsuji, and V. Ryzhii “Plasmonic field-effect transistors (TeraFETs) for 6G communications,” Sensors, 21, 7907 (2021).
  * [4] T. Otsuji, S.A. Boubanga-Tombet, A. Satou, D. Yadav, H. Fukidome, T. Watanabe, T. Suemitsu, A.A. Dubinov, V.V. Popov, W. Knap, V. Kachorovskii, K. Narahara, M. Ryzhii, V. Mitin, M.S. Shur, and V. Ryzhii, “Graphene-based plasmonic metamaterial for terahertz laser transistors,” Nanophoton. 11,1677 (2022).
  * [5] V. Ryzhii,M. Ryzhii, A. Satou, T. Otsuji, M. Ryzhii, V. Mitin, M. S. Shur, and T. Otsuji, “Ballistic injection terahertz plasma instability in graphene n+-i-n-n+ field-effect transistors and lateral diodes, ” Phys. Stat. Solidi A 219,2100694 (2022).
  * [6] L. Wang, N. An, S. Gong, X. Sheng, Y. Li, B. Yao, C. Yu, Z. He, Q. Liu, Z. Feng, T. Otsuji, and Y. Zhang “Ultrafast terahertz transparency boosting in graphene meta-cavities,” Nanophoton., 11, 4899 (2022).
  * [7] M. Ryzhii, V. Ryzhii, M. S. Shur, V. Mitin, C. Tang, and T. Otsuji, “Terahertz bolometric detectors based on graphene field-effect transistors with the composite h-BN/black-P/h-BN gate layers using plasmonic resonances,” J. Appl. Phys., vol. 134, no. 8, 084501, 2023.
  * [8] V. Ryzhii, C. Tang, T. Otsuji, M. Ryzhii, V. Mitin, and M. S. Shur “Effect of electron thermal conductivity on resonant plasmonic detection in terahertz hot-electron bolometers based on metal/black-AsP/graphene FETs,” Phys. Rev. Appl. 19, 064033( 2023).
  * [9] V. Ryzhii, C. Tang, T. Otsuji, M. Ryzhii, V. Mitin, and M.S. Shur “Hot-electron resonant terahertz bolometric detection in the graphene/black-AsP field-effect transistors with a floating gate,” J. Appl. Phys. 133, 174501 (2023).
  * [10] K. Tamura, C. Tang, D. Ogiura, K. Suwa, H. Fukidome, Y. Takida, H. Minamide, T. Suemitsu, T. Otsuji, and A. Satou, “Fast and sensitive terahertz detection in a current-driven epitaxial-graphene asymmetric dual-grating-gate FET structure,” APL Photonics 7, 16101 (2022).
  * [11] X. Guo, J. An, H. Wu, Z. Cai, and P. Wang, “Graphene acoustic transducers based on electromagnetic interactions,” Ultrasonics 114, 106420 (2021).
  * [12] J. He, Z. Fang, C. Gao, W. Zhang, X. Guan, and J. Lin, “Graphene enhanced flexible piezoelectric transducers for dynamic strain measurement: from material preparation to application,” Smart Mater. Struct. 32, 025012 (2023).
  * [13] X. Guo, Y. Zhang, J. An, Q. Zhang, R. Wang, and X. Yu, “Experimental investigation on characteristics of graphene acoustic transducers driven by electrostatic and electromagnetic forces,” Ultrasonics 127, 1064857 (2023)
  * [14] M. Šiškins, M. Lee, D. Wehenkel, R. van Rijn, T. W. de Jong, J. R. Renshof, B. C. Hopman, W. S. J. M. Peters, D. Davidovikj, H. S. J. van der Zant, and P. G. Steeneken, “Sensitive capacitive pressure sensors based on graphene membrane arrays,” Microsyst. Nanoeng. 6, 102 (2020).
  * [15] V. Ryzhii, A. Satou, and T. Otsuji, “Plasma waves in two-dimensional electron-hole system in gated graphene heterostructures,” J. Appl. Phys. 101, 024509 (2007).
  * [16] H. Yan, X. Li, B. Chandra, G. Tulevski, Y. Wu, M. Freitag, W. Zhu, P. Avouris, and F. Xia, “Tunable infrared plasmonic devices using graphene/insulator stacks,” Nat. Nanotech. 7, 330 (2012).
  * [17] N. K. Emani, D. Wang, T.-F. Chung, L. J. Prokopeva, A. V. Kildishev, V. M. Shalaev, Y. P. Chen, and A. Boltasseva, “Plasmon resonance in multilayer graphene nanoribbons,” Laser $\&$ Photonics Reviews 9, 650 (2015). $2L+2L_{G}=150$ nm.
  * [18] D. Svintsov, Zh. Devizorova, T. Otsuji, and V. Ryzhii, “Plasmons in tunnel-coupled graphene layers: backward waves with quantum cascade gain,” Phys. Rev. B. 94, 115301 (2016).
  * [19] D. Rodrigo, A. Tittl, O. Limaj, F. J. G. de Abajo, V. Pruneri, and H. Altug, “Double-layer graphene for enhanced tunable infrared plasmonics,” Light: Science $\&$ Applications 6, e16277 (2017).
  * [20] G. Li, V. Semenenko, V. Perebeinos, and P. Liu, “Multilayer graphene terahertz plasmonic structures for enhanced frequency tuning range,” ACS Photonics 6, 3180 (2019).
  * [21] S. Ogawa, S. Fukushima, and M. Shimatani, “Graphene plasmonics in sensor applications: A Review,” Sensors (Basel) 20, 3563 (2020).
  * [22] L. Wang, N. An, Xu-s. He, X. Zhang A. Zhu, B. Yao,and Y. Zhang, “Dynamic and active THz graphene metamaterial devices,” Nanomaterials (Basel) 17, 2097 (2022)
  * [23] S. Gong, L. Wang, Y. Zhang, Z. Yang, X. Li, Q. Wen, Z. He, S. Liang, L. Yuan, C. Yu, Z. Feng, Z. Yang, and X. Zhang, “Ultra-extraordinary optical transmission induced by cascade coupling of surface plasmon polaritons in composite graphene–dielectric stack,” Opt.Express 28, 30502 (2020)
  * [24] V. V. Popov, O. V. Polischuk, A. R. Davoyan, V. Ryzhii, T. Otsuji, and M. S. Shur, “Plasmonic terahertz lasing in an array of graphene nanocavities,” Phys. Rev. B 86, 1955437 (2012).
  * [25] P. A. D. Goncalves, E. J. C. Dias, Y. V. Bludov, and N. M. R. Peres, “Modeling the excitation of graphene plasmons in periodic grids of graphene ribbons: An analytical approach,” Phys. Rev. B 94, 195421 (2016).
  * [26] L. Xiong, C. Forsythe, M. Jung, A. S. McLeod, S. S. Sunku, Y. M. Shao, G. X. Ni, A. J. Sternbach, S. Liu, J. H. Edgar, E. J. Mele, M. M. Fogler, G. Shvets, C. R. Dean, and D. N. Basov, “Photonic crystal for graphene plasmons,” Nat. Comm. 10, 4780 (2019).
  * [27] V. Semenenko, M. Liu, and V. Perebeinos, “Tunable graphene plasmons in nanoribbon arrays: the role of interactions,” Opt. Express 11, 1390 (2021).
  * [28] L. Cui, J. Wang, and M. Sun, “Graphene plasmon for optoelectronics,” Reviews in Physics 6, 100054 (2021).
  * [29] Z. Ahmad, S. Oh, and E. Muljarov, “Transverse-electric surface plasmon polaritons in periodically modulated graphene,” arXiv:2301.1027v2 [cond-mat.mes-hall], 2023
  * [30] A. Sh. Achoyan, A. É. Yesayan, É. M. Kazaryan, and S. G. Petrosyan, “Two-dimensional p–n junction under equilibrium conditions,” Semiconductors,36, 9023 (2002).
  * [31] S. G. Petrosyan and A. Ya. Shik, “Contact phenomena in low-dimensional electron systems,” Sov. Phys. JETP 69, 1261 (1989).
  * [32] B. Gelmont,M. Shur, and C. Moglestue, “Theory of junction between two-dimensional electron gas and p-type semiconductor,” IEEE Trans. Electron Devices 39, 1216 (1992).
  * [33] D. B. Chklovskii, B. I. Shklovskii, and L. I. Glazman, “Electrostatics of edge channels,” Phys. Rev. B 46, 4026 (1992).
  * [34] Ginés Lifante Pedrola, Beam Propagation Method for Design of Optical Waveguide Devices, Wiley, 2015
  * [35] V. Ryzhii, M. Ryzhii,V. Mitin, M.S. Shur, and T. Otsuji, “S-shaped current-voltage characteristics of n+-i-n-n+ graphene field-effect transistors due to the Coulomb drag of quasiequilibrium electrons by ballistic electrons,” Phys. Rev. Appl. 16, 014001 (2021).
  * [36] L. Britnell, R. V. Gorbachev, R. Jalil, B. D. Belle, F. Shedin, A. Mishenko, T. Georgiou, M. I. Katsnelson, L. Eaves, S. V. Morozov, N. M. R. Peres, J. Leist, A. K. Geim, K. S. Novoselov, and L. A. Ponomarenko, “Field-effect tunneling transistor based on vertical graphene heterostructures,” Science 335, 947 (2012).
  * [37] T. Georgiou, R. Jalil, B. D. Bellee, L. Britnell, R. V. Gorbachev, S. V. Morozov, Y.-J. Kim, A. Cholinia, S. J. Haigh, O. Makarovsky, and L. Eaves, “Vertical field-effect transistor based on graphene-WS2 heterostructures for flexible and transparent electronics,” Nat. Nanotechnol. 8, 100 (2013).
  * [38] L. Britnell, R. V. Gorbachev, A. K. Geim, L. A. Ponomarenko, A. Mishchenko, M. T. Greenaway, T. M. Fromhold, K. S. Novoselov, and L. Eaves, “Resonant tunnelling and negative differential conductance in graphene transistors,” Nat. Commun. 4, 1794 (2013).
  * [39] M. Liu , X. Yin, and X. Zhang, “Double-layer graphene optical modulator,” Nano. Lett. 12, 1482 (2012).
  * [40] V. Ryzhii, T. Otsuji, M. Ryzhii, V. G. Leiman, S. O. Yurchenko, V. Mitin, M.S. and Shur, “Effect of plasma resonances on dynamic characteristics of double graphene-layer optical modulators,” J. Appl. Phys. 112, 104507 (2012). .
  * [41] C. T. Phare, Y.-Ho D. Lee1, J. Cardenas, and M. Lipson, “Graphene electro-optic modulator with 30 GHz bandwidth,” Nat. Photon. 9, 511 (2015).
  * [42] G. Kovacevic and S. Yamashita, “Design optimization for a high-speed two-layer graphene optical modulator on silicon,” IEICE Electronic Express, 13, 1 (2016).
  * [43] M. A. Giambra, V. Sorianello, V. Miseikis, S. Marconi, A. Montanaro, P. Galli, S. Pezzini, C. Coletti, and M. Romagnoli, “High-speed double layer graphene electro-absorption modulator on SOI waveguide,” Opt. Express 27, 20145 (2019).

  * [44] V. Ryzhii, M. Ryzhii, V. Mitin, M.S. Shur, A.Satou, and T. Otsuji, “Terahertz photomixing using plasma resonances in double-graphene layer structure,” J. Appl. Phys. 113, 174506 (2013).
  * [45] V. Ryzhii, T. Otsuji, V. Ya. Aleshkin, A. A. Dubinov, M. Ryzhii, V. Mitin, and M. S. Shur, “Voltage-tunable terahertz and infrared photodetectors based on double-graphene-layer structures,” Appl. Phys. Lett. 104, 163505 (2014).
  * [46] E. G. Mishchenko, A. V. Shytov, and P.G. Silvestrov, “Guided plasmons in Graphene p-n junctions,” Phys. Rev. Lett. 104, 156806 (2010).