  *

# Comprehensive studies on the universality of BKT transitions — Machine-
learning study, Monte Carlo simulation, and Level-spectroscopy method

Hiromi Otsuka1, Kenta Shiina2,1, and Yutaka Okabe1 1Department of Physics,
Tokyo Metropolitan University, Hachioji, Tokyo 192-0397, Japan
2AiCAN Inc., Kanagawa Science Park, 3-2-1 Sakado, Takatsu-ku, Kawasaki
213-0012, Japan<EMAIL_ADDRESS>
###### Abstract

Comprehensive studies are made on the six-state clock universality of two
models using several approaches. We apply the machine-learning technique of
phase classification to the antiferromagnetic (AF) three-state Potts model on
the square lattice with ferromagnetic next-nearest-neighbor (NNN) coupling and
the triangular AF Ising model with anisotropic NNN coupling to study two
Berezinskii-Kosterlitz-Thouless transitions. We also use the Monte Carlo
simulation paying attention to the ratio of correlation functions of different
distances for these two models. The obtained results are compared with those
of the previous studies using the level-spectroscopy method. We directly show
the six-state clock universality for totally different systems with the
machine-learning study.

  * Keywords: BKT transitions, Machine-learning, Monte Carlo, Level-spectroscopy

##  1 Introduction

The concepts of scaling and universality have played an essential role in
statistical physics, especially in the study of phase transitions [2, 3].
There are sometimes implicit symmetries in the models of physics. Universality
appears in totally different systems. One example is a six-state clock
universality.

The two-dimensional (2D) spin systems with a continuous XY symmetry exhibit a
unique phase transition called the Berezinskii-Kosterlitz-Thouless (BKT)
transition [4, 5, 6, 7]. A BKT phase of a quasi-long-range order, a fixed
line, exists, wherein the correlation function decays as a power law. The
$q$-state clock model, which is a discrete version of the XY model,
experiences two BKT transitions for $q>4$ due to the discreteness [8, 9]. It
is noteworthy that in the clock model with a modified interaction of Villain
type [10], there is an exact duality relation between the two (higher and
lower temperatures) transitions.

Antiferromagnetic (AF) spin systems may have a variety of phase transitions
compared to ferromagnetic (F) counterparts. The F three-state Potts model on
the square lattice exhibits a second-order phase transition with $T_{c}$ of
$1/\ln(1+\sqrt{3})$ in units of the coupling $J$. On the contrary, the AF
three-state Potts model on the square lattice does not possess the order,
because there are many possibilities to choose different states for the
nearest-neighbor (NN) pairs. We note that these features are not only present
in the three-state AF Potts model on the square lattice but are also present
in the three-state AF Potts model on a large class of plane quadrangulations
[11]. There is a macroscopic degeneracy in the ground state. However, if the F
next-nearest-neighbor (NNN) coupling is added [12], diagonal spins take the
same state. That is, all the spins in the sublattice ”a” take the same state,
one of three, and all the spins in the sublattice ”b” take one of the other
two states. Then, the degeneracy of the ground state becomes six ($3\times
2$). This model is thought to be in the same universality class as the F six-
state clock model in the analysis of essentially the same model with three-
fold symmetry-breaking fields [13]. Another example is the AF Ising model on
the triangular lattice with an anisotropic NNN coupling, which was introduced
by Kitatani and Oguchi [14]. A detailed description of the anisotropic NNN
coupling will be given in Sec. II. There is no order in the AF Ising model on
the triangular lattice without NNN couplings due to frustration, and a
macroscopic degeneracy appears in the ground state. With the introduction of
the anisotropic NNN coupling, spins are classified into interpenetrating three
sublattices, and the system is considered to belong to the six-state clock
universality class [14].

There is a difficulty in the numerical analysis of BKT transitions. The
divergence of the correlation length for the BKT transition is more rapid than
any power law, and there are logarithmic corrections. These pathological
features make it difficult to determine the BKT transition point from finite-
size calculations. Nomura proposed a method to determine BKT points observed
in quantum spin chains with high precision [15], the so-called level
spectroscopy (LS) method [16]. It was pointed out that an enhancement of the
U(1) symmetry of the 2D Gaussian critical model to an SU(2) symmetry of a
Wess-Zumino-Witten model brings about a typical degeneracy in an excitation
spectrum, which then can be detected as a crossing of excitation levels of
constituent physical quantities to an SU(2) multiplet [15]. Although the
levels to analyze depend on the types of BKT transitions, we can specify them
in numerical calculations of finite-size systems with the help of the symmetry
properties of physical quantities and thus can determine BKT transition
points. Nakamura et al. extended the LS method to apply to the one-dimensional
electron systems; for example, see [17]. Alternatively, Otsuka et al. studied
a 2D classical spin model with frustration, the AF three-state Potts model on
the square lattice with F NNN coupling [18]. Combining the exact
diagonalization calculations and the LS analysis, they showed two BKT
transitions and the universality with the F six-state clock model (see Sec.
II). The precise phase diagram in the space of the NN coupling and the NNN
coupling was obtained. The triangular AF Ising model with anisotropic NNN
coupling was also studied by Otsuka et al. [19] by the LS method. The six-
state clock universality was discussed. Although some numerical works on these
two models had been reported previously, the LS studies [18, 19] gave reliable
global phase diagrams.

The Monte Carlo (MC) simulation has served as a standard method for the
numerical analysis of many particle systems [20]. The finite-size scaling
(FSS) study [21, 22] of the Binder ratio [23], essentially the moment ratio,
is a powerful tool for the second-order phase transition. However, for the BKT
transition, the Binder ratio is not so effective [24] because of the
multiplicative logarithmic corrections [7, 25]. There are other quantities
that have the same FSS behavior with a single variable as the Binder ratio.
The ratio of the correlation functions of different distances [26] is one
example. The size-dependent second-moment correlation length for a finite
system, $\xi(L)$, has been frequently used in spin glass problems [27]. The
FSS of $\xi(L)/L$ also has the same form of a single scaling variable.
Surungan et al. [28] made the MC study of the correlation ratio and the size-
dependent second-moment correlation length for the $q$-state clock model of
both the standard cosine-type interaction and the Villain-type interaction.
They observed the collapsing curves of different sizes at intermediate
temperature regime and the spray out at lower and higher temperatures, which
demonstrates two BKT transitions.

In addition to conventional approaches, new numerical strategies are
developing. Recent advances in machine-learning-based techniques have been
applied to fundamental research, such as statistical physics [29].
Carrasquilla and Melko [30] used a technique of supervised learning for image
classification, which is complementary to the conventional approach of
studying interacting spin systems. They classified and identified a high-
temperature paramagnetic phase and a low-temperature ferromagnetic phase of
the 2D Ising model by using data sets of spin configurations. Shiina et al.
[31] extended and generalized this machine-learning approach to studying
various spin models including the multi-component systems and the systems with
a vector order parameter. The configuration of a long-range spatial
correlation was treated instead of the spin configuration itself. Not only the
second-order and the first-order transitions but also the BKT transition was
studied. They detected two BKT transitions for the six-state clock model.
Tomita et al. [32] made further progress in the machine-learning study of
phase classification. When the cluster update is possible in the MC
simulation, the Fortuin-Kasteleyn (FK) [33, 34] representation-based improved
estimators [35, 36] for the configuration of two-spin correlation were
employed as an alternative to the ordinary spin configuration. This method of
improved estimators was applied not only to the classical spin models but also
to the quantum MC simulation. Miyajima et al. also used a machine-learning
approach to study two BKT transitions of the eight-state clock model [37]. We
note that the principal component analysis was used in the machine-learning
study on the XY model [38, 39]. As another approach to machine learning, the
image super-resolution using the deep convolutional neural networks was
applied to the study of the inverse renormalization group of spin systems
[40].

Quite recently, much attention has been paid to the study of the $q$-state
clock models [41, 42, 43, 44, 45]. Some works used newly developed methods,
such as the tensor renormalization methods [41, 42, 44], the corner-transfer
matrix renormalization group [43], and the worm-type simulation [45].

In the present paper, we perform the comprehensive studies on two spin models,
the AF three-state Potts model on the square lattice with F NNN coupling and
the triangular AF Ising model with anisotropic NNN coupling. We first apply
the machine-learning approach to these models. Using the training data of the
F six-state clock model, we demonstrate the classification of the disordered,
the BKT, and the ordered phase. We then perform the MC studies on two models.
Analyzing the correlation ratios, we examine two BKT transitions. Finally, we
compare the results with the LS method by Otsuka et al. [18, 19].

The rest of the paper is organized as follows: In Sec. II, we explain the
models and simulation models. Sections III and IV are devoted to the results
of the AF three-state Potts model on the square lattice with F NNN coupling
and the triangular AF Ising model with anisotropic NNN coupling, respectively.
The summary and discussion are given in Sec. V. The calculations of the XY and
six-state clock models on the triangular lattice are presented in the
Appendix.

##  2 Models and numerical methods

We treat two spin models in the present study. The first one is the AF three-
state Potts model on the square lattice with F NNN coupling [12], whose
Hamiltonian is given by

$H=J_{1}\sum_{\langle
ij\rangle}\delta_{s_{i}s_{j}}-J_{2}\sum_{[ij]}\delta_{s_{i}s_{j}},\quad
s_{i}=1,2,3,$ (1)

where $\delta_{ab}$ is the Kronecker delta. The first and the second sums run
over all NN and NNN pairs, respectively. The periodic boundary conditions are
imposed in numerical simulations. The couplings $J_{1}$ and $J_{2}$ are
positive, and hereafter we use $r=J_{2}/J_{1}$. The temperature will be
measured in units of $J_{1}$.

The second model is the AF Ising model on the triangular lattice with
anisotropic NNN coupling, which was studied by Kitatani and Oguchi [14]. The
Hamiltonian is given by

$H=J_{1}\sum_{\langle
ij\rangle}\delta_{s_{i}s_{j}}-J_{2}\sum_{[ij]^{\prime}}\delta_{s_{i}s_{j}},\quad
s_{i}=1,2,$ (2)

where the first sum runs over all NN pairs, whereas the second sum runs over
anisotropic NNN pairs in two of the three directions. The schematic
representation of Hamiltonian is given in figure 1. The periodic boundary
conditions are imposed in numerical simulations. The couplings $J_{1}$ and
$J_{2}$ are positive, and we use $r=J_{2}/J_{1}$. The temperature will be
measured in units of $J_{1}$ again.

Figure 1:  The schematic representation of Hamiltonian (2) on the triangular
lattice [19]. The NN couplings are denoted by straight lines, whereas the
anisotropic NNN couplings in two of the three directions are denoted by dotted
lines.

The correlation function with the distance $r$ for the $q$-state Potts model
including the Ising model ($q=2$) is given by

$g_{i}(r)=\frac{q\delta_{s_{i}s_{i+r}}-1}{q-1}.$ (3)

It assumes a value of $+1$ or $-1/(q-1)$.

We use three approaches to study the BKT transitions. The first approach is
the machine-learning study which was developed by Shiina et al. [31], to
classify the ordered, the BKT, and the disordered phases. A fully connected
neural network is implemented with a standard library of TensorFlow of the
100-hidden unit model. For the spin data, we use the long-range correlation
configuration with the distance of $L/2$ instead of the spin configuration
itself. The average of the correlation along the $x$-axis and the $y$-axis is
considered. For the training data, we use the results of the F six-state clock
model [31]. For the test data, we investigate the correlation configurations
of the AF three-state Potts model with F NNN coupling and that of the
triangular AF Ising model with anisotropic NNN coupling. We classify the
phases of these two models using the training data of the six-state clock
model.

As the second approach, we perform the MC simulation of the ratio of the
correlation functions of different distances, $R(T)=\langle
g(L/2)\rangle/\langle g(L/4)\rangle$; for the two distances, we chose $L/2$
and $L/4$. Surungan et al. [28] made MC studies of the BKT transitions by
investigating the correlation ratio and the size-dependent second-moment
correlation length. Although there are other candidate quantities, we take the
correlation ratio here. The correlation ratio has a single scaling variable
for FSS

$R(T)=\frac{\langle g(L/2)\rangle}{\langle g(L/4)\rangle}=\tilde{f}(L/\xi),$
(4)

as in the Binder ratio [23], which is essentially the ratio of different
moments. Here, $\xi$ stands for the correlation length. At the critical
region, where the correlation length $\xi$ is infinite, the correlation ratio
does not depend on the system size $L$. Thus, we expect the data collapsing of
different sizes in the BKT phase. Above the upper BKT temperature $T_{2}$ and
below the lower BKT temperature $T_{1}$, the data of different sizes start to
separate. In the case of the BKT transition, the correlation length diverges
as

$\xi\propto\exp(c/\sqrt{|t|}),$ (5)

where $t=T-T_{1,2}$. We employed the Metropolis local-update algorithm for the
two models in combination with the replica-exchange algorithm [46] that
overcomes slow dynamics. For the clock model we used the Swendsen-Wang
cluster-update algorithm [47] in combination with the replica-exchange
algorithm. The typical number of Monte Carlo steps was 200000 after discarding
the first 10000 steps. Since the temperature exchange was typically performed
for 40 temperatures, the effective Monte Carlo step number is much larger.

The third approach is the LS method employed in the previous publications [18,
19]. To make the explanations concrete, we shall borrow some notations there.
Based on a dual sine-Gordon model, an effective field theory of the lattice
models (see (2) in [18] and (4) in [19]), we analyzed the scaling dimensions
of some local operators. For the upper BKT transition, the so-called ${\cal
M}$ operator (a part of the Gaussian term) hybridizes with an operator
$\sqrt{2}\cos\\!\sqrt{2}\theta$ given in the disorder field and yields the
${\cal M}$-like operator with the dimension ${\bar{x}}_{0}$ [16]. We also
focused on an operator $\sqrt{2}\sin 3\sqrt{2}\phi$ given in the order field
with the dimension $x_{2}$. The numerical analysis of the eigenvalues of
transfer matrices supposes a lattice model with a $L\times\infty$ cylindrical
geometry. As explained in Refs. [18, 19], we can evaluate the renormalized
scaling dimensions ${\bar{x}}_{0}(l)$ and $x_{2}(l)$ ($l=\ln L$) using
discrete symmetry properties of these operators. We then obtained finite-size
estimates of $T_{2}$ by the level-crossing condition

${\bar{x}}_{0}(l)=\frac{16}{9}x_{2}(l).$ (6)

For the lower BKT transition, the ${\cal M}$ operator hybridizes with
$\sqrt{2}\cos 6\sqrt{2}\phi$ and yields another ${\cal M}$-like operator with
$x_{0}$ [16]. We thus calculated finite-size estimates of $T_{1}$ by the
level-crossing condition

$x_{0}(l)=4x_{2}(l).$ (7)

We compare the estimated BKT temperatures $T_{1,2}$ by three approaches.

##  3 AF three-state Potts model with ferromagnetic NNN coupling

### Machine-learning study

We first show the results of the AF three-state Potts model with F NNN
coupling. As a machine-learning study, the output layer averaged over a test
set as a function of $T$ is plotted in figure 2. The data for $r$ = 0.2, 0.4,
and 0.8 are given in (a), (b), and (c), respectively. The system sizes are 32,
48, and 64. For the training data, we use the F six-state clock model. We show
the result of the previous study [31] for the F six-state clock model in
figure 3, for convenience. The probabilities of predicting the phases are
plotted at each temperature. The samples of $T$ within the ranges $0.4\leq
T\leq 0.64$, $0.78\leq T\leq 0.82$ and $0.96\leq T\leq 1.2$ are used for the
low-temperature, middle-temperature, and high-temperature training data,
respectively. For a whole temperature range, around 35000 training data sets
are used, and we use 500 test data sets for each temperature.

Figure 2:  The output layer for the AF three-state Potts model with F NNN
coupling for $r=0.2$ (a), $r=0.4$ (b), and $r=0.8$ (c), using the training
data of the F six-state clock model. The system sizes are $L=32$, 48, and 64.
Figure 3:  The output layer averaged over a test set as a function of $T$ for
the 2D six-state clock model on the square lattice. The system sizes are $L$ =
32, 48, and 64. The samples of $T$ within the ranges $0.4\leq T\leq 0.64$,
$0.78\leq T\leq 0.82$ and $0.96\leq T\leq 1.2$ are used for training data.

Figure 2 clearly shows the existence of three phases, the ordered phase, the
intermediate phase, and the disordered phase. We estimate the size-dependent
$T_{1,2}(L)$ from the point that the probabilities of predicting two phases
are 50%. Systematic size dependence is not observed. The rough estimates of
the transition temperatures, $T_{1}$ and $T_{2}$, are tabulated in table 1. We
emphasize that the BKT transitions of a complex model, AF three-state Potts
model on the square lattice with F NNN coupling can be studied by the training
data of the six-state clock model.

Figure 4:  The plot of the correlation ratio $R(T)$ for the AF three-state
Potts model on the square lattice with F NNN coupling. In the inset, the FSS
plots are given, where $X(c,T)=L/\exp(c_{1,2}/\sqrt{|T-T_{1,2}|})$.

### Monte Carlo study

We here show the MC results of the AF three-state Potts model on the square
lattice with F NNN coupling. Following the study by Surungan et al. [28], we
calculate the correlation ratio, $R(T)=\langle g(L/2)\rangle/\langle
g(L/4)\rangle$, to investigate the BKT transitions. The results for $r$ = 0.2,
0.4, and 0.8 are plotted in figure 4. The system sizes are $L=32$, 48, 64, 96,
and 128. We observe the collapsing curves of different sizes at intermediate
temperature regimes and the spray out at lower and higher temperatures, which
indicates the existence of two BKT transitions. The behavior of the collapsing
is not good enough, compared to the regular six-state clock model [28], which
may be attributed to the large corrections to FSS.

In the case of the BKT transition, the correlation length diverges as

$\xi\propto\exp(c/\sqrt{|t|}),$ (8)

where $t=T-T_{1,2}$. Thus, we can get the FSS plot if the temperature $T$ is
scaled as $L/\exp(c_{1,2}/\sqrt{|T-T_{1,2}|})$, where $c$ is a fitting
parameter. The FSS plots are given in the inset of figure 4.

The rough estimates of $T_{1}$ and $T_{2}$ are also tabulated in table 1.

### Level-spectroscopy study

This model was studied using the LS method [18]. The phase diagram, in the
space of $u=e^{-K_{1}}=e^{-J_{1}/T}$ and $v=e^{K_{2}}=e^{rJ_{1}/T}$, was given
in figure 2 of [18]. The phase diagram is reproduced in figure 5, and the
curves $v=e^{-r\ln(u)}$ with $r=0.2$, 0.4, and 0.8 are added. From the
crossing point with these curves, we can read off the estimates of $T_{1}$ and
$T_{2}$. The estimates are tabulated in table 1.

Figure 5:  The phase diagram of the AF three-state Potts model on the square lattice with F NNN coupling using the LS method [18] in the space of $(u,v)$ = $(e^{-K_{1}},e^{K_{2}})$. The curves $v=e^{-r\ln(u)}$ with $r=0.2$, 0.4, and 0.8 are also plotted. From the crossing point with these curves, we can obtain the estimates of $T_{1}$ and $T_{2}$.  Table 1:  The estimates of $T_{1}$ and $T_{2}$ for the AF three-state Potts model on the square lattice with F NNN coupling. The estimates of the LS method are based on the calculation in [18].  | $r=0.2$ | $r=0.4$ | $r=0.8$
---|---|---|---
| $T_{1}$ | $T_{2}$ | $T_{1}$ | $T_{2}$ | $T_{1}$ | $T_{2}$
machine-learning | 0.25 | 0.44 | 0.50 | 0.74 | 0.98 | 1.23
Monte Carlo | 0.30 | 0.51 | 0.57 | 0.775 | 1.075 | 1.255
level-spectroscopy | 0.277 | 0.490 | 0.554 | 0.755 | 1.056 | 1.240

By comparing various estimates of $T_{1}$ and $T_{2}$ tabulated in table 1, we
see the consistency of various methods. The machine-learning study using the
training data of the six-state clock model yields reasonable estimates.

##  4 Triangular AF Ising model with anisotropic NNN coupling

### Machine-learning study

Next, we consider the triangular AF Ising model with anisotropic NNN coupling.
Because the $q$-state clock models (XY model in the limit of $q\to\infty$) on
the triangular lattice were not studied previously in the literature, the
detailed calculation of the clock models on the triangular lattice is given in
the Appendix.

As the result of the machine-learning study, we show the output layer averaged
over a test set as a function of $T$ for the triangular AF Ising model with
anisotropic NNN coupling in figure 6. For the training data, we use the F six-
state clock model on the triangular lattice, which will be shown in the
Appendix. The data for $r=0.2$, 0.4, and 0.8 are given in (a), (b), and (c),
respectively. The system sizes are 24, 36, and 48. The probabilities of
predicting the phases are plotted at each temperature. Figure 6 clearly
demonstrates the existence of three phases, the ordered phase, the
intermediate BKT phase, and the disordered phase. The rough estimates of the
transition temperatures, $T_{1}$ and $T_{2}$, are tabulated in table 2.

Figure 6:  The output layer for the triangular AF Ising model with anisotropic
NNN coupling for $r=0.2$ (a), $r=0.4$ (b), and $r=0.8$ (c), using the training
data of the F six-state clock model. The system sizes are $L$ = 24, 36, and
48.

Figure 7:  The plot of the correlation ratio $R(T)$ for the triangular AF
Ising model with anisotropic NNN coupling. In the inset, the FSS plots are
given, where $X(c,T)=L/\exp(c_{1,2}/\sqrt{|T-T_{1,2}|})$.

### Monte Carlo study

We show the MC results of the correlation ratio, $R(T)$, of the triangular AF
Ising model with anisotropic NNN coupling. The results for $r=0.2$, 0.4, and
0.8 are plotted in figure 7. The system sizes are $L=24$, 48, 72, and 96. We
again observe the collapsing curves of different sizes at intermediate
temperature regimes and the spray out at lower and higher temperatures. The
behavior of the collapsing in the BKT phase is better than the case of the AF
three-state Potts model on the square lattice with F NNN coupling (figure 4).

The FSS plots based on the exponential divergence of the correlation length
are given in the inset of figure 7, where $T$ is scaled as
$L/\exp(c_{1,2}/\sqrt{|T-T_{1,2}|})$.

The rough estimates of $T_{1}$ and $T_{2}$ are also tabulated in table 2.

### Level-spectroscopy study

Otsuka et al. [19] studied this model using the LS method. The phase diagram,
in the space of $u=e^{-K_{1}}=e^{-J_{1}/T}$ and $v=e^{K_{2}}=e^{rJ_{1}/T}$, is
reproduced in figure 8, where the curves $v=e^{-r\ln(u)}$ with $r=0.2$, 0.4,
and 0.8 are added. From the crossing point with the curve $v=e^{-r\ln(u)}$, we
can obtain the estimates of $T_{1}$ and $T_{2}$. The estimates are tabulated
in table 2.

Figure 8:  The phase diagram of the triangular AF Ising model with anisotropic NNN coupling using the level-spectroscopy method [19] in the space of $(u,v)$ = $(e^{-K_{1}},e^{K_{2}})$. The curves $v=e^{-r\ln(u)}$ with $r$ = 0.2, 0.4, and 0.8 are also plotted. From the crossing point with these curves, we can obtain the estimates of $T_{1}$ and $T_{2}$.  Table 2:  The estimates of $T_{1}$ and $T_{2}$ for the triangular AF Ising model with anisotropic NNN coupling. The estimates of the level-spectroscopy method are based on the calculation in [19].  | $r=0.2$ | $r=0.4$ | $r=0.8$
---|---|---|---
| $T_{1}$ | $T_{2}$ | $T_{1}$ | $T_{2}$ | $T_{1}$ | $T_{2}$
machine-learning | 0.25 | 0.48 | 0.51 | 0.82 | 1.01 | 1.39
Monte Carlo | 0.29 | 0.465 | 0.57 | 0.795 | 1.08 | 1.36
level-spectroscopy | 0.271 | 0.494 | 0.547 | 0.801 | 1.070 | 1.330

We see the consistency of the estimates of $T_{1}$ and $T_{2}$ by various
methods tabulated in table 2. We may say that the machine-learning study of
this model yields reasonable estimates.

##  5 Effective exponents

In this section, we make a comment on the calculation of the decay exponent
$\eta$ in the MC study, and discuss the self-dual point.

The correlation function with the distance $r$ decays as $g(r)\sim
r^{-\eta(T)}$ at the critical point (line). We consider the correlation
function with the distance $L/2$ for the finite system of the size $L$

$g_{2}(L)=\langle\bm{s}(i)\cdot\bm{s}(i+L/2)\rangle.$ (9)

We also consider

$g_{4}(L)=\langle{\bm{s}}(i)\cdot{\bm{s}}(i+L/4)\rangle.$ (10)

At the critical point (line), $g_{2}(L)$ behaves as $L^{-\eta(T)}$ and for
$T<T_{1}$ and $T>T_{2}$, we may have the following relation:

$g_{2}(L)\sim L^{-\eta(T)}f(L/\xi).$ (11)

We examine the ratio of $g_{2}$ of the two sizes $L$ and $2L$:

$\frac{g_{2}(L)}{g_{2}(2L)}=\frac{L^{-\eta}f(L/\xi)}{(2L)^{-\eta}f(2L/\xi)}=2^{-\eta}\frac{f(L/\xi)}{f(2L/\xi)}.$
(12)

Then, we define the effective exponent $\eta_{\rm eff}$ by

$\eta_{\rm eff}=\frac{\ln(g_{2}(L)/g_{2}(2L))}{\ln(2)}.$ (13)

Theoretically, the critical exponent $\eta$ becomes 1/4 at $T_{2}$ and
$4/q^{2}$ at $T_{1}$. We may estimate $T_{2}$ and $T_{1}$ from the condition
that $\eta_{\rm eff}$ takes the theoretical value. The exponent $\eta_{\rm
eff}$ is expected to take the temperature-dependent value for $T_{1}\leq T\leq
T_{2}$; it does not depend on $L$, and there is no difference between $g_{2}$
and $g_{4}$. For $T<T_{1}$ and $T>T_{2}$, it becomes $L$-dependent, and
$g_{2}$ and $g_{4}$ take different values. We show the calculated results of
$\eta_{\rm eff}$ using (13) for the AF three-state Potts model with
ferromagnetic NNN coupling in figure 9 and for the AF Ising model with
anisotropic NNN coupling on the triangular lattice in figure 10, respectively.
There, the theoretical values of $\eta_{2}=1/4$ and $\eta_{1}=1/9$ are shown
by the dashed and dotted lines, respectively. Not only the data for $g_{2}$
but those for $g_{4}$ are given. We see an expected behavior of the
temperature dependence of the effective exponent $\eta_{\rm eff}$ in figures 9
and 10. The temperatures at which $\eta_{1,2}$ take the theoretical values are
consistent with the estimates of $T_{1}$ and $T_{2}$ in tables 1 and 2.

In the Villain model, where the exact duality transformation holds, the
exponent $\eta$ becomes $1/q$ at the self-dual (sd) point. Moreover, there is
an exact relation between the temperatures:

$T_{1}T_{2}=4\pi^{2}/q^{2},\quad{\rm therefore},\quad T_{\rm sd}=2\pi/q.$ (14)

Although this self-dual point is particular to the Villain model, we give the
value of $\eta_{\rm eff}$, $1/q$, in figures 9 and 10, by the dash-dotted
lines. Let us denote the temperature $T_{\rm sd}$ by the condition that
$\eta_{\rm eff}$ takes the self-dual value $\eta_{\rm sd}$. We can see that
$T_{\rm sd}$ is close to the temperature at which the probability of
predicting the BKT phase takes the maximum in Figs. 2 and 6 in the machine-
learning studies. We note that for the triangular AF Ising model with
anisotropic NNN coupling, the self-dual temperature was discussed using the LS
method [19].

Figure 9:  The effective exponents $\eta_{\rm eff}$, calculated by (13), for
AF three-state Potts model with F NNN coupling (square lattice). The
theoretical values of $\eta_{2}=1/4$, $\eta_{\rm sd}=1/6$, and $\eta_{1}=1/9$,
are shown by the dashed, dash-dotted, and dotted lines, respectively.

Figure 10:  The effective exponents $\eta_{\rm eff}$, calculated by (13), for
the AF Ising model with anisotropic NNN coupling (triangular lattice). The
theoretical values of $\eta_{2}=1/4$, $\eta_{\rm sd}=1/6$, and $\eta_{1}=1/9$
are shown by the dashed, dash-dotted, and dotted lines, respectively.

##  6 Summary and discussion

We used the machine-learning study for the AF three-state Potts model on the
square lattice with F NNN coupling and the triangular AF Ising model with
anisotropic NNN coupling. We classified the ordered, the BKT, and the
disordered states using the training data of the F six-state clock model. We
explicitly showed the universal behavior of implicit symmetries in totally
different models.

We also used the MC studies paying attention to the correlation ratio, which
is appropriate in studying the BKT transitions. In the recent publication on
the antiferromagnetic Potts model [48], the authors pointed out that the
Binder ratio is not necessarily the best method in determining the critical
point as in Refs. [27] and [24]. The two BKT transitions were determined with
precision using the valid metric of correlation ratios, confirming the six-
state clock universality.

Comparing the results of the machine-learning study and the MC simulation with
the previous results of the LS method [18, 19], we confirmed the six-state
clock universality for the AF three-state Potts model on the square lattice
with F NNN coupling and the triangular AF Ising model with anisotropic NNN
coupling.

We make a brief comment on a future study. Surface effects are important
subjects of critical phenomena, and the extraordinary-log phase has captured
recent attention. With this regard, recently the three-dimensional six-state
clock model [49] and the three-dimensional three-state antiferromagnetic Potts
model [48] with a free surface were studied. The application of the present
approaches to these models will be interesting problems.

Universality is an essential concept in phase transitions. We performed a
comprehensive study on the six-state clock universality in completely
different models. In future research, we would like to explore the six-state
clock universality in other models.

The authors thank Hiroyuki Mori and Hwee Kuan Lee for valuable discussions.
This work was supported by JSPS KAKENHI Grant Number JP22K03472.

## Appendix. $\bm{XY}$ and six-state clock models on the triangular lattice

The reports on the XY or clock models on the triangular lattice have not been
found in the literature. Only the exception is the MC study of the XY model by
Sorokin [50]. By using the MC study of the helicity modulus, the BKT
temperature was estimated as 1.418(2). We here show our MC results on the XY
model and the six-state clock model on the triangular lattice. We also perform
a machine-learning study.

Figure 11:  Helicity modulus (a) and correlation ratio (b) for XY model on the
triangular lattice. The system sizes are 48, 72, 96, 144, and 192. In the plot
of helicity modulus, we give the straight line $(2/\pi)T$. In the inset the
FSS plots are given, where $X(c,T)=L/\exp(c_{2}/\sqrt{|T-T_{2}|})$.

The MC results on the XY model and the six-state clock model are plotted in
Figs. 11 and 12, respectively. The helicity modulus and correlation ratio
$R(T)$ are plotted in both figures. The system sizes are 48, 72, 96, 144, and
192. In the plot of helicity modulus, we give the straight line $(2/\pi)T$.
The crossing point gives a universal jump. The $q$-state clock model, a
discrete version of the XY model, experiences two BKT transitions because of
the discreteness. We note that there is no anomaly in helicity modulus for the
lower transition, $T_{1}$. The BKT phase is indicated by the collapsing curves
of different sizes in the plot of $R(T)$. It is a direct consequence of the
power law behavior of the correlation function at the BKT phase. The spray out
of the curves at lower and higher temperatures signifies BKT transitions. The
rough estimates of the BKT temperatures are $T_{2}$ = 1.43 for the XY model,
and $T_{2}$=1.44 and $T_{1}$ = 1.12 for the six-state clock model. The result
of the XY model is consistent with the report by Sorokin [50]. The values of
the triangular lattice are roughly 3/2 of those of the square lattice. It is
because the coordination number of the triangular lattice is six, whereas that
of the square lattice is four.

Figure 12:  Helicity modulus (a) and correlation ratio (b) for the six-state
clock model on the triangular lattice. The system sizes are 48, 72, 96, 144,
and 192. In the plot of helicity modulus, we give the straight line
$(2/\pi)T$. In the inset the FSS plots are given, where
$X(c,T)=L/\exp(c_{1,2}/\sqrt{|T-T_{1,2}|})$. For the helicity modulus, only
the FSS for the high-temperature side is given.  Figure 13:  The output layer
averaged over a test set as a function of $T$ for the triangular lattice six-
state clock model. The system sizes are $L=24$, 36, and 48. The samples of $T$
within the ranges $0.76\leq T\leq 1.04$, $1.21\leq T\leq 1.31$ and $1.52\leq
T\leq 1.8$ are used for training data.

We also show the results of the machine-learning study. The output layer
averaged over a test set as a function of $T$ for the triangular F six-state
model is shown in figure 13. The system sizes are 24, 36, and 48. The samples
of $T$ within the ranges $0.76\leq T\leq 1.04$, $1.21\leq T\leq 1.31$ and
$1.52\leq T\leq 1.8$ are used for the low-temperature, middle-temperature, and
high-temperature training data, respectively. The probabilities of predicting
the phases are plotted at each temperature. We clearly observe three phases,
the ordered, the BKT, and the disordered phases. We estimate the size-
dependent $T_{1,2}(L)$ from the point that the probabilities of predicting two
phases are 50%. Systematic size dependence is not observed. The rough
estimates of $T_{1}$ and $T_{2}$ are 1.07 and 1.48, respectively. The
estimates are compatible to those of the MC study. However, the BKT phase is
slightly wider, which is the same situation as the case of the square lattice.
It comes from the finite-size effects.

## References

  * [1]
  * [2] L. P. Kadanoff, W. Götze, D. Hamblen, R. Hecht, E. A. S. Lewis, V. V. Palciauskas, M. Rayl, J. Swift, D. Aspnes, and J. Kane 1967 Static Phenomena Near Critical Points: Theory and Experiment Rev. Mod. Phys. 39 395\.
  * [3] L. P. Kadanoff 1990 Scaling and universality in statistical physics Physica A 163 1\.
  * [4] V. L. Berezinskii 1970 Destruction of Long-range Order in One-dimensional and Two-dimensional Systems having a Continuous Symmetry Group I. Classical Systems Sov. Phys. JETP 32 493\.
  * [5] V. L. Berezinskii 1972 Destruction of Long-range Order in One-dimensional and Two-dimensional Systems Possessing a Continuous Symmetry Group. II. Quantum Systems Sov. Phys. JETP 34 610\.
  * [6] J. M. Kosterlitz and D. Thouless 1973 Ordering, metastability and phase transitions in two-dimensional systems J. Phys. C: Solid State Phys. 6 1181\.
  * [7] J. M. Kosterlitz 1974 The critical properties of the two-dimensional xy model J. Phys. C: Solid State Phys. 7 1046\.
  * [8] J. V. José, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson 1977 Renormalization, vortices, and symmetry-breaking perturbations in the two-dimensional planar model Phys. Rev. B 16 1217\.
  * [9] S. Elitzur, R. B. Pearson, and J. Shigemitsu 1979 Phase structure of discrete Abelian spin and gauge systems Phys. Rev. D 19 3698\.
  * [10] J. Villain 1975 Theory of one- and two-dimensional magnets with an easy magnetization plane. II. The planar, classical, two-dimensional magnet J. Phys. (Paris) 36 581\.
  * [11] J.-P. Lv, Y. Deng, J. L. Jacobsen, and J. Salas 2018 The three-state Potts antiferromagnet on plane quadrangulations J. Phys. A: Math. Theor. 51 365001\.
  * [12] M. P. M. den Nijs, M. P. Nightingale, and M. Schick 1982 Critical fan in the antiferromagnetic three-state Potts model Phys. Rev. B 26 2490\.
  * [13] J. L. Cardy 1981 Antiferromagnetic clock models in two dimensions Phys. Rev. B 24 5128\.
  * [14] H. Kitatani and T. Oguchi 1988 Antiferromagnetic Triangular Ising Model with Ferromagnetic Next Nearest Neighbor Interactions-Transfer Matrix Method- J. Phys. Soc. Jpn. 57 1344\.
  * [15] K. Nomura and K. Okamoto 1994 Critical properties of $S=1/2$ antiferromagnetic XXZ chain with next-nearest-neighbour interactions J. Phys. A: Math, Gen. 27 5773\.
  * [16] K. Nomura 1995 Correlation functions of the 2D sine-Gordon model J. Phys. A: Math, Gen. 28 5451\.
  * [17] M. Nakamura, K. Nomura, and, A. Kitazawa 1997 Renormalization Group Analysis of the Spin-Gap Phase in the One-Dimensional $t$-$J$ Model Phys. Rev. Lett. 79 3214\.
  * [18] H. Otsuka, K. Mori, Y. Okabe, and K. Nomura 2005 Level spectroscopy of the square-lattice three-state Potts model with a ferromagnetic next-nearest-neighbor coupling Phys. Rev. E 72 046103\.
  * [19] H. Otsuka, Y. Okabe, and K. Nomura 2006 Global phase diagram and six-state clock universality behavior in the triangular antiferromagnetic Ising model with anisotropic next-nearest-neighbor coupling: Level-spectroscopy approach Phys. Rev. E 74 011104\.
  * [20] D. P. Landau and K. Binder 2014 A guide to MC Simulations in Statistical Physics (Cambridge: Cambridge University Press).
  * [21] M. N. Barber 1983 in Phase Transitions and Critical Phenomena, vol. 8 ed C. Domb and J. Lebowitz (New York: Academic Press).
  * [22] J. L. Cardy 1988 Finite Size Scaling ed J. L. Cardy (Amsterdam: North-Holland).
  * [23] K. Binder 1981 Finite size scaling analysis of ising model block distribution functions Z. Phys. B 43 119\.
  * [24] M. Hasenbusch 2008 The Binder cumulant at the Kosterlitz-Thouless transition J. Stat. Mech. P08003.
  * [25] W. Janke 1997 Logarithmic corrections in the two-dimensional XY model Phys. Rev. B 55 3580\.
  * [26] Y. Tomita and Y. Okabe 2002 Finite-size scaling of correlation ratio and generalized scheme for the probability-changing cluster algorithm Phys. Rev. B 66 180401(R).
  * [27] H. G. Katzgraber, M. Körner, and A. P. Young 2006 Universality in three-dimensional Ising spin glasses: A Monte Carlo study Phys. Rev. B 73 224432\.
  * [28] T. Surungan, S. Masuda, Y. Komura, and Y. Okabe 2019 Berezinskii-Kosterlitz-Thouless transition on regular and Villain types of $q$-state clock models J. Phys. A: Math. Theor. 52 275002\.
  * [29] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová 2019 Rev. Mod. Phys. 91 045002\.
  * [30] J. Carrasquilla and R. G. Melko 2017 Machine learning phases of matter Nat. Phys. 13 431\.
  * [31] K. Shiina, H. Mori, Y. Okabe, and H. K. Lee 2020 Machine-Learning Studies on Spin Models Sci. Rep. 10 2177\.
  * [32] Y. Tomita, K. Shiina, Y. Okabe, and H. K. Lee 2020 Machine-learning study using improved correlation configuration and application to quantum Monte Carlo simulation Phys. Rev. E 102 021302(R).
  * [33] P. W. Kasteleyn and C. M. Fortuin 1969 Phase Transitions in Lattice Systems with Random Local Properties J. Phys. Soc. Jpn. Suppl. 26 11\.
  * [34] C. M. Fortuin and P. W. Kasteleyn 1972 On the random-cluster model: I. Introduction and relation to other models Physica 57 536\.
  * [35] U. Wolff 1990 Asymptotic freedom and mass generation in the O(3) nonlinear $\sigma$-model Nucl. Phys. B 334 581-610.
  * [36] H. G. Evertz, G. Lana, and M. Marcu 1993 Cluster algorithm for vertex models Phys. Rev. Lett. 70 875\.
  * [37] Y. Miyajima, Y. Murata, Y. Tanaka, and M. Mochizuki 2021 Machine-learning detection of the Berezinskii-Kosterlitz-Thouless transitions in the $q$-state clock models Phys. Rev. B 104 075114\.
  * [38] W. J. Hu, R. Singh, and R. Scalettar 2017 Discovering phases, phase transitions, and crossovers through unsupervised machine learning: A critical examination Phys. Rev. E 95 062122\.
  * [39] C. Wang and H. Zhai 2017 Machine learning of frustrated classical spin models. I. Principal component analysis Phys. Rev. E 96 144432\.
  * [40] K. Shiina, H. Mori, Y. Tomita, H. K. Lee, and Y. Okabe 2021 Inverse renormalization group based on image super-resolution using deep convolutional networks Sci. Rep. 11 9617\.
  * [41] Z.-Q. Li, L.-P. Yang, Z. Y. Xie, H.-H. Tu, H.-J. Liao, and T. Xiang 2020 Critical properties of the two-dimensional $q$-state clock model Phys. Rev. E 101 060105(R).
  * [42] S. Hong and D.-H. Kim 2020 Logarithmic finite-size scaling correction to the leading Fisher zeros in the $p$-state clock model: A higher-order tensor renormalization group study Phys. Rev. E 101 012124\.
  * [43] H. Ueda, K. Okunishi, K. Harada, R. Krčmár, A. Gendiar, S. Yunoki, and T. Nishino 2020 Finite-$m$ scaling analysis of Berezinskii-Kosterlitz-Thouless phase transitions and entanglement spectrum for the six-state clock model Phys. Rev. E 101 062111\.
  * [44] G. Li, K. H. Pai, and Z.-C. Gu 2022 Tensor-network renormalization approach to the $q$-state clock model Phys. Rev. Research 4 023159\.
  * [45] H. Chen, P. Hou, S. Fang, and Y. Deng 2022 Monte Carlo study of duality and the Berezinskii-Kosterlitz-Thouless phase transitions of the two-dimensional $q$-state clock model in flow representations Phys. Rev. E 106 024106\.
  * [46] K. Hukushima and K. Nemoto 1996 Exchange Monte Carlo Method and Application to Spin Glass Simulations J. Phys. Soc. Jpn. 65 1604\.
  * [47] R. H. Swendsen and J. S. Wang 1987 Nonuniversal critical dynamics in Monte Carlo simulations Phys. Rev. Lett. 58 86\.
  * [48] L.-R. Zhang, C. Ding, Y. Deng, and L. Zhang 2022 Surface criticality of the antiferromagnetic Potts model Phys. Rev. B 105 224415\.
  * [49] X. Zou, S. Liu, and W. Guo 2022 Surface critical properties of the three-dimensional clock model Phys. Rev. B 106 064420\.
  * [50] A. O. Sorokin 2019 Critical density of topological defects upon a continuous phase transition Ann. Phys. (Amsterdam) 411 167952\.


# Light scattering by a subwavelength plasmonic array: anisotropic model

Anton Nemykin, Leonid Frumin, David Shapiro
Institute of Automation and Electrometry
Russian Academy of Sciences, Siberian Branch
Novosibirsk 630090, 1 Koptjug Ave, Russia
e-mail<EMAIL_ADDRESS>
###### Abstract

We calculate the light transmission by a subwavelength plasmonic array using
the boundary element method for parallel cylinders with different cross-
sections: circular or elliptic with axes ratio 4:1. We demonstrate that the
plasmonic resonance is sharper for the case of ellipses with the large
horizontal axis. This structure is susceptible to the refractive index
variations of media since the high derivative of reflection and transmission
coefficients are near the angle of total internal reflection. To obtain an
approximate analytical expression, we use the model of a metallic layer. We
explore the ”sandwich” structure with an anisotropic film between two
dielectrics and demonstrate its quantitative agreement with numerical results.

##  1 Introduction

Plasmons are surface waves of conduction electrons inside the metallic film
near the border with a dielectric. The dispersion relation for plasmons
differs from the bulk plasma due to an interaction with the dielectric medium.
The plasmon electric field decays when one penetrates deeper both into the
metal and dielectric areas. At plasmon resonance, there is a notable increase
in the electric field, and intensity [1].

Applications of surface plasmons are diverse: from optical biosensors [2, 3,
4] to the acceleration of relativistic electrons [5] and space jet engines
[6]. Let us consider the lower dielectric half-space as the substrate with a
planar periodic structure on top. The refractive index of the upper half-space
strongly affects the angular and wavelength spectral characteristics of the
scattering layer due to crucial changes in layer effective permittivity. These
resonance responses underlie plasmon biosensors. The scattering on the layer
is also essential for research in silicon photonics, such as the study,
development, and manufacture of optical microcircuits in which photons
propagate instead of electrons [7]. The all-optical circuits can significantly
enhance the density of communication channels and lead to considerable energy
savings at once.

Recently the study of the plasmon-enhanced local field in the array of
parallel circular metallic cylinders is proposed to improve the sensitivity to
refractive index variations [8]. The matter is rapid changes in the Fresnel
coefficients near the angle of total internal reflection (TIR). Earlier, we
analyzed the grid of parallel metallic cylinders between two dielectrics; see
[9] and references therein. The boundary element method (BEM) idea is to
transform the Maxwell equations to the boundary integral ones applying the
Green theorem. In the case of a grid with parallel cylinders, the Floquet
theorem reduces the problem of infinite grating to that with one cylinder
within one elementary cell. To get analytical formulas, we exploited
simplified models. The sequence of nanowires can be approximately replaced by
a thin layer with the averaged permittivity [10]. In the present paper, we
treat cylinders of the elliptic cross-section with different orientations and
compare the results to more general models of an anisotropic layer.

##  2 FEF calculation

Figure 1: The periodic set of circular cylinders near the interface between
free space and the glass. The inset describes the shapes of cross-section
considered: circle (i), vertical ellipse (ii), horizontal (iii).

We carry out the numerical modeling by the BEM [9] for cylinders near the
interface between upper free half-space ($\varepsilon_{2}=1$) and glass lower
half-space $\varepsilon_{1}=2.25$. We analyze an infinite grating with period
$d$ consisting of gold ($\varepsilon_{3}=-23.6+i1.27$ at
$\lambda=0.7749~{}\mu$m [11]). We consider three shapes of cylinder cross-
sections: circle of radius $a=0.05~{}\mu$m and ellipses with axial ratio 1:4
and 4:1, as shown by an inset in fig. 1. The cross-section area is fixed, and
the gap dimension is $0.01~{}\mu$m for all the cases; meanwhile, the period is
different. It is $d=0.11~{}\mu$m for circular cylinders, $0.06~{}\mu$m for
vertical ellipse, and $0.21~{}\mu$m for horizontal. The distance between their
centers and glass half-space is fixed $0.11~{}\mu$m for all the cases. The
incidence field is $p$-wave near angle $\theta=\theta_{0}=41.81^{\circ}$,
where $\theta_{0}$ is the TIR angle. The grid sampling in calculations is
nonuniform near $\theta_{0}$ by a power law with exponents 3/2.

Figure 2: Field enhancement factor (a), absorption (b), and transmission (c)
as a function of the incidence angle $\theta$ (degrees): circles (short
dashes), vertical ellipses (long dashes), horizontal ellipses (solid line).

Formulas for the extinction and scattering are

$C_{ext}=-\frac{1}{I_{0}}\int(\mathbf{S}_{ext}\cdot\mathbf{e}_{r})\,dA,\quad
C_{s}=\frac{1}{I_{0}}\int(\mathbf{S}_{s}\cdot\mathbf{e}_{r})\,dA,$ (1)

where $\mathbf{S}_{s}=[\mathbf{E}_{z}\times\mathbf{H}_{s}^{*}]$,
$\mathbf{S}_{ext}=[\mathbf{E}_{0}\times\mathbf{H}_{s}^{*}+\mathbf{E}_{s}\times\mathbf{H}_{0}^{*}]$
are the Pointing vectors of scattered radiation and the energy flux of
interaction between scattered and incident fields, responsible for the
extinction. We take the integral over the surface of a cylinder $A$. The
absorption coefficient was calculated as the difference $C_{a}=C_{ext}-C_{s}$
[12, 13]. We see that FEF and absorption have close qualitative behavior. The
transmission coefficient is calculated in the far-field domain, neglecting all
the evanescent modes. The transmission has different angular dependence with a
sharp edge at the angle of TIR.

Figure 2 shows a comparison of the angular dependencies near the TIR angle of
the field enhancement factor (FEF) (a), absorption (b), and transmission
coefficients (c) for three versions of cross-sections. Figure 2(a)
demonstrates the $\mbox{FEF}=|E/E_{0}|^{2}$ in the middle between neighbor
cylinders as a function of the incidence angle. Here $E$ is the electric field
in the middle of the gap, and $E_{0}$ is the incident field. A plasmon
resonance induces a peak, exceptionally sharp for the horizontal ellipses. We
observe the highest peak since the horizontal ellipse has a minimum radius of
curvature in the gap.

As previously mentioned [8], the resonance of the slit is excited by the
refracted into the upper half-space (Fresnel) field. There is only an
evanescent wave when $\theta$ exceeds the TIR angle. As a result, the Fresnel
field achieves its maximum near the TIR angle; its $x$-component vanishes. The
absolute maximum of the intensity is created only by the $y$-component of the
Fresnel field. Component $E_{x}$, and not $E_{y}$, experiences resonance in
the gaps. A natural question arises: how is this possible if the Fresnel
$x$-component of the field vanishes. It turns out that the field scattered on
the cylinders has a significant $x$-component. Since the cylinders are not
limiting subwavelength ($kd=0.9,0.5,1.7$ for (i)–(iii) cases), then the
scattering differs from the dipole case. The phase difference between neighbor
wires originates the additional contribution to the $x$-component. We
emphasize that the gaps between the cylinders are subwavelength: $k\Delta<1$,
where $\Delta=0.01~{}\mu$m.

Absorption, fig. 2(b), correlate well with FEF in fig. 2(a). It cannot be
surprising since the most intensive electric field and the largest damping is
achieved precisely in plasmon resonance. The absorption characteristics, to a
certain extent, reproduce the angular dependence of FEF. As for the
transmittance presented in fig 2(c), it vanishes at a greater angle than the
TIR. The angular distribution of the horizontal ellipses has a shape similar
to a resonance curve. For other cross-sections, we see a substantial
broadening. The curve is wider for the vertical ellipses, as well as in FEF
and absorption.

##  3 Anisotropic layer

Figure 3: The anisotropic Fresnel’s model of a layer between two half-spaces.
Permittivity is $\varepsilon$ for the parallel electric field and
$\varepsilon^{\bot}$ for transverse field with respect to the boundary.

Helmholz equation gives dispersion relations in each part of the area, and
boundary conditions match tangential components of electric and magnetic
vectors $\mathbf{E},\mathbf{H}$, respectively:

$\displaystyle-E_{0}\cos\theta+E_{\mathrm{R}}\cos\theta$ $\displaystyle=$
$\displaystyle-\frac{D}{\varepsilon}\cos\psi-\frac{D_{\mathrm{R}}}{\varepsilon}e^{\delta}\cos\psi,$
$\displaystyle-\frac{D}{\varepsilon}e^{\delta}\cos\psi-\frac{D_{\mathrm{R}}}{\varepsilon}\cos\psi$
$\displaystyle=$ $\displaystyle-E_{\mathrm{T}}\cos\phi\;,$ (2)
$\displaystyle\varepsilon_{1}E_{0}\sin\theta+\varepsilon_{1}E_{\mathrm{R}}\sin\theta$
$\displaystyle=$ $\displaystyle D\sin\psi-D_{\mathrm{R}}e^{\delta}\sin\psi,$
$\displaystyle De^{\delta}\sin\psi-D_{\mathrm{R}}\sin\psi$ $\displaystyle=$
$\displaystyle\varepsilon_{2}E_{\mathrm{T}}\sin\phi\;.$ (3)

We consider the scattering of incident wave with electric field
$\mathbf{E}_{0}$ by an isotropic or anisotropic layer of thickness $h$.

For the isotropic material, we substitute cylinders with dielectric constant
$\varepsilon_{3}$ and surrounding free space with a dielectric or metallic
thin film $\varepsilon$. In anisotropic case the permittivity tensor has two
diagonal components
$\varepsilon_{3}^{xx}=\varepsilon,\varepsilon_{3}^{yy}=\varepsilon^{\bot}$.
The incidence angle is $\theta$, the inside layer refraction angle is $\psi$,
and the transmittance angle is $\phi$. Electric fields of reflected and
transmitted waves are $\mathbf{E}_{\mathrm{R}}$ and $\mathbf{E}_{\mathrm{T}}$,
respectively. Vectors of electric displacement $\mathbf{D}$ and
$\mathbf{D}_{\mathrm{R}}$ are orthogonal to their wavevectors inside the
layer.

Figure 3 illustrates the electric field notation: we calculate tangential
components by multiplying vector lengths by cosines and sines of corresponding
refraction angles. The imaginary part of cosine or sine turns to zero only in
non-absorptive dielectric domains at refraction angles up to TIR. In general,
they are complex and should be defined carefully. Furthermore, due to nonzero
layer thickness, we take the phase shift into account:

$\delta={i}\frac{2\pi}{\lambda}h\,\sqrt{\varepsilon-\varepsilon_{1}\frac{\varepsilon}{\varepsilon^{\perp}}\sin^{2}\theta}\;.$
(4)

The complex angle $\psi$ inside anisotropic layer satisfies relations

$\sin\psi=\frac{\sqrt{\varepsilon_{1}}\sin\theta}{\sqrt{\varepsilon+\varepsilon_{1}\left(1-\frac{\varepsilon}{\varepsilon^{\perp}}\right)\sin^{2}\theta}}\;,\quad\cos\psi=\frac{\sqrt{\varepsilon-\varepsilon_{1}\frac{\varepsilon}{\varepsilon^{\perp}}\sin^{2}\theta}}{\sqrt{\varepsilon+\varepsilon_{1}\left(1-\frac{\varepsilon}{\varepsilon^{\perp}}\right)\sin^{2}\theta}}\;,$
(5)

here the cut in complex plane of the square root function goes along the
negative real semi-axis. The transmit refraction angle is

$\sin\phi=\frac{\sqrt{\varepsilon_{1}}}{\sqrt{\varepsilon_{2}}}\,\sin\theta\;,\quad\cos\phi=\frac{\sqrt{\varepsilon_{2}-\varepsilon_{1}\sin^{2}\theta}}{\sqrt{\varepsilon_{2}}}\;.$
(6)

Eq. (3) and (3) are sufficient to find transmittance

$\displaystyle\mathrm{T}=4\,\frac{\textrm{Re}\,(\sqrt{\varepsilon_{2}}\cos\phi)}{\textrm{Re}\,(\sqrt{\varepsilon_{1}}\cos\theta)}$
$\displaystyle\times\left|\left(\frac{\sqrt{\varepsilon_{2}}}{\sqrt{\varepsilon_{1}}}+\frac{\cos\phi}{\cos\theta}\right)\cosh\delta-\left(\frac{\varepsilon}{\varepsilon_{1}}\frac{\sin\psi}{\sin\theta}\frac{\cos\phi}{\cos\psi}+\frac{\varepsilon_{2}}{\varepsilon}\frac{\sin\phi}{\sin\psi}\frac{\cos\psi}{\cos\theta}\right)\sinh\delta\right|^{-2}\;,$
(7)

and absolute reflectance

$\mathrm{R}=\left|\frac{\left(\frac{\sqrt{\varepsilon_{2}}}{\sqrt{\varepsilon_{1}}}-\frac{\cos\phi}{\cos\theta}\right)\cosh\delta-\left(\frac{\varepsilon}{\varepsilon_{1}}\frac{\sin\psi}{\sin\theta}\frac{\cos\phi}{\cos\psi}-\frac{\varepsilon_{2}}{\varepsilon}\frac{\sin\phi}{\sin\psi}\frac{\cos\psi}{\cos\theta}\right)\sinh\delta}{\left(\frac{\sqrt{\varepsilon_{2}}}{\sqrt{\varepsilon_{1}}}+\frac{\cos\phi}{\cos\theta}\right)\cosh\delta-\left(\frac{\varepsilon}{\varepsilon_{1}}\frac{\sin\psi}{\sin\theta}\frac{\cos\phi}{\cos\psi}+\frac{\varepsilon_{2}}{\varepsilon}\frac{\sin\phi}{\sin\psi}\frac{\cos\psi}{\cos\theta}\right)\sinh\delta}\right|^{2}\;.$
(8)

These coefficients mean the ratio of transmitted through the layer and
reflected energy flux to incident one. So the sum of transmittance
$\mathrm{T}$, reflectance $\mathrm{R}$ and absorption $\mathrm{A}$ add up to
one. Then to get the absorption coefficient we apply relation
$\mathrm{A}=1-\mathrm{T}-\mathrm{R}$. Here the absorption is specific, i.e.,
normalized by period $d$ and $\cos\theta$. If phase shift $\delta$ tends to
zero Eq. (7) and (8) become simpler and turn to Fresnel’s ones for two half-
spaces with common boundary plane. If anisotropic layer becomes isotropic
$\varepsilon=\varepsilon^{\bot}$, then the transmittance (7) coincides with
found one for magnetic field [14, 15]. Note that the formulas for isotropic
medium are applicable also for reflectometry sensor based on metallic-
dielectric structure [16, 17].

Figure 4: Absorption (a) and transmission (b) of circular cylinders as a
function of the incidence angle $\theta$ (degrees): BEM calculation (solid
line), isotropic layer (dashes), non-isotropic medium (circles).

Figure 5: The same as in fig. 4, but for elliptic cylinders with vertical
orientation of large axis.

Figure 6: The same as in fig. 4, but for elliptic cylinders with horizontal
orientation of large axis.

Let us plot the angular distribution within the isotropic and anisotropic
layer models and compare them with numerical results. We fit the model
parameters (thickness and permittivity) to minimize the offsets by the least-
squares method. The functional includes both the residues of transmission and
specific absorption coefficients with equal weights:

$F=\sum_{n=1}^{N}\left|T(\theta_{i})-T_{i}\right|^{2}+\sum_{n=1}^{N}\left|A(\theta_{i})-A_{i}\right|^{2}\;,$
(9)

here $T(\theta),A(\theta)$ are the model transmittance or specific absorption
depending on refraction angle in glass, $T_{i},A_{i}$ are the transmittance or
absorption calculated by BEM numerically at the incidence angle $\theta_{i}$
in glass. Number $N=201$ corresponds to numerical samples of incidence angle
glass in our calculations. The best fit parameters of the least-squares
procedure are gathered in table 1.

How accurate the fit for absorption and transmission is, one can see from fig.
4, 5, 6. The isotropic model describes the BEM curves qualitatively. When the
real part of permittivity occurs negative (the third line and third column of
the table), the fitting is good in transmission but loses its accuracy in
absorption. Otherwise, the curve has a huge mismatch in transmission and
better fitting in absorption, though the distortion function remains large.
For example, it is negative for the isotropic (scalar) model of vertical
ellipses, fig. 5(a). Meanwhile, the anisotropic model describes all the curves
quantitatively with high accuracy. It becomes evident from the bottom line in
table 1. The better fitting by the anisotropic model is at least due to more
fitting parameters (5 vs. 3).

The anisotropy turns out to be sufficiently high in this structure. Different
signs of the real part in the second and third lines from the bottom indicate
that the media behaves as hyperbolic metamaterials, i.e., a medium with the
hyperbolic dispersion relation [18]. One could explain this observation within
the concatenated capacitor model in the circuit diagram. For circular
cylinders (nanowires) at a small fill-fraction of the metal compared to
dielectric media in a unit cell, the Maxwell — Garnett approach is exploited
to calculate the dielectric tensor components for metamaterials [19].

Table 1: Resulting parameters in isotropic and anisotropic layer models. Parameter |  Circle |  Vertical ellipse |  Horizontal ellipse
---|---|---|---
Isotropic/anisotropic layer thickness, nm $(h)$ |  55/40 |  9.2/103 |  24/19
Isotropic permittivity $(\varepsilon)$ |  $20.7+0.7\rm{i}$ |  $-27.6+1.5\rm{i}$ |  $93.7+8.6{i}$
Longitudinal permittivity $(\varepsilon)$ |  $27.0+1.0{i}$ |  $-0.95+0.08{i}$ |  $\;117+19.1{i}$
Transversal permittivity $(\varepsilon^{\bot})$ |  $-2.55+0.0{i}$ |  $1.34+0.04{i}$ |  $-2.96-0.55{i}$
Isotropic/anisotropic layer distortion functional $(F\cdot 10^{3})$ |  65.3/4.7 |  705/0.67 |  149/11.5

##  4 Conclusions

To study the angular dependence of scattering parameters near the TIR angle,
we perform numerical calculations by the BEM. We consider circular and
elliptic cylinders with vertical and horizontal large-axis orientations and
treat the changes in absorption, transmittance, and the field enhancement
factor. We examine two possible directions of the large ellipse axis:
horizontal and vertical. The plasmonic resonance is sharper for horizontal
ellipses, where the curvature radius at the slit is the least.

We compare the calculated absorption and transmission coefficient with the
‘sandwich’ model. In the model, we place a medium with some effective average
permittivity and fixed thickness between two dielectric half-spaces instead of
the parallel wires. This model allows the analytical expressions of
transmission and absorption coefficients. The least-square minimization shows
that both models reproduce well the behavior of curves. Moreover, the curves
practically coincide with BEM computation when one chooses the model of the
anisotropic layer. Therefore, the anisotropic model is more appropriate to
describe artificial nanostructured media known as metamaterials.

## Acknowledgments

We are grateful to O. V. Belai for helpful discussions. This work is supported
by the Russian Foundation of Basic Research, grant # 20-02-00211.

## References

  * [1] Zayats, A.V.; Smolyaninov, I.I.; Maradudin, A.A.  Nano-optics of surface plasmon polaritons.  Physics Reports 2005, 408, 131–314.
  * [2] Gupta, B.D.; Kant, R.  Recent advances in surface plasmon resonance based fiber optic chemical and biosensors utilizing bulk and nanostructures.  Optics & Laser Technology 2018, 101, 144–161.
  * [3] Xu, Y.; Bai, P.; Zhou, X.; Akimov, Y.; Png, C.E.; Ang, L.K.; Knoll, W.; Wu, L.  Optical refractive index sensors with plasmonic and photonic structures: promising and inconvenient truth.  Advanced Optical Materials 2019, 7, 1801433.
  * [4] Zhao, Y.; Tong, R.j.; Xia, F.; Peng, Y.  Current status of optical fiber biosensor based on surface plasmon resonance.  Biosensors and Bioelectronics 2019, 142, 111505.
  * [5] Peralta, E.A.; Soong, K.; England, R.J.; Colby, E.R.; Wu, Z.; Montazeri, B.; McGuinness, C.; McNeur, J.; Leedle, K.J.; Walz, D.; others.  Demonstration of electron acceleration in a laser-driven dielectric microstructure.  Nature 2013, 503, 91–94.
  * [6] Rovey, J.L.; Friz, P.D.; Hu, C.; Glascock, M.S.; Yang, X.  Plasmonic force space propulsion.  Journal of Spacecraft and Rockets 2015, 52, 1163–1168.
  * [7] Bogaerts, W.; Chrostowski, L.  Silicon photonics circuit design: methods, tools and challenges.  Laser & Photonics Reviews 2018, 12, 1700237.
  * [8] Frumin, L.; Shapiro, D.  Sensitivity enhancement of plasmonic grating in the local field.  Optics Express 2020, 28, 26143–26150.
  * [9] Frumin, L.; Nemykin, A.; Perminov, S.; Shapiro, D.  Plasmons excited by an evanescent wave in a periodic array of nanowires.  Journal of Optics 2013, 15, 085002.
  * [10] Nemykin, A.; Perminov, S.; Frumin, L.; Shapiro, D.  Excitation of a plasmon resonance in metal cylinders by an evanescent wave.  Quantum Electronics 2015, 45, 240.
  * [11] Palik, E.D., Ed.  Handbook of Optical Constants of Solids. Vol 1,2.; Academic Press: London, 1998.
  * [12] Bohren, C.F.; Huffman, D.R., Eds.  Absorption and scattering of light by small particles; Wiley: Weinheim, 2004.
  * [13] Zymovetz, S.V.; Geshev, P.I.  Boundary Integral Equation Method for Analysis of Light Scattering by 2D Nanoparticles.  Technical Physics 2006, 51, 291–296.
  * [14] Kotkin, G.L.  On excitation of surface wave.  Phys. Metals and Metallography 1966, 21, 479–480.
  * [15] Raether, H.  Surface Plasmons on Smooth and Rough Surfaces and on Gratings; Springer: Berlin, 1988.
  * [16] Terent’ev, V.S.; Simonov, V.A.  Experimental method of fabrication of a matched metal–dielectric structure for a sensor based on the effect of frustrated total internal reflection.  Optoelectronics, Instrumentation and Data Processing 2015, 51, 609–616.
  * [17] Goldina, N.D.; Terent’ev, V.S.; Simonov, V.A.  Spectral properties of a metal–dielectric sensor structure.  Optics and Spectroscopy 2016, 120, 796–802.
  * [18] Poddubny, A.; Iorsh, I.; Belov, P.; Kivshar, Y.  Hyperbolic metamaterials.  Nature photonics 2013, 7, 948–957.
  * [19] Cortes, C.L.; Newman, W.; Molesky, S.; Jacob, Z.  Quantum nanophotonics using hyperbolic metamaterials.  Journal of Optics 2012, 14, 063001.