# ScGAN: A Generative Adversarial Network to Predict Hypothetical
Superconductors

Evan Kim<EMAIL_ADDRESS>Tesla STEM High School, Redmond, WA 98053, USA
S.V. Dordevic<EMAIL_ADDRESS>Department of Physics, The University of
Akron, Akron, OH 44325, USA

###### Abstract

Despite having been discovered more than three decades ago, High Temperature
Superconductors (HTSs) lack both an explanation for their mechanisms and a
systematic way to search for them. To aid this search, this project proposes
ScGAN, a Generative Adversarial Network (GAN) to efficiently predict new
superconductors. ScGAN was trained on compounds in OQMD and then transfer
learned onto the SuperCon database or a subset of it. Once trained, the GAN
was used to predict superconducting candidates, and approximately 70% of them
were determined to be superconducting by a classification model–a 23-fold
increase in discovery rate compared to manual search methods. Furthermore,
more than 99% of predictions were novel materials, demonstrating that ScGAN
was able to potentially predict completely new superconductors, including
several promising HTS candidates. This project presents a novel, efficient way
to search for new superconductors, which may be used in technological
applications or provide insight into the unsolved problem of high temperature
superconductivity.

##  I Introduction

In recent years, superconductors have been applied in a variety of important
technologies such as power transmission lines, MRI magnets, Maglev trains, and
quantum computers. Quantum computers are especially important as they are
expected to solve problems that are too computationally expensive for current
classical computers. However, to be used in these applications the
superconductors must be cooled below their critical temperatures ($T_{c}$),
which for most current superconductors are very low. For instance, Google’s
quantum computer must be maintained at $0.02\;\mathrm{K}$, which severely
limits its general use [1]. This points to a growing need for superconductors
with higher $T_{c}$, which has made them an active research topic for the last
couple of decades.

The mechanisms of superconductivity in most materials with relatively high
$T_{c}$ is not fully understood, which means there is no systematic way to
search for new materials or to predict their critical temperatures [2]. Thus,
the current procedure for finding HTSs is essentially trial–and–error, which
is extremely inefficient. This was exemplified in a recent study, which found
that only about $3\%$ of the approximately $1000$ materials surveyed were
superconducting [3]. Furthermore, the study failed to find any superconductors
with $T_{c}>60\;\mathrm{K}$. This extreme inefficiency means that the
likelihood of manually finding new superconductors, especially HTSs, is
extremely low.

To address this difficulty, the use of computational tools in
superconductivity research has become popular in recent years [4]. In
particular, there have been several studies utilizing machine learning to
predict whether a given chemical compound will be superconducting or not. The
earliest such study was done by Stanev  _et al._ in which two random forest-
based models were built: a classification model for predicting
superconductivity and a regression model for predicting superconducting
transition temperature [5]. The models were successfully applied on the
SuperCon database, achieving a 92% accuracy in classification and a
$R^{2}=0.88$ for regression. They did run across one limitation of their
machine learning model, however. When trained on a certain class of
superconductors (e.g. cuprates), the model was unable to make good predictions
on other classes of superconductors (e.g. pnictides).

Following this pioneering study, there have been several other studies
applying machine learning methods, such as a K-nearest neighbors algorithm
(Ref. [6]) and a deep learning model (Ref. [7]). The K-nearest neighbors
algorithm reported improvements on the previous study from Ref. [5] in terms
of overall performance: an $R^{2}$ of $0.93$ and a classification accuracy of
$96.5\%$. The deep learning model, on the other hand, showed that it might be
possible to overcome the limitation that Stanev  _et al._ faced, as they were
able to make predictions about pnictide superconductors from the training data
that did not contain them.

The way each of these previous studies predicted new superconducting materials
was by running their trained model on a database of known existing chemical
compounds and finding which compounds the model indicated could be
superconducting. This procedure has several notable limitations. First, these
studies miss out on the vast chemical composition space that is not already
contained in existing databases. Second, commonly used databases (such as ICSD
and OQMD) contain mostly stoichiometric compounds, whereas many
superconductors are non-stoichiometric (cuprates and pnictides, in
particular), and so many possibilities were missed in that manner as well (see
for example Table 3 in Ref. [5] and Table 1 in Ref. [6]). Finally, the
discovery of superconductors in the lab usually does not happen that way; they
are discovered by synthesizing new materials and testing them, rather than
checking the known ones.

In order to overcome these limitations, in this work we employ Generative
Adversarial Networks (GANs) [8]. Generative models refer to a general class of
machine learning models which are able to generate things that resemble the
input dataset. They have proven to be extremely powerful, and in recent years
have found numerous applications such as science, engineering, medicine, art,
video games, deepfakes, etc. [9]. In most cases GANs performed better than
other generative models, such as variational autoencoders, because they are
able to learn more hidden rules in the input data set. For example, Dan  _et
al._ [10] reported a GAN model which generated new chemical compounds with
92.53% novelty and 84.5% validity. Another work, Hu  _et al._ [11] applied
GANs to the SuperCon dataset, but for the purpose of characterization rather
than prediction.

In this work, we combine the general idea of GANs with the previous
superconductor models, and propose the first GAN to predict new
superconductors (ScGAN). Our ScGANs are based on chemical composition only,
and are able to generate new superconducting materials with high novelty and
validity. The paper is organized as follows. In Section II we present the
details of the creation of the ScGANs. In section III the main results of the
study are discussed. In particular, we present a list of hypothetical
superconducting materials generated by our ScGANs, as well as their predicted
critical temperatures. Finally, in Section IV we summarize the most important
findings.

##  II Methodology

As stated in the introduction, we chose the GAN as our generative model. Its
structure is shown in Fig. 1. A GAN is composed of two competing neural
networks, the generator and discriminator. The generator takes in random noise
and generates a “fake” compound, which the discriminator attempts to determine
if it is real or not using existing data of real compounds. Each of them then
updates their parameters based on the performance of the discriminator. The
two networks improve their performance iteratively until a generator can
generate realistic looking compounds, while the discriminator can detect
unrealistic looking ones [8].

Figure 1: The architecture of the GAN model used. Y1.2Ba0.8CuO4 and
O7.25Ca0.62Cu2.13Sr1.91Y0.37Bi1.71 are examples of real (from SuperCon) and
“fake” superconductors, respectively.

Fig. 2 depicts the training/testing process of our GANs and it has three main
stages. The main idea is to have the model first learn general chemical
composition rules by training off of a larger dataset (OQMD database) of
general chemical compounds and then transfer learning it onto a (much smaller)
dataset of known superconducting materials (SuperCon). The idea of transfer
learning is to allow a model to learn even with a limited amount of data,
which is the case here as the general compounds dataset (OQMD) has on the
order of $\sim 10^{6}$ data points, while the superconductor datasets will
only have on the order of $\sim 10^{4}$ data points.

Figure 2: The overall training / testing process from the data to the final
GAN model. It is composed of three main stages: (a) data processing, (b)
training on the OQMD dataset, and finally (c) transfer learning onto the
superconductor dataset (the whole SuperCon, or a part of it).

###  II.1 Data Collection

Data were sourced from two open-source databases: SuperCon [12] and the Open
Quantum Materials Database (OQMD) [13, 14]. SuperCon is the largest database
for superconducting materials with around $30,000$ superconductors before
filtering. Similar to what was done in some previous studies [6, 15], we only
used the chemical compositions of the materials extracted from SuperCon. On
the other hand, OQMD is a much larger database with around $10^{6}$ DFT-
calculated compounds, most of which are not superconductors.

###  II.2 Data Processing

Figure 3: Chemical composition data represented as a matrix [6] in $\mathbb{R}^{96\times m}$, where $m$ is the number of datapoints. Each column is a single compound, with each entry representing the number of each element present in the compound. Note that that the numbers in the matrix are for illustration purposes only; they do not represent any real compounds or superconductors. Figure 4: An “adjusted one-hot” encoding of the chemical composition as a matrix in $\mathbb{R}^{96\times 8}$. The idea here is that the amount of each element is encoded through both the vertical location of the yellow box and the value in it, which allows the GAN to learn the chemical compositions better. Class | Quantity | Percentage
---|---|---
Cuprates | 7,304 | 44.4%
Pnictides | 1,436 | 8.7%
Others | 7,749 | 47.0%
Everything | 16489 | 100%
Table 1: Distribution of the superconductors in our filtered version of
SuperCon by class. Each of these served as different training sets. We notice
that the number of pnictide entries is much smaller compared with cuprates and
others.

Before filtering the data, all compounds had to be expressed in a common
format so that unwanted datapoints could be detected regardless of small
differences in formatting across the databases, such as the ordering of
elements. First, all datapoints were formatted as $1\times 96$ matrices [6]
(the $96$ is the maximum atomic number present across all the compounds in the
database) and then combined so that each dataset was a matrix in
$\mathbb{R}^{96\times m}$, where $m$ is the number of data points (see Fig.
3). Then, both datasets were filtered for duplicates, which reduced OQMD to
around $800,000$ datapoints. The SuperCon dataset required further processing
as a number of entries were incomplete, and they were either corrected or
removed, leaving around 16,000 datapoints. Lastly, the SuperCon dataset was
split into three different groups (classes): cuprates, pnictides (iron-based)
and others (anything else that is neither a curpate nor a pnictide). Table 1
has the quantitative counts of these different groups. The purpose of this was
to test the model’s ability to learn the different classes of superconductors,
especially the high temperature classes of cuprates and pnictides.

Once the data was filtered, the chemical compositions were then transformed
into a form for the GAN to train on. A previous study, Ref. [10], used a one-
hot encoding for general chemical compounds. However, that encoding was
designed for stoichiometric compounds, i.e. for integer values of parameters,
and many superconductors are non-stoichiometric, i.e. have decimal
compositions. Instead, we propose the use of a “adjusted one-hot” encoding
that works for decimals, in which each compound is represented by a 96
$\times$ 8 matrix of real numbers (Fig. 4). As shown in the figure, from the
$1\times 96$ vectors in the columns of Fig. 3, each nonzero component of that
vector $\mathbf{v}_{i}$ was expanded into an 8-dimensional vector with the
following process. First the nearest integer $k$ between $1$ and $7$ inclusive
to $\mathbf{v}_{i}$ was found. Then, the matrix values were set as

$A_{mi}=\delta_{mk}\cdot\mathbf{v}_{i}/k,$ (1)

where $\delta$ here is the Kronecker delta, we index from $0$ (top left is
$A_{00}$), and $m$ ranges from $0$ to $7$. For zero components
($\mathbf{v}_{i}=0$), a $1$ was simply placed at $A_{0i}$. We point out that
we tested other encodings, such as the one from Ref. [7], but these were
susceptible to mode collapse. The encoding proposed in this work successfully
encodes decimal values both through the actual matrix entry and its location.

###  II.3 Model

A GAN is a type of generative model that has two competing neural networks, a
generator and a discriminator, as shown in Fig. 1. Traditional GANs, however,
can suffer from issues such as mode collapse and gradient vanishing, so the
Wasserstein GAN with Gradient Penalty [16] was used instead. In the
Wasserstein GAN with gradient penalty, the loss functions are

$\displaystyle\mathrm{Loss}_{D}$
$\displaystyle=\underset{\bm{\tilde{x}}\sim\mathbb{P}_{g}}{\mathbb{E}}[D(\bm{\tilde{x}})]-\underset{\bm{x}\sim\mathbb{P}_{r}}{\mathbb{E}}[D(\bm{x})]$
(2)
$\displaystyle\quad+\lambda\underset{\bm{\hat{x}}\sim\mathbb{P}_{\bm{\hat{x}}}}{\mathbb{E}}[(\|\nabla_{\bm{\hat{x}}}D(\hat{\bm{x}})\|-1)^{2}],$
$\displaystyle\mathrm{Loss}_{G}$
$\displaystyle=-\underset{\bm{\tilde{x}}\sim\mathbb{P}_{g}}{\mathbb{E}}[D(\bm{\tilde{x}})].$
(3)

Here $D(x)$ represents the output of the discriminator and $\mathbb{E}$ is the
expectation value (average). Then the parameters are updated with an
optimizer,

$\displaystyle w$ $\displaystyle\leftarrow
w+\alpha\cdot\operatorname{Optimizer}(w,\nabla_{w}\mathrm{Loss}_{D}),$ (4)
$\displaystyle\theta$
$\displaystyle\leftarrow\theta-\alpha\cdot\operatorname{Optimizer}(\theta,\nabla_{\theta}\mathrm{Loss}_{G}),$
(5)

where $w$ are the discriminator’s parameters, $\theta$ is the generator’s
parameters, and $\alpha$ is the learning rate. RMSProp was chosen as the
optimizer [17], after testing out several other options, such as Adam [18].

####  II.3.1 Training

The model was first trained on the OQMD dataset for 400 epochs. It was then
transfer learned onto the SuperCon dataset or a subset of it, on which it
would train for another 500 epochs. Transfer learning onto four different
datasets (cuprates, pnictides, others, and everything together) resulted in
four different versions of the GAN. Afterwards, the testing procedure and the
data analysis were conducted, and then the hyperparameters were updated based
on the results. The training curves for the final model on each of these sets
are displayed in Fig. 5. Notably, they were all able to converge and stabilize
over the 500 epochs.

Figure 5: The generator loss against training epoch for each of the four
datasets the GAN trained on: (a) All of SuperCon; (b) Others (i.e. not
cuprates or pnictides); (c) Cuprates; and (d) Pnictides.

####  II.3.2 Testing

After each training process, $5,000$ hypothetical compounds were generated
from the versions of the GAN trained on the smaller sets and $30,000$ from the
version trained on everything. The generated predictions were then inspected
with various quality checks. Each compound was first tested for validity using
the charge neutrality and electronegativity check features of the SMACT
package [19]. The package tests the compound for electronegativity and charge
balance to determine whether it is a valid compound or not. Each prediction
was then checked for uniqueness—whether the compound showed up earlier in the
generated list—and novelty—whether the compound was in the training dataset.
These three checks looked at the general quality of the model. Then, more
specific to superconductivity, each compound was run through the model from
Ref. [6] to check whether it is a superconductor or not, as well as to predict
its critical temperature. Of course, to be sure of superconductivity, the
compounds must be synthesized and tested. Lastly, the formation energy of each
generated compound was calculated using the model from Ref. [20], which
indicates the stability of the compound. These tests will be discussed further
in Section III, along with the actual results.

###  II.4 Clustering

In order to further assess the quality of predictions, clustering analysis was
performed. Clustering is an unsupervised machine learning technique whose main
goal is to unveil hidden patterns in the data. It was recently applied to
superconductors from SuperCon database [15]. Depending on the data set,
different clustering algorithms were used, such as k-means, hierarchical,
Gaussian mixtures, self-organizing maps, etc. The results showed that in case
of superconductors, clustering methods can achieve, and in some cases exceed,
human level performance.

In order to visualize clustering results, different techniques can be used. It
was shown that for superconductors the so-called t-SNE produces the best
results [15]. t-SNE is a non-linear dimensionality reduction technique which
allows higher dimensional data to be represented in 2D or 3D [21]. In case of
superconductors, the data points are $96$-dimensional (each compound is
represented by a 1 $\times$ 96 matrix), as discussed in Section II.2. t-SNE
reduces these dimensions down to either two or three, which allows easy
visualization. We point out, however, that these reduced dimensions do not
have any physical meaning.

##  III Results

After training, from the four different versions of the GAN, we generated four
superconductor candidate lists with either 5,000 or 30,000 chemical
compositions. We then ran the predicted superconductors through a series of
tests to evaluate their quality. The first few were general tests, and the
rest were in the context of superconductivity using existing computational
models, as experimentally testing all of them is unfeasible.

###  III.1 Duplicates and Novelty

We first screened the output for duplicates within the generated sets
themselves and then for duplicates between the generated set and the dataset
of known superconductors that it trained on. The results are tabulated in
Table 2. As seen in the table, the number of repeats within the generated
samples were relatively low (high uniqueness), with the exception of
pnictides, which had more duplicates than the rest. This is likely due to the
fact that the dataset of pnictides was significantly smaller than the rest
(see Table 1). We speculate that this overall low rate of duplicates stems
from the fact that the model is able to handle decimals (see Section IIB and
Fig.4), which opens up a large composition space for it to explore.

The percent of predicted compounds that were novel, i.e. not already known to
exist, listed in Table 2 is also very high across all four GANs. These two
results demonstrate that all versions of ScGAN can generate a diverse array of
new compositions.

GAN Version | Novel % | Unique %
---|---|---
Entirety of SuperCon | 99.69% | 96.78%
Cuprate | 99.74% | 92.98%
Pnictides | 99.32% | 58.74%
Others | 98.89% | 91.58 %
Table 2: The percentage of generated predictions that were novel (not in the
training set) and unique (distinct from others in the given generated set) for
each of the versions of ScGAN trained on the given training datasets on the
left.

###  III.2 Formation Energy

As mentioned in the previous section, we found the formation energies of the
predicted compounds from the GANs using ElemNet [20]. It was indicated in Ref.
[20] that a negative value of formation energy is a good indicator of
stability, i.e. the possibility of being synthesized in the lab. In Fig. 6 we
display the values of formation energy for all predictions. We see from the
distributions of formation energies that most of predicted compounds have
calculated formation energies less than zero. Even though this does not
provide definitive proof of stability, it is a general indication that most of
the predicted compounds are stable and can be synthesized in the lab.

Figure 6: Distribution of the formation energies of the predicted compounds
from the four versions of the GAN: (a) everything, (b) others, (c) cuprates,
(d) pnictides.

###  III.3 Superconductivity

As a next test, we ran the predicted compounds through the $K$-Nearest
Neighbors (KNN) classification model from Ref. [6] for predicting
superconductivity based on elemental composition, in order to to check if the
predictions of two machine learning models (GAN and KNN) would agree. However,
the probabilistic nature of the machine learning model had to be taken into
account. If $p_{sc}$ is the proportion of predictions that came up
superconducting according to the model, then we can estimate $\rho_{sc}$, the
true proportion of superconducting entries, using Bayesian statistics.
Denoting tp and fp as the true positive and false positive rates of the
classification model, respectively, we can write

$\rho_{sc}\cdot\textit{{tp}}+(1-\rho)\cdot\textit{{fp}}\approx p_{sc}.$ (6)

Solving for $\rho_{sc}$ gives

$\rho_{sc}\approx\frac{p_{sc}-\textit{{fp}}}{\textit{{tp}}-\textit{{fp}}}.$
(7)

The true positive and false positive rates here are reported from Ref. [6]:
$\textit{{tp}}=98.69\%$ and $\textit{{fp}}=16.94\%$. The output percentages
along with the estimates of the true proportions calculated from equation 7
are tabulated in Table 3.

GAN Version | Output % | True % Estimate
---|---|---
Entirety of SuperCon | 74.50 % | 70.42%
Cuprates | 75.76% | 71.95%
Pnictides | 72.44% | 67.89%
Others | 69.58% | 64.39%
Table 3: The percentages of generated predictions that were determined by the
KNN model to be superconducting for different training sets along with the
estimated real percentage of the predictions that were superconducting. The
true percentages were estimated according to Eq. 7.

All four GANs achieved very high percentages of superconducting material
according to the KNN model, especially when compared to the 3% figure from the
manual search in Ref. [3]. However, the only definite test of
superconductivity can come from experimental measurements.

###  III.4 Critical Temperature Estimates

Figure 7: Distributions of the critical temperatures of the predictions of the
four different versions of ScGAN: trained on (a) everything, (b) others, (c)
cuprates, (d) pnictides.

We also calculated the critical temperatures of our predictions using the
regression model from Ref. [6]. Similar to the superconductivity tests in the
previous subsection, these calculated values can only be taken as
approximations. However, the regression model can still provide us with a
general understanding of the capabilities of ScGAN. The critical temperature
outputs from the model are summarized in Table 4 and the distributions are in
shown in Fig. 7. While the distributions are somewhat broad, the GANs were
still able to find several superconductors with predicted critical
temperatures higher than $100\,\mathrm{K}$, which exceeds the manual search
maximum of $58\,\mathrm{K}$ (though that search was mostly restricted to
pnictides) and several of the previous machine learning approaches [6, 5].
This is not surprising, as those previous searches were limited to existing
databases of stoichiometric compounds, which meant that these forward design
approaches could only produce a limited number of candidates, mostly with low
critical temperatures.

Training Data | Average $T_{c}$ | Standard Dev. | Max $T_{c}$
---|---|---|---
Entirety of SuperCon | $6.53\,\mathrm{K}$ | $11.76\,\mathrm{K}$ | $123.25\,\mathrm{K}$
Cuprates | $59.34\,\mathrm{K}$ | $24.78\,\mathrm{K}$ | $133\,\mathrm{K}$
Pnictides | $20.41\,\mathrm{K}$ | $13.69\,\mathrm{K}$ | $51.98\,\mathrm{K}$
Others | $72.68\,\mathrm{K}$ | $21.24\,\mathrm{K}$ | $116.55\,\mathrm{K}$
Table 4: Summary statistics of the predicted critical temperatures of the
generated predictions that were determined by the regression model for
different training sets.

###  III.5 Ability to Learn Features

We then looked at the types of superconductors that were generated by the
different versions of ScGAN. The distributions are given in Table 5, and we
can see that each versions of ScGAN generated mostly superconductors that
matched their training data. This indicates the GAN was able to detect the
different underlying features behind these different major classes of
superconductors.

Training Data | Cuprate % | Pnictide % | Other %
---|---|---|---
Cuprate | 92.76% | 0.06% | 7.18%
Pnictides | 0.02 % | 99.84 % | 0.14%
Others | 0.14% | 0.6 % | 99.26%
Table 5: The distribution of the predicted superconductors across the
different classes of superconductors for the different versions of the GAN.

###  III.6 Clustering results

In Fig. 8 we display the results of clustering analysis on three sets of
predictions: panel (a) for cuprates, panel (b) for pnictides and panel (c) for
others. The results are visualized with the help of t-SNE. As pointed out
above, the two t-SNE dimensions, Y1 and Y2, do not have any physical meaning.
Full circles of different colors represent different families of
superconductors from SuperCon database, whereas purple open circles represent
GAN predictions. As can be seen from all three panels, GANs were able to
generate new superconductors from all known families of cuprates, pnictides
and other superconductors. However, GANs did not predict any new families of
superconductors.

Figure 8: Clustering of the predicted compounds from various versions of the
GAN: (a) cuprates, (b) pnictides and (c) others. Full circles represent the
data points from SuperCon and purple open circles are GAN predictions.

###  III.7 Promising Candidates

After running through the $T_{c}$ prediction model, we manually identified the
ones that looked the most promising, including some with very high critical
temperatures. It turns out that most of these were cuprates, which is not
surprising, considering that superconductors with highest critical
temperatures in the training set are cuprates. We then checked the
Crystallography Open Database (COD) [22] to see if those compounds were in the
database. The ones listed in Table 6 could not be found neither in COD nor in
SuperCon, showing that this model overcomes the limitations of the previous
forward design models and finds completely new superconductors. A more
comprehensive list of predictions is available upon reasonable request.

Compound | Predicted $T_{c}$ | Class
---|---|---
$\mathrm{PrCaBiSr_{2}Cu_{2}O_{7.46}}$ | $104.6\,\mathrm{K}$ | Cuprates
$\mathrm{YTiSr_{2}Cu_{2.74}O_{6.76}}$ | $91.7\,\mathrm{K}$ | Cuprates
$\mathrm{TeYSr_{2}Cu_{2}O_{7.75}}$ | $89.8\,\mathrm{K}$ | Cuprates
$\mathrm{TlCaSr_{2}Cu_{2}O_{7.82}}$ | $73.9\,\mathrm{K}$ | Cuprates
$\mathrm{YCaBa_{2}ZnCu_{2.36}O_{7.54}}$ | $71.5\,\mathrm{K}$ | Cuprates
$\mathrm{HgCsSrCa_{2}Cu_{2.56}O_{8.66}}$ | $69.8\,\mathrm{K}$ | Cuprates
$\mathrm{GdCaRuSr_{1.83}Cu_{2}O_{8.71}}$ | $40.8\,\mathrm{K}$ | Cuprates
$\mathrm{C_{2.52}Ni_{0.92}Y_{0.71}Th}$ | $85.3\,\mathrm{K}$ | Others
$\mathrm{Si_{0.62}V_{0.91}Zr_{0.83}}$ | $84.7\,\mathrm{K}$ | Others
$\mathrm{Al_{2.34}Te_{0.64}Ir_{1.07}}$ | $84.7\,\mathrm{K}$ | Others
$\mathrm{Be_{0.16}Si_{1.09}V_{2.67}Y_{1.72}}$ | $62.4\,\mathrm{K}$ | Others
$\mathrm{Cu_{1.13}Nb_{3.0}Sb_{0.72}Ir_{1.05}}$ | $59.4\,\mathrm{K}$ | Others
$\mathrm{Ga_{0.62}Nb_{2.88}Sn_{0.65}Te_{0.79}}$ | $40.8\,\mathrm{K}$ | Others
$\mathrm{B_{1.73}C_{1.03}Ni_{1.12}Y_{0.66}Pt_{0.64}}$ | $40.8\,\mathrm{K}$ | Others
$\mathrm{RuTeSeFe}$ | $35.6\,\mathrm{K}$ | Pnictides
$\mathrm{TeSSeFe_{1.05}}$ | $31.0\,\mathrm{K}$ | Pnictides
$\mathrm{CeCoAs_{2.15}Fe_{1.39}}$ | $23.3\,\mathrm{K}$ | Pnictides
$\mathrm{CeThPAsFe_{1.59}}$ | $12.2\,\mathrm{K}$ | Pnictides
$\mathrm{GaPrCa_{2.58}As_{12.44}Fe_{6.34}}$ | $11.9\,\mathrm{K}$ | Pnictides
$\mathrm{NdOAsFe}$ | $4.5\,\mathrm{K}$ | Pnictides
Table 6: Promising superconductor candidates generated by ScGANs, that do not
exist in current databases. Also shown are their predicted critical
temperatures [6].

##  IV Conclusion

For decades the search for new superconductors has relied on the serendipity
of material scientists to synthesize a new material with superconducting
proprieties. This paper introduced a novel method to search for
superconductors—discovering candidates with a generative adversarial network.
In contrast to previous computational methods which attempted to predict new
superconductors off of existing datasets, this model predicted compounds
directly. This “inverse design” approach proved to be far more powerful than
manual search methods and previous computational methods, with the model being
able to generate thousands of candidates with a wide-range of critical
temperatures that lied outside of existing databases (both superconductor and
general inorganic compound databases). Even while only training on chemical
compositions, more than $70\%$ of the GANs predictions were cross-checked with
a separate model to be potentially superconducting (the only way to know for
sure, however, would be to synthesize and check these compounds in a lab). Of
these, several were promising HTS candidates listed in Table 6. We point out
that previous models would have been unable to find such candidates as they
were outside of current databases.

While the compounds generated were new, our clustering showed that the GAN did
not generate any new families of superconductors. However, it was still able
to generate non-stoichiometric compounds, widening the scope of the
computational search. Future studies should look into some improvements that
can be made, such as being able to account for charge neutrality and crystal
structure in the compound encodings. Chemical checks on the predicted
compounds (using SMACT as detailed earlier) revealed that while there was
electronegativity balance, charge balance was not always exact for the
superconductor candidates. Furthermore, crystal structure has been known to
play a significant role in superconductivity, so including it in calculations
would most likely result in improvements. However, this may prove to be a
difficult endeavor as crystal structure is not well-documented in existing
databases. Active transfer learning could also be attempted to narrow down the
predictions to high temperature superconductors only [23]. However, even
without these possible improvements, the model at its current version is still
very promising and can be applied to search for superconductors, starting with
the candidates identified in this paper.

## References

  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, and J. M. Martinis, Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Hirsch  _et al._ [2015] J. Hirsch, M. Maple, and F. Marsiglio, Superconducting materials classes: Introduction and overview, Physica C: Superconductivity and its Applications 514, 1 (2015), superconducting Materials: Conventional, Unconventional and Undetermined.
  * Hosono  _et al._ [2015] H. Hosono, K. Tanabe, E. Takayama-Muromachi, H. Kageyama, S. Yamanaka, H. Kumakura, M. Nohara, H. Hiramatsu, and S. Fujitsu, Exploration of new superconductors and functional materials, and fabrication of superconducting tapes and wires of iron pnictides, Science and Technology of Advanced Materials 16, 033503 (2015).
  * Bedolla  _et al._ [2020] E. Bedolla, L. C. Padierna, and R. Castañeda-Priego, Machine learning for condensed matter physics, Journal of Physics: Condensed Matter 33, 053001 (2020).
  * Stanev  _et al._ [2018] V. Stanev, C. Oses, A. G. Kusne, E. Rodriguez, J. Paglione, S. Curtarolo, and I. Takeuchi, Machine learning modeling of superconducting critical temperature, npj Computational Materials 4, 29 (2018).
  * Roter and Dordevic [2020] B. Roter and S. Dordevic, Predicting new superconductors and their critical temperatures using machine learning, Physica C: Superconductivity and its Applications 575, 1353689 (2020).
  * Konno  _et al._ [2021] T. Konno, H. Kurokawa, F. Nabeshima, Y. Sakishita, R. Ogawa, I. Hosako, and A. Maeda, Deep learning model for finding new superconductors, Phys. Rev. B 103, 014509 (2021).
  * Goodfellow  _et al._ [2014] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, Generative adversarial nets, in  _Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 2_ , NIPS’14 (MIT Press, Cambridge, MA, USA, 2014) p. 2672–2680.
  * Dash  _et al._ [2021] A. Dash, J. Ye, and G. Wang, A review of generative adversarial networks (gans) and its applications in a wide variety of disciplines–from medical to remote sensing, arXiv preprint arXiv:2110.01442 (2021).
  * Dan  _et al._ [2020] Y. Dan, Y. Zhao, X. Li, S. Li, M. Hu, and J. Hu, Generative adversarial networks (gan) based efficient sampling of chemical composition space for inverse design of inorganic materials, npj Computational Materials 6, 84 (2020).
  * Hu  _et al._ [2020] T. Hu, H. Song, T. Jiang, and S. Li, Learning representations of inorganic materials from generative adversarial networks, Symmetry 12, 10.3390/sym12111889 (2020).
  * [12] National Institute for Materials Science, .
  * Saal  _et al._ [2013] J. E. Saal, S. Kirklin, M. Aykol, B. Meredig, and C. Wolverton, Materials design and discovery with high-throughput density functional theory: The open quantum materials database (oqmd), JOM 65, 1501 (2013).
  * Kirklin  _et al._ [2015] S. Kirklin, J. E. Saal, B. Meredig, A. Thompson, J. W. Doak, M. Aykol, S. Rühl, and C. Wolverton, The open quantum materials database (oqmd): assessing the accuracy of dft formation energies, npj Computational Materials 1, 15010 (2015).
  * Roter  _et al._ [2022] B. Roter, N. Ninkovic, and S. Dordevic, Clustering superconductors using unsupervised machine learning, Physica C: Superconductivity and its Applications 598, 1354078 (2022).
  * Gulrajani  _et al._ [2017] I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, and A. C. Courville, Improved training of wasserstein gans, Advances in neural information processing systems 30 (2017).
  * Hinton [2012] G. Hinton, Neural networks for machine learning: Lecture 6a overview of min-batch gradient descent (2012).
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * Davies  _et al._ [2016] D. Davies, K. Butler, A. Jackson, A. Morris, J. Frost, J. Skelton, and A. Walsh, Computational screening of all stoichiometric inorganic materials, Chem 1, 617 (2016).
  * Jha  _et al._ [2018] D. Jha, L. Ward, A. Paul, W.-k. Liao, A. Choudhary, C. Wolverton, and A. Agrawal, Elemnet: Deep learning the chemistry of materials from only elemental composition, Scientific Reports 8, 17593 (2018).
  * van der Maaten and Hinton [2008] L. van der Maaten and G. Hinton, Visualizing data using t-SNE, Journal of Machine Learning Research 9, 2579 (2008).
  * Vaitkus  _et al._ [2021] A. Vaitkus, A. Merkys, and S. Gražulis, Validation of the Crystallography Open Database using the Crystallographic Information Framework, Journal of Applied Crystallography 54, 661 (2021).
  * Kim  _et al._ [2021] Y. Kim, Y. Kim, C. Yang, K. Park, G. X. Gu, and S. Ryu, Deep learning framework for material design space exploration using active transfer learning and data augmentation, npj Computational Materials 7, 140 (2021).


# Screened plasmons of graphene near a perfect electric conductor

Afshin Moradi1<EMAIL_ADDRESS>Nurhan Türker Tokan2<EMAIL_ADDRESS>_1 Department of Engineering Physics, Kermanshah University of Technology,
Kermanshah, Iran
2Department of Electronic and Communications Engineering, Yildiz Technical
University, Istanbul, Turkey _

###### Abstract

Screened plasmon properties of graphene near a perfect electric conductor are
investigated using classical electrodynamics and a linearized hydrodynamic
model that includes Fermi correction. A general expression for the dispersion
relation of the mentioned screened plasmonic waves is given and illustrated
graphically. The result indicates that for realistic wavenumbers, the
dispersion relation of plasmonic waves of isolated graphene is almost
unaffected by the Fermi correction, while this correction is an important
factor for the screened plasmons of graphene near a perfect electric
conductor, where it increases the frequency of surface waves. The results show
that near the graphene neutrality point, the surface wave has a linear
dispersion with a universal speed close to $v_{\mathrm{F}}/\sqrt{2}$. Such
linear dispersion for surface waves (also known as energy waves) appears to be
a common occurrence when a splitting of plasma frequencies occurs, e.g. in the
electron-hole plasma of graphene [W. Zhao et al., Nature 614, 688 (2023)].
Furthermore, analytical expressions for the energy parameters (the power flow,
energy density, and energy velocity) of screened plasmons of the system are
derived. Also, the analytical expressions are derived and analyzed for the
damping function and surface plasmon and electromagnetic field strength
functions of surface waves of the system with small intrinsic damping.

##  I Introduction

Graphene is a two-dimensional (2D) material consisting of carbon atoms
arranged in a hexagonal lattice, which was discovered by Novoselov et al.
K.S.N666  in 2004. Graphene has electrons that behave like massless Dirac
fluid B.W318 ; E.H.H205418 , and therefore extraordinary properties can be
observed in this 2D material. For example, graphene has carriers (i.e.,
electrons and holes) with extremely high mobility. Also, graphene supports the
propagation of surface plasmon polariton (SPP) in the region from infrared to
THz frequencies S.A.M016803 ; X.L351 ; P.A.D.G . Furthermore, the low loss of
SPPs of graphene up to mid-infrared frequencies also makes it a promising
alternative for future applications X.L351 .

The most important advantage of graphene in the plasmonics world is the
tunability of surface plasmons because the density of carriers in graphene can
be easily adjusted with doping and an electric gate. The conductivity
characteristic of graphene M.M1052  and graphene’s plasmonic properties
A.N.G749  can be well explained by the hydrodynamic model derived by Müller et
al. M.M025301  in the long-wavelength limit, i.e., $k\ll k_{\mathrm{F}}$,
where $k_{\mathrm{F}}$ is the Fermi wavenumber in doped graphene and $k$ is
the wavenumber of the plasmonic wave. Chaves et al. A.J.C195438  investigated
the excitation of plasmonic waves of graphene in the presence of a fast-moving
charge using the hydrodynamic model in the electrostatic approximation.
Ferreira et al. B.A.F033817  performed the quantization of graphene plasmons
using the hydrodynamic model in the absence of losses for three graphene-based
structures, i.e., a monolayer graphene, a bilayer graphene, and a graphene
near a perfect electric conductor (PEC). The hydrodynamic Dirac fluid also
shows interesting collective excitations, such as hydrodynamic bipolar plasmon
polaritons that exhibit a coupled collective excitation of electromagnetic and
electron-hole oscillations with the opposite motion of electrons and holes,
and energy wave (also known as the demon mode), which is a quasi-acoustic mode
in which the motion of relativistic electrons and holes are in the same
direction D.S083715 ; A.L245153 ; Z.S3285 ; A.L053001 ; D.S121405 ; A.L115449
; I.T144307 ; B.N.N167979 ; E.I.K245434 ; J.D023036 ; D.F941 ; B.N115402 .
Importantly, more recently Zhao et al. W.Z688  observed both hydrodynamic
plasmons and the hydrodynamic energy waves of Dirac fluid based on new on-chip
THz spectroscopy techniques, where the report by Zhao et al. may reveal new
opportunities to study the collective hydrodynamic excitations in graphene-
based materials.

Figure 1: (a) Side view of the system under study. A monolayer graphene near a
PEC. Monolayer graphene and PEC are separated from each other by a dielectric
medium of thickness $d$ and dielectric constant $\varepsilon_{1}$, while the
region $z>d$ is a semi-infinite dielectric with the dielectric constant
$\varepsilon_{2}$. (b) Snapshot of the electric-field pattern of screened
$p$-polarized plasmon of graphene near a PEC in the $xz$ plane that shows the
tangential component of the electric field vanishes at $z=0$. The magnitude of
the electric field is oscillating in the $x$-direction but decreases away from
the boundary. For the attenuation of the electric field, it can be found that
the rate of hyperbolic sine decaying to zero is much faster than that of the
exponential function. (c) Snapshot of the Poynting-vector pattern close to the
boundary. The power flows in the $0<z<d$ region and in the $z>d$ region are in
the same direction, giving a net power flow to the right for a positive value
of wavenumber $k$. One can see that $x$-component of the Poynting vector
$S_{x}$ is maximum at the dielectric-PEC boundary and the Poynting vector is
everywhere orthogonal to the electric field, as expected on physical grounds.

One of the main advantages of the hydrodynamic model for the study of
plasmonic waves of graphene is the possibility to include nonlocal and quantum
effects in its plasmonic response without a high computational burden
I.S.E133104 . Note that the condition $kc/(k_{\mathrm{F}}v_{\mathrm{F}})\gg 1$
must be satisfied for nonlocality to play an important role in the optical
spectrum of graphene, where $v_{\mathrm{F}}$ is the Fermi speed in doped
graphene, and $c$ is the speed of light in free space. Using the hydrodynamic
model, we analyzed the characteristics of energy density and power flow of the
$p$-polarized plasmonic waves of monolayer A.M63 ; A.M043103 ; A.M072114 ,
bilayer graphene A.M135 , and graphene on a conducting substrate A.M353 .
However, to the best of our knowledge, no explicit calculation can be found
for the energy behaviors of screened plasmons of graphene near a PEC. We note
that the plasmonic properties of graphene near a PEC may present some new
behaviors that make it appropriate for applications in the mid-infrared to a
few THz range of frequencies. For instance, the results by Gu et al. X.G071103
show that graphene near a PEC leads to additional field localization near
graphene and may increase the amplification factor of plasmons, as shown by
Morozov et al. M.Y.M40 .

Therefore, in the present work, we wish to investigate the screened plasmon
properties of graphene near a PEC for future applications. In this way, by
using the linearized hydrodynamic model that includes Fermi correction
A.J.C195438  and classical electrodynamic formulations, we derive the general
expressions for the dispersion relation, power flow, energy density, energy
(group) velocity, damping function, and surface plasmon and electromagnetic
field strength functions of the $p$-polarized surface waves of the system
under investigation.

##  II Theory

The side view of the system under study, i.e., a monolayer graphene near a PEC
is shown in panel (a) of Fig. 1 in a Cartesian coordinate system with
coordinates $(x,y,z)$. Note that the metal, typically gold or silver, may be
modeled as a PEC for the frequency range of interest, i.e., from GHz to a few
THz A.M2023 ; M.M.B494 . The monolayer graphene and PEC are separated by a
dielectric of thickness $d$ and dielectric constant $\varepsilon_{1}$, whereas
the region $z>d$ is assumed to be a semi-infinite dielectric with the
dielectric constant $\varepsilon_{2}$. The electronic behavior of the graphene
layer is modeled as a 2D massless Dirac electron fluid and we assume that the
equilibrium doping density of electrons in graphene is $n$. Now, let us
consider the propagation of a plasmon polariton that is a $p$-polarized
surface electromagnetic wave having $E_{x}$, $H_{y}$, and $E_{z}$ components
along the $x$-direction at the boundary between two dielectric media. In this
way, the homogeneous 2D massless Dirac electron fluid will be perturbed and
can be regarded as a charged fluid with the first-order perturbed values of
the electron fluid density per unit area $n_{\mathrm{g}}(x,t)$ and electric
current density flowing on the graphene surface
$\textbf{J}(x,t)=j_{x}(x,t)\textbf{e}_{x}$, where $\textbf{e}_{x}$ being the
unit vector along the $x$-axis. We assume that all physical quantities vary as
$e^{i(kx-\omega t)}$, where $k$ is the wavenumber (propagation constant) of
the wave.

Based on the linear Drude model with Fermi correction in the limit of large
wavelengths, covering a range of frequencies from the mid-infrared to the THz,
the electronic excitations on a doped graphene surface can be described by the
following set of hydrodynamic equations A.J.C195438 ; A.S.P195437

Figure 2:  The dimensionless screened plasmon-polariton frequencies
$\Omega=\omega/\omega_{0}$ of graphene near a PEC, as a function of
$K=k/k_{0}$, for different values of the parameter $D=k_{0}d$, when
$\varepsilon_{1}=\varepsilon_{\mathrm{SiO_{2}}}=3.9$, $\varepsilon_{2}=1$, and
$n=n_{0}$ that means $k_{0}=k_{\mathrm{F}}$, and
$\omega_{0}=\omega_{\mathrm{F}}$. The different panels refer to (a)
$d\rightarrow\infty$, (b) $d=1/k_{\mathrm{F}}$, (c) $d=0.5/k_{\mathrm{F}}$,
and (d) $d=0.1/k_{\mathrm{F}}$. In each panel, the dashed blue line, and red
line correspond to $s=0$ and $s\neq 0$, respectively. The dashed black line
marks a linear dispersion with a velocity of
$V_{\mathrm{s}}=v_{\mathrm{s}}/v_{\mathrm{F}}=1/\sqrt{2}$.

Figure 3:  The dimensionless screened plasmon-polariton frequencies
$\Omega=\omega/\omega_{0}$ of graphene near a PEC, as a function of
$K=k/k_{0}$, for different values of the parameter $k_{\mathrm{F}}=\sqrt{\pi
n}$, when $\varepsilon_{1}=\varepsilon_{\mathrm{SiO_{2}}}=3.9$,
$\varepsilon_{2}=1$, and $d=0.8/k_{0}$. The different panels refer to (a)
$k_{\mathrm{F}}=k_{0}$, (b) $k_{\mathrm{F}}=0.7k_{0}$, (c)
$k_{\mathrm{F}}=0.4k_{0}$, and (d) $k_{\mathrm{F}}=0.1k_{0}$. In each panel,
the dashed blue line, and red line correspond to $s=0$ and $s\neq 0$,
respectively. The dashed black line marks a linear dispersion with a velocity
of $V_{\mathrm{s}}=v_{\mathrm{s}}/v_{\mathrm{F}}=1/\sqrt{2}$.

$-e\partial_{t}n_{\mathrm{g}}(x,t)+\partial_{x}j_{x}(x,t)=0\;,$ (1)
$\partial_{t}j_{x}(x,t)=\frac{D_{\mathrm{g}}}{\pi}E_{x}\big{|}_{z=d}+ev_{\mathrm{s}}^{2}\partial_{x}n_{\mathrm{g}}(x,t)\;,$
(2)

where $e$ is the electron charge,
$D_{\mathrm{g}}=(e^{2}/\hbar^{2})E_{\mathrm{F}}$ is the Drude weight of
graphene, $E_{\mathrm{F}}=\hbar\omega_{\mathrm{F}}$ (with
$\omega_{\mathrm{F}}=v_{\mathrm{F}}k_{\mathrm{F}}$, and $\hbar=h/2\pi$, where
$h$ is the Planck constant) is the Fermi energy, $k_{\mathrm{F}}=\sqrt{\pi n}$
is the Fermi wavenumber and $v_{\mathrm{F}}\approx c/300$ is the Fermi speed
in doped graphene, as mentioned before. In the right-hand side of Eq. (2), the
first term is the force on electrons due to the tangential component of the
electric field, evaluated at the graphene surface $z=d$, and the second term
shows the Fermi pressure in the 2D electron gas with
$v_{\mathrm{s}}=\upsilon_{\mathrm{F}}/\sqrt{2}$ that is the energy sound speed
of Dirac electrons in graphene W.Z688 . Note that Eqs. (1) and (2) provides an
adequate description of the low-energy intraband electronic transitions in
doped graphene in the optical limit, specifically, for the wavenumbers $k\ll
k_{\mathrm{F}}$.

Now, by eliminating the induced density $n$ from Eqs. (1) and (2), and
applying $j_{x}=\sigma_{\mathrm{g}}E_{x}$ (where $\sigma_{\mathrm{g}}$ is the
conductivity of graphene), we find G.A.M1150

$\sigma_{\mathrm{g}}(k,\omega)=\dfrac{i}{\pi}\dfrac{\omega
D_{\mathrm{g}}}{\omega^{2}-v_{\mathrm{s}}^{2}k^{2}}\;.$ (3)

In order to determine the plasmonic properties of the system, we look for an
evanescent $p$-polarized wave described by an electric field of the form A.M
$\textbf{E}(x,z)=\left[\textbf{e}_{x}E_{x}+\textbf{e}_{z}E_{z}\right]e^{ikx}$.
Note that we have considered $\partial\textbf{E}/\partial y=0$, since the
$x$\- and $y$-directions are equivalent. The associated magnetic field has the
form $\textbf{H}(y)=\textbf{e}_{y}H_{y}e^{ikx}$, where we have omitted writing
explicitly a factor $\exp(-i\omega t)$ describing the time-dependence of the
wave. Note that $H_{y}$ and $E_{z}$ can be determined if the non-zero
longitudinal component $E_{x}$ is known. The Helmholtz equation for the
x-component electric field $E_{x}$ of the $p$-polarized surface wave can be
given by

$\left[\dfrac{d^{2}}{dz^{2}}-\kappa_{\ell}^{2}\right]E_{x}(z)=0\;,$ (4)

where $\kappa_{\ell}=\left[k^{2}-\varepsilon_{\ell}k_{0}^{2}\right]^{1/2}$
with $\ell=1,2$, denotes the attenuation constant in the regions $0<z<d$ and
$z>d$, and $k_{0}=\omega/c$ is free space wavenumber. To solve Eq. (4), we
have to provide appropriate boundary conditions. With the electric
conductivity of graphene, these boundary conditions at the surface $z=d$ can
be written as

$E_{x}\big{|}_{z=d+}=E_{x}\big{|}_{z=d-}\;,$ (5)
$H_{y}(z)|_{z=d+}-H_{y}(z)|_{z=d-}=-\sigma_{\mathrm{g}}E_{x}\big{|}_{z=d}\;,$
(6)

which express the continuity and discontinuity of the tangential components of
the electric and magnetic fields, respectively, across the surface $z=d$.
Also, the boundary condition satisfied by $E_{x}(z)$ at the surface $z=0$ is

$E_{x}\big{|}_{z=0}=0\;,$ (7)

that implies the tangential component of the electric field should vanish at
$z=0$ as can be seen in panel (b) of Fig. 1. With the above equations and also
using $E_{z}=-(ik/\kappa_{\ell}^{2})\partial E_{x}/\partial z$, and
$H_{y}=(i\omega\varepsilon_{0}\varepsilon_{\ell}/\kappa_{\ell}^{2})\partial
E_{x}/\partial z$, we will investigate the screened plasmon-polariton
properties of graphene near a PEC in the following sections.

###  II.1 Dispersion relation

We note that the surface field should decay for $z\rightarrow\infty$. Also,
the presence of a PEC at $z=0$ implies that $E_{x}(z=0)=0$. Therefore, the
appropriate expressions for $E_{x}(z)$, are as follows:

$E_{x}(z)=\left\\{\begin{array}[]{clcr}A_{-}\sinh\kappa_{1}z\;,&\mbox{$0\leq
z\leq d$\;.}\\\ A_{+}e^{-\kappa_{2}z}\;,&\mbox{$d\leq z$\;,}\\\
\end{array}\right.$ (8)

where the relations between the coefficients $A_{+}$ and $A_{-}$ can be
determined from the matching boundary conditions at the graphene surface,
i.e., $z=d$. Use of Eqs. (4) and (8) in the boundary conditions (5) and (6)
yields the condition that

$\dfrac{\varepsilon_{1}}{\kappa_{1}}\coth\kappa_{1}d+\dfrac{\varepsilon_{2}}{\kappa_{2}}+\dfrac{i\sigma_{\mathrm{g}}}{\omega\varepsilon_{0}}=0\;.$
(9)

The roots of this transcendental equation, which can only be solved
numerically, provide the dispersion relation of the screened plasmon-
polaritons of graphene near a PEC. This dispersion relation has two tuning
parameters: the graphene-PEC distance $d$, and the graphene sheet carrier
density $n$, which controls the conductivity of graphene. When $d$ is very
large such that $\coth\kappa_{1}d\approx 1$, Eq. (9) reduces to the dispersion
relation for the surface plasmon-polaritons supported by isolated monolayer
graphene P.A.D.G ; A.M .

Let us note that the dispersion relation for the screened plasmon-polaritons
of graphene near a PEC coincides with the quasi-acoustic plasmon-polaritons in
symmetric bilayer graphene, provided $d=d_{\mathrm{bilayer}}/2$, where
$d_{\mathrm{bilayer}}$ is the interlayer distance in the bilayer graphene.
This fact can be understood in terms of image charges. That is why the
screened plasmon-polariton introduced here is also called quasi-acoustic
plasmon-polariton. Note that such a quasi-acoustic mode seems to be a common
occurrence when a splitting of plasma frequencies happens due to the
electrostatic interaction, e.g., in the electron-hole plasma of graphene
D.S083715 ; W.Z688 , or in the coupling between the interlayer in the bilayer
graphene A.M135 .

If we neglect the retardation effects, i.e., $k\gg k_{0}$, from Eq. (9) we
find the dispersion relation of the screened plasmons of graphene near a PEC,
as

$\omega=\left[v_{\mathrm{s}}^{2}k^{2}+\dfrac{D_{\mathrm{g}}}{\pi\varepsilon_{0}}\dfrac{k}{\varepsilon_{1}\coth
kd+\varepsilon_{2}}\right]^{1/2}\;.$ (10)

From Eq. (10), we can distinguish two different dimensionality regimes
depending on two cases of $kd\gg 1$ and $kd\ll 1$. For $kd\gg 1$, where
graphene and PEC decouple, we may use the asymptotic expression $\coth
kd\approx 1$. Thus, the dispersion relation can be written as

$\omega=\left[v_{\mathrm{s}}^{2}k^{2}+\dfrac{D_{\mathrm{g}}}{\pi\varepsilon_{0}\left(\varepsilon_{1}+\varepsilon_{2}\right)}k\right]^{1/2}\approx\sqrt{\dfrac{D_{\mathrm{g}}}{\pi\varepsilon_{0}\left(\varepsilon_{1}+\varepsilon_{2}\right)}k}\;.$
(11)

which is exactly the same as the well-known dispersion relation of the surface
plasmons of isolated graphene in free space, when
$\varepsilon_{1}=1=\varepsilon_{2}$ A.J.C195438 , with the approximate result
valid for realistic ($k\ll k_{\mathrm{F}}$ ) wavenumbers. This result means
that the dispersion relation of an isolated graphene is almost unaffected by
the Fermi correction, i.e., the internal pressure force of the electron [the
term with $v_{\mathrm{s}}^{2}$ in Eq. (3)].

On the other hand, for $kd\ll 1$ we may use the asymptotic expression $\coth
kd\approx 1/kd$. Thus, the dispersion relation can be written as

$\omega=\left[v_{\mathrm{s}}^{2}+\dfrac{D_{\mathrm{g}}}{\pi\varepsilon_{0}}\dfrac{d}{\varepsilon_{1}+\varepsilon_{2}kd}\right]^{1/2}k\approx\left[v_{\mathrm{s}}^{2}+\dfrac{D_{\mathrm{g}}d}{\pi\varepsilon_{0}\varepsilon_{1}}\right]^{1/2}k\;.$
(12)

As a new interesting result, it is clear that in this case, the internal
pressure force of the electron is an important term and increases the
frequency of surface plasmons. In fact, for graphene near a PEC, the
dispersion is strongly dependent on $d$. We note that Chaves et al.
A.J.C195438  showed that for $d$ about $1.5$nm, the screened plasmons of
graphene near a PEC can appear in the mid-infrared with a wavenumber of the
order of $200\mu$m-1 (corresponding to a $\lambda_{\mathrm{spp}}=2\pi/k\approx
30$nm). For a Fermi energy of graphene about $E_{\mathrm{F}}=0.4$eV, we find
$kc/(k_{\mathrm{F}}v_{\mathrm{F}})\sim 100$, which places graphene in the
strong nonlocal regime.

Also, the phase and group velocities of screened plasmons of the system can be
obtained from Eq. (10). For the phase velocity, we have

$v_{\mathrm{phase}}=\left[v_{\mathrm{s}}^{2}+\dfrac{D_{\mathrm{g}}}{\pi\varepsilon_{0}k}\dfrac{1}{\varepsilon_{1}\coth
kd+\varepsilon_{2}}\right]^{1/2}\;,$ (13)

while for the group velocity by derivation of Eq. (10) with respect to
$\omega$, we find

$v_{\mathrm{group}}=\dfrac{v_{\mathrm{s}}^{2}k+\dfrac{D_{\mathrm{g}}}{2\pi\varepsilon_{0}}\dfrac{\varepsilon_{1}\coth
kd+\varepsilon_{2}+\varepsilon_{1}\dfrac{kd}{\sinh^{2}kd}}{\left[\varepsilon_{1}\coth
kd+\varepsilon_{2}\right]^{2}}}{\left[v_{\mathrm{s}}^{2}k^{2}+\dfrac{D_{\mathrm{g}}}{\pi\varepsilon_{0}}\dfrac{k}{\varepsilon_{1}\coth
kd+\varepsilon_{2}}\right]^{1/2}}\;.$ (14)

To see clearly the character of the dispersion relation for the screened
plasmon-polariton of graphene near a PEC, first let us introduce the
dimensionless variables $K=k/k_{0},\Omega=\omega/\omega_{0}$,
$V_{\mathrm{s}}=v_{\mathrm{s}}/v_{\mathrm{F}}=1/\sqrt{2}$,
$V_{\mathrm{c}}=c/v_{\mathrm{F}}\approx 300$, $D=k_{0}d$, where
$k_{0}=\sqrt{\pi n_{0}}$, and $\omega_{0}=v_{\mathrm{F}}k_{0}$. Note that for
$n_{0}=n$, we have $k_{0}=k_{\mathrm{F}}$ and
$\omega_{0}=\omega_{\mathrm{F}}$. Now, in Fig. 2, we show the dependence of
the dimensionless frequency $\Omega=\omega/\omega_{0}$ on the dimensionless
variable $K=k/k_{0}$, for different values of the parameter $D=k_{0}d$, when
$\varepsilon_{1}=\varepsilon_{\mathrm{SiO_{2}}}=3.9$, $\varepsilon_{2}=1$, and
$n=n_{0}$ that means $k_{0}=k_{\mathrm{F}}$, and
$\omega_{0}=\omega_{\mathrm{F}}$. One can see that the behavior of the
screened plasmon-polariton depends on the value of $d$, where the decreasing
thickness of the spacer layer $d$, red-shifts the frequency of the surface
wave. More importantly, it can be seen that with the decreasing thickness of
the spacer layer $d$, internal interaction force plays an important role in
the dispersion relation of the surface wave for $k\ll k_{\mathrm{F}}$. We
observe that in the presence of the nonlocal effects, screened plasmon-
polariton of the system has a phase velocity that can be made arbitrarily
close to the energy sound speed of Dirac electrons in graphene by tuning the
spacer layer $d$. Furthermore, from Fig. 2, we observe for the constant values
of frequency and graphene sheet carrier density, the wavelength of screened
plasmon-polariton considerably decreases when the PEC draws near to graphene,
as can be easily concluded from Eq. (12).

The effect of the graphene sheet carrier density $n$ on the dispersion
relation for the screened plasmon-polariton of graphene near a PEC is shown in
Fig. 3. It can be seen that as the carrier density decreases to lower values,
the surface wave frequency decreases. It is clear that at extremely low charge
density, the surface wave has a velocity close to $v_{\mathrm{F}}/\sqrt{2}$.
This means that near the graphene neutrality point, the surface wave has a
linear dispersion with a universal speed close to $v_{\mathrm{F}}/\sqrt{2}$
W.Z688 . One may conclude from panel (d) of Fig. 3 that near the neutral point
of graphene, the local model shows an incorrect result. Also, from Fig. 3, it
is clear that for the constant values of $\omega$ and $d$, the wavelength of
screened plasmon-polariton considerably decreases with a decrease in the
graphene sheet carrier density, as can be easily seen from Eq. (12).

###  II.2 Power flow

For the power flow density associated with a surface wave of graphene near a
PEC, we have, in the three media,

$\textbf{S}=\left\\{\begin{array}[]{clcr}\textbf{E}_{-}\times\textbf{H}_{-}\;,&\mbox{$0<z<d$\;,}\\\
\textbf{E}_{+}\times\textbf{H}_{+}\;,&\mbox{$d<z$\;,}\\\ \end{array}\right.$
(15)

where subscripts $-$ and $+$ denote the regions below and above the graphene
layer, and for $z=0$ using Eq. (A-4) in Appendix, we have
$S_{\mathrm{g}x}=-\frac{\pi
e}{D_{\mathrm{g}}}v_{\mathrm{s}}^{2}n_{\mathrm{g}}j_{x}$. After the
elimination of $j_{x}=\sigma_{\mathrm{g}}E_{x}$ and
$n_{\mathrm{g}}=-(k/e\omega)j_{x}$, and also using
$A_{+}=A_{-}\sinh(\kappa_{1}d)\exp(\kappa_{2}d)$, the cycle-averaged
$x$-components of S in Eq. (15), and also on the surface of graphene can be
written as

$S_{x}=\dfrac{\varepsilon_{0}k\omega}{2}A_{-}^{2}\left\\{\begin{array}[]{clcr}\dfrac{\varepsilon_{1}}{\kappa_{1}^{2}}\cosh^{2}\kappa_{1}z\;,&\mbox{$0<z<d$\;,}\\\
\dfrac{\varepsilon_{2}}{\kappa_{2}^{2}}\sinh^{2}\kappa_{1}d\;e^{-2\kappa_{2}(z-d)}\;,&\mbox{$d<z$\;,}\\\
\end{array}\right.$ (16)
$S_{\mathrm{g}x}=\dfrac{\varepsilon_{0}k\omega}{2}A_{-}^{2}\frac{\pi}{\varepsilon_{0}D_{\mathrm{g}}\omega^{2}}v_{\mathrm{s}}^{2}|\sigma_{\mathrm{g}}|^{2}\sinh^{2}\kappa_{1}d\;.\
$ (17)

The total power flow density associated with the $p$-polarized surface wave is
determined by integration over $z$. The power flow through an area in the $yz$
plane of infinite length in the $z$-direction and unit width in the
$y$-direction is

$\left\langle
S_{x}\right\rangle=\dfrac{\varepsilon_{0}k\omega}{4}\sinh^{2}\kappa_{1}d\;A_{-}^{2}\left[\dfrac{\varepsilon_{1}}{\kappa_{1}^{2}}\dfrac{d}{\sinh^{2}\kappa_{1}d}\right.\\\
\qquad\left.+\dfrac{\varepsilon_{1}}{\kappa_{1}^{3}}\coth\kappa_{1}d+\dfrac{\varepsilon_{2}}{\kappa_{2}^{3}}+\dfrac{2\pi}{\varepsilon_{0}D_{\mathrm{g}}\omega^{2}}v_{\mathrm{s}}^{2}|\sigma_{\mathrm{g}}|^{2}\right]\;,$
(18)

where $\left\langle\cdots\right\rangle\equiv\int_{0}^{+\infty}\cdots dz$. This
total power flow density (per unit width) is positive, as can be seen in panel
(c) of Fig. 1.

###  II.3 Energy distribution

For the cycle-averaged energy density distribution associated with a surface
wave of graphene near a PEC, we have,

$U=\dfrac{1}{4}\left\\{\begin{array}[]{clcr}\varepsilon_{0}\varepsilon_{1}|\textbf{E}_{-}|^{2}+\mu_{0}|\textbf{H}_{-}|^{2}\;,&\mbox{$0<z<d$\;,}\\\
\varepsilon_{0}\varepsilon_{2}|\textbf{E}_{+}|^{2}+\mu_{0}|\textbf{H}_{+}|^{2}\;,&\mbox{$d<z$\;,}\\\
\end{array}\right.$ (19)

where for $z=0$ using Eq. (A-3) in Appendix, we have

$U_{\mathrm{g}}=\frac{1}{4}\frac{\pi}{D_{\mathrm{g}}}\left(|j_{x}|^{2}+e^{2}v_{\mathrm{s}}^{2}|n_{\mathrm{g}}|^{2}\right)\;.$
(20)

Then Eqs. (19) and (20) yield,

$U=\dfrac{\varepsilon_{0}}{4}A_{-}^{2}\left\\{\begin{array}[]{clcr}\varepsilon_{1}\left[\sinh^{2}\kappa_{1}z+\dfrac{k^{2}+\varepsilon_{1}\frac{\omega^{2}}{c^{2}}}{\kappa_{1}^{2}}\cosh^{2}\kappa_{1}z\right]\;,&\mbox{$z<d$\;,}\\\
2\varepsilon_{2}\dfrac{k^{2}}{\kappa_{2}^{2}}\sinh^{2}\kappa_{1}d\;e^{-2\kappa_{2}(z-d)}\;,&\mbox{$z>d$\;,}\\\
\end{array}\right.$ (21)
$U_{\mathrm{g}}=\dfrac{1}{4}A_{-}^{2}\frac{\pi}{D_{\mathrm{g}}}\frac{\omega^{2}+v_{\mathrm{s}}^{2}k^{2}}{\omega^{2}}|\sigma_{\mathrm{g}}|^{2}\sinh^{2}\kappa_{1}d\;,$
(22)

where all contributions to the energy density are positive. The total energy
density associated with the screened plasmon-polaritons of graphene near a PEC
is again determined by integration over $z$, the energy per unit surface area
is

Figure 4:  The wavenumber dependence of dimensionless relaxation rate of the
long-wavelength screened plasmon-polaritons of graphene near a PEC as a
function of $K=k/k_{0}$, when
$\varepsilon_{1}=\varepsilon_{\mathrm{SiO_{2}}}=3.9$, $\varepsilon_{2}=1$. (a)
For different values of the parameter $D=k_{0}d$, when $n=n_{0}$. (b) For
different values of the parameter $k_{\mathrm{F}}=\sqrt{\pi n}$, when $D=5$.

Figure 5: The wavenumber dependence of the strength functions
$\Theta_{\mathrm{sp}}$ and $\Theta_{\mathrm{ph}}$ of the screened plasmon-
polaritons of graphene near a PEC, as a function of $K=k/k_{0}$, when
$\varepsilon_{1}=\varepsilon_{\mathrm{SiO_{2}}}=3.9$, $\varepsilon_{2}=1$ for
different values of the parameter $D=k_{0}d$, when $n=n_{0}$.

$\left\langle
U\right\rangle=\dfrac{\varepsilon_{0}}{4}A_{-}^{2}\sinh^{2}\kappa_{1}d\left[\frac{\varepsilon_{1}^{2}}{\kappa_{1}^{2}}\dfrac{\omega^{2}}{c^{2}}\dfrac{d}{\sinh^{2}\kappa_{1}d}+k^{2}\dfrac{\varepsilon_{1}}{\kappa_{1}^{3}}\coth\kappa_{1}d\right.\\\
\qquad\left.+k^{2}\dfrac{\varepsilon_{2}}{\kappa_{2}^{3}}+\frac{\pi}{\varepsilon_{0}D_{\mathrm{g}}}\frac{\omega^{2}+v_{\mathrm{s}}^{2}k^{2}}{\omega^{2}}|\sigma_{\mathrm{g}}|^{2}\right]\;.$
(23)

Let us note that in a recent work, Morozov and Popov M.Y.M22209  prepared a
concept of a terahertz waveguide plasmon amplifier based on a metal groove
with active graphene. In this way, they used concepts of plasmon energy to
gain insight into the physical origins of terahertz waveguide plasmon
amplification. However, in their results, we cannot see the contribution of
the energy density on the graphene surface, as shown by Eq. (22). Fortunately,
there is an easy way to check such results. In fact, we should find the group
and energy velocities of the wave under consideration, when losses are
neglected. If the group velocity is the same as the energy velocity, then the
obtained energy formulas are correct, as we are going to check this important
issue for our results in the following subsection.

###  II.4 Energy velocity

The energy velocity of a surface wave of graphene near a PEC is given as the
ratio of the total power flow density (per unit width) and total energy
density (per unit area), such as

$v_{\mathrm{energy}}=\omega k\\\
\dfrac{\frac{\varepsilon_{1}}{\kappa_{1}^{2}}\frac{d}{\sinh^{2}\kappa_{1}d}+\frac{\varepsilon_{1}}{\kappa_{1}^{3}}\coth\kappa_{1}d+\frac{\varepsilon_{2}}{\kappa_{2}^{3}}+\frac{2\pi}{\varepsilon_{0}D_{\mathrm{g}}\omega^{2}}v_{\mathrm{s}}^{2}|\sigma_{\mathrm{g}}|^{2}}{\frac{\varepsilon_{1}^{2}}{\kappa_{1}^{2}}\frac{\omega^{2}}{c^{2}}\frac{d}{\sinh^{2}\kappa_{1}d}+k^{2}\frac{\varepsilon_{1}}{\kappa_{1}^{3}}\coth\kappa_{1}d+k^{2}\frac{\varepsilon_{2}}{\kappa_{2}^{3}}+\frac{\pi}{\varepsilon_{0}D_{\mathrm{g}}}\frac{\omega^{2}+v_{\mathrm{s}}^{2}k^{2}}{\omega^{2}}|\sigma_{\mathrm{g}}|^{2}}\;.$
(24)

If we neglect the retardation effects, i.e., $k\gg k_{0}$, from Eq. (25) we
find

$v_{\mathrm{energy}}=\dfrac{v_{\mathrm{s}}^{2}k+\dfrac{D_{\mathrm{g}}}{2\pi\varepsilon_{0}}\dfrac{\varepsilon_{1}\coth
kd+\varepsilon_{2}+\varepsilon_{1}\dfrac{kd}{\sinh^{2}kd}}{\left[\varepsilon_{1}\coth
kd+\varepsilon_{2}\right]^{2}}}{\left[v_{\mathrm{s}}^{2}k^{2}+\dfrac{D_{\mathrm{g}}}{\pi\varepsilon_{0}}\dfrac{k}{\varepsilon_{1}\coth
kd+\varepsilon_{2}}\right]^{1/2}}\;,$ (25)

which is identical to the group velocity obtained in Eq. (14). This equality
confirms the correctness of the presented results. In fact, in general, the
group velocity is equal to the energy velocity in the absence of damping
A.M18373 ; A.M143901 ; A.M10760 .

###  II.5 Damping property

Now, we study the damping function of surface waves of the system. To obtain
an analytical expression for the damping function of the surface waves of
graphene near a PEC we use the perturbative method proposed by Loudon R.L233
and Nkoma et al. J.N3547 . Such a procedure enables us to calculate the true
surface wave damping rate to the first order in the damping parameter
$\gamma$, introduced to describe the intrinsic damping of crystal
oscillations. Also, this theory enables us to discuss both the propagation
length and the lifetime of a surface wave. The advantage of the perturbative
method is that the damping properties result from the calculation of real
dispersion relations. The plasmonic damping parameter or relaxation rate
$\Gamma(k,\omega)$ of the present case may be determined by the following
procedure. The kinetic and total energy densities (per unit area)
$U_{\mathrm{gk}}$, and $\left\langle U\right\rangle$ are first calculated in
the absence of damping. If a small amount of damping is now reintroduced, the
surface energy relaxation rate to the lowest order in $\gamma$ is
$\Gamma(k,\omega)=2\gamma U_{\mathrm{gk}}/\left\langle U\right\rangle$. where
from Eq. (A-3) we have
$U_{\mathrm{gk}}=\dfrac{\varepsilon_{0}}{4}A_{-}^{2}\frac{\pi}{\varepsilon_{0}D_{\mathrm{g}}}|\sigma_{\mathrm{g}}|^{2}\sinh^{2}\kappa_{1}d$.
Therefore, we get

$\Gamma(k,\omega)\\\
=\dfrac{2\gamma\frac{\pi}{\varepsilon_{0}D_{\mathrm{g}}}|\sigma_{\mathrm{g}}|^{2}}{\dfrac{d\frac{\varepsilon_{1}^{2}}{\kappa_{1}^{2}}\frac{\omega^{2}}{c^{2}}}{\sinh^{2}\kappa_{1}d}+k^{2}\frac{\varepsilon_{1}}{\kappa_{1}^{3}}\coth\kappa_{1}d+k^{2}\frac{\varepsilon_{2}}{\kappa_{2}^{3}}+\frac{\pi|\sigma_{\mathrm{g}}|^{2}}{\varepsilon_{0}D_{\mathrm{g}}}\frac{\omega^{2}+v_{\mathrm{s}}^{2}k^{2}}{\omega^{2}}}\;.$
(26)

Let us note that the frequency and wavenumber dependence of the damping
function comes from the retarded part of the plasmonic waves and it is easy to
find that, in the nonretarded limit, the total energy density (per unit area)
becomes twice as large as the kinetic energy density (per unit area) of the
system. As a consequence, the damping function of surface waves of the system
equals $\gamma$, i.e., it becomes a constant. Also, let us note that the
surface wave lifetime $T$ is simply the inverse of Eq. (26), i.e.,
$T(k,\omega)=\Gamma^{-1}(k,\omega)$, while their propagation length is given
by $L(k,\omega)=v_{\mathrm{group}}\Gamma^{-1}(k,\omega)$, where
$v_{\mathrm{group}}$ can be found from Eq. (24).

By using (26), the damping rate of the long-wavelength screened plasmon-
polaritons of graphene near a PEC, in terms of the dimensionless variables are
presented in Fig. 4 when $\varepsilon_{1}=\varepsilon_{\mathrm{SiO_{2}}}=3.9$,
$\varepsilon_{2}=1$. It is clear that the damping function of long-wavelength
screened plasmon-polaritons is approximately equal to $\gamma$ for a large
value of $k$. From panel (a) one can see that by decreasing values of $D$ for
a fixed value of $k_{\mathrm{F}}$, the screened plasmon-polaritons relaxation
rate increases sharply for a low value of $k$. On the other hand, from panel
(b) it is obvious that for a fixed value of $D$, by decreasing values of
$k_{\mathrm{F}}$, the screened plasmon-polaritons relaxation rate increases.

###  II.6 Surface plasmon and electromagnetic field strength functions

Since a screened plasmon-polariton is a coupled optical plasmon-photon wave,
we can introduce strength functions that characterize the quantitative
compositions of the mixed wave. We have $\left\langle
U_{\mathrm{sp}}\right\rangle/\left\langle U\right\rangle+\left\langle
U_{\mathrm{ph}}\right\rangle/\left\langle U\right\rangle=1$, where
$\left\langle U_{\mathrm{ph}}\right\rangle$ is the total photon energy density
and $\left\langle U_{\mathrm{sp}}\right\rangle$ is the total surface plasmon
energy density associated with the screened plasmon-polariton. Note that
$\left\langle U\right\rangle$ is the sum of the integrated energy densities,
i.e., Eq. (23), and also

$\left\langle
U_{\mathrm{ph}}\right\rangle=\dfrac{\varepsilon_{0}}{4}A_{-}^{2}\sinh^{2}\kappa_{1}d\left[\frac{\varepsilon_{1}^{2}}{\kappa_{1}^{2}}\dfrac{\omega^{2}}{c^{2}}\dfrac{d}{\sinh^{2}\kappa_{1}d}\right.\\\
\qquad\left.+k^{2}\dfrac{\varepsilon_{1}}{\kappa_{1}^{3}}\coth\kappa_{1}d+k^{2}\dfrac{\varepsilon_{2}}{\kappa_{2}^{3}}\right]\;,$
(27) $\left\langle
U_{\mathrm{sp}}\right\rangle=\dfrac{\varepsilon_{0}}{4}A_{-}^{2}\sinh^{2}\kappa_{1}d\frac{\pi}{\varepsilon_{0}D_{\mathrm{g}}}\frac{\omega^{2}+v_{\mathrm{s}}^{2}k^{2}}{\omega^{2}}|\sigma_{\mathrm{g}}|^{2}\;.$
(28)

Thus, for the surface plasmon strength function $\Theta_{\mathrm{sp}}$ and the
electromagnetic strength function $\Theta_{\mathrm{ph}}$ we obtain

$\Theta_{\mathrm{ph}}=\dfrac{\left\langle
U_{\mathrm{ph}}\right\rangle}{\left\langle U\right\rangle}\\\
=\dfrac{\frac{\varepsilon_{1}^{2}}{\kappa_{1}^{2}}\frac{\omega^{2}}{c^{2}}\frac{d}{\sinh^{2}\kappa_{1}d}+k^{2}\frac{\varepsilon_{1}}{\kappa_{1}^{3}}\coth\kappa_{1}d+k^{2}\frac{\varepsilon_{2}}{\kappa_{2}^{3}}}{\frac{d\frac{\varepsilon_{1}^{2}}{\kappa_{1}^{2}}\frac{\omega^{2}}{c^{2}}}{\sinh^{2}\kappa_{1}d}+k^{2}\frac{\varepsilon_{1}}{\kappa_{1}^{3}}\coth\kappa_{1}d+k^{2}\frac{\varepsilon_{2}}{\kappa_{2}^{3}}+\frac{\pi|\sigma_{\mathrm{g}}|^{2}}{\varepsilon_{0}D_{\mathrm{g}}}\frac{\omega^{2}+v_{\mathrm{s}}^{2}k^{2}}{\omega^{2}}}\;.$
(29) $\Theta_{\mathrm{sp}}=\dfrac{\left\langle
U_{\mathrm{sp}}\right\rangle}{\left\langle U\right\rangle}\\\
=\dfrac{\frac{\pi}{\varepsilon_{0}D_{\mathrm{g}}}\frac{\omega^{2}+v_{\mathrm{s}}^{2}k^{2}}{\omega^{2}}|\sigma_{\mathrm{g}}|^{2}}{\frac{d\frac{\varepsilon_{1}^{2}}{\kappa_{1}^{2}}\frac{\omega^{2}}{c^{2}}}{\sinh^{2}\kappa_{1}d}+k^{2}\frac{\varepsilon_{1}}{\kappa_{1}^{3}}\coth\kappa_{1}d+k^{2}\frac{\varepsilon_{2}}{\kappa_{2}^{3}}+\frac{\pi|\sigma_{\mathrm{g}}|^{2}}{\varepsilon_{0}D_{\mathrm{g}}}\frac{\omega^{2}+v_{\mathrm{s}}^{2}k^{2}}{\omega^{2}}}\;.$
(30)

Fig. 5 shows the variation of $\Theta_{\mathrm{ph}}$ and
$\Theta_{\mathrm{sp}}$ with respect to the dimensionless wavenumber of the
screened plasmon-polaritons of graphene near a PEC. One can see that for low
values of $D$, a plasmonic wave of the system is largely plasmon-like and the
role of the surface electromagnetic wave is small compared with that of the
surface plasmon wave. Only for high values of $D$ the strength functions of
surface plasmon and surface electromagnetic wave are comparable.

##  III Conclusion

In summary, we have investigated the properties of screened plasmons of
graphene near a PEC based on classical electrodynamics and the linearized
hydrodynamic model with Fermi correction. In this way, at first, we have
derived the dispersion relation of the mentioned screened plasmonic waves. We
have studied numerically the effects of graphene distance from the PEC and
graphene sheet carrier density on the surface wave properties of the system.
Numerical results for the present system indicate alterations in the physical
behavior of the surface waves, in comparison with those obtained for an
isolated monolayer graphene. Also, we have derived the analytical expressions
for the power flow, energy density, and energy velocity of screened plasmons
of the system. Furthermore, we have obtained the analytical expressions for
the damping function, and surface plasmon and electromagnetic field strength
functions of surface waves of the system with small intrinsic damping. The
results show that the plasmonic properties of graphene near the PEC may
present new behaviors that make it suitable for applications in the mid-
infrared to a few THz range of frequencies.

## ACKNOWLEDGMENTS

The first-named author would like to thank the Department of Electronic and
Communications Engineering at the Yildiz Technical University for its
hospitality during his visit. Also, A.M. would like to acknowledge the
financial support of the Kermanshah University of Technology for this research
opportunity under grant number S/P/F/6.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

Afshin Moradi: Project administration (equal); Conceptualization (lead);
Investigation (lead); Methodology (lead); Validation (lead); Formal analysis
(lead); Software (lead); Writing - original draft (lead); Writing - review and
editing (equal). Nurhan Türker Tokan: Project administration (equal);
Conceptualization (supporting); Investigation (supporting); Methodology
(supporting); Validation (supporting); Formal analysis (supporting); Writing -
review and editing (equal).

### DATA AVAILABILITY

The data that supports the findings of this study are available within the
article.

## Appendix: The energy density and power flow on graphene with Fermi
correction

From the Poynting theorem for energy in standard electrodynamics, we have:

$\nabla\cdot\textbf{S}+\partial_{t}u=-\textbf{E}\cdot\textbf{J}\;,$ (A-1)

where $\textbf{S}=\textbf{E}\times\textbf{H}$ is known as the Poynting vector,
which is a power density vector associated with an electromagnetic field and
$u=\frac{1}{2}\left(\textbf{E}\cdot\textbf{D}+\textbf{H}\cdot\textbf{B}\right)$
is the energy density of electromagnetic waves. For the present nonmagnetic
system, we have $\textbf{D}=\varepsilon_{0}\varepsilon\textbf{E}$ and
$\textbf{B}=\mu_{0}\textbf{H}$ (where $\mu_{0}$ is the permeability of free
space). At this stage, by employing Eqs. (1) and (2) and after doing some
algebra, we rewrite the right-hand side of Eq. (A-1) on graphene surface as:

$\textbf{E}\cdot\textbf{J}=\frac{\pi}{2D_{\mathrm{g}}}\partial_{t}j_{x}^{2}+\frac{\pi
e^{2}}{2D_{\mathrm{g}}}v_{\mathrm{s}}^{2}\partial_{t}n_{\mathrm{g}}^{2}-\frac{\pi
e}{D_{\mathrm{g}}}v_{\mathrm{s}}^{2}\partial_{x}\textbf{e}_{x}\cdot
n_{\mathrm{g}}j_{x}\textbf{e}_{x}\;.$ (A-2)

Then, this equation together with Eq. (A-1) yields the energy density
$U_{\mathrm{g}}$ and the power flow $S_{\mathrm{g}}$ on the graphene surface
in the forms, as:

$U_{\mathrm{g}}=U_{\mathrm{gk}}+U_{\mathrm{gp}}\;,$ (A-3)
$S_{\mathrm{g}}=-\frac{\pi
e}{D_{\mathrm{g}}}v_{\mathrm{s}}^{2}n_{\mathrm{g}}j_{x}\;.$ (A-4)

On the right-hand side of Eq. (A-3), the first term is the kinetic-energy
density $U_{\mathrm{gk}}=\pi j_{x}^{2}/2D_{\mathrm{g}}$ and the second term
represents the the potential-energy density $U_{\mathrm{gp}}=\pi
e^{2}v_{\mathrm{s}}^{2}n_{\mathrm{g}}^{2}/2D_{\mathrm{g}}$.

## References

  * (1) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Electric field effect in atomically thin carbon films, Science 306, 666 (2004)
  * (2) B. Wunsch, T. Stauber, F. Sols, and F. Guinea, Dynamical polarization of graphene at finite doping, New J. Phys. 8, 318 (2006)
  * (3) E. H. Hwang and S. D. Sarma, Dielectric function, screening, and plasmons in two-dimensional graphene, Phys. Rev. B 75, 205418 (2007)
  * (4) S. A. Mikhailov and K. Ziegler, New electromagnetic mode in graphene, Phys. Rev. Lett. 99, 016803 (2007)
  * (5) X. Luo, T. Qiu, W. Lu, and Z. Ni, Plasmons in graphene: recent progress and applications, Mater. Sci. Eng. R 74, 351 (2013)
  * (6) P.A.D. Goncalves, N.M.R. Peres, An Introduction to Graphene Plasmonic (World Scientific, Singapore, 2016)
  * (7) M. Mendoza, H. J. Herrmann, and S. Succi, Hydrodynamic model for conductivity in graphene, Sci. Rep. 3, 1052 (2013)
  * (8) A. N. Grigorenko, M. Polini, and K. S. Novoselov, Graphene plasmonics, Nat. Photon. 6, 749 (2012)
  * (9) M. Müller, J. Schmalian, and L. Fritz, Graphene: A Nearly Perfect Fluid, Phys. Rev. Lett. 103, 025301 (2009)
  * (10) A. J. Chaves, N. M. R. Peres, G. Smirnov, and N. Asger Mortensen, Hydrodynamic model approach to the formation of plasmonic wakes in graphene, Phys. Rev. B 96, 195438 (2017)
  * (11) B. A. Ferreira, B. Amorim, A. J. Chaves, and N. M. R. Peres, Quantization of graphene plasmons, Phys. Rev. A 101, 033817 (2020)
  * (12) D. Svintsov, V. Vyurkov, S. Yurchenko, T. Otsuji, and V. Ryzhii, Hydrodynamic model for electron-hole plasma in graphene, J. Appl. Phys. 111, 083715 (2012)
  * (13) A. Lucas, Sound waves and resonances in electron-hole plasma, Phys. Rev. B 93, 245153 (2016)
  * (14) Z. Sun, D. N. Basov, and M. M. Fogler, Universal linear and nonlinear electrodynamics of a Dirac fluid, Proc. Natl Acad. Sci. USA 115, 3285 (2018)
  * (15) A. Lucas, and K. C. Fong, Hydrodynamics of electrons in graphene, J. Phys. Condens. Matter 30, 053001 (2018)
  * (16) D. Svintsov, Hydrodynamic-to-ballistic crossover in Dirac materials,Phys. Rev. B 97, 121405 (2018)
  * (17) A. Lucas, and S. D. Sarma, Electronic sound modes and plasmons in hydrodynamic two-dimensional metals, Phys. Rev. B 97, 115449 (2018)
  * (18) I. Torre, L. V. de Castro, B. V. Duppen, D. B. Ruiz, F. M. Peeters, F. H. L. Koppens, and M. Polini, Acoustic plasmons at the crossover between the collisionless and hydrodynamic regimes in two-dimensional electron liquids, Phys. Rev. B 99, 144307 (2019)
  * (19) B. N. Narozhny, Electronic hydrodynamics in graphene, Ann. Phys. 411, 167979 (2019)
  * (20) E. I. Kiselev, and J. Schmalian, Nonlocal hydrodynamic transport and collective excitations in Dirac fluids, Phys, Rev. B 102, 245434 (2020)
  * (21) J. Dufty, K. Luo, and J. Wrighton, Generalized hydrodynamics revisited, Phys. Rev. Res. 2, 023036 (2020)
  * (22) D. Fateev, and V. Popov, Hydrodynamic terahertz plasmons and electron sound in graphene with spatial dispersion, Semiconductors 54, 941 (2020)
  * (23) B. Narozhny, I. Gornyi, and M. Titov, Hydrodynamic collective modes in graphene, Phys. Rev. B 103, 115402 (2021)
  * (24) W. Zhao, S. Wang, S. Chen, Z. Zhang, K. Watanabe, T. Taniguchi, A. Zettl, and F. Wang, Observation of hydrodynamic plasmons and energy waves in graphene, Nature 614, 688 (2023)
  * (25) I. S. Eid, B. F. Mohamed, and B. Guo, Electron exchange effect on surface magnetoplasmon polaritons dynamics in a graphene-plasmonic structure, J. Appl. Phys. 133, 133104 (2023)
  * (26) A. Moradi, Energy density and energy flow of magnetoplasmonic waves on graphene, Solid State Commun. 253, 63 (2017)
  * (27) A. Moradi, Damping properties of plasmonic waves on graphene, Phys. Plasmas 24, 072114 (2017)
  * (28) A. Moradi, Energy density and energy flow of surface waves in a strongly magnetized graphene, J. Appl. Phys. 123, 043103 (2018)
  * (29) A. Moradi, Energy density and energy flow of plasmonic waves in bilayer graphene, Opt. Commun. 394, 135 (2017)
  * (30) A. Moradi, Plasmonic waves of graphene on a conducting substrate, J. Mod. Opt. 66, 353 (2019)
  * (31) X. Gu, I-T. Lin, and J.-M. Liu, Extremely confined terahertz surface plasmon-polaritons in graphene-metal structures, Appl. Phys. Lett. 103, 071103 (2013)
  * (32) M. Y. Morozov, I. M. Moiseenko, and V. V. Popov, Amplification of plasma waves in shielded active graphene. Tech. Phys. Lett. 42, 40 (2016)
  * (33) A. Moradi, Comment on: Tunable surface waves supported by graphene-covered left-handed material structures, Opt. Commun. 545, 129735 (2023)
  * (34) M. M. Bait Suwailam, Z. Chen, Surface waves on a grounded double-negative (DNG) slab waveguide, Microw. Opt. Technol. Lett. 44, 494 (2005)
  * (35) A. S. Petrov, and D. Svintsov, Perturbation theory for two-dimensional hydrodynamic plasmons, Phys. Rev. B 99, 195437 (2019)
  * (36) G. A. Marks, D. Blankespoor, Z. L. Miskovic, Launching plasmons in a two-dimensional material traversed by a fast charged particle, Materials 16, 1150 (2023)
  * (37) A. Moradi, Canonical Problems in the Theory of Plasmonics: From 3D to 2D Systems (Switzerland, Springer, 2020)
  * (38) A. Moradi, M. Wubs, Strongly direction-dependent magnetoplasmons in mixed Faraday-Voigt configurations, Sci Rep 11, 18373 (2021)
  * (39) A. Moradi, N. T. Tokan, Magnetostatic microwaves in circular metallic waveguides filled with uniaxial negative permeability media, J. Appl. Phys. 132, 143901 (2022)
  * (40) A. Moradi, P.-G. Luan, Electromagnetic energy density in hyperbolic metamaterials, Sci. Rep. 12, 10760 (2022)
  * (41) M. Y. Morozov, and V. V. Popov, Concept of terahertz waveguide plasmon amplifier based on a metal groove with active graphene, Sci Rep 12, 22209 (2022)
  * (42) R. Loudon, The propagation of electromagnetic energy through an absorbing dielectric, J. Phys. A: Gen. Phys. 3, 233 (1970)
  * (43) J. Nkoma, R. Loudon, D.R. Tille, Elementary properties of surface polaritons. J. Phys. C: Solid State Phys. 7, 3547 (1974)