#  Temperature Controlled Open Quantum System Dynamics using Time-dependent
Variational Method

Mantas Jakučionis1, Darius Abramavičius1 1Institute of Chemical Physics,
Vilnius University, Sauletekio Ave. 9-III, LT-10222 Vilnius, Lithuania

###### Abstract

Dirac-Frenkel variational method with Davydov $\text{D}_{2}$ trial
wavefunction is extended by introducing a thermalization algorithm and applied
to simulate dynamics of a general open quantum system. The algorithm allows to
control temperature variations of a harmonic finite size bath, when in contact
with the quantum system. Thermalization of the bath vibrational modes is
realised via stochastic scatterings, implemented as a discrete-time Bernoulli
process with Poisson statistics. It controls bath temperature by steering
vibrational modes’ evolution towards their canonical thermal equilibrium.
Numerical analysis of the exciton relaxation dynamics in a small molecular
cluster reveals that thermalization additionally provides significant
calculation speed up due to reduced number of vibrational modes needed to
obtain the convergence.

##  I Introduction

Obtaining dynamics of open quantum systems, i.e., quantum systems that are
identified as separate from its environment, yet remain in thermal contact
with it, is one of the most general quantum mechanical problems. Its
applicability range from excited state relaxation in optical response (Dorfman
_et al._ , 2016; Mukamel, 1995), energy transport in molecular aggregates
(Jakučionis  _et al._ , 2018; Chorošajev  _et al._ , 2017; Schröter  _et al._
, 2015; Chorošajev  _et al._ , 2014; Valkunas  _et al._ , 2013; May and Kühn,
2011), photosynthetic complexes (Jang and Mennucci, 2018; Thyrhaug  _et al._ ,
2018; Malý  _et al._ , 2016; Chenu and Scholes, 2015) and others (Sjakste  _et
al._ , 2018; Flesch  _et al._ , 2008; Lombardo and Turiaci, 2013; Yu and
Zhang, 2008; Ruostekoski and Walls, 1998). Prevalent theoretical description
is given in terms of system- bath model in constant temperature bath
conditions (Breuer and Petruccione, 2007; Weiss, 2012), where the system
degrees of freedom are coupled to the bath-induced thermal fluctuations
representing the environment, e.g., phonons or vibrational motion of
surrounding molecules. Fluctuations are modeled by an infinite number of
quantum harmonic oscillators _constituting_ the quantum bath at thermal
equilibrium.

These conditions can be fulfilled using the reduced density matrix approach
(Valkunas  _et al._ , 2013; Mukamel, 1995). First perturbation order, with
respect to the system-bath coupling, leads to the reduced equations of motion
of the system-only variables, while the bath is averaged out. Then the system
variables indirectly depend on the bath degrees of freedom via fluctuation
correlation functions, which are well-behaved analytical functions. At the
second perturbation order (Breuer and Petruccione, 2007; Valkunas  _et al._ ,
2013), equations of motion are reminiscent of the Pauli master equation with
relaxation coefficients calculated with respect to the thermal equilibrium.
However, now the resulting equations can lead to unphysical results, e.g.,
negative probabilities (Montoya-Castillo  _et al._ , 2015). The more
complicated fourth order equations of motion include divergent parameters and
are often avoided (Jang  _et al._ , 2002). Non-perturbative, numerically exact
approach of hierarchical equations of motion for the exponential fluctuation
correlation functions is available to obtain exact dynamics (Tanimura, 1990,
2006; Xu and Yan, 2007), chain-mapping techniques together with the time-
dependent density matrix renormalization group are alternatively possible for
structured environments (Tamascelli  _et al._ , 2019; Prior  _et al._ , 2010).
However, computational costs limit these methods to models with just few
degrees of freedom. A well known method of stochastic Schrödinger equation
requires averaging over many entangled trajectories to obtain dynamics at
finite temperature (Abramavicius and Abramavicius, 2014; Appel and Di Ventra,
2009; Biele and D’Agosta, 2012; Diósi and Strunz, 1997; de Vega  _et al._ ,
2005; Link and Strunz, 2017). Its hierarchical realisation (Hartmann and
Strunz, 2017) improves convergence, meanwhile, thermofield dynamics approach
tries to directly compute thermally averaged dynamics by mapping the initial
thermal density matrix onto a fictitious bath vacuum state and then coupling
system to it (Reddy and Prasad, 2015; Ritschel  _et al._ , 2015; Borrelli and
Gelin, 2016; Chen and Zhao, 2017). Alternativelly, dissipative dynamics can be
obtained by straightforward addition of a linear friction coefficient to the
model Hamiltonian (Martinazzo  _et al._ , 2006), however, it only applies at
zero temperature. Yet, in all these cases, thermal state of the nearest
surrounding is not under control.

An important aspect of the bath, more explicitly, of the finite number of bath
oscillators, is its heat capacity. For a single quantum harmonic oscillator
the heat capacity in the limit of weak system-bath coupling is given by

$\displaystyle c\left(\beta^{-1}\right)=$
$\displaystyle(\beta\omega)^{2}\frac{\exp(\beta\omega)}{(\exp(\beta\omega)-1)^{2}}\
,$ (1)

here throughout the paper $\hbar=1$, $\beta=\left(k_{\text{B}}T\right)^{-1}$
is an inverse temperature, $\omega$ is the oscillator frequency. When the
system exchanges energy with a bath made of such oscillators, its temperature
may be affected. If the system-bath energy exchange is excessively large, the
thermal energy can accumulate in the bath oscillators and this will
effectively change thermostat temperature (Abramavicius  _et al._ , 2018). In
most cases, the bath heating effect is undesirable as, in the system-bath
models, the bath is generally supposed to represent a constant temperature
thermostat.

On the other hand, the bath heating effect could be related to the natural
phenomenon of molecular local heating (Chen and Sorbello, 1993; Ichikawa  _et
al._ , 2007), i.e., if a molecule quickly dissipates a large amount of thermal
energy to its environment, e.g., due to exciton-exciton annihilation (Gulbinas
_et al._ , 1996; Valkunas and Gulbinas, 1997; Gulbinas  _et al._ , 2006) or
ultrafast molecular internal-conversion (Balevičius  _et al._ , 2018;
Balevičius Jr  _et al._ , 2019), the local heating of the molecule nearest
surrounding takes place and the further cooling process, the quantum
thermalization (Choi  _et al._ , 2019; Scarani  _et al._ , 2002), becomes an
important ingredient to consider when describing the corresponding
experiments.

In this work, we introduce the thermalization algorithm to the time dependent
variational theory that allows explicit control over the bath temperature. By
varying the bath size and the thermalization rate, both the degree of bath
heating and the cooling time can be adjusted. These properties allow to mimic
realistic physical conditions, making presented approach superior to the
density operator based approaches, where the bath heating is excluded, and to
the explicit bath models, where the bath temperature is not controlled.

##  II Fluctuating exciton model

We consider a molecular aggregate made of $N$ coupled chromophores at specific
sites. In the simplest case, the sites represent distinct molecules that can
be electronically excited by, e.g., laser or sunlight irradiation in visible
spectral region. Vibrational normal modes of molecules and of the surrounding
medium will be treated as the baths of harmonic oscillators. Each chromophore
is directly affected only by its own intramolecular vibrations and of its
closest environment, therefore, a separate and uncorrelated (local) manifold
of vibrational modes $q=1,2,\ldots,Q$ is associated with each chromophore.
Such model is characterized by a Hamiltonian
$\hat{H}=\hat{H}_{\text{S}}+\hat{H}_{\text{B}}+\hat{H}_{\text{SB}}$, with
system, bath and system-bath coupling terms being

$\displaystyle\hat{H}_{\text{S}}=$
$\displaystyle\sum_{n}\varepsilon_{n}\hat{a}_{n}^{\dagger}\hat{a}_{n}+\sum_{n,m}^{n\neq
m}J_{nm}\hat{a}_{n}^{\dagger}\hat{a}_{m}\ ,$ (2)
$\displaystyle\hat{H}_{\text{B}}=$
$\displaystyle\sum_{n,q}\omega_{nq}\hat{b}_{nq}^{\dagger}\hat{b}_{nq}\ ,$ (3)
$\displaystyle\hat{H}_{\text{SB}}=$
$\displaystyle-\sum_{n}\hat{a}_{n}^{\dagger}\hat{a}_{n}\sum_{q}\omega_{nq}g_{nq}\left(\hat{b}_{nq}^{\dagger}+\hat{b}_{nq}\right)\
,$ (4)

Here $\varepsilon_{n}$ denotes $n$th chromophore electronic excitation energy
$J_{nm}$ is the resonant coupling between $n$th and $m$th chromophores, while
$\hat{a}_{n}^{\dagger}$ and $\hat{a}_{n}$ are the corresponding electronic
excitation creation and annihilation bosonic operators. Frequency of the $q$th
vibrational mode in the $n$th bath is $\omega_{nq}$, the electron-vibrational
coupling is characterized by $g_{nq}$, while $\hat{b}_{nq}^{\dagger}$ and
$\hat{b}_{nq}$ are the creation and annihilation bosonic operators of the
$q$th mode in the $n$th bath.

In the following we consider only single electronic excitation in the
aggregate. Time evolution of a non-equilibrium state is described by Davydov
$\text{D}_{2}$ wavefunction (Frenkel and Frenkel, 1936; Sun  _et al._ , 2010)

$|\Psi_{\text{D}_{2}}\rangle=\sum_{n}^{N}\alpha_{n}\left(t\right)\hat{a}_{n}^{\dagger}\hat{a}_{n}|0\rangle_{\text{el}}\times\prod_{m,q}^{N,Q}|\lambda_{mq}\left(t\right)\rangle\
,$ (5)

where $\alpha_{n}\left(t\right)$ is the electronic excitation amplitude,
$|0\rangle_{\text{el}}=\prod_{n}|0\rangle_{n}$ is the global ground state,
when all sites are in their electronic ground states.
$|\lambda_{mq}\left(t\right)\rangle$ is the coherent state of the $q$th mode
in the $m$th bath (Zhang  _et al._ , 1990; Kais and Levine, 1990). It is fully
described by the time-dependent complex displacements,
$\lambda_{mq}\left(t\right)$. The time dependent Dirac-Frenkel variational
method allows to obtain equations of motion for parameters $\alpha_{n}$,
$\lambda_{mq}$ (Chorošajev  _et al._ , 2014; Jakučionis  _et al._ , 2018; Chen
_et al._ , 2018; Yan  _et al._ , 2020)

$\displaystyle\frac{\text{d}\alpha_{n}\left(t\right)}{\text{dt}}=-\text{i}\alpha_{n}\left(t\right)\varepsilon_{n}-\text{i}\sum_{m}^{m\neq
n}\alpha_{m}\left(t\right)J_{nm}$
$\displaystyle\qquad+\text{i}\alpha_{n}\left(t\right)\sum_{q}\omega_{nq}\left(2g_{nq}-h_{q}\right)\text{Re}\left(\lambda_{nq}\left(t\right)\right)\
,$ (6)
$\displaystyle\frac{\text{d}\lambda_{mq}\left(t\right)}{\text{dt}}=-\text{i}\omega_{mq}\left(\lambda_{mq}\left(t\right)-h_{q}\left(t\right)\right)\
.$ (7)

Here
$h_{q}\left(t\right)=\sum_{i}^{N}g_{iq}\left|\alpha_{i}\left(t\right)\right|^{2}$
is the site population-weighted electron-vibrational coupling strength. First
line in Eq. (6) describes dynamics of an isolated system. Accordingly, the
first term on the right hand side of Eq. (7) describes isolated oscillators.
Other terms are due to the system-bath interaction.

Description of the model at the given temperature $T$ requires creation of
statistical ensemble. This is achieved by Monte Carlo sampling over a
statistical thermal ensemble, i. e., over initial bath oscillator
displacements $\lambda_{mq}\left(0\right)$, sampled from the Glauber-Sudarshan
probability distribution (Glauber, 1963)

$\mathcal{P}\left(\lambda_{mq}\right)=\mathcal{Z}^{-1}\exp\left(-\left|\lambda_{mq}\right|^{2}\left[\text{e}^{\beta\omega_{mq}}-1\right]\right)\
.$ (8)

The ensemble describes canonical statistics of quantum harmonic oscillators,
which applies to our model prior to external perturbations. The ensemble
averaged quantities will be denoted by
$\left\langle\cdots\right\rangle_{\text{th}}$. The ensemble of exciton
trajectories allows to describe irreversible excitation energy relaxation.
While the initial thermal state before excitation can be properly defined, the
bath accepts energy during exciton relaxation and the state of the bath
_after_ relaxation is away from equilibrium. Equations of motion guarantee
energy conservation, hence the combined system-bath cannot thermalize. In
order to thermalize the bath, we extend the original model by introducing the
secondary bath (we will refer to the local baths as the primary baths).
Effective heat capacity of the secondary bath is infinite, hence, the bath can
be characterized by a constant temperature, $T_{\infty}$. The secondary bath
will not be treated explicitly: modes of the primary baths interact with the
secondary bath via stochastic scattering events, or quantum jumps (Plenio and
Knight, 1998; Luoma  _et al._ , 2020), which affect the kinetic energy of
primary baths modes.

The scattering statistics follows Poisson distribution
$P_{mq}\left(\theta,\tau\right)=\frac{1}{\theta!}\left(\tau\nu_{mq}\right)^{\theta}\text{e}^{-\tau\nu_{mq}}$,
which defines the probability of observing $\theta$ scattering events per time
interval $\tau$ with individual event scattering rate $\nu_{mq}$. Poisson
statistics is obtained by simulating a discrete-time Bernoulli process
(Kampen, 2007; Bertsekas and Tsitsiklis, 2008) in a limit of $\tau\rightarrow
0$ and $\nu_{mq}\tau\ll 1$. This is realized in simulations by dividing the
total evolution time $t_{\text{total}}$ into equidistant length $\tau$
intervals. At the end of each interval, for each mode in the primary bath, we
flip a biased coin with probability $\nu_{mq}\tau$ of landing “heads”. If the
coin lands heads, we shift momentum of the mode
$p_{mq}\left(kt\right)=\sqrt{2}\text{Im}\lambda_{mq}\left(kt\right)$ to a
value drawn from the Glauber-Sudarshan distribution, see Eq. (8), while the
coordinate remains unchanged. Otherwise, if coin lands tails, no changes are
done. To obtain converged statistics, we apply the thermalization algorithm to
every trajectory of the thermal ensemble.

##  III Simulation results

We first demonstrate control of the primary bath temperature of the simplest
possible system, a single, $N=1$, chromophore unit. For demonstration we set
up artificial conditions: the initial primary bath temperature
$T_{1}\left(0\right)=300\ \text{K}$, the secondary bath is at $T_{\infty}=200\
\text{K}$. The primary bath consists of $Q=750$ vibrational modes with
frequencies, $\omega_{q}=\omega_{0}+\left(q-1\right)\Delta\omega$. An offset
by $\omega_{0}=0.01\ \text{cm}^{-1}$ is introduced for stability and a step
size $\Delta\omega=1\ \text{cm}^{-1}$. The coupling parameters $g_{nq}$ follow
the super-Ohmic spectral density function
$C^{"}\left(\omega\right)=\omega^{s}\exp\left(-\omega/\omega_{c}\right)$ with
$s=2$, $\omega_{\text{c}}=100\ \text{cm}^{-1}$ (Kell  _et al._ , 2013;
Jakučionis  _et al._ , 2018). The number of modes and discretization
parameters are sufficient to obtain convergent model dynamics. For
thermalization, we consider scattering rates of all modes to be equal
$\nu_{mq}\rightarrow\nu$, and the scattering step size $\tau=0.01\
\text{ps}^{-1}$. Thermal ensemble consists of $5000$ trajectories.

Figure 1: Phase space trajectory of one specific bath mode $\omega_{q}=100\
\text{cm}^{-1}$ for a single excited chromophore calculated with various
scattering rates $\nu$. Initial temperature of the primary bath is
$T_{1}\left(0\right)=300\ \text{K}$ and the secondary bath is at a constant
temperature $T_{\infty}=200\ \text{K}$. Scattering step size is $\tau=0.01\
\text{ps}$. Wiggles in dynamics are due to finite size ensemble averaging
(5000 trajectories).

In Fig. (1) the coordinate $\left\langle
x_{1q}\right\rangle_{\text{th}}=\left\langle\sqrt{2}\text{Re}\lambda_{1q}\right\rangle_{\text{th}}$
and momentum $\left\langle
p_{1q}\right\rangle_{\text{th}}=\left\langle\sqrt{2}\text{Im}\lambda_{mq}\right\rangle_{\text{th}}$
phase space trajectory of a single $100\ \text{cm}^{-1}$ frequency vibrational
mode, calculated with various scattering rates $\nu$ is presented. The
oscillator, in the absence of thermalization, evolves along a closed
trajectory around $x_{1q}^{\text{min}}=\sqrt{2}g_{1q}$. Applying
thermalization procedure, a dissipative type trajectory is observed. The
coordinate $\left\langle x_{1q}\right\rangle_{\text{th}}$ equilibrates to
$x_{1q}^{\text{min}}$ (equilibrium is shifted from zero due to coupling with
the system), while momentum $\left\langle p_{1q}\right\rangle_{\text{th}}$
approaches zero. The thermalization time can be adjusted by changing the
scattering rate, $\nu$. Both weakly damped and overdamped regimes become
available.

Figure 2: The primary bath temperature $T_{1}\left(t\right)$ calculated with
various scattering rates $\nu$. Initial temperature of the primary bath is
$T_{1}\left(0\right)=300\ \text{K}$ and the secondary bath is at a constant
temperature $T_{\infty}=200\ \text{K}$. Scattering step size is $\tau=0.01\
\text{ps}$.

Transient temperature of the primary bath can be estimated (Abramavicius  _et
al._ , 2018) by computing the average kinetic energy $\left\langle
K_{mq}\left(t,\epsilon\right)\right\rangle_{\text{th}}$ over time interval
$\epsilon$. The parameter $\epsilon$ then implies the resolution. For the
whole primary bath the transient temperature is then given by

$T_{m}\left(t\right)=\frac{1}{k_{\text{B}}Q}\sum_{q=1}^{Q}\omega_{mq}\ln\left(1+\frac{\omega_{mq}}{2\left\langle
K_{mq}\left(t,\epsilon\right)\right\rangle_{\text{th}}}\right)^{-1}\ .$ (9)

In Fig. (2) we present the primary bath temperature calculated with
$\epsilon=50\ \text{fs}$ and various scattering rates, $\nu$. In the absence
of thermalization, the primary bath temperature remains at the initial value
of $T_{1}\left(0\right)=300\ \text{K}$. Meanwhile, thermalization introduces
cooling of the primary bath down to the temperature of the secondary bath. The
scattering rate, $\nu$, allows to control the thermalization time.

The temperature control and stability considerably affects the electronic
excitation dynamics. To demonstrate the sensitivity of the excitation
evolution to the thermalization we consider a linear $N=3$ chromophore
aggregate, with chromophore transition energies $0,\ 250,\ 500\
\text{cm}^{-1}$, and nearest neighbor coupling $J=100\ \text{cm}^{-1}$.
_Excited states of such chromophore aggregate are_ excitons (Valkunas  _et
al._ , 2013; van Amerongen  _et al._ , 2000). They represent electronic
excitations delocalized over several sites with time dependent delocalization
length (Chorošajev  _et al._ , 2016). Hence, we switch to the eigenstate basis
(exciton representation, defined by
$\hat{H}_{\text{S}}\psi^{(\mathrm{exc})}=\varepsilon\psi^{(\mathrm{exc})}$):
$\rho_{e}^{\text{(exc)}}\left(t\right)=\sum_{n,m}\left(\psi_{ne}^{(\mathrm{exc})}\right)^{\star}\left\langle\alpha_{n}^{\ast}\left(t\right)\alpha_{m}\left(t\right)\right\rangle_{th}\psi_{me}^{(\mathrm{exc})}$.
The initial electronic state correspond to the optically excited highest
energy exciton eigen state. The parameters of the primary baths of
chromophores are the same as above, however, now the initial primary bath
temperature and the secondary bath temperature is the same
$T_{m}\left(0\right)=T_{\infty}=77\ \text{K}$. The thermal ensemble consists
of $240$ trajectories. In Fig. (3) we present exciton state populations
$\rho_{e}^{\text{exc}}\left(t\right)$ and the primary bath temperatures
$T_{m}\left(t\right)$ calculated in (i) the dense primary bath discretization
regime without thermalization (the bath discretization step size is
$\Delta\omega=1\ \text{cm}^{-1}$, $Q=750$ vibrational modes per site), in (ii)
the sparse discretization regime without thermalization ($\Delta\omega=50\
\text{cm}^{-1}$, $Q=15$) and (iii) the sparce discretized bath with
thermalization ($\nu=2.5\ \text{ps}^{-1}$).

Figure 3: Multi site-bath model exciton state populations
$\rho_{n}^{\text{exc}}\left(t\right)$ and local bath temperatures
$T_{m}\left(t\right)$ calculated in (i) the dense primary bath discretization
regime without thermalization, (ii) the sparse discretization regime without
thermalization and (iii) the sparcely discretized bath with thermalization
($\nu=2.5\ \text{ps}^{-1}$).

Consider the excitation dynamics without thermalization. In models (i) and
(ii) exciton populations sequentially relax to lower energy exciton states,
eventually reaching the lowest energy state. (Moix  _et al._ , 2012; Subasi
_et al._ , 2012; Gelzinis and Valkunas, 2020). The final population
distribution in the sparse regime, model (ii), significantly differs from the
dense case. Origin of the discrepancy is two-fold: the bath recursion time
$t_{\text{rec}}=2\pi/\Delta\omega$ for model (ii) is shorter than the
calculation time $t_{\text{rec}}<t_{\text{total}}$, and the sparse primary
bath shows significant growth of the bath temperature, compare Fig. (3B and
3D). Both of these drawbacks are addressed by introducing the bath
thermalization in model (iii). Looking at Fig. (3E), we see that the exciton
population dynamics and steady state values for model (iii) become
quantitatively comparable to case of model (i).

##  IV Discussion

A single quantum harmonic oscillator is characterized by a specific heat
$c(\beta^{-1})<k_{\text{B}}$, which depends on temperature as given by Eq.
(1). For a given set of bath oscillators the specific heat at a given
temperature can be estimated, however, the harmonic oscillators of the bath as
defined by Eq. (3) do not exchange energy. Accordingly, as the system relaxes,
only a few in-resonant oscillators accept the energy and diverges away from
equilibrium (Chin  _et al._ , 2013). Hence, the temperature at which
excitation dynamics occur no longer match the initial bath temperature - local
heating takes place.

A straightforward approach to avoid heating is to increase the bath density of
states until dynamics of interest converges (in our model, this is achieved by
increasing the number of bath oscillators). However, this is acceptable only
for small systems, since computation effort scales quadratically with both
number of sites and bath oscillators. Thermalization can be utilized to steer
the bath to the required temperature. Additional merit of thermalization is
the significant reduction of the the number of vibrational modes needed per
bath. Our simulations show convergence with just $15$ modes per bath while
maintaining comparable exciton relaxation dynamics (Fig. 3).

In effort to reduce the computational effort, Wang et al. have used (Wang  _et
al._ , 2016) a logarithmic bath discretization. However, high frequency
representation of the continuous spectral density becomes poor. Our model is
in line with explicit surrogate Hamiltonian (Baer and Kosloff, 1997) and its
stochastic realization (Katz  _et al._ , 2008; Torrontegui and Kosloff, 2016;
Habecker  _et al._ , 2019), while our approach does not require the explicit
modeling of the secondary bath, it still maintaining proper quantum dynamics
in the system.

The time-dependent variational approach with Davydov $\text{D}_{2}$ ansatz can
be improved by considering more complex Davydov ans$\ddot{a}$tze family
members, e.g., multitude of $\text{D}_{1}$ ansatz (multi-$\text{D}_{1}$) and
multi-$\text{D}_{2}$ (Zhou  _et al._ , 2015; Wang  _et al._ , 2016; Zhou  _et
al._ , 2016; Chen  _et al._ , 2018) or its Born-Oppenheimer approximated
variant (Jakučionis  _et al._ , 2020), $\text{sD}_{2}$. In either way, they
all suffer from finite bath heating capacity, in most cases, even stronger
than the $\text{D}_{2}$ ansatz, because of significantly increased
computational effort needed to propagate numerous bath oscillators. Work is in
progress on adapting the presented thermalization algorithm to these more
intricate ans$\ddot{a}$tze.

In conclusion, we present a system-bath model with stochastic bath
thermalization using the time-dependent variational approach with Davydov
$\text{D}_{2}$ ansatz. Thermalization allows to steer bath vibrational modes
evolution towards equilibrium thermal state of selected temperature in a
controlled way, and at the same time, for the bath to still maintain an aspect
of being coupled to the system. In addition, by analysing exciton relaxation
dynamics of a chromophore aggregate with thermalization, we found the exciton
dynamics to converge with much smaller number of bath modes, significantly
speeding up numerical computation.

###### Acknowledgements.

We thank the Research Council of Lithuania for financial support (grant No:
SMIP-20-47). Computations were performed on resources at the High Performance
Computing Center, “HPC Sauletekis” in Vilnius University Faculty of Physics.

## References

  * Dorfman  _et al._ (2016) K. E. Dorfman, F. Schlawin, and S. Mukamel, Rev. Mod. Phys. 88, 45008 (2016).
  * Mukamel (1995) S. S. Mukamel,  _Principles of nonlinear optical spectroscopy_ (Oxford University Press, 1995) p. 543.
  * Jakučionis  _et al._ (2018) M. Jakučionis, V. Chorošajev, and D. Abramavičius, Chem. Phys. 515, 193 (2018).
  * Chorošajev  _et al._ (2017) V. Chorošajev, T. Marčiulionis, and D. Abramavicius, J. Chem. Phys. 147, 74114 (2017).
  * Schröter  _et al._ (2015) M. Schröter, S. Ivanov, J. Schulze, S. Polyutov, Y. Yan, T. Pullerits, and O. Kühn, Phys. Rep. 567, 1 (2015).
  * Chorošajev  _et al._ (2014) V. Chorošajev, A. Gelzinis, L. Valkunas, and D. Abramavicius, J. Chem. Phys. 140, 244108 (2014).
  * Valkunas  _et al._ (2013) L. Valkunas, D. Abramavicius, and T. Mančal,  _Molecular Excitation Dynamics and Relaxation_ (Wiley-VCH, 2013).
  * May and Kühn (2011) V. May and O. Kühn,  _Charge and Energy Transfer Dynamics in Molecular Systems_ (Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2011) p. 562.
  * Jang and Mennucci (2018) S. J. Jang and B. Mennucci, Rev. Mod. Phys. 90, 35003 (2018).
  * Thyrhaug  _et al._ (2018) E. Thyrhaug, C. N. Lincoln, F. Branchi, G. Cerullo, V. Perlík, F. Šanda, H. Lokstein, and J. Hauer, Photosynth. Res. 135, 45 (2018).
  * Malý  _et al._ (2016) P. Malý, J. M. Gruber, R. J. Cogdell, T. Mančal, and R. Van Grondelle, Proc. Natl. Acad. Sci. U.S.A. 113, 2934 (2016).
  * Chenu and Scholes (2015) A. Chenu and G. D. Scholes, Annu. Rev. Phys. Chem. 66, 69 (2015).
  * Sjakste  _et al._ (2018) J. Sjakste, K. Tanimura, G. Barbarino, L. Perfetti, and N. Vast, J. Phys. Condens. Matter 30, 353001 (2018).
  * Flesch  _et al._ (2008) A. Flesch, M. Cramer, I. P. McCulloch, U. Schollwöck, and J. Eisert, Phys. Rev. A 78, 033608 (2008).
  * Lombardo and Turiaci (2013) F. C. Lombardo and G. J. Turiaci, Phys. Rev. D 87, 84028 (2013).
  * Yu and Zhang (2008) H. Yu and J. Zhang, Phys. Rev. D 77, 24031 (2008).
  * Ruostekoski and Walls (1998) J. Ruostekoski and D. F. Walls, Phys. Rev. A 58, R50 (1998).
  * Breuer and Petruccione (2007) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007) p. 636.
  * Weiss (2012) U. Weiss,  _Quantum Dissipative Systems_ (World Scientific, 2012).
  * Montoya-Castillo  _et al._ (2015) A. Montoya-Castillo, T. C. Berkelbach, and D. R. Reichman, J. Chem. Phys. 143, 194108 (2015).
  * Jang  _et al._ (2002) S. Jang, J. Cao, and R. J. Silbey, J. Chem. Phys. 116, 2705 (2002).
  * Tanimura (1990) Y. Tanimura, Phys. Rev. A 41, 6676 (1990).
  * Tanimura (2006) Y. Tanimura, J. Phys. Soc. Japan 75, 82001 (2006).
  * Xu and Yan (2007) R.-X. Xu and Y.J. Yan, Phys. Rev. E 75, 031107 (2007).
  * Tamascelli  _et al._ (2019) D. Tamascelli, A. Smirne, J. Lim, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 123, 090402 (2019).
  * Prior  _et al._ (2010) J. Prior, A. W. Chin, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 105, 050404 (2010).
  * Abramavicius and Abramavicius (2014) V. Abramavicius and D. Abramavicius, J. Chem. Phys 140, 065103 (2014).
  * Appel and Di Ventra (2009) H. Appel and M. Di Ventra, Phys. Rev. B 80, 212303 (2009).
  * Biele and D’Agosta (2012) R. Biele and R. D’Agosta, J. Phys. Condens. Matter 24, 273201 (2012).
  * Diósi and Strunz (1997) L. Diósi and W. T. Strunz, Phys. Lett. Sect. A Gen. At. Solid State Phys. 235, 569 (1997).
  * de Vega  _et al._ (2005) I. de Vega, D. Alonso, and P. Gaspard, Phys. Rev. A - At. Mol. Opt. Phys. 71, 023812 (2005).
  * Link and Strunz (2017) V. Link and W. T. Strunz, Phys. Rev. Lett. 119, 180401 (2017).
  * Hartmann and Strunz (2017) R. Hartmann and W. T. Strunz, J. Chem. Theory Comput. 13, 5834 (2017).
  * Reddy and Prasad (2015) C. S. Reddy and M. D. Prasad, Mol. Phys. 113, 3023 (2015).
  * Ritschel  _et al._ (2015) G. Ritschel, D. Suess, S. Möbius, W. T. Strunz, and A. Eisfeld, J. Chem. Phys. 142, 034115 (2015).
  * Borrelli and Gelin (2016) R. Borrelli and M. F. Gelin, J. Chem. Phys. 145, 224101 (2016).
  * Chen and Zhao (2017) L. Chen and Y. Zhao, J. Chem. Phys. 147, 214102 (2017).
  * Martinazzo  _et al._ (2006) R. Martinazzo, M. Nest, P. Saalfrank, and G. F. Tantardini, J. Chem. Phys. 125, 194102 (2006).
  * Abramavicius  _et al._ (2018) D. Abramavicius, V. Chorošajev, and L. Valkunas, Phys. Chem. Chem. Phys. 20, 21225 (2018).
  * Chen and Sorbello (1993) Z. Chen and R. S. Sorbello, Phys. Rev. B 47, 13527 (1993).
  * Ichikawa  _et al._ (2007) M. Ichikawa, H. Ichikawa, K. Yoshikawa, and Y. Kimura, Phys. Rev. Lett. 99, 148104 (2007).
  * Gulbinas  _et al._ (1996) V. Gulbinas, L. Valkunas, D. Kuciauskas, E. Katilius, V. Liuolia, W. Zhou, and R. E. Blankenship, J. Phys. Chem. 100, 17950 (1996).
  * Valkunas and Gulbinas (1997) L. Valkunas and V. Gulbinas, Photochem. Photobiol. 66, 628 (1997).
  * Gulbinas  _et al._ (2006) V. Gulbinas, R. Karpicz, G. Garab, and L. Valkunas, Biochemistry 45, 9559 (2006).
  * Balevičius  _et al._ (2018) V. Balevičius, C. N. Lincoln, D. Viola, G. Cerullo, J. Hauer, and D. Abramavicius, Photosynth. Res. 135, 55 (2018).
  * Balevičius Jr  _et al._ (2019) V. Balevičius Jr, T. Wei, D. Di Tommaso, D. Abramavicius, J. Hauer, T. Polívka, and C. D. P. Duffy, Chem. Sci. 10, 4792 (2019).
  * Choi  _et al._ (2019) J. Choi, H. Zhou, S. Choi, R. Landig, W. W. Ho, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, D. A. Abanin, and M. D. Lukin, Phys. Rev. Lett. 122, 043603 (2019).
  * Scarani  _et al._ (2002) V. Scarani, M. Ziman, P. Štelmachovič, N. Gisin, and V. Bužek, Phys. Rev. Lett. 88, 097905 (2002).
  * Frenkel and Frenkel (1936) I. Frenkel and J. Frenkel,  _Wave Mechanics: Elementary Theory_ (Oxford University Press, 1936).
  * Sun  _et al._ (2010) J. Sun, B. Luo, and Y. Zhao, Phys. Rev. B - Condens. Matter Mater. Phys. 82, 014305 (2010).
  * Zhang  _et al._ (1990) W.-M. Zhang, D. H. Feng, and R. Gilmore, Rev. Mod. Phys 62, 867 (1990).
  * Kais and Levine (1990) S. Kais and R. D. Levine, Phys. Rev. A 41, 2301 (1990).
  * Chen  _et al._ (2018) L. Chen, M. Gelin, and Y. Zhao, Chem. Phys. 515, 108 (2018).
  * Yan  _et al._ (2020) Y. Yan, L. Chen, J. Y. Luo, and Y. Zhao, Phys. Rev. A 102, 023714 (2020).
  * Glauber (1963) R. J. Glauber, Phys. Rev. 131, 2766 (1963).
  * Plenio and Knight (1998) M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998).
  * Luoma  _et al._ (2020) K. Luoma, W. T. Strunz, and J. Piilo, Phys. Rev. Lett. 125, 150403 (2020).
  * Kampen (2007) V. N. Kampen,  _Stochastic Processes in Physics and Chemistry_ (Elsevier, 2007).
  * Bertsekas and Tsitsiklis (2008) D. Bertsekas and J. Tsitsiklis,  _Introduction to probability_ (Athena Scientific, 2008).
  * Kell  _et al._ (2013) A. Kell, X. Feng, M. Reppert, and R. Jankowiak, J. Phys. Chem. B 117, 7317 (2013).
  * van Amerongen  _et al._ (2000) H. van Amerongen, R. van Grondelle, and L. Valkunas,  _Photosynthetic Excitons_ (World Scientific, 2000).
  * Chorošajev  _et al._ (2016) V. Chorošajev, O. Rancova, and D. Abramavicius, Phys. Chem. Chem. Phys. 18, 7966 (2016).
  * Moix  _et al._ (2012) J. M. Moix, Y. Zhao, and J. Cao, Phys. Rev. B 85, 115412 (2012).
  * Subasi  _et al._ (2012) Y. Subasi, C. H. Fleming, J. M. Taylor, and B. L. Hu, Phys. Rev. E 86, 61132 (2012).
  * Gelzinis and Valkunas (2020) A. Gelzinis and L. Valkunas, J. Chem. Phys. 152, 51103 (2020).
  * Chin  _et al._ (2013) A. W. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. F. Huelga, and M. B. Plenio, Nat. Phys. 9, 113 (2013).
  * Wang  _et al._ (2016) L. Wang, L. Chen, N. Zhou, and Y. Zhao, J. Chem. Phys. 144, 024101 (2016).
  * Baer and Kosloff (1997) R. Baer and R. Kosloff, J. Chem. Phys. 106, 8862 (1997).
  * Katz  _et al._ (2008) G. Katz, D. Gelman, M. A. Ratner, and R. Kosloff, J. Chem. Phys. 129, 034108 (2008).
  * Torrontegui and Kosloff (2016) E. Torrontegui and R. Kosloff, New J. Phys. 18, 093001 (2016).
  * Habecker  _et al._ (2019) F. Habecker, R. Röhse, and T. Klüner, J. Chem. Phys. 151, 134113 (2019).
  * Zhou  _et al._ (2015) N. Zhou, Z. Huang, J. Zhu, V. Chernyak, and Y. Zhao, J. Chem. Phys. 143, 014113 (2015).
  * Zhou  _et al._ (2016) N. Zhou, L. Chen, Z. Huang, K. Sun, Y. Tanimura, and Y. Zhao, J. Phys. Chem. A 120, 1562 (2016), arXiv:1908.09243 .
  * Jakučionis  _et al._ (2020) M. Jakučionis, T. Mancal, and D. Abramavičius, Phys. Chem. Chem. Phys. 22, 8952 (2020).


# Time-reversal symmetric topological metal

L. C. Xie  H. C. Wu  L. Jin<EMAIL_ADDRESS>Z. Song  School of
Physics, Nankai University, Tianjin 300071, China

###### Abstract

Topological metals possess gapless band structures accompanied with nontrivial
edge states. Topological metals are created from the topological insulators by
adjusting the magnetic flux. Here we propose the time-reversal symmetry
protected topological metal without employing the magnetic flux. Topological
metallic phase presents although the Chern number vanishes. In the
topologically nontrivial phase, three different cases in the metallic phase
are distinguished from the band-touching and in-gap features of the
topological edge states. These findings shed light on the time-reversal
symmetric topological metals and greatly simplify the realization of
topological metals.

##  I Introduction

The introduction of topology in the condensed matter physics greatly changes
the traditional way of classifying phases of matter, novel topological states
including insulators, superconductors, semimetals, and metals are discovered
E.Mele2007 ; Kane2010 ; XGWen2017 ; T.Das2016 ; SCZhang2011 ; L.Balents2011 ;
AMRappe2012 ; A.Vish2018 ; AB16 ; BYan17 ; CKChiu16 ; ABansil16 ; LJin19 ;
Kawabata19 ; Ashida20 . These topological phases all support robust edge
states, which are counted by the corresponding topological invariant that
capturing the topology of the bulk band DJTh1982 ; DXiao2010 ; Kane2010b . The
tremendous progresses in various platforms of photonics, acoustics, and
electronic circuits allow us to study these topological phases in experiments,
which provide crucial enlightenments for topological materials YRan2009 ;
LJLang2012 ; XJLiu13 ; MAide13 ; LLu14 ; SYYu18 ; YBYang18 ; SLZhu19 ;
XDChen19 ; YXu19 ; SImh18 ; THel19 ; YLu19 .

Topological insulators have a fully gapped bulk band structure like ordinary
insulators, but topological insulators support conductive edge states on the
surface of the sample, where electrons move along the surface of materials
Kane2010 . The Chern insulator is a typical two-dimensional (2D) topological
insulator, which supports robust chiral edge states that are predicted from
the nonzero Chern number. A landmark development of topological insulators is
the discovery of time-reversal invariant topological insulators S.M04 ;
CLKane05a ; BAB06 ; XLQi06 , where the Chern number vanishes and becomes
invalid for topological characterization; alternatively, the $Z_{2}$
topological number classifies the time-reversal invariant insulators CLKane05b
and these topological insulators are featured by their helical edge states,
being confirmed in nonlocal transport experiment ARoth09 . The helical edge
states have opposite spins and counter-propagate on the boundaries in the
opposite directions. The spin Chern number LSheng05 ; EProdan09 ; YYang11 ,
valley Chern number FZhang11 ; FZhang13 ; MEzawa13 , and wave polarization
RResta94 ; PDelp11 ; FLiu17  are developed to characterize the quantum spin
Hall effect.

Topological semimetals and metals are metallic states with nontrivial gapless
band structure and unique quasi-particle excitations A.Vish2018 , represented
by Weyl and Dirac semimetals. In the Weyl (Dirac) semimetals, energy bands are
double (fourfold) degenerately connected by the Weyl (Dirac) points SMura17 ;
XWan11 ; ZWang13 . In addition, the existence of other multiple degeneracy
points has been discovered, like threefold and sixfold ZZhu16 ; BJW16 ;
Hweng16 ; SSingh18 ; GShan19 ; RChapai19 . For systems protected by certain
symmetries, the band connection is a closed node-line instead of nodes. The
corresponding topological node-line semimetals have been reported in the
graphene network AAB11 ; RYu15 ; HMW15 .

Topological metals (TMs) are gapless phases featured by their separable
conduction and valence bands. Thus, the topological characterization for the
metals is similar as that for the insulators. TMs can be divided into nodal-
point, nodal-line, and nodal-surface TMs according to the dimension of band
crossing WWu18 ; or alternatively be divided into type I and type II TMs
according to the dispersion close to the band crossing AAS15 ; YXu15 ; SLi17 .
TMs in different spatial dimensions have been extensively studied MBahari19 ;
SYXu15 ; ZJWang16 . To create 2D topological metal, the traditional method is
through adjusting the magnetic flux to tune the energy bands of the
topological insulators XZYing18 . However, introducing and tuning the magnetic
flux is difficult in many experimental platforms. Here we demonstrate a time-
reversal symmetric topological metal without the magnetic flux.

In this paper, we report the time-reversal symmetric topological metal in a 2D
square lattice, which belongs to class AI in the tenfold way of topological
classification. The 2D square lattice possesses nontrivial topology although
the Chern number is zero ensured by the time-reversal symmetry. The proposed
topological metal is free from the magnetic flux and significantly differs
from the traditional topological metal created through using the magnetic
flux. The reflection symmetry ensures that the quantized polarization; the
zero polarization predicts the trivial topology in the vertical direction. The
Zak phase at the phase transition point characterizes the topology in the
horizontal direction. We demonstrate three different cases in the
topologically nontrivial phases that are distinguished from their different
types of topological edge states. The topological phases as well as the
topological edge states can be deformed into each other without undergoing
topological phase transition. Our findings provide novel insights into time-
reversal symmetric topological metals.

The remainder of the paper is organized as follows. In Sec. II, the time-
reversal symmetric 2D square lattice is introduced. Section III presents the
phase diagram and the energy band of the topological metallic phase. Three
different topological edge states are elucidated. Section IV provides the
topological characterization. The discussion and conclusion are summarized in
Secs. V and VI.

##  II Time-reversal symmetric 2D square lattice

Figure 1: Schematic of the 2D square lattice. The unit cell is consist of two
sublattices $A$ and $B$ in cyan and magenta as indicated inside the dashed
blue rectangle. $\lambda$, $\kappa$, and $m$ are the couplings between the
sublattices.

We consider a time-reversal symmetric 2D square lattice. The schematic of the
lattice is illustrated in Fig. 1. The Hamiltonian of the 2D square lattice in
the real space is written in the form of

$\displaystyle H$ $\displaystyle=$ $\displaystyle\sum\limits_{l,j}(\lambda
a_{l,j}^{\dagger}b_{l,j}+\kappa
a_{l,j}^{\dagger}b_{l-1,j}+ma_{l,j}^{\dagger}b_{l,j+1}+$ (1) $\displaystyle
ma_{l,j}^{\dagger}b_{l,j-1}+\text{H.c.})+\lambda
a_{l,j}^{\dagger}(a_{l,j+1}+a_{l,j-1})$ $\displaystyle+\kappa
b_{l,j}^{\dagger}(b_{l,j+1}+b_{l,j-1}),$

where $a^{\dagger}$, $b^{\dagger}$ $(a$, $b)$ represent the creation
(annihilation) operators for the sublattices $A$ and $B$. The dashed box in
blue indicates the unit cell of the system. The couplings $\lambda$ and
$\kappa$ indicate the nearest neighbor coupling strengths, the coupling $m$ is
the next nearest neighbor coupling strength. The 2D square lattice does not
enclose any magnetic flux and can be realized in many versatile experimental
platforms for the topological insulators including the ultracold atoms in the
optical lattice, the photonic crystals of coupled optical
waveguides/resonators, the acoustic metamaterial dubbed sonic crystals, and
the electric circuits MAide15 ; SMi18 ; THofma19 ; CHe16 ; ZJYang15 ; YGPeng16
; JYLu18 .

Applying the Fourier transformation
$\rho_{l,j}=M^{-1}\sum_{\mathbf{k}}e^{i\mathbf{k\cdot r}}\rho_{\mathbf{k}}$
$(\rho=a,b)$ to the two sublattices, we rewrite the lattice Hamiltonian $H$ in
the momentum space in the form of

$H=\sum_{\mathbf{k}}H\left(\mathbf{k}\right)=\sum_{\mathbf{k}}\psi_{\mathbf{k}}^{\dagger}h\left(\mathbf{k}\right)\psi_{\mathbf{k}},$
(2)

where the basis is $\psi_{\mathbf{k}}=(a_{\mathbf{k}},b_{\mathbf{k}})^{T}$ and
the Bloch Hamiltonian $h\left(\mathbf{k}\right)$ is a $2\times 2$ matrix

$h\left(\mathbf{k}\right)=d_{0}\sigma_{0}+\mathbf{d\cdot\sigma,}$ (3)

where $\sigma_{0}$ is the $2\times 2$ identical matrix,
$d_{0}=(\lambda+\kappa)\cos k_{y}$.
$\mathbf{\sigma=(\sigma}_{x},\mathbf{\sigma}_{y},\mathbf{\sigma}_{z}\mathbf{)}$
is the Pauli matrices, and the effective magnetic field $\mathbf{d}$ is

$\displaystyle d_{x}$ $\displaystyle=$ $\displaystyle\lambda+2m\cos
k_{y}+\kappa\cos k_{x},$ $\displaystyle d_{y}$ $\displaystyle=$
$\displaystyle\kappa\sin k_{x},$ (4) $\displaystyle d_{z}$ $\displaystyle=$
$\displaystyle(\lambda-\kappa)\cos k_{y}.$

The energies of the Bloch bands are

$E_{\pm}=d_{0}\pm\sqrt{d_{x}^{2}+d_{y}^{2}+d_{z}^{2}},$ (5)

where the subscripts $+$ and $-$ of energy represent the upper and lower
bands, respectively. Without loss of generality, we fix the coupling $m$ as
unity. Notably, the band gap closes at $\lambda=\pm\kappa$, associated with
the appearance of degenerate points (DPs) or degenerate lines.

The Bloch Hamiltonian $h\left(\mathbf{k}\right)$ of the square lattice is
protected by the time-reversal symmetry. Thus, the band energies satisfy
$E_{\pm,\mathbf{k}}=E_{\pm,-\mathbf{k}}^{\ast}$ as reflected from Eq. (5). The
symmetries are important CKChiu16 ; Kawabata19 ; LJin21 . Notably, the time-
reversal symmetry ensures a zero Chern number for the energy band; however,
the topologically nontrivial phase still exists in the 2D lattice FLiu17 ;
HCWu20 . $d_{0}$ adjusts the energy band and changes the band structure, and
plays a crucial role in creating the metallic phase when the strengths of the
couplings $\lambda$ and $\kappa$ are not equal $\lambda^{2}\neq\kappa^{2}$. In
the absence of $d_{0}$, the bands are fully gapped except for the band
touching and the system is in an insulator phase.

The energy bands have highly degenerate points, where the energy bands become
fully flat. The highly degenerate points appear at
$k_{y}=\pm\arccos(-\lambda/2)$, but disappear as $\left|\lambda\right|>2$.
Time-reversal symmetry ensures that the two highly degenerate points are
symmetric about $k_{y}=0$. The upper and lower bands reduce to two momentum
independent flat bands

$\varepsilon_{\pm}=-(\lambda^{2}+\kappa\lambda\pm\sqrt{\Delta})/2,$ (6)

where
$\Delta=4\kappa^{2}+\kappa^{2}\lambda^{2}-2\kappa\lambda^{3}+\lambda^{4}$.
Notice that the presence of $d_{0}$ breaks the symmetry of the upper and lower
bands about zero energy.

We emphasize that the proposed topological metal significantly differs from
the traditional topological metal created through altering the magnetic flux.
Their symmetry protections, topological classifications, and topological
invariants are different. The proposed topological metal without magnetic flux
is time-reversal symmetric and belongs to the symmetry class AI SRu10 . In
contrast, the topological metal with magnetic flux breaks the time-reversal
symmetry; for example, the topological metal in Ref. XZYing18  created via
tuning the magnetic flux from the ordinary metal is particle-hole symmetric
and belongs to the symmetry class D. Consequently, the topological invariants
are different. In addition, the energy bands are symmetric about $k_{y}=0$ in
the time-reversal symmetric topological metal. Thus, the edge states cross
twice at zero energy as a consequence of zero Chern number $C=0$; and the
dispersions and velocities near the zero edge states on the two sides of
$k_{y}=0$ are opposite. We can observe these features in the following
section. However, the edge states cross once at zero energy in the traditional
topological metal without the time-reversal symmetry protection as a
consequence of nonzero Chern number $C=\pm 1$.

##  III Metal-insulator transition

In this section, we discuss the topological phases. The prominent feature of
the considered 2D square lattice is the presence of topological metallic
phase. For a two-band model, if the bulk bands are separable and the Fermi
level lies in the band gap, the system is an insulator. Once the highest
energy level of the valence band exceeds the lowest energy level of the
conduction band, the two band energies partly overlap each other and the Fermi
level crosses the two bands. Then, the system completes the transition from
insulator to metal, known as the Wilson transition. Figure 2 is the phase
diagram, being inversion invariant about the origin $(\kappa,\lambda)=(0,0)$.
The dashed lines represent the topological phase transition, and the solid
black lines represent the Wilson transition. The purple and cyan regions are
the metallic and insulator phases, respectively; being the topologically
nontrivial phases. Similarly, the green and yellow regions are the metallic
and insulator phases, respectively; being the trivial phases.

The topological insulator phase is indicated by the cyan region in the phase
diagram Fig. 2. The two bands are fully gapped as shown in the energy spectra
Figs. 3(a)-3(c). The topological metallic phase is indicated by the purple
region in the phase diagram Fig. 2. The lowest energy of the conduction band
is lower than the highest energy of the valence band as shown in the energy
spectra Figs. 3(d)-3(f). When the system in the insulator phase crosses the
solid black lines in the phase diagram Fig. 2, the system undergoes a Wilson
phase transition and enters the metallic phase. However, the band gap is
always open, which ensures that the topological characteristics remain
unchanged. This indicates that the topological metallic phase and the
topological insulator phase have identical topological properties.

Figure 2: Phase diagram of $h\left(\mathbf{k}\right)$ in the
$\kappa$-$\lambda$ parameter space. Regions with different colors indicate
different phases, topological metal (TM) in purple; ordinary metal (OM) in
green; topological insulator (TI) in cyan; ordinary insulator (OI) in yellow.
The black solid lines indicate the Wilson transitions. The red dashed lines
indicate the topological phase transition with two DPs, while the blue dashed
line indicates the topological phase transition with four DPs. The brown solid
lines are the boundary that characterizes the positional relationship between
the edge states and bulk bands, indicating three different cases in the
topological metallic or insulator phase. Figure 3: Band structures of three
cases in the insulator phase at (a) $\kappa=1,\lambda=-0.5$, (b)
$\kappa=2,\lambda=-1$, (c) $\kappa=5,\lambda=-1$ and three cases in the
metallic phase at (d) $\kappa=1,\lambda=0.5$, (e) $\kappa=2,\lambda=1$, (f)
$\kappa=5,\lambda=1$ as indicated by the black solid dots in the phase diagram
of Fig. 2. The energy band is depicted for the OBC in the $x$ direction and
the PBC in the $y$ direction. The other system parameter is $m=1$. The red
line indicates the edge state. Figure 4: Plot of the total inverse
participation ratio along the momentum $k_{y}$ as the counterpart of Fig. 3.
The red line represents the TIPR including the localized edge state, and the
black line represents the TIPR of extended states in bulk. (a)
$\kappa=1,\lambda=-0.5$, (b) $\kappa=2,\lambda=-1$, (c) $\kappa=5,\lambda=-1$,
(d) $\kappa=1,\lambda=0.5$, (e) $\kappa=2,\lambda=1$, (f)
$\kappa=5,\lambda=1$. Figure 5: Plot of the local density of state near the
zero energy as the counterpart of Fig. 3. (a) $\kappa=1,\lambda=-0.5$, (b)
$\kappa=2,\lambda=-1$, (c) $\kappa=5,\lambda=-1$, (d) $\kappa=1,\lambda=0.5$,
(e) $\kappa=2,\lambda=1$, (f) $\kappa=5,\lambda=1$. Here we take $\alpha=2$
and lattice size is $30\times 30$. The color bar indicates the probability of
wave function and the maximum is renormalized to 1.

Interestingly, the system has three different types of edge states featured
from the relation between the edge states and the bulk bands in both the
topological metallic phase and the topological insulator phase. The edge
states always intersect at the momentum $k_{y}=\pm\pi/2$, which is also the
band closing point. These are observed from the energy spectra of Figs.
3(a)-3(c) and Figs. 3(d)-3(f), respectively. The different types of edge
states can be deformed into each other without the band gap closing; thus,
they possess identical topological properties. These edge states are
distinguished by the lines $\lambda+\kappa\pm 2=0$, $\lambda-\kappa\pm 2=0$ as
the solid brown lines in the phase diagram of Fig. 2. The topological
insulator phase in the $\kappa>0$ is divided into three regions; the spectra
for points $a$, $b$, and $c$ are representative cases. In Fig. 3(a), the
system has two pairs of gapless edge states, connecting the upper and the
lower bulk bands. In Fig. 3(b), the gapless edge states near the center of the
Brillouin zone (BZ) detach the bulk bands, forming two edge states cross twice
at $k_{y}=\pm\pi/2$ with each edge state connected to one bulk band,
respectively. In Fig. 3(c), the edge states are completely detached from the
bulk bands and become the in-gap edge states. The edge states in the
topological metallic phase inherit the characteristics of the edge states in
the topological insulator phase. Figures 3(d) to 3(f) show the three typical
energy spectra for the topological metallic phase. The topological metal still
has separable energy bands. In Fig. 3(d), the gapless edge states intact the
bulk bands. In Fig. 3(e), the edge states detach the bulk bands near the
boundary of the BZ. In Fig. 3(f), the edge states are completely detached from
the bulk band and become in-gap edge states.

For the three cases of insulator or metal as separated by the brown solid
lines in Fig. 2, we introduce the total inverse participation ratio (TIPR) as
a criterion to distinguish their different features. The existence of edge
states can be measured by the quantity of TIPR, which is defined as

$\text{TIPR}(k_{y})=\sum_{n,l}\left|\langle
l\left|\psi_{n,k_{y}}\right\rangle\right|^{4},$ (7)

where $n$ is the band index and $\left|\psi_{n,k_{y}}\right\rangle$ is
eigenstate. In the large size limit, the value of TIPR$(k_{y})$ is stable. The
distribution of TIPR has three different configurations. The TIPR$(k_{y})$ has
relative large values in two separated regions, in single incomplete region,
and in the whole region of the Brillouin zone of $k_{y}$. The TIPR$(k_{y})$ in
the topological insulator phase is shown in Figs. 4(a)-4(c) and the
TIPR$(k_{y})$ in the topological metal phase is shown in Figs. 4(d)-4(f). The
distribution can be characterized by the scattering process, where the mid-gap
levels are detected. Each region opens a window in the momentum space (or
energy) for the channel of electron transport, and then can be detected by
measuring the edge electrical conductivity and thermal conductivity XWen19 .

The insulator phase and metallic phase can be distinguished from the local
density of states (LDOS), which describes a space-resolved density of states.
The LDOS at the position $x$ and the energy $E$ is defined as
$D(E,x)=\sum_{n}\left|\psi_{n}(x)\right|^{2}\delta(E-\varepsilon_{n})$ in a
continuous system, and the LDOS reduces to the form

$D(E,l)=\frac{1}{\chi}\sum_{n=1}^{N}e^{-\alpha^{2}\left(E_{n}-E\right)^{2}}\left|\langle
l\left|\psi_{n}\right\rangle\right|^{2},$ (8)

at the lattice site ${l}$ and the energy ${E}$ in a discrete system. In Eq.
(8), the delta-function is replaced by a Gaussian function and $\chi$ is the
renormalization factor. Here, we consider the half-filling situation. The LDOS
at the Fermi surface $E_{F}=0$ is calculated for the finite size systems with
several representative parameters; the results are shown in Fig. 5. $D(E,l)$
has vanishing distributions in all the insulating phases as shown in Figs.
5(a)-5(c), but has finite distributions in all the conducting phases as shown
in Figs. 5(d)-5(f). In addition, ${D(E,l)}$ has evidently large distributions
along only two sides in all the phases, resulting the relatively large
electrical and thermal conductivities. The LDOS is detected by a scanning
tunneling microscope (STM), which is capable of imaging electron densities of
states with atomic resolution NHLe20 . The TIPR and LDOS are measurable
quantities in experiment. Furthermore, the nonzero bulk current of the
topological metal differs from that of the insulator and the relatively larger
edge current of the topological metal differs from that of the normal
conductor.

In the topological phases, the edge states exist for the system under the
periodic boundary condition (PBC) in the $y$ direction and the open boundary
condition (OBC) in the $x$ direction. The energies of the edge states that
localized on the left and right boundaries are

$E_{L}(k_{y})=2\lambda\cos k_{y},E_{R}(k_{y})=2\kappa\cos k_{y}.$ (9)

The wave function for the edge state with $E_{L}\left(k_{y}\right)$
[$E_{R}\left(k_{y}\right)$] localized at the left (right) boundary is denoted
as $\left|\psi_{L}\right\rangle$ ($\left|\psi_{R}\right\rangle$). We denote
the expression of the edge states as

$|\psi_{L,R}\rangle=\left(\psi_{1A},\psi_{1B},...,\psi_{NA},\psi_{NB}\right),$
(10)

where $N$ is the total number of the unit cells. The edge states at the
limitation of infinity large size system ($N\rightarrow\infty$) are
analytically obtained. For the edge state $|\psi_{L}\rangle$, the components
of $|\psi_{L}\rangle$ in the $n$-th unit cell are
$(\psi_{nA},\psi_{nB})=(\rho^{n-1},0)$ with the decay factor
$\rho=-(\lambda+2\cos k_{y})/\kappa$. For the edge state
$\left|\psi_{R}\right\rangle$, the components of $\left|\psi_{R}\right\rangle$
in the $n$-th unit cell are $(\psi_{nA},\psi_{nB})=(0,\rho^{N-n})$. The
expression is valid for all the three types of topological edge states.

##  IV Topological characterization

Nonspatial (internal) symmetry plays a key role in the topological
classification and determines the topological characterization CKChiu16 .
Time-reversal symmetry (TRS), particle-hole symmetry (PHS), and chiral
symmetry (CS) classify the topological phases in a ten-fold way SRu10 . The
time-reversal symmetry is defined as

$\mathrm{TRS}:\mathcal{T}h\left(\mathbf{k}\right)\mathcal{T}^{-1}=h\left(-\mathbf{k}\right),\mathcal{T}^{2}=\pm
1.$ (11)

The time-reversal operator $\mathcal{T}=UK$ is an anti-unitary operator, where
$U$ is a linear operator and the complex conjugation $K$ is an anti-linear
operator.

The TRS ensures that the Berry curvature $F$ satisfies
$F(\mathbf{k})=-F(-\mathbf{k})$, the integral of $F$ in the entire BZ results
in a vanishing Chern number. How to characterize the band topology and
determine the topologically nontrivial phase of a time-reversal symmetric
system? For the time-reversal symmetric system with $\mathcal{T}^{2}=-1$, the
topological properties can be characterized by a $Z_{2}$ topological
invariant. $Z_{2}=0$ represents the trivial insulator and $Z_{2}=1$ represents
the topological insulator. The Pfaffian method, time-reversal polarization,
and non-Abelian Berry connection have been developed as the $Z_{2}$ invariant,
but these are unable to describe the topological systems with
$\mathcal{T}^{2}=+1$. For the time-reversal symmetric system with
$\mathcal{T}^{2}=+1$, the 2D Zak phase, also known as the wave polarization,
is used for the topological characterization. The wave polarization is
quantized to $0$ or $1/2$ under the symmetry protection and is used to
describe the band topology of the time-reversal symmetric system with
$\mathcal{T}^{2}=+1$ FLiu17 . The topological edge states appear inside the
band gap HCWu20 .

For the Bloch Hamiltonian of the 2D square lattice $h\left(\mathbf{k}\right)$,
the linear operator $U=\sigma_{0}$ is the $2\times 2$ identical matrix; thus,
the 2D square lattice possesses the time-reversal symmetry of
$\mathcal{T}^{2}=+1$. The term $d_{0}\sigma_{0}$ in the Bloch Hamiltonian
$h\left(\mathbf{k}\right)$ is crucial for the appearance of the metallic
phase; however, $d_{0}\sigma_{0}$ only adjusts the band energy, but not the
eigenstates of $h\left(\mathbf{k}\right)$. The topology of
$h\left(\mathbf{k}\right)$ only relates to the eigenstate and thus is
determined from $\mathbf{d\cdot\sigma}$. The eigenstate of
$h\left(\mathbf{k}\right)$ is in the form of

$\left|\Psi_{\pm}(\mathbf{k})\right\rangle=\frac{1}{\Omega_{\pm}}\dbinom{d_{z}\pm
d}{d_{x}+id_{y}},$ (12)

where the subscript $+$ represents the upper band and the subscript $-$
represents the lower band; and $\Omega_{\pm}=\sqrt{2d(d\pm d_{z})}$ is the
normalization factor.

The wave polarization is defined as

$\mathbf{P=}\frac{1}{2\pi}\int_{BZ}\mathrm{Tr}[\mathbf{A}(k_{x},k_{y})]dk_{x}dk_{y},$
(13)

where
$\mathbf{A}(k_{x},k_{y})=\left\langle\Psi(\mathbf{k})\right|i\nabla_{\mathbf{k}}\left|\Psi(\mathbf{k})\right\rangle$
is the Berry connection. $\mathbf{P=(}P_{x},P_{y}\mathbf{)}$ includes two
components, where $P_{x}$ characterizes the topology in the $x$ direction and
predicts the edge states for the 2D square lattice under the OBC in the $x$
direction. $P_{y}$ characterizes the topology in the $y$ direction. The wave
polarization defines the 2D Zak phase $\Theta_{\nu}=2\pi P_{\nu}$,
($\nu=x,y$). Protected by the reflection symmetry or inversion symmetry, the
value of wave polarization exhibits quantized, being $0$ or $\pm 1/2$. For the
Bloch Hamiltonian $h(\mathbf{k})$, the reflection symmetry in the $y$
direction guarantees $P_{y}=0$. The Bloch Hamiltonian $h(\mathbf{k})$ is
reflection symmetric in the $y$ direction $h(k_{x},k_{y})=h(k_{x},-k_{y})$.
$P_{y}=0$ reflects the lack of nontrivial edge state for the 2D square lattice
under the OBC in the $y$ direction and the PBC in the $x$ direction.

Figure 6: Topological phase diagram of the time-reversal symmetric 2D square
lattice. The $|\lambda|<|\kappa|$ region in purple is the topologically
nontrivial phase; in contrast, the $|\lambda|>|\kappa|$ region in green is the
trivial phase. The Zak phase $\Theta_{x}$ distinguishes different topological
phases.

The polarization $P_{x}$ is not quantized without the symmetry protection.
Therefore, $P_{x}$ is invalid for the topological characterization. We apply a
Fourier transformation to decompose the 2D lattice Hamiltonian into decoupled
1D lattice Hamiltonian parameterized by the momentum $k_{y}$. The band gap
closes at $k_{y}=\pm\pi/2$, where the topological phase transition occurs and
the edge states intersect. For $\lambda=-\kappa$ $(\kappa\neq 0)$, the band
degenerate points locate at $(k_{x},k_{y})=(0,\pm\pi/2)$. For $\lambda=\kappa$
$(|\kappa|\leq 1)$, the band gap closes at
$(k_{x},k_{y})=(0,\pm\arccos\lambda)$ and $(\pi,\pm\pi/2)$. For
$\lambda=\kappa$ $(|\kappa|>1)$, the band degenerate points locate at
$(k_{x},k_{y})=(\pi,\pm\pi/2)$. In particular, when $(\lambda,\kappa)=(0,0)$,
the band degenerate points become two degenerate lines at $k_{y}=\pm\pi/2$.

The 1D Zak phase is defined as

$\Theta_{x}\left(k_{y}\right)=\int\left\langle\psi(k_{x})\right|i\partial_{k_{x}}\left|\psi(k_{x})\right\rangle
dk_{x},$ (14)

which reflects the winding of the Bloch Hamiltonian when $k_{x}$ varies
throughout the entire BZ of the decoupled 1D lattice. The 1D Zak phase
$\Theta_{x}(\pm\pi/2)$ is quantized to $\pi$ in the topologically nontrivial
phases (the purple and cyan regions in Fig. 2), whereas $\Theta_{x}(\pm\pi/2)$
is quantized to $0$ in the trivial phase (the green and yellow regions in Fig.
2). The topological characterization is shown in Fig. 6 and is consistent with
the phase diagram in Fig. 2. $\Theta_{x}(\pm\pi/2)=\pi$ characterizes the
topological metal and topological insulator. $\Theta_{x}(\pm\pi/2)=0$
characterizes the ordinary metal and ordinary insulator.

##  V Discussion

The vertical coupling $\lambda\neq\pm\kappa$ is crucial for the topological
metallic phase as indicated in Fig. 2. In addition, we consider the horizontal
couplings $\lambda$ and $\kappa$ are replaced by $\lambda^{\prime}$ and
$\kappa^{\prime}$. Then, when $\lambda^{\prime}=\kappa^{\prime}$, the energy
band must have degenerate points located at $(k_{x},k_{y})=(\pm\pi,\pm\pi/2)$;
when $\lambda^{\prime}=-\kappa^{\prime}$, the energy band must have degenerate
points located at $(k_{x},k_{y})=(0,\pm\pi/2)$. Notably,
$\lambda^{\prime}\neq\pm\kappa^{\prime}$ is necessary for the topological
metallic phase, but $\lambda=\lambda^{\prime}$ and $\kappa=\kappa^{\prime}$
are not necessary. The nearest neighbor coupling $m\neq 0$ is important for
the different edge states in the topological metallic phases. For $m=0$, only
the in-gap edge states exist in the topological phase and only one type of
topological phase presents in this situation.

The spin-orbit coupling is used to create the magnetic flux and the
traditional topological metal. However, the proposed 2D topological metal is
time-reversal symmetric without the magnetic flux. This requires the spin-
orbit coupling is negligible in the material. The organic materials composed
of light elements have weak spin-orbit coupling LJin15 . Therefore, the
organic materials may be suitable candidates for the class AI topological
phase ZSLiao21 .

If the time-reversal symmetry breaks when the square lattice is threaded by
the magnetic flux. The topological metallic phase still exists through tuning
the magnetic flux and the coupling strengths. In this situation, the
topological phase is characterized by the nonzero Chern number. Under the PBC
in the $y$ direction, a pair of anti-chiral edge states exist after imposing
the OBC in the $x$ direction. The antichiral edge states also connect the
upper and lower energy bands to form the conductive surface states, but the
two antichiral edge states on different boundaries propagate in the same
direction in contrast to the chiral edge states EColo18 .

Furthermore, the non-Hermiticity compresses the band gap between the valence
band and the conduction band. The non-Hermitian topological metallic phase
exits in the proposed 2D square lattice when the additional gain and loss are
presented in different sublattices. This is an interesting topic that
deserving further investigation.

##  VI Conclusion

In conclusion, we provide an alternative way of generating the topological
metals. The proposed topological metallic phases are protected by the time-
reversal symmetry and are free from magnetic flux. The time-reversal symmetry
ensures a zero Chern number; however, the 2D square lattice supports
topologically nontrivial phases, which are distinguished by the different
types of edge states. The 2D square lattice supports three different cases in
the topological metallic phase that significantly differ from the traditional
topological metallic phases created through tuning the magnetic flux. These
different cases can be deformed into each other without topological phase
transition. Our findings provide insights into the time-reversal symmetric
topological metals.

## Acknowledgement

This work was supported by the National Natural Science Foundation of China
(Grants No. 11975128 and No. 11874225).

## References

  * (1) L. Fu, C. L. Kane, and E. J. Mele, Topological Insulators in Three Dimensions, Phys. Rev. Lett. 98, 106803 (2007).
  * (2) M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * (3) X.-G. Wen, Colloquium: Zoo of quantum-topological phases of matter, Rev. Mod. Phys. 89, 041004 (2017).
  * (4) A. Bansil, H. Lin, and T. Das, Colloquium: Topological band theory, Rev. Mod. Phys. 88, 021004 (2016).
  * (5) X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * (6) A. A. Burkov and L. Balents, Weyl semimetal in a topological insulator multilayer, Phys. Rev. Lett. 107, 127205 (2011).
  * (7) S. M. Young, S. Zaheer, J. C. Y. Teo, C. L. Kane, E. J. Mele, and A. M. Rappe, Dirac Semimetal in Three Dimensions, Phys. Rev. Lett. 108, 140405 (2012).
  * (8) N. P. Armitage, E. J. Mele, and A. Vishwanath, Weyl and Dirac semimetals in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018).
  * (9) A. A. Burkov, Topological semimetals, Nat. Mater. 15, 1145 (2016).
  * (10) B. Yan and C. Felser, Topological materials: Weyl semimetals, Ann. Rev. Cond. Mat. Phys. 8, 337 (2017).
  * (11) C. K. Chiu, J. C. Teo, A. P. Schnyder, and S. Ryu, Classification of topological quantum matter with symmetries, Rev. Mod. Phys. 88, 035005 (2016).
  * (12) A. Bansil, H. Lin, and T. Das, Topological band theory, Rev. Mod. Phys. 88, 021004 (2016).
  * (13) L. Jin and Z. Song, Bulk-boundary correspondence in a non-Hermitian system in one dimension with chiral inversion symmetry, Phys. Rev. B 99, 081103(R) (2019).
  * (14) K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and Topology in Non-Hermitian Physics, Phys. Rev. X 9, 041015 (2019).
  * (15) Y. Ashida, Z. Gong, M. Ueda, Non-Hermitian Physics, Adv. Phys. 69, 249 (2020).
  * (16) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Quantized Hall Conductance in a Two-Dimensional Periodic Potential, Phys. Rev. Lett. 49, 405 (1982).
  * (17) D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010).
  * (18) J. C. Y. Teo and C. L. Kane, Topological Defects and Gapless Modes in Insulators and Superconductors, Phys. Rev. B 82, 115120 (2010).
  * (19) Y. Ran, Y. Zhang, and A. Vishwanath, One-dimensional topologically protected modes in topological insulators with lattice dislocations. Nat. Phys. 5, 298 (2009).
  * (20) L. J. Lang, X. Cai, and S. Chen, Edge states and topological phases in one-dimensional optical superlattices, Phys. Rev. Lett. 108, 220401 (2012).
  * (21) X.-J. Liu, K. T. Law, T. K. Ng, and P. A. Lee, Detecting Topological Phases in Cold Atoms, Phys. Rev. Lett. 111, 120402 (2013).
  * (22) M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Realization of the Hofstadter Hamiltonian with Ultracold Atoms in Optical Lattices, Phys. Rev. Lett. 111, 185301 (2013).
  * (23) L. Lu, J. D. Joannopoulos, and M. Soljačić, Topological photonics, Nat. Photon. 8, 821 (2014).
  * (24) S.-Y. Yu, C. He, Z. Wang, F.-K. Liu, X.-C. Sun, Z. Li, H.-Z. Lu, M.-H. Lu, X.-P. Liu, and Y.-F. Chen, Elastic pseudospin transport for integratable topological phononic circuits, Nat. Commun. 9, 3072 (2018).
  * (25) Y.-B. Yang, L.-M. Duan, and Y. Xu, Continuously tunable topological pump in high-dimensional cold atomic gases, Phys. Rev. B 98, 165128 (2018).
  * (26) D.-W. Zhang, Y.-Q. Zhu, Y. X. Zhao, H. Yan, and S.-L. Zhu, Topological quantum matter with cold atoms, Adv. Phys. 67, 253 (2019).
  * (27) X.-D. Chen, W.-M. Deng, F.-L. Shi, F.-L. Zhao, M. Chen, and J.-W. Dong, Direct Observation of Corner States in Second-Order Topological Photonic Crystal Slabs, Phys. Rev. Lett. 122, 233902 (2019).
  * (28) Y. Xu, Topological gapless matters in three-dimensional ultracold atomic gases, Front. Phys. 14, 43402 (2019).
  * (29) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, and R. Thomale, Topolectrical-circuit realization of topological corner modes, Nat. Phys. 14, 925 (2018).
  * (30) T. Helbig, T. Hofmann, C. H. Lee, R. Thomale, S. Imhof, L. W. Molenkamp, and T. Kiessling, Band structure engineering and reconstruction in electric circuit networks, Phys. Rev. B 99, 161114 (2019).
  * (31) Y. Lu, N. Jia, L. Su, C. Owens, G. Juzeliunas, D. I. Schuster, and J. Simon, Probing the Berry curvature and Fermi arcs of a Weyl circuit, Phys. Rev. B 99, 020302 (2019).
  * (32) S. Murakami, N. Nagaosa, and S. C. Zhang, Spin-hall insulator, Phys. Rev. Lett. 93, 156804 (2004).
  * (33) C. L. Kane and E. J. Mele, Quantum Spin Hall Effect in Graphene, Phys. Rev. Lett. 95, 226801 (2005).
  * (34) B. A. Bernevig, T. L. Hughes, and S. C. Zhang, Quantum spin Hall effect and topological phase transition in HgTe quantum wells, Science, 314, 1757 (2006).
  * (35) X.-L. Qi, Y.-S. Wu, and S.-C. Zhang, Topological quantization of the spin hall effect in two-dimensional paramagnetic semiconductors, Phys. Rev. B 74, 085308 (2006).
  * (36) C. L. Kane and E. J. Mele, $\mathbb{Z}_{2}$ Topological Order and the Quantum Spin Hall Effect, Phys. Rev. Lett. 95, 146802 (2005).
  * (37) A. Roth, C. Brüne, H. Buhmann, L. W. Molenkamp, J. Maciejko, X. L. Qi, and S. C. Zhang, Nonlocal transport in the quantum spin Hall state, Science 325, 294 (2009).
  * (38) L. Sheng, D. N. Sheng, C. S. Ting, and F. D. M. Haldane, Nondissipative spin Hall effect via quantized edge transport, Phys. Rev. Lett. 95, 136602 (2005).
  * (39) E. Prodan, Robustness of the spin-Chern number, Phys. Rev. B 80, 125327 (2009).
  * (40) Y. Yang, Z. Xu, L. Sheng, B. Wang, D. Y. Xing, and D. N. Sheng, Time-reversal-symmetry-broken quantum spin Hall effect, Phys. Rev. Lett. 107, 066602 (2011).
  * (41) F. Zhang, J. Jung, G. A. Fiete, Q. Niu, and A. H. MacDonald, Spontaneous quantum Hall states in chirally stacked few-layer graphene systems, Phys. Rev. Lett. 106, 156801 (2011).
  * (42) F. Zhang, A. H. MacDonald, and E. J. Mele, Valley Chern numbers and boundary modes in gapped bilayer graphene, Proc. Natl. Acad. Sci. USA 110, 10546 (2013).
  * (43) M. Ezawa, Topological Kirchhoff law and bulk-edge correspondence for valley Chern and spin-valley Chern numbers, Phys. Rev. B 88, 161406(R) (2013).
  * (44) R. Resta, Macroscopic polarization in crystalline dielectrics: the geometric phase approach, Rev. Mod. Phys. 66, 899 (1994).
  * (45) P. Delplace, D. Ullmo, and G. Montambaux, The Zak phase and the existence of edge states in graphene, Phys. Rev. B 84, 195452 (2011).
  * (46) F. Liu and K. Wakabayashi, Novel Topological Phase with a Zero Berry Curvature, Phys. Rev. Lett. 118, 076803 (2017).
  * (47) S. Murakami, Phase transition between the quantum spin Hall and insulator phases in 3D: Emergence of a topological gapless phase, New J. Phys. 9, 356 (2007).
  * (48) X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Topological semimetal and Fermi-arc surface states in the electronic structure of pyrochloreiridates, Phys. Rev. B 83, 205101 (2011).
  * (49) Z. Wang, H. Weng, Q. Wu, X. Dai, and Z. Fang, Three-dimensional Dirac semimetal and quantum transport in Cd3As2, Phys. Rev. B 88, 125427 (2013).
  * (50) Z. Zhu, G. W. Winkler, Q. Wu, J. Li, and A. A. Soluyanov, Triple point topological metals, Phys. Rev. X 6, 031003 (2016).
  * (51) B. J. Wieder, Y. Kim, A. Rappe, and C. Kane, Double Dirac semimetals in three dimensions, Phys. Rev. Lett. 116, 186402 (2016).
  * (52) H. Weng, C. Fang, Z. Fang, and X. Dai, Topological semimetals with triply degenerate nodal points in $\theta$-phase tantalum nitride, Phys. Rev. B 93, 241202 (2016).
  * (53) S. Singh, Q. Wu, C. Yue, A. H. Romero, and A. A. Soluyanov, Topological phonons and thermoelectricity in triple-point metals, Phys. Rev. Mater. 2, 114204 (2018).
  * (54) G. Shan and H. B. Gao, New topological semimetal candidate of nonsymmorphic PdSb2 with unique six-fold degenerate point, Front. Phys. 14, 43201 (2019).
  * (55) R. Chapai, Y. Jia, W. Shelton, R. Nepal, M. Saghayezhian, J. DiTusa, E. Plummer, C. Jin, and R. Jin, Fermions and bosons in nonsymmorphic PdSb2 with six-fold degeneracy, Phys. Rev. B 99, 161110 (2019).
  * (56) A. A. Burkov, M. D. Hook, and L. Balents, Topological nodal semimetals, Phys. Rev. B 84, 235126 (2011).
  * (57) R. Yu, H. M. Weng, Z. Fang, X. Dai, and X. Hu, Topological Node-line Semimetal and Dirac Semimetal State in Antiperovskite Cu3PdN, Phys. Rev. Lett. 115, 036807 (2015).
  * (58) H. M. Weng, Y. Y Liang, Q. N. Xu, R. Yu, Z. Fang, X. Dai, and Y. Kawazoe, Topological node-line semimetal in three-dimensional graphene networks, Phys. Rev. B 92, 045108 (2015).
  * (59) W. Wu, Y. Liu, S. Li, C. Zhong, Z.-M. Yu, X.-L. Sheng, Y. X. Zhao, and S. A. Yang, Nodal surface semimetals: Theory and material realization, Phys. Rev. B 97, 115125 (2018).
  * (60) A. A. Soluyanov, D. Gresch, Z. Wang, Q. Wu, M. Troyer, X. Dai, and B. A. Bernevig, Type-II Weyl semimetals, Nature 527, 495 (2015).
  * (61) Y. Xu, F. Zhang, and C. Zhang, Structured Weyl points in spin-orbit coupled fermionic superfluids, Phys. Rev. Lett. 115, 265304 (2015).
  * (62) S. Li, Z. M. Yu, Y. Liu, S. Guan, S. S. Wang, X. Zhang, Y. Yao, and S. A. Yang, Type-II nodal loops: Theory and material realization, Phys. Rev. B 96, 081106 (2017).
  * (63) M. Bahari and M. V. Hosseini, One-dimensional topological metal, Phys. Rev. B 99, 155128 (2019).
  * (64) S.-Y. Xu, C. Liu, S. K. Kushwaha, R. Sankar, J. W. Krizan, I. Belopolski, M. Neupane, G. Bian, N. Alidoust, T.-R. Chang, et al., Observation of Fermi arc surface states in a topological metal: A New Type of 2D Electron Gas, Science 347, 294 (2015).
  * (65) Z. J. Wang, D. Gresch, A. A. Soluyanov, et al., MoTe2: A Type-II Weyl Topological Metal, Phys. Rev. Lett. 117, 056805 (2016).
  * (66) X. Z. Ying and A. Kamenev, Symmetry-Protected Topological Metals, Phys. Rev. Lett. 121, 086810 (2018).
  * (67) M. Aidelsburger, M. Lohse, C. Schweizer, M. Atala, J. T. Barreiro, S. Nascimbène, N. R. Cooper, I. Bloch, and N. Goldman, Measuring the Chern number of Hofstadter bands with ultracold bosonic atom, Nat. Phys. 11, 162 (2015).
  * (68) S. Mittal, E. A. Goldschmidt, and M. Hafezi, A topological source of quantum light, Nature 561, 502 (2018).
  * (69) T. Hofmann, T. Helbig, C. H. Lee, M. Greiter, and R. Thomale, Chiral voltage propagation and calibration in a topolectrical Chern circuit, Phys. Rev. Lett. 122, 247702 (2019).
  * (70) Z. J. Yang, F. Gao, X. H. Shi, X. Lin, Z. Gao, Y. D. Chong, and B. L. Zhang, Topological acoustics, Phys. Rev. Lett. 114, 114301 (2015).
  * (71) C. He, X. Ni, H. Ge, X.-C. Sun, Y.-B. Chen, M.-H. Lu, X.-P. Liu, and Y.-F. Chen, Acoustic topological insulator and robust one-way sound transport, Nat. Phys. 12, 1124 (2016).
  * (72) Y. G. Peng, C. Z. Qin, D. G. Zhao, Y. X. Shen, X. Y. Xu, M. Bao, H. Jia, and X. F. Zhu, Experimental demonstration of anomalous Floquet topological insulator for sound, Nat. Commun. 7, 13368 (2016).
  * (73) J. Y. Lu, C. Y. Qiu, W. Y. Deng, X. Q. Huang, F. Li, F. Zhang, S. Q. Chen, and Z. Y. Liu, Valley topological phases in bilayer sonic crystals. Phys. Rev. Lett. 120, 116802 (2018).
  * (74) L. Jin and Z. Song, Symmetry-Protected Scattering in Non-Hermitian Linear Systems, Chin. Phys. Lett. 38, 024202 (2021).
  * (75) H. C. Wu, L. Jin, and Z. Song, Nontrivial topological phase with a zero Chern number, Phys. Rev. B 102, 035145 (2020).
  * (76) S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, Topological insulators and superconductors: tenfold way and dimensional hierarchy, New J. Phys, 12, 065010 (2010).
  * (77) X. Wen, C. Qiu, Y. Qi, L. Ye, M. Ke, F. Zhang, and Z. Liu, Acoustic Landau quantization and quantum-Hall-like edge states, Nat. Phys. 15, 352 (2019).
  * (78) N. H. Le, A. J. Fisher, N. J. Curson, and E. Ginossar, Topological phases of a dimerized fermi–hubbard model for semiconductor nano-lattices, npj Quantum Inf. 6, 24 (2020).
  * (79) L. Jin, M. Pfender, N. Aslam, P. Neumann, S. Yang, J. Wrachtrup, and R.-B. Liu, Proposal for a room-temperature diamond maser, Nat. Commun. 6, 8251 (2015).
  * (80) Z. S. Liao, H. H. Zhang, and Z. B. Yan, Nonlinear Hall effect in two-dimensional class AI metals, arXiv:2104.08477 (2021).
  * (81) E. Colomés and M. Franz, Antichiral edge states in a modified Haldane nanoribbon, Phys. Rev. Lett. 120, 086603 (2018).