#  Axion Electrodynamics in Topological Insulators for beginners.
Part I: Introduction and basic equations
Part II: Images of a charge close to an interface ordinary insulator-
topological insulator
Part III: Appendices on boundary conditions, the electrostatic image method
and on electrodynamic atomic units

Josep Planelles, Dept. Química-Física i Analítica, Universitat Jaume I

##  1 Introduction

According to band theory, insulators are defined as those materials that have
a gap between the occupied valence and the vacant conduction bands, while
conductors do not have a gap (either, because the last band is semi-occupied
or because the conduction and valence bands overlap).

In 2006 unusual electromagnetic properties were observed in a $CdTe|HgTe|CdTe$
quantum well:[1] this quantum well behaves in bulk as an insulator but it was
observed electric current across the interface, i.e., it behaves like a
conductor in surface. The most striking feature is that this behavior cannot
be accounted by the Maxwell’s equations.

Earlier, in 1987, Frank Wilczek[2] suggested the possible behavior of this
kind of materials and pointed out that it could be explained by means the
axion electrodynamics that himself[3] and Steven Weinberg[4] developed to
understand the violation of combined symmetries of charge conjugation and
parity in the strong interactions. The name axion is related to the name of
the particle associated to this peculiar field.

In a relatively recent paper, Qi and Zhang[5] present and discuss experimental
results on the $CdTe|HgTe|CdTe$ quantum well that shows an almost infinite
resistance (i.e. behaves like an insulator) if the $HgTe$ thickness is smaller
than a critical distance $d\sim 6.5nm$, whereas its resistance is small and
displays the typical plateaus of the quantum Hall effect for $d>6.5nm$.

Axion electrodynamics can account for this behavior. To this end, a $\theta$
parameter, related to the Berry phase and the Chern number, is introduced.
This parameter, called magneto-electric polarizability, is a piece-wise
constant function. Its value is $\theta=0$ for ordinary and $\theta=\pi$ for
time-reversal symmetry topological insulators, as e.g. $HgTe$. Axion
electrodynamics with $\theta(x,t)=2{\bf b\cdot x}-2b_{0}t$ also describes Weyl
semimetals[6] or, in general, the electrodynamics of magneto-electric
media.[7] A constant axion angle, $\theta(x,t)=k$, implies spatial and
temporal translation symmetry conservation.[8]

##  2 $\theta$-electrodynamics

Electromagnetism in material media is described by Maxwell equations. In
differential form and a.u. they can be written as (see a comment on
electrodynamics atomic units in appendix 3):

$\begin{array}[]{lll}{\boldsymbol{\nabla}}\cdot{\bf
D}=4\pi\rho\;{\rm\;\;\;(Gauss\,law)}&&{\boldsymbol{\nabla}}\times{\bf
E}=-\frac{1}{c}\frac{\partial{\bf B}}{\partial
t}{\rm\;\;\;\;\;\;\;\;\;(Faraday\,law)}\\\ \\\ {\boldsymbol{\nabla}}\cdot{\bf
B}=0\;\;\;\;\;\;\;\;{\rm(Gauss\,law\,for\,magnetism)}&&{\bf\nabla}\times{\bf
H}=\frac{1}{c}\frac{\partial{\bf D}}{\partial t}+\frac{4\pi}{c}{\bf
J}\;{\rm(Ampere\,law)}\\\ \end{array}$

We should add to these equations the ${\bf D},{\bf H}$ constitutive relations
in terms of ${\bf E},{\bf B}$. For linear material media these relations are
simple: ${\bf D}=\epsilon{\bf E}$, ${\bf H}=\frac{{\bf B}}{\mu}$, with
$\epsilon$, $\mu$ the dielectric constant and magnetic permeability. For
isotropic materials $\epsilon$ and $\mu$ are constants, while in anisotropic
materials are tensors, eventually coordinate-dependent.

In topological media, where magneto-electric effects take place, the electric
displacement vector ${\bf D}$ is modified by the magnetic induction ${\bf B}$
and the magnetic field intensity ${\bf H}$ of is in turn influenced by the
electric field. Then, the relations ${\bf D}=\epsilon{\bf E}$, ${\bf
B}=\mu{\bf H}$ must be modified:[9]

${\bf D}=\epsilon{\bf E}-\frac{\theta\alpha}{\pi}{\bf B}\;\;\;\;\;\;\;\;\;{\bf
H}=\frac{{\bf B}}{\mu}+\frac{\theta\alpha}{\pi}{\bf E}$

where $\alpha=1/137$ is the fine-structure constant and $\theta$ an additional
parameter that can be considered at the same level as the permittivity
$\epsilon$ or the permeability $\mu$. In ordinary media $\theta=0$ and we
recover Maxwell equations.

These modifications in the ${\bf D},{\bf H}$ constitutive relations entail
changes of two Maxwell equations, as we will show later. But first we will try
to outline on the natural appearance of this extra term in electromagnetism.
To this end, lets recall the electromagnetic energy expression:

$W=\frac{1}{2}({\bf E}\cdot{\bf D}+{\bf B}\cdot{\bf H})=\frac{1}{2}\epsilon
E^{2}+\frac{1}{2}\frac{1}{\mu}B^{2}$

On the other hand, the Lagrangian is:[10]

${\cal L}=\frac{1}{2}\epsilon E^{2}-\frac{1}{2}\frac{1}{\mu}B^{2}$

From this Lagrangian, Maxwell equations can be obtained by means the Euler-
Lagrange variational calculus (see e.g. Civelek et al.[11]).

We observe that the Lagrangian is quadratic with respect to the
electromagnetic field. F. Wilczek[2] pointed out that there exists an
additional quadratic term missing in the previous Lagrangian:

$\Delta{\cal L}=\kappa\,\theta\,{\bf E}\cdot{\bf B}$

By including this term, Wilczek obtains the axion electrodynamics equations.
Here, however, we pursue a less elegant but simpler derivation: we
incorporate, as already said, the ${\bf D},{\bf H}$ constitutive relations in
a $\theta\neq 0$ medium. By looking to Maxwell equations we can see:

$\displaystyle{\boldsymbol{\nabla}}\cdot{\bf D}$ $\displaystyle=4\pi\rho$
$\displaystyle{\rm(changes\,as\,D\,changes)}$ (1)
$\displaystyle{\boldsymbol{\nabla}}\times{\bf E}$
$\displaystyle=-\frac{1}{c}\frac{\partial{\bf B}}{\partial t}$
$\displaystyle{\rm(does\,not\,change)}$ (2)
$\displaystyle{\boldsymbol{\nabla}}\cdot{\bf B}$ $\displaystyle=0$
$\displaystyle{\rm(does\,not\,change)}$ (3)
$\displaystyle{\bf\nabla}\times{\bf H}$
$\displaystyle=\frac{1}{c}\frac{\partial{\bf D}}{\partial
t}+\frac{4\pi}{c}{\bf J}$
$\displaystyle{\rm(changes\,as\,H\,and\,D\,changes)}$ (4)

We add now the constitutive relations:

$\displaystyle{\bf D}$ $\displaystyle=\epsilon{\bf
E}-\frac{\theta\alpha}{\pi}{\bf B}$ (5) $\displaystyle{\bf H}$
$\displaystyle=\frac{{\bf B}}{\mu}+\frac{\theta\alpha}{\pi}{\bf E}$ (6)

From equations (1), (5) we can write ${\boldsymbol{\nabla}}\cdot(\epsilon{\bf
E})=4\pi\rho+\frac{\alpha}{\pi}\;{\boldsymbol{\nabla}}\cdot(\theta{\bf B})$.
Since ${\boldsymbol{\nabla}}\cdot{\bf B}=0$, eq. (2), we finally rewrite eq.
(1) as:

${\boldsymbol{\nabla}}\cdot(\epsilon{\bf
E})=4\pi\rho+\frac{\alpha{\boldsymbol{\nabla}}\theta}{\pi}\cdot{\bf B}$ (7)

In a similar way, from equations (4), (5) and (6) we have:

${\bf\nabla}\times(\frac{1}{\mu}\,{\bf
B})-\frac{1}{c}\frac{\partial{(\epsilon{\bf E})}}{\partial
t}=\frac{4\pi}{c}{\bf
J}-\frac{1}{c}\frac{\alpha}{\pi}\;\frac{\partial(\theta{\bf B})}{\partial
t}-\frac{\alpha}{\pi}\;{\boldsymbol{\nabla}}\times(\theta\,{\bf E})$

Now we carry out the two last terms derivatives:

$\displaystyle-\frac{1}{c}\frac{\alpha}{\pi}\;\frac{\partial(\theta{\bf
B})}{\partial t}-\frac{\alpha}{\pi}\;{\boldsymbol{\nabla}}\times(\theta{\bf
E})$
$\displaystyle=-\frac{1}{c}\frac{\alpha}{\pi}\;\left(\frac{\partial\theta}{\partial
t}\,{\bf B}+\theta\,\frac{\partial{\bf B}}{\partial
t}\right)-\frac{\alpha}{\pi}\;\left({\boldsymbol{\nabla}}\theta\times{\bf
E}+\theta\;{\boldsymbol{\nabla}}\times{\bf E}\right)$
$\displaystyle=\frac{\alpha}{\pi}\,\theta\,\left(-\frac{1}{c}\frac{\partial{\bf
B}}{\partial t}-{\boldsymbol{\nabla}}\times{\bf
E}\right)-\frac{\alpha}{\pi}\,\left(\frac{1}{c}\frac{\partial\theta}{\partial
t}\,{\bf B}+{\boldsymbol{\nabla}}\theta\times{\bf E}\right)$

Taking into account eq (2), the first bracket in the last equation must be
zero. Then, eq (4) turns into:

${\bf\nabla}\times(\frac{1}{\mu}\,{\bf
B})-\frac{1}{c}\frac{\partial{(\epsilon{\bf E})}}{\partial
t}=\frac{4\pi}{c}{\bf
J}-\frac{1}{c}\frac{\alpha}{\pi}\frac{\partial\theta}{\partial t}\,{\bf
B}-\frac{\alpha}{\pi}\,{\boldsymbol{\nabla}}\theta\times{\bf E}$ (8)

The new equations (7) and (8), replacing eqs. (1) and (4), suggest the
definition of effective charges and current densities. Thus, from (7) we
define the effective charge $\rho_{\theta}$:

$4\pi\rho_{\theta}=\frac{\alpha{\boldsymbol{\nabla}}\theta}{\pi}\cdot{\bf
B}\to\boxed{\rho_{\theta}=\frac{\alpha}{4\pi^{2}}\,{\boldsymbol{\nabla}}\theta\cdot{\bf
B}}\\\ $ (9)

and from (8) the effective current density ${\bf J}_{\theta}$:

$\frac{4\pi}{c}{\bf
J}_{\theta}=\frac{\alpha}{\pi}\left(-\frac{1}{c}\frac{\partial\theta}{\partial
t}\,{\bf B}-{\boldsymbol{\nabla}}\theta\times{\bf E}\right)\to\boxed{{\bf
J}_{\theta}=\frac{\alpha}{4\pi^{2}}\,\left(-\frac{\partial\theta}{\partial
t}\,{\bf B}-c\,{\boldsymbol{\nabla}}\theta\times{\bf E}\right)}\\\ $ (10)

Should $\theta=0$ then, equations (7) and (8) goes back to (1) and (4). Should
be $\theta$ coordinates and time independent, i.e. should be $\theta$ a
constant, then $\rho_{\theta}={\bf J}_{\theta}=0$ and we return again to the
Maxwell equations.

As said above, $\theta=0$ in a vacuum or in an ordinary insulator, while
$\theta=\pi$ for time reversal topological insulators. Therefore, $\theta$ is
a piecewise constant function (see figure).

If $\theta(z)$ is a step function, then
${\boldsymbol{\nabla}}\theta=\pi\delta(z_{0}){\bf n}$, where ${\bf n}$ is a
unit vector in the $z$ direction and $\delta(z_{0})$ is the Dirac delta (see
note[12]). Therefore, since $\theta$ is time-independent, eq. (10) leads to:

${\bf J}_{\theta}=-\frac{c\,\alpha}{4\pi}\,\delta(z_{0})\,{\bf n}\times{\bf
E}$ (11)

In a similar way, from l’eq. (9) we have:

$\rho_{\theta}=\frac{\alpha}{4\pi}\,\delta(z_{0})\,{\bf n}\cdot{\bf B}$ (12)

The set of equations (7), (2), (3), (8) are those introduced by F. Wilczek[2]
to define axion electrodynamics and are broadly used (see e.g. equation 39 in
[13], or equation 68 in [8] etc.). However, Luca Visinelli[14] points out that
Maxwell equations for an electromagnetic field show an internal symmetry,
known as the duality transformation,

$\left(\begin{matrix}E^{\prime}\\\
B^{\prime}\end{matrix}\right)=\left(\begin{matrix}\cos\xi&\sin\xi\\\
-\sin\xi&\cos\xi\end{matrix}\right)\left(\begin{matrix}E\\\
B\end{matrix}\right)$

and that whenever a pseudoscalar axion-like field $\theta=\theta(x)$ is
introduced in the theory, the dual symmetry is spontaneously and explicitly
broken. He relates this broken symmetry to the fact that the introduction of
an axion-like interaction with the electromagnetic field only modifies two of
the four Maxwell equations (Gauss and Ampere laws), but not the remaining two
equations (Faraday and Gauss law for B). The requirement that the electric and
magnetic fields must satisfy also the above duality relation along with Gauss
and Ampere laws for axion electrodynamics leads him we obtain new terms that
also modify Faraday law and Gauss law for B (see eq. 26 in [14]. See also
[15]).

Part II:

Images of a charge close to an interface ordinary insulator-topological
insulator

##  3 Image charge in topological insulators

###  3.1 Boundary conditions at the interface with a topological insulator

In Appendix 1 section 2 we obtain the boundary conditions (BCs) for the
electric displacement ${\bf D}$ and the magnetic field intensity ${\bf H}$ in
electrostatics and magnetostatics (i.e., for time-independent fields). In
particular, ${\bf D}_{2\perp}={\bf D}_{1\perp}+\sigma$, where we implicitly
assume that the border hosting the free surface charge $\sigma$ belongs to
medium 1. At an interface between ordinary insulators there is no free charge
and therefore $\sigma=0$. If medium 1 is ordinary and medium 2 topological,
the border surface containing the charge density
$\sigma=\frac{\alpha}{4\pi}\,{\bf n}\cdot{\bf B}$, eq. (12), belongs to the
topological medium 2. Incidentally, the sign of $\sigma$ is determined by the
direction of ${\bf B}_{z}$. With $B_{z}\equiv{\bf B}_{z}(z=0)$, the boundary
condition can be rewritten as:

$\epsilon_{1}{\bf E}_{1\perp}=\epsilon_{2}{\bf E}_{2\perp}+\alpha B_{z}$ (13)

The second BC is ${\bf H}_{1\parallel}={\bf H}_{2\parallel}+\frac{4\pi}{c}{\bf
J}_{s}$. Since ${\bf J}_{s}=-\frac{c\,\alpha}{4\pi}\,{\bf n}\times{\bf E}$,
see eq. (11), it yields:

$\frac{1}{\mu_{1}}{\bf B}_{1\parallel}=\frac{1}{\mu_{2}}{\bf
B}_{2\parallel}-\alpha{\bf E}_{\parallel}$ (14)

where ${\bf E}_{\parallel}={\bf E}_{\parallel}(x,y,0)$ (we assume the border
located at $z=0$).

###  3.2 Images of an electric charge close to the interface between ordinary
and topological insulator

We assume a static problem i.e., without magnetic induction ${\bf B}$ or
electric field ${\bf E}$ temporary dependence, and that there are no charges
or currents except at the interface. Then, except at the interface, Maxwell
equations (1), (2), (3), (4) become ${\boldsymbol{\nabla}}\cdot{\bf D}=0$,
${\boldsymbol{\nabla}}\cdot{\bf B}=0$, ${\boldsymbol{\nabla}}\times{\bf E}=0$
and ${\boldsymbol{\nabla}}\times{\bf H}=0$. Should the curl of a vector field
be zero, then it can be written as the gradient of a scalar field. Therefore,
from ${\boldsymbol{\nabla}}\times{\bf E}=0$ we conclude that ${\bf
E}=-{\boldsymbol{\nabla}}V$. On the other hand, from
${\boldsymbol{\nabla}}\times{\bf H}=0$ and the constituent equation (6), it
follows that ${\boldsymbol{\nabla}}\times\left(\frac{{\bf
B}}{\mu}+\frac{\theta\alpha}{\pi}{\bf E}\right)=0$ and then $\frac{{\bf
B}}{\mu}+\frac{\theta\alpha}{\pi}{\bf E}=-{\boldsymbol{\nabla}}W$ or, in an
equivalent way, ${\bf
B}=-{\boldsymbol{\nabla}}\left[\mu\,(W-\frac{\theta\alpha}{\pi}V)\right]=-{\boldsymbol{\nabla}}U$.
Therefore, we can define electric and magnetic potentials whose opposite sign
gradients yield the field.

Let $Q$ be an electric charge located at $(0,0,a)$ in an
($\epsilon_{1},\mu_{1}$) ordinary medium in contact with an
($\epsilon_{2},\mu_{2},\alpha$) topological insulator ( $\theta=0,\pi$ for
ordinary and topological insulators, respectively). Assume the boundary at
$z=0$. According to the image method (see Appendix 1), in order to calculate
the potential in the ordinary insulator (zone 1, $z>0$) we must add a
fictitious charge at $(0,0,-a)$ while to do it in the topological insulator
(zone 2, $z<0$) we must add the fictitious charge at $(0,0,a)$. We have
learned in the above section that a magnetic field ${\bf B}$ induces a surface
charge $\sigma=\frac{\alpha}{4\pi}\,{\bf n}\cdot{\bf B}$ (not relevant now,
since no magnetic field is produced by an electric $Q$ source) while an
electric field (which is actually the case with the electric source $Q$)
generates surface currents ${\bf J}_{s}=-\frac{c\,\alpha}{4\pi}\,{\bf
n}\times{\bf E}$. The axial symmetry of the $Q$-generated electric field leads
us to conclude from ${\bf J}_{s}$ formulae that circular currents are
generated around the axis joining $Q$ with the interface (see Figure).

This surface current, generated by $Q$ located at $z>0$, originates a magnetic
field on the other side of the interface, $z<0$, proportional to the
$Q$-electric field. Therefore proportional to $\frac{\bf r}{r^{3}}$. This
magnetic field is equivalent to the magnetic field generated by a magnetic
monopole $p_{1}$ located at $(0,0,b)$. Additionally, the surface current
generates another magnetic field in the zone $z>0$, also proportional to the
electric field, equivalent to that generated by a magnetic monopole $p_{2}$
located, for sake of symmetry, at $(0,0,-b)$.

The electric and magnetic potentials in both regions, according to the image
method (Appendix 1) are:111The lack of symmetry of the magnetic field (axial
vector) with respect to a horizontal plane lead us to write $p_{1}$ and
$p_{2}$ as unknowns, with the expectation to find out $p_{1}=-p_{2}$.

$\displaystyle V(x,z>0)$
$\displaystyle=\frac{Q/\epsilon_{1}}{[x^{2}+(z-a)^{2}]^{1/2}}+\frac{q}{[x^{2}+(z+a)^{2}]^{1/2}}$
(15) $\displaystyle V(x,z<0)$
$\displaystyle=\frac{Q/\epsilon_{1}}{[x^{2}+(z-a)^{2}]^{1/2}}+\frac{q}{[x^{2}+(z-a)^{2}]^{1/2}}$
(16) $\displaystyle U(x,z>0)$
$\displaystyle=\frac{p_{1}}{[x^{2}+(z+b)^{2}]^{1/2}}$ (17) $\displaystyle
U(x,z<0)$ $\displaystyle=\frac{p_{2}}{[x^{2}+(z-b)^{2}]^{1/2}}$ (18)

The ${\bf B}_{\perp}=-\frac{\partial U}{\partial z}$ continuity across the
boundary means that $(\frac{\partial U(z>0)}{\partial z})_{0}=(\frac{\partial
U(z<0)}{\partial z})_{0}$. Then:

$\displaystyle
p_{1}\,\frac{(\cancelto{0}{z}+b)}{[x^{2}+(\cancelto{0}{z}+b)^{2}]^{3/2}}=p_{2}\,\frac{(\cancelto{0}{z}-b)}{[x^{2}+(\cancelto{0}{z}-b)^{2}]^{3/2}}\to\boxed{p_{1}=-p_{2}}$

Continuity at $z=0$ of ${\bf D}_{\perp}$, eq. (13), $\epsilon_{1}{\bf
E}_{1\perp}=\epsilon_{2}{\bf E}_{2\perp}+\alpha B_{z}$ with
$B_{z}=-(\frac{\partial U}{\partial z})_{0}$ and ${\bf
E}_{\perp}=-(\frac{\partial V}{\partial z})_{0}$, leads to:222In the
derivation we have calculated $B_{z}$ as the $U(z>0)$ derivative at $z=0$. The
same result is achieved with $U(z<0)$ since a double change of sign (monopole
charge and position) compensate.

$\displaystyle\epsilon_{1}\,\left(\frac{Q}{\epsilon_{1}}\frac{(\cancelto{0}{z}-a)}{[x^{2}+(\cancelto{0}{z}-a)^{2}]^{3/2}}+q\,\frac{(\cancelto{0}{z}+a)}{[x^{2}+(\cancelto{0}{z}+a)^{2}]^{3/2}}\right)$
$\displaystyle=$
$\displaystyle\epsilon_{2}\,\left(\frac{Q}{\epsilon_{1}}\frac{(\cancelto{0}{z}-a)}{[x^{2}+(\cancelto{0}{z}-a)^{2}]^{3/2}}+q\,\frac{(\cancelto{0}{z}-a)}{[x^{2}+(\cancelto{0}{z}-a)^{2}]^{3/2}}\right)$
$\displaystyle+\alpha\,p_{1}\,\frac{(\cancelto{0}{z}+b)}{[x^{2}+(\cancelto{0}{z}+b)^{2}]^{3/2}}$
$\displaystyle\to
Q\,(1-\frac{\epsilon_{2}}{\epsilon_{1}})\frac{(-a)}{[x^{2}+a^{2}]^{3/2}}+\frac{a\,q}{[x^{2}+a^{2}]^{3/2}}\,(\epsilon_{1}+\epsilon_{2})$
$\displaystyle=$
$\displaystyle\frac{b}{[x^{2}+b^{2}]^{3/2}}\,\alpha\,p_{1}\;\;\forall x$
$\displaystyle\to$ $\displaystyle\boxed{a=b}$ $\displaystyle\to
Q\,(\frac{\epsilon_{2}-\epsilon_{1}}{\epsilon_{1}})+(\epsilon_{1}+\epsilon_{2})\,q$
$\displaystyle=$
$\displaystyle\alpha\,p_{1}\to\boxed{p_{1}=\frac{Q}{\alpha}\frac{\epsilon_{2}-\epsilon_{1}}{\epsilon_{1}}+\frac{q}{\alpha}\;(\epsilon_{1}+\epsilon_{2})}$
(20)

The ${\bf H}_{\parallel}$ continuity at $z=0$, eq. (14),
$\frac{1}{\mu_{1}}{\bf B}_{1\parallel}=\frac{1}{\mu_{2}}{\bf
B}_{2\parallel}-\alpha{\bf E}_{\parallel}$, with $p_{1}=-p_{2}$, ${\bf
B}_{\parallel}=-(\frac{\partial U}{\partial x})_{z=0}$, ${\bf
E}_{\parallel}=-(\frac{\partial V}{\partial x})_{z=0}$, and $(\frac{\partial
V(z>0)}{\partial x})_{z=0}=(\frac{\partial V(z<0)}{\partial x})_{z=0}$, yields
(at $z=0$):

$\displaystyle\frac{1}{\mu_{1}}\frac{x\,p_{1}}{[x^{2}+a^{2}]^{3/2}}$
$\displaystyle=$
$\displaystyle\frac{1}{\mu_{2}}\frac{x\,(-p_{1})}{[x^{2}+a^{2}]^{3/2}}+\alpha\;\left(\frac{Q}{\epsilon_{1}}\frac{(-x)}{[x^{2}+a^{2}]^{3/2}}+q\;\frac{(-x)}{[x^{2}+a^{2}]^{3/2}}\right)$
$\displaystyle\to(\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}})p_{1}$ $\displaystyle=$
$\displaystyle-\alpha\;(\frac{Q}{\epsilon_{1}}+q)$ (21) $\displaystyle\to$
$\displaystyle\boxed{p_{1}=-\frac{\alpha}{\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}}}\;(\frac{Q}{\epsilon_{1}}+q)}$

From equations (20) and (21) it follows:

$\begin{array}[]{l}\frac{Q}{\alpha}\frac{\epsilon_{2}-\epsilon_{1}}{\epsilon_{1}}+\frac{q}{\alpha}(\epsilon_{1}+\epsilon_{2})=-\frac{\alpha}{\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}}}\;(\frac{Q}{\epsilon_{1}}+q)\\\
\\\
\to\frac{Q}{\epsilon_{1}}\left(\frac{\epsilon_{2}-\epsilon_{1}}{\alpha}+\frac{\alpha}{\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}}}\right)+q\;\left(\frac{\epsilon_{1}+\epsilon_{2}}{\alpha}+\frac{\alpha}{\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}}}\right)=0\\\
\\\ \to
q=-\frac{Q}{\epsilon_{1}}\frac{\frac{\epsilon_{2}-\epsilon_{1}}{\alpha}+\frac{\alpha}{\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}}}}{\frac{\epsilon_{2}+\epsilon_{1}}{\alpha}+\frac{\alpha}{\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}}}}\\\
\end{array}$
$\begin{array}[]{l}\to\boxed{q=\frac{Q}{\epsilon_{1}}\frac{(\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}})(\epsilon_{1}-\epsilon_{2})-\alpha^{2}}{(\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}})(\epsilon_{1}+\epsilon_{2})+\alpha^{2}}}\end{array}$
(22)

By replacing $q$, eq. (22), in eq. (21) and calling
$M=\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}}$, we find:

$\begin{array}[]{l}p_{1}=-\frac{\alpha}{M}\;(\frac{Q}{\epsilon_{1}}+\frac{Q}{\epsilon_{1}}\frac{M\,(\epsilon_{1}-\epsilon_{2})-\alpha^{2}}{M\,(\epsilon_{1}+\epsilon_{2})-\alpha^{2}})=-\frac{\alpha}{M}\;\frac{Q}{\epsilon_{1}}(1+\frac{(\epsilon_{1}-\epsilon_{2})-\alpha^{2}/M}{(\epsilon_{1}+\epsilon_{2})+\alpha^{2}/M})\\\
\\\
=-\frac{\alpha}{M}\;\frac{Q}{\epsilon_{1}}\frac{(\epsilon_{1}+\epsilon_{2})+\alpha^{2}/M+(\epsilon_{1}-\epsilon_{2})-\alpha^{2}/M}{(\epsilon_{1}+\epsilon_{2})+\alpha^{2}/M}=-\frac{\alpha
Q}{\epsilon_{1}}\frac{2\epsilon_{1}}{M(\epsilon_{1}+\epsilon_{2})+\alpha^{2}}\end{array}$
$\begin{array}[]{l}\to\boxed{p_{1}=-\alpha\frac{2Q}{(\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}})(\epsilon_{1}+\epsilon_{2})+\alpha^{2}}}\end{array}$
(23)

Eqs. (22) i (23) provide the value of the $Q$-induced electric charge and
magnetic monopole originated as images across the ordinary-topological
insulator interface.333By identifying $\alpha^{2}$ with $4P_{3}^{2}\alpha^{2}$
(for $P_{3}=\pm\frac{1}{2}$ i.e. $4P_{3}^{2}=1$, see [16] page 1185) we can
check that Eq. (22) matches eq. 3.2 and eq. (23) matches eq. 3b in [16].

###  3.3 Images of a magnetic charge close to the interface between ordinary
and topological insulator

Beyond the existence or not of a magnetic monopole, it is worth getting its
electric and magnetic images across an ordinary-topological interface. On the
one hand because two monopole at a given distance constitute a dipole, and
doubtless magnetic dipoles do exist, so that monopole image and the
superposition principle can be used to calculate dipole images. Additionally,
an electric charge between two interfaces yields a cascade of image
charges,[17] so it can be useful in this and other kind of more elaborate
calculations.

We have already pointed out that the electric field of a charge in the
vicinity of an interface with a topological insulator generates a surface
current, eq. (11), ${\bf J}_{\theta}=-\frac{c\alpha}{4\pi}\delta(z_{0}){\bf
n}\times{\bf E}$ (it is implied that the interface is located at $z_{0}$).
This current is the same as that generated by a magnetic monopole,
symmetrically located with respect to the source electric charge $Q$, having a
magnetic charge, eq. (23),
$p=-\frac{2\alpha}{(\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}})(\epsilon_{1}+\epsilon_{2})+\alpha^{2}}\,Q$.

Similarly, a magnetic monopole in the vicinity of a $z_{0}$-located interface
with a topological insulator creates a magnetic field $\bf B$ which in turn
generates a surface charge density, eq. (12),
$\rho_{\theta}=\frac{\alpha}{4\pi}\,\delta(z_{0})\,{\bf n}\cdot{\bf B}$.

As the monopole generates a radial field (like that of a charge or electric
monopole, see Figure), the charge density induced by the magnetic monopole
will generate an electric potential like that of an electric charge $q$
located at the position determined the image method.

Like images of an electric charge, a magnetic monopole $P$ generates an
electric ($V$) and magnetic ($U$) potentials given by:

$\displaystyle U_{1}(x,z>0)$
$\displaystyle=\frac{\mu_{1}\,P}{[x^{2}+(z-a)^{2}]^{1/2}}+\frac{p_{1}}{[x^{2}+(z+a)^{2}]^{1/2}}$
(24) $\displaystyle U_{2}(x,z<0)$
$\displaystyle=\frac{\mu_{1}\,P}{[x^{2}+(z-a)^{2}]^{1/2}}+\frac{p_{2}}{[x^{2}+(z-a)^{2}]^{1/2}}$
(25) $\displaystyle V_{1}(x,z>0)$
$\displaystyle=\frac{q}{[x^{2}+(z+b)^{2}]^{1/2}}$ (26) $\displaystyle
V_{2}(x,z<0)$ $\displaystyle=\frac{q}{[x^{2}+(z-b)^{2}]^{1/2}}$ (27)

with boundary conditions (BCs):

$\displaystyle\epsilon_{1}\left(\frac{\partial V_{1}}{\partial z}\right)_{0}$
$\displaystyle=\epsilon_{2}\left(\frac{\partial V_{2}}{\partial
z}\right)_{0}+\alpha\,\left(\frac{\partial U}{\partial z}\right)_{0}$ (29)
$\displaystyle\frac{1}{\mu_{1}}\left(\frac{\partial U_{1}}{\partial
x}\right)_{0}$ $\displaystyle=\frac{1}{\mu_{2}}\left(\frac{\partial
U_{2}}{\partial x}\right)_{0}-\alpha\,\left(\frac{\partial V}{\partial
x}\right)_{0}$ (30)

where since $\left(\frac{\partial U_{1}}{\partial
z}\right)_{0}=\left(\frac{\partial U_{2}}{\partial z}\right)_{0}$ and
$\left(\frac{\partial V_{1}}{\partial x}\right)_{0}=\left(\frac{\partial
V_{2}}{\partial x}\right)_{0}$, we can employ $U$ and $V$ of either region.

From the first BC, eq. (29), we have:

$\frac{\epsilon_{1}\,q\,b}{[x^{2}+b^{2}]^{3/2}}=\frac{-\epsilon_{2}\,q\,b}{[x^{2}+b^{2}]^{3/2}}+\alpha\,\left(\frac{\mu_{1}\,P\,a}{[x^{2}+a^{2}]^{3/2}}+\frac{\mu_{1}\,p_{1}\,a}{[x^{2}+a^{2}]^{3/2}}\right)\;\forall
x\;\to\boxed{a=b}$ (31)

and coming back to this equation with $a=b$ it follows:

$\epsilon_{1}\,q=-\epsilon_{2}\,q+\,\alpha\,\left(\mu_{1}\,P+p_{1}\right)\;\to\boxed{q=\frac{\alpha}{\epsilon_{1}+\epsilon_{2}}\,(p_{1}-\mu_{1}\,P)}$
(32)

From the second BC, eq. (30), we find out:

$\displaystyle\frac{1}{\mu_{1}}\left(\frac{\mu_{1}\,P\,x}{[x^{2}+a^{2}]^{3/2}}+\frac{p_{1}\,x}{[x^{2}+a^{2}]^{3/2}}\right)=\frac{1}{\mu_{2}}\left(\frac{\mu_{1}\,P\,x}{[x^{2}+a^{2}]^{3/2}}+\frac{p_{2}\,x}{[x^{2}+a^{2}]^{3/2}}\right)-\alpha\,\frac{q\,x}{[x^{2}+a^{2}]^{3/2}}\;\;\;\;\forall
x$ $\displaystyle\to
P+\frac{p_{1}}{\mu_{1}}=\frac{\mu_{1}}{\mu_{2}}\,P+\frac{1}{\mu_{2}}\,p_{2}-\alpha\,q\;\to\boxed{P\;\left(1-\frac{\mu_{1}}{\mu_{2}}\right)+\frac{p_{1}}{\mu_{1}}-\frac{p_{2}}{\mu_{2}}-\alpha\;q=0}$
(33)

The boundary condition ${\bf E}_{2\parallel}={\bf E}_{1\parallel}$ turns to be
an identity which brings nothing. Finally the boundary condition ${\bf
B}_{1\perp}={\bf B}_{2\perp}$ with ${\bf B}_{\perp}=-\frac{\partial
U}{\partial z}$ leads to:

$\displaystyle\frac{\mu_{1}\,P\,(\cancelto{0}{z}-a)}{[x^{2}+a^{2}]^{3/2}}+\frac{p_{1}\,(\cancelto{0}{z}+a)}{[x^{2}+a^{2}]^{3/2}}=\frac{\mu_{1}\,P\,(\cancelto{0}{z}-a)}{[x^{2}+a^{2}]^{3/2}}+\frac{p_{2}\,(\cancelto{0}{z}-a)}{[x^{2}+a^{2}]^{3/2}}$
$\displaystyle\to\boxed{p_{1}=-p_{2}}$ (34)

Now, by combining eqs. (3.3) and (34), taking into account (32), we find:

$\displaystyle
P\,\left(1-\frac{\mu_{1}}{\mu_{2}}\right)+p_{1}\left(\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}}\right)+\frac{\alpha^{2}}{\epsilon_{1}+\epsilon_{2}}\,\left(p_{1}-\mu_{1}\,P\right)=0$
$\displaystyle\to
P\,\left(1-\frac{\mu_{1}}{\mu_{2}}-\frac{\alpha^{2}\,\mu_{1}}{\epsilon_{1}+\epsilon_{2}}\right)+p_{1}\left(\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}}+\frac{\alpha^{2}}{\epsilon_{1}+\epsilon_{2}}\right)=0$
$\displaystyle\to
p_{1}=-P\;\frac{1-\frac{\mu_{1}}{\mu_{2}}-\mu_{1}\,\frac{\alpha^{2}}{\epsilon_{1}+\epsilon_{2}}}{\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}}+\frac{\alpha^{2}}{\epsilon_{1}+\epsilon_{2}}}\to\boxed{p_{1}=-\mu_{1}\,P\;\frac{\frac{1}{\mu_{1}}-\frac{1}{\mu_{2}}-\frac{\alpha^{2}}{\epsilon_{1}+\epsilon_{2}}}{\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}}+\frac{\alpha^{2}}{\epsilon_{1}+\epsilon_{2}}}}$
(35)

Note that if the insulator is not topological but ordinary, then we should
take $\alpha=0$ and obtain
$p_{1}=-\mu_{1}\,P\;\frac{\mu_{1}-\mu_{2}}{\mu_{1}+\mu_{2}}$ which is the
image of a monopole in front of an interface between two different ordinary
media.

By injecting $p_{1}$, eq. (3.3), into $q$ given by equation (32) we get:

$\displaystyle q=\frac{\alpha}{\epsilon_{1}+\epsilon_{2}}\,\mu_{1}\,P\
\left(\frac{\frac{1}{\mu_{2}}-\frac{1}{\mu_{1}}+\frac{\alpha^{2}}{\epsilon_{1}+\epsilon_{2}}}{\frac{1}{\mu_{2}}+\frac{1}{\mu_{1}}+\frac{\alpha^{2}}{\epsilon_{1}+\epsilon_{2}}}-1\right)$
(36)

We can simplify eqs. (3.3) i (36) as follow:

$\displaystyle\boxed{p_{1}=-\mu_{1}\,P\;\frac{(\frac{1}{\mu_{1}}-\frac{1}{\mu_{2}})(\epsilon_{1}+\epsilon_{2})-\alpha^{2}}{(\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}})(\epsilon_{1}+\epsilon_{2})+\alpha^{2}}}$
(37)
$\displaystyle\boxed{q=-\alpha\;\frac{2\,P}{(\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}})(\epsilon_{1}+\epsilon_{2})+\alpha^{2}}}$
(38)

The result is a perfect analogy of the image of an electric charge, with the
detail that the image $q$ of the charge $Q$, eq. (22), is proportional to
$(\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}})(\epsilon_{1}-\epsilon_{2})$ while the
image $p$ of the monopole $P$, eq. (38), is proportional to
$(\frac{1}{\mu_{1}}-\frac{1}{\mu_{2}})(\epsilon_{1}+\epsilon_{2})$. Dielectric
constants and magnetic permeabilities exchange their role.

Exercise
An electric charge is located in the central region of a quantum well build up
by a topological insulator surrounded by an ordinary insulator. (a) Calculate
the electrical potential generated in the central region, assuming that both
the insulator of the central region and that of the surrounding barriers, are
ordinary. Check that the result agrees with equations (2.8) and (2.14) by
Kumagai and Takagahara.[17] (b) Calculate the electrical and magnetic
potentials under the assumption that the insulator of the central region is
topological. Check that by forcing $\alpha=0$ the magnetic potential goes to
zero and the electric potential matches that obtained in the previous section.

## References

  * [1] B. Andrei Bernevig, Taylor L. Hughes, and Shou-Cheng Zhang, Science 314 (2006) 1757.
  * [2] F. Wilczek, Phys. Rev. Lett. 58 (1987) 1799.
  * [3] F. Wilczek, Phys. Rev. Lett. 40 (1978) 279.
  * [4] S. Weinberg, Phys. Rev. Lett. 40 (1978) 223.
  * [5] X-L Qi and S-C Zhang, Physics Today 63 (2010) 33.
  * [6] M. M. Vazifeh and M. Franz, Phys. Rev. Lett. 111 (2013) 027201.
  * [7] T. H. O’Dell, The electrodynamics of magneto-electric media, American Elsevier Pub. Co., New York 1970.
  * [8] See e.g. section 3, after eq. 62, in X-L Qi and S-C Zhang, Rev. Mod. Phys. 83 (2011) 1057.
  * [9] A. Martín-Ruiz, M. Cambiaso and L. F. Urrutia, Int. J. Mod. Phys. A 34 (2019) 1941002.
  * [10] The electromagnetic energy density $W=\frac{1}{2\mu}B^{2}+\frac{1}{2}\epsilon E^{2}$ can formally be identified with a momentum $B$, mass $\mu$, coordinate $E$ and force constant $\epsilon$ harmonic oscil$\cdot$lator, so that energy $W$ is a sum of kinetic $T$ and potential $V$ energies. By convention, Lagrangian is defined as ${\cal L}=T-V$, but can also be defined with opposite sign ${\cal L}=V-T=\frac{1}{2}\epsilon E^{2}-\frac{1}{2}\frac{1}{\mu}B^{2}$ for both definitions lead to the same Euler-Lagrange equation $\frac{d}{dt}\left(\frac{\partial{\cal L}}{\partial\dot{q}}\right)-\left(\frac{\partial{\cal L}}{\partial q}\right)=0$.
  * [11] C. Civelek and T. F. Bechteler, Int. J. of Eng. Sci. 46 (2008) 1218.
  * [12] Dirac delta is defined as the derivative of the Heaviside function $\\{H(z>0)=1,\;H(z<0)=0\\}$, $\delta(z)=\frac{d}{dz}H(z)$. Since $\theta=\pi\,H(z_{0})$, then $\frac{d\theta}{dz}=\pi\,\frac{dH}{dz}=\pi\,\delta(z_{0})$ and therefore ${\boldsymbol{\nabla}}\theta=\pi\delta(z_{0}){\bf n}$.
  * [13] A. Sekine and K. Nomura, J. Appl. Phys. 129 (2021) 141101.
  * [14] L. Visinelli, Mod. Phys. Lett. A 28 (2013) 135062.
  * [15] Since for image charges calculation only equations (7) i (8) are employed, Visinelli proposal[14] does not affect this kind of calculation.
  * [16] X-L Qi, R. Li, J. Zang and S-C Zhang, Science 329 (2009) 1184.
  * [17] M. Kumagai and T. Takagahara, Phys. Rev. B 40 (1989) 12359.

Part III:

##  4 Appendix 1: Boundary conditions at the interface of ordinary insulators

Maxwell’s equations in integral form (and MKS rational system) are:

$\begin{array}[]{lll}\iint{\bf D}\cdot d{\bf
S}=\iiint\rho\;dv&&\iint{\bf\nabla}\times{\bf E}\cdot d{\bf S}=\oint{\bf
E}\cdot d{\boldsymbol{\ell}}=-\iint\frac{\partial{\bf B}}{\partial t}\cdot
d{\bf S}\\\ \iint{\bf B}\cdot d{\bf S}=\iiint{\bf\nabla}\cdot{\bf
B}\;dv=0&&\iint{\bf\nabla}\times{\bf H}\;d{\bf S}=\oint{\bf H}\cdot
d{\boldsymbol{\ell}}=\iint{\bf J}\cdot d{\bf S}+\iint\frac{\partial{\bf
D}}{\partial t}\cdot d{\bf S}\\\ \end{array}$ (39)

Should we consider time independence, the last two equations including time
derivatives are simplified. Then, Maxwell equations become:

$\begin{array}[]{lll}\iint{\bf D}\cdot d{\bf S}=\iiint\rho\;dv&&\oint{\bf
E}\cdot d{\boldsymbol{\ell}}=0\\\ \iint{\bf B}\cdot d{\bf S}=0&&\oint{\bf
H}\cdot d{\boldsymbol{\ell}}=\iint{\bf J}\cdot d{\bf S}\\\ \end{array}$ (40)

On left in the Figure we show a slab of area $\Delta S$ and infinitesimal
height $dx$. We rewrite the integral $\iint{\bf B}\cdot d{\bf S}=0$ for this
slab as:

$0=\iint{\bf B}\cdot d{\bf S}=\iint{\bf B}\cdot{\bf n}\,dS=({\bf
B}_{2}\cdot{\bf n}-{\bf B}_{1}\cdot{\bf n})\,\Delta S=({\bf B}_{2}-{\bf
B}_{1})\cdot{\bf n}\,\Delta S\to\boxed{{\bf B}_{2\perp}={\bf B}_{1\perp}}$

In a similar way, $\iint{\bf D}\cdot d{\bf S}=({\bf D}_{2}-{\bf
D}_{1})\cdot{\bf n}\,\Delta S$. On the other hand
$\iiint\rho\;dv=\iiint\rho\,dx\,\Delta S=\sigma\Delta S$. Then,

$({\bf D}_{2}-{\bf D}_{1})\cdot{\bf n}=\sigma\to\boxed{{\bf D}_{2\perp}={\bf
D}_{1\perp}+\sigma}$

On the right hand side of the figure it is drawn a circuit in the shape of a
rectangle with a long side $\Delta{\boldsymbol{\ell}}$ and a narrow
infinitesimal side $dx$. The circulation of the electric field $\oint{\bf
E}\cdot d{\boldsymbol{\ell}}=0$ in this circuit is:

$0=\oint{\bf E}\cdot d{\boldsymbol{\ell}}=({\bf E}_{2}-{\bf
E}_{1})\cdot\Delta{\boldsymbol{\ell}}\to\boxed{{\bf E}_{2\parallel}={\bf
E}_{1\parallel}}$

In a similar way, $\oint{\bf H}\cdot d{\boldsymbol{\ell}}=({\bf H}_{2}-{\bf
H}_{1})\cdot\Delta\ell$. Also, $\iint{\bf J}\cdot d{\bf S}=\int{\bf
J}\cdot(dx\,\Delta\ell){\boldsymbol{\tau}}=(\int{\bf
J}dx)\cdot{\boldsymbol{\tau}}\Delta\ell={\bf J}_{s}\Delta\ell$, where $J_{s}$
represents surface current at the interface. Therefore:

$\boxed{{\bf H}_{2\parallel}={\bf H}_{1\parallel}+{\bf J}_{s}}$

###  4.1 Boundary conditions at the interface of ordinary insulators in a.u.

Maxwell equations in differential form and a.u. reads:

$\begin{array}[]{lll}{\boldsymbol{\nabla}}\cdot{\bf
D}=4\pi\rho\;{\rm\;\;\;(Gauss\,law)}&&{\boldsymbol{\nabla}}\times{\bf
E}=-\frac{1}{c}\frac{\partial{\bf B}}{\partial
t}{\rm\;\;\;\;\;\;\;\;\;(Faraday\,law)}\\\ \\\ {\boldsymbol{\nabla}}\cdot{\bf
B}=0\;\;\;\;\;\;\;\;{\rm(Gauss\,law\,for\,magnetism)}&&{\bf\nabla}\times{\bf
H}=\frac{1}{c}\frac{\partial{\bf D}}{\partial t}+\frac{4\pi}{c}{\bf
J}\;{\rm(Ampere\,law)}\\\ \end{array}$

Therefore, $4\pi$ must be added to the ${\bf D}$ boundary condition and
$\frac{4\pi}{c}$ to the ${\bf H}$ boundary condition:

$\boxed{{\bf D}_{2\perp}={\bf
D}_{1\perp}+4\pi\sigma}\;\;\;\;\;\;\;\;\boxed{{\bf H}_{2\parallel}={\bf
H}_{1\parallel}+\frac{4\pi}{c}\,{\bf J}_{s}}$

##  5 Appendix 2: Electrostatic Image Method

Let $Q$ be a charge in a dielectric constant $\epsilon_{1}$ medium. $Q$
polarizes this medium. Should the medium be infinite then the result is that
$Q$ plus the medium generate the same electric field as that of an effective
charge $Q/\epsilon_{1}$ in a vacuum. If there is an interface separating two
different polarizability media, i.e., with different dielectric constant (see
figure), on the surface of the upper medium appears the lower ends of the last
dipole layer induced in this medium which are not compensated with the first
dipole layer induced in the lower medium. The result is that the charge $Q$
located at a distance $a$ from the interface plus the two media generate the
same field as that generated in a vacuum by an effective charge
$Q/\epsilon_{1}$ located where the charge $Q$ is located plus a surface
density $\sigma$ located at the interface between the media. The charge
distribution $\sigma$ generates in the upper medium the same field as that
generated by a charge $q_{1}$ located symmetrically in the lower medium at a
distance $b$, while in the lower medium it generates the same field than that
generated by a charge $q_{2}$ located in the same vertical as $Q$ (for
symmetry reasons) at a distance $c$ from the interface.

The electric potential in the upper medium $V(z>0)$ is generated by
$Q/\epsilon_{1}$ located at $z_{Q}=a$ plus $q_{1}$ located at $z_{q_{1}}=-b$:

$V(z>0)=\frac{Q/\epsilon_{1}}{[x^{2}+(z-a)^{2}]^{1/2}}+\frac{q_{1}}{[x^{2}+(z+b)^{2}]^{1/2}}$
(41)

The electric potential in the lower medium $V(z<0)$ is generated by
$Q/\epsilon_{1}$ located at $z_{Q}=a$ and $q_{2}$ located a $z_{q_{2}}=c$:

$V(z<0)=\frac{Q/\epsilon_{1}}{[x^{2}+(z-a)^{2}]^{1/2}}+\frac{q_{2}}{[x^{2}+(z-c)^{2}]^{1/2}}$
(42)

To determine the values $b,c,q_{1},q_{2}$ we apply the BCs obtained in
Appendix 1. From BC ${\bf E}_{2\parallel}={\bf E}_{1\parallel}$, with ${\bf
E}_{\parallel}=-\frac{\partial V}{\partial x}$ we find:

${\bf
E}_{1\parallel}(z>0)=\frac{Q}{\epsilon_{1}}\frac{x}{[x^{2}+(z-a)^{2}]^{3/2}}+q_{1}\;\frac{x}{[x^{2}+(z+b)^{2}]^{3/2}}$
(43) ${\bf
E}_{2\parallel}(z<0)=\frac{Q}{\epsilon_{1}}\frac{x}{[x^{2}+(z-a)^{2}]^{3/2}}+q_{2}\;\frac{x}{[x^{2}+(z-c)^{2}]^{3/2}}$
(44)

Equating the two equations at $z=0$ we find that:

$\frac{q_{1}}{[x^{2}+b^{2}]^{3/2}}=\frac{q_{2}}{[x^{2}+c^{2}]^{3/2}}\to\frac{q_{1}}{q_{2}}=\frac{[x^{2}+b^{2}]^{3/2}}{[x^{2}+c^{2}]^{3/2}}$
(45)

As the ratio $\frac{q_{1}}{q_{2}}$ is constant, so must be the fraction
$\frac{[x^{2}+b^{2}]^{3/2}}{[x^{2}+c^{2}]^{3/2}}$. Then, $b=c$ and therefore
$q_{1}=q_{2}$.

From $b=c$, $q_{1}=q_{2}=q$, BC ${\bf D}_{2\perp}={\bf D}_{1\perp}+\sigma$,
taking into account that there is no free surface charge ($\sigma=0$), and
that ${\bf D}_{\perp}=\epsilon{\bf E}_{\perp}$, with ${\bf
E}_{\perp}=-\frac{\partial V}{\partial z}$, we have:

${\bf
E}_{1\perp}(z>0)=\frac{Q}{\epsilon_{1}}\frac{z-a}{[x^{2}+(z-a)^{2}]^{3/2}}+q\;\frac{z+b}{[x^{2}+(z+b)^{2}]^{3/2}}$
(46) ${\bf
E}_{2\perp}(z<0)=\frac{Q}{\epsilon_{1}}\frac{z-a}{[x^{2}+(z-a)^{2}]^{3/2}}+q\;\frac{z-b}{[x^{2}+(z-b)^{2}]^{3/2}}$
(47)

Equating $\epsilon_{1}{\bf E}_{1\perp}=\epsilon_{2}{\bf E}_{2\perp}$ at $z=0$
we find:

$-Q\;\frac{a}{[x^{2}+a^{2}]^{3/2}}+\epsilon_{1}\;q\;\frac{b}{[x^{2}+b^{2}]^{3/2}}=-\frac{\epsilon_{2}}{\epsilon_{1}}\;Q\;\frac{a}{[x^{2}+a^{2}]^{3/2}}-\epsilon_{2}\;q\;\frac{b}{[x^{2}+b^{2}]^{3/2}}$
(48)

Grouping factors:

$Q\;(1-\frac{\epsilon_{2}}{\epsilon_{1}})\frac{a}{[x^{2}+a^{2}]^{3/2}}=q(\epsilon_{1}+\epsilon_{2})\;\frac{b}{[x^{2}+b^{2}]^{3/2}}\to\frac{\frac{1}{\epsilon_{1}}\;Q\;(\epsilon_{1}-\epsilon_{2})}{q\;(\epsilon_{1}+\epsilon_{2})}\frac{a}{b}=\frac{[x^{2}+a^{2}]^{3/2}}{[x^{2}+b^{2}]^{3/2}}.$
(49)

As the first term of the last equation is a constant, so must be the second,
which implies $a=b$, a result that entails:

$q=\frac{1}{\epsilon_{1}}\;\frac{(\epsilon_{1}-\epsilon_{2})}{(\epsilon_{1}+\epsilon_{2})}\;Q$
(50)

In summary, the electric potential generated by a charge $Q$ in a medium of
dielectric constant $\epsilon_{1}$ located at a distance $a$ from the
interface with another medium of dielectric constant $\epsilon_{2}$ is the
same to that generated in a vacuum by an effective charge
$Q^{*}=Q/\epsilon_{1}$ located at $Q$ site plus that generated by a image
charge
$q=\frac{1}{\epsilon_{1}}\;\frac{(\epsilon_{1}-\epsilon_{2})}{(\epsilon_{1}+\epsilon_{2})}\;Q$
located symmetrically in the other medium at the same distance from the
interface. In addition, $Q$ also generates in the other medium a potential
like the one generated in a vacuum by an effective charge
$Q^{*}=Q/\epsilon_{1}$ plus an effective charge
$q=\frac{1}{\epsilon_{1}}\;\frac{(\epsilon_{1}-\epsilon_{2})}{(\epsilon_{1}+\epsilon_{2})}\;Q$
both located at the $Q$ site.

##  6 Appendix 3: On the electrodynamics atomic units

There is some ambiguity in defining the atomic units of magnetic field and the
own Maxwell equations. Within the Lorentz force convention, the Maxwell
equations read:

$\displaystyle{\boldsymbol{\nabla}}\cdot{\bf D}$ $\displaystyle=4\pi\rho$
$\displaystyle{\boldsymbol{\nabla}}\cdot{\bf B}$ $\displaystyle=0$ (51)
$\displaystyle{\boldsymbol{\nabla}}\times{\bf E}$
$\displaystyle=-\frac{\partial{\bf B}}{\partial t}$
$\displaystyle{\bf\nabla}\times{\bf H}$ $\displaystyle=\frac{\partial{\bf
D}}{\partial t}+4\pi\,{\bf J}$ (52)

with the constitutive relations ${\bf D}={\bf E}+4\pi{\bf P}$, ${\bf
H}=c^{2}\,{\bf B}-4\pi{\bf M}$. In the last equation the value $c$ of the
speed of light in a.u. is numerically equal to $\frac{1}{\alpha}$, with
${\alpha}=1/137.036$ being the dimensionless fine-structure constant. These
equations are complemented by the Lorentz force law: ${\bf F}=e\,({\bf E}+{\bf
v}\times{\bf B})$.

Within the Gaussian convention, the Maxwell equations are:

$\displaystyle{\boldsymbol{\nabla}}\cdot{\bf D}$ $\displaystyle=4\pi\rho$
$\displaystyle{\boldsymbol{\nabla}}\cdot{\bf B}$ $\displaystyle=0$ (53)
$\displaystyle{\boldsymbol{\nabla}}\times{\bf E}$
$\displaystyle=-\alpha\frac{\partial{\bf B}}{\partial t}$
$\displaystyle{\bf\nabla}\times{\bf H}$
$\displaystyle=\alpha\frac{\partial{\bf D}}{\partial t}+4\pi\alpha{\bf J}$
(54)

with the constitutive relations ${\bf D}={\bf E}+4\pi{\bf P}$, ${\bf H}={\bf
B}-4\pi{\bf M}$, which look like those of the cgs-Gaussian system. In
addition, we have the Lorentz force law: ${\bf F}=e\,({\bf E}+\frac{{\bf
v}}{c}\times{\bf B})$.

Please note that the formal replacement $\alpha{\bf B}\to{\bf B}$, $\frac{{\bf
H}}{\alpha}\to{\bf H}$, $\frac{{\bf M}}{\alpha}\to{\bf M}$ in the Gaussian
convention equations retrieves those of the Lorentz convention.

Care should we have in using either convention. For example, the atomic unit
for magnetic field in the Lorentz force convention is $1\,a.u._{LF}\approx
2.35\cdot 10^{5}\,T$ while $1\,a.u._{G}\approx 1.72\cdot 10^{7}\,G$. Since
$1\,T=10^{4}\,G$ we can check that
$\frac{1\,a.u._{LF}}{1\,a.u._{G}}=\frac{1}{\alpha}$.

It would be instructive to have a look to the units guide provided by Andrea
Dal Corso at the address:
https:$\backslash\backslash$people.sissa.it$\backslash\sim$dalcorso$\backslash$notes$\backslash$units.pdf
The text is motivated by the conversion factors implemented in the QUANTUM
ESPRESSO (https:$\backslash\backslash$www.quantum-espresso.org/) but is of
general interest.

###  6.1 A bit more on electrodynamics units

In order to dig on the different atomic units we start from two fundamental
laws, the Coulomb force between electrical charges,
$F=k_{c}\frac{q_{1}\,q_{2}}{r^{2}}$ and the Biot-Savart force between
infinitesimal $dl_{1}$, $dl_{2}$ fragments of wires supporting current of
intensity $i_{1}$ and $i_{2}$, respectively,
$F=k_{b}\frac{i_{1}\,i_{2}}{r^{2}}dl_{1}\,dl_{2}$. The constants $k_{c}$ and
$k_{b}$ are related by the equation,

$k_{c}/k_{b}=c^{2},$ (55)

with $c$ the speed of light.

We can assume this relationship as an experimental fact. All the same, it can
be derived since, according to theory of relativity, magnetism is just
electric interaction between moving charges (see e.g. chap. 5 in E.M. Purcell
and J. Morin, Electricity and Magnetism, Cambridge University Press 2013).444
We can provide a simple, rather oversimplified, approach that can help to have
a taste of this. To this end we consider $S_{0}$, a lab frame of reference
attached to the positive cores of a wire supporting a current density $i_{1}$,
carried by negative electrons. In this inertial frame we observe a neutral
wire, i.e., the linear charge densities of cores $\rho_{+}$ and electrons
$\rho_{-}$ are equal: $\rho_{+}=\rho_{-}$. We employ the notation $L_{0}^{+}$
for the observed proper length between consecutive cores and $L^{-}$ for the
observed length between consecutive electrons (remember: the proper length
$L_{0}$ is the largest, so that the length $L$ observed from another inertial
frame $S$, moving at the speed $v$ –where the system is not seen at rest–, is
related to $L_{0}$ by $L=L_{0}\,\sqrt{1-v^{2}/c^{2}}$). We can write
$\rho_{+}=\frac{e}{L_{0}^{+}}$ and $\rho_{-}=-\frac{e}{L^{-}}$, with $e$
representing the absolute value of the electron charge. From
$\rho_{+}=-\rho_{-}$ we conclude,  $L_{0}^{+}=L^{-}$ (56) Next, we consider a
second inertial frame $S^{\prime}$ attached to electrons. Then, relative to
$S_{0}$, it has a speed $v$. In this new frame the distance between
consecutive moving cores is $L^{+}=L_{0}\sqrt{1-v^{2}/c^{2}}$ and that of
electrons $L_{0}^{-}$. Interestingly, from this inertial frame we see a net
density charge in the wire:
$\rho^{\prime}=\frac{e}{L^{+}}-\frac{e}{L_{0}^{-}}=\frac{e}{L_{0}^{+}}\frac{1}{\sqrt{1-v^{2}/c^{2}}}-\frac{e}{L^{-}}\,\sqrt{1-v^{2}/c^{2}}$
(57) With eq. (56) we have:
$\rho^{\prime}=\frac{e}{L_{0}^{+}}\frac{1}{\sqrt{1-v^{2}/c^{2}}}\,(1-1+\frac{v^{2}}{c^{2}})=\frac{e}{L^{+}}\frac{v^{2}}{c^{2}}$
(58) Since the density of cores observed in $S^{\prime}$ is
$\rho=\frac{e}{L^{+}}$, we can write:
$\rho^{\prime}=\rho\frac{v^{2}}{c^{2}}.$ (59) The cores are observed to move
in $S^{\prime}$ at speed $v$ then, the current intensity observed is
$i=\frac{dq}{dt}=\rho\,\frac{dl}{dt}=\rho\,v$, so that
$\rho^{\prime}=i\,\frac{v}{c^{2}}$.
A wire-length $dl$ with a non-zero density charge $\rho^{\prime}$ yields an
electrical field,
$E=\frac{k_{c}}{r^{2}}\rho^{\prime}dl=\frac{k_{c}}{r^{2}}\,i\,dl\,\frac{v}{c^{2}}.$
(60) Let’s write $i_{1}$ and $dl_{1}$ instead of $i$ and $dl$ to indicate that
they correspond to a given wire that we refer to as wire 1.
Consider next that $S^{\prime}$ (attached to electrons of wire 1) sees a net
charge $dq_{2}=\rho_{2}dl_{2}$ moving along wire 2 at speed $v$. Then
$S^{\prime}$ sees a current intensity $i_{2}=\rho_{2}v$ in wire 2 and exert on
it a force:
$F=E\,dq_{2}=\frac{k_{c}}{r^{2}}\,i_{1}\,dl_{1}\,\frac{v}{c^{2}}\,\rho_{2}\,dl_{2}=\frac{k_{c}}{c^{2}}\,\frac{i_{1}\,i_{2}}{r^{2}}\,dl_{1}\,dl_{2},$
(61) that, if compared to the Biot-Savart force,
$F=k_{b}\frac{i_{1}\,i_{2}}{r^{2}}dl_{1}\,dl_{2}$, lead us to say:
$\frac{k_{c}}{c^{2}}=k_{b}$.  Once we have settled eq. (55), there are several
options to define these constants, yielding the different unit systems. The
cgs Gaussian unit system assumes $k_{c}=1$ and, according to (55),
$k_{b}=1/c^{2}$. The same does the atomic unit system in the so-called
Gaussian convention. However, while $k_{c}$ has the same numerical value in
both unit systems, taking into account that $c=2.89\cdot 10^{10}\;cm/s$ and
$c=137.036\;a.u.$, $k_{b}$ has a different numerical value in either system.

On the other hand, the international SI or MKS unit system assumes
$k_{c}=\frac{1}{4\pi\epsilon_{0}}$ and $k_{b}=\frac{\mu_{0}}{4\pi}$ with
$\epsilon_{0}=8.854\cdot 10^{-12}\;C^{2}/N\cdot m^{2}$, $\mu_{0}=4\pi\cdot
10^{-7}\;N/A^{2}$. Again,
$\frac{k_{c}}{k_{b}}=\frac{1}{\epsilon_{0}\mu_{0}}=8.988\cdot
10^{16}\;m^{2}/s^{2}=c^{2}$.

Another challenge is the transformation of electromagnetic equations between
both unit systems. To this end we first invoke the Coulomb equation:

$\frac{q^{G}_{1}\,q^{G}_{2}}{r^{2}}\hskip 28.45274ptvs.\hskip
28.45274pt\frac{1}{4\pi\epsilon_{0}}\frac{q^{SI}_{1}\,q^{SI}_{2}}{r^{2}}=\frac{\frac{q^{SI}_{1}}{\sqrt{4\pi\epsilon_{0}}}\frac{q^{SI}_{2}}{\sqrt{4\pi\epsilon_{0}}}}{r^{2}}$

where $r$ must be expressed either, in $cm$ (Gaussian formula) or in $m$ (SI
formula). Then, we see that

$q^{SI}=\sqrt{4\pi\epsilon_{0}}\,q^{G}.$ (62)

We consider next the electric field:

$E^{SI}=\frac{1}{4\pi\epsilon_{0}}\,\frac{q^{SI}}{r^{2}}=\frac{1}{\sqrt{4\pi\epsilon_{0}}}\,\frac{q^{G}}{r^{2}}=\frac{1}{\sqrt{4\pi\epsilon_{0}}}\,E^{G}$
(63)

In a similar way, we derive the potential:

$V^{SI}=\frac{V^{G}}{\sqrt{4\pi\epsilon_{0}}}$ (64)

In dielectric media, the Coulomb law turns into
$F=\frac{1}{4\pi\epsilon_{0}}\,\frac{1}{\epsilon^{\prime}}\frac{q^{SI}_{1}\,q^{SI}_{2}}{r^{2}}=\frac{1}{\epsilon^{\prime}}\frac{q^{G}_{1}\,q^{G}_{2}}{r^{2}}$.
Then,

$\epsilon^{SI}=\epsilon_{0}\epsilon^{G}=\epsilon_{0}\epsilon^{\prime}\hskip
28.45274pti.e.\hskip 28.45274pt\epsilon^{G}=\epsilon^{\prime}$ (65)

Finally,

$D^{SI}=\epsilon^{SI}E^{SI}=\epsilon_{0}\epsilon^{G}\,\frac{1}{\sqrt{4\pi\epsilon_{0}}}\,E^{G}=\sqrt{\frac{\epsilon_{0}}{4\pi}}\;D^{G}.$
(66)

In the magnetic equations we have $\frac{\mu_{0}}{4\pi}$ instead of
$\frac{1}{4\pi\epsilon_{0}}$ and poles instead of charges. Then, proceeding in
a similar way, we find out:

$\displaystyle p^{SI}=\sqrt{\frac{4\pi}{\mu_{0}}}\,p^{G}$ (67) $\displaystyle
B^{SI}=\sqrt{\frac{\mu_{0}}{4\pi}}\,B^{G}$ (68) $\displaystyle
U^{SI}=\sqrt{\frac{\mu_{0}}{4\pi}}\,U^{G}$ (69)
$\displaystyle\mu^{SI}=\mu_{0}\mu^{G}$ (70) $\displaystyle
H^{SI}=\frac{1}{\sqrt{4\pi\mu_{0}}}\,H^{G}$ (71)

A complete table with Gaussian-SI equivalences and constant numerical values
is reported by wikipedia
(https:$\backslash\backslash$en.wikipedia.org$\backslash$wiki$\backslash{\rm
Gaussian\\_units}$).

With these equation we can translate eqs. in section 3.2 from Gaussian (or
atomic unit in Gaussian convention) to SI units. I enclose, next, some
examples.

  1. 1.

Equation (17):
$U^{G}(x,z>0)=\frac{p^{G}_{1}}{[x^{2}+(z+b)^{2}]^{1/2}}$ $\to$
$U^{SI}=\sqrt{\frac{\mu_{0}}{4\pi}}U^{G}=\sqrt{\frac{\mu_{0}}{4\pi}}\frac{\sqrt{\frac{\mu_{0}}{4\pi}}\,p^{SI}_{1}}{[x^{2}+(z+b)^{2}]^{1/2}}=\frac{\mu_{0}}{4\pi}\frac{p^{SI}_{1}}{[x^{2}+(z+b)^{2}]^{1/2}}$

  2. 2.

Boundary condition ${\bf D}^{G}_{1\perp}={\bf D}^{G}_{2\perp}+\alpha
B^{G}_{z}$ yields ${\bf D}^{SI}_{1\perp}={\bf
D}^{SI}_{2\perp}+\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\alpha B^{SI}_{z}$
i.e. $\epsilon^{SI}_{1}{\bf E}^{SI}_{1\perp}=\epsilon^{SI}_{2}{\bf
E}^{SI}_{2\perp}+\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\alpha B^{SI}_{z}$.

  3. 3.

Boundary condition ${\bf H}^{G}_{1\parallel}={\bf
H}^{G}_{2\parallel}-\alpha{\bf E}^{G}_{2\parallel}$ yields ${\bf
H}^{SI}_{1\parallel}={\bf
H}^{SI}_{2\parallel}-\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\alpha{\bf
E}^{SI}_{2\parallel}$
i.e. $\frac{1}{\mu^{SI}_{1}}{\bf
B}^{SI}_{1\parallel}=\frac{1}{\mu^{SI}_{2}}{\bf
B}^{SI}_{2\parallel}-\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\alpha{\bf
E}^{SI}_{2\parallel}$

  4. 4.

Equation (21):
$p^{G}_{1}=-\frac{\alpha}{\frac{1}{\mu^{G}_{1}}+\frac{1}{\mu^{G}_{2}}}\;(\frac{Q^{G}}{\epsilon^{G}_{1}}+q^{G})$
with $p^{SI}_{1}=\sqrt{\frac{4\pi}{\mu_{0}}}\,p^{G}_{1}$ yields,
$p^{SI}_{1}=-\sqrt{\frac{4\pi}{\mu_{0}}}\frac{\alpha}{\frac{\mu_{0}}{\mu^{SI}_{1}}+\frac{\mu_{0}}{\mu^{SI}_{2}}}\;\left(\frac{1}{\sqrt{4\pi\epsilon_{0}}}\frac{Q^{SI}}{\epsilon^{\prime}_{1}}-\frac{1}{\sqrt{4\pi\epsilon_{0}}}q^{SI}\right)=-\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\frac{\alpha/\mu_{0}}{\frac{1}{\mu^{SI}_{1}}+\frac{1}{\mu^{SI}_{2}}}\left(\frac{Q^{SI}}{\epsilon^{SI}_{1}}-\frac{q^{SI}}{\epsilon_{0}}\right)$

  5. 5.

Equation (22):
$q^{G}=\frac{1}{\sqrt{4\pi\epsilon_{0}}}q^{SI}=\frac{Q^{G}}{\epsilon^{G}_{1}}\frac{(\frac{1}{\mu^{G}_{1}}+\frac{1}{\mu^{G}_{2}})(\epsilon^{G}_{1}-\epsilon^{G}_{2})-\alpha^{2}}{(\frac{1}{\mu^{G}_{1}}+\frac{1}{\mu^{G}_{2}})(\epsilon^{G}_{1}+\epsilon^{G}_{2})+\alpha^{2}}$
yields,
$q^{SI}=\epsilon_{0}\;\frac{Q^{SI}}{\epsilon^{SI}_{1}}\frac{(\frac{1}{\mu^{SI}_{1}}+\frac{1}{\mu^{SI}_{2}})(\epsilon^{SI}_{1}-\epsilon^{SI}_{2})-\frac{\epsilon_{0}}{\mu_{0}}\,\alpha^{2}}{(\frac{1}{\mu^{SI}_{1}}+\frac{1}{\mu^{SI}_{2}})(\epsilon^{SI}_{1}+\epsilon^{SI}_{2})+\frac{\epsilon_{0}}{\mu_{0}}\,\alpha^{2}}$

  6. 6.

Finally, Equation (23):
$p^{G}_{1}=-\alpha\,\frac{2Q^{G}}{(\frac{1}{\mu^{G}_{1}}+\frac{1}{\mu^{G}_{2}})(\epsilon^{G}_{1}+\epsilon^{G}_{2})+\alpha^{2}}$
yields,
$p^{SI}_{1}=-\frac{1}{\mu_{0}}\left(\sqrt{\frac{\epsilon_{0}}{\mu_{0}}}\,\alpha\right)\,\frac{2Q^{SI}}{(\frac{1}{\mu^{SI}_{1}}+\frac{1}{\mu^{SI}_{2}})(\epsilon^{SI}_{1}+\epsilon^{SI}_{2})+\frac{\epsilon_{0}}{\mu_{0}}\,\alpha^{2}}$

As for the dimensionless fine-structure constant $\alpha$, it has the same
numerical value in all unit systems ($\alpha=1/137.036$), but presents
different formulas in different systems. Thus,
$\alpha^{SI}=\frac{e^{2}}{4\pi\epsilon_{0}\hbar c}=7.297\cdot 10^{-3}$. Since
$\epsilon^{G}_{0}=\frac{1}{4\pi}$, then $\alpha^{G}=\frac{e^{2}}{\hbar
c}=7.297\cdot 10^{-3}$. Finally, from $e=\hbar=1\,a.u.$,
$\alpha^{a.u.}=\frac{1}{c}=7.297\cdot 10^{-3}$.


* [143] S. James, P. Wohlhart, M. Kalakrishnan, D. Kalashnikov, A. Irpan, J. Ibarz, S. Levine, R. Hadsell, and K. Bousmalis, “Sim-to-Real via Sim-to-Sim: Data-efficient Robotic Grasping via Randomized-to-Canonical Adaptation Networks,” in _CVPR_ , 2019.
  * [144] Y. Chebotar, A. Handa, V. Makoviychuk, M. Macklin, J. Issac, N. Ratliff, and D. Fox, “Closing the Sim-to-Real Loop: Adapting Simulation Randomization with Real World Experience,” in _ICRA_ , 2018.
  * [145] J. Tan, T. Zhang, E. Coumans, A. Iscen, Y. Bai, D. Hafner, S. Bohez, and V. Vanhoucke, “Sim-to-Real: Learning Agile Locomotion For Quadruped Robots,” in _RSS_ , 2018.
  * [146] X. B. Peng, M. Andrychowicz, W. Zaremba, and P. Abbeel, “Sim-to-real transfer of robotic control with dynamics randomization,” in _ICRA_ , 2018.
  * [147] M. Andrychowicz, F. Wolski, A. Ray, J. Schneider, R. Fong, P. Welinder, B. McGrew, J. Tobin, O. P. Abbeel, and W. Zaremba, “Hindsight experience replay,” in _NIPS_ , 2017.
  * [148] M. Feurer, J. T. Springenberg, and F. Hutter, “Initializing bayesian hyperparameter optimization via meta-learning,” in _AAAI_ , 2015.
  * [149] C. Finn, P. Abbeel, and S. Levine, “Model-agnostic meta-learning for fast adaptation of deep networks,” in _ICML_ , 2017.
  * [150] I. Clavera, A. Nagabandi, R. S. Fearing, P. Abbeel, S. Levine, and C. Finn, “Learning to Adapt in Dynamic, Real-World Environments through Meta-Reinforcement Learning,” in _ICLR_ , 2019.
  * [151] S. Sæmundsson, K. Hofmann, and M. P. Deisenroth, “Meta Reinforcement Learning with Latent Variable Gaussian Processes,” in _UAI_ , 2018.
  * [152] J. Hwangbo, J. Lee, A. Dosovitskiy, D. Bellicoso, V. Tsounis, V. Koltun, and M. Hutter, “Learning agile and dynamic motor skills for legged robots,” _Science Robotics_ , vol. 4, no. 26, 2019.
  * [153] N. G. Tsagarakis, G. Metta, _et al._ , “icub: the design and realization of an open humanoid platform for cognitive and neuroscience research,” _Advanced Robotics_ , vol. 21, no. 10, pp. 1151–1175, 2007.
  * [154] P. Maiolino, M. Maggiali, G. Cannata, G. Metta, and L. Natale, “A flexible and robust large scale capacitive tactile system for robots,” _IEEE Sensors Journal_ , vol. 13, no. 10, pp. 3910–3917, 2013.
  * [155] T. Dean and K. Kanazawa, “A model for reasoning about persistence and causation,” _Comput. Intell._ , vol. 5, pp. 142–150, 1989.
  * [156] C. Boutilier, R. Dearden, and M. Goldszmidt, “Stochastic dynamic programming with factored representations,” _Artif. Intel._ , vol. 121, pp. 49–107, 2000\.
  * [157] A. J. Ijspeert, “Central pattern generators for locomotion control in animals and robots: a review,” _Neural Netw._ , vol. 21, pp. 642–653, 2008.
  * [158] V. C. Kumar, S. Ha, and K. Yamane, “Improving Model-Based Balance Controllers using Reinforcement Learning and Adaptive Sampling,” in _ICRA_ , 2018.
  * [159] R. Jonschkowski and O. Brock, “Learning state representations with robotic priors,” _Autonomous Robots_ , vol. 39, no. 3, pp. 407–428, 2015.
  * [160] T. Lesort, N. Díaz-Rodríguez, J.-F. Goudou, and D. Filliat, “State representation learning for control: An overview,” _Neural Netw._ , vol. 108, pp. 379–392, 2018.
  * [161] J. Oh, X. Guo, H. Lee, R. L. Lewis, and S. Singh, “Action-conditional video prediction using deep networks in atari games,” in _NIPS_ , 2015.
  * [162] D. Ha and J. Schmidhuber, “World models,” _arXiv preprint arXiv:1803.10122_ , 2018.
  * [163] J.-A. M. Assael, N. Wahlström, T. B. Schön, and M. P. Deisenroth, “Data-efficient learning of feedback policies from image pixels using deep dynamical models,” _NIPS Deep RL Workshop_ , 2015.
  * [164] S. A. Eslami, D. J. Rezende, F. Besse, F. Viola, A. S. Morcos, M. Garnelo, A. Ruderman, A. A. Rusu, I. Danihelka, K. Gregor, _et al._ , “Neural scene representation and rendering,” _Science_ , vol. 360, no. 6394, pp. 1204–1210, 2018.
  * [165] J. A. Musick and C. J. Limpus, “Habitat utilization and migration in juvenile sea turtles,” _The biology of sea turtles_ , vol. 1, pp. 137–163, 1997.
  * [166] T. Lesort, M. Seurin, X. Li, N. D. Rodríguez, and D. Filliat, “Unsupervised state representation learning with robotic priors: a robustness benchmark,” _arXiv preprint arXiv:1709.05185_ , 2017.
  * [167] T.-H. Pham, G. De Magistris, and R. Tachibana, “OptLayer-Practical Constrained Optimization for Deep Reinforcement Learning in the Real World,” in _ICRA_ , 2018.
  * [168] T. Haarnoja, V. Pong, A. Zhou, M. Dalal, P. Abbeel, and S. Levine, “Composable Deep Reinforcement Learning for Robotic Manipulation,” in _ICRA_ , 2018.
  * [169] T. Pinville, S. Koos, J.-B. Mouret, and S. Doncieux, “How to promote generalisation in evolutionary robotics: the progab approach,” in _GECCO_ , 2011.
  * [170] B. C. Da Silva, G. Konidaris, and A. G. Barto, “Learning parameterized skills,” in _ICML_ , 2012.
  * [171] J. Kober, A. Wilhelm, E. Oztop, and J. Peters, “Reinforcement learning to adjust parametrized motor primitives to new situations,” _Autonomous Robots_ , vol. 33, no. 4, pp. 361–379, 2012.
  * [172] A. Fabisch and J. H. Metzen, “Active contextual policy search,” _JMLR_ , vol. 15, no. 1, pp. 3371–3399, 2014.
  * [173] T. Schaul, D. Horgan, K. Gregor, and D. Silver, “Universal value function approximators,” in _ICML_ , 2015.
  * [174] P. Karkus, A. Kupcsik, D. Hsu, and W. S. Lee, “Factored Contextual Policy Search with Bayesian Optimization,” in _BayesOpt’16: Proceedings of the International Workshop “Bayesian Optimization: Black-box Optimization and Beyond” at NIPS_ , 2016.
  * [175] S. Ha and C. K. Liu, “Evolutionary optimization for parameterized whole-body dynamic motor skills,” in _ICRA_ , 2016.
  * [176] Y. Zhu, R. Mottaghi, E. Kolve, J. J. Lim, A. Gupta, L. Fei-Fei, and A. Farhadi, “Target-driven visual navigation in indoor scenes using deep reinforcement learning,” in _ICRA_ , 2017.
  * [177] P. Rauber, A. Ummadisingu, F. Mutz, and J. Schmidhuber, “Hindsight policy gradients,” _arXiv preprint arXiv:1711.06006_ , 2017.
  * [178] D. Ghosh, A. Singh, A. Rajeswaran, V. Kumar, and S. Levine, “Divide-and-conquer reinforcement learning,” in _ICLR_ , 2018.
  * [179] D. J. Mankowitz, A. Žídek, A. Barreto, D. Horgan, M. Hessel, J. Quan, J. Oh, H. van Hasselt, D. Silver, and T. Schaul, “Unicorn: Continual learning with a universal, off-policy agent,” _arXiv preprint arXiv:1802.08294_ , 2018.
  * [180] M. P. Deisenroth, P. Englert, J. Peters, and D. Fox, “Multi-task policy search for robotics,” in _ICRA_ , 2014.
  * [181] S. Paul, K. Chatzilygeroudis, K. Ciosek, J.-B. Mouret, M. A. Osborne, and S. Whiteson, “Alternating Optimisation and Quadrature for Robust Control,” in _AAAI_ , 2018.
  * [182] V. Vassiliades and C. Christodoulou, “Toward nonlinear local reinforcement learning rules through neuroevolution,” _Neural Computation_ , vol. 25, no. 11, pp. 3020–3043, 2013.
  * [183] J. X. Wang, Z. Kurth-Nelson, D. Tirumala, H. Soyer, J. Z. Leibo, R. Munos, C. Blundell, D. Kumaran, and M. Botvinick, “Learning to reinforcement learn,” _arXiv preprint arXiv:1611.05763_ , 2016.
  * [184] J. Harrison, A. Sharma, R. Calandra, and M. Pavone, “Control Adaptation via Meta-Learning Dynamics,” in _Workshop on Meta-Learning at NeurIPS 2018_ , 2018.
  * [185] W. Yu, C. K. Liu, and G. Turk, “Preparing for the unknown: Learning a universal policy with online system identification,” in _RSS_ , 2017.
  * [186] A. Rajeswaran, S. Ghotra, B. Ravindran, and S. Levine, “Epopt: Learning robust neural network policies using model ensembles,” _arXiv preprint arXiv:1610.01283_ , 2016.
  * [187] R. E. Kalman, “A new approach to linear filtering and prediction problems,” _J. Basic. Eng._ , vol. 82, pp. 35–45, 1960.
  * [188] D. Mayne, “A second-order gradient method for determining optimal trajectories of non-linear discrete-time systems,” _International Journal of Control_ , vol. 3, no. 1, pp. 85–95, 1966.
  * [189] D. H. Jacobson and D. Q. Mayne, _Differential dynamic programming_. Elsevier, 1970.
  * [190] E. Todorov and W. Li, “A generalized iterative lqg method for locally-optimal feedback control of constrained nonlinear stochastic systems,” in _American Control Conference_ , 2005.
  * [191] J. Koenemann, A. Del Prete, Y. Tassa, E. Todorov, O. Stasse, M. Bennewitz, and N. Mansard, “Whole-body model-predictive control applied to the HRP-2 humanoid,” in _IROS_ , 2015.
  * [192] K. L. Moore, M. Dahleh, and S. Bhattacharyya, “Iterative learning control: A survey and new results,” _Journal of Field Robotics_ , vol. 9, no. 5, pp. 563–594, 1992.
  * [193] D. A. Bristow, M. Tharayil, and A. G. Alleyne, “A survey of iterative learning control,” _IEEE Control Systems_ , vol. 26, pp. 96–114, 2006.
  * [194] K. S. Lee, I.-S. Chin, H. J. Lee, and J. H. Lee, “Model predictive control technique combined with iterative learning for batch processes,” _AIChE Journal_ , vol. 45, no. 10, pp. 2175–2187, 1999.
  * [195] J. H. Lee, K. S. Lee, and W. C. Kim, “Model-based iterative learning control with a quadratic criterion for time-varying linear systems,” _Automatica_ , vol. 36, no. 5, pp. 641–657, 2000.
  * [196] Y. Wang, D. Zhou, and F. Gao, “Iterative learning model predictive control for multi-phase batch processes,” _Journal of Process Control_ , vol. 18, no. 6, pp. 543–557, 2008.
  * [197] T. Zhang, G. Kahn, S. Levine, and P. Abbeel, “Learning deep control policies for autonomous aerial vehicles with mpc-guided policy search,” in _ICRA_ , 2016.
  * [198] S. Karaman and E. Frazzoli, “Sampling-based algorithms for optimal motion planning,” _IJRR_ , vol. 30, no. 7, pp. 846–894, 2011.
  * [199] C. B. Browne, E. Powley, D. Whitehouse, S. M. Lucas, P. I. Cowling, P. Rohlfshagen, S. Tavener, D. Perez, S. Samothrakis, and S. Colton, “A survey of monte carlo tree search methods,” _IEEE Transactions on Computational Intelligence and AI in Games_ , vol. 4, pp. 1–43, 2012.
  * [200] M. Duarte, J. Gomes, S. M. Oliveira, and A. L. Christensen, “Evolution of repertoire-based control for robots with complex locomotor systems,” _IEEE Trans. Evol. Comput._ , 2017.
  * [201] D. Clever, M. Harant, K. Mombaur, M. Naveau, O. Stasse, and D. Endres, “Cocomopl: A novel approach for humanoid walking generation combining optimal control, movement primitives and learning and its transfer to the real robot hrp-2,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 2, pp. 977–984, 2017.

|  Konstantinos Chatzilygeroudis is currently a post-doctoral fellow at the
LASA team at EPFL. He obtained a B.Sc. and M.Sc. in Computer Science and
Engineering from the University of Patras in 2014, and a Ph.D. in Robotics and
Machine Learning from Inria Nancy-Grand Est (France) and the University of
Lorraine. His research interests lie in the area of artificial intelligence
and focus on reinforcement learning and fast robot adaptation.
Website: http://costashatz.github.io
---|---
|  Vassilis Vassiliades received the Ph.D. degree from the University of
Cyprus (2015). He is currently a team leader at the Research Centre on
Interactive Media, Smart Systems and Emerging Technologies (RISE) in Cyprus.
He held post-doctoral and research engineer positions at Inria Nancy, France
(2015-2018), and research associate positions at the University of Cyprus
(2015-2019) and RISE (2019). His research focuses on reinforcement learning,
neural networks and evolutionary computation.
---|---
|  Freek Stulp received his doctorate degree in Computer Science from the
Technische Universität München in 2007. He is currently the head of the
department of Cognitive Robotics at the Institute of Robotics and Mechatronics
at the German Aerospace Center (DLR). Previously, he was an assistant
professor at the École Nationale Supérieure de Techniques Avancées (ENSTA-
ParisTech). He currently serves as an Associate Editor in IEEE Transactions on
Robotics.
---|---
|  Sylvain Calinon received the Ph.D. degree from the Ecole Polytechnique
Fédérale de Lausanne (EPFL) in 2007. He is a Senior Researcher at the Idiap
Research Institute, and a Lecturer at the EPFL. From 2009 to 2014, he was a
Team Leader at the Department of Advanced Robotics, Italian Institute of
Technology. From 2007 to 2009, he was a Postdoc at EPFL. He currently serves
as an Associate Editor in IEEE Transactions on Robotics and IEEE Robotics and
Automation Letters. Website: http://calinon.ch
---|---
|  Jean-Baptiste Mouret received the Ph.D. degree in 2008 from the Pierre and
Marie Curie University (Paris, France). He is currently a senior researcher
(“Directeur de recherche”) at Inria, the French research institute dedicated
to computer science and mathematics; from 2009 to 2015, he was an assistant
professor (“maître de conférences”) at the Pierre and Marie Curie University.
His work was recently featured on the cover of Nature (Cully et al., 2015) and
it received several national and international scientific awards, including
the “Prix La Recherche 2016” and the “Distinguished Young Investigator in
Artificial Life 2017”. Website: http://members.loria.fr/jbmouret
---|---


# Influence of Co and Mn on Electronic and Magnetic properties of
$\mathrm{Ni_{2}MnGa}$ Heusler alloy

Karunakaran M  Rudra Banerjee<EMAIL_ADDRESS>Department of Physics and
Nanotechnology, SRM Institute of Science and Technology, Kattankuhur, Tamil
Nadu, 603203, India

###### Abstract

The ferromagnetic Heusler alloy $\mathrm{Ni_{2}MnGa}$ had been of major
interest in past few years because of its magnetic properties which can be
easily tuned. $\mathrm{Ni_{2}MnGa}$ Heusler alloys are intermetallic alloy
with $L2_{1}$ structure. Here we report a detailed investigation of effect of
doping of Co and Mn in $\mathrm{Ni_{2}MnGa}$. Magnetic properties and
electronic structure of $\mathrm{Ni_{2-x}Co_{x}MnGa_{1-z}Mn_{z}}$ Heusler
alloys have been studied by using Green’s function based KKR-CPA method based
DFT calcuions. We will show the magnetisation can be tuned depending on the
cite Mn occupies. We will also discuss the magnetic interactions and magnetic
stability of the systems.

##  I Introduction

Heusler alloys (HAs) has gained a lot of attention from the community due to
their possible usage in magnetocaloric(MC) [1], magnetic shape memory[2],
spintronics[3] and GMR devices[4, 5, 6]. These ternary alloys, with generic
symbol $X_{2}YZ$, where $X$ and $Y$ are $d$-block and $Z$ is main group
elements, shows various structures[7]. The completely ordered phase
demonstrates $L2_{1}$ structure. $Y-Z$ disordered phase demonstrates $B_{2}$
phase and the complete $X-Y-Z$ disorder demonstrates $A_{2}$ structure. Wycoff
positions 4a(0, 0, 0), 4b(1/2, 1/2, 1/2) and 8c(1/4, 1/4, 1/4) are occupied by
$Z$, $Y$ and $X$ respectively[8, 9, 10] shown in figure (1). As we are having
more lab-scale results, both theoretical and experimental, the commercial
implementation is getting nearer. The downsize of the material with various
applications like topological insuors, magnetic cooling and spintroncis
becomes imperative[11, 12].

Among the HAs, Ni-Mn-Ga system has got much attention due to their
demonstration of both magnetic and structural phase transition[13, 14]. Ni-Mn-
Ga systems also show properties like shape memory effect (SME) and magnetic
field induced strain (MFIS) favourable to be used as actuators and sensors[15,
16]. They have also shown very favourable MC properties, suitable to replace
the century old Joule-Thomson cooling[17, 18]. Ni-Mn-Ga system has a
magnetostructural phase transition with huge change in isothermal magnetic
entropy ($\Delta S_{m}$), but reively narrow range of working temperature
($T_{fwhw}$), typical to the materials with first order phase
transition(FOPT)[19]. However, the thermal hysteresis of FOPT and the Curie
temperature is detrimental for magnetic refrigeration.

One of the best characteristics of HAs are the extreme tunability of their
various properties due to doping. After Groot’s[20] theoretical observation of
half-metalicity in half-heusler $\mathrm{NiMnSb}$ alloy, the same was
demonstrated very soon in full-heusler $L2_{1}$ structure[21, 22, 23] as well.
Half-metallic properties has also been proposed in Co doped
$\mathrm{Ni_{2-x}Co_{x}MnGa}$[24] and $\mathrm{Ni_{2-x}Co_{x}MnAl}$[25]
systems. Off-stoichiometric $\mathrm{Ni_{2-x}Q_{x}MnGa}$ and
$\mathrm{Ni_{2-x}Q_{x}MnGa_{1-z}R_{z}}$, where Q and R is any suitable
elements has used to improve the various properties of $\mathrm{Ni_{2}MnGa}$.
$\mathrm{Ni_{2-x}Q_{x}MnGa}$/Al is under study for long time for their
martensite phase transition ($M_{s}$) and magnetic critical
temperature($T_{C}$). It is observed that MC effect is highest around
$0.18\leqslant x\leqslant 0.27$ as magnetic and structural phase change occurs
in this range[6, 26]. Partially substituting Ni with Co is known to enhance
the ferromagnetic coupling and hence the $T_{C}$, and decreases the
martensitic temperature, $T_{M}$. This increases the possibility of that the
quaternary system undergoes a martensitic transition together with a meta
magnetic phase transition[27].

In General, both structural and magnetic properties are highly depends on off-
stoichiometric combination of $Z$. There are also a few study to observe the
effect of $Y$ replacing $Z$, mostly $Ni_{2}MnZ_{1-y}Mn$, $Z$=Ga,In etc[28, 27,
5, 29]. Substituting the Mn element in $Z$ position, can stabilize the
ferromagnetic cubic state[5] and these MnZ atom interact antiferromagnetically
between the surrounding MnZ and normal MnY because the distance between the
MnY-MnZ is less compared to the MnY-MnY and MnZ-MnZ[30]. Here, MnZ is denotes
the Mn atom in Ga site (Z), and MnY is the Mn’s normal position (Y). We
expected to enhance the magnetic moments of the material for doping of Mn in
Ga site.

In the present article, we have reported our systematic study of
$\mathrm{Ni_{2-x}Co_{x}MnGa_{1-z}Mn_{z}}$ (which will be referred as
NiMnGa($x,z$), where, $x,z$ is the concentration of Co and Mn in Ga site,
respectively, for brevity. For example,
$\mathrm{Ni_{1.88}Co_{0.12}MnGa_{0.74}Mn_{0.26}}$ will be denoted as
NiMnGa($0.12,0.26$)). We have studied the variation of electronic structure,
magnetic exchange and magnetic critical temperature of NiMnGa($x,z$) with $x$
ranging from 0 to 0.26 and $z$ ranging from 0 to 0.5. This is technologically
relevant as the $T_{C}$ is one of the detrimental factor of working
temperature range.

Figure 1: Crystal structure of the full Heusler alloys $X_{2}YZ$ (  $X$ atom,
$Y$ atom,  $Z$ atom)

##  II Methods

The off-stoichiometric composition are best handled like
$\mathrm{Ni_{2-x}Q_{x}MnGa_{1-z}R_{z}}$are best handled using greens function
based formalism[31]. We have performed ab-initio calcuions using multiple
scattering greens function formalism as implemented in SPRKKR code[32, 33,
34]. We have used the Perdew-Burke-Ernzerhof(PBE)[35] method is used as GGA
functional. First Brillouin zone integrations were performed on the 250 grid
of k-points and energy convergence criteria was set at the calcuion is
$10^{-5}$Ry. We have implemented full potential spin polarised scaler
relativistic implimentation of SPRKKR with angular momentum cut-off
$\ell_{max}=2$ as suitable for our system.

The Lattice parameter with minimum energy of
$\mathrm{Ni_{2}MnGa}$(NiMnGa($0,0$)) is calculated using the following
procedure:  (i) obtained the lattice parameter from materials project
database[36];  (ii) calcue the scf of the system, with varying lattice
parameter ranging from 94% to 106%, with identical calculation for each
lattice parameters.  (iii) fit the lattice parameter vs energy plot obtained
in last step using a 4th order polynomial.  The minima of the curve is the
optimized lattice parameter. For NiMnGa($x,z$), we have taken the optimized
lattice parameter of previous calculations that has a minimum $x,z$ change as
the starting point and followed the steps above. We have shown the
optimization curve of $\mathrm{Ni_{2}MnGa}$in figure (2(a)). Other
minimization energy curves are not shown here due to brevity.

The magnetic exchange energy ($\mathcal{J}^{\nu\mu}_{ij}$) was calcued to
understand the properties of magnetic interactions. The Heisenberg model,
defined as

$H=\sum_{\nu,\mu}\sum_{i,j}\mathcal{J}^{\nu\mu}_{ij}\mathbf{e}_{i}^{\nu}\cdot\mathbf{e}_{j}^{\mu}$
(1)

where $\nu,\mu$ represent atoms in different subtices, $i,j$ is different
lattice point, $\mathbf{e^{\nu}_{i}}$ is the magnetic orientation of $i$th
atom at $\nu$ subtice. The $\mathcal{J}_{i,j}^{\nu\mu}$ is calculated via the
energy difference due to infinitesimal change of magnetic direction, as
formued by Lichtenstein[37].

The magnetic relaxation, as revealed by Gilbert damping parameter ($\alpha$)
is calcued using Landau-Lifshitz-Gilbert equation[38],

$\partial_{t}\mathbf{m}=-\gamma\mathbf{m}\times\mathbf{H}_{\mathrm{eff}}+\mathbf{m(r)}\times\int
d\mathbf{r}^{\prime}\mathbf{\alpha(r,r^{\prime})}\partial_{t}\mathbf{m(r^{\prime})}$
(2)

where $\mathbf{m}\equiv\mathbf{m(r},t)$ is the unit vector along the direction
of magnetization, $\gamma$ is the gyromagnetic ratio, $H_{\mathbf{eff}}$ is
the effective field, and $\alpha$, the magnetic damping, is the 2nd rank
tensor with $\alpha^{ij}=\alpha^{ij}(\mathbf{m,\Delta m})$.

Finally, the $\mathrm{T_{C}}$ is estimated using mean field theory, yielding

$k_{B}T_{C}=\frac{3}{2}\mathcal{J}_{l}$ (3)

where $\mathcal{J}_{l}$ is the largest eigen value of the determinant, as
described in [30]. It must be remembered that mean field calcuions generally
overestimate the $T_{C}$, which is true for our calculations, as well.

(a) Lattice parameter optimization of $\mathrm{Ni_{2}MnGa}$.

(b) Variation of lattice parameter of NiMnGa($x,z$) with $x,z$

Figure 2: Lattice parameter calculations for NiMnGa($x,z$). In fig. (2(a)),
the blue dots are the actual calculation points; orange dot is the lattice
parameter for minimal energy.

##  III Results

We have calcued the electronic and magnetic properties of NiMnGa($x,z$) for
$x=0.0,0.12,0.24$ and $z=0.0,0.15,0.26,0.35,0.5$. The optimised lattice
parameter of each sample calculated using the method described in the sec.
(II) are tabued in Table (4). The Figure (2(b)) shows the variation of lattice
parameters with $x$ and $z$, which is mostly linear.

###  III.1 Pure $\mathrm{Ni_{2}MnGa}$ system

Figure (3) shows the density of states (DOS) (Figure (3(a))) and magnetic
exchange interactions ($\mathcal{J}_{ij}$) (Figure (3(b))) with Mn atom at the
center for $\mathrm{Ni_{2}MnGa}$. The $\mathcal{J}_{ij}$ is highest for Mn-Ni
interactions, with around 5meV. All the interactions are predominantly
ferromagnetic in this case. The complete table of total and individual
magnetic moment per atom is tabued in Table (4).

(a) DOS of $\mathrm{Ni_{2}MnGa}$

(b) $\mathcal{J}_{ij}$ of $\mathrm{Ni_{2}MnGa}$

Figure 3: The DOS and $\mathcal{J}_{ij}$ of $\mathrm{Ni_{2}MnGa}$.

###  III.2 NiMnGa($x,0$) systems

The NiMnGa($x,0$) system, i.e. doping in Ni site only, though heavily studied,
we have included them for completeness. The electronic and magnetic structures
of NiMnGa($0.12,0$) is shown in Figure (4) as representative of the series.
Other NiMnGa($x,z$) plots are not shown due to brevity. The Co atom has much
higher magnetic moment ($\approx 1\mu B$) compared to Ni ($\approx 0.3\mu B$),
and the $\mathcal{J}_{ij}$ is much higher compared to $\mathrm{Ni_{2}MnGa}$.

(a) DOS of NiMnGa($0.12,0$)

(b) $\mathcal{J}_{ij}$ of NiMnGa($0.12,0$)

Figure 4: The DOS and Exchange interaction of NiMnGa($0.12,0$).

Due to very fragile magnetic moment of Ni,it is seen that moment of Ni
increases rapidly with increase of doping concentrations of Co, as shown in
Table (1). Off course, in current case, MnZ =0.

Table 1: Magnetic moments on individual atoms for NiMnGa($x,0$) system Concentrations | Magnetic Moment($\mu B$)
---|---
$x$ | $z$ | NiMnGa($x,0$) | Ni | Co | MnZ | Ga | MnY
0 | 0 | 4.0805 | 0.2997 | 0 | 0 | -0.0807 | 3.5617
0.12 | 0 | 4.2151 | 0.3134 | 1.0107 | 0 | -0.0859 | 3.5068
0.24 | 0 | 4.3452 | 0.3316 | 0.9985 | 0 | -0.0907 | 3.4527

###  III.3 NiMnGa($0,z$) systems

The opposite systems of the one shown in sec. (LABEL:sec:x,0_systems) is
NiMnGa($0,z$) system, where $X$ site is fully ordered, but Ga is replaced by
Mn, yielding $\mathrm{Ni_{2-x}Co_{x}MnGa_{1-z}Mn_{z}}$ system. Doping Mn in
$Z$ position will decreases the lattice parameter of the material figure
(2(b)) and table (4). As before, only one representative figure for electronic
and magnetic interaction is plotted in Figure (5). The electronic structure
doesn’t change much due to doping(Figure (5(a))), but the absence of Co in Ni
site decreases the magnetic interactions of Mn-Ni by half w.r.t. the case of
NiMnGa($x,0$)(Figure (5(b))).

The most dominant interaction is MnY-MnZ interaction, which are
antiferromagnetically coupled.

(a) DOS of NiMnGa($0,0.26$)

(b) $\mathcal{J}_{ij}$ of NiMnGa($0,0.26$)

Figure 5: The DOS and Exchange interaction of NiMnGa($0,0.26$). Table 2: Magnetic moments on individual atoms for NiMnGa($0,z$) system. Concentrations | Magnetic Moment($\mu B$)
---|---
$x$ | $z$ | NiMnGa($0,z$) | Ni | Co | MnZ | Ga | MnY
0 | 0.15 | 4.7603 | 0.3573 | 0 | 3.6542 | -0.0861 | 3.5708
0 | 0.26 | 5.2739 | 0.4043 | 0 | 3.6567 | -0.0891 | 3.5805
0 | 0.35 | 5.6955 | 0.4433 | 0 | 3.6567 | -0.0913 | 3.5885

###  III.4 Complete disorder: NiMnGa($x,z$) systems

Finally, we study the systems with disorder both at X and Z cite, i.e.
$\mathrm{Ni_{2-x}Co_{x}MnGa_{1-z}Mn_{z}}$. Figure (6) and Figure (7)
represents the electronic and magnetic structures of NiMnGa($0.12,0.26$) and
NiMnGa($0.24,0.5$) respectively.

(a) DOS of NiMnGa($0.12,0.26$)

(b) $\mathcal{J}_{ij}$ of NiMnGa($0.12,0.26$)

Figure 6: The DOS and Exchange interaction for NiMnGa($0.12,0.26$)

(a) DOS of NiMnGa($0.24,0.50$)

(b) $\mathcal{J}_{ij}$ of NiMnGa($0.24,0.50$)

Figure 7: The DOS and Exchange interaction for NiMnGa($0.24,0.50$)

The variation of magnetic moments of NiMnGa($x,z$) are tabued in Table (3).
This shows the presence of Co at Ni induces magnetism in a negative way to all
other atoms except the Ni, and their total magnetic moment is changed
accordingly.

Table 3: Magnetic moments on individual atoms for NiMnGa($0,z$) system. Concentrations | Magnetic Moment($\mu B$)
---|---
$x$ | $z$ | NiMnGa($x,z$) | Ni | Co | MnZ | Ga | MnY
0.12 | 0.15 | 4.9006 | 0.3745 | 1.0861 | 3.6002 | -0.0912 | 3.5184
0.26 | 5.4133 | 0.4239 | 1.1445 | 3.6014 | -0.0942 | 3.5260
0.35 | 5.8264 | 0.4629 | 1.1893 | 3.5986 | -0.0968 | 3.5297
0.50 | 6.5204 | 0.5260 | 1.2626 | 3.5980 | -0.1019 | 3.5435
0.24 | 0.15 | 5.0373 | 0.3981 | 1.0806 | 3.5444 | -0.0956 | 3.4631
0.26 | 5.5422 | 0.4460 | 1.1368 | 3.5457 | -0.0994 | 3.4703
0.35 | 5.9575 | 0.4856 | 1.1826 | 3.5460 | -0.1025 | 3.4773
0.50 | 6.6384 | 0.5485 | 1.2545 | 3.5415 | -0.1085 | 3.4860

Replacing Ga with MnZ on the other hand has a positive effect on magnetic
moment. Specially, Mn${}_{Z}=50\%$ gives a comparatively large magnetic moment
compare to others.

Table 4: We reported the Magnetic moment, Curie temperature($T_{c}$) and Gilbert value of the different X and Z configuration of NiMnGa($x,z$) Combination | Lattice | Megnetic | $T_{C}$ | $\alpha$
---|---|---|---|---
$x$ | $z$ | Parameter (au) | Moment($\mu_{B}$) | (K) | ($10^{-3}$)
0 | 0 | 11.0608 | 4.0805 | 392.975 | 0.000
0.15 | 11.0745 | 4.7603 | 741.170 | 8.645
0.26 | 11.0840 | 5.2739 | 690.785 | 12.097
0.35 | 11.0909 | 5.6955 | 664.351 | 14.390
0.50 | 11.1005 | 6.3938 | 622.667 | 16.710
0.12 | 0 | 11.0472 | 4.2151 | 1138.093 | 6.460
0.15 | 11.0615 | 4.9006 | 1374.932 | 10.396
0.26 | 11.0663 | 5.4133 | 1464.969 | 12.448
0.35 | 11.0731 | 5.8264 | 1538.180 | 13.820
0.50 | 11.0819 | 6.5204 | 1653.155 | 15.422
0.24 | 0 | 11.0322 | 4.3452 | 1147.953 | 9.178
0.15 | 11.0419 | 5.0373 | 1388.860 | 11.985
0.26 | 11.0479 | 5.5422 | 1476.459 | 13.394
0.35 | 11.0543 | 5.9575 | 1548.840 | 14.103
0.50 | 11.0627 | 6.6384 | 1661.804 | 14.430

The variation of magnetic moment of NiMnGa($x,z$) as a function of doping
concentration is shown in Figure (8(a)). The change in magnetic moment in the
system changes very linearly within the doping range.

(a) Variation of magnetic moment with doping concentration.

(b) Variation of $T_{C}$ with doping concentration.

Figure 8: Variation of magnetic moment and $T_{C}$ of NiMnGa($x,z$)

This well behaved nature is good for tuning and applicability of this system.

Finally, we have done mean field based calcuions to find the $T_{C}$ of the
systems. Our result shows that doping both $X$ and $Z$ sites increases the
$T_{C}$. But, for NiMnGa($0,z$), $T_{C}$ is highest around $z=0.12$, then it
starts decreasing. On the other hand, doping a little in $X$ site, increases
the $T_{C}$ and increases steadily. Our $T_{C}$ for NiMnGa($0,0$)
commensurates with [39]. Unfortunately we are unable to find any literature on
variation of $T_{C}$ with doping. The variation of $T_{C}$ is shown in Figure
(8(b)) and tabued in Table (4).

##  IV Conclusion

We have investigated the electronic and magnetic properties of
$\mathrm{Ni_{2-x}Co_{x}MnGa_{1-z}Mn_{z}}$, with $0\leqslant x\leqslant 0.24$,
$0\leqslant z\leqslant 0.5$. The study of replacing both X and Z is
significantly sparse compared to any one of them. We have used DFT and Mean
Field approximations to study the compound effect of dual-doping. Electronic
structures, DOS and magnetic properties including magnetic exchange and
moments are calcued using DFT as implemented in SPR-KKR package. The critical
temperature, $T_{C}$ is calculated using MFA. Our calcuion shows the existence
of strong Mn-Mn antiferromagnetic interaction between Mn at Ga and Mn at its
own subtice. It is also noticed that changing concentration at MnZ does not
change the magnetic exchange significantly. This is an interesting result, as
opposed to the findings in the case of $Ni_{2}Mn_{1+x}Sn_{1-x}$[30].

The $T_{C}$s are obtained using MFA simuions using the magnetic exchange
values obtained from ab-initio calculations. The numerical results for pure
system is very close to that reported by experiments.

We would like to point out the calcuions are done for martensetic phases.
Nevertheless, the dual-doping case, we have shown the variation of magnetic
exchange parameters and magnetic critical temperatures. We believe, these
findings will be helpful for designing the functional properties like MCE and
shape memory alloys by alloying $\mathrm{Ni_{2}MnGa}$ suitably.

##  V Acknowledgement

We Acknowledge the High Performance Computing Center (HPCC), SRM IST for
providing the computational facility to carry out this research work
effectively.

## References

  * Krenke  _et al._ [2005] T. Krenke, E. Duman, M. Acet, E. F. Wassermann, X. Moya, L. Mañosa, and A. Planes, Inverse magnetocaloric effect in ferromagnetic ni–mn–sn alloys, Nature Materials 4, 450 (2005).
  * Sozinov  _et al._ [2002] A. Sozinov, A. A. Likhachev, N. Lanska, and K. Ullakko, Giant magnetic-field-induced strain in NiMnGa seven-layered martensitic phase, Applied Physics Letters 80, 1746 (2002).
  * Bainsla  _et al._ [2015] L. Bainsla, A. I. Mallick, M. M. Raja, A. A. Coelho, A. K. Nigam, D. D. Johnson, A. Alam, and K. G. Suresh, Origin of spin gapless semiconductor behavior in CoFeCrGa: Theory and experiment, Physical Review B 92, 10.1103/physrevb.92.045201 (2015).
  * Kundu  _et al._ [2017] A. Kundu, S. Ghosh, R. Banerjee, S. Ghosh, and B. Sanyal, New quaternary half-metallic ferromagnets with large curie temperatures, Scientific Reports 7, 10.1038/s41598-017-01782-5 (2017).
  * Kihara  _et al._ [2021] T. Kihara, T. Roy, X. Xu, A. Miyake, M. Tsujikawa, H. Mitamura, M. Tokunaga, Y. Adachi, T. Eto, and T. Kanomata, Observation of inverse magnetocaloric effect in magnetic-field-induced austenite phase of heusler alloys ni50-xCoxMn31.5ga18.5(x=9 and 9.7), Physical Review Materials 5, 10.1103/physrevmaterials.5.034416 (2021).
  * Liu  _et al._ [2012] J. Liu, N. Scheerbaum, S. Kauffmann-Weiss, and O. Gutfleisch, NiMn-based alloys and composites for magnetically controlled dampers and actuators, Advanced Engineering Materials 14, 653 (2012).
  * Manna  _et al._ [2018] K. Manna, Y. Sun, L. Muechler, J. Kübler, and C. Felser, Heusler, weyl and berry, Nature Reviews Materials 3, https://doi.org/10.1038/s41578-018-0036-5 (2018).
  * Yang  _et al._ [2019] X. Yang, Y. Wang, M. Du, and Y. Xue, First-principles study of pt doping effects on ni2mnga and ni2fega ferromagnetic shape memory alloys, J. Appl. Phys. 126, https://doi.org/10.1063/1.5116247 (2019).
  * Heusler [1903] F. Heusler, Ueber magnetische manganlegierungen, Phys. Ges 5 (1903).
  * Heusler  _et al._ [1903] F. Heusler, W. Starck, and E. Haupt, Magnetisch-chemische studien, Phys. Ges 5 (1903).
  * Wang  _et al._ [2014] C. Wang, J. Meyer, N. Teichert, A. Auge, E. Rausch, B. Balke, A. HÃ¼tten, G. H. Fecher, and C. Felser, Heusler nanoparticles for spintronics and ferromagnetic shape memory alloys, Journal of Vacuum Science & Technology B, Nanotechnology and Microelectronics: Materials, Processing, Measurement, and Phenomena 32, 020802 (2014).
  * Zhang  _et al._ [2020] Y. Zhang, F. Qin, D. Estevez, V. Franco, and H. Peng, Structure, magnetic and magnetocaloric properties of ni2mnga heusler alloy nanowires, Journal of Magnetism and Magnetic Materials 513, 167100 (2020).
  * Webster  _et al._ [1984] P. J. Webster, K. R. A. Ziebeck, S. L. Town, and M. S. Peak, Magnetic order and phase transformation in ni2mnga, Philosophical Magazine B 49, 10.1080/13642817408246515 (1984).
  * Dilmieva  _et al._ [2017] E. T. Dilmieva, Y. S. Koshkid’ko, V. V. Koledov, A. P. Kamantsev, A. V. Mashirov, J. Cwik, V. V. Khovaylo, and V. G. Shavrov, Formation of a martensitic twins structure in ni2.16mn0.84ga heusler alloy by high magnetic fields under adiabatic and isothermal conditions, Bulletin of the Russian Academy of Sciences: Physics 81, 1283 (2017).
  * Wuttig  _et al._ [2000] M. Wuttig, L. Liu, K. Tsuchiya, and R. D. James, Occurrence of ferromagnetic shape memory alloys (invited), Journal of Applied Physics 87, 4707 (2000).
  * Xin  _et al._ [2007] Y. Xin, Y. Li, L. Chai, and H. Xu, Shape memory characteristics of dual-phase ni–mn–ga based high temperature shape memory alloys, Scripta Materialia 57, 599 (2007).
  * Kamantsev  _et al._ [2015] A. P. Kamantsev, V. V. Koledov, A. V. Mashirov, E. T. Dilmieva, V. G. Shavrov, J. Cwik, A. S. Los, V. I. Nizhankovskii, K. Rogacki, I. S. Tereshina, Y. S. Koshkid'ko, M. V. Lyange, V. V. Khovaylo, and P. Ari-Gur, Magnetocaloric and thermomagnetic properties of ni2.18mn0.82ga heusler alloy in high magnetic fields up to 140 kOe, Journal of Applied Physics 117, 163903 (2015).
  * Khovaylo  _et al._ [2008] V. V. Khovaylo, K. P. Skokov, Y. S. Koshkid’ko, V. V. Koledov, V. G. Shavrov, V. D. Buchelnikov, S. V. Taskaev, H. Miki, T. Takagi, and A. N. Vasiliev, Adiabatic temperature change at first-order magnetic phase transitions:ni2.19mn0.81gaas a case study, Physical Review B 78, 10.1103/physrevb.78.060403 (2008).
  * Umetsu  _et al._ [2016] R. Y. Umetsu, X. Xu, and R. Kainuma, NiMn-based metamagnetic shape memory alloys, Scripta Materialia 116, 1 (2016).
  * de Groot  _et al._ [1983] R. A. de Groot, F. M. Mueller, P. G. van Engen, and K. H. J. Buschow, New class of materials: Half-metallic ferromagnets, Physical Review Letters 50, 2024 (1983).
  * Ishida  _et al._ [1982] S. Ishida, S. Akazawa, Y. Kubo, and J. Ishida, Band theory of Co2MnSn, Co2TiSn and Co2TiAl, Journal of Physics F: Metal Physics 12, 1111 (1982).
  * Fujii  _et al._ [1990] S. Fujii, S. Sugimura, Ishida, and S. Asano, Hyperfine fields and electronic structures of the heusler alloys Co2MnX (X=Al, Ga, Si, Ge, Sn), Journal of Physics: Condensed Matter 2, 8583 (1990).
  * Jourdan  _et al._ [2014] M. Jourdan, J. Minár, J. Braun, A. Kronenberg, S. Chadov, B. Balke, A. Gloskovskii, M. Kolbe, H. Elmers, G. Schönhense, H. Ebert, C. Felser, and M. Kläui, Direct observation of half-metallicity in the heusler compound co2mnsi, Nature Communications 5, 10.1038/ncomms4974 (2014).
  * Okubo  _et al._ [2011] A. Okubo, X. Xu, R. Y. Umetsu, T. Kanomata, K. Ishida, and R. Kainuma, Magnetic properties of co50-xNixMn25al25 alloys with b2 structure, Journal of Applied Physics 109, 07B114 (2011).
  * Kanomata  _et al._ [2009] T. Kanomata, Y. Kitsunai, K. Sano, Y. Furutani, H. Nishihara, R. Y. Umetsu, R. Kainuma, Y. Miura, and M. Shirai, Magnetic properties of quaternary heusler alloysNi2-xCoxMnGa, Physical Review B 80, 10.1103/physrevb.80.214402 (2009).
  * Buchelnikov and Sokolovskiy [2011] V. D. Buchelnikov and V. V. Sokolovskiy, Magnetocaloric effect in ni-mn-x (x = ga, in, sn, sb) heusler alloys, The Physics of Metals and Metallography 112, 633 (2011).
  * Kainuma  _et al._ [2006] R. Kainuma, Y. Imano, W. Ito, H. Morito, Y. Sutou, K. Oikawa, A. Fujita, K. Ishida, S. Okamoto, O. Kitakami, and T. Kanomata, Metamagnetic shape memory effect in a heusler-type ni43co7mn39sn11 polycrystalline alloy, Applied Physics Letters 88, 192513 (2006).
  * Orlandi  _et al._ [2020] F. Orlandi, A. Çakır, P. Manuel, D. D. Khalyavin, M. Acet, and L. Righi, Neutron diffraction and symmetry analysis of the martensitic transformation in co-doped ni2mnga, Phys. Rev. B 101, 10.1103/physrevb.101.094105 (2020).
  * Sokolovskiy  _et al._ [2014] V. Sokolovskiy, A. GrÃ¼nebohm, V. Buchelnikov, and P. Entel, Ab initio and monte carlo approaches for the magnetocaloric effect in co- and in-doped ni-mn-ga heusler alloys, Entropy 16, 4992 (2014).
  * Sokolovskiy  _et al._ [2012] V. V. Sokolovskiy, V. D. Buchelnikov, M. A. Zagrebin, P. Entel, S. Sahoo, and M. Ogura, First-principles investigation of chemical and structural disorder in magnetic ni2mn1+xsn1-xheusler alloys, Physical Review B 86, 10.1103/physrevb.86.134418 (2012).
  * Banerjee and Mookerjee [2010] R. Banerjee and A. Mookerjee, Augmented space recursion code and application in simple binary metallic alloy, International Journal of Modern Physics C 21, 205 (2010).
  * Ebert  _et al._ [2011] H. Ebert, D. Ködderitzsch, and J. Minár, Calculating condensed matter properties using the KKR-green's function method—recent developments and applications, Reports on Progress in Physics 74, 10.1088/0034-4885/74/9/096501 (2011).
  * Ebert [2000] H. Ebert, Fully relativistic band structure calculations for magnetic solids - formalism and application, in  _Electronic Structure and Physical Properies of Solids_ , edited by H. Dreyssé (Springer Berlin Heidelberg, Berlin, Heidelberg, 2000) pp. 191–246.
  * Mavropoulos Phivos [2006] P. N. Mavropoulos Phivos, The korringa-kohn-rostoker (kkr) green function method i. electronic structure of periodic systems, in  _Computational Nanoscience: Do It Yourself!_ , edited by D. M. J. Grotendorst, S. Blügel (2006).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Physical Review Letters 77, 3865 (1996).
  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, Commentary: The materials project: A materials genome approach to accelerating materials innovation, APL Materials 1, 011002 (2013).
  * Liechtenstein  _et al._ [1987] A. Liechtenstein, M. Katsnelson, V. Antropov, and V. Gubanov, Local spin density functional approach to the theory of exchange interactions in ferromagnetic metals and alloys, Journal of Magnetism and Magnetic Materials 67, https://doi.org/10.1016/0304-8853(87)90721-9 (1987).
  * Skadsem  _et al._ [2007] H. J. Skadsem, Y. Tserkovnyak, A. Brataas, and G. E. W. Bauer, Magnetization damping in a local-density approximation, Physical Review B 75, 10.1103/physrevb.75.094416 (2007).
  * Brown  _et al._ [2002] P. J. Brown, J. Crangle, T. Kanomata, M. Matsumoto, K.-U. Neumann, B. Ouladdiaf, and K. R. A. Ziebeck, The crystal structure and phase transitions of the magnetic shape memory compound Ni2MnGa, Journal of Physics: Condensed Matter 14, 10159 (2002).