# Efficient calculation of three-dimensional tensor networks

Li-Ping Yang  Department of Physics, Chongqing University, Chongqing 401331,
China  Y. F. Fu  Department of Physics, Renmin University of China, Beijing
100872, China  Z. Y. Xie<EMAIL_ADDRESS>Department of Physics, Renmin
University of China, Beijing 100872, China  T. Xiang<EMAIL_ADDRESS>Institute of Physics, Chinese Academy of Sciences, Beijing, 100190, China
Beijing Academy of Quantum Information Sciences, Beijing 100193, China
University of Chinese Academy of Sciences, Department of Physics, Beijing
100190, China

###### Abstract

We have proposed an efficient algorithm to calculate physical quantities in
the translational invariant three-dimensional tensor networks, which is
particularly relevant to the study of the three-dimensional classical
statistical models and the (2+1)-dimensional quantum lattice models. In the
context of a classical model, we determine the partition function by solving
the dominant eigenvalue problem of the transfer matrix, whose left and right
dominant eigenvectors are represented by two projected entangled simplex
states. These two projected entangled simplex states are not Hermitian
conjugate to each other but are appropriately arranged so that their inner
product can be computed much more efficiently than in the usual prescription.
For the three-dimensional Ising model, the calculated internal energy and
spontaneous magnetization agree with the published results in the literature.
The possible improvement and extension to other models are also discussed.

##  I Introduction

The many-body problem is one of the central problems in physics, and
developing accurate and efficient numerical methods that can effectively
handle the exponential growth of the corresponding Hilbert space has always
been a great challenge, especially for quantum systems. Based on the idea of
renormalization group and the tensor-network representation of partition
functions and wave functions, tensor-network methods have evolved
progressively to be an important member of many-body computational methods in
recent years [1, 2, 3]. In fact, due to the absence of sign problem and the
ability to deal with two-dimensional systems, tensor networks are drawing
increasing attention, and have been applied successfully to strongly-
correlated electron systems [4, 5], frustrated spin systems [6, 7, 8, 9, 10],
statistical models [11, 12, 13], topological order [14, 15, 16, 17], quantum
field theory [18, 19, 20], machine learning [21, 22], and quantum circuit
simulation [23], etc. Among the various tensor-network methods, imaginary time
evolution is a highly efficient method to determine the ground state of low-
dimensional quantum systems [3, 24, 25].

A central task in the application of tensor network states to 2+1 dimensional
quantum lattice models is to contract a two-dimensional tensor network, which
has a double-layer structure and bond dimension $D^{2}$, with $D$ the maximal
virtual bond dimension of the tensor network representation of a quantum state
[3]. This is extremely costly, and thus although $D$ is expected to be larger
to produce a more accurate result, it is limited to about $13$ in practical
calculations [26, 9, 17]. There have been some efforts in recent years to
solve this problem, such as employing symmetries [27, 28], combining Monte
Carlo sampling [29], and the nested tensor network method [30]. However,
apparently the problem is not completely solved, and developing efficient
algorithms to resolve this issue is still an important topic for tensor-
network communities.

Meanwhile, it is known that, in the formalism of path integral, the
equilibrium quantum many-body problem in $d$ dimensions is similar to the
classical many-body problem in $d$+1 dimensions [31, 32, 33, 34, 35, 11]. Thus
inspired by the success in one- and two-dimensional quantum systems, there are
also some efforts to apply tensor-network methods to three-dimensional (3D)
classical models. Actually, both the coarse-graining tensor renormalization
group algorithms [11, 36, 37, 38] and the transfer-matrix-based tensor-network
state methods [39, 40, 41, 42, 43, 44, 45] have been applied to these models
or related previously. Although the transfer-matrix-based methods can work
excellently in two-dimensional networks, they seem not so efficient in three-
dimensions [46]. In the literature, in order to obtain reliable results, the
variational optimization procedure is indispensable, thus due to the high
computational cost, as will be discussed later, the bond dimension of the
involved tensor-network states and the accuracy are rather limited, and some
auxiliary techniques are usually employed further to analyze the data.
Especially, to our best knowledge, in practice, whether a simple yet efficient
algorithm analogous to the simple-update algorithm [47] in two-dimensional
quantum lattice models can be developed for 3D classical models is unclear up
to now.

In this paper, we are trying to address the above issues. To be specific,
following the usual prescription of the transfer-matrix-based methods, first
we express the partition function of a 3D classical model in terms of some
special two-dimensional transfer matrices, and reduce the problem to the
dominant eigenvalue problem of the matrices. Then we solve the dominant
eigenvalue problem by representing the dominant eigenvector as a special
tensor-network state, namely the projected entangled simplex state (PESS)
[26], which is efficiently determined through a power iteration procedure
analogous to imaginary time evolution, as done for two-dimensional classical
models similarly [48]. Furthermore, a simple nesting technique is proposed in
which the PESS representations of the left and right dominant eigenvectors are
designed appropriately so that their inner-product can be expressed as a
tensor network with bond dimension $D$ instead of $D^{2}$, and thus the
contraction can be carried out much more efficiently. Combining this nesting
technique with the corner transfer matrix renormalization group (CTMRG)
algorithm [49, 50, 4], we are able to push the bond dimension $D$ to 20 in
this work. For the 3D Ising model, it shows that even if the tensor network
state is renormalized by the so-called simple update technique [25, 47, 26]
after each evolution step, the obtained local quantities, such as energy
density and spontaneous magnetization, are in good consistent with the
previous studies, and the estimated critical temperature $T_{c}$ is about
$4.50984(2)$ which has only a relative deviation of about $10^{-4}$ from the
best Monte Carlo estimations.

The rest of the paper is organized as follows. In Sec. II, we introduce some
details of the algorithm employed in this paper, including the nesting
technique. The numerical results for statistical averages, such as energy
density $E$ and spontaneous magnetization $M$, as well as the convergence
analysis, are presented in Sec. III. In Sec. IV, we summarize our paper, and
discuss the possible improvement and promising extensions briefly.

##  II Method

###  II.1 Tensor-network representation of the partition function

For concreteness, hereinafter, we will focus on the 3D Ising model. The
partition function can be written as

$Z=\sum_{\\{s\\}}\prod_{\langle ij\rangle}e^{\beta s_{i}s_{j}},$ (1)

where $s_{i}=\pm 1$ is the spin variable located on the $i$th lattice site,
$\beta$ is the inverse temperature, $\langle ij\rangle$ means the product is
performed over all the nearest-neighboring bonds, and the summation is over
all the spin configurations.

Figure 1: Special transfer matrices in the expression of the partition
function of 3D Ising model. (Left) The green and blue cubes constitute
transfer matrices $T_{2}$ and $T_{1}$, as expressed in Eqs. (4) and (5).
(Right) The definition of local tensor $T$ located in each colored cube, as
expressed in Eq. (3). As mentioned in the main text, the vertical direction is
referred to as $z$ direction for convenience.

.

For later use, we regroup the product in Eq. (1) in the unit of cubes that are
the building block of a cubic lattice, i.e.,

$Z=\sum_{\\{s\\}}\prod_{\alpha}T^{(\alpha)},$ (2)

where $T^{(\alpha)}$ is a rank-8 tensor defined at the $\alpha$th cube. For
example, as shown in the right panel of Fig. 1, for a cube where the spin
variables residing on the eight vertices are denoted by $s_{1}$-$s_{8}$, the
local tensor $T$ can be expressed as the following product of twelve Boltzmann
weights corresponding to each edge of the cube, respectively,

$\displaystyle T_{s_{1}s_{2}s_{3}s_{4}s_{5}s_{6}s_{7}s_{8}}$ $\displaystyle=$
$\displaystyle\exp\left[\beta(s_{1}s_{2}+s_{2}s_{3}+s_{3}s_{4}+s_{4}s_{1}+s_{5}s_{6}+s_{6}s_{7}+\right.$
$\displaystyle\left.s_{7}s_{8}+s_{8}s_{5}+s_{1}s_{5}+s_{2}s_{6}+s_{3}s_{7}+s_{4}s_{8})\right].$
(3)

Manifestly, to make Eq. (2) and (3) consistent, $T$ should be defined only in
two kinds of inequivalent cubes, as denoted as green and blue, respectively,
in the left panel of Fig. 1. Thus the two kinds of cubes form an alternative
or staggered structure in the vertical direction. This is very similar to the
case of the imaginary time evolution in quantum lattice models, where the
Trotter-Suzuki decomposition of the evolution operator $e^{-\tau H}$ always
leads to an alternative structure in imaginary time $\tau$ direction. In the
following, we will use this special structure extensively. For convenience,
the vertical direction will be referred to as $z$ direction hereinafter.

###  II.2 Determination of the tensor-network representation of the dominant
eigenvector

Following the prescription of the transfer-matrix-based method, one needs to
express the partition function in terms of some transfer matrices. To this
end, firstly, we introduce two tensor-network operators

$\displaystyle T_{1}=\bigotimes_{\alpha\in b}T^{(\alpha)},\quad
T_{2}=\bigotimes_{\alpha\in g}T^{(\alpha)},$ (4)

where the direct products are performed over $T$s defined at blue cubes and
green cubes, respectively, as illustrated in Fig. 1. And then we can identify
the following equality:

$Z={\rm{Tr}}\left(T_{2}T_{1}\right)^{2n},$ (5)

where the length in the z direction is denoted as $2n$ for convenience. It is
worth noting that $T_{1}$ and $T_{2}$ should be understood as matrices by
grouping indices properly in Eqs. (4) and (5), in order to make the operations
for matrices therein meaningful.

Once Eq. (5) is established, the calculation of partition function is
immediately reduced to the dominant eigenvalue problem of the transfer matrix
$T_{2}T_{1}$, which can be solved by the power iteration method. In this
paper, we represent the corresponding dominant eigenvector $|\Psi\rangle$ as a
PESS form [26, 51], which can be written as

$|\Psi\rangle=\sum_{\\{s\\}}\mathrm{Tr}(...A^{(\alpha)}_{a_{\alpha}b_{\alpha}c_{\alpha}d_{\alpha}}B^{(\beta)}_{a_{\beta}b_{\beta}c_{\beta}d_{\beta}}P^{(i)}_{a_{i}b_{i}}[s_{i}]...)|...s_{i}...\rangle,$
(6)

as illustrated in Fig. 2. Here $\alpha$ and $\beta$ denote the coordinates of
two inequivalent squares, at the center of which a rank-4 simplex tensor $A$
or $B$ is introduced to characterize the four-spin entanglement in that
square. $i$ denotes the coordinates of lattice sites, where a rank-3
projection tensor $P$ is defined with two virtual indices labeled as $a,b...$
and a single physical index labeled as $s$. Every two virtual indices
associated with the same bond take the same values. $\mathrm{Tr}$ is over all
the repeated virtual indices and $\sum$ is over all the spin configurations
$\\{s\\}$. In this paper, we employed the $C_{4}$ rotational symmetry of $A$
and $B$ that corresponds to the symmetry of local tensor $T$ defined in Eq.
(3), and the translational symmetry of the 3D lattice, thus there are only one
independent $P$, $A$, and $B$ present in Eq. (6). More background about the
PESS wave function can be found in the Appendix.

Figure 2: Tensor network representation of the dominant eigenvector of the
transfer matrix $T_{2}T_{1}$ appeared in Eq. (5), as formulated in Eq. (6).
The black vertical lines denote the physical configurations $\\{s\\}$ appeared
in Eq. (6). In the right panel, the indices of the local tensors are
explicitly shown for clarity.

As done in time evolution in quantum lattice systems, using the alternative
structure of $T_{1}$ and $T_{2}$, we can take a similar simple update strategy
[25, 47, 26] to determine the variational parameters $A$, $B$, and $P$ in
$|\Psi\rangle$. To be specific, starting from a random state
$|\Psi_{0}\rangle$ which has the same structure as $|\Psi\rangle$, we apply
$T_{1}$ and $T_{2}$ alternatively to $|\Psi_{0}\rangle$ and update the
parameters accordingly by local decompositions of the related clusters after
each projection. This procedure is applied repeatedly until the convergence is
reached, and the obtained tensor network state provides an approximate
representation of $|\Psi\rangle$ [48]. Taking $A$ as an example, Fig. 3
illustrates how a single projection step is performed. For more details, we
suggest referring to Ref. [26].

Figure 3: A single projection step of $T$ on the PESS ansatz (take the update
of simplex tensor $A$ for example), as described in the main text. To
illustrate the operation, the indices of the local tensor $T$ are labeled
explicitly. (a) The indices $\\{s_{i}s_{j}s_{k}s_{l}\\}$ in $T$ and $A$ are
summed over. A cluster tensor $\bar{S}$ is obtained by further absorbing four
mean-field entanglements, i.e., the positive definite vector $\lambda_{B}$
obtained from the higher-order singular value decomposition (HOSVD) of the $B$
tensor. (b) The HOSVD of $\bar{S}$ is performed to obtain four unitary
matrices $U$s and the core tensor $\bar{A}$. (c) Truncate $\bar{A}$ to keep
bond dimension $D$, recover the original mean-field entanglement structure,
and obtain the updated local tensors $A$ and $P$. For more details, one can
refer to the elaboration in Ref. [26].

###  II.3 Inner product between the left and right dominant eigenvectors, and
its nesting structure

Having described how to obtain the dominant eigenvectors, in this subsection,
firstly, we briefly describe how the statistical averages are obtained in the
traditional reduced tensor network method (RTN) [3, 30]. As long as the
dominant eigenvector $|\Psi\rangle$ is obtained, we can use the fundamental
formula to calculate the statistical averages of local physical quantities.
For example, the average magnetization located at the $i$th site can be
determined by

$\displaystyle
M_{i}=\frac{\sum_{\\{s\\}}s_{i}\prod_{\alpha}T^{(\alpha)}}{\sum_{\\{s\\}}\prod_{\alpha}T^{(\alpha)}}=\frac{{\rm{Tr}}(T_{2}T_{1})^{n}s_{i}(T_{2}T_{1})^{n}}{{\rm{Tr}}(T_{2}T_{1})^{n}(T_{2}T_{1})^{n}},$
(7)

where we have used the transfer-matrix expression of $Z$ and assumed that
$s_{i}$ sits in the middle in the $z$-direction. In the thermodynamic limit,
$n\rightarrow\infty$, we reach

$M_{i}=\frac{\langle\Psi^{\prime}|s_{i}|\Psi\rangle}{\langle\Psi^{\prime}|\Psi\rangle},\quad
M=\frac{\sum_{i\in\alpha}M_{i}}{8},$ (8)

where $\langle\Psi^{\prime}|$ is the left dominant eigenvector of $T_{2}T_{1}$
and can be derived easily due to the symmetry between $T_{1}$ and $T_{2}$. The
spontaneous magnetization $M$ is obtained eventually as above by averaging
over the eight spins in the same cube due to the translational invariance of
the tensor networks. Internal energy $E$ of the bonds can be obtained
similarly. To evaluate the ratios like Eq. (8), we can follow the widely-used
impurity method as discussed in Refs. [47, 27, 52], and some details can be
found in the Appendix.

As indicated explicitly in Eq. (8), the statistical average of a local
physical quantity needs to contract a two-dimensional tensor network, which is
identical to the expectation value calculation for quantum lattice models in
essence. Suppose both $|\Psi\rangle$ and $\bra{\Psi^{\prime}}$ have bond
dimension $D$, then the generated two-dimensional tensor network has dimension
$D^{2}$, as shown in Fig. 4. When $D$ is small, this tensor network can be
contracted by the CTMRG algorithm [49, 50, 4] effectively. More introduction
of the CTMRG algorithm used in this work is provided in the Appendix.

Figure 4: Sketch of $\langle\Psi^{\prime}|\Psi\rangle$ appeared in Eq. (8).
Here $\langle\Psi^{\prime}|$ is the left dominant eigenvector of $T_{2}T_{1}$.
The relevance to $|\Psi\rangle$ and thus the structure of the local tensors
$T^{a}$ and $T^{b}$, come from the symmetry between $T_{1}$ and $T_{2}$, i.e.,
they are transformed to each other by a finite translation in the diagonal
direction.

The simple-update-like algorithm described in Sec. II.2 is highly efficient,
and the leading computational cost scales as $D^{5}$ in our case. Thus PESS
eigenvector with large bond dimensions can be obtained without too much cost.
However, contracting a two-dimensional tensor network with bond dimension
$D^{2}$, as shown in Fig. 4, has a leading cost $D^{12}$ empirically, as
discussed in Ref. [30], and is extremely costly when $D$ is large. A possible
solution is the nested tensor network technique, which can reduce the cost to
$D^{9}$ [30]. Nevertheless, in this paper, we did not take this strategy.
Instead, we proposed a much simpler nesting technique by taking advantage of
the special cubic lattice structure. This nesting technique can reduce the
computational cost similarly as achieved by the nested tensor network method
(also scales as $D^{9}$), keeps the symmetry of the ground state properly, and
is especially suitable for cubic systems.

To see how it works, we first divide the cubic lattice into three parts,
namely two bulks whose contribution to the full partition function can be
expressed in terms of transfer matrices, and the surface part which connects
the two bulks and combines as the whole cubic lattice. This means we have
rewritten the partition function in another manner, as illustrated in Fig.
5(a),

$Z={\rm{Tr}}(T_{4}T_{3})^{n}YVX(T_{2}T_{1})^{n},$ (9)

where $T_{1}$ and $T_{2}$ are the transfer matrices corresponding to the lower
bulk, $T_{3}$ and $T_{4}$ correspond to the upper bulk similarly, and $X$,
$Y$, $V$ are auxiliary matrices to describe the surface part. Here we
intentionally decompose the lower and upper parts in different manners, and
their relative locations projected in the $xy$ plane are shown in Fig. 5(b).
Mathematically,

$X=\prod_{\langle ij\rangle_{t}\not\in g}e^{\beta s_{i}s_{j}},\quad
Y=\prod_{\langle ij\rangle_{b}\not\in o}e^{\beta s_{i}s_{j}},\quad
V=\prod_{\langle ij\rangle_{bt}}e^{\beta s_{i}s_{j}},$ (10)

which are explicitly shown in Figs. 5(a), and 5(c)–5(e). Here $X$ and $Y$ can
be understood as onsite diagonal matrices, $\langle ij\rangle_{t}\not\in g$
means the product is over the nearest spin pairs between green cubes on the
top surface of the lower bulk, and $\langle ij\rangle_{b}\not\in o$ means the
product is over the nearest spin pairs between orange cubes at the bottom
surface of the upper bulk. $\langle ij\rangle_{bt}$ means $V$ is the product
of Boltzmann weights corresponding to all the bonds connecting the two bulks
in the $z$ direction. In Fig. 5, the dashed bonds corresponding to $X$, $Y$,
and $V$, are denoted as black, blue, and red, respectively. Note in Eq. (9),
the size in the $z$ direction is assumed as $2n$+1 for convenience.

Figure 5: Decomposition of the cubic lattice into three parts, as described in
the nesting technique in Sec. II. (a) The cubes with colors blue, green,
orange and grey, constitute the transfer matrices $T_{1}$, $T_{2}$, $T_{3}$,
and $T_{4}$, respectively. The dashed lines with color black, blue, and red,
constitute the transfer matrices $X$, $Y$, and $V$, respectively. The full
partition function can be represented in terms of these transfer matrices, as
expressed in Eq. (9). (b) The spatial distribution of the cubes in $xy$ plane,
for the sixteen cubes shown in (a). (c) Relative location of X (black dashed
lines) to $T_{1}$ (blue) and $T_{2}$ (green), and this explains how the
operation $X|\Psi\rangle$ is performed. (d) Relative location of Y (blue
dashed lines) to $T_{3}$ (orange) and $T_{4}$ (grey), and this explains how
the operation $\langle\Phi|Y$ is performed. (e) Relative location of V (on-
site red dots) to all four $T$s, and this explains how the operation
$\langle\tilde{\Phi}|V|\tilde{\Psi}\rangle$ is performed. Note that (b)–(e)
display the projection of (a) onto the horizontal plane, therefore the simplex
tensors appeared in the PESS ansatz Eq. (6) and updated in Fig. 3, reside
exactly in the center of the colored squares.

In the thermodynamic limit, the evaluation of Eq. (9) is again reduced to the
dominant eigenvalue problem as before, and the statistical averages can be
determined similarly. For example, the bond energy in the $z$ direction can be
obtained by

$E_{p}=\frac{\langle\Phi|YV^{\prime}X|\Psi\rangle}{\langle\Phi|YVX|\Psi\rangle},\quad
V^{\prime}=(-s_{a}s_{b}e^{\beta s_{a}s_{b}})\prod_{\langle ij\rangle_{bt}\neq
p}e^{\beta
s_{i}s_{j}},{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{}}$
(11)

where $\langle\Phi|$ is the left dominant eigenvector of $T_{4}T_{3}$, and
$V^{\prime}$ differs from $V$ by only a single vertical bond denoted as $p$
with $s_{a}$ and $s_{b}$ located at its ends. Finally we have

$E_{p}=\frac{\langle\tilde{\Phi}|V^{\prime}|\tilde{\Psi}\rangle}{\langle\tilde{\Phi}|V|\tilde{\Psi}\rangle},\quad\langle\tilde{\Phi}|\equiv\langle\Phi|Y,\quad|\tilde{\Psi}\rangle\equiv
X|\Psi\rangle,$ (12)

The key observation is that $\langle\tilde{\Phi}|$ and $|\tilde{\Psi}\rangle$
have a similar structure , but the parameters are distributed separately in
space, i.e., their simplex tensors are defined separately in space at
different sets of squares, which can be seen from Figs. 5(b)–5(d) and Figs.
6(a)–6(c) straightforwardly. And this leads to a great advantage, i.e., if the
wavefunctions $\langle\tilde{\Phi}|$ and $|\tilde{\Psi}\rangle$ are
approximated by PESS with bond dimension $D$, then the generated two-
dimensional tensor network also has bond dimension $D$, instead of $D^{2}$. To
be specific, let us take the denominator as an example. Exploring the symmetry
between $\langle\tilde{\Phi}|$ and $|\tilde{\Psi}\rangle$, the resulting
tensor network are composed of three parameters, i.e., the local tensors in
$|\tilde{\Psi}\rangle$ and $|\tilde{\Phi}\rangle$, and a new local tensor
$T^{p}$ derived from $P$ therein, i.e.,

$\displaystyle
T^{p}_{ijkl}=\sum_{s_{a}s_{b}}P_{lj}[s_{a}]P_{ki}[s_{b}]e^{\beta s_{a}s_{b}},$
(13)

as illustrated in Figs. 6(c)–6(d). Therefore, compared with the reduced
method, which generates a tensor network with a squared bond dimension, this
nesting technique can greatly reduce the cost and leads to much more efficient
contraction, as achieved in Ref. [30] similarly. Considering the fact that the
environment dimension $\chi$ used in the CTMRG iterations has the order
$D^{2}$ empirically, the nesting technique can reduce the computational cost
from $D^{12}$ to $D^{9}$. In this paper, we push $D$ to 20 with the help of
this technique.

Figure 6: Illustration of the nesting structure of
$\langle\tilde{\Phi}|V|\tilde{\Psi}\rangle$ appeared in Eq. (12). The spatial
location corresponds to Fig. 5 exactly. (a) Structure of the ket state
$|\tilde{\Psi}\rangle$, whose virtual bonds are indicated as solid black
lines. (b) Structure of the bra state $\bra{\tilde{\Phi}}$, whose virtual
bonds are indicated as dashed green lines. (c) Structure of the inner product
$\langle\tilde{\Phi}|V|\tilde{\Psi}\rangle$, where new rank-4 tensor $T^{p}$
is defined in Eq. (13) and illustrated in (d). Here we have used the symmetry
between the two wavefunctions, and both definitions of $T^{p}$ hold, although
it is unnecessary for the nesting method. The local tensors are the parameters
in $|\tilde{\Psi}\rangle$ in this context, not to be confused with parameters
in $|\Psi\rangle$ in Fig. 2. Clearly, the bond dimension of the resulting
tensor network is not squared, which is different from the case in Fig. 4.

##  III Results

In this paper, we focus on the 3D Ising model, which is of long-standing
interest in statistical physics and condensed matter physics. Although there
is no analytical solution as in two-dimensional case, the higher-order tensor
renormalization group (HOTRG) [11, 53] and Monte Carlo methods [54, 55, 56,
57] have provided very accurate numerical data, which confirms a second-order
finite temperature phase transition. Therefore, this model provides a suitable
touchstone to test new numerical algorithms in higher dimensions.

In this paper, as described in Sec. II, we use the developed evolution method
combined with the simple update technique to determine the tensor-network
representation of the dominant eigenvectors of the transfer matrices.
Furthermore, to push $D$ to a larger value (up to 20 in this paper), in the
expectation value calculations, we always employ the nesting technique
described in Sec. II.3 to improve the efficiency of CTMRG algorithm. Note that
when $D$ is small, and the conventional reduced method can be efficiently
performed, the results obtained by the CTMRG with and without nesting
technique are consistent, as illustrated and discussed in the Appendix.

The result of internal energy $E$ is shown in Fig. 7. The reference curve
denotes high-precision Monte Carlo data obtained in Ref. [55]. It shows that
the proposed method can give consistent results at the off-critical region,
while near the critical point, enlarging the bond dimension $D$ can produce
more accurate results as expected. The spontaneous magnetization $M$ shows
similar behavior, as shown in Fig. 8, where the reference curve is obtained
from Monte Carlo [54]. For both quantities, our result coincides well with
previous studies, and the singular behavior can be seen clearly.

Figure 7: Energy estimation obtained from the simple update and the nesting
method. The Monte Carlo data from Ref. [55] are plotted as a reference. Figure
8: Magnetization estimation obtained from the simple update and the nesting
method. The Monte Carlo data from Ref. [54] are plotted as a reference.

From $E$ and $M$, and the critical temperatures derived as shown in Fig. 9, it
seems that the proposed method tends to underestimate the critical temperature
a little, which is different from the coarse-graining tensor renormalization
group methods as shown, e.g., in Refs. [58, 11]. This is probably related to
the simple update technique used in the evolution process, since it
essentially provides a Bethe lattice approximation of the dominant eigenvector
[59, 60] and thus tends to be disordered at finite temperature, especially at
the temperature close to but lower than critical point where the approximation
cannot give a good estimate of the correlation length. When $D=20$, the
estimated critical temperature $T_{c}$ is about $4.50984(2)$, which has only
about $10^{-4}$ relative deviation from the best estimation $T_{c}\sim
4.51152$ [11, 53, 57]. The comparison with other variational calculations is
summarized in Table 1. Although obtaining a benchmark result is not our
motivation, it shows that our method performs rather well. Compared to the
most recent PEPS calculations [44, 45], though our result seems slightly less
accurate, our method does not need other auxiliary techniques, such as
ignoring the long tails of magnetization with careful fitting [44], scaling
hypothesis and data collapse [45], bootstrap exponent fixing [45], etc. The
accuracy of the critical temperature and the critical exponent can be further
improved by further refining the simple-updated wave functions through
variational calculations and performing these advanced techniques.

Figure 9: Critical temperature obtained from the magnetization $M$, with respect to the bond dimension $D$ of the tensor network state representation of the transfer matrices’ dominant eigenvectors. The critical temperature is estimated as the lowest temperature which has almost zero magnetization, as shown in Fig. 8. Method | $T_{c}$
---|---
Monte Carlo [57] | 4.51152326(11)
HOTRG [11, 53] | 4.51152469(1)
KWA-based TNS∗ [40] | 4.5788
Vertex-type TNS∗ [41] | 4.5392
TPVA∗ [40, 61] | 4.5704, 4.554
TNS∗ [41] | 4.504
Algebraic variation∗ [42] | 4.547
TNS data fitting∗ [44] | 4.5118057(41)
TNS data collapse∗ [45] | 4.5104, 4.51170
Simple update (this paper) | 4.50984(2)
Table 1: Comparison of the critical point $T_{c}$ for the 3D Ising model
obtained from different tensor-network-state-based methods. The Monte Carlo
and HOTRG results are included as references. Note that for the methods
denoted by $(^{\ast})$, variational calculations have been employed and the
bond dimensions are rather small ($D\sim 4$). Our results are derived from the
magnetization data with $D=20$ simple-updated wave function.

Furthermore, as shown in Fig. 10, the critical exponent can be obtained by
fitting the magnetization data near critical temperature according to the
scaling relation

$\displaystyle M\sim\left(1-T/T_{c}\right)^{\beta}.$ (14)

We find that the exponent $\beta=0.335(16)$, close to the Monte Carlo [62]
(0.3262), conformal bootstrap [63] (0.326419), and HOTRG [11] (0.3295)
results.

Figure 10: Critical exponent $\beta$ for the magnetization obtained from
$D=20$ simple-updated wave function and $\chi=120$ nesting CTMRG method.
$\beta$ is estimated as 0.335(16), and the data obtained at temperature
ranging from 4.5050 to 4.5096 are used for fitting.

As to the nesting technique proposed in this paper, its convergence with
respect to the bond dimension $\chi$ of the environment tensors used in CTMRG
is shown in Fig. 11, for $D=20$. It shows clearly that for all the three
typical temperatures ranging from symmetry-breaking phase to paramagnetic
phase, the convergence is at least equally satisfying, compared with the
nested tensor network method proposed in Ref. [30]. Particularly, even at the
temperature that is very close to the critical value, the convergence is
already very acceptable when $\chi=180$ for $D=20$, with an error about
$10^{-6}$ for $E$ and $10^{-4}$ for $M$. This is a very nice feature for
tensor-network methods, especially for large $D$, as discussed in Ref. [30] in
detail.

Figure 11: Convergence analysis about energy and magnetization for $D=20$.
Three typical temperatures are chosen, i.e., temperature in symmetry-breaking
phase (4.45), close to critical point (4.508), and in disordered phase (4.6).

##  IV Summary

In this paper, we propose an efficient numerical method to contract the 3D
tensor networks with translational invariance. The result of the contraction
is expressed in terms of some transfer matrices, whose dominant eigenvectors
are represented as the PESS form and determined by power iterations analogous
to the imaginary time evolution algorithm [24, 48]. Especially, the PESS
representations of the left and right dominant eigenvectors are designed
appropriately so that their innerproduct has a nesting structure and can be
expressed as a tensor network with bond dimension $D$ instead of $D^{2}$, and
thus the contraction can be carried out much more efficiently. As to the 3D
Ising model, it can give very consistent results for both energy and
magnetization with previous studies, even though only the simple update
strategy is employed to update the tensor network states. The convergency of
the nesting technique is shown to be at least equally good to that of the
nested tensor network method [30]. When $D=20$, we obtain a critical
temperature about 4.50984(2), which is very close to the best-known
estimation. As far as we know, this is probably the successful trial of
applying the imaginary-time-evlution-like method to 3D classical models
without variational update procedure, thus this extends the application scope
of tensor network states. The application to other interesting but unsolvable
classical models in three dimensions, such as the Potts model [53], clock
model [64], dimer model [44], and lattice gauge models [19], is
straightforward.

As mentioned in Sec. III, the simple update strategy employed in this paper is
probably a reason why the obtained $T_{c}$ tends to be underestimated. Besides
the scaling hypothesis and data collapse technique [45], this tendency could
be eased in two different manners. One is to resort to the more involved but
more accurate update methods, such as the cluster update [65] and the full
update strategies [66, 67], which consider the renormalization effect of the
environment better. The other possible way is to utilize more sophisticated
ansatz and evolution techniques, such as PESS ansatz with stronger simplex
entanglement [26], the recently proposed regularized scheme of time evolution
[68], and the so-called minimal canonical form of tensor network states [69],
all of which are also expected to produce more accurate representation of the
dominant eigenvectors. In both cases, the nesting technique proposed in this
paper can still be applied without changing too much, and the improvement in
performance can be expected.

At last, it is worth mentioning that the proposed nesting technique in this
paper can also be extended to the two-dimensional quantum lattice models. In
order to achieve this, one needs to represent the partition function in a
different manner, e.g., similarly as we have done in Eq. (9), so that two
representations of the ground state with separated distributions in space can
be used simultaneously, e.g., as illustrated in Fig. 5(b) and expressed in Eq.
(12). This can be advantageous in the study of models which have square
structure and separatable block interactions, such as Shastry-Sutherland model
[70, 6], checkboard system [71], and so on. Exploring its potential and limits
is an interesting and promising project, and we would like to leave it as a
future pursuit.

## Acknowledgement

We thank Zhong-Cao Wei and Zhi-Yuan Liu for their contribution in the early
stage of this work. We are supported by the National R$\&$D Program of China
(Grants No. 2017YFA0302900 and No. 2016YFA0300503), the National Natural
Science Foundation of China (Grants No. 12274458, No. 11888101, No. 11874095,
and No. 11774420), and by the Research Funds of Renmin University of China
(Grant No. 20XNLG19). L.-P.Y. and Y.F.F contribute equally to this work.

## Appendix:FURTHER BACKGROUND AND COMPARISONS

In this Appendix, some further background on the concepts and techniques
involved in this paper, and some relevant comparisons are provided.

### 1\. PESS wave function ansatz

The PESS ansatz is particularly important in the nesting technique developed
in this paper. It was proposed in Ref. [26], and can be regarded as a
generalization of the projected entangled pair state (PEPS) ansatz proposed in
Ref. [1]. It extends pair entanglement well characterized by PEPS to simplex
entanglement that might be more important in some situations, and has been
successfully applied to a series of frustrated spin systems, such as kagome
spin liquid systems [8, 72], chiral spin liquid systems [17, 73], quantum
spin-orbital liquid systems [9, 74], etc.

As a wavefunction ansatz, the defining feature of a PESS is the existence of
simplex tensors in it, each of which connects more than two physical degrees
of freedom and captures the many-body entanglement among them. For example, in
Eq.(6) and Fig. 2, simplex tensors $A^{\alpha}$ and $B^{\beta}$ are introduced
to characterize the entanglement among the four spins defined on the vertices
of the $\alpha$th and $\beta$th squares respectively. For any given base ket
$|\\{s\\}\rangle$, the corresponding coefficient in the superposition is given
by the $\mathrm{Tr}$ operation, which means the contraction of a two-
dimensional tensor network and can be accomplished by, e.g., the CTMRG method.
For general discussions about tensor-network state, Refs. [1, 2, 3] are
excellent references covering almost all important topics.

In this paper, the translational symmetry and the $C_{4}$ rotational symmetry
of the simplex tensors are employed, thus we have only one independent $A$,
$B$, and $P$ as variational parameters that need to be determined by the
simple-update-aided imaginary-time evolution method, as discussed in the main
text and Ref. [26].

### 2\. CTMRG algorithm for arbitrary supercell

The CTMRG is a highly-efficient and extensively-used algorithm to contract
two-dimensional tensor networks, and in this paper it has been used to
calculate the expectation values, since the numerator and denominator in both
Eqs. (8) and (12), as well as other terms similar to the inner product of two
tensor-network states, can be represented as two-dimensional tensor networks
exactly. In this paper, for a given tensor-network state with bond dimension
$D$, the difference in the expectation value calculation between the
contraction of a normal network and a network with nesting structure, lies in
the fact that the normal one has bond dimension $D^{2}$ while the nested one
has dimension $D$. Therefore, though both the two kinds of tensor networks can
be contracted by CTMRG algorithm, the leading cost differs greatly, i.e.,
$D^{12}$ for the normal one and $D^{9}$ for the nested one (can be reduced to
$D^{10}$ and $D^{8}$ respectively by some advanced techniques, e.g., partial
SVD.)

In this paper, we used the CTMRG algorithm for arbitrary supercell that was
proposed in Ref. [4]. For the normal network depicted in Fig.4, the supercell
size is $2\times 2$ (one can always absorb the Q tensors into $T^{a}$ and/or
$T^{b}$), while for the nested network depicted in Fig. 6, the supercell size
is $4\times 4$. In the CTMRG algorithm, each independent tensor in the
supercell is associated with four corner tensors and four edge tensors, which
mimic its effective environment and are used to extract the expectation
values. For a supercell with size $m\times n$, we have stored $4mn$ corners
and $4mn$ edges in total, which are initialized as random numbers and are
updated iteratively until convergence is reached. The detailed iterative
relations among these corners and edges are determined according to their
relative spatial locations, and have been given clearly in Ref. [4]. For
general discussions about the CTMRG algorithm, Ref. [3] provide excellent and
friendly reviews. For the most recent development, Refs. [75, 76] provide
promising trials.

### 3\. Impurity method in tensor-network states

The impurity method is frequently used to calculate the expectation values
without resorting to the derivatives. It has been successfully applied early
in Refs. [47, 58], and was elaborated later in Refs. [27, 52]. In essence, it
aims to calculate the expectation values of some local operators, e.g., the
bond energy or the magnetization, which can be expressed as a ratio of two
quantities represented as tensor networks and differing from with each other
by only very few tensors (dubbed impurity tensors). For simplicity, let us
take the on-site magnetization operator $M_{i}$ as an example, and repeat Eq.
(8) in the following

$\displaystyle
M_{i}=\frac{\langle\Psi^{\prime}|s_{i}|\Psi\rangle}{\langle\Psi^{\prime}|\Psi\rangle}$
(15)

where $\langle\Psi^{\prime}|$ and $|\Psi\rangle$ are the left and right
dominant eigenvectors of $T_{2}T_{1}$ respectively. Suppose the two
eigenvectors have been represented as PESS form, as shown in Fig. 2, then the
numerator and the denominator can be represented as Figs. 12(a) and 12(b)
respectively. It shows that there is only one special impurity tensor
$Q^{\prime}$ in Fig. 12(a), which constitutes the only difference from Fig.
12(b). Thus $Q$ and $Q^{\prime}$ share exactly the same effective environment
$N^{(e)}$, and one can evaluate the ratio $M_{i}$ by determining $N^{(e)}$
through well-developed methods, e.g., the CTMRG method as described above, and
obtain the following simple ratio form of $M_{i}$,

$\displaystyle
M_{i}=\frac{\mathrm{Tr}\left(N^{e}Q^{\prime}\right)}{\mathrm{Tr}\left(N^{e}Q\right)}$
(16)

finally. For other local operators like bond energy, the process is very
similar but involves two impurity tensors. More details and examples can be
found in Ref. [52].

Figure 12: Sketch of the impurity method and the explanation of Eq. (15). (a)
The numerator $\langle\Psi^{\prime}|s_{i}|\Psi\rangle$ represented as a 2D
tensor network. (b) The denominator $\langle\Psi^{\prime}|\Psi\rangle$
represented as a 2D tensor network. Note that almost all the tensors in (a)
and (b) are the same, as defined in Fig. 4, and the only difference lies in
two tensors defined on the $i$th site, i.e., $Q^{\prime}$ (impurity tensor,
denoted as bright green) as defined in (c) and $Q$ as defined in (d). In the
definition of $Q^{\prime}$, the operator $s_{i}$ is included, and here its
matrix representation $\sigma^{z}$ is a diagonal matrix with nonzero elements
$\pm 1$.

### 4\. The result comparison between the CTMRG with and without nesting
technique

In this paper, the data shown in the figures in the main text are all obtained
by the CTMRG algorithms combined with the nesting technique. When $D$ is small
and the conventional reduced tensor network method, i.e., RTN, can be
efficiently performed, we can compare the results obtained by the two methods.
The results are shown in Figs. 13 and 14. In both figures, complete
convergence of the CTMRG calculations concerning the environment dimensions
$\chi$ has been reached for both methods. It shows that both the internal
energy $E$ and $M$ are very consistent for the two methods. Even though it
seems that near critical region, the energy result obtained by the nesting
technique is more closer to the Monte Carlo data, they are essentially the
same. Actually the tiny difference comes not from the two methods themselves
[i.e., whether use Fig. 4 or Fig. 6(c) to represent the overlap of two PESSs,
both of which are exact for two given representations], but from the fact that
the methods employ different decompositions of the partition function [i.e.,
Fig. 1 and Eq. (2) in RTN, while Fig. 5 and Eq. (9) in the nesting technique],
and that accordingly the two expectation value calculations involve different
tensor networks. This tiny difference is expected to vanish gradually as $D$
becomes larger, and should vanish in quantum lattice models where exact Fig. 5
and Eq. (9) are unnecessary when employing the nesting technique.

Figure 13: Comparison between the internal energy $E$ obtained from the
conventional RTN method and the nesting technique, for $D=10$. Note that the
PESS wave functions used in these two methods are the same, and in both
methods, the network is contracted by the CTMRG algorithm. Figure 14:
Comparison between the magnetization $M$ obtained from the conventional RTN
method and the nesting technique, for $D=10$. Note that the PESS wave
functions used in these two methods are the same, and in both methods, the
network is contracted by the CTMRG algorithm.

## References

  * [1] F. Verstraete and J. I. Cirac, arXiv:cond-mat/0407066 (2004); J. I. Cirac, D. Perez-Garcia, N. Schuch, and F. Verstraete, Rev. Mod. Phys. 93, 045003 (2021).
  * [2] S. Montangero, Introduction to Tensor Network Methods, Springer (2018).
  * [3] R. Orus, Ann. Phys. 349, 117 (2014); R. Orus, Eur. Phys. J. B 87, 280 (2014); R. Orus, Nat. Rev. Phys. 1, 538 (2019).
  * [4] P. Corboz, T. M. Rice, and M. Troyer, Phys. Rev. Lett. 113, 046402 (2014).
  * [5] Simons Collaboration on the Many-Electron Problem, Phys. Rev. X 5, 041041 (2015), and Phys. Rev. X 10, 031016 (2020); B. X. Zheng, et al., Science 358, 1155 (2017).
  * [6] P. Corboz and F. Mila, Phys. Rev. Lett. 112, 147203 (2014); P. Corboz and F. Mila, Phys. Rev. B 87, 115144 (2013).
  * [7] L. Wang, Z. C. Gu, F. Verstraete, and X. G. Wen, Phys. Rev. B 94, 075143 (2016).
  * [8] H. J. Liao, Z. Y. Xie, J. Chen, Z. Y. Liu, H. D. Xie, R. Z. Huang, B. Normand, and T. Xiang, Phys. Rev. Lett. 118, 137202 (2017)
  * [9] Q. Li, H. Li, J. Z. Zhao, H. G. Luo, and Z. Y. Xie, Phys. Rev. B 105, 184418 (2022).
  * [10] G. Ferrari, G. Magnifico, and S. Montangero, Phys. Rev. B105, 214201 (2022).
  * [11] Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, and T. Xiang, Phys. Rev. B 86, 045139 (2012).
  * [12] C. Wang, S. M. Qin, and H. J. Zhou, Phys. Rev. B 90, 174201 (2014).
  * [13] J. G. Liu, L. Wang, and P. Zhang, Phys. Rev. Lett. 126, 090506 (2021).
  * [14] C. Wille, O. Buerschaper, and J. Eisert, Phys. Rev. B 95, 245127 (2017).
  * [15] W. T. Xu, Q. Zhang, and G. M. Zhang, Phys. Rev. Lett. 124, 130603 (2020).
  * [16] D. J. Williamson, N. Bultinck, and F. Verstraete, arXiv:1711.07982.
  * [17] R. Wang, Z. Y. Xie, B. G. Wang, and T. Sedrakyan, Phys. Rev. B 106, L121117 (2022).
  * [18] F. Verstraete and J. I. Cirac, Phys. Rev. Lett. 104, 190405 (2010).
  * [19] Y. Z. Liu, Y. Meurice, M. P. Qin, J. Unmuth-Yockey, T. Xiang, Z. Y. Xie, J. F. Yu, and H. Y. Zou, Phys. Rev. D 88, 056005 (2013).
  * [20] A. Tilloy and J. I. Cirac, Phys. Rev. X 9, 021040 (2019).
  * [21] Z. Y. Han, J. Wang, H. Fan, L. Wang, and P. Zhang, Phys. Rev. X 8, 031012 (2018).
  * [22] Z. F. Gao, S. Cheng, R. Q. He, Z. Y. Xie, H. H. Zhao, Z. Y. Lu, and T. Xiang, Phys. Rev. Research 2, 023300 (2020).
  * [23] Y. Meurice, R. Sakai, and J. Unmuth-Yockey, Rev. Mod. Phys. 94, 025005 (2022).
  * [24] G. Vidal, Phys. Rev. Lett. 91, 147902 (2003), Phys. Rev. Lett. 93, 040502 (2004).
  * [25] G. Vidal, Phys. Rev. Lett. 98, 070201 (2007).
  * [26] Z. Y. Xie, J. Chen, J. F. Yu, X. Kong, B. Normand, and T. Xiang, Phys. Rev. X 4, 011025 (2014).
  * [27] H. H. Zhao, Z. Y. Xie, Q. N. Chen, Z. C. Wei, J. W. Cai, and T. Xiang, Phys. Rev. B 81, 174411 (2010).
  * [28] P. Schmoll, S. Singh, M. Rizzi, and R. Orus, Ann. Phys. 419, 168232 (2020); M. Mambrini, R. Orus, and D. Poilblanc, Phys. Rev. B 94, 205124 (2016).
  * [29] L. Wang, I. Pizorn, and F. Verstraete, Phys. Rev. B 83, 134421 (2011); W.-Y. Liu, S.-J. Dong, Y.-J. Han, G.-C. Guo, and L. X. He, Phys. Rev. B 95, 195154 (2017); H.-H. Zhao, K. Ido, S. Morita, and M. Imada, Phys. Rev. B 96, 085103 (2017).
  * [30] Z. Y. Xie, H. J. Liao, R. Z. Huang, H. D. Xie, J. Chen, Z. Y. Liu, and T. Xiang, Phys. Rev. B 96, 045128 (2017).
  * [31] S. Sachdev, Quantum Phase Transitions, 2nd Edition, Cambridge University Press (2011).
  * [32] M. Suzuki, Prog. Theo. Phys. 56, 1454 (1976).
  * [33] T. Nishino, J. Phys. Soc. Jpn. 64, 3598 (1995).
  * [34] R. J. Bursill, T. Xiang, and G. A. Gehring. J. Phys.: Conden. Matt., 8, L583 (1996).
  * [35] H. Niggemann, A. Klumper, and J. Zittartz, Z. Phys. B: Condens. Matter 104, 103 (1997).
  * [36] P. Y. Teng, Physica A 472, 117 (2017).
  * [37] D. Kadoh, and K. Nakayama, arXiv:1912.02414.
  * [38] D. Adachi, T. Okubo, and S. Todo, Phys. Rev. B 102, 054432 (2020).
  * [39] T. Nishino and K. Okunishi, J. Phys. Soc. Jpn. 67, 3066 (1998).
  * [40] K. Okunishi and T. Nishino, Prog. Theor. Phys. 103, 541 (2000); T. Nishino, K. Okunishi, Y. Hieida, N. Maeshima, and Y. Akutsu, Nucl. Phys. B 575, 504 (2000); A. Gendiar and T. Nishino, Phys. Rev. E 65, 046702 (2002).
  * [41] T. Nishino, Y. Hieida, K. Okunishi, N. Maeshima, Y. Akutsu, and A. Gendiar, Prog. Theor. Phys. 105, 409 (2001); A. Gendiar, N. Maeshima, and T. Nishino, Prog. Theor. Phys. 110, 691 (2003).
  * [42] S. G. Chung, Phys. Lett. A 359, 707 (2006).
  * [43] R. Orus, Phys. Rev. B 85, 205117 (2012).
  * [44] L. Vanderstraeten, B. Vanhecke, and F. Verstraete, Phys. Rev. E 98, 042145 (2018).
  * [45] B. Vanhecke, J. Hasik, F. Verstraete, and L. Vanderstraeten, Phys. Rev. Lett. 129, 200601 (2022).
  * [46] K. Okunishi, T. Noshinio, and H. Ueda, J. Phys. Soc. Jpn. 91, 062001 (2022).
  * [47] H. C. Jiang, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. 101, 090603 (2008).
  * [48] Note that, essentially the similar idea has been applied to two-dimensional classical models under the help of the canonical form of matrix product state. For example, see R. Orus, and G. Vidal, Phys. Rev. B 78, 155117 (2008).
  * [49] T. Nishino and K. Okunishi, J. Phys. Soc. Jpn. 65, 891 (1996).
  * [50] R. Orús and G. Vidal, Phys. Rev. B 80, 094403 (2009).
  * [51] N. Schuch, D. Poilblanc, J. I. Cirac, and D. P. Garcia, Phys. Rev. B 86, 115108 (2012).
  * [52] H. H. Zhao, Z. Y. Xie, T. Xiang, and M. Imada, Phys. Rev. B 93, 125115 (2016).
  * [53] S. Wang, Z. Y. Xie, J. Chen, B. Normand, and T. Xiang, Chin. Phys. Lett. 31, 070503 (2014).
  * [54] A. L. Talapov and H. W. J. Blote, J. Phys. A: Math. Gen. 29, 5727 (1996).
  * [55] Y. J. Deng and H. W. J. Blote, Phys. Rev. E 68, 036125 (2003); X. M. Feng and H. W. J. Blote, Phys. Rev. E 81, 031103 (2010).
  * [56] M. Hasenbusch, Phys. Rev. B 82, 174433 (2010).
  * [57] A. M. Ferrenberg, J. Xu, and D. P. Landau, Phys. Rev. E 97, 043301 (2018).
  * [58] Z. Y. Xie, H. C. Jiang, Q. N. Chen, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. 103, 160601 (2009).
  * [59] W. Li, J. von Delft, and T. Xiang, Phys. Rev. B, 86, 195137 (2012).
  * [60] H. J. Liao, Z. Y. Xie, J. Chen, X. J. Han, H. D. Xie, B. Normand, and T. Xiang, Phys. Rev. B, 93, 075154 (2016).
  * [61] A. Gendiar, and T. Nishino, Phys. Rev. B 71, 024404 (2005).
  * [62] M. Hasenbusch, Int. J. Mod. Phys. C 12, 911 (2001).
  * [63] F. Kos, D. Poland, D. S. Duffin, and A. Vichi, J. High Energ. Phys. 08, 036 (2016)
  * [64] H. Shao, W. Guo, and A. W. Sandvik, Phys. Rev. Lett. 124, 080602 (2020).
  * [65] L. Wang, and F. Verstraete, arXiv:1110.4362.
  * [66] J. Jordan, R. Orus, G. Vidal, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008);
  * [67] M. Lubasch, J. I. Cirac, and M. C. Banuls, Phys. Rev. B 90, 064425 (2014); H. N. Phien, J. A. Bengua, H. D. Tuan, P. Corboz, and R. Orus, Phys. Rev. B 92, 035142 (2015).
  * [68] L. X. Cen, arXiv:2208.03436.
  * [69] A. Acuaviva, V. Makam, H. Nieuwboer, D. P. Garcia, F. Sittner, M. Walter, and F. Witteveen, arXiv:2209.14358.
  * [70] B. S. Shastry and B. Sutherland, Physica 108B, 1069 (1981); J. Yang, A. W. Sandvik, and L. Wang, Phys. Rev. B 105, L060409 (2022).
  * [71] H. Y. Zou, F. Yang, and W. Ku, arXiv:2011.06520.
  * [72] C. Y. Lee, B. Normand, and Y. J. Kao, Phys. Rev. B 98, 224414 (2018).
  * [73] S. Niu, J. Hasik, J. Y. Chen, and D. Poilblanc, Phys. Rev. B 106, 245119 (2022).
  * [74] Y. Liu, Z. Y. Xie, H. G. Luo, and J. Z. Zhao, Phys. Rev. B 107, L041106 (2023).
  * [75] X. F. Liu, Y. F. Fu, W. Q. Yu, J. F. Yu, and Z. Y. Xie, Chin. Phys. Lett. 39, 067502 (2022).
  * [76] I. V. Lukin, and A. G. Sotnikov, Phys. Rev. B 107, 054424 (2023).


# Sensing Strain-induced Symmetry Breaking by Reflectance Anisotropy
Spectroscopy

M. Volpi<EMAIL_ADDRESS>ETH Zurich, Department of Materials,
Laboratory for Nanometallurgy, Vladimir-Prelog-Weg 5, Zurich, 8093,
Switzerland  S. Beck  ETH Zurich, Department of Materials, Materials Theory,
Wolfgang-Pauli-Str. 27, Zurich, 8093, Switzerland  Center for Computational
Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
A. Hampel  Center for Computational Quantum Physics, Flatiron Institute, 162
5th Avenue, New York, NY 10010, USA  H. Galinski  ETH Zurich, Department of
Materials, Laboratory for Nanometallurgy, Vladimir-Prelog-Weg 5, Zurich, 8093,
Switzerland  A. Sologubenko  ETH Zurich, Department of Materials, Laboratory
for Nanometallurgy, Vladimir-Prelog-Weg 5, Zurich, 8093, Switzerland  ETH
Zurich, Scientific Center for Optical and Electron Microscopy, ScopeM,
Vladimir-Prelog-Weg 5, Zürich, 8093, Switzerland  R. Spolenak  ETH Zurich,
Department of Materials, Laboratory for Nanometallurgy, Vladimir-Prelog-Weg 5,
Zurich, 8093, Switzerland

###### Abstract

Intentional breaking of the lattice symmetry in solids is a key concept to
alter the properties of materials by modifying their electronic band
structure. However, the correlation of strain-induced effects and breaking of
the lattice symmetry is often indirect, resorting to vibrational spectroscopic
techniques such as Raman scattering. Here, we demonstrate that reflectance
anisotropy spectroscopy (RAS), which directly depends on the complex
dielectric function, enables the direct observation of electronic band
structure modulation. Studying the strain-induced symmetry breaking in copper,
we show how uniaxial strain lifts the degeneracy of states in the proximity of
the both L and X symmetry points, thus altering the matrix element for
interband optical transitions, directly observable in RAS. We corroborate our
experimental results by analysing the strain-induced changes in the electronic
structure based on ab-initio density functional theory calculations. The
versatility to study breaking of the lattice symmetry by simple reflectance
measurements opens up the possibility to gain a direct insight on the band-
structure of other strain-engineered materials, such as graphene and two-
dimensional (2D) transition metal dichalcogenides (TMDCs).

Strain Engineering, Density functional Theory, Reflectance Anisotropy
Spectroscopy

††preprint: AIP/123-QED

Strain engineering, i.e. the controlled deformation of solid state materials,
is an elegant and versatile method to unlock new material properties. Thus,
strain-induced symmetry breaking has gained considerable attention in
condensed matter physics, enabling the observation of the anomalous Hall
effect (AHE) in magnetic thin films [1], breaking the absorption limit of
silicon [2] and band-gap engineering of halide perovskites [3]. Moreover, the
concept proves most helpful to tune the band structure of 2D materials [4, 5,
6], where strain engineering combined with the outstanding mechanical
flexibility of 2D materials enables tailorable enhanced electronic and optical
properties.

For instance, single layer MoS2 subjected to uniform uniaxial and biaxial
tensile strain shows its bandgap linearly decreasing with increasing tensile
strain, first undergoing a direct-to-indirect bandgap transition and then a
semiconductor-to-metal transition [7, 8]. Also, strain abruptly modifies the
electronic band structure of graphene, leading to the appearance of an optical
band gap [9] accompanied by a change in the optical conductivity [10].

However, the measurement of the strain state and the correlation of the
reported strain-induced effects with the breaking of the lattice symmetry is
often indirect. Commonly, Raman scattering is employed, which allows for an
indirect characterization of strain due to shifting vibrational modes [11]. A
determination of the single strain components is not straightforward and
requires additional information about the stress state.
In this Letter, we demonstrate that reflectance anisotropy spectroscopy (RAS),
thanks to its direct correlation to the dielectric function, enables the
direct observations of electronic band structure modulation due to strain-
induced symmetry breaking. We employ density functional theory (DFT)
calculations to explain and characterize the emergence of the RAS signal as
result of the symmetry breaking. RAS, also known as reflectance difference
spectroscopy (RDS) [12, 13, 14], is a non-destructive phase modulating
technique capable of detecting materials optical anisotropy (Fig. 1a). RAS
measures the complex reflection anisotropy $\Delta r/r$, defined as the
normalized difference of complex frequency dependent Fresnel reflection
amplitudes $r_{x}$ and $r_{y}$ of two orthogonal directions $x$ and $y$:

$\frac{\Delta r}{r}=\frac{r_{x}-r_{y}}{r_{x}+r_{y}}\quad.$ (1)

Figure 1: RAS setup and spectra. RAS is composed of a Xe-lamp, a collimator
and a polarizer along the incident beam direction that polarize the light beam
at 45° with respect to the x and y specimen axes. Along the reflected beam
direction, it features a photoelastic modulator and an analyser to modulate
the phase of the light, which is eventually measured at different wavelengths
by means of a monochromator and a photomultiplier. The annealed Cu
polycrystalline specimen is reported in an EBSD fashion. b) reports the
sample’s orientation map, where colours represent out-of-plane
crystallographic orientations as in the legend, and the relative pole figures
along the 3 main crystallographic axes. c) and d) are the real and imaginary
part, respectively, of the RAS signal for a set of uniaxial tensile and
compressive strains. The absence of signal at 0 % strain is due to an overall
in-plane optical isotropy. The appearance of the elastic feature for uniaxial
strains higher than 0 % – indicated by the black dashed line – is due to
symmetry break. The inset in d) shows the linear relation between the
$\Delta{r}/{r}$ at the elastic feature energy ( 4.2 eV) with the strain,
ensuring that yield has not been reached.

RAS, developed by Aspnes et al. [15], is a versatile method to study surface
reconstruction [16], dimer orientations [17], epitaxial growth [16], piezo-
optical properties [13], and thin film growth [14]. To date, it also extends
to surface physics and chemistry, enabling to resolve phenomena such as
surface molecular absorption [14].
Intuitively, cubic crystallographic orientations like [1 1 1] and [1 0 0] do
not contribute to $\Delta r/r$ as $r_{x}=r_{y}$ nor do polycrystalline
specimens with random grain orientations whose in-plane anisotropy averages
out [12]. Therefore, even though largely employed in industry, polycrystalline
specimens have been disregarded for long. The RA signal is often linked to
band structure modulation [14–17] and first ab-initio DFT simulations confirm
this [18], focusing on surface states in single-crystalline materials.
However, a comprehensive study linking the strain-induced reflection
anisotropy spectra to symmetry breaking in the band structure and a general
understanding of its physical origin is still lacking.

Here, we probe the strain-induced symmetry breaking using the specific example
of Cu thin films. 200 nm thick polycrystalline Cu films were deposited using
magnetron assisted physical vapour deposition (PVD) on strain-free polyimide
substrate (Kapton E® Du Pont) (see also supplementary material [19]). The
choice of pure Cu, with a plasma frequency at 7.2 eV, excludes absorptions
caused by non-metal-like behavior within the explored energy range (1.5-5.5
eV) [20], and facilitates an undisturbed study of interband transitions. The
same accounts for free electron contributions (Drude model) which only play a
role for energies below 1.5 eV and hence can be neglected. Moreover, these
films are easy to produce and simpler to strain than e.g. 2D materials [11];
also, they guarantee a high surface quality and flatness that do not cause
topography related contributions to the RA signal [Fig. S1]. Prior to
analysis, thin films were annealed in vacuum at 300°C for 24 h. This thermal
treatment ensures the desired (1 1 1) out-of-plane texture shown in Figure 1b,
which prevents parasitic optical contribution from inherent anisotropic
surface energy states [12].

Figure 2: DFT simulations of RA signal. The figure reports the band structure
obtained by DFT for (a) the unstrained Cu unit cell (space group #225) and (b)
the elongated cell along $b$ [2 -1 -1] (space group #12). A schematics of a 2D
representation (1 1 1) cell is represented next to the band structure. (c) and
(d) display the real and imaginary components, respectively, of the simulated
RA signal (solid lines) superimposed to the experimental (dashed lines) for
both pristine and 0.8 % elongated unit cell and thin film, respectively. The
experimental signal is offset with respect to the $y$ axis in both (c) and (d)
to facilitate the comparison with simulations. The black dashed lines
highlight the features of interest at 2.9 and 4.2 eV.

The pole figures observed in Figure 1b, indicate an in-plane random
orientation. Thus, the unstrained specimens considered in this work do not
carry in-plane inherent asymmetry and $r_{x}$ and $r_{y}$ coincide. So, the
numerator in Eq. 1 vanishes and no RA signal is seen. This situation changes
once the films are uniaxially strained. Strain is induced on 3 x 25 mm2 Cu
films by using a tensile stage (Kammrath & Weiss GmbH, Germany). Considering
the enhanced mechanical properties of thin films’ fine structure defined by
the Hall-Petch theory [21], we ensure elastic behavior up to 0.8 % strain,
avoiding optical contributions from dislocations. Figure 1c and 1d display the
real and imaginary reflectance anisotropy spectra of a Cu thin film at varying
uniaxial strain within the elastic regime. The observed lineshape resembles a
complex Lorentzian oscillator centered at a resonant frequency (energy).
Indeed, similarly to the frequency dependence of the Lorentzian response [22],
the imaginary part Im($\Delta{r}/{r}$) shows a peak, while the real part
exhibits a typical 2${}^{\text{nd}}$ derivative trend, featuring a positive or
negative slope depending on the strain state. Figure 1c and 1d also show that
the feature disappears and reverses in sign the strain state changes from
tension to compression. In accordance with synchrotron measurements [23], the
inset in Fig. 1d reveals that the peak of Im($\Delta{r}/{r}$) scales linearly
with strain. Owing to the linear elastic response, we will refer to it as
”elastic feature” throughout the text.

Figure 3: Azimuthal RAS signal. The schematics (a) shows the experimental
approach. The polycrystalline Cu elongates along the strain direction and
contracts in the perpendicular one. This induces optical anisotropy which is
probed over a set of azimuthal angles $\vartheta$. The resulting imaginary
component of the experimental RAS signal upon 0.8 % strain over a set of
angles is reported in (b), e.g. 0∘ equals $2*(r_{x}-r_{y})/(r_{x}+r_{y})$, 15°
equals
$2*(r_{x+15^{\circ}}-r_{y+15^{\circ}})/(r_{x+15^{\circ}}+r_{y+15^{\circ}})$,
30° equals
$2*(r_{x+30^{\circ}}-r_{y+30^{\circ}})/(r_{x+30^{\circ}}+r_{y+30^{\circ}})$,
etc. The imaginary component of the simulated RAS signal upon 0.8 % elongation
in b [1 -1 0] over the same set of angles is reported in (c), e.g. 0∘ equals
$2*(r_{b}-r_{c})/(r_{b}+r_{c})$, 15° equals
$2*(r_{b+15^{\circ}}-r_{c+15^{\circ}})/(r_{b+15^{\circ}}+r_{c+15^{\circ}})$,
30∘ equals
$2*(r_{b+30^{\circ}}-r_{c+30^{\circ}})/(r_{b+30^{\circ}}+r_{c+30^{\circ}})$,
etc. The y-axis scale bars between b and c are adapted to the features’
amplitude.

To simulate the dielectric response of uniaxially strained Cu, we carry out a
set of ab-initio DFT calculations using the Quantum ESPRESSO software package
[24], and the exchange correlation functional according to Perdew, Burke, and
Ernzerhof [25]. We model three structures: two strained ones with the strain
applied perpendicular to the [1 1 1] direction, i.e. along [2 -1 -1] and [0 1
-1] (from now on referred to as $b$ and $c$ direction, respectively), to
replicate the experimental (1 1 1) texture, and the unstrained one for
comparison. Details on the construction of the supercells and computational
parameters are given in the supplementary information [19]. Tensile strain of
a magnitude of 0.8 % is applied in agreement to the experimental setup. It is
to note that uniaxial strain perpendicular to the [1 1 1]-direction reduces
the space group from #225 of the unstrained to #12 for the strained lattices
(Figure 2a-b), which results in a breaking of symmetry visible in the
bandstructure. To obtain the reflectance anisotropy, i.e. $\Delta{r}/{r}$, at
the (1 1 1) surface, the diagonal elements of the dielectric tensor
$\epsilon_{\alpha\alpha}$ in Cartesian coordinates are calculated with the
epsilon.x post-processing tool of Quantum ESPRESSO. We only consider interband
transitions in the calculation of $\epsilon_{\alpha\beta}(\omega)$, since the
Drude term is not relevant for the RAS signal in the energy range considered
in this work.

The contribution of the interband transitions to the total dielectric function
are directly related to the joint density of states, i.e. the sum over all
vertical transitions across the Fermi level across the Brillouin zone.
Consequently, its variation and in turn the RAS modulation stem from the
interband transitions between occupied and unoccupied energy levels.

To understand the origin of the calculated signal, we turn to the
bandstructure in Figure 2a, which shows the projection of the Kohn-Sham states
on atomic wavefunctions of bulk Cu. The path of high-symmetry points is chosen
such that it allows for a simple comparison with the strained structures. With
the occupied 3$d$ states in an energy range [-6,-1] eV, the main contribution
of the bands crossing the Fermi level is given by the 4$p$ states, followed by
the unoccupied 4$s$ states.

In the unstrained configuration, the high symmetry of the structure results in
an isotropic dielectric tensor, such that the calculated RA signal (Fig. 2c-d
blue solid lines) vanishes in agreement with experimental results (Fig. 2c-d
red dashed lines). Figure 2b shows the band structure of the strained
structure along the crystallographic $b$-axis. At first glance the overall
result seems unchanged in comparison to the unstrained structure, however, due
to the symmetry breaking, previously symmetry-related Kohn-Sham eigenvalues
are no longer degenerate along the chosen $k$-point path. This is emphasized
in the inset of both bandstructure plots, focusing on the eigenstates in an
energy range of [3.9-4.3] eV along L-W-L and M2-D$||$D2-A in unstrained and
strained structure, respectively. While the unoccupied 4$s$ states are
degenerate in the former case, the uniaxial strain leads to an increase
(decrease) in energy of the corresponding state at M2 (A), leading to an
energy difference of 200 meV.

Thus, when summing over all $k$-points in the calculation of the dielectric
tensor, the transitions do not cancel out exactly. The relevant transition
channel resulting in the signal are excitations from occupied 4$p$ to
unoccupied 4$s$ states, which we crosschecked by performing calculations
excluding these transitions in the calculation of the dielectric function (see
supplementary information and Fig. S2 for details [19]). If the channel is not
suppressed, it results in the observed oscillation in the real part and a dip
in the imaginary part of the RA signal in an energy range of [3.9-4.6] eV, as
shown in Fig. 2c-d (blue solid lines). The corresponding analysis for strain
along the $c$-direction is in agreement with these results, see Fig. S4 in
supplementary [19]. This reveals the exceptional sensitivity of the RAS to the
asymmetries of the energy bands, enabling to probe band energy shifts on the
meV scale.

A similar trend is also found at the high-symmetry point X in the unstrained
structure, and Y2 and L2 in the strained one. Here, the signal results from
transitions from the occupied 3$d$ to the unoccupied 4$p$ states. However,
since the energy gap is smaller the resulting signal in the RA spectrum
appears at smaller energies. Precisely, Figure 2c-d reports it in a range
[2.7-3.2] eV with a dip at 2.9 eV. The sign is opposite in this case, as the
band at Y2 is shifted down, while at L2 it is shifted up.

Our experimental spectra are in excellent agreement with data from literature
[26, 27, 20, 28], affirming the special role of interband absoprtion due to a
van Hove singularity [29] close to the L symmetry point. Also, another
transition at a joint density of states critical point in the vicinity of the
X appears in the measured RA signal. Quite interestingly, as reported in
Figure 2d, the latter exhibits sign inversion with respect to simulations.

To analyze the origin of the flipped response at 2.9 eV, we test the
robustness of the signal when the polarization direction is not aligned with
the strain direction. Figure 3 reports a comparison between experiments and
DFT of the ImRAS signal for a set of azimuthal angles $\theta$ between
$0^{\circ}$ and $90^{\circ}$ (schematics in Fig. 3a). The polycrystalline
specimens elongated in the experimental reference system $x$ direction (Fig.
3a) experience optical anisotropy which is detectable at different in-plane
orientations (Fig. 3b). This coincides with the trend of simulated ImRAS
signal for the unit cell elongated along the $b$-direction (Fig. 3c). The sign
agreement and inversion at 4.2 eV and 2.9 eV, respectively, is conserved over
the in-plane rotation. This consistency indicates that, while the experimental
band shift in the vicinity of the L symmetry point equals the simulated one,
the one at the X symmetry point is opposite. This suggests that, although the
same symmetries are broken in both cases, the experimental strain state
differs from the calculated one. This is intuitive considering that the L
symmetry point lies along the strained direction, whereas the X point does
not. Thus, the uniaxial strain applied in this work directly breaks the
symmetry in the vicinity of L, whereas only subsequent relaxation along one of
the other directions breaks the symmetry in the vicinity of X.

Because of this, we attribute the flip in sign of the 2.9 eV feature to the
non-ideal experimental conditions. Indeed, contrary to the simulated single
crystal, polycrystalline thin films are bound to a polymeric substrate [19],
which constrains the natural transversal relaxation that would arise as a
result of a Poisson ratio mismatch. This prevents the atoms distribution in
specimens to reach the equilibrium, leaving residual strains on the Cu film
that are not present in the calculations. Experimental confirmation using free
Cu single crystals is challenging considering that the elastic limit
drastically drops compared to polycrystals and a 0.8 % strain would exceed the
yield point, thus, leading to the appearance of plasticity deformation RAS
contributions in the same spectral region [30, 27]. Alternatively, follow-up
DFT calculations without cell relaxation according to the Poisson ratio could
shed light into the sign flip. However, a clear instruction of the particular
state of relaxation of the experimental film remains elusive to date.
Furthermore, we note that the RAS sensitivity to the matrix element for
interband optical transitions spans over 360° in-plane orientations, as
emphasized in Figure 3. Therefore, although the (1 1 1) crystallographic plane
exhibits a three-fold symmetry, the signal disappears at $45^{\circ}$ and it
changes in sign (for constant magnitude) at a $90^{\circ}$ rotation (visible
in Fig. 3b and Fig. 3c). This result not only indicates a 2-fold symmetry in
the system (supplementary III [19]) but it also suggests that one can in
principle determine the relative orientation of different crystallographic
domains within a given material. Therefore, DFT calculations for surface [18]
and bulk could be used to generate RA-spectra of any desired material,
orientation and strain direction, paving the way for new applications.

To conclude, we have demonstrated that reflectance anisotropy spectroscopy
(RAS) can serve as robust method to study symmetry breaking, as it establishes
a direct link between strain state and electronic band structure in strain-
engineered materials. Using the specific example of Cu we show that uniaxial
strain breaks the crystal symmetry and lifts the degeneracy of the Kohn-Sham
eigenvalues. Consequence of the lifted degenerancy is an anisotropic interband
absorption from 4$p$ to 4$s$ states and 3$d$ and 4$p$, located in the vicinity
of van Hove singularities at the X and L points, respectively. We conclude
that RAS can be used to probe and resolve these strain-induced modulations of
interband transitions down to the meV range, allowing for a large sensitivity
to the strain direction and sign. Since RAS is versatile and can be easily
adapted to study strain-engineered phenomena of other materials, this opens up
possibilities for bandstructure-tuning in 2D materials [31] and new insights
into Weyl semimetals [32] and strain-controlled photoluminescence of quantum
dots [33].

###### Acknowledgements.

M.V. and S.B. contributed equally to this work. The authors would like to
acknowledge the support from Scientific Center for Optical and Electron
Microscopy (ScopeM) of ETH Zurich. This study is part of the strategic focus
area advanced manufacturing project PREcision Additive Manufacturing of
Precious metals Alloys (PREAMPA). The Flatiron Institute is a division of the
Simons Foundation.

## References

  * Kim  _et al._ [2020] W. J. Kim, T. Oh, J. Song, E. K. Ko, Y. Li, J. Mun, B. Kim, J. Son, Z. Yang, Y. Kohama, M. Kim, B.-j. Yang, and T. W. Noh, Science Advances 6 (2020).
  * Katiyar  _et al._ [2020] A. K. Katiyar, K. Y. Thai, W. S. Yun, J. Lee, and J.-h. Ahn, Science Advances 6, 10.1126/sciadv.abb0576 (2020).
  * Chen  _et al._ [2020] Y. Chen, Y. Lei, Y. Li, Y. Yu, J. Cai, M.-h. Chiu, R. Rao, Y. Gu, C. Wang, W. Choi, H. Hu, C. Wang, Y. Li, J. Song, J. Zhang, B. Qi, M. Lin, Z. Zhang, A. E. Islam, and B. Maruyama, Nature 577, 209â€“215 (2020).
  * Novoselov  _et al._ [2014] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2014).
  * Bai  _et al._ [2020] Y. Bai, L. Zhou, J. Wang, W. Wu, L. J. Mcgilly, D. Halbertal, C. Fan, B. Lo, F. Liu, J. Ardelean, P. Rivera, N. R. Finney, X.-c. Yang, D. N. Basov, W. Yao, X. Xu, J. Hone, A. N. Pasupathy, and X. Zhu, Nature Materials 19, 1068â€“1073 (2020).
  * Hsu  _et al._ [2020] C. Hsu, M. L. Teague, J. Wang, and N. Yeh, Science Advances 6, 10.1126/sciadv.aat9488 (2020).
  * Scalise  _et al._ [2012] E. Scalise, M. Houssa, G. Pourtois, V. Afanas’ev, and A. Stesmans, Nano Research 5, 43 (2012).
  * Ghorbani-Asl  _et al._ [2013] M. Ghorbani-Asl, S. Borini, A. Kuc, and T. Heine, Physical Review B - Condensed Matter and Materials Physics 87, 1 (2013).
  * Guinea  _et al._ [2010] F. Guinea, M. I. Katsnelson, and A. K. Geim, Nature Physics 6, 30 (2010).
  * Pellegrino  _et al._ [2010] F. M. Pellegrino, G. G. Angilella, and R. Pucci, Physical Review B - Condensed Matter and Materials Physics 81, 1 (2010).
  * Dai  _et al._ [2019] Z. Dai, L. Liu, and Z. Zhang, Advanced Materials 31, 1 (2019).
  * Weightman  _et al._ [2005] P. Weightman, D. S. Martin, R. J. Cole, and T. Farrell, Reports on Progress in Physics 68, 1251 (2005).
  * Rönnow  _et al._ [1999] D. Rönnow, L. F. Lastras-Martinez, M. Cardona, and P. V. Santos, Journal of the Optical Society of America A 16, 568 (1999).
  * Martin and Weightman [2001] D. S. Martin and P. Weightman, Surface and Interface Analysis 31, 915 (2001).
  * Aspnes and Studna [1985] D. Aspnes and A. Studna, Physical review letters 54, 1956 (1985).
  * Aspnes  _et al._ [1988] D. Aspnes, J. Harbison, A. Studna, and L. Florez, J. Vac. Sci. Technol. A 6, 1327 (1988).
  * Power  _et al._ [1998] J. R. Power, P. Weightman, S. Bose, A. I. Shkrebtii, and R. D. Sole, Physical Review Letters 80, 3133 (1998).
  * Harl  _et al._ [2007] J. Harl, G. Kresse, L. D. Sun, M. Hohage, and P. Zeppenfeld, Physical Review B - Condensed Matter and Materials Physics 76, 20 (2007).
  * [19] See supplemental material [URL to be inserted by publisher] for technical details of the DFT calculations,etc. .
  * Ehrenreich and Philipp [1962] H. Ehrenreich and H. Philipp, Physical Review 128, 1622 (1962).
  * Hall [1951] E. O. Hall, Proceedings of the Physical Society. Section B 64, 747 (1951).
  * Hummel [2012] R. E. Hummel, Electronic properties of materials (Springer, New York, 2012) Chap. 11, pp. 238–243, 4th ed.
  * Wyss  _et al._ [2017] A. Wyss, A. S. Sologubenko, N. Mishra, P. A. Gruber, and R. Spolenak, Journal of Materials Science 52, 6741 (2017).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. De Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics Condensed Matter 21, 10.1088/0953-8984/21/39/395502 (2009).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
  * Cole  _et al._ [2003] R. J. Cole, S. Kheradmand, D. D. Higgins, F. Madani, B. F. Macdonald, V. Koutsos, and J. R. Blackford, Journal of Physics D: Applied Physics 36, L115 (2003).
  * Wyss  _et al._ [2015] A. Wyss, M. Schamel, A. S. Sologubenko, R. Denk, M. Hohage, P. Zeppenfeld, and R. Spolenak, Journal of Physics D: Applied Physics 48, 10.1088/0022-3727/48/41/415303 (2015).
  * B.R. Cooper; H. Ehrenreich; H.R. Philipp [1964] B.R. Cooper; H. Ehrenreich; H.R. Philipp, Physical Review 138, 494 (1964).
  * Segall [1962] B. Segall, Physical Review 125, 109 (1962).
  * Blackford  _et al._ [2005] J. R. Blackford, V. Koutsos, D. S. Martin, D. S. Roseburgh, F. Madani-Grasset, J. V. Schwerdtfeger, C. Kratiras, M. Zaiser, and R. J. Cole, Physica Status Solidi C: Conferences 2, 3997 (2005).
  * Peng  _et al._ [2020] Z. Peng, X. Chen, Y. Fan, D. J. Srolovitz, and D. Lei, Light: Science and Applications 9, 10.1038/s41377-020-00421-5 (2020).
  * Toudert and Serna [2017] J. Toudert and R. Serna, Opt. Mater. Express 7, 2299 (2017).
  * Shimomura and Kamiya [2015] K. Shimomura and I. Kamiya, Applied Physics Letters 106, 10.1063/1.4913443 (2015).