# Mu-Metal Enhancement of Effects in Electromagnetic Fields Over Single
Emitters Near Topological Insulators

Eitan Dvorqueza, Benjamín Paveza, Qiang Sunb, Felipe Pintoa, Andrew D.
Greentreeb, Brant C. Gibsonb, and Jerónimo R. Mazea aInstitute of Physics,
Pontificia Universidad Católica de Chile, Santiago, Chile.  bARC Centre of
Excellence for Nanoscale Biophotonics, RMIT University, Melbourne, VIC 3001,
Australia.

###### Abstract

We focus on the transmission and reflection coefficients of light in systems
involving of topological insulators (TI). Due to the electro-magnetic coupling
in TIs, new mixing coefficients emerge leading to new components of the
electromagnetic fields of propagating waves. We have discovered a simple
heterostructure that provides a 100-fold enhancement of the mixing
coefficients for TI materials. Such effect increases with the TI’s wave
impedance. We also predict a transverse deviation of the Poynting vector due
to these mixed coefficients contributing to the radiative electromagnetic
field of an electric dipole. Given an optimal configuration of the dipole-TI
system, this deviation could amount to $0.18\%$ of the Poynting vector due to
emission near not topological materials, making this effect detectable.

##  I Introduction

Topological insulators (TIs) are materials with non-trivial topological order.
These materials present an insulating bulk, while still having conductive edge
states on their surface Qi and Zhang (2011); Hasan and Moore (2011); Ando
(2013). These properties allow interesting electronic effects such as the
quantum spin Hall effect in 2D TIs, which is a quantized Hall current that is
present despite the absence of an external magnetic field.

TIs exhibit unusual optical properties, such as the coupling between the
electric and magnetic fields, which produces magnetic monopole-like behavior
when a point charge is near a TI’s surface Qi  _et al._ (2009), varying Hall
conductivity Essin  _et al._ (2009), and the magneto-electric effectMartín-
Ruiz  _et al._ (2019). Remarkably, the magneto-electric effect that TIs
exhibits has been measured to respond faster than multiferroic’s magneto-
electric effectKitagawa  _et al._ (2010); Velev  _et al._ (2011) and does so
without material fatigueMoore (2010).

In recent years, these properties have been at the center of the development
of topological devices such as nanowiresLegg  _et al._ (2022); Fischer  _et
al._ (2022); Münning  _et al._ (2021), memory storageWu  _et al._ (2021),
nanoparticlesCastro-Enríquez  _et al._ (2022); Lin  _et al._ (2015), among
others. Application of these devices have been proposed for optoelectronics,
spintronics and quantum computingHe  _et al._ (2019); Pandey  _et al._ (2021);
Breunig and Ando (2021). To make such devices, it is crucial to develop
reliable methods to characterize TIs.

Efforts to characterize the optical effects exhibited by TIs usually seek to
measure the Faraday rotation of a transmitted ray, as it should be quantized
by the fine structure constant $\alpha$Dziom  _et al._ (2017); C and F (2017).
However, this effect is usually very small (1-10 mrad), and the conditions for
it to be purely topological are such that it requires a very strong magnetic
field (on the Tesla range), low temperatures, and low-frequency emissions (THz
range) Crosse (2017); Shuvaev  _et al._ (2022); Tse and MacDonald (2010,
2011). The Faraday effect can be enhanced by a wide range of material
configurations. Composite TI multi-layered media materialsCrosse (2016) and
defects on the TI’s surfaceArdakani and Zare (2021) are examples of such
systems. The geometry of the system can also have an impact on the enhancement
of these effectsMartín-Ruiz  _et al._ (2016); Ge  _et al._ (2014); Gangaraj
_et al._ (2020). It is also possible to observe the rotation on semi-magnetic
TIs, without a magnetic fieldMogi  _et al._ (2022).

Here we show an optical method to characterize TIs. We want to find new, fast,
and spatially resolved techniques to characterize TIs optically. With this in
mind, we will start by setting the foundations for the optical
characterization of a TI’s planar surface via the Axion Lagrangian, remarking
on the fundamental aspects of the derivation for new boundary conditions,
Fresnel coefficients, and the electric field’s Green’s function. Then, we will
study the inclusion of a Mu-Metal layer into the system to enhance these
topological effects, in a setup like the one in figure 1.

Figure 1: Sketch of the studied dipole-TI-Mu metal system for an $\hat{x}$
oriented dipole $\vec{d}$, where $z_{p}$ is the distance of the dipole
relative to the TI’s surface, and $z_{l}$ is the distance from the surface to
the observation plane.

##  II Theoretical Background

###  1 Axion Formalism for Topological Insulators

To describe the effect of an interface over any incident emission, we start by
obtaining the Maxwell equations and boundary conditions.

Due to this electromagnetic coupling, we can describe these materials via
axion theoryMartín-Ruiz  _et al._ (2021), adding an extra term to the usual
Maxwell Lagrangian, such that
$\mathcal{L}=\mathcal{L}_{M}+\mathcal{L}_{\Theta}$, where

$\displaystyle\mathcal{L}_{M}=$
$\displaystyle\frac{1}{8\pi}\left(\epsilon|\bf{E}|^{2}-\frac{1}{\mu}|\bf{B}|^{2}\right),$
(1) $\displaystyle\mathcal{L}_{\Theta}=$
$\displaystyle\frac{\alpha}{4\pi^{2}}\frac{\Theta(\mathbf{r})}{\mu_{0}c^{2}}\mathbf{E}\cdot\mathbf{B},$
(2)

where $\bf{E}$ and $\bf{B}$ are the electric and magnetic fields,
respectively, $\mu_{0}$ is the vacuum’s magnetic permeability, and $c$ is the
speed of light. $\alpha$ is the fine structure constant and $\Theta(\bf{r})$
is the axionic coupling parameter as a function of position in space $\bf{r}$,
which is 0 for a non-topological material and an odd multiple of $\pi$ for
TIsQi and Zhang (2011). The Euler-Lagrange equation outputs a new set of
modified time-dependent Maxwell equations;

$\displaystyle\nabla\cdot\mathbf{E}=\mu_{0}c^{2}\rho-c\frac{\alpha}{4\pi^{2}}\nabla\Theta\cdot\mathbf{B},$
(3)
$\displaystyle\nabla\times\textbf{B}=\mu_{0}\textbf{J}+\frac{1}{c^{2}}\frac{\partial\textbf{E}}{\partial
t}+\frac{1}{c}\frac{\alpha}{4\pi^{2}}\nabla\Theta\times\mathbf{E}+\frac{1}{c}\frac{\partial\Theta}{\partial
t}\mathbf{B},$ (4)
$\displaystyle\nabla\times\textbf{E}=-\tfrac{\partial\textbf{B}}{\partial t},$
(5) $\displaystyle\nabla\cdot\textbf{B}=0,$ (6)

and, therefore, new terms on the boundary conditions to the half-plane space
problem of an interface Air-TI, are directly dependent on the difference
between the axion parameters of the two materials.

Taking advantage of the linear relations between $\mathbf{E}$ and $\mathbf{B}$
upon converting the time $t$ to frequency $\omega$ via Fourier Transform, we
define the $\mathbf{D}$ and $\mathbf{H}$ fields as

$\displaystyle\mathbf{D}(\mathbf{r},\omega)$
$\displaystyle=\epsilon_{0}\epsilon_{\mathrm{r}}(\mathbf{r},\omega)\mathbf{E}(\mathbf{r},\omega)+\frac{\alpha}{\pi}\frac{\Theta(\mathbf{r})}{\mu_{0}c}\mathbf{B}(\mathbf{r},\omega),$
(7) $\displaystyle\mathbf{H}(\mathbf{r},\omega)$
$\displaystyle=\frac{1}{\mu_{0}\mu_{\mathrm{r}}(\mathbf{r},\omega)}\mathbf{B}(\mathbf{r},\omega)-\frac{\alpha}{\pi}\frac{\Theta(\mathbf{r})}{\mu_{0}c}\mathbf{E}(\mathbf{r},\omega).$
(8)

where $\epsilon_{0}$ is the vacuum’s electric permittivity, and
$\epsilon_{\mathrm{r}}(\mathbf{r},\omega)$ and
$\mu_{\mathrm{r}}(\mathbf{r},\omega)$ are the relative electric permittivity
and magnetic permeability of the TI, respectively. With these definitions, we
can retain the well-known boundary equations
$\hat{z}\times\mathbf{E}_{1}=\hat{z}\times\mathbf{E}_{2}$ and
$\hat{z}\times\mathbf{H}_{1}=\hat{z}\times\mathbf{H}_{2}$, where $\hat{z}$ is
the direction on the axis perpendicular to the TI’s plane and the subscripts
$1$ and $2$ indicate fields over and under the TI’s surface, respectively.
Also, $\mathbf{D}$ and $\mathbf{H}$ follow the usual Maxwell equations.

The new Helmholtz wave equation is obtained by combining the modified Maxwell
equations:

$\textstyle{}\nabla\times\nabla\times\mathbf{E}-\frac{\omega^{2}}{c^{2}}\mathbf{E}-i\omega\mu_{0}\left(\frac{\alpha}{4\pi^{2}\mu_{0}c}\right)\nabla\Theta\times\mathbf{E}=-i\omega\mu_{0}\mathbf{J},$
(9)

where $\mathbf{J}$ stands for the usual current, constituted by the external
free current and the current generated by the non-topological
magnetodielectric responses of the material.

The extra term in equation (9) can be understood as an extra superficial
current due to the topological properties of the material, given by the
coupling parameter $\Theta(\mathbf{r})$. This current is exclusively located
on the surface of the interface due to the step function like the behavior of
$\Theta(\mathbf{r})$, which implies that $\nabla\Theta(\mathbf{r})$ behaves as
a Dirac delta on the interface of the materials.

###  2  Two-Layered Half-Plane System

To characterize the reflected or transmitted waves, we explore the Fresnel
coefficients that follow these new boundary conditions. Comparing these new
coefficients with the regular ones for a non-topological material gives us
insight into the effect we seek to measure. We start by considering a planar
wave ansatz.

When considering an incident planar wave, we must take into account its
polarization; Transverse Electric (TE) or Transverse Magnetic (TM), referring
to which field, $\mathbf{E}$ or $\mathbf{H}$, is initially pointing in the
$\hat{y}$ axis. For magnetodielectric media, reflected and transmitted waves
maintain the initial axis of their respective fields. On the other hand, due
to the electromagnetic coupling, the waves interacting with TIs present new
fields components on all axes. For example, given a TE polarized wave, with
$\mathbf{E}=E_{0}\hat{y}$, we should expect the reflected wave to present a
resulting electric field component not only on the $\hat{y}$ axis, but on the
$\hat{x}$ and $\hat{z}$ axis as well.

Using a planar wave ansatz for the electric and magnetic field of initially TE
and TM polarized waves, and combining them with the new boundary conditions
previously discussed yields the following new reflection and transmission
matricesCrosse  _et al._ (2015),

$\displaystyle\mathbf{R}$
$\displaystyle=\begin{pmatrix}R_{\textit{TE,TE}}&R_{\textit{TE,TM}}\\\
R_{\textit{TM,TE}}&R_{\textit{TM,TM}}\end{pmatrix}=\frac{1}{\gamma}\left[\begin{matrix}(\mu_{2}k_{z1}-\mu_{1}k_{z2})\Omega_{\epsilon}-k_{z1}k_{z2}\Delta^{2}&-2\mu_{2}n_{1}k_{z1}k_{z2}\Delta\\\
-2\mu_{2}n_{1}k_{z1}k_{z2}\Delta&(\epsilon_{2}k_{z1}-\epsilon_{1}k_{z2})\Omega_{\mu}-k_{z1}k_{z2}\Delta^{2}\end{matrix}\right],$
(10) $\displaystyle\mathbf{T}$
$\displaystyle=\begin{pmatrix}T_{\textit{TE,TE}}&T_{\textit{TE,TM}}\\\
T_{\textit{TM,TE}}&T_{\textit{TM,TM}}\end{pmatrix}=\frac{2k_{z1}}{\gamma}\left[\begin{matrix}\mu_{2}\Omega_{\epsilon}&-\mu_{2}n_{1}k_{z2}\Delta\\\
\mu_{2}n_{2}k_{z1}\Delta&\tfrac{n_{2}}{n_{1}}\epsilon_{1}\Omega_{\mu}\end{matrix}\right],$
(11)

where $\epsilon_{i}$ and $\mu_{i}$ are the relative electric permittivity and
magnetic permeability, respectively, of the $i$-layer.
$\mathbf{k}_{i}=(k_{xi},k_{yi},k_{zi})$ is the wave vector in each layer,
$n_{i}=\sqrt{\mu_{i}\epsilon_{i}}$ is the refractive index of the
corresponding layer, $\Delta=\alpha\mu_{1}\mu_{2}(\Theta_{2}-\Theta_{1})/\pi$,
$\Omega_{\mu}=\mu_{1}\mu_{2}(k_{z1}\mu_{2}+k_{z2}\mu_{1})$,
$\Omega_{\epsilon}=\mu_{1}\mu_{2}(k_{z1}\epsilon_{2}+k_{z2}\epsilon_{1})$ and
$\gamma=\frac{\Omega_{\mu}\Omega_{\epsilon}}{\mu_{1}\mu_{2}}+k_{z1}k_{z2}\Delta^{2}$.

It is worth remarking that, when $\Theta_{2}-\Theta_{1}\rightarrow 0$ the
mixed coupling coefficients $R_{\textit{TM,TE}}$, $R_{\textit{TE,TM}}$,
$T_{\textit{TE,TM}}$, $T_{\textit{TM,TE}}\rightarrow 0$ and
$R_{\textit{TE,TE}}$, $R_{\textit{TM,TM}}$, $T_{\textit{TE,TE}}$, and
$T_{\textit{TM,TM}}$ reduce to the conventional coefficients for
magnetodielectric media Chew (1995); Novotny and Hecht (2012). Further limit
cases are found in table 1.

|  | Air/Perfect
---
Conductor
Air/Mu-Metal | Mu-Metal/Air |  | TI/Perfect
---
Conductor
TI/Mu-Metal | Mu-Metal/TI
| Reflection
---
Matrix
$\left[\begin{matrix}-1&0\\\ 0&1\end{matrix}\right]$ | $\left[\begin{matrix}1&0\\\ 0&-1\end{matrix}\right]$ | $\left[\begin{matrix}-1&0\\\ 0&1\end{matrix}\right]$ | $\left[\begin{matrix}-1&0\\\ 0&1\end{matrix}\right]$ | $\left[\begin{matrix}1&0\\\ 0&-1\end{matrix}\right]$ | $\left[\begin{matrix}-1&0\\\ 0&1\end{matrix}\right]$
| Transmission
---
Matrix
$\left[\begin{matrix}0&0\\\ 0&0\end{matrix}\right]$ | $\left[\begin{matrix}2&0\\\ 0&2\frac{\cos{\phi_{r}}}{\cos{\phi_{t}}}\end{matrix}\right]$ | $\left[\begin{matrix}0&0\\\ 0&0\end{matrix}\right]$ | $\left[\begin{matrix}0&0\\\ 0&0\end{matrix}\right]$ | $\frac{\displaystyle 1}{\displaystyle\epsilon_{1}+\tilde{\alpha}^{2}\mu_{1}}\left[\begin{matrix}2\epsilon_{1}&\frac{n_{2}k_{z1}}{k_{z2}}2\tilde{\alpha}\\\ -2n_{1}\tilde{\alpha}&\frac{n_{2}k_{z1}}{n_{1}k_{z2}}2\epsilon_{1}\end{matrix}\right]$ | $\left[\begin{matrix}0&0\\\ 0&0\end{matrix}\right]$
Table 1: Limit case results of the Fresnel matrices for different 2 layered-
Systems involving Perfect Conductors ($\epsilon\rightarrow\infty$) and Mu-
Metals ($\mu\rightarrow\infty$). The perfect conductor retains its reflecting
behavior when paired with TIs. The Mu-metal and TI pairing retains the Mu-
Metal shielding phenomena, irrespective of the new off-diagonal coefficients
in the TI/Mu-Metal configuration, given that the transmission coefficients of
the Mu-Metal/TI configuration are null.

The Fresnel coefficients are given by

$r_{ij}=\left|R_{ij}\right|^{2},\qquad
t_{ij}=\frac{k_{z2}}{k_{z1}}\frac{\mu_{1}}{\mu_{2}}\left|T_{ij}\right|^{2}.$
(12)

New mixed Fresnel coefficients, derived from the new off-diagonal terms on the
$\mathbf{R}$ and $\mathbf{T}$ matrices, are about $10^{-5}$ times smaller
compared to the usual non-mixed onesCrosse  _et al._ (2015) because they
scaled as $\alpha^{2}$. This indicates that any effects from these
coefficients are very small with respect to the usual fields.

The reflected mixed coefficients tend to increase with the TI’s impedance
relative to the air. However, since the impedance is an intrinsic property of
the material and the increase is relatively minor, we will search for another
way to maximize these coefficients. Topologically non-trivial heterostructures
have been widely proposed to reach said goal, primarily because proximity
coupling introduces magnetism to the TI without modifying the structure of the
material, both enabling and enhancing topological responseLiu and Hesjedal
(2021); Mogi  _et al._ (2017); An  _et al._ (2020); Zuo  _et al._ (2013);
Crosse (2016). Similarly, we propose the addition of a single Mu-Metal layer
below the TI.

###  3 Electric Green’s Function

We define the Green’s function such that it satisfies

$\mathbf{E}(\mathbf{r},\omega)=i\omega\mu_{0}\int
d^{3}r^{\prime}\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)\cdot\mathbf{J}(\mathbf{r}^{\prime}),$
(13)

then, knowing the behavior of $\mathbf{G}$ gives us the behavior of
$\mathbf{E}$. The Green’s function also satisfies the following Helmholtz
equation,

$\nabla\times\frac{1}{\mu_{\mathrm{r}}(\mathbf{r},\omega)}\nabla\times\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)-\frac{\omega^{2}}{c^{2}}\epsilon_{\mathrm{r}}(\mathbf{r},\omega)\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)=\delta(\mathbf{r}-\mathbf{r}^{\prime}).$
(14)

This equation is equal to the one for a traditional magnetodielectric, as we
are interested in the Green’s function on the bulk of the upper layer and not
on the interface.

The corresponding Green’s function, in its singularity extracted form, for
non-diagonal Fresnel matrices, isCrosse  _et al._ (2015)

$\displaystyle\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)=$
$\displaystyle\frac{i}{8\pi^{2}}\int
d^{2}k_{p}\,\mu_{\mathrm{r}}(\omega)\times$
$\displaystyle\left[\frac{\mathbf{C}(\mathbf{k},\mathbf{k}^{\prime}):\mathbf{F}_{z\lessgtr
z^{\prime}}(z,z^{\prime})}{k_{z}k^{2}_{p}}\right]e^{i\mathbf{k}_{p}\cdot(\mathbf{r}_{p}-\mathbf{r}^{\prime}_{p})},$
(15)

where $\mu_{\mathrm{r}}(\omega)$ is the magnetic permeability in the layer
with the emitting source, $k^{2}=k_{z}^{2}+k_{p}^{2}$, the operation $:$
corresponds to the Frobenius inner product, $\mathbf{F}_{z\lessgtr
z^{\prime}}(z,z^{\prime})$ corresponds to a matrix that describes the waves
traveling upward or downward in a given layer $n$, and

$\mathbf{C}(\mathbf{k},\mathbf{k}^{\prime})=\begin{pmatrix}\mathbf{m}^{*}(\mathbf{k})\otimes\mathbf{m}(\mathbf{k^{\prime}})&\mathbf{n}^{*}(\mathbf{k})\otimes\mathbf{m}(\mathbf{k^{\prime}})\\\
\mathbf{m}^{*}(\mathbf{k})\otimes\mathbf{n}(\mathbf{k^{\prime}})&\mathbf{n}^{*}(\mathbf{k})\otimes\mathbf{n}(\mathbf{k^{\prime}})\end{pmatrix},$
(16)

where $\mathbf{m}(\mathbf{r})=i\nabla_{\mathbf{r}}\times\hat{z}$ and
$\mathbf{n}(\mathbf{r})=\frac{1}{k}\nabla_{\mathbf{r}}\times\nabla_{\mathbf{r}}\times\hat{z}$
are the dyadic operators that generate the s and p polarization vectors
Novotny and Hecht (2012) and $\otimes$ is the outer product.

For only one interface below the radiation emitter, located at $z^{\prime}$,
we have that

$\mathbf{F}_{z<z^{\prime}}(z,z^{\prime})=\mathbf{F}_{z>z^{\prime}}(z,z^{\prime})=e^{ik_{z}|z-z^{\prime}|}\mathbf{I}+e^{ik_{z}(z+z^{\prime})}\mathbf{R}.$
(17)

Note that if we introduce $\mathbf{F}(z,z^{\prime})$ into equation (3) we
recover the electric free-space dyadic Green’s function $\mathbf{G}_{0}$, plus
another extra term corresponding to the reflected dyadic Green’s function
$\mathbf{G}_{R}$,

$\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)=\mathbf{G}_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega)+\mathbf{G}_{R}(\mathbf{r},\mathbf{r}^{\prime},\omega).$
(18)

##  III Adding a Mu-Metal Layer

As the intrinsic parameters of the TI are difficult to change without altering
its properties, it is far more convenient to find a way to incorporate another
material with a very high value for $\mu$ that enhances the mixed reflection
coefficients similarly, as $\Delta\sim\mu_{1}\mu_{2}$. This addition increases
the relative scale of these coefficients with respect to the usual
coefficients. Thus, we propose to examine a three-layered system with a Mu-
metal with $\mu\sim 10^{5}\mu_{0}$, as the one seen in figure 2. It remains to
be seen where to set this layer: over or under the topological insulator.

Figure 2: Sketch of a TE polarized ray impacting the three-layered Air-TI-Mu-
Metal system, where the reflected wave has electric and magnetic field
components in all directions due to the Electromagnetic coupling.

To determine the ordering of the layers, and how the reflection matrix will
behave in each case, we study the reflection model for a 3-layered system Chew
(1995). By taking into account not only the first reflected wave but also the
transmitted portion which is reflected on the second interface and transmitted
once again through the first one, we can define a new reflection matrix, such
that

$\mathbf{\tilde{R}}_{12}=\mathbf{R}_{12}+e^{2ik_{z2}(d_{12}-d_{23})}\mathbf{T}_{21}\cdot\mathbf{R}_{23}\cdot\mathbf{M}_{2123}^{-1}\cdot\mathbf{T}_{12},$
(19)

where

$\mathbf{M}_{2123}=\mathbf{I}-e^{2ik_{z2}(d_{12}-d_{23})}\mathbf{R}_{21}\cdot\mathbf{R}_{23}.$
(20)

If the emission source remains on the first layer of our configuration, it is
only necessary to replace the reflection matrix in $\mathbf{F}_{z\lessgtr
z^{\prime}}(z,z^{\prime})$ with (19).

###  1 3-Layers Fresnel Coefficients

For the following numerical results, we used the optical properties of room-
temperature TI TlBiSe2Sato  _et al._ (2010); Mitsas  _et al._ ; Castro-
Enriquez  _et al._ (2020); Phutela  _et al._ (2022), as its relatively high
impedance magnifies the topological effects previously discussed. These
optical properties are $\epsilon\approx 4,\mu=1$. We also consider another TI
with $\epsilon=1,\mu=2$, which can be achieved by increasing the magnetic
permeability through dopingChen  _et al._ (2019); Kim  _et al._ (2019); Choi
_et al._ (2011). Magnetically doping TIs is a prominent technique for
enhancing the anomalous magnetic effects displayed, despite the risk of losing
the topological properties due to anisotropy in the BulkTeng  _et al._ (2019).
We use the latter optical properties in contrast to TlBiSe2 to show relevant
dependencies that arise with the TI’s impedance.

In figure 3 we can see how, as we increase the mu-metal $\mu$ value, the mixed
reflective coefficients tend asymptotically to a maximum value, so after a
certain $\mu$ the gain from increasing this value will diminish. As for the
thickness of the topological insulator, due to the oscillatory nature of the
exponential terms on equations (19) and (20), for a given angle of incidence
over the surface, we can select an appropriate thickness on which the mixed
reflection coefficients are at their maximum. In figure 4(a), we can see this
behavior mirrored, and we can use these peaks to determine how the TI’s
parameters influence these increases.

Figure 3: Mixed reflective coefficients as a function of the $\mu$ of the mu-
metal, for different incident angles, with
$\mu_{1}=\epsilon_{1}=\epsilon_{\Theta}=1$, $\mu_{\Theta}=2$,
$\lambda=600\;\text{nm}$ and $d_{12}-d_{23}=100\;\text{nm}$.

Here it is also clear that the Air/TI/Mu-Metal configuration’s peaks deviate
above from the baseline, while the Air/Mu-Metal/TI configuration peaks are at
most at the base value. This last phenomenon can be attributed to Mu-Metal
shieldingYamazaki  _et al._ (2005); ter Brake  _et al._ (1991); Gao  _et al._
(2023); Dubbers (1986), as starting to increase the thickness of the Mu-Metal
in this configuration results in the total reflection of the light before it
reaches the TI layer, causing any mixed reflection coefficients to nullify.
This argument is sustained by the reflection and transmission matrices in
table 1, for the Mu-Metal cases.

As seen in figure 4(b), for TI with lower impedance, we are able to enhance by
approximately half an order of magnitude the mixed reflection coefficients
introducing the Air/TI/Mu-Metal configuration, taking into account the optimal
thickness for these parameters obtained numerically. The Air/Mu-Metal/TI
configuration is found to peak at the base configuration of Air/TI. In figure
4(b), we also see that the increase is more significant as the TI’s impedance
is raised. Meanwhile, the maximum of the Air/Mu-Metal/TI does not improve with
higher impedance. We predict mixed coefficients around $10^{2}$ times bigger
than those without the Mu-Metal sub-layer.

(a) |
---|---
(b) |
|
Figure 4: (a) Mixed reflective coefficients as a function of the TI layer’s
thickness, in units of the optical wavelength in the medium ($\lambda/n_{2}$),
to visualize the same amount of peaks in both configurations. Here, the angle
of incidence $\phi_{i}=0^{\circ}$, $\mu_{\Theta}=2$ and $\epsilon_{\Theta}=1$.
(b) Mixed reflective coefficients as a function of the incident angle for the
three different possible configurations for two different TIs, with
$\mu_{1}=\epsilon_{1}=\epsilon_{\mu-\text{metal}}=1$,
$\mu_{\mu-\text{metal}}=10^{5}$, $\lambda=600\;\text{nm}$ and middle layer
thickness is strategically chosen to maximize them. The TI’s parameters of
TlBiSe2 are $\mu_{\Theta}=1$ and $\epsilon_{\Theta}=4$, and for the doped TI
are $\mu_{\Theta}=2$ and $\epsilon_{\Theta}=1$.

We did not consider absorption in our models. In this case, absorption would
add an exponential decay to the reflection coefficients, which depends on the
middle layer’s thickness. Thus, selecting the thickness corresponding to the
first peak in the coefficients is advised.

We make the approximation
$\mu_{3}\gg\mu_{2},\mu_{1},\epsilon_{1},\epsilon_{2},\epsilon_{3}$ in order to
show the algebraic dependencies on the parameters. From (20), we obtain

$\displaystyle\mathbf{M}_{2123}=$
$\displaystyle\left(\begin{matrix}1-\mathcal{F}(R_{\textit{EE}}^{21}R_{\textit{EE}}^{23}+R_{\textit{EM}}^{21}R_{\textit{ME}}^{23})&-\mathcal{F}(R_{\textit{EE}}^{21}R_{\textit{EM}}^{23}+R_{\textit{EM}}^{21}R_{\textit{MM}}^{23})\\\
-\mathcal{F}(R_{\textit{ME}}^{21}R_{\textit{EE}}^{23}+R_{\textit{MM}}^{21}R_{\textit{ME}}^{23})&1-\mathcal{F}(R_{\textit{ME}}^{21}R_{\textit{EM}}^{23}+R_{\textit{MM}}^{21}R_{\textit{MM}}^{23})\end{matrix}\right),$
(21)

where $\mathcal{F}=e^{2ik_{z2}(d_{12}-d_{23})}$ and, for example,
$R_{\textit{EE}}^{21}$ is the $R_{\textit{TE,TE}}$ entry of the
$\mathbf{R}_{21}$ reflection matrix. Now, using our approximation, we can
reduce the $\mathbf{R}_{23}$ such that,

$\mathbf{R}_{23}\approx\left(\begin{matrix}1&-2Z_{2}\alpha(\Theta_{3}-\Theta_{2})\\\
-2Z_{2}\alpha(\Theta_{3}-\Theta_{2})&-1\end{matrix}\right),$ (22)

$Z_{2}=\sqrt{\mu_{2}/\epsilon_{2}}$ being the TI’s impedance. Multiplying (1)
with this approximation, we obtain the following matrix,

$\displaystyle\mathbf{R}_{23}\cdot\mathbf{M}_{2123}^{-1}=$
$\displaystyle\frac{1}{\det(\mathbf{M}_{2123})}\left(\begin{matrix}1&-\beta\\\
-\beta&-1\end{matrix}\right)+\frac{(1+\beta^{2})\mathcal{F}}{\det(\mathbf{M}_{2123})}\left(\begin{matrix}R_{\textit{MM}}^{21}&R_{\textit{EM}}^{21}\\\
R_{\textit{ME}}^{21}&R_{\textit{EE}}^{21}\end{matrix}\right),$ (23)

where $\beta=2Z_{2}\alpha(\Theta_{3}-\Theta_{2})/\pi$. It is worth remarking
that new terms appear on the off-diagonal components of our matrix that are
dependant on $\beta/\det(\mathbf{M}_{2123})$. Given that $\beta\sim 10^{-2}$
in order, and $\det(\mathbf{M}_{2123})\sim 1$, these new terms should dominate
over the mixed reflection coefficients that are of order $10^{-3}$. Given
this, we roughly estimate the order of the new reflection coefficients with
equations (19) and (12). We see that the maximum of the new reflection
coefficients is of order $10^{-4}$, which is shown in figure 4.

Therefore, we can confirm that the increase in the mixed reflection
coefficients, due to the addition of a Mu-Metal layer, directly correlates to
the off-diagonal terms in the first term of equation (1). The increase scales
linearly with the TI’s impedance, as long as $\mu_{3}\gg\mu_{2}$ is satisfied.

###  2 3-Layers Green’s Function/Electric Field

In this section, we will display the results obtained for the Green’s function
over the TI. From equation (13), we can understand each column of
$\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)$ as a scaled representation
of the directional real components of the electric field generated by an
oscillating dipole on $\mathbf{r}^{\prime}$, oriented along the
$\hat{x},\hat{y}$ and $\hat{z}$ directions, as for that case

$\mathbf{J}(\mathbf{r}^{\prime})=-i\omega
d(\omega)\delta(\mathbf{r}^{\prime})\hat{e}_{\mu},$ (24)

with $\hat{e}_{\mu}$ the unitary dipole orientation and $d(\omega)$ the dipole
strength. We obtain that the electric field is proportional by a factor to
$\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)\cdot\hat{e}_{\mu}$.

For free-space and Air/Magnetodielectric configurations, for an $\hat{x}$
oriented dipole, the electric field $\hat{y}$ component in the $x-z$ plane is
null, and the $\hat{x}$ and $\hat{z}$ components are not. This can be seen in
figures 8 and 9 in appendix A and is crucial to characterize the difference
between magnetodielectric materials and topological insulators, as this will
not continue to be true for topologically complex materials.

Now, analyzing the 2-layer Air/TI system configuration, as we can see in
figure 5 (a), the mixed reflected coefficients caused the y-combined matrix
components to appear, which can be physically seen by measuring, for example,
the $y$-oriented component of the electric field of an $x$ or $z$-oriented
dipole near the surface. Given the relatively small value of the mixed
reflection coefficients, the mixed components of the matrix have been scaled
by a $\times 1000$ factor to be able to appreciate them in the same magnitude
as the non-mixed entries.

Figure 5: Electric field components $\mathbf{E}_{i}\times\lambda$, graphed in
the $y=0$ plane, as a function of $z/\lambda$ and $x/\lambda$ for (a) an
Air/TI configuration, and (b) an Air/TI/Mu-Metal configuration. $E_{y}$ is
solely due to the topological response of the configuration and had to be
enhanced to be observed with the non-topological $E_{x}$. Note that the
necessary enhancement factor is reduced by an order of magnitude when adding a
Mu-Metal to the system. The corresponding parameters, other components for
each figure, and the components for the free-space and Air/Magnetodielectric
configurations can be found in appendix A.

As previously discussed, the mixed Fresnel coefficients can be enhanced by
implementing an Air/TI/Mu-metal 3-layered configuration. In figure 5 (b), we
can see how the mixed entries of the electric field are now scaled by a
$\times 100$ factor, which means we have reduced by an order of magnitude the
threshold necessary to be able to measure this effect via the electric field
in an experimental environment. We can assure that this enhancement of the
effect can be greater for different values of $\mu$ for the topological
insulator layer, for which we can determine the optimal thickness of the TI,
as discussed in the previous section.

##  IV New Component on Poynting Vector

During previous sections, we have studied how the electric and magnetic
fields, upon reflection or transmission on a TI, present components not only
in their original polarization but also in the complementary polarization. In
this section, we will study if this is also the case for the Poynting vector.
If there is a substantial deviation, this might be the basis for a new way to
characterize and measure this effect.

So, we are interested in the time-averaged Poynting vector, given by

$\mathbf{S}=\frac{1}{2}\text{Re}\\{\mathbf{E}\times\mathbf{H}^{*}\\}.$ (25)

Figure 6: Poynting vector (first row) and the difference between the Poynting
vector (second row) generated by an $\hat{x}$ oriented dipole located at
$z^{\prime}=0.05\times\lambda$ over a TI’s surface $\mathbf{S_{\theta}}$, and
over the equivalent magnetodielectric surface $\mathbf{S_{0}}$, normalized by
the average of the non-topological Poynting vector over the plane. Simulating
an 0.9 numerical aperture (NA) at $z=0.1\times\lambda$, $0.35\times\lambda$,
and $0.75\times\lambda$ over the TI’s surface. All plots consider the
parameters $\mu_{1}=\epsilon_{1}=\mu_{\Theta}=1$, $\epsilon_{\Theta}=4$,
$\lambda=600\;\text{nm}$, and $\Theta_{2}=\pi$ for the TI.

Upon taking the ansatz for a single TE polarized planar wave and calculating
the corresponding Poyinting vector for the top layer of an Air-TI
configuration, we obtain the following solution,

$\mathbf{S}=\frac{E_{0}^{2}}{2\mu_{0}\mu_{1}\omega}\left(\begin{matrix}k_{p}\left[1+R_{\textit{TE,TE}}^{2}+2\cos{(2k_{z})}R_{\textit{TE,TE}}-R_{\textit{TM,TE}}^{2}\right]\\\
\frac{k_{z1}k_{p}}{k_{1}}R_{\textit{TE,TE}}\;R_{\textit{TM,TE}}\\\
k_{z1}\left[1+R_{\textit{TE,TE}}^{2}+R_{\textit{TM,TE}}^{2}\right]\end{matrix}\right).$
(26)

The important feature of equation (26) is its $\hat{y}$ component, which is
nullified for regular magnetodielectric media. This new component suggests a
relatively small transversal deviation of the Poynting vector.

For a dipole over a TI’s surface, we can expand the analysis for the Green’s
function to also obtain the magnetic field, given byNovotny and Hecht (2012)

$\mathbf{H}(\mathbf{r},\omega)=\int
d^{3}r^{\prime}\left[\nabla\times\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)\right]\cdot\mathbf{J}(\mathbf{r}^{\prime}).$
(27)

Expanding on the result for the electric Green’s functionCrosse  _et al._
(2015), computing its rotor, and using the corresponding Hankel
transformGaskill (1978), we are able to add the reflected part of the magnetic
Green’s function for a TI half-plane space. Here are the results for an
$\hat{x}$ oriented dipole:

$\displaystyle M^{\,xx}_{R}$
$\displaystyle(\mathbf{r},\mathbf{r}^{\prime},\omega)=\frac{i\mu(\omega)}{8\pi}\int_{0}^{k_{1}}dk_{p}k_{p}\;e^{ik_{z}(|z|+|z^{\prime}|)}\left[J_{0}(k_{p}R_{p})\left(\frac{k_{1}}{k_{z}}+\frac{k_{z}}{k_{1}}\right)-J_{2}(k_{p}R_{p})\frac{k_{p}^{2}}{k_{1}\,k_{z}}\right]R_{\text{TE,TM}},$
(28) $\displaystyle M^{\,yx}_{R}$
$\displaystyle(\mathbf{r},\mathbf{r}^{\prime},\omega)=\frac{i\mu(\omega)}{8\pi}\int_{0}^{k_{1}}dk_{p}k_{p}\;e^{ik_{z}(|z|+|z^{\prime}|)}\Big{[}(J_{0}(k_{p}R_{p})+J_{2}(k_{p}R_{p}))(R_{\text{TE,TE}}+R_{\text{TM,TM}})+J_{2}(k_{p}R_{p})\frac{2k_{p}k_{z}^{2}}{k_{1}^{2}}R_{\text{TM,TM}}\Big{]},$
(29) $\displaystyle M^{\,zx}_{R}$
$\displaystyle(\mathbf{r},\mathbf{r}^{\prime},\omega)=\frac{\mu(\omega)}{4\pi}\int_{0}^{k_{1}}dk_{p}e^{ik_{z}(|z|+|z^{\prime}|)}J_{1}(k_{p}R_{p})\;\frac{k_{p}^{2}}{k_{1}^{2}}R_{\text{TE,TM}},$
(30)

where $k_{z}=\sqrt{k_{1}^{2}-k_{p}^{2}}$, and $J_{0}$, $J_{1}$ and $J_{2}$ are
first kind Bessel functions. This can also be extended to the rest of the
dipole orientations. Substituting these results and the electrical field’s
Green’s function into equation (26) and computing the difference between this
vector field for TI and its non-topological counterpart
$\mathbf{S}_{\Theta}-\mathbf{S}_{0}$, we obtain the graph on figure 6. Here we
see a clear new transverse component on the $\hat{y}$ direction, due to the
topological effect.

Figure 7: Topological Poynting vector ($\mathbf{S}_{\Theta}$) and the
topological Poynting vector’s deviation ($\mathbf{S}_{\Theta}-\mathbf{S}_{0}$)
normalized by $\mathbf{S}_{\Theta}$ on the far-field. Note that both have
components that decay as $1/r^{2}$. The oscillatory behavior on the deviation
is due to the Bessel functions dependency of $\mathbf{S}_{\Theta}$.

We also notice that the relative difference tends to decrease the further we
go from the emission source. Looking at equations (28), (29) and (30),
specifically their dependence on the Bessel functions $J_{n}(k_{p}\cdot
R_{p})$ and the exponential, we expect that the results for the Poynting
vector should scale by $\sim 1/r^{2}$, as evidenced in figure 7. This
indicates that the relative difference is indeed a far-field effect, as the
integral over the solid angle should be constant, as in the case of the non-
topological Poynting vector.

Moreover, for a room temperature TI with a lower impedance, we should expect
this effect to be lower, although we are able to find the optimal
configuration of the dipole-TI system to find the maximum possible percentual
deviation. Figure 6 shows a maximum transversal deviation of $0.18\%$ over the
expected Poynting vector of the equivalent magnetodielectric.

##  V Conclusions

We considered the inclusion of a Mu-Metal layer in the system Air/Topological
Insulator to enhance the optical effects of the mixed reflection Fresnel
coefficients. We found that it is most beneficial to situate it under the
Topological Insulator layer and that the increase is more substantial for TIs
with higher impedance. We also discussed how to find the optimal thickness for
the TI middle layer, given the angle of incidence and the TI’s properties. In
our simulations, we predict an increase of $\sim\times 100$ for a TI with
impedance $Z=\sqrt{2}$.

Finally, we extended the electric Green’s function analysis to obtain the
magnetic Green’s function and the Poynting vector field for a dipole over a
TI. We found that the Poynting vector presents a new transverse component due
to the new mixed Fresnel coefficients, which scale like $1/r^{2}$, where $r$
is the distance to the source, indicating this effect will be present on the
far field. Moreover, for a room-temperature TI, we predict a maximum
percentual deviation of $0.18\%$, which indicates that this effect is
measurable given the right experimental setup.

Since the result showed that the new magnetic components that appear in TIs
are enhanced, we can more confidently expect that magnetic probing, using
single emitters, might prove a reliable way to characterize TIs. Moreover, the
enhancement shown by adding Mu-Metal materials may be aided with additional
meta-materials in a variety of geometrical forms being introduced into the
layered configuration to enhance these effects further.

##  VI ACKNOWLEDGEMENTS

The authors acknowledge the support from the Asian Office of Aerospace
Research and Development (AOARD) FA2386-21-1-4125, Fondecyt Regular No
1221512, and Anillo ACT192023 for this work.

## References

  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Hasan and Moore (2011) M. Z. Hasan and J. E. Moore, Annual Review of Condensed Matter Physics 2, 55 (2011).
  * Ando (2013) Y. Ando, Journal of the Physical Society of Japan 82, 102001 (2013).
  * Qi  _et al._ (2009) X.-L. Qi, R. Li, J. Zang, and S.-C. Zhang, Science 323, 1184 (2009).
  * Essin  _et al._ (2009) A. M. Essin, J. E. Moore, and D. Vanderbilt, Physical Review Letters 102 (2009), 10.1103/physrevlett.102.146805.
  * Martín-Ruiz  _et al._ (2019) A. Martín-Ruiz, O. Rodríguez-Tzompantzi, J. R. Maze, and L. F. Urrutia, Phys. Rev. A 100, 042124 (2019).
  * Kitagawa  _et al._ (2010) Y. Kitagawa, Y. Hiraoka, T. Honda, T. Ishikura, H. Nakamura, and T. Kimura, Nature Materials 9, 797 (2010).
  * Velev  _et al._ (2011) J. P. Velev, S. S. Jaswal, and E. Y. Tsymbal, Philosophical Transactions: Mathematical, Physical and Engineering Sciences 369, 3069 (2011).
  * Moore (2010) J. E. Moore, Nature 464, 194 (2010).
  * Legg  _et al._ (2022) H. F. Legg, M. Rößler, F. Münning, D. Fan, O. Breunig, A. Bliesener, G. Lippertz, A. Uday, A. A. Taskin, D. Loss, J. Klinovaja, and Y. Ando, Nature Nanotechnology 17, 696 (2022).
  * Fischer  _et al._ (2022) R. Fischer, J. Picó-Cortés, W. Himmler, G. Platero, M. Grifoni, D. A. Kozlov, N. N. Mikhailov, S. A. Dvoretsky, C. Strunk, and D. Weiss, Phys. Rev. Res. 4, 013087 (2022).
  * Münning  _et al._ (2021) F. Münning, O. Breunig, H. F. Legg, S. Roitsch, D. Fan, M. Rößler, A. Rosch, and Y. Ando, Nature Communications 12 (2021), 10.1038/s41467-021-21230-3.
  * Wu  _et al._ (2021) H. Wu, A. Chen, P. Zhang, H. He, J. Nance, C. Guo, J. Sasaki, T. Shirokura, P. N. Hai, B. Fang, S. A. Razavi, K. Wong, Y. Wen, Y. Ma, G. Yu, G. P. Carman, X. Han, X. Zhang, and K. L. Wang, Nature Communications 12 (2021), 10.1038/s41467-021-26478-3.
  * Castro-Enríquez  _et al._ (2022) L. A. Castro-Enríquez, A. Martín-Ruiz, and M. Cambiaso, Scientific Reports 12 (2022), 10.1038/s41598-022-24939-3.
  * Lin  _et al._ (2015) Y.-H. Lin, S.-F. Lin, Y.-C. Chi, C.-L. Wu, C.-H. Cheng, W.-H. Tseng, J.-H. He, C.-I. Wu, C.-K. Lee, and G.-R. Lin, ACS Photonics 2, 481 (2015).
  * He  _et al._ (2019) M. He, H. Sun, and Q. L. He, Frontiers of Physics 14 (2019), 10.1007/s11467-019-0893-4.
  * Pandey  _et al._ (2021) A. Pandey, R. Yadav, M. Kaur, P. Singh, A. Gupta, and S. Husale, Scientific Reports 11 (2021), 10.1038/s41598-020-80738-8.
  * Breunig and Ando (2021) O. Breunig and Y. Ando, Nature Reviews Physics 4, 184 (2021).
  * Dziom  _et al._ (2017) V. Dziom, A. Shuvaev, A. Pimenov, G. V. Astakhov, C. Ames, K. Bendias, J. Böttcher, G. Tkachov, E. M. Hankiewicz, C. Brüne, H. Buhmann, and L. W. Molenkamp, Nature Communications 8 (2017), 10.1038/ncomms15197.
  * C and F (2017) G. E. J. C and R. D. F, Journal of Physics: Conference Series 850, 012024 (2017).
  * Crosse (2017) J. A. Crosse, Scientific Reports 7 (2017), 10.1038/srep43115.
  * Shuvaev  _et al._ (2022) A. Shuvaev, L. Pan, L. Tai, P. Zhang, K. L. Wang, and A. Pimenov, Applied Physics Letters 121, 193101 (2022).
  * Tse and MacDonald (2010) W.-K. Tse and A. H. MacDonald, Physical Review Letters 105 (2010), 10.1103/physrevlett.105.057401.
  * Tse and MacDonald (2011) W.-K. Tse and A. H. MacDonald, Physical Review B 84 (2011), 10.1103/physrevb.84.205327.
  * Crosse (2016) J. A. Crosse, Physical Review A 94 (2016), 10.1103/physreva.94.033816.
  * Ardakani and Zare (2021) A. G. Ardakani and Z. Zare, J. Opt. Soc. Am. B 38, 2562 (2021).
  * Martín-Ruiz  _et al._ (2016) A. Martín-Ruiz, M. Cambiaso, and L. Urrutia, Physical Review D 94 (2016), 10.1103/physrevd.94.085019.
  * Ge  _et al._ (2014) L. Ge, T. Zhan, D. Han, X. Liu, and J. Zi, Optics Express 22, 30833 (2014).
  * Gangaraj  _et al._ (2020) S. A. H. Gangaraj, C. Valagiannopoulos, and F. Monticone, Physical Review Research 2 (2020), 10.1103/physrevresearch.2.023180.
  * Mogi  _et al._ (2022) M. Mogi, Y. Okamura, M. Kawamura, R. Yoshimi, K. Yasuda, A. Tsukazaki, K. S. Takahashi, T. Morimoto, N. Nagaosa, M. Kawasaki, Y. Takahashi, and Y. Tokura, Nature Physics 18, 390 (2022).
  * Martín-Ruiz  _et al._ (2021) A. Martín-Ruiz, M. Cambiaso, and L. F. Urrutia, in  _Topics in Applied Physics_ (Springer International Publishing, 2021) pp. 459–492.
  * Crosse  _et al._ (2015) J. A. Crosse, S. Fuchs, and S. Y. Buhmann, Phys. Rev. A 92, 063831 (2015).
  * Chew (1995) W. C. Chew,  _Waves and fields in inhomogeneous media_ (IEEE Press, New York, 1995).
  * Novotny and Hecht (2012) L. Novotny and B. Hecht,  _Principles of Nano-Optics_, 2nd ed. (Cambridge University Press, 2012).
  * Liu and Hesjedal (2021) J. Liu and T. Hesjedal, Advanced Materials , 2102427 (2021).
  * Mogi  _et al._ (2017) M. Mogi, M. Kawamura, R. Yoshimi, A. Tsukazaki, Y. Kozuka, N. Shirakawa, K. S. Takahashi, M. Kawasaki, and Y. Tokura, Nature Materials 16, 516 (2017).
  * An  _et al._ (2020) H. An, R. Zeng, M. Zhang, H. Li, M. Hu, Q. Li, and X. Zeng, Optics Communications 477, 126335 (2020).
  * Zuo  _et al._ (2013) Z.-W. Zuo, D.-B. Ling, L. Sheng, and D. Xing, Physics Letters A 377, 2909 (2013).
  * Sato  _et al._ (2010) T. Sato, K. Segawa, H. Guo, K. Sugawara, S. Souma, T. Takahashi, and Y. Ando, Phys. Rev. Lett. 105, 136802 (2010).
  * (40) C. Mitsas, D. I. Siapkas, E. K. Polychroniadis, O. Valassiades, and K. M. Paraskevopoulos, Physica status solidi 136, 483.
  * Castro-Enriquez  _et al._ (2020) L. A. Castro-Enriquez, L. F. Quezada, and A. Martín-Ruiz, Phys. Rev. A 102, 013720 (2020).
  * Phutela  _et al._ (2022) A. Phutela, P. Bhumla, M. Jain, and S. Bhattacharya, Scientific Reports 12 (2022), 10.1038/s41598-022-26445-y.
  * Chen  _et al._ (2019) B. Chen, F. Fei, D. Zhang, B. Zhang, W. Liu, S. Zhang, P. Wang, B. Wei, Y. Zhang, Z. Zuo, J. Guo, Q. Liu, Z. Wang, X. Wu, J. Zong, X. Xie, W. Chen, Z. Sun, S. Wang, Y. Zhang, M. Zhang, X. Wang, F. Song, H. Zhang, D. Shen, and B. Wang, Nature Communications 10 (2019), 10.1038/s41467-019-12485-y.
  * Kim  _et al._ (2019) J. Kim, E.-H. Shin, M. K. Sharma, K. Ihm, O. Dugerjav, C. Hwang, H. Lee, K.-T. Ko, J.-H. Park, M. Kim, H. Kim, and M.-H. Jung, Scientific Reports 9 (2019), 10.1038/s41598-018-37663-8.
  * Choi  _et al._ (2011) Y. H. Choi, N. H. Jo, K. J. Lee, J. B. Yoon, C. Y. You, and M. H. Jung, Journal of Applied Physics 109, 07E312 (2011).
  * Teng  _et al._ (2019) J. Teng, N. Liu, and Y. Li, Journal of Semiconductors 40, 081507 (2019).
  * Yamazaki  _et al._ (2005) K. Yamazaki, K. Kato, K. Muramatsu, A. Haga, K. Kobayashi, K. Kamata, K. Fujiwara, and T. Yamaguchi, IEEE transactions on magnetics 41, 4087 (2005).
  * ter Brake  _et al._ (1991) H. J. ter Brake, H. Wieringa, and H. Rogalla, Measurement Science and Technology 2, 596 (1991).
  * Gao  _et al._ (2023) Y. Gao, D. Ma, K. Wang, X. Xu, S. Li, Y. Dou, and J. Li, Sensors and Actuators A: Physical 352, 114207 (2023).
  * Dubbers (1986) D. Dubbers, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 243, 511 (1986).
  * Gaskill (1978) J. D. Gaskill,  _Linear systems, Fourier transforms, and optics_ , Wiley Series in Pure and Applied Optics (John Wiley & Sons, Nashville, TN, 1978).

##  Appendix A Green’s Function components plots for different configurations

Figure 8: Respective matrix components of
$\mathbf{G}_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega)\times\lambda$, graphed
in the $y=0$ plane, as a function of $z/\lambda$ and $x/\lambda$ for free-
space, with parameters $\mu_{1}=\epsilon_{1}=1$ and $\lambda=600\;\text{nm}$.
Here, the emitting dipole is arbitrarily located at
$\mathbf{r}^{\prime}=1.5\times\lambda$. Here $\mathbf{G}_{ij}$ can be
understood as the $E_{i}$ electric field component generated by a dipole
oriented in $j$. Note that the off-diagonal components are null, regardless of
the enhancement factor. Figure 9: Respective matrix components of
$\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)\times\lambda$, graphed in
the $y=0$ plane, as a function of $z/\lambda$ and $x/\lambda$ for the
Air/Magnetodielectric configuration, with parameters
$\mu_{1}=\epsilon_{1}=\epsilon_{2}=1$, $\mu_{2}=2$, $\Theta_{2}=0$ and
$\lambda=600\;\text{nm}$, and the interface located at $z=0$. Here, the
emitting dipole is arbitrarily located at
$\mathbf{r}^{\prime}=1.5\times\lambda$. Here $\mathbf{G}_{ij}$ can be
understood as the $E_{i}$ electric field component generated by a dipole
oriented in $j$. Note that the off-diagonal components remain null, regardless
of the enhancement factor. Figure 10: Respective matrix components of
$\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)\times\lambda$, graphed in
the $y=0$ plane, as a function of $z/\lambda$ and $x/\lambda$ for the Air/TI
configuration, with parameters $\mu_{1}=\epsilon_{1}=\epsilon_{2}=1$,
$\mu_{2}=2$, $\Theta_{2}=\pi$ and $\lambda=600\;\text{nm}$, and the interface
located at $z=0$. Here, the emitting dipole is arbitrarily located at
$\mathbf{r}^{\prime}=1.5\times\lambda$. Here $\mathbf{G}_{ij}$ can be
understood as the $E_{i}$ electric field component generated by a dipole
oriented in $j$. Note that here we had to enhance the off-diagonal components
by a factor of $1000$ and $10000$, respectively, to see them with the normal
emissions. These components appear due to topological contributions. Figure
11: Respective matrix components of
$\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime},\omega)\times\lambda$, graphed in
the $y=0$ plane, as a function of $z/\lambda$ and $x/\lambda$ for the
Air/TI/Mu-metal configuration, with parameters
$\mu_{1}=\epsilon_{1}=\epsilon_{2}=\epsilon_{3}=1$, $\mu_{2}=2$,
$\mu_{3}=10^{5}$, $\Theta_{2}=\pi$, $\Theta_{3}=0$, $d=100\;\text{nm}$
(thickness of the TI-layer) and $\lambda=600\;\text{nm}$. Here, the emitting
dipole is arbitrarily located at $\mathbf{r}^{\prime}=1.5\times\lambda$. Here
$\mathbf{G}_{ij}$ can be understood as the $E_{i}$ electric field component
generated by a dipole oriented in $j$. Note that here we had to enhance the
off-diagonal components by a factor of $100$ to see them with the normal
emissions. This implies an enhancement of an order of magnitude for $G_{xy}$
and $G_{yx}$, and of two orders of magnitude for $G_{yz}$ and $G_{zy}$,
compared with the Air/TI configuration.


# Trigonal warping effects on optical properties of anomalous Hall materials

Jyesta M. Adhidewata<EMAIL_ADDRESS>Research Center for Quantum Physics,
National Research and Innovation Agency (BRIN), South Tangerang 15314,
Indonesia  Theoretical High Energy Physics Research Division, Faculty of
Mathematics and Natural Sciences, Institut Teknologi Bandung, Bandung 40132,
Indonesia  Ravanny W. M. Komalig  Research Center for Quantum Physics,
National Research and Innovation Agency (BRIN), South Tangerang 15314,
Indonesia  Advanced Functional Materials Research Group, Faculty of Industrial
Technology, Institut Teknologi Bandung, Bandung 40132, Indonesia  M Shoufie
Ukhtary  Research Center for Quantum Physics, National Research and Innovation
Agency (BRIN), South Tangerang 15314, Indonesia  Ahmad R. T. Nugraha  Research
Center for Quantum Physics, National Research and Innovation Agency (BRIN),
South Tangerang 15314, Indonesia  Department of Engineering Physics, School of
Electrical Engineering, Telkom University, Bandung 40257, Indonesia  Research
Collaboration Center for Quantum Technology 2.0, Bandung 40132, Indonesia
Bobby E. Gunara<EMAIL_ADDRESS>Theoretical High Energy Physics Research
Division, Faculty of Mathematics and Natural Sciences, Institut Teknologi
Bandung, Bandung 40132, Indonesia  Eddwi H. Hasdeo
<EMAIL_ADDRESS>Research Center for Quantum Physics, National
Research and Innovation Agency (BRIN), South Tangerang 15314, Indonesia
Department of Physics and Materials Science, University of Luxembourg, L-1511
Luxembourg, Luxembourg

###### Abstract

The topological nature of topological insulators are related to the symmetries
present in the material, for example, quantum spin Hall effect can be observed
in topological insulators with time reversal symmetry, while broken time
reversal symmetry may give rise to the presence of anomalous quantum Hall
effect (AHE). Here we consider the effects of broken rotational symmetry on
the Dirac cone of an AHE material by adding trigonal warping terms to the
Dirac Hamiltonian. We calculate the linear optical conductivity semi-
analytically to show how by breaking the rotational symmetry we can obtain a
topologically distinct phase. The addition of trigonal warping terms causes
the emergence of additional Dirac cones, which when combined has a total Chern
number of $\mp 1$ instead of $\pm 1/2$. This results in drastic changes in the
anomalous Hall and longitudinal conductivity. The trigonal warping terms also
activates the higher order Hall responses which does not exist in a
$\mathcal{R}$ symmetric conventional Dirac material. We found the presence of
a non-zero second order Hall current even in the absence of Berry curvature
dipole. This shift current is also unaffected by the chirality of the Dirac
cone, which should lead to a non-zero Hall current in time reversal symmetric
systems.

##  I Introduction

The quantum anomalous Hall effect (QAHE) is the presence of a quantized Hall
conductance in a material without an external magnetic field due to time-
reversal (TR) symmetry breaking [1, 2, 3]. Similar to the quantum Hall effect,
QAHE is a manifestation of the topological properties of a material [4, 5, 3].
In QAHE, the Hall conductivity is determined by the Chern number, which is a
topological invariant given by the integral of the Berry curvature over the
Brillouin zone [6, 7, 8]. As a topological invariant, the Chern number is
insensitive to small changes applied to the bulk form of the material.
However, several ways have been shown to control the winding number of the
wave functions allowing the emergence of higher Chern insulators [9].
Therefore, it is intriguing for us to present a particular case where the
Chern number can change due to the breaking of rotational symmetry
$\mathcal{R}$. One possible model of QAHE material is a gapped Dirac material
with broken TR symmetry in the two-band model [10], which can be used to
describe the Chern insulators or magnetic topological insulators [11].

In this work, we examine the optical conductivity of a broken-TR-symmetry
Dirac material with trigonal warping terms. The trigonal warping, as occurs in
trigonal or honeycomb lattices, breaks the rotational symmetry of the band
structure. We find that the breaking of rotational symmetry can add new three
Dirac cones carrying half integer Chern number with opposite sign from the one
at K point giving the total values of Chern number equals $-1$, in contrast to
the regular Dirac dispersion which carries a $+1/2$ Chern number. This opens a
new possibility to find higher Chern insulators by tracking the strength of
trigonal warping and the location of the Dirac cones. Another advantage of
breaking this symmetry is the activation of higher-harmonic red and
longitudinal currents arising as a nonlinear response to external fields since
those currents are normally suppressed in the $\mathcal{R}$ symmetric system.
It is well-known that a nonlinear anomalous Hall current can be observed in
materials without inversion symmetry [12, 13], such as transition metal
dichalcogenides [14] and strained graphene [15] (even in the presence of TR
symmetry). This nonlinear current is usually attributed to the presence of a
Berry curvature dipole, although similar nonlinear transverse current might
arise from other sources, such as the semiclassical Jerk term [16]. Here we
show that, even with zero Berry curvature dipole, we still have a nonzero
transverse current contribution from the nonlinear shift current in our Dirac
material. The shift current arising from a single Dirac point is unaffected by
the chirality of that point, so that it warrants the presence of a transverse
current even in the TR-symmetric systems.

We organize this paper as follows. In Sec. II, we describe our model
Hamiltonian and the general band structure of our Hamiltonian. In Sec. III we
describe the method used for calculation of the linear conductivity and its
results. We discuss the effects of the trigonal warping such as a non-trivial
Chern number and a van Hove singularity. In Sec. IV we described our
calculation of nonlinear conductivity, and we discuss the results of our
calculation. Finally, we give conclusion of this work in Sec. V.

##  II Model and methods

The Hamiltonian for a two-band material can be written in terms of Pauli
matrices $\boldsymbol{\tau}=(\tau_{x},\tau_{y},\tau_{z})$ as follows:

$\displaystyle H=\boldsymbol{\tau}\cdot\mathbf{d},$ (1)

where the matrices work on the vector $|\psi\rangle=(\psi_{A},\psi_{B})$
consisting of the wave function of each site in a lattice. Around the $K$ and
$K^{\prime}$ points of a gapped graphene-like materials , $\mathbf{d}$ is
given by:

$\displaystyle\mathbf{d}=(\chi
vp_{x}+\lambda(p_{y}^{2}-p_{x}^{2}),vp_{y}+2\chi\lambda p_{x}p_{y},\Delta),$
(2)

where $\Delta$ is half of the band gap, $\chi=+1$ ($\chi=-1$) is the chirality
around the $K$-point ($K^{\prime}$-point), and $v$ along with $\lambda$ are
constants determined by the material’s structure. Expressed in the tight-
binding hopping parameter $t$ and the lattice constant $a$, one can write $v$
and $\lambda$ as $v=\sqrt{3}ta/2\hbar$ and $\lambda=ta^{2}/8\hbar^{2}$,
respectively. The term quadratic in $p$ generates a trigonal warping of the
rotationally-symmetric Dirac energy dispersion, in which the $\lambda$
parameter determines the strength or degree of the warping.

By diagonalizing Eq. (1), we obtain the energy dispersion of our system:

$\displaystyle\epsilon(\mathbf{p})=$ $\displaystyle\pm|\mathbf{d}|$
$\displaystyle=$
$\displaystyle\sqrt{\Delta^{2}+v^{2}p^{2}+2\chi\lambda(3p_{y}^{2}-p_{x}^{2})p_{x}+\lambda^{2}p^{4}}.$
(3)

In Fig. 1, we plot the energy bands and each corresponding energy contour. In
Fig. 1.(a) we plot the energy bands of a conventional Dirac cone ($\lambda=0$)
while in Fig. 1.(b) we include the trigonal warping terms with
$\lambda\Delta/v^{2}=0.05$. We observe in Fig.1.(b) that with nonzero
$\lambda$ the bands no longer form a single symmetric (gapped) Dirac cones,
but are warped into a triangular pattern with multiple extrema. From the
contour, we can observe the presence of three additional extrema with angular
spacing at $2\pi/3$ around the central point. The distance (in $p$) between
the central point and the secondary Dirac points can be found to be
$v^{2}/\lambda\Delta$.

Figure 1: Surface plot of the energy dispersion for the (a) conventional Dirac
material and (b) Dirac material with trigonal warping, expressed in the
normalized quantities $P=vp/\Delta$ and $d/\Delta$. We take $L=0.05$ for the
trigonally warped Dirac material.

##  III Linear conductivity

To calculate the linear conductivity of our system, we use the pseudospin
approach previously described in [10]. The dynamics of the electron
wavefunction can be described in terms of
$\mathbf{m}=\langle\psi|\boldsymbol{\tau}|\psi\rangle$. The time evolution of
$\mathbf{m}$ can be derived from the Bloch equation of motion,

$\displaystyle\frac{d\mathbf{m}}{dt}=\left\langle\psi\bigg{|}\frac{1}{i\hbar}[\boldsymbol{\tau},H]\bigg{|}\psi\right\rangle=\frac{2}{\hbar}\mathbf{d}\times\mathbf{m}.$
(4)

We can consider an external electric field $\mathbf{E}=E_{0}e^{-i\omega
t}\mathbf{\hat{x}}$ as a perturbation on $\mathbf{d}$ and expand $\mathbf{d}$
up to the first order as $\mathbf{d}^{(0)}+\delta\mathbf{d}$ where
$\delta\mathbf{d}=\frac{ieE_{0}}{\omega}(\chi v-2\lambda p_{x},2\chi\lambda
p_{y},0)$. Then the perturbation in $\mathbf{m}$ is given by

$\displaystyle\frac{d}{dt}\delta\mathbf{m}=-\frac{2}{\hbar}(\delta\mathbf{m}\times\mathbf{d}+\mathbf{m}^{(0)}\times\delta\mathbf{d}),$
(5)

where $\mathbf{m}^{(0)}=-\frac{|\mathbf{d}|}{d}$ is the equilibrium value of
$\mathbf{m}$.

Let us take $\delta\mathbf{m}=(\delta m_{x},\delta m_{y},\delta
m_{z})e^{-i\omega t}$. Thus, we can solve Eq. (5) for $(\delta m_{x},\delta
m_{y},\delta m_{z})$ as follows:

$\displaystyle\delta m_{x}$
$\displaystyle=\frac{1}{id((\hbar\omega/2)^{2}-d^{2})}\left(-i(d_{y}^{2}+d_{z}^{2})\delta
d_{x}-((\hbar\omega/2)d_{z}-id_{x}d_{y})\delta
d_{y}+((\hbar\omega/2)d_{y}+id_{x}d_{z})\delta d_{z}\right),$ (6)
$\displaystyle\delta m_{y}$
$\displaystyle=\frac{1}{id((\hbar\omega/2)^{2}-d^{2})}\left(((\hbar\omega/2)d_{z}+id_{x}d_{y})\delta
d_{x}-i(d_{x}^{2}+d_{z}^{2})\delta
d_{y}-((\hbar\omega/2)d_{x}-id_{y}d_{z})\delta d_{z}\right),$ (7)
$\displaystyle\delta m_{z}$
$\displaystyle=\frac{1}{id((\hbar\omega/2)^{2}-d^{2})}\left(((\hbar\omega/2)d_{y}-id_{x}d_{z})\delta
d_{x}+((\hbar\omega/2)d_{x}+id_{y}d_{z})\delta
d_{y}-i(d_{x}^{2}+d_{y}^{2})\delta d_{z}\right),$ (8)

The electron current can then be calculated as follows:

$\mathbf{j}=\sum_{\mathbf{p}}e\left\langle\frac{\partial
H}{\partial\mathbf{p}}\right\rangle=\sum_{\mathbf{p}}e\frac{\partial
d_{x}}{\partial\mathbf{p}}\delta m_{x}+\frac{\partial
d_{y}}{\partial\mathbf{p}}\delta m_{y}.$ (9)

By substituting Eqs. (6)) - (8) to Eq. (9) and using $j_{i}=\sigma_{ij}E_{j}$,
we obtain the Hall conductivity as follows:

$\displaystyle\sigma_{xy}=\frac{\chi e^{2}}{4\pi^{2}\hbar^{2}}\int
d^{2}p\frac{(\hbar\omega/2)\Delta}{d\omega((\hbar\omega/2)^{2}-d^{2})}(-4\lambda^{2}p^{2}+v^{2}),$
(10)

where we change the summation on $\mathbf{p}$ to integration and we define
$d\equiv\sqrt{\Delta^{2}+v^{2}p^{2}+2\chi\lambda
v(3p_{y}^{2}-p_{x}^{2})p_{x}+\lambda^{2}p^{4}}$. Furthermore, we define
$P\equiv vp/\Delta$ and $L\equiv\lambda\Delta/v$ to normalize the integral and
we obtain,

$\displaystyle\sigma_{xy}=\frac{\chi e^{2}}{8\pi^{2}\hbar}\int
d^{2}\mathbf{P}\frac{(-4L^{2}P^{2}+1)}{D(\Omega^{2}-D^{2})},$ (11)

where $D\equiv\sqrt{1+P^{2}+2\chi L(3P_{y}^{2}-P_{x}^{2})P_{x}+L^{2}P^{4}}$
and $\Omega=\hbar\omega/(2\Delta)$. We note that

$\mathcal{F}_{xy}=\frac{\chi\hbar^{2}\Delta(-4\lambda^{2}p^{2}+v_{F}^{2})}{d^{3}},$
(12)

is the Berry curvature, such that for $\omega=0$ we see that the Hall
conductivity is given by $e^{2}C/h$, where
$C=(1/4\pi\hbar^{2})\int\mathcal{F}_{xy}d^{2}\mathbf{p}$ is the Chern number.

By using similar method, we obtain the longitudinal conductivity $\sigma_{xx}$
as follows:

$\displaystyle\sigma_{xx}$
$\displaystyle=\frac{-ie^{2}}{8\pi^{2}\hbar}\int\frac{(1+P_{y}^{2}+4L^{2}P^{2}+4L^{2}P_{y}^{2}(4P_{x}^{2}+P_{y}^{2})+4L^{4}P_{y}^{2}P^{4})-4\chi
LP_{x}(1+P_{y}^{2}-2L^{2}P_{y}^{2}P^{2})}{D\Omega(\Omega^{2}-D^{2})}d^{2}P.$
(13)

It is noted that due to the symmetry of the band structure under a
$120^{\circ}$ rotation, the longitudinal conductivity is isotropic, in
particular $\sigma_{xx}=\sigma_{yy}$. The integrands in the $\sigma_{xx}$ and
$\sigma_{yx}$ have a singularity at $D=\pm\Omega$ whenever $|\hbar\omega|\geq
2\Delta$. To remove the singularity during the calculation of conductivity we
replace $\Omega$ to $\Omega-i\eta$ and taking the limit $\eta\rightarrow 0$.
Then, by using the Sokhostski-Plemelj relation, we obtain the imaginary (real)
part of $\sigma_{xy}(\Omega)$ ($\sigma_{xx}$ and $\sigma_{yy}$) . To obtain
the real (imaginary) part of $\sigma_{xy}$ ($\sigma_{xx}$ and $\sigma_{yy}$)
we can use the Kramers-Kronig relations

$\displaystyle\textrm{Re}[\sigma_{xy}(\Omega)]$
$\displaystyle=\frac{2}{\pi}\int_{0}^{\infty}\frac{\Omega^{\prime}\textrm{Im}[\sigma_{xy}(\Omega^{\prime})]}{\Omega^{\prime
2}-\Omega^{2}}d\Omega^{\prime}$ (14)
$\displaystyle\textrm{Im}[\sigma_{xx}(\Omega)]$
$\displaystyle=-\frac{2}{\pi}\int_{0}^{\infty}\frac{\Omega\textrm{Re}[\sigma_{xx}(\Omega^{\prime})]}{\Omega^{\prime
2}-\Omega^{2}}d\Omega^{\prime}$ (15)

The delta function integrals, and the Kramers-Kronig integrals, can then be
calculated numerically. We used Python’s SciPy library [17] for the numerical
integrations.

For our numerical calculation, we set $L=\Delta/6t$ to $L=0.05$, and assume
the bandgap $\Delta$ to be nonzero. We plot the linear conductivity of our
model as a function of normalized frequency $\Omega=\hbar\omega/2\Delta$ in
Fig. 2. For comparison, we plot the conductivity of a conventional Dirac
material with $L=0$ in Fig. 2(a) and the conductivity of a Dirac material with
trigonal warping material in Figs. 2(b)–(c). We take the warping coefficient
$L=0.05$, which roughly corresponds to a band gap of
$\Delta=0.5~{}\mathrm{eV}$ and a first neighbour hopping integral similar to
graphene, $t=2.76~{}\mathrm{eV}$. We note from Eq. (10) that in case of zero
bandgap, the Hall conductivity will vanish.

Figure 2: (a) The linear optical conductivity for conventional Dirac material
with $L=0$. (b) The Hall conductivity and (c) The longitudinal conductivity
for Dirac material with trigonal warping ($L=0.05$). The $x$ axis shows the
normalized frequency $\Omega=\hbar\omega/2\Delta$.

From these results, we see that the DC Hall conductivity for our model is
twice higher than the DC Hall conductivity for a Dirac material, meaning the
Chern number of our trigonal warping band structure is $C=\pm 1$, in contrast
to the conventional $C=\pm 1/2$ Chern number for a single Dirac cone. This is
due to the presence of three extra Dirac cone around the central Dirac point
created by the trigonal warping term (see Fig. 1). By integrating the Berry
curvature around each of these Dirac cones, it can be seen that the central
cone carries a $C=\pm 1/2$ Chern number while the secondary cones carry a
$C^{\prime}=\mp 1/2$ Chern number, yielding a total of $\mp 1$ overall Chern
number. This presence of additional Dirac points is reminiscent of the model
examined in Ref.[18] concerning a Haldane-model like material with third and
higher-neighbor hoppings. However, we must note that here we only consider a
single Dirac point, while generally in materials with intact TR symmetry, each
Dirac cone is paired with another Dirac point with opposite chirality and
opposite Chern number. As a result, when accounting for all the Dirac cones in
the material, the total Chern number is zero, unless we break its TR symmetry
in which we can obtain a non-zero Chern number. Still, our result demonstrates
that trigonal warping allows the presence of a nontrivial Chern number under
the appropriate parameters. For example, in the third neighbor hopping model
of [18] we have three additional Dirac cone in a trigonal configuration, which
give a total Chern number of $\pm 2$ for a certain range of the third-neighbor
hopping, while with $\mathrm{LaN_{2}}$, there are six Dirac cones in the
Brillouin zone, which results in a total Chern number of $C=3$ [19].

Another notable feature from Fig. 2 is the presence of a singularity at
$\Omega\approx 5$ for the longitudinal conductivity, in addition to the
singularity at resonant frequency $\hbar\omega=2\Delta$. This singularity is
related to the van Hove singularity, where the density of states diverge due
to a saddle point in the band structure. By taking the derivative of Eq. (3),
we can verify that a saddle point of $|\mathbf{d}(p)|$ occurs when
$|\mathbf{d}|=\sqrt{1+1/(16L^{2})}=5.09$ for $L=0.05$. This saddle point
agrees with previous calculations of optical conductivity [20, 21], which
shows the presence of a singularity for materials with honeycomb lattice.
Additional results regarding Kerr and Faraday angles related to this material
are given in the Appendix.

##  IV Nonlinear conductivity

To further explore the consequences of our band structure’s asymmetry, we next
calculate the nonlinear dc response of the trigonally warped Dirac material.
Following the method given by [14], we can write the current as the sum of
intraband $\mathbf{j}_{\textrm{intra}}$ and interband
$\mathbf{j}_{\textrm{inter}}$ contributions. We first note that the intraband
contribution from the anomalous photocurrent term [22] is zero for our system.
For a 2D system with normally-incident light, the anomalous photocurrent term
is reduced to the curl of the Berry curvature
$\mathbf{\mathcal{F}}_{xy}=(\partial_{k_{x}}\mathcal{A}_{y}-\partial_{k_{y}}\mathcal{A}_{x})\mathbf{\hat{z}}$,

$\mathbf{j}^{(2)}_{ap}(0)=-\frac{ie^{3}}{2h^{2}\omega}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}d^{2}\mathbf{k}(\mathbf{E}\times\mathbf{E}^{*})\times(\nabla\times\mathbf{\mathcal{F}}_{xy}).$
(17)

Expanding the integrand into component form, we see that the anomalous
photocurrent can be expressed in terms of the Berry curvature dipole (see
[12]) $D_{i}$ as follows:

$j^{(2)}_{\mathrm{ap},i}(0)=-\frac{ie^{3}(E_{x}E_{y}^{*}-E_{x}^{*}E_{y})D_{i}}{2\hbar^{2}\omega},$
(18)

where

$D_{i}=\int d^{2}\mathbf{k}~{}\partial_{i}\mathcal{F}_{xy}.$ (19)

The Berry curvature for our system is

$\mathcal{F}_{xy}=\frac{\chi\hbar^{2}\Delta(-4\lambda^{2}p^{2}+v^{2})}{(\Delta^{2}+v^{2}p^{2}-2\chi\lambda
vp^{3}\cos{3\theta}+\lambda^{2}p^{4})^{3/2}},$ (20)

which has three mirror symmetry planes at $\theta=0^{\circ},120^{\circ},$ and
$240^{\circ}$. As the derivative of the Berry curvature behaves as a
pseudovector under reflection symmetry, the Berry curvature dipole will always
be zero, guaranteeing zero photocurrent unless this symmetry is broken (such
as by an uniaxial strain as in [15]). However, unlike the Berry curvature, the
Berry connections $\mathcal{A}_{cc},\mathcal{A}_{vv}$ does not possess this
threefold symmetry. This lack of symmetry allows the presence of a nonlinear
transverse current arising from the interband shift current. Therefore, here
we are most interested in the second-order contribution from the shift
current,

$\mathbf{j}_{\textrm{inter}}^{(2)}=-\frac{e^{2}}{\hbar}\sum_{n\neq
m,\mathbf{k}}\rho_{nm}\mathbf{\mathcal{D}_{mn}}(\mathbf{E}(t)\cdot\mathbf{\mathcal{A}_{mn}}).$
(21)

Here, $\rho_{nm}$ is the density matrix element,
$\mathbf{\mathcal{A}_{mn}}=i\langle\psi_{m}|\frac{\partial}{\partial\mathbf{k}}|\psi_{n}\rangle$
is the Berry connection vector, and
$\mathbf{\mathcal{D}_{mn}}=\frac{\partial}{\partial\mathbf{k}}+i(\mathbf{\mathcal{A}_{mm}}-\mathbf{\mathcal{A}_{nn}})$.

The second-order contribution of the shift current term can be written in term
of the shift vector
$\mathbf{r}_{\textrm{shift}}=\mathbf{\mathcal{A}}_{cc}-\mathbf{\mathcal{A}}_{vv}-\partial_{\mathbf{k}}\mathrm{arg}[\mathbf{E}\cdot\mathbf{\mathcal{A}}_{vc}]$[23,
24, 25, 26],

$\mathbf{j}^{(2)}_{\textrm{{shift}}}=\frac{e^{3}}{2h}\int
d^{2}\mathbf{k}\delta(\hbar\omega-2d)|\mathbf{E}\cdot\mathbf{\mathcal{A}}_{vc}^{*}|^{2}\mathbf{r}_{\textrm{shift}}.$
(22)

For linear electric fields in the $x$ or $y$ direction, the conductivities are
then

$\sigma_{ijj}(\omega)=\frac{e^{3}}{2h}\int
d^{2}\mathbf{k}~{}\delta(\hbar\omega-2d)\mathcal{A}_{vc,j}^{*}\mathcal{A}_{vc,j}r_{\textrm{shift,i}}$
(23)

As the shift vector $\mathbf{r}_{shift}$ is always real, this equation shows
the real part of $\sigma_{ijj}$. The imaginary parts of $\sigma_{ijj}$ can be
calculated from the Kramers-Kronig relations.

We now consider the dc current generated by a linearly polarized light with
normal incident. The resulting shift current is directed to the $+y$ direction
when $\mathbf{E}$ is polarized in the $x$ direction, while if $\mathbf{E}$ is
polarized in the $y$ direction the resulting shift current is directed to the
$-y$ direction (Fig. 3).

Figure 3: The direction of the current $\mathbf{j}^{\textrm{shift}}$ and the
electric field $\mathbf{E}$ polarization superimposed on the energy contour
plot for (a) $x$-polarized electric field (b) $y$-polarized electric field.

This anisotropy in current direction is related to the underlying symmetry of
our Hamiltonian. The energy dispersion of our system is invariant with respect
to reflection over the $p_{x}$ axis ($d(p_{x},p_{y})=d(p_{x},-p_{y})$) while
the Hamiltonian itself obeys $H(p_{x},-p_{y})=H^{*}(p_{x},p_{y})$ as can be
seen by simple substitution on Eq. (2). Then, this symmetry imposes the
following condition on the eigenvalues $|u_{m}(k_{x},k_{y})\rangle$

$\displaystyle|u_{m}(k_{x},-k_{y})\rangle$
$\displaystyle=|u_{m}(k_{x},k_{y})\rangle^{*},$ (24)

where $m=c$ for the conduction band and $m=v$ for the valence band. The Berry
connections then obey the following conditions,

$\displaystyle(\mathcal{A}_{mn})_{x}(k_{x},-k_{y})$
$\displaystyle=-(\mathcal{A}_{mn})^{*}_{x}(k_{x},k_{y}),$ (25)
$\displaystyle(\mathcal{A}_{mn})_{y}(k_{x},-k_{y})$
$\displaystyle=(\mathcal{A}_{mn})^{*}_{y}(k_{x},k_{y}).$ (26)

As the diagonal Berry connections $\mathbf{\mathcal{A}}_{mm}$ are real,

$\displaystyle(\mathcal{A}_{mm})_{x}(k_{x},-k_{y})$
$\displaystyle=-(\mathcal{A}_{mm})_{x}(k_{x},k_{y}),$ (27)
$\displaystyle(\mathcal{A}_{mm})_{y}(k_{x},-k_{y})$
$\displaystyle=(\mathcal{A}_{mm})_{y}(k_{x},k_{y}),$ (28)

while $|\mathcal{A}_{vc}(k_{x},-k_{y})|=|\mathcal{A}_{vc}(k_{x},k_{y})|$ and

$\displaystyle\partial_{k_{x}}\mathrm{arg}[\mathbf{E}\cdot\mathcal{A}_{vc}](k_{x},-k_{y})$
$\displaystyle=-\partial_{k_{x}}\mathrm{arg}[\mathbf{E}\cdot\mathcal{A}_{vc}](k_{x},k_{y}),$
(29)
$\displaystyle\partial_{k_{y}}\mathrm{arg}[\mathbf{E}\cdot\mathcal{A}_{vc}](k_{x},-k_{y})$
$\displaystyle=\partial_{k_{y}}\mathrm{arg}[\mathbf{E}\cdot\mathcal{A}_{vc}](k_{x},k_{y}).$
(30)

Thus we can conclude that
$(\mathbf{r}_{\textrm{shift}}(k_{x},-k_{y}))_{x}=-(\mathbf{r}_{\textrm{shift}}(k_{x},k_{y}))_{x}$
while
$(\mathbf{r}_{\textrm{shift}}(k_{x},-k_{y}))_{y}=(\mathbf{r}_{\textrm{shift}}(k_{x},k_{y}))_{y}$.
As a result, with integration over $k_{y}$, the $x$ component of the shift
current vanishes. Note that the addition of the trigonal warping terms to the
Dirac Hamiltonian breaks a similar symmetry with respect to $k_{x}\rightarrow-
k_{x}$, allowing a nonzero shift current in the $y$ direction.

We plot the shift conductivity in Fig. 4. There are several notable qualities
from this plot. First, we see that $\sigma_{yxx}=-\sigma_{yyy}$, as required
by the $C_{3}$ symmetry of the system [20]. Also, since the shift current is
created by interband transitions only, we see that for frequencies below the
band gap the real part of the conductivity is zero, both for the longitudinal
($\sigma_{yyy}$) and transversal ($\sigma_{yxx}$) current. A singular peak is
observed when $\hbar\omega=2\Delta$, as in the linear case. We can also check
the effects of different chirality $\chi=\pm 1$ on the nonlinear conductivity.
If we consider the remaining nonzero shift conductivity,
$\sigma_{yyy}^{\textrm{shift}}=-\sigma_{yxx}^{\textrm{shift}}$, all the Berry
connections on the integrand,
$\mathcal{A}_{ij,y}=i\langle\psi_{i}|\partial_{k_{y}}|\psi_{j}\rangle$, depend
on $\chi$ only through $|\psi_{i}\rangle$, which can be written in terms of
the two-band Hamiltonian $\mathbf{d}$ vector as [7]

$\displaystyle|\psi_{v}\rangle$
$\displaystyle=\frac{1}{\sqrt{2d(d+d_{z})}}\begin{bmatrix}d_{z}+d\\\
d_{x}-id_{y}\end{bmatrix}$ (31) $\displaystyle|\psi_{c}\rangle$
$\displaystyle=\frac{1}{\sqrt{2d(d-d_{z})}}\begin{bmatrix}d_{z}-d\\\
d_{x}-id_{y}\end{bmatrix}$ (32)

We see from the Hamiltonian that the chirality $\chi$ is always coupled with
$p_{x}$, thus flipping the chirality is equivalent to a substitution
$p_{x}\rightarrow-p_{x}$ to the Hamiltonian. So, we have
$\mathcal{A}_{ij,y}^{\chi=-1}(k_{x},k_{y})=\mathcal{A}_{ij,y}^{\chi=+1}(-k_{x},k_{y})$.
As flipping the chirality is equivalent to the substitution $p_{x}\rightarrow-
p_{x}$, while Eq. (22) is already integrated over $p_{x}\in(-\infty,\infty)$,
we have
$\sigma_{ijj}^{\textrm{shift},\chi=+1}=\sigma_{ijj}^{\textrm{shift},\chi=-1}$,
which has also been checked manually through numerical integration. As a
result, this shift current is nonzero even with TR symmetry, in contrast to
the linear Hall conductivity [27], which vanishes when summed over
$\chi=+1,-1$ unless TR symmetry is broken.

Figure 4: The nonlinear transverse conductivity (in dimensionless units) as a
function of normalized light frequency $\hbar\omega/2\Delta$. We plot the
second order contribution from the shift current, stimulated by linearly
polarized light in the $x$ direction (red) and in the $y$ direction (blue)

To discover possible realizations for our system, we calculated the shift
current of a monolayer $\mathrm{MoS}_{2}$. $\mathrm{MoS}_{2}$ possesses a
hexagonal lattice, and it has been noted that the band structure around the K
and K’ points exhibits trigonal warping [28]. By fitting our Hamiltonian (Eq.
(1)) to the band structure of $\mathrm{MoS}_{2}$ around the K point, we obtain
the following value of the band parameters: $\Delta=0.89~{}\mathrm{eV}$,
$v=0.22~{}\mathrm{nm/fs}$ and $L=\lambda\Delta/v^{2}=0.065$. We then calculate
the shift current of $\mathrm{MoS}_{2}$ using these parameters and Eq. (22)
and also from a first principle approach using Wannier90. We plot the
calculated shift current of $\mathrm{MoS}_{2}$ in Fig.5. In Fig.5.(a) we
compare our theoretical result from Eq.(23) with the DFT results, at
frequencies around the band gap, where the effects of multiband interactions
are negligible, and in Fig.5.(b) we show the DFT results for the full
frequency range. We note that as the DFT software calculates the 3D
conductivity, to fit our theoretical results with the DFT calculations’ units
we divide Eq.(23) with the confinement length, which we take to be 12 nm. The
first principle results diverges from the theoretical results for frequencies
further from the band gap, with the shift current taking negative values for
$\hbar\omega>2.3~{}\mathrm{eV}$. However, we emphasize that in our
calculations we assumed a two band model, while we believe this divergence to
be a result of multiband interaction, which are activated at higher
frequencies. We note that Eq. (22) does allow a negative value of the shift
current as $r_{{\rm shift},i}$ might take negative values, depending on
$\mathcal{A}_{cc}-\mathcal{A}_{vv}$.

Figure 5: The shift current in $\mathrm{MoS_{2}}$ as a function of frequency.
(a) We compare the results of the first-principles approach (solid lines) with
our theoretical approach (dotted lines) at frequencies near the band gap,
where interband effects are negligible (b) The first principle calculation of
the shift current for the full frequency range.

##  V Conclusions

We have shown that around a Dirac point with trigonal warping, the appearance
of three additional Dirac valleys with Chern number opposite to the central
Dirac point unlocks the possibility of an anomalous Hall material with higher
Chern number from multi-valley contributions.. Furthermore, by breaking
rotational symmetry we show the existence of a nonlinear transverse current
with a zero Berry curvature dipole even in the absence of broken TR symmetry.
These results provide a new perspective and directions for exploring
topological materials.

###### Acknowledgements.

We acknowledge HPC BRIN for their high-performance computing facilities.
J.M.A. and B.E.G. are supported by P2MI FMIPA ITB and P2MI KK FTETi ITB.
E.H.H. acknowledges supports from the National Research Fund Luxembourg under
grants CORE C20/MS/14764976/TOPREL and C21/MS/15752388/NAVSQM.

##  Appendix A Kerr and Faraday angles

In addition to our main result above, to illuminate the possible methods for
the characterization of warping effects, we also plot the resulting Kerr and
Faraday angles produced by a linearly polarized light normally incident on a
thin film of our anomalous Hall material. For a 2D material on a substrate
with refraction index $n$ expression for Kerr angle and ellipticity is given
by [21]

$\displaystyle\tan{2\theta_{K}}$ $\displaystyle\sim
2\mathrm{Re}\left[\frac{8\pi\sigma_{yx}}{c(n^{2}-1)}\right]$ (33)
$\displaystyle\sin{2\eta_{K}}$ $\displaystyle\sim
2\mathrm{Im}\left[\frac{8\pi\sigma_{yx}}{c(n^{2}-1)}\right]$ (34)

while the Faraday angles is approximately one half of the Kerr angles. We plot
the results in Fig. 6.

Figure 6: (a) Kerr angle, (b) Kerr ellipticity, (c) Faraday angle, and (d)
Faraday ellipticity of a linearly polarized light that is incident in the
normal direction to the material. All properties are plotted as a function of
normalized frequency. We take the refractive index of the substrate to be
$n=2$.

##  Appendix B First-principles calculations

We carried out supporting first-principles, density functional theory (DFT)
calculations using the Quantum ESPRESSO package [29]. The monolayer 2H-MoS2
lattice parameters are adopted from previous calculations [30] with vacuum
space set to 20 Å in $z$-direction. We utilized generalized gradient
approximation exchange-correlation functional as implemented by Perdew, Burke,
and Ernzerhof (PBE-GGA)[31] and full-relativistic norm-conserving
pseudopotentials [32]. The self-consistent calculations are found optimally
converged with basis set energy cutoff of 35 Ry and 16 $\times$ 16 $\times$ 1
Monkorst pack k-mesh grid. We constructed maximally localized Wannier
functions (MLWF) [33, 34] using WANNIER90 package[35], projected from both
Mo-d and S-p orbitals. Shift current properties were computed using Wannier
interpolation scheme [36] with 500 $\times$ 500 $\times$ 1 integration k-grid
tested for convergence.

## References

  * Laughlin [1981] R. B. Laughlin, Quantized hall conductivity in two dimensions, Phys. Rev. B 23, 5632 (1981).
  * Haldane [1988] F. D. M. Haldane, Model for a quantum hall effect without landau levels: Condensed-matter realization of the ”parity anomaly”, Phys. Rev. Lett. 61, 2015 (1988).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * Thouless  _et al._ [1982] 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).
  * Qi  _et al._ [2008] X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Topological field theory of time-reversal invariant insulators, Phys. Rev. B 78, 195424 (2008).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010).
  * Bernevig and Hughes [2013] B. Bernevig and T. Hughes,  _Topological Insulators and Topological Semiconductors_ , 1st ed. (Princeton University Press, 2013).
  * Yang  _et al._ [2012] S. Yang, Z.-C. Gu, K. Sun, and S. Das Sarma, Topological flat band models with arbitrary Chern numbers, Phys. Rev. B 86, 241112(R) (2012).
  * Zhao  _et al._ [2020] Y.-F. Zhao, R. Zhang, R. Mei, L.-J. Zhou, H. Yi, Y.-Q. Zhang, J. Yu, R. Xiao, K. Wang, N. Samarth, M. H. W. Chan, C.-X. Liu, and C.-Z. Chang, Tuning the Chern number in quantum anomalous Hall insulators, Nature 588, 419 (2020).
  * Hasdeo  _et al._ [2019] E. H. Hasdeo, A. J. Frenzel, and J. C. W. Song, Cyclotron motion without magnetic field, New J. Phys. 21, 083026 (2019).
  * Tokura  _et al._ [2019] Y. Tokura, K. Yasuda, and A. Tsukazaki, Magnetic topological insulators, Nat. Rev. Phys. 1, 126 (2019).
  * Sodemann and Fu [2015] I. Sodemann and L. Fu, Quantum nonlinear hall effect induced by berry curvature dipole in time-reversal invariant materials, Phys. Rev. Lett. 115, 216806 (2015).
  * Nandy and Sodemann [2019] S. Nandy and I. Sodemann, Symmetry and quantum kinetics of the nonlinear hall effect, Phys. Rev. B 100, 195117 10.1103/PhysRevB.100.195117 (2019), cited By 43.
  * Taghizadeh and Pedersen [2019] A. Taghizadeh and T. G. Pedersen, Nonlinear excitonic spin Hall effect in monolayer transition metal dichalcogenides, 2D Mater. 7, 015003 (2019).
  * Battilomo  _et al._ [2019] R. Battilomo, N. Scopigno, and C. Ortix, Berry curvature dipole in strained graphene: A fermi surface warping effect, Phys. Rev. Lett. 123, 196403 (2019).
  * Matsyshyn and Sodemann [2019] O. Matsyshyn and I. Sodemann, Nonlinear hall acceleration and the quantum rectification sum rule, Phys. Rev. Lett. 123, 246602 (2019).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nat. Methods 17, 261 (2020).
  * Sticlet and Piéchon [2013] D. Sticlet and F. Piéchon, Distant-neighbor hopping in graphene and haldane models, Phys. Rev. B 87, 115402 (2013).
  * Li  _et al._ [2020] L. Li, X. Kong, X. Chen, J. Li, B. Sanyal, and F. M. Peeters, Monolayer 1t-lan2: Dirac spin-gapless semiconductor of p-state and chern insulator with a high chern number, Appl. Phys. Lett. 117, 143101 (2020).
  * Hipolito  _et al._ [2016] F. Hipolito, T. G. Pedersen, and V. M. Pereira, Nonlinear photocurrents in two-dimensional systems based on graphene and boron nitride, Phys. Rev. B 94, 045434 (2016).
  * Pratama  _et al._ [2020] F. R. Pratama, M. S. Ukhtary, and R. Saito, Circular dichroism and faraday and kerr rotation in two-dimensional materials with intrinsic hall conductivities, Phys. Rev. B 101, 045426 (2020).
  * Rostami and Polini [2018] H. Rostami and M. Polini, Nonlinear anomalous photocurrents in Weyl semimetals, Phys. Rev. B 97, 195151 (2018).
  * Sipe and Shkrebtii [2000] J. E. Sipe and A. I. Shkrebtii, Second-order optical response in semiconductors, Phys. Rev. B 61, 5337 (2000).
  * Sinitsyn  _et al._ [2006] N. A. Sinitsyn, Q. Niu, and A. H. MacDonald, Coordinate shift in the semiclassical boltzmann equation and the anomalous hall effect, Phys. Rev. B 73, 075318 (2006).
  * Ahn  _et al._ [2020] J. Ahn, G.-Y. Guo, and N. Nagaosa, Low-frequency divergence and quantum geometry of the bulk photovoltaic effect in topological semimetals, Phys. Rev. X 10, 041041 (2020).
  * Shi  _et al._ [2021] L.-K. Shi, D. Zhang, K. Chang, and J. C. W. Song, Geometric photon-drag effect and nonlinear shift current in centrosymmetric crystals, Phys. Rev. Lett. 126, 197402 (2021).
  * Watanabe and Yanase [2021] H. Watanabe and Y. Yanase, Chiral photocurrent in parity-violating magnet and enhanced response in topological antiferromagnet, Phys. Rev. X 11, 011001 (2021).
  * Kormányos  _et al._ [2013] A. Kormányos, V. Zólyomi, N. D. Drummond, P. Rakyta, G. Burkard, and V. I. Fal’ko, Monolayer mos2: Trigonal warping, the $\Gamma$ valley, and spin-orbit coupling effects, Phys. Rev. B 88, 045416 (2013).
  * 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. D. 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, Quantum espresso: a modular and open-source software project for quantum simulations of materials, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Schankler  _et al._ [2021] A. M. Schankler, L. Gao, and A. M. Rappe, Large bulk piezophotovoltaic effect of monolayer 2h-mos2, The Journal of Physical Chemistry Letters 12, 1244 (2021), pMID: 33497221.
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Hamann [2013] D. R. Hamann, Optimized norm-conserving vanderbilt pseudopotentials, Phys. Rev. B 88, 085117 (2013).
  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Maximally localized generalized wannier functions for composite energy bands, Phys. Rev. B 56, 12847 (1997).
  * Souza  _et al._ [2001] I. Souza, N. Marzari, and D. Vanderbilt, Maximally localized wannier functions for entangled energy bands, Phys. Rev. B 65, 035109 (2001).
  * Mostofi  _et al._ [2008] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, wannier90: A tool for obtaining maximally-localised wannier functions, Computer Physics Communications 178, 685 (2008).
  * Ibañez Azpiroz  _et al._ [2018] J. Ibañez-Azpiroz, S. S. Tsirkin, and I. Souza, Ab initio calculation of the shift photocurrent by wannier interpolation, Phys. Rev. B 97, 245143 (2018).