# Quantum probabilities for the causal ordering of events

Charis<EMAIL_ADDRESS>
Department of Physics, University of Patras, 26500 Greece
and
Maria-Electra Plakitsi
Department of Electrical and Computing Engineering, University of Patras,
26500 Greece

###### Abstract

We develop a new formalism for constructing probabilities associated to the
causal ordering of events in quantum theory, where by an event we mean the
emergence of a measurement record on a detector. We start with constructing
probabilities for the causal ordering events in classical physics, where
events are defined in terms of worldline coincidences. Then, we show how these
notions generalize to quantum systems, where there exists no fundamental
notion of trajectory. The probabilities constructed here are experimentally
accessible, at least in principle. Our analysis here clarifies that the
existence of quantum orderings of events does not require quantum gravity
effects: it is a consequence of the quantum dynamics of matter, and it appears
in presence of a fixed background spacetime.

##  1 Introduction

A bet on any type of race (with humans, horses, chariots, or cars) is
equivalent to the assignment of probabilities to a causal ordering of events.
The relevant events are the crossings of the finish line by the runners, and
the causal ordering of such events is the results of the race. In this sense,
assigning probabilities to causal orderings is both one of the oldest
applications of probabilistic thinking, dating at least to the ancient
Olympics, and one of the most common uses of probability theory today. In this
paper, we describe causal ordering of events (COoE) for quantum systems, where
by an event we mean the emergence of a macroscopic measurement record that is
localized in space and in time [1, 2]. We construct the probabilities for such
causal orderings, and we suggest physical set-ups where such probabilities can
be measured, at least in principle.

This work is partially motivated by the recent studies of indefinite causal
ordering of events in quantum computing [3, 4]. In this context, the word
“event” is used to denote an operation on a quantum system, for example, a
step in an algorithm. An indefinite sequence of operations can arguably lead
to significant advances in quantum computation and other technologies [5, 6,
7, 8]. The most common set-up to witness such phenomena involves the quantum
switch, that is, a quantum operation in which two or more quantum channels act
on a quantum system with the order of application determined by the state of
another quantum system. Systems that manifest indefinite causal order in this
sense have been realized in the laboratory [9, 10, 11].

This quantum informational notion of causal ordering differs and may even
conflict with relativistic causality—for a detailed analysis of this issue,
see [12, 13]. The meaning of the term “event” is crucial in this context. In
this paper, we employ a notion of event that is similar to the crossing of the
finishing line by a racer. This notion closely reflects the notion of an event
in relativity, where physical events are invariantly defined in terms of
world-line coincidences222Einstein emphasized this perspective in his first
review of General Relativity [14], see also [15] and [16].. For example, a
particle-detection event is defined as the intersection of the particle’s and
the detector’s worldlines. With this definition, we use the lightcone
structure of spacetime, in order to define causal relations between events.
Since events are defined in terms of trajectories, the dynamical behavior of
trajectories directly influences the properties of causal ordering. If the
trajectories are stochastic, then the ordering of events is also a stochastic
variable.

Quantum theory does not admit trajectories as physical observables, so in
quantum systems, we have to define events and their causal ordering
differently. We define an event as the occurrence of a measurement outcome,
that is, the emergence of a measurement record on a detector. This is the most
conservative definition of a quantum event. It is also the most natural one in
the Copenhagen interpretation [17], and the standard use of the term in
particle physics. If the quantum events can be embedded in spacetime, i.e., if
we can associate a spacetime point or region to the emergence of a measurement
record, then we can define quantum probabilities for the COoE that are natural
analogues of the classical ones.

Such a definition requires the treatment of the time associated to an event as
a quantum observable. It is an old result by Pauli [18] that time cannot be
treated as a self-adjoint operator. The only way to have time as an observable
is to represent it by a Positive-Operator-Valued measure (POVM). To this end,
we use the Quantum Temporal Probabilities (QTP) approach that has been
developed for constructing probabilities for temporal observables [1, 19, 20,
21, 2].

The key idea in QTP is to distinguish between the time parameter of
Schrödinger’s equation and the time variable associated to particle detection
[23, 24]. The latter is then treated as a macroscopic quasi-classical variable
associated to the detector degrees of freedom. A quasi-classical variable is a
coarse-grained quantum variable that satisfies appropriate decoherence
conditions, so that its time evolution can be well approximated by classical
equations [25, 26, 27]. Hence, the detector admits a dual description: in
microscopic scales it is described by quantum theory, but its macroscopic
records are expressed in terms of classical spacetime coordinates. The key
point here is that the spacetime coordinates of a quantum event are random
variables, and they can be used to define quantum probabilities for the causal
order of events.

It is important to emphasize that there is no relation between the COoE
considered here and a quantum causal structure of spacetime, as commonly
postulated in quantum gravity research. The quantum behavior of the COoE,
considered here, is due to the quantum nature of matter, and it coexists
peacefully with a fixed background spacetime. In fact, the background
spacetime structure is essential for defining quantum probabilities for the
COoE.

We also show that probabilities for quantum COoEs can also be defined even if
there is no macroscopic record about the time at which the events occur. To
this end, we construct a simple detection model, in which different orderings
of events correspond to different measurement records. Hence, the
probabilities of such records coincide with probabilities for different causal
orders.

The structure of this paper is the following. In Sec. 2, we provide a general
definition of the notion of an event, and explain how we can construct
probabilities for the causal order of events. In Sec. 3, we apply these
definitions to classical physics, including Hamiltonian mechanics and
stochastic processes. In Sec. 4, we define probabilities for the causal
ordering of quantum events, using temporal observables. In Sec. 5, we present
a simple model for the quantum order of events in absence of records about
temporal observables. In Sec. 6, we summarize our results.

##  2 Main concepts

In this section, we present a general characterization of events, and we
identify the mathematical properties that are satisfied by a causal order of
events.

By “event” we mean a uniquely identifiable occurrence with definite
characteristics. In classical mechanics, events are typically defined as the
intersection of two world-lines. For example, one worldline may correspond to
a particle and the other to an observer with a particle detector; their
coincidence is a particle-detection event. We can improve on this description,
by defining an event as the first intersection of a worldline with a specific
time-like surface. We can consider, for example, the worldline of a runner
crossing the world tube of the finish line in a marathon.

However, in quantum theory, definite characteristics are attributed only to
measurement outcomes; trajectories are not observables. For this reason, we
will define events in terms of measurement records. For example, a detection
event is the “click” of a particle detector. Hence, we identify event with
changes in a macroscopic apparatus that denote the occurrence of a
measurement.

Let $E$ be a set of events in the physical system under study. Events are
discrete occurrences, so $E$ is a discrete set. We will denote events by Greek
letters, $\alpha,\beta$, and so on.

The events in the set $E$ may be ordered causally. We say that
$\alpha\prec\beta$, if event $\alpha$ occurs prior to an event $\beta$. A
causal order on $E$ is the consistent assignment of the order relation $\prec$
to pairs of elements of $E$. A causal order satisfies the properties of a
partial-order relation, namely,

  1. 1.

Irreflexivity: It is never true that $\alpha\prec\alpha$.

  2. 2.

Asymmetry: If $\alpha\prec\beta$, then $\beta\prec\alpha$ is false.

  3. 3.

Transitivity: If $\alpha\prec\beta$ and $\beta\prec\gamma$, then
$\alpha\prec\gamma$.

In a partial order, it is not necessary that all pairs of elements are related
with the order relation. Physically, we can distinguish two cases. Some
elements may be simultaneous, in which case we write $\alpha\sim\beta$. Or
they may be uncomparable, in which case we write $\alpha|\beta$. We therefore
define a causal order as a partial order that also include the distinction
between simultaneous and incomparable pairs of elements.

We will denote the set of all possible causal orders on $E$ by $CO(E)$. We
will denote elements of $CO(E)$ by capital Greek letters. For example, in a
set $E$ that consists of two distinct elements $\alpha$ and $\beta$, there are
four possible causal orders

  * •

$M_{1}=\\{\alpha\prec\beta\\}$

  * •

$M_{2}=\\{\beta\prec\alpha,\\}$

  * •

$M_{3}=\\{\alpha|\beta\\}$

  * •

$M_{4}=\\{\alpha\sim\beta\\}$

We say that a causal order defines a time order on $E$, if there exist no pair
$\alpha,\beta\in E$ such that $\alpha|\beta$. We will denote the set of all
time orders on $E$ by $TO(E)$. Clearly $TO(E)\subset CO(E)$.

Ever since Newton, we define the causal ordering of physical events in terms
of the spacetime causal structure. That is, we consider a four-dimensional
manifold $M$ with points $(x^{0},x^{1},x^{2},x^{3})$ that is equipped with a
partial ordering relation $<$ that defines the causal structure of spacetime.

  * •

In non-relativistic physics, $x<y$ if $x^{0}<y^{0}$, and $x\sim y$ if
$x^{0}=y^{0}$. There are no incomparable elements.

  * •

In relativistic physics, $x<y$, if $y$ lies in the future lightcone of $x$,
and $x|y$ if $x$ is spacelike separated from $y$. Spacelike separated events
are incomparable; there are no simultaneous events.

Since all physical events occur in spacetime, we consider embeddings $X$ of
the set of events $E$ into spacetime, that is, onto maps $X:E\rightarrow M$.
Then, the pullback of the spacetime causal structure with respect to $X$
defines a causal order on $E$, that is,

$\displaystyle\alpha\prec\beta,\ \ \text{if}\ \ X(\alpha)<X(\beta)$ (1)

Hence, the physical COoE reflects the causal structure on spacetime333Here, we
associated events with spacetime points. A more general analysis should
associate events with spacetime regions, but it will not be needed in this
paper. .

It is imperative to distinguish COoEs from the causal structure of spacetime.
As long as we ignore gravitational interactions, the latter is fixed an
unchanging. It is defined by the lightcone structure of Minkowski spacetime,
or of any other background spacetime. However, COoEs are not fixed: they can
be stochastic or quantum variables. This is because they depend on the
embeddings $X$, which are themselves stochastic or quantum variables. The
quantum behavior of the COoEs does not require a quantum behavior of
spacetime, as postulated in quantum gravity theories. As a matter of fact, in
quantum gravity proper, we expect to have no external spacetime causal
structure, hence, the definitions of the COoEs given here do not work. The
difficulties that arise from this fact are known as the problem of time in
quantum gravity [28, 29, 30].

##  3 COoE in classical physics

In this section, we construct probabilities for the COoEs for classical
systems, namely, for Hamiltonian systems and for systems described by
stochastic processes.

###  3.1 Classical mechanics

Let $\Gamma$ be the state space of a classical system; we will denote the
elements of $\Gamma$ by $\xi$. By Hamilton’s equation, a system found at $\xi$
at time $t=0$, will evolve to a point $\sigma_{t}(\xi)$ at time $t$. The map
$\sigma_{t}$ is a diffeomorphism.

An event in classical mechanics corresponds to the first intersection of a
state space trajectory with a surface on $\Gamma$. Surfaces of codimension $s$
are locally determined by the vanishing of $s$ functions on $\Gamma$, hence,
we can represent an event $\alpha$ by a set of $s$ functions $F^{i}_{\alpha}$,
where $i=1,2,\ldots s$. A set $E$ of $n$ events consists $n$ such families.
For simplicity, we will consider only surfaces of codimension one in this
paper, so that one event corresponds to a single function on $\Gamma$.

The causal ordering in the set of events is defined through the parameter of
time evolution $t$, which is assumed to coincide with the Newtonian absolute
time. For any event $\alpha$, we define the null set $N_{\alpha}$ of $\alpha$,
as the set of all $\xi\in\Gamma$, such that the equation
$F_{\alpha}[\sigma_{t}(\xi)]=0$ has no solution for all $t\geq 0$. Then, we
define the time $T_{\alpha}$ of the event $\alpha$, as a function
$T_{\alpha}:\Gamma-N_{\alpha}\rightarrow\mathbb{R}^{+}$, such that
$T_{\alpha}(\xi)$ is the smallest positive value of $t$ that solves the
equation $F_{\alpha}[\sigma_{t}(\xi)]=0$. This means that $T_{\alpha}(\xi)$ is
the time it takes a trajectory that starts at $\xi$ to cross the surface
$F_{\alpha}=0$ for the first time.

Suppose that the initial state of the system corresponds to a probability
distribution $\rho(\xi)$. Then, we can construct joint probability
distributions for the times of events

$\displaystyle p(t_{1},t_{2},\ldots,t_{n})=\int
d\xi\rho(\xi)\delta[T_{1}(\xi)-t_{1}]\delta[T_{2}(\xi)-t_{2}]\ldots\delta[T_{n}(\xi)-t_{n}].$
(2)

These probability densities are not normalized to unity. For proper
normalization, we have to include the probability densities for no events,
which corresponds to the null sets $N_{\alpha}$. For example, for $n=2$, we
have the probability densities $p(t_{1},t_{2})$ as above, together with the
probability densities

$\displaystyle p(N_{1},t_{2})=\int
d\xi\chi_{N_{1}}(\xi)\delta[T_{2}(\xi)-t_{2}]$ (3) $\displaystyle
p(t_{1},N_{2})=\int d\xi\chi_{N_{2}}(\xi)\delta[T_{1}(\xi)-t_{1}]$ (4)
$\displaystyle p(N_{1},N_{2})=\int d\xi\chi_{N_{1}}(\xi)\chi_{N_{2}}(\xi).$
(5)

Here, $\chi_{C}$ is the characteristic function of a set $C$.

We can define the following four causal orders for the two events,

  * •

$M_{1}=\\{1\prec 2\\}$ corresponds to $t_{1}<t_{2}$, or $N_{2}$ together with
finite $t_{1}$.

  * •

$M_{2}=\\{2\prec 1\\}$ corresponds to $t_{2}<t_{1}$, or $N_{1}$ together with
finite $t_{2}$.

  * •

$M_{3}=\\{1|2\\}$ corresponds to $N_{1}$ and $N_{2}$.

  * •

$M_{4}=\\{1\sim 2\\}$ corresponds to $t_{1}=t_{2}$.

Then, we obtain the associated probabilities

$\displaystyle p(M_{1})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}dt_{1}\int_{0}^{t_{1}}dt_{2}p(t_{1},t_{2})+\int_{0}^{\infty}dt_{1}p(t_{1},N_{2})$
$\displaystyle p(M_{2})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}dt_{2}\int_{0}^{t_{2}}dt_{1}p(t_{1},t_{2})+\int_{0}^{\infty}dt_{2}p(N_{1},t_{2})$
$\displaystyle p(M_{3})$ $\displaystyle=$ $\displaystyle p(N_{1},N_{2})$
$\displaystyle p(M_{4})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}dtp(t,t).$ (6)

This procedure is straightforwardly generalized to $n$ events.

As an illustration, consider a system of two free particles of mass $m$ in one
dimension, with state space $\Gamma=\\{x_{1},x_{2},p_{1},p_{2}\\}$. We
restrict to $x_{1}\leq 0$ and $x_{2}\leq 0$, and we consider the pair of
events that correspond to either of the two particles crossing the line $x=0$.
Hence, the two functions that define events are $F_{1}=x_{1}$ and
$F_{2}=x_{2}$. The equations of motion are $x_{1}(t)=x_{1}+p_{1}t/m$ and
$x_{2}(t)=x_{2}+p_{2}t/m$. We straightforwardly find that
$N_{1}=\\{(x_{1},x_{2},p_{1},p_{2})|p_{1}\leq 0\\}$, and
$N_{2}=\\{(x_{1},x_{2},p_{1},p_{2})|p_{2}\leq 0\\}$, that is, the particles
never cross the line $x=0$ if they have non positive momentum. Similarly, we
compute the time functions $T_{1}=-mx_{1}/p_{1}$ and $T_{2}=-mx_{2}/p_{2}$.

It is simple to identify the subsets of $\Gamma$ that correspond to the
different causal orders,

$\displaystyle M_{1}$ $\displaystyle=$
$\displaystyle\\{(x_{1},x_{2},p_{1},p_{2})|p_{1}>0,p_{2}>0,x_{1}p_{2}>x_{2}p_{1}\\}\cup\\{(x_{1},x_{2},p_{1},p_{2})|p_{1}>0,p_{2}\leq
0\\},$ $\displaystyle M_{2}$ $\displaystyle=$
$\displaystyle\\{(x_{1},x_{2},p_{1},p_{2})|p_{1}>0,p_{2}>0,x_{2}p_{1}>x_{1}p_{2}\\}\cup\\{(x_{1},x_{2},p_{1},p_{2})|p_{1}\leq
0,p_{1}>0\\},$ $\displaystyle M_{3}$ $\displaystyle=$
$\displaystyle\\{(x_{1},x_{2},p_{1},p_{2})|p_{1}\leq 0,p_{2}\leq 0\\},$
$\displaystyle M_{4}$ $\displaystyle=$
$\displaystyle\\{(x_{1},x_{2},p_{1},p_{2})|p_{1}>0,p_{2}>0,x_{2}p_{1}=x_{1}p_{2}\\}.$
(7)

The associated probabilities are simply $p(M_{i})=\int
d\xi\chi_{M_{i}}(\xi)\rho(\xi)$. Note that $M_{4}$ is a set of measure zero,
so, in general, the associated probability vanishes.

Suppose, for example, that both particles start from $x_{0}<0$, and that they
have the same momentum distribution $f(p)$, so that

$\displaystyle\rho(x_{1},x_{2},p_{1},p_{2})=\delta(x_{1}-x_{0})\delta(x_{2}-x_{0})f(p_{1})f(p_{2}).$
(8)

Then, we compute, $p(M_{1})=p(M_{2})=w_{+}-\frac{1}{2}w_{+}^{2}$, and
$p(M_{3})=(1-w_{+})^{2}$, where $w_{+}=\int_{0}^{\infty}dpf(p)$ is the
fraction of particles with positive momentum.

###  3.2 Stochastic processes

The analysis of Sec. 3.1 passes with little change to classical stochastic
systems. Consider a system characterized by a state space $\Gamma$ with
elements $\xi$. Let us denote by $P(\Gamma)$ the space of paths on $\Gamma$,
that is, of continuous maps from the time interval $[0,T]$ to $\Gamma$. Here,
we are restricting to paths between an initial time $t=0$, and a final time
$t=T$. A stochastic system is described by a probability measure $\mu$ on
$P(\Gamma)$, such that the expectation of any function $A$ of $P(\Gamma)$ is
given by

$\displaystyle\langle A\rangle=\int d\mu[\xi(\cdot)]A[\xi(\cdot)]]$ (9)

Again, an event $\alpha$ is defined by the first intersection of a path with a
surface, and it can be represented by a function $F_{\alpha}$ on $\Gamma$. We
can still define a null space $N_{\alpha}$, and a time function $T_{\alpha}$,
however, in absence of a deterministic evolution law, these objects are
defined on the space of paths $P(\Gamma)$, and not on $\Gamma$. In particular,
we define by $N_{\alpha}$ the subset of $P(\Gamma)$ that consists of paths
$\xi(\cdot)$ for which the equation $F_{\alpha}(\xi(t))=0$ has no solution for
any $t\in[0,T]$; we will denote the complement of $N_{\alpha}$ by
$\bar{N}_{\alpha}$. We also define the time function $T_{\alpha}$ for any path
$\xi(\cdot)\in\bar{N}_{\alpha}$ by setting the value $T_{\alpha}[\xi(\cdot)]$
on a path $\xi(\cdot)$ equal to the smallest value of $t$ such that
$F_{\alpha}(\xi(t))=0$.

The definition of joint probabilities for the times of events proceeds in a
similar way to Sec. 2. For example, the joint probability distribution for $n$
events is

$\displaystyle p(t_{1},t_{2},\ldots,t_{n})=\int
d\mu[\xi(\cdot)]\delta(t_{1}-T_{1}(\xi(\cdot))\delta(t_{2}-T_{2}(\xi(\cdot))\ldots\delta(t_{n}-T_{n}(\xi(\cdot)).$
(10)

The space of paths $P(\Gamma)$ is split into mutually exclusive and exhaustive
subsets, each corresponding to an element of $CO(E)$. For two events, we have
four elements of $CO(E)$, which correspond to the following subsets.

$\displaystyle M_{1}$ $\displaystyle=$ $\displaystyle\\{\gamma\in
P(\Gamma)|T_{1}[\gamma]<T_{2}[\gamma]\\}\cup\left(N_{2}\cap\bar{N}_{1}\right),$
$\displaystyle M_{2}$ $\displaystyle=$ $\displaystyle\\{\gamma\in
P(\Gamma)|T_{2}[\gamma]<T_{1}[\gamma]\\}\cup\left(N_{1}\cap\bar{N}_{2}\right),$
$\displaystyle M_{3}$ $\displaystyle=$ $\displaystyle N_{1}\cap N_{2},$
$\displaystyle M_{4}$ $\displaystyle=$ $\displaystyle\\{\gamma\in
P(\Gamma)|T_{1}[\gamma]=T_{2}[\gamma]\\}.$ (11)

As an example, consider the case of a Wiener process. We have two particles
undergoing Brownian motion, that is, each particle is described by the
evolution of a single-time probability density $\rho$ on $\mathbb{R}$, by

$\displaystyle\frac{\partial\rho}{\partial
t}=\frac{D}{2}\frac{\partial^{2}\rho}{\partial x^{2}},$ (12)

where $D$ is the diffusion constant. For a particle that starts at $x_{0}=-L$,
the probability density of crossing the line $x=0$ is

$\displaystyle f(t)=\frac{1}{\sqrt{2\pi
Dt}}\frac{L}{2t}e^{-\frac{L^{2}}{2Dt}},$ (13)

with the probability of not crossing $x=0$ for any $t\in[0,\infty)$ equal to
$\frac{1}{2}$.

We assume that the two particles move independently, and that the associated
diffusion constants are different, $D_{1}$ and $D_{2}$. (This is possible, for
example, if the particle masses are different.) The joint probability density
that the first crosses $x=0$ at time $t_{1}$ and the second at time $t_{2}$ is
simply $f_{1}(t_{1})f_{2}(t_{2})$, where $f_{i}$ is the probability density
(13), with diffusion constant $D_{i}$. Then, we evaluate

$\displaystyle p(M_{1})$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\arctan\left(\sqrt{D_{1}/D_{2}}\right)+\frac{1}{4},$
$\displaystyle p(M_{2})$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\arctan\left(\sqrt{D_{2}/D_{1}}\right)+\frac{1}{4},$
$\displaystyle p(M_{3})$ $\displaystyle=$ $\displaystyle\frac{1}{4},$ (14)

where we ignored $M_{4}$, as it is of measure zero.

##  4 COoE in quantum systems

In this section, we define probabilities for the COoE in quantum systems.

###  4.1 Probability assignment

For quantum systems, the definition of events in terms of paths crossing a
surface does not work, because paths are not physical observables in quantum
theory. The only meaningful observables are measurement outcomes. In the most
common measurement scheme, namely, von Neumann measurements, the timing of the
measurement events is fixed a priori. Hence, the causal order of events is
also fixed.

We need a measurement scheme that treats the time of an event as a random
variable, if we are to treat the causal order of events as a random variable
quantum mechanically. This is achieved by the QTP approach that was described
in the introduction.

Suppose that we have a particle detector located at a fixed region in space.
Then, via QTP, we can construct a set of positive $\hat{\Pi}(t)$, such that
the probability density of detection at time $t>0$ is
$p(t)=Tr(\hat{\rho}_{0}\hat{\Pi}(t))$, where $\hat{\rho}_{0}$ is the initial
state if the particle. Together with the positive operator $\hat{\Pi}(N)$ of
no detection, the operators $\hat{\Pi}_{t}$ define a POVM.

For example, we can construct a POVM for the time of arrival of a particle of
mass $m$. We assume that the particle moves at a line and that the particle
detector is located at $x=L$. In the momentum basis,

$\displaystyle\langle
k|\hat{\Pi}(t)|k^{\prime}\rangle=\int\frac{dkdk}{2\pi}S(k,k^{\prime})\sqrt{v_{k}v_{k^{\prime}}}e^{i(k-k^{\prime})L-i(\epsilon_{k}-\epsilon_{k^{\prime}})t},$
(15)

where $\epsilon_{k}=\sqrt{m^{2}+k^{2}}$ is the particle’s energy, $v_{k}$ is
the particle’s velocity, and $S(k,k^{\prime})$ is the localization operator,
that is an operator that determines the irreducible spread of the detection
record. The sharpest localization is achieved for $S(k,k^{\prime})=1$. The
operators $\hat{\Pi}(t)$ are not normalized to unity for $t\in[0,\infty)$.
However, if we restrict to quantum states with strictly positive momentum
content, the contribution to the total probability from negative values of $t$
is negligible, and we can consider the normalization of $\hat{\Pi}(t)$ in the
full real line. In this case,

$\displaystyle\int_{-\infty}^{\infty}dt\hat{\Pi}(t)=\hat{I}.$ (16)

For $n$ independent detectors, each detecting a different particle, we can
identify a POVM
$\hat{\Pi}(t_{1},t_{2},\ldots,t_{n})=\hat{\Pi}_{1}(t_{1})\otimes\hat{\Pi}_{2}(t_{2})\otimes\ldots\otimes\hat{\Pi}_{n}(t_{n})$,
where $\hat{\Pi}_{i}(t_{i})$ corresponds to the POVM for the $i$-th detector.
Thus, we can define probability densities $p(t_{1},t_{2},\ldots,t_{n})$ for
the $n$ measurement events, and we can follow the same procedure as in Sec. 3,
in order to obtain probabilities for different causal orders of $n$
measurement events.

For example, for two events, $1$ and $2$, we have the three COoEs
$M_{1}=\\{1\prec 2\\}$, $M_{2}=\\{2\prec 1\\}$, and $M_{3}=\\{1||2\\}$. The
positive operators

$\displaystyle\hat{E}(M_{1})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}dt_{2}\int_{0}^{t_{2}}dt_{1}\hat{\Pi}_{1}(t_{1})\otimes\hat{\Pi}_{2}(t_{2})+[\hat{I}-\hat{\Pi}_{1}(N)]\otimes\hat{\Pi}_{2}(N),$
$\displaystyle\hat{E}(M_{2})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}dt_{1}\int_{0}^{t_{1}}dt_{2}\hat{\Pi}_{1}(t_{1})\otimes\hat{\Pi}_{2}(t_{2})+\hat{\Pi}_{1}(N)[\hat{I}-\hat{\Pi}_{2}(N)],$
$\displaystyle\hat{E}(M_{3})$ $\displaystyle=$
$\displaystyle\hat{\Pi}_{1}(N)\otimes\hat{\Pi}_{2}(N),$ (17)

define a POVM for the causal orders.

Assume that the two events correspond to the detection of identical particles
with the two detectors located at $x_{1}=L_{1}$ and $x_{2}=L_{2}$ from the
source—see Fig. 1. We use the POVM (15) for both $\hat{\Pi}_{1}(t)$ and
$\hat{\Pi}_{2}(t)$. Taking $-\infty$ for the lower bound in the time integral,
we find

$\displaystyle\hat{E}(M_{1})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\hat{I}+\hat{B}$ (18) $\displaystyle\hat{E}(M_{2})$
$\displaystyle=$ $\displaystyle\frac{1}{2}\hat{I}-\hat{B}$ (19)
$\displaystyle\hat{E}(M_{3})$ $\displaystyle=$ $\displaystyle 0$ (20)

where the operator $\hat{B}$ is defined by matrix elements

$\displaystyle\langle
k_{1},k_{2}|\hat{B}|k_{1}^{\prime},k_{2}^{\prime}\rangle$ $\displaystyle=$
$\displaystyle
iS_{1}(k_{1},k_{1}^{\prime})S_{2}(k_{2},k_{2}^{\prime})\sqrt{v_{k_{1}}v_{k_{2}}v_{k_{1}^{\prime}}v_{k_{2}^{\prime}}}$
(21) $\displaystyle\times$ $\displaystyle
e^{i(k_{1}-k_{1}^{\prime})L_{1}-i(k_{2}-k_{2}^{\prime})L_{2}}\delta(\epsilon_{k_{1}}+\epsilon_{k_{2}}-\epsilon_{k_{1}^{\prime}}-\epsilon_{k_{2}^{\prime}})\;\text{PV}\frac{1}{\epsilon_{k_{2}}-\epsilon_{k_{2}^{\prime}}}.$

Here, PV stands for Cauchy principal value.

Figure 1:  A set-up by which to measure the causal order of two events that
correspond to detections at detectors 1 and 2.

Consider a general initial state of the form
$|\psi\rangle=\sum_{i}c_{i}|\psi_{1i}\rangle\otimes|\psi_{2i}\rangle$. The
probability densities associated to the three orders are

$\displaystyle p(M_{1})=\frac{1}{2}+w,\ \ p(M_{2})=\frac{1}{2}-w,\ \
p(M_{3})=0,$ (22)

where the asymmetry

$\displaystyle w=\sum_{ij}c_{i}c^{*}_{j}\int_{-\infty}^{\infty}ds\
[\alpha^{(1)}_{ij}(s)\dot{\alpha}^{(2)}_{ij}(s)-\dot{\alpha}^{(1)}_{ij}(s)\alpha^{(2)}_{ij}(s)],$
(23)

is expressed in terms of the quantities

$\displaystyle\alpha_{ij}^{(a)}(s)=\int\frac{dkdk^{\prime}}{2\pi}\psi_{ai}(k)\psi_{aj}^{*}(k^{\prime})S_{a}(k,k^{\prime})\sqrt{v_{k}v_{k^{\prime}}}e^{i(k-k^{\prime})L_{a}-i(\epsilon_{k}-\epsilon_{k^{\prime}})s}\
\text{PV}\frac{1}{\epsilon_{k}-\epsilon_{k^{\prime}}},$ (24)

where $a=1,2$.

###  4.2 Examples

We analyze the case of massless particles, $m=0$, and ideal detector,
$S(k,k^{\prime})=1$. For two particles prepared in the same state
$\psi_{0}(k)$, that is centered around $x=0$. However, the distances traveled
by the two particles may be different, $L_{1}\neq L_{2}$. We take for
$\psi_{0}$ a Gaussian centered around $k_{0}$,

$\displaystyle\psi_{0}(k)=(2\pi\sigma)^{-1/4}\exp\left[(k-k_{0})^{2}/(4\sigma^{2})\right].$
(25)

Then, we find that $w$ in Eq. (23) equals $Q_{1}[\sigma(L_{1}-L_{2})]$, where

$\displaystyle
Q_{1}(\delta)=\int_{-\infty}^{\infty}dxe^{-2(x-\delta)^{2}}\mbox{erf}\left(\sqrt{2}x\right).$
(26)

The dependence of the function $Q$ on $\delta$ is plotted in Fig. 2. As
expected $w$ vanishes for $L_{1}=L_{2}$ and tends to $\pm\frac{1}{2}$ for
large differences between $L_{1}$ and $L_{2}$, in which case the ordering of
the events is almost certain.

Figure 2:  The function $Q_{1}$ of Eq. (26) (solid) and the function $Q_{2}$
of Eq. (29) (dashed) as a function of $\delta$.

The probabilities in the example above could have been derived from a
classical theory. To see quantum behavior, we consider an superposition state
for the first particle

$\displaystyle\psi(k_{1},k_{2})=\frac{1}{\sqrt{2(1+\nu)}}\psi_{0}(k_{1})\left[1+e^{ik_{1}\ell}\right]\psi_{0}(k_{2}),$
(27)

Here $\ell$ is a path difference for the first particle in one component of
the superposition and $\nu=\int dk_{1}dk_{2}|\psi_{0}(k)|^{2}\cos(k\ell)$. For
the Gaussian (25), $\nu=e^{-\sigma^{2}\ell^{2}/2}\cos(k_{0}\ell)$. For this
initial state

$\displaystyle
w=\frac{Q_{1}(\delta)+2Q_{2}(\delta)\cos\left(\frac{k_{0}}{\sigma}\delta\right)}{2\left[1+e^{-\delta^{2}/2}\cos\left(\frac{k_{0}}{\sigma}\delta\right)\right]},$
(28)

where $\delta=\sigma\ell$, and the function

$\displaystyle
Q_{2}(\delta)=\int_{-\infty}^{\infty}dxe^{-x^{2}-(x-\delta)^{2}}\mbox{erf}\left(\sqrt{2}x\right),$
(29)

is plotted in Fig. 2.

In Fig. (3), we plot $w$ as a function of $\delta$ and of $k_{0}/\sigma$. The
quantum nature of the system is manifested in the oscillatory behavior of the
probabilities.

Figure 3:  The asymmetry $w$ of Eq. (28) as a function of $\delta$ for
constant $k_{0}/\sigma=10$ (left) and as a function of $k_{0}/\sigma$ for
constant $\delta=1$ (right).

##  5 Probabilities for the COoE via a detection model

In the examples of Sec. 4, the measurements of the causal ordering of events
are coarse-grained. This means that the measuring apparatuses record the time
of detection events, and the probabilities for the causal order of events are
obtained by integrating over probabilities with respect to detection times. In
this sense, the construction is formally similar to the one of classical
physics, even if there are no paths at the fundamental level.

However, in quantum theory, it may be possible to define probabilities for
causal ordering of events as fine-grained observables, even if we cannot
distinguish between the times of the individual events. In this section, we
will present a simple model that provides such probabilities for the case of
two potential events, and which can be straightforwardly generalized for $n$
events.

###  5.1 The model

The key idea is to direct a pair of particles towards a detector that can
record only one of them, but not both. As an example, we consider a three-
level system (3LS), with states $|0\rangle,|1\rangle$, and $|2\rangle$.
Suppose that particle 1 can excite only the transition $0\rightarrow 1$, and
particle 2 only the transition $0\rightarrow 2$. If after the interaction of
the particles with the three-level system, we find the system in state
$|1\rangle$, we can surmise that particle 1 was detected first and particle 2
was not detected, and vice versa.

To conform with our definition of an event with a measurement record, we must
place an identical 3LS after the first, in which the particle not absorbed by
the first can be detected. However, this is superfluous for identifying the
COoE in this system, so we will consider the case of a single 3LS.

This set-up is straightforwardly generalized for determining the probabilities
for $n$ events. We require $n$ particles that can be sharply distinguished by
their energies and $n-1$ systems with $n+1$ energy levels, so that the
detection of each particle can be associated to a single transition.

To implement our model, we assume that the incoming particles are described by
a free scalar field $\hat{\phi}(x)$ with mass $m$. The two particles are
distinguished by their initial energies; we can assume that they are prepared
from different sources. The particles interact with one 3LS, which we take to
be located at ${\bf x}=0$.

The total Hamiltonian is a sum of three terms
$\hat{H}_{\phi}+\hat{H}_{3LS}+\hat{H}_{int}$, where $\hat{H}_{\phi}$ is the
field Hamiltonian, expressed in terms of field creation and annihilation
operators,

$\displaystyle\hat{H}_{\phi}=\int d{\bf k}\epsilon_{\bf
k}\hat{a}^{\dagger}_{\bf k}\hat{a}_{bfk},$ (30)

where $d{\bf k}$ stands for $d^{3}k/(2\pi)^{3}$;

$\displaystyle\hat{H}_{3LS}=\Omega_{1}|1\rangle\langle
1|+\Omega_{2}|2\rangle\langle 2|$ (31)

is the 3LS Hamiltonian, and the interaction Hamiltonian

$\displaystyle\hat{H}_{int}=\sum_{a=1}^{2}\lambda_{a}\int\frac{dk}{\sqrt{2\omega_{\bf
k}}}(\hat{a}_{\bf k}\hat{u}_{a+}+\hat{a}^{\dagger}_{\bf k}\hat{u}_{a-}),$ (32)

describes a dipole coupling between the field and the 3LS. Here, $\lambda_{a}$
is the coupling constants associated to the transition $0\rightarrow a$,
$\hat{u}_{a+}=|a\rangle\langle 0|$, and $\hat{u}_{a-}=|0\rangle\langle a|$;
$a=1,2$. This Hamiltonian is a variation of Lee’s Hamiltonian that is commonly
employed in the study of spontaneous decay [27].

###  5.2 Time evolution

To derive the time evolution law for this model, we work in the interaction
picture. Then, the quantum state satisfies the equation

$\displaystyle i\frac{\partial}{\partial
t}|\psi(t)\rangle=\sum_{a=1}^{2}\lambda_{a}\int\frac{d{\bf
k}}{\sqrt{2\omega_{\bf k}}}(\hat{a}_{\bf k}\hat{u}_{a+}e^{-i(\epsilon_{\bf
k}-\Omega_{a})t}+\hat{a}^{\dagger}_{\bf k}\hat{u}_{a-}e^{i(\epsilon_{\bf
k}-\Omega_{a})t})|\psi(t)\rangle.$ (33)

We assume an initial two-particle state for the field and the ground state for
the 3LS. The Hamiltonian employed here causes transitions only to one-
particles states with an excited state for the 3LS. Hence, the state is of the
form

$\displaystyle|\psi(t)\rangle=\int d{\bf k}d{\bf k}^{\prime}c({\bf k},{\bf
k}^{\prime};t)|{\bf k},{\bf k}^{\prime},0\rangle+\sum_{a}\int d{\bf
k}d_{a}({\bf k};t)|k,a\rangle,$ (34)

Substituting into Eq. (33), we obtain

$\displaystyle i\dot{c}({\bf k},{\bf k}^{\prime};t)$ $\displaystyle=$
$\displaystyle\sum_{a}\lambda_{a}\left[\frac{d_{a}({\bf
k};t)}{\sqrt{2\epsilon_{{\bf k}^{\prime}}}}e^{i(\epsilon_{{\bf
k}}-\Omega_{a})t}+\frac{d_{a}({\bf k}^{\prime};t)}{\sqrt{2\epsilon_{\bf
k}}}e^{i(\epsilon_{{\bf k}^{\prime}}-\Omega_{a})t}\right]$ (35) $\displaystyle
i\dot{d}_{a}({\bf k};t)$ $\displaystyle=$ $\displaystyle
2\lambda_{a}\int\frac{d{\bf k}^{\prime}}{\sqrt{2\epsilon_{{\bf
k}^{\prime}}}}c({\bf k},{\bf k}^{\prime};t)e^{-i(\epsilon_{{\bf
k}^{\prime}}-\Omega_{a})t}$ (36)

These equations are to be solved subject to the initial conditions $d_{a}({\bf
k};0)=0$ and $c({\bf k},{\bf k}^{\prime};0)=c_{0}({\bf k},{\bf k}^{\prime})$,
where $c_{0}({\bf k},{\bf k}^{\prime})$ is the initial state of the two
particles. We integrate both sides of Eq. (35) and substitute $c({\bf k},{\bf
k}^{\prime};t)$ to Eq. (36). We obtain

$\displaystyle\dot{d}_{a}({\bf k};t)$ $\displaystyle=$ $\displaystyle
2\lambda_{a}\int\frac{d{\bf k}^{\prime}}{\sqrt{2\epsilon_{{\bf
k}^{\prime}}}}c_{0}({\bf k},{\bf k}^{\prime})e^{-i(\epsilon_{{\bf
k}^{\prime}}-\Omega_{a})t}$ (37) $\displaystyle-$ $\displaystyle
2\lambda_{a}\sum_{b}\lambda_{b}\int\frac{d{\bf k}^{\prime}}{2\epsilon_{{\bf
k}^{\prime}}}e^{-i(\epsilon_{{\bf
k}^{\prime}}-\Omega_{a})t}\int_{0}^{t}dsd_{a}({\bf k};s)e^{i(\epsilon_{{\bf
k}}-\Omega_{b})s}$ $\displaystyle-$ $\displaystyle
2\frac{\lambda_{a}}{\sqrt{\epsilon}_{\bf k}}\sum_{b}\lambda_{b}\int\frac{d{\bf
k}^{\prime}}{2\sqrt{\epsilon_{{\bf k}^{\prime}}}}e^{-i(\epsilon_{{\bf
k}^{\prime}}-\Omega_{a})t}\int_{0}^{t}dsd_{a}({\bf
k}^{\prime};s)e^{i(\epsilon_{{\bf k}^{\prime}}-\Omega_{b})s}.$

Eq. (37) is exact. The term in the second line is proportional to the vacuum
Wightman function $W(t)=\int\frac{d{\bf k}^{\prime}}{2\epsilon_{{\bf
k}^{\prime}}}e^{-i\epsilon_{{\bf k}^{\prime}}t}$, which drops at least with
$e^{-mt}$ for $m\neq 0$ and as $t^{-2}$ for $m=0$. Assuming that the particle
starts sufficiently far from the detector, $d_{a}({\bf k};t)$ becomes
appreciable different from zero at times such that the term proportional to
$W(t)$ is strongly suppressed.

The third-line term in Eq. (37) is of a structure that commonly appears in
elementary treatments of spontaneous decay [31, 27]. It can be calculated by
invoking a version of the Wigner-Weisskopf approximation. For
$\Omega_{a}t>>1$, this expression is strongly dominated by the term with
$b=a$. By carrying out the integration over ${\bf k}^{\prime}$, we obtain

$\displaystyle-\frac{\lambda_{a}^{2}}{2\pi^{2}\sqrt{\epsilon_{\bf
k}}}\int_{m}^{\infty}d\epsilon\frac{(\epsilon^{2}-m^{2})^{3/2}}{\sqrt{\epsilon}}\int_{0}^{t}dse^{-i(\epsilon-\Omega_{a})(t-s)}d_{a}({\bf
k},s).$ (38)

The time integral is negligible except for values of $\epsilon$ around
$\Omega_{a}$. Hence, we are justified in substituting
$(\epsilon^{2}-m^{2})^{3/2}/\sqrt{\epsilon}$ with
$(\Omega_{a}^{2}-m^{2})^{3/2}/\sqrt{\Omega_{a}}$, and then, to extend
integration over $\epsilon$ to $(-\infty,\infty)$. Then, the term (38)
simplifies to $-\frac{1}{2}\eta_{a}\epsilon_{\bf k}^{-1/2}d_{a}({\bf k},t)$,
where

$\displaystyle\eta_{a}=-\frac{\lambda_{a}^{2}(\Omega_{a}^{2}-m^{2})^{3/2}}{\pi\sqrt{\Omega_{a}}}.$
(39)

Eq. (37) becomes

$\displaystyle\dot{d}_{a}({\bf k};t)+\frac{1}{2}\eta_{a}\epsilon_{\bf
k}^{-1/2}d_{a}({\bf k};t)=2\lambda_{a}\int\frac{d{\bf
k}^{\prime}}{\sqrt{2\epsilon_{{\bf k}^{\prime}}}}c_{0}({\bf k},{\bf
k}^{\prime})e^{-i(\epsilon_{{\bf k}^{\prime}}-\Omega_{a})t}.$ (40)

This is a linear inhomogenous equation of first order. The Green function for
the corresponding homogeneous equation is simply
$\theta(t-t^{\prime})e^{-\frac{1}{2}\eta_{a}\epsilon_{\bf
k}^{-1/2}(t-t^{\prime})}$. Hence, we obtain

$\displaystyle d_{a}({\bf k};t)=2\lambda_{a}\int\frac{d{\bf
k}^{\prime}}{\sqrt{2\epsilon_{{\bf k}^{\prime}}}}c_{0}({\bf k},{\bf
k}^{\prime})\int_{0}^{t}dse^{-\frac{1}{2}\eta_{a}\epsilon_{\bf
k}^{-1/2}(t-s)}e^{-i(\epsilon_{{\bf
k}^{\prime}}-\Omega_{a})s}=2\lambda_{a}\int d{\bf k}^{\prime}c_{0}({\bf
k},{\bf k}^{\prime})h_{a}(\epsilon_{\bf k},\epsilon_{{\bf k}^{\prime}};t),$

where

$\displaystyle
h_{a}(\epsilon,\epsilon^{\prime};t)=\frac{e^{-\frac{1}{2}\eta_{a}\epsilon^{-1/2}t}-e^{-i(\epsilon^{\prime}-\Omega_{a})t}}{\sqrt{2\epsilon^{\prime}}\left[\frac{1}{2}\eta_{a}\epsilon^{-1/2}-i(\epsilon^{\prime}-\Omega_{a})\right]}.$
(41)

The detection probability is non-negligible only if ${\bf k}$ is along the
axis that connects the source to the detector. Hence, the problem is
effectively one-dimensional. Therefore, we can substitute the initial state
with $c_{0}(k,k^{\prime})$, where $k,k^{\prime}>0$, and write
$d_{a}(k;t)=2\lambda_{a}\int\frac{dk^{\prime}}{2\pi}c_{0}(k,k^{\prime})h_{a}(\epsilon_{k},\epsilon_{k^{\prime}};t)$.

###  5.3 An example

Consider an initial state

$\displaystyle
c_{0}(k,k^{\prime})=\frac{1}{\sqrt{2}}\left[\psi_{1}(k)\psi_{2}(k^{\prime})+\psi_{1}(k^{\prime})\psi_{2}(k)\right],$
(42)

where $\psi_{i}$, for $i=1,2$, is centered around momentum $k_{i}$, or,
equivalently, on energy $\epsilon_{i}=\sqrt{k_{i}^{2}+m^{2}}$. We assume that
there is no overlap between $\psi_{1}$ and $\psi_{2}$. Then, we can
approximate

$\displaystyle
d_{a}(k;t)=\lambda_{a}\left[\psi_{1}(k)F_{2a}(t)+\psi_{2}(k)F_{1a}(t)\right],$
(43)

where

$\displaystyle
F_{ia}(t)=\int\frac{dk}{2\pi}\psi_{i}(k)\frac{e^{-\Gamma_{ia}t}-e^{-i(\epsilon_{k}-\Omega_{a})t}}{\sqrt{2\epsilon_{k}}\left[\Gamma_{ia}-i(\epsilon_{k}-\Omega_{a})\right]},$
(44)

where $\Gamma_{1a}=\frac{1}{2}\eta_{a}\epsilon_{2}^{-1/2}$ and
$\Gamma_{2a}=\frac{1}{2}\eta_{a}\epsilon_{1}^{-1/2}$. Then, the probability
$p_{a}(t)$ that the 3LS is found in an excited state is given by

$\displaystyle p_{a}(t)=\int
dk|d_{a}(k;t)|^{2}=\lambda_{a}^{2}\left(|F_{1a}(t)|^{2}+|F_{2a}(t)|^{2}\right).$
(46)

Let the states $\psi_{i}(k)$ be well localized around $x=-L$, so that they can
be written as $\chi(k-k_{i})e^{ikL}$, where $\chi$ is a positive function
peaked around $k=0$, for example, a Gaussian. Then, the typical behavior of
$|F_{ia}(t)|^{2}$ is given in Fig. 4. The function is negligible prior to the
arrival time $t_{a}=mL/k_{i}$ of the particle to the locus of 3LS. Then, it
jumps to a finite value, which then decays with a rate given by $\Gamma_{ia}$.

Figure 4: Typical plot of the functions $|F_{ia}(t)|^{2}$ as a function of
time $t$, for a Gaussian function $\chi(k)$. The function jumps to a finite
value when the particle arrives at the 3LS, and then it decays with a rate of
$\Gamma_{ia}$.

The peak value of $|F_{ia}(t)|^{2}$ is approximately proportional to the
Breit-Wigner term $[(\Gamma_{ia}^{2}+(\epsilon_{i}-\Omega_{a})^{2}]^{-1}$.
Supposing that we choose $\epsilon_{a}\simeq\Omega_{a}$, and that
$\Gamma_{ia}<<|\Omega_{1}-\Omega_{2}|$, for all $i,a=1,2$, then, the terms
$|F_{11}|^{2}$ and $|F_{22}|^{2}$ dominate in the probability assignment, and

$\displaystyle p_{a}(t)=\lambda_{a}^{2}|F_{aa}(t)|^{2}.$ (47)

The behavior of the probabilities is characteristic of resonant fluorescence.
The 3LS absorbs one of the two particles, and after a time of order
$\Gamma_{aa}^{-1}$ it re-emits the particle, albeit in a different direction.
Hence, the energy of the fluorescent particle determines whether the ordering
$M_{1}$ or $M_{2}$ was realized.

##  6 Conclusions

We provided a general definition of events in quantum theory, and showed how
to construct probabilities associated to the causal ordering such events. Our
notion of events is very different from that of Refs. [3, 4], and it is
naturally related to the relativistic notion of events. Our analysis clarifies
that the existence of an indefinite quantum causal order of events has no
relation to quantum gravity, as this causal order is a dynamical consequence
of the quantum nature of the matter degrees of freedom. The COoE should not be
conflated with the causal structure of spacetime, which we take to be fixed
and unchanged in absence of gravity. Further work is needed, in order to
explore how the quantum probabilities for causal order defined here differ
from the corresponding classical one, for example, if they violate Bell’ type
inequalities.

The model systems considered in this paper are experimentally accessible. The
set-ups considered in Sec. 5 are essentially quantum races, that is the causal
order of events coincides with the order that a number of distinguishable
particles arrive in a specific finish line. The set-up of Sec. 6, when applied
to photons, involves a variation of resonant fluorescence with specially
engineered multi-level atoms that play the role of detectors for the causal
ordering that is being realized.

## References

  * [1] C. Anastopoulos and N. Savvidou, Time-of-Arrival Probabilities for General Particle Detectors, Phys. Rev. A86, 012111 (2012).
  * [2] C. Anastopoulos, B.L. Hu and K. Savvidou, Quantum field theory based quantum information: Measurements and correlations, Ann. Phys. 450, 169239 (2023).
  * [3] O. Oreshkov, F. Costa, and C. Brukner, Quantum correlations with no causal order, Nat. Commun. 3, 1092 (2012).
  * [4] J. Chiribella, G. M. D’Ariano, P. Perinotti, and B. Valiron, Quantum computations without definite causal structure, Phys. Rev. A 88, 022318 (2013).
  * [5] M. Araújo, F. Costa, and C. Brukner, Computational Advantage from Quantum-Controlled Ordering of Gates, Phys. Rev. Lett. 113, 250402 (2014).
  * [6] M. Araújo, P. A. Guérin, and Ä. Baumeler, Quantum computation with indefinite causal structures, Phys. Rev. A96, 052315 (2017).
  * [7] P. A. Guérin, A. Feix, M. Araújo, and C. Brukner, Exponential communication complexity advantage from quantum superposition of the direction of communication Phys. Rev. Lett. 117, 100502 (2016).
  * [8] M. M. Taddei, J. Cariñe, D. Martínez, T. García, N. Guerrero, A. A. Abbott, M. Araújo, C. Branciard, E. S. Gómez, S. P. Walborn, L. Aolita, and G. Lima, Computational advantage from the quantum superposi- tion of multiple temporal orders of photonic gates, PRX Quantum 2, 0100320 (2021).
  * [9] G. Rubino, L. A. Rozema, A. Feix, M. Araújo, J. M. Zeuner, L. M. Procopio, C. Brukner, and P. Walther, Experimental verification of an indefinite causal order, Sci. Adv. 3, e1602589 (2017).
  * [10] K. Goswami, C. Giarmatzi, M. Kewming, F. Costa, C. Branciard, J. Romero, and A. White, Indefinite causal order in a quantum switch, Phys. Rev. Lett. 121, 090503 (2018)
  * [11] K. Goswami and J. Romero, Experiments on quantum causality, AVS Quantum Sci. 2, 037101 (2020).
  * [12] V. Vilasini and R. Colbeck, A general framework for cyclic and fine-tuned causal models and their compatibility with space-time, Phys. Rev. A106, 032204 (2022).
  * [13] V. Vilasini and R. Renner, Embedding cyclic causal structures in acyclic spacetimes: no-go results for process matrices, arXiv:2203.11245.
  * [14] A. Einstein, Die Grundlage der allgemeinen Relativitätstheorie, Annalen der Physik 49, 769–822 (1916).
  * [15] P. G. Bergmann, Observables in General Relativity, Rev. Mod. Phys. 33, 510 (1961).
  * [16] N. D. Mermin, It’s about time: Understanding Einstein’s Relativity (Princeton, Princeton University Press, 2005).
  * [17] A. Peres, Quantum Theory: Concepts and Method (Kluwer, 2002).
  * [18] W. Pauli, in Handbuch der Physik, Vol. 5, Part 1, edited by S. Flügge, S. 1–168 (Springer, Berlin, 1958).
  * [19] C. Anastopoulos and N. Savvidou, Time-of-Arrival Correlations, Phys. Rev. A95, 032105 (2017).
  * [20] C. Anastopoulos and N. Savvidou, Time of arrival and Localization of Relativistic Particles, J. Math. Phys. 60, 0323301 (2019).
  * [21] C. Anastopoulos and N. Savvidou, Coherences of Accelerated Detectors and the Local Character of the Unruh Effect, J. Math. Phys. 53, 012107 (2012).
  * [22] C. Anastopoulos, Time-of-arrival probabilities and quantum measurements. III. Decay of unstable states, J. Math. Phys. 49, 022103 (2008).
  * [23] K. Savvidou, The Action Operator for Continuous-time Histories J. Math. Phys. 40, 5657 (1999); Continuous Time in Consistent Histories, gr-qc/9912076.
  * [24] N. Savvidou, Space-time Symmetries in Histories Canonical Gravity, in ”Approaches to Quantum Gravity”, edited by D. Oriti (Cambridge University Press, Cambridge 2009).
  * [25] M. Gell-Mann and J. B. Hartle, Classical Equations for Quantum Systems, Phys. Rev. D47, 3345 (1993).
  * [26] R. Omnés, The Interpretation of Quantum Mechanics, (Princeton University Press, Princeton 1994).
  * [27] C. Anastopoulos,Quantum Theory: A Foundational Approach (Cambridge: Cambridge University Press, 2023).
  * [28] C. J. Isham, Canonical Quantum Gravity and the Problem of Time, Lecture Notes at Salamanca Summer School (1992), gr-qc/9210011.
  * [29] K. Kuchar, Time and interpretations of quantum gravity ,in ‘Proceedings of the 4th Canadian Conference on General Relativity and Relativistic Astrophysics’,World Scientific, Singapore 1992.
  * [30] E. Anderson, The Problem of Time: Quantum Mechanics Versus General Relativity (Springer, 2017)
  * [31] M. O. Scully and M. S. Zubairy, Quantum Optics (Cambridge : Cambridge University Press 2012).


# First principles feasibility assessment of a topological insulator at the
InAs/GaSb interface

Shuyang Yang  Department of Materials Science and Engineering, Carnegie Mellon
University, Pittsburgh, PA 15213, USA  Derek Dardzinski  Department of
Materials Science and Engineering, Carnegie Mellon University, Pittsburgh, PA
15213, USA  Andrea Hwang  Department of Materials Science and Engineering,
Carnegie Mellon University, Pittsburgh, PA 15213, USA  Dmitry I. Pikulin
Microsoft Quantum, Redmond, WA 98052, USA  Microsoft Quantum, Microsoft
Station Q, University of California, Santa Barbara, California 93106-6105, USA
Georg W. Winkler  Microsoft Quantum, Microsoft Station Q, University of
California, Santa Barbara, California 93106-6105, USA  Noa Marom
<EMAIL_ADDRESS>Department of Materials Science and Engineering,
Carnegie Mellon University, Pittsburgh, PA 15213, USA  Department of
Chemistry, Carnegie Mellon University, Pittsburgh, PA 15213, USA  Department
of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA

###### Abstract

First principles simulations are conducted to shed light on the question of
whether a two-dimensional topological insulator (2DTI) phase may be obtained
at the interface between InAs and GaSb. To this end, the InAs/GaSb interface
is compared and contrasted with the HgTe/CdTe interface. Density functional
theory (DFT) simulations of these interfaces are performed using a machine-
learned Hubbard U correction [npj Comput. Mater. 6, 180 (2020)]. For the
HgTe/CdTe interface our simulations show that band crossing is achieved and an
inverted gap is obtained at a critical thickness of 5.1 nm of HgTe, in
agreement with experiment and previous DFT calculations. In contrast, for
InAs/GaSb the gap narrows with increasing thickness of InAs; however the gap
does not close for interfaces with up to 50 layers (about 15 nm) of each
material. When an external electric field is applied across the InAs/GaSb
interface, the GaSb-derived valence band maximum is shifted up in energy with
respect to the InAs-derived conduction band minimum until eventually the bands
cross and an inverted gap opens. Our results show that it may be possible to
reach the topological regime at the InAs/GaSb interface under the right
conditions. However, it may be challenging to realize these conditions
experimentally, which explains the difficulty of experimentally demonstrating
an inverted gap in InAs/GaSb.

††preprint: APS/123-QED

##  I INTRODUCTION

Two-dimensional topological insulators (2DTIs) have attracted increasing
attention in recent years owing to the emergence of helical edge states and
backscattering-free edge currents relevant for applications in spintronics and
quantum computing [1, 2, 3]. 2DTIs were first proposed based on a theoretical
model of graphene incorporating spin-orbit interactions [4]. However, the
required type of spin-orbit coupling in graphene is too weak to observe the
quantum spin Hall effect (QSHE) experimentally [5]. Later, a proposal for 2DTI
was made based on a HgTe/CdTe quantum well (QW) [6] and the signatures of the
QSHE were experimentally demonstrated [7]. When the thickness of the HgTe in
the QWs is varied, the band structure changes from a trivial insulator to a
2DTI with an inverted gap when a critical thickness is reached [6, 7, 8, 9].
2DTIs have been proposed in additional materials systems, some of which have
shown promising signs [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22].

In the present work we focus on another QW structure, InAs/GaSb. It has been
proposed that a 2DTI may be realized in InAs/GaSb QWs because the band lineup
of coupled InAs/GaSb QWs could lead to the coexistence of electrons and holes
at the charge neutrality point [23, 24]. The topological insulator phase would
arise if the band ordering were inverted and coupling between electron and
hole states opened a hybridization gap which is necessarily topological due to
the orbital structure of the hybridized bands [25]. Such band ordering could
potentially be achieved by choosing appropriate QW thickness and by applying
an external electric field [5]. InAs/GaSb QWs are in the family of well-
studied III-V compounds and have thus attracted considerable experimental
interest [26, 27, 28, 29, 30, 31, 32, 33, 34]. The experiments have provided
some encouraging signs of edge conductance in the material. However, a phase
diagram showing a clear topological transition accompanied by edge state
formation has yet to be demonstrated. Here, we use first principles
simulations to investigate whether it would be possible to realize a 2DTI at
the InAs/GaSb interface and under what conditions.

The HgTe/CdTe and InAs/GaSb interfaces have been studied theoretically using a
variety of methods. This includes the k$\cdot$p method [35, 9, 36, 37, 38,
39], pseudopotential models [40, 41], and tight-binding [42, 43]. The drawback
of these semi-empirical methods is that the fitting to experimental data
largely determines the extent of their predictive capability. Atomistic ab
initio simulations may provide a more accurate representation of the
electronic properties and their dependence on the structure of the interface.
First principles studies based on density functional theory (DFT) have
investigated the influence of thickness on the edge states of HgTe/CdTe(100)
[44, 45]. Using different exchange-correlation functionals and different
thicknesses of CdTe, Ref. [44] predicted a critical thickness of 4.6 nm of
HgTe, whilst Ref. [45] predicted a critical thickness 6.5 nm of HgTe. Both
results are close to the experimental critical thickness of 6.3 nm [6].

For InAs and GaSb, local and semi-local exchange-correlation functionals
severely underestimate the band gaps to the point that they reduce to zero
[46], due to the self-interaction error (SIE). Some DFT studies of InAs/GaSb
have applied an empirical correction to the DFT band gaps [47, 48]. Others
have used hybrid functionals, which mitigate the effect of SIE by including a
fraction of exact exchange [49]. An alternative approach, which has been used
to obtain more accurate band gaps for InAs/GaSb is many-body perturbation
theory within the $GW$ approximation, where $G$ stands for the one-particle
Green’s function and $W$ stand for the screened Coulomb interaction [50].
Although hybrid DFT functionals and the $GW$ approximation produce
significantly improved band gaps, their high computational cost limits their
applicability to relatively small system sizes. Therefore, these methods have
been used only for periodic heterostructures of InAs/GaSb with very few layers
[49, 50]. DFT studies of large interface slab models with vacuum regions have
not been conducted. All previous ab initio studies of InAs/GaSb have not
reported the band structure and band alignment at the interface and have not
shown an inverted band gap. Furthermore, previous studies have not considered
the effect of applying an electric field, which plays an important role in
experiments, and therefore should be considered computationally.

Recently, we have introduced a new method of DFT with a machine-learned
Hubbard U correction, which can provide a solution for accurate and efficient
simulations of InAs and GaSb [51]. Within the Dudarev formulation of DFT+U
[52] the effective Hubbard U is defined as $U_{eff}=U-J$, where $U$ represents
the on-site Coulomb repulsion, and $J$ represents the exchange interaction.
For a given material, the $U_{eff}$ parameters of each element are machine-
learned using Bayesian optimization (BO). The BO algorithm finds the optimal
$U_{eff}$ values that maximize an objective function formulated to reproduce
as closely as possible the band gap and the qualitative features of the band
structure obtained with a hybrid functional. The DFT+U(BO) method allows for
negative $U_{eff}$ values. Negative $U_{eff}$ values are theoretically
permissible when the exchange term, $J$, is larger than the on-site Coulomb
repulsion, $U$ [53, 54, 55, 56, 57]. We have found that negative $U_{eff}$
values are necessary to produce band gaps for narrow-gap semiconductors, such
as InAs and GaSb. Because the reference hybrid functional calculation is
performed only once for the bulk material to determine the optimal $U_{eff}$
values, the computational cost of DFT+U(BO) calculations for interfaces is
comparable to semi-local DFT.

In this work, we use the DFT+U(BO) method to study the HgTe/CdTe and InAs/GaSb
interfaces. For the HgTe/CdTe interface, we obtain band crossing at a critical
thickness of 5.1 nm of HgTe, and subsequently an inverted gap is observed. Our
results are in agreement with experiment and previous DFT studies, thus
validating the DFT+U(BO) method. For the InAs/GaSb interface, we find that
increasing the thickness of InAs leads to gap narrowing. However, band
crossing is not obtained up to the largest number of layers calculated here.
When an external electric field is applied across the InAs/GaSb interface, the
GaSb-derived valence band maximum is shifted up in energy compared to the
InAs-derived conduction band minimum. Band crossing is achieved at a critical
field, followed by an inverted gap which widens and shifts higher above the
Fermi level as the field is increased. Our results indicate that it may be
possible to reach the topological regime in InAs/GaSb QWs. However, doing so
would require a combination of careful interface engineering, a considerable
electric field across the interface, and gating to tune the position of the
Fermi level. This explains the difficulty of experimentally demonstrating an
inverted gap in InAs/GaSb.

##  II METHODS

###  II.1 Computational details

DFT calculations were performed using the Vienna ab initio simulation package
(VASP) [58] with the projector augmented wave method (PAW) [59, 60]. The
generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE)
[61, 62] was used with a Hubbard U correction [52] determined by Bayesian
optimization [51], as detailed below. Spin-orbit coupling (SOC) [63] was
included throughout and the energy cutoff was set to 400 eV. For bulk band
structure calculations a 8$\times$8$\times$8 k-point grid was used to sample
the Brillouin zone. For interface calculations a 8$\times$8$\times$1 k-point
grid was used to sample the interface Brillouin zone and dipole corrections
[64] were included. Bulk unfolding [65] was applied to project interface band
structures onto the primitive cell, as described in the SI.

###  II.2 Performance of PBE+U(BO)

Figure 1:  Performance of different DFT functionals for CdTe and HgTe: PBE
band structures of (a) CdTe and (b) HgTe; HSE band structures of (c) CdTe and
(d) HgTe; PBE+U(BO) band structures of (e) CdTe and (f) HgTe; The
contributions of the Cd/Hg $s$, Cd/Hg $d$, and Te $p$ states are indicated by
the red, green, and yellow dots, respectively. $A$ is the point along
$X-\Gamma$ in the bulk’s Brillouin zone with the coordinates (0.1, 0.1, 0).
$A-\Gamma-A$ is mapped to $\overline{A}-\overline{\Gamma}-\overline{A}$ in the
(001) direction.

The PBE functional fails to provide an adequate description of the band
structures of the materials studied here. The cases of InAs and GaSb have been
discussed in detail in [51]. For CdTe, Fig. 1a shows that PBE severely
underestimates the bad gap compared to the experimental value of 1.60 eV [66].
This is because the Cd 4$d$ states, which contribute significantly to the top
of the valence band, are pushed up in energy due to the SIE [67]. For HgTe,
Fig. 1b shows that PBE produces an incorrect band shape and band ordering at
the $\Gamma$ point with the Hg $s$ orbitals and Te $p$ orbitals inverted
around 1 eV below the Fermi level [45]. These issues are rectified by the
Heyd-Scuseria-Ernzerhof (HSE) [68, 69] hybrid functional, as shown in Fig. 1c
for CdTe and Fig. 1d for HgTe. However, the computational cost of HSE is too
high for simulations of large interface models.

To achieve a balance between accuracy and efficiency, a Hubbard $U$ correction
was applied to the $p$ orbitals of In, As, Ga and Sb and the $d$ orbitals of
Hg and Cd within the Dudarev approach [52]. For each orbital, the optimal
value of Ueff was machine learned by Bayesian optimization [51]. The objective
function was formulated to reproduce as closely as possible the band structure
produced by HSE:

${f(\vec{U})}=-{\alpha_{1}}(\text{E}_{\text{g}}^{\text{HSE}}-\text{E}_{\text{g}}^{\text{PBE+U}})^{2}-{\alpha_{2}}(\Delta\text{Band})^{2}$
(1)

Here, $\vec{U}$ = [$U^{1}$, $U^{2}$,…,$U^{n}$] is the vector of $U_{eff}$
values applied to different atomic species and $U^{i}\in[-10,10]$ eV.
$\Delta\text{Band}$ is defined as the mean squared error of the PBE+U band
structure with respect to HSE:

${\Delta\text{Band}}=\sqrt{\frac{1}{N_{E}}\sum^{N_{k}}_{i=1}\sum^{N_{b}}_{j=1}(\epsilon_{HSE}^{j}[k_{i}]-\epsilon_{PBE+U}^{j}[k_{i}])^{2}}$
(2)

$N_{E}$ represents the total number of eigenvalues, $\epsilon$, included in
the comparison, $N_{k}$ is the number of $k$-points, and $N_{b}$ is the number
of bands selected for comparison. To avoid double counting the band gap
difference in the calculation of $\Delta\text{Band}$, the valence band maximum
(VBM) and conduction band minimum (CBM) are shifted to zero for both the PBE+U
and HSE band structures. Hence, $\Delta\text{Band}$ captures differences in
the qualitative features of the band structures produced by PBE+U vs. HSE,
independently of the difference in the band gap. The coefficients $\alpha_{1}$
and $\alpha_{2}$ may be used to assign different weights to the band gap vs.
the band structure. The default values are 0.25 and 0.75, respectively. For
CdTe, we set $\alpha_{1}=\alpha_{2}=0.5$ to assign a higher weight to the band
gap term. For HgTe, which is a metal, we set $\alpha_{1}=0$ and
$\alpha_{2}=1$.

For InAs and GaSb the optimal values of Ueff have been found to be:
U${}_{eff}^{In,p}$ = -0.5 eV, U${}_{eff}^{As,p}$ = -7.5 eV, U${}_{eff}^{Ga,p}$
= 0.8 eV, U${}_{eff}^{Sb,p}$ = -6.9 eV, as reported in [51]. With these
parameters, DFT+U(BO) yields a band gap of 0.31 eV for InAs, in good agreement
with the experimental value of 0.41 eV [70], and a band gap of 0.45 eV for
GaSb, which is somewhat underestimated compared to the experimental value of
0.81 eV [70]. For CdTe, BO produces an optimal value of U${}_{eff}^{Cd,d}$ =
8.3 eV, somewhat higher than the value of 7 eV used in [67]. This results in a
band gap of 0.87 eV, which is closer to experiment than previous ab initio
calculations [71, 72]. The qualitative features of the PBE+U(BO) band
structure are in agreement with HSE, as shown in Fig. 1e, however the gap and
the band width are still somewhat underestimated. For HgTe, BO produces a
value of U${}_{eff}^{Hg,d}$ = 8.4 eV, somewhat lower than the value of 9.4 eV
used in Ref. [45]. The band structure, shown in Fig. 1d, has the correct band
shape, comparable to the HSE band structure, and is in agreement with Ref.
[45]. To demonstrate the transferability of the U${}_{e}ff$ values obtained by
BO from bulk materials to interfaces, we compare the band structures produced
by PBE+(BO) and HSE for an InAs/GaSb interface with 5 layers of InAs and 5
layers of GaSb, constructed as detailed below. Fig 2 shows that overall good
agreement is obtained between PBE+U(BO) and HSE, however PBE+U(BO) somewhat
underestimates the band gap and the band width. We note that the Ueff values
obtained here are based on the implementation of the Dudarev formalism in
VASP. Different DFT+U implementations may yield different results [73, 74].

Figure 2:  The band structure of an InAs/GaSb interface with 5 layers of InAs
and 5 layers of GaSb obtained with (a) HSE and (b) DFT+U(BO). Orange and green
dots indicate the contributions of InAs and GaSb, respectively.

###  II.3 Interface model construction

For the HgTe/CdTe(100) interface, we constructed periodic heterostructures,
similar to Ref. [45]. However, we used a larger number of CdTe layers to
ensure convergence, as detailed below. The thickness of HgTe was varied to
study the evolution of the electronic structure. The experimental lattice
constants of 6.45 $\mathrm{\SIUnitSymbolAngstrom}$ for HgTe and 6.48
$\mathrm{\SIUnitSymbolAngstrom}$ for CdTe are closely matched [75]. We assumed
that an epitaxially matched HgTe film would grow on top of a CdTe substrate
with the experimental lattice constant of 6.48
$\mathrm{\SIUnitSymbolAngstrom}$.

For the InAs/GaSb interface, we constructed two types of interface slab
models: The InSb-type interface has In and Sb as the terminal atoms at the
surfaces and interface. The GaAs-type interface has Ga and As as the terminal
atoms. The experimental lattice constants of 6.058
$\mathrm{\SIUnitSymbolAngstrom}$ for InAs and 6.096
$\mathrm{\SIUnitSymbolAngstrom}$ for GaSb [70] are closely matched. We assumed
that an epitaxially matched InAs film would grow on top of GaSb with the
lattice constant of 6.096 $\mathrm{\SIUnitSymbolAngstrom}$, based on the
experiment in Ref. [34]. To study the effect of the InAs and GaSb thickness,
interface models were constructed with the number of layers of each material
varying from 10 to 50. The notation ”A/B” is used to describe an InAs/GaSb
interface with A layers of InAs and B layers of GaSb. A vacuum region of about
40 $\mathrm{\SIUnitSymbolAngstrom}$ was added to the interface model to
prevent spurious interactions between periodic images (for the purpose of band
unfolding the closest integer number of primitive cells to 40
$\mathrm{\SIUnitSymbolAngstrom}$ was used [65]). In order to terminate
dangling bonds, In and Ga atoms on the surface were passivated by pseudo
hydrogen atoms with 1.25 fractional electrons, whereas As and Sb atoms on the
surface were passivated by pseudo hydrogen atoms with 0.75 fractional
electrons. Structural relaxation was performed for the surface atoms and
passivating pseudo-hydrogen atoms until the change of the all forces was below
10-3 eV/$\mathrm{\SIUnitSymbolAngstrom}$.

The number of layers included in slab models needs to be converged to the bulk
limit to avoid quantum size effects. For semiconductors the band gap is
typically used as a the convergence criterion [76, 65]. Fig. 3 shows the band
gap as a function of the number of layers for InAs(100), GaSb(100), and
CdTe(100). We note that here ”layer” is defined as one atomic layer. In each
iteration, the number of layers was increased by 8 for InAs and GaSb and by 6
for CdTe. If the band gap difference between the current iteration and the
previous iteration was within $1\times 10^{-2}$ eV, the current number of
layers was regarded as converged. For InAs and GaSb surfaces, 50 layers are
required, whereas for CdTe 40 layers are required to converge the band gap.
The converged band gap values are close to the bulk values. The size of the
interface models used to simulate the effect of an electric field was limited
to 10 layers of InAs with 10 layers of GaSb due to convergence issues, as
detailed in the SI.

Figure 3:  The band gap obtained with PBE+U(BO) as a function of the number of
layers for InAs(100), GaSb(100), and CdTe(100) surface slabs.

##  III RESULTS AND DISCUSSION

###  III.1 HgTe/CdTe

To validate the DFT+U(BO) method, we begin by applying it to the well-studied
HgTe/CdTe interface. Bulk-unfolded band structures of HgTe/CdTe
heterostructures with 40 layers of CdTe and a varying number of HgTe layers
are shown in Fig. 4. The red dots indicate the contributions from Hg $s$
orbitals and the blue dots indicate the contributions from Te $p$ orbitals.
The band gap value as a function of the number of HgTe layers is shown in Fig.
5. Negative values indicate an inverted band gap. A drastic change is observed
with the thickness of HgTe. When the number of layers is below 16, the
interface behaves as a trivial insulator, with the Hg $s$ orbitals forming the
bottom of the conduction band and the Te $p$ orbitals forming the top of the
valence band. When the number of HgTe layers reaches 16, a transition point
from a trivial insulator to a topological insulator occurs. At this transition
point, both the CBM and VBM show a hybridized $sp$ character. When the number
of HgTe layers exceeds 16, an inverted gap opens, leading to the occurrence of
a topologically nontrivial phase, in which the VBM is dominated by Hg $s$
states and the CBM is dominated by Te $p$ states. The critical thickness of 16
layers, corresponds to 5.1 nm in good agreement with the experimental result
of 6.3 nm (around 19 layers) [6]. Our result is comparable to previous DFT
calculations, which used different functionals and considered structures with
fewer layers of CdTe. Ref. [44] obtained a critical thickness of 4.6 nm of
HgTe on top of 4 layers of CdTe using the modified Becke-Johnson (MBJ)
functional.Ref. [45] obtained a critical thickness of 6.5 nm of HgTe on top of
10 layers of CdTe using GGA+U for HgTe and GGA for CdTe. Thus, the DFT+U(BO)
method successfully describes the electronic structure of the HgTe/CdTe
interface and captures the transition from trivial to topological behavior.

Figure 4:  Band structures of a HgTe/CdTe interface with 40 layers of CdTe and
(a) 4 layers, (b) 16 layers, and (c) 20 layers of HgTe. The red dots indicate
the contributions of Hg $s$ states and the blue dots indicate the
contributions of the Te $p$ states. Figure 5:  The band gap of a HgTe/CdTe
interface with 40 layers of CdTe as a function of the number of HgTe layers.
Negative values indicate an inverted band gap.

###  III.2 InAs/GaSb

####  III.2.1 Effect of layer thickness

Figure 6:  Band gap values as a function of number of layers for 50-layer
InAs/X-layer GaSb and X-layer InAs/50-layer GaSb of InSb-type and GaAs-type
interface.

To investigate the influence of the thickness of InAs and GaSb on the band
gap, we conducted two series of calculations for InSb-type and GaAs-type
interfaces. In one series, the thickness of InAs was fixed at 50 layers InAs
and the number of GaSb layers (X) was varied. In the other series, the
thickness of GaSb was fixed at 50 layers and the number of InAs layers (X) was
varied. The results are shown in Fig. 6. For the InAs(50)/GaSb(X) series, the
band gap of the InSb-type interface increases with increasing GaSb thickness,
whereas the band gap of the GaAs-type interface does not change significantly.
For the InAs(X)/GaSb(50) series, the band gap decreases with increasing InAs
thickness for both interface types, although the gap of the InSb-type
interface remains smaller than that of the GaAs-type interface throughout. The
trend of the gap decreasing with the increase in InAs thickness is in
agreement with experimental observations [77]. The thickest interface we were
able to calculate comprises 50 layers, which corresponds to about 15 nm of
each material. The bulk-unfolded band structure of a 50/50 InSb-type interface
is shown in the SI. Because this interface still has a gap of over 0.2 eV, and
the rate of the gap narrowing decreases with increasing InAs thickness, as
shown in Fig. 6, we estimate that it would either require a significantly
thicker film of InAs for the gap to completely close or the gap would approach
a finite asymptotic limit rather than close. In addition to increasing the QW
thickness, strain engineering, which is not taken into account here, may also
help modulate the gap. [78, 79, 80].

We note that an analysis based on the empirical 8-band Kane model found band
inversion and the quantum spin Hall phase for an InAs thickness above 9 nm at
fixed 10 nm GaSb thickness [25]. However, this analysis was based on empirical
parameters for the material and interface properties and did not take the
atomic details of the interface structure into account. For example, in Ref.
[25] the band alignment at the interface was chosen such that the GaSb valence
band is 150 meV higher than the InAs conduction band leading to a band
inversion even for relatively thin layers. In contrast, within our first
principles approach, we find that band inversion is not achieved up to an InAs
thickness of 15 nm for a range of GaSb thicknesses including 10 nm.
Furthermore, we found that the atomic details of the interface, like the type
of bonds formed at the interface (InSb or GaAs), are relevant, which was
neglected in the effective theory of Ref. [25]. Finally, it should be noted
that experiments seem to indicate that an electric field is required to
achieve an inverted regime in InAs/GaSb heterostructures [81].

####  III.2.2 Effect of electric field

Figure 7:  Electronic structure of a 10/10 InSb-type interface with different
external electric fields. a-c) Bulk unfolded band structures with the
contributions of the interface layers of InAs and GaSb colored in orange and
green, respectively.

The band alignment at the interface of InAs/GaSb can be manipulated by
applying external gate voltages. Ref.[81] has presented strong experimental
evidence that the gap closes when the external gate voltages reach a critical
value. Therefore, we performed DFT simulations for interface slabs in presence
of electric field. In the VASP code, an external electric field is simulated
by adding an artificial dipole sheet in the vacuum region of the unit cell
[64]. Due to screening effects and the electric susceptibility inside the
materials, the effective electric field at the interface may be significantly
smaller than the input electric field [64, 82]. To estimate the effective
electric field, we calculated the gradient of the potential in the InAs and in
GaSb, based on the electrostatic potential averaged over xy plane. The
averaged gradient is taken as the effective electric field. The full account
of the effective field estimation is provided in the SI. The electric field is
applied perpendicular to the plane of the interface and points from the GaSb
side to the InAs side. We note that in VASP only an external electric field
can be set, whereas in experiments the position of the Fermi level can be
independently controlled by applying front-gate and back-gate voltages. Owing
to convergence issues in DFT calculations with external electric fields (see
SI), the largest interfaces we were able to calculate comprise 10 layers of
InAs and 10 layers of GaSb.

Fig. 7 shows the band structure of a 10/10 InSb-type interface. When no
electric field is applied, the interface is in the trivial insulator state.
The CBM is dominated by the interface InAs layer (orange), whereas the
interface GaSb layer (green) contributes predominantly to the VBM. As the
electric field increases, the bands contributed by the GaSb shift upwards with
respect to the bands contributed by the InAs and the gap narrows. When the
input electric field reaches 0.25 V/$\mathrm{\SIUnitSymbolAngstrom}$, which
corresponds to an effective field of 0.014 V/$\mathrm{\SIUnitSymbolAngstrom}$,
the GaSb VBM overlaps with the InAs CBM, the gap closes, and band crossing
occurs. Our results are qualitatively in agreement with previous studies [83,
25, 81], which indicated that the band gap in InAs/GaSb could be closed via an
external electric field. When the electric field is increased further, an
inverted gap opens. As the electric field is increased, the inverted gap
expands, but also shifts higher above the Fermi level. Fig 8 shows the
position of the inverted gap above the Fermi level at the $\Gamma$ point as a
function of the electric field. With an input electric field of 0.35
V/$\mathrm{\SIUnitSymbolAngstrom}$, which corresponds to an effective field of
0.018 V/$\mathrm{\SIUnitSymbolAngstrom}$, the gap at the $\Gamma$ point is 65
meV and the bottom of the inverted gap is found 66 meV above the Fermi level.
For the GaAs-type interface, shown in the SI, the band gap also decreases as
the electric field increases. However, because the GaAs-type interface has a
larger band gap and the effect of the electric field is weaker than for the
InSb-type interface, the gap does not close even for an input electric field
as high as 0.55 V/$\mathrm{\SIUnitSymbolAngstrom}$.

Fig 9 shows the change in the band gap, $\Delta$, as a function of the input
electric field, $E_{in}$, for 10/10 GaAs-type and InSb-type interfaces:

$\Delta=Gap(E_{in})-Gap(E_{in}=0)$ (3)

The blue and orange dashed lines indicate the band gaps of the 10/10 GaAs-type
and InSb-type interfaces, respectively. The gap closes when the dashed line is
crossed. To estimate the input electric field that would be required for the
gap to close for a 50/50 interface, we assume that the change in the gap would
behave similarly to a 10/10 interface. The green and red dashed lines indicate
the band gaps of the 50/50 GaAs-type and InSb-type interfaces, respectively.
Based on this, we estimate that an input electric field of 0.19
V/$\mathrm{\SIUnitSymbolAngstrom}$, which corresponds to an effective electric
field of 0.012 V/$\mathrm{\SIUnitSymbolAngstrom}$, would be needed to close
the gap for a 50/50 InSb-type interface, as indicated by the red solid line.
For the GaAs-type interface an input electric field of 0.55
V/$\mathrm{\SIUnitSymbolAngstrom}$, which corresponds to an effective electric
field of 0.017 V/$\mathrm{\SIUnitSymbolAngstrom}$, would be needed to close
the gap, as indicated by the green solid line. We highlight that the effective
electric field of 0.017 V/$\mathrm{\SIUnitSymbolAngstrom}$ corresponds to a
potential drop of $2.55$ V over the 15 nm thickness of the QW in this case,
which is likely to make the material conducting well before the topological
transition. Our results indicate that while it may be possible to tune the
InAs/GaSb interface into the topological regime, it would not be trivial.

Figure 8:  The inverted band gap at $\Gamma$ and its position above the Fermi
level as a function of the input electric field, where the shift is defined as
the energy difference between the position of the bottom of the inverted band
gap and the Fermi level at the $\Gamma$ point. Figure 9:  The band gap
reduction, $\Delta$, as a function of electric field for 10/10 GaAs-type and
InSb-type interfaces.

##  IV Conclusion

In summary, we have studied the HgTe/CdTe and InAs/GaSb quantum wells using
DFT with a Hubbard U correction determined by Bayesian optimization. DFT+U(BO)
produces band structures of comparable accuracy to a hybrid functional at the
computational cost of a semi-local functional. This enables us to conduct
simulations of large interface models with hundreds of atoms.

For the HgTe/CdTe interface we find that an inverted gap opens at a critical
thickness of 5.1 nm of HgTe, in agreement with experimental observations and
previous theoretical studies. For InAs/GaSb QWs with 50 layers (about 15 nm)
of GaSb we find that the gap narrows with increasing thickness of InAs in
agreement with the previous theory estimations. However, the gap does not
completely close with up to 50 layers (about 15 nm) of InAs. Based on the rate
of gap narrowing, we estimate that it would either require a significantly
thicker InAs film to close the gap or the gap would decay to a finite
asymptotic limit.

Simulations with an external electric field applied perpendicular to the
interface, pointing from GaSb to InAs, have been conducted for models with 10
layers of each material. We find that with increasing field strength the GaSb
VBM shifts upwards relative to the InAs CBM, leading to narrowing of the gap
at the interface. For the InSb-type interface, band crossing is observed at a
critical field and subsequently an inverted gap opens. As the electric field
increases the gap increases but also shifts higher in energy above the Fermi
level. Because the 10/10 interface has a larger gap due to the quantum size
effect, we estimate the reduced critical field that would be required to
achieve band inversion and reach the topological regime for thicker QWs
comprising 50 layers of each material.

Our results explain the difficulty of experimentally reaching the topological
regime in InAs/GaSb QWs. In principle, under the right conditions, an inverted
gap could be produced in this system. However, achieving this requires a
delicate balance between several parameters. To tune the initial gap, the
structure of the QWs must be precisely controlled, including the layer
thickness, the bonding configuration at the interface, and possibly also the
lattice strain. Even if a smaller zero-field gap is obtained by interface
engineering, a considerable electric field may still be required to obtain
band crossing and drive the system into the topological regime. Finally,
gating or doping may be required to tune the Fermi level position inside the
inverted gap.

The HgTe/CdTe QW does not suffer from this difficulty because HgTe has an
inverted band structure intrinsically. Therefore, no electric field is
necessary to achieve band inversion at the HgTe/CdTe interface and it is
easier to reach the topological regime. Thus, our results make a case for
limited applicability of InAs/GaSb quantum wells for 2DTI production and
suggest that alternative, more promising materials should be sought.

###### Acknowledgements.

We would like to thank Sergey Frolov from the University of Pittsburgh, Chris
Palmstrøm from the University of California, Santa Barbara, Vlad Pribiag from
the University of Minnesota, and Michael Wimmer from TU Delft for helpful
discussions. Work at CMU was funded by the National Science Foundation (NSF)
through grant OISE-1743717. This research used resources of the National
Energy Research Scientific Computing Center (NERSC), a DOE Office of Science
User Facility supported by the Office of Science of the U.S. Department of
Energy under contract no. DE-AC02-05CH11231.

## References

  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: topological insulators, Reviews of modern physics 82, 3045 (2010).
  * Ando [2013] Y. Ando, Topological insulator materials, Journal of the Physical Society of Japan 82, 102001 (2013).
  * Kou  _et al._ [2017] L. Kou, Y. Ma, Z. Sun, T. Heine, and C. Chen, Two-dimensional topological insulators: Progress and prospects, The Journal of Physical Chemistry Letters 8, 1905 (2017), pMID: 28394616, https://doi.org/10.1021/acs.jpclett.7b00222 .
  * Kane and Mele [2005] C. L. Kane and E. J. Mele, Quantum spin hall effect in graphene, Phys. Rev. Lett. 95, 226801 (2005).
  * Bernevig and Zhang [2006] B. A. Bernevig and S.-C. Zhang, Quantum spin hall effect, Phys. Rev. Lett. 96, 106802 (2006).
  * Bernevig  _et al._ [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Quantum spin hall effect and topological phase transition in hgte quantum wells, science 314, 1757 (2006).
  * König  _et al._ [2007] M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Quantum spin hall insulator state in hgte quantum wells, Science 318, 766 (2007).
  * Vidal  _et al._ [2012] J. Vidal, X. Zhang, V. Stevanović, J.-W. Luo, and A. Zunger, Large insulating gap in topological insulators induced by negative spin-orbit splitting, Phys. Rev. B 86, 075316 (2012).
  * Küfner and Bechstedt [2014] S. Küfner and F. Bechstedt, Topological transition and edge states in hgte quantum wells from first principles, Physical Review B 89, 195312 (2014).
  * Yang  _et al._ [2012] F. Yang, L. Miao, Z. Wang, M.-Y. Yao, F. Zhu, Y. Song, M.-X. Wang, J.-P. Xu, A. V. Fedorov, Z. Sun, _et al._ , Spatial and energy distribution of topological edge states in single bi (111) bilayer, Physical review letters 109, 016801 (2012).
  * Li  _et al._ [2016] X.-B. Li, W.-K. Huang, Y.-Y. Lv, K.-W. Zhang, C.-L. Yang, B.-B. Zhang, Y. Chen, S.-H. Yao, J. Zhou, M.-H. Lu, _et al._ , Experimental observation of topological edge states at the surface step edge of the topological insulator zrte 5, Physical review letters 116, 176803 (2016).
  * Wu  _et al._ [2016] R. Wu, J.-Z. Ma, S.-M. Nie, L.-X. Zhao, X. Huang, J.-X. Yin, B.-B. Fu, P. Richard, G.-F. Chen, Z. Fang, _et al._ , Evidence for topological edge states in a large energy gap near the step edges on the surface of zrte 5, Physical Review X 6, 021017 (2016).
  * Pauly  _et al._ [2015] C. Pauly, B. Rasche, K. Koepernik, M. Liebmann, M. Pratzer, M. Richter, J. Kellner, M. Eschbach, B. Kaufmann, L. Plucinski, _et al._ , Subnanometre-wide electron channels protected by topology, Nature Physics 11, 338 (2015).
  * Reis  _et al._ [2017] F. Reis, G. Li, L. Dudy, M. Bauernfeind, S. Glass, W. Hanke, R. Thomale, J. Schäfer, and R. Claessen, Bismuthene on a sic substrate: A candidate for a high-temperature quantum spin hall material, Science 357, 287 (2017).
  * Zhu  _et al._ [2019] S.-Y. Zhu, Y. Shao, E. Wang, L. Cao, X.-Y. Li, Z.-L. Liu, C. Liu, L.-W. Liu, J.-O. Wang, K. Ibrahim, _et al._ , Evidence of topological edge states in buckled antimonene monolayers, Nano letters 19, 6323 (2019).
  * Deng  _et al._ [2018] J. Deng, B. Xia, X. Ma, H. Chen, H. Shan, X. Zhai, B. Li, A. Zhao, Y. Xu, W. Duan, _et al._ , Epitaxial growth of ultraflat stanene with topological band inversion, Nature materials 17, 1081 (2018).
  * Wang  _et al._ [2016] Z. Wang, H. Zhang, D. Liu, C. Liu, C. Tang, C. Song, Y. Zhong, J. Peng, F. Li, C. Nie, _et al._ , Topological edge states in a high-temperature superconductor fese/srtio 3 (001) film, Nature materials 15, 968 (2016).
  * Tang  _et al._ [2017] S. Tang, C. Zhang, D. Wong, Z. Pedramrazi, H.-Z. Tsai, C. Jia, B. Moritz, M. Claassen, H. Ryu, S. Kahn, _et al._ , Quantum spin hall state in monolayer 1t’-wte 2, Nature Physics 13, 683 (2017).
  * Jia  _et al._ [2017] Z.-Y. Jia, Y.-H. Song, X.-B. Li, K. Ran, P. Lu, H.-J. Zheng, X.-Y. Zhu, Z.-Q. Shi, J. Sun, J. Wen, _et al._ , Direct visualization of a two-dimensional topological insulator in the single-layer 1 t’- wt e 2, Physical Review B 96, 041108 (2017).
  * Fei  _et al._ [2017] Z. Fei, T. Palomaki, S. Wu, W. Zhao, X. Cai, B. Sun, P. Nguyen, J. Finney, X. Xu, and D. H. Cobden, Edge conduction in monolayer wte 2, Nature Physics 13, 677 (2017).
  * Wu  _et al._ [2018] S. Wu, V. Fatemi, Q. D. Gibson, K. Watanabe, T. Taniguchi, R. J. Cava, and P. Jarillo-Herrero, Observation of the quantum spin hall effect up to 100 kelvin in a monolayer crystal, Science 359, 76 (2018).
  * Sajadi  _et al._ [2018] E. Sajadi, T. Palomaki, Z. Fei, W. Zhao, P. Bement, C. Olsen, S. Luescher, X. Xu, J. A. Folk, and D. H. Cobden, Gate-induced superconductivity in a monolayer topological insulator, Science 362, 922 (2018).
  * Lakrimi  _et al._ [1997] M. Lakrimi, S. Khym, R. J. Nicholas, D. M. Symons, F. M. Peeters, N. J. Mason, and P. J. Walker, Minigaps and novel giant negative magnetoresistance in inas/gasb semimetallic superlattices, Phys. Rev. Lett. 79, 3034 (1997).
  * Cooper  _et al._ [1998] L. J. Cooper, N. K. Patel, V. Drouot, E. H. Linfield, D. A. Ritchie, and M. Pepper, Resistance resonance induced by electron-hole hybridization in a strongly coupled inas/gasb/alsb heterostructure, Phys. Rev. B 57, 11915 (1998).
  * Liu  _et al._ [2008] C. Liu, T. L. Hughes, X.-L. Qi, K. Wang, and S.-C. Zhang, Quantum spin hall effect in inverted type-ii semiconductors, Physical review letters 100, 236601 (2008).
  * Suzuki  _et al._ [2013] K. Suzuki, Y. Harada, K. Onomitsu, and K. Muraki, Edge channel transport in the inas/gasb topological insulating phase, Phys. Rev. B 87, 235311 (2013).
  * Suzuki  _et al._ [2015] K. Suzuki, Y. Harada, K. Onomitsu, and K. Muraki, Gate-controlled semimetal-topological insulator transition in an inas/gasb heterostructure, Phys. Rev. B 91, 245309 (2015).
  * Knez  _et al._ [2011] I. Knez, R.-R. Du, and G. Sullivan, Evidence for helical edge modes in inverted $\mathrm{InAs}/\mathrm{GaSb}$ quantum wells, Phys. Rev. Lett. 107, 136603 (2011).
  * Pribiag  _et al._ [2015] V. S. Pribiag, A. J. Beukman, F. Qu, M. C. Cassidy, C. Charpentier, W. Wegscheider, and L. P. Kouwenhoven, Edge-mode superconductivity in a two-dimensional topological insulator, Nature nanotechnology 10, 593 (2015).
  * Knez  _et al._ [2014] I. Knez, C. T. Rettner, S.-H. Yang, S. S. P. Parkin, L. Du, R.-R. Du, and G. Sullivan, Observation of edge transport in the disordered regime of topologically insulating $\mathrm{InAs}/\mathrm{GaSb}$ quantum wells, Phys. Rev. Lett. 112, 026602 (2014).
  * Nichele  _et al._ [2016] F. Nichele, H. J. Suominen, M. Kjaergaard, C. M. Marcus, E. Sajadi, J. A. Folk, F. Qu, A. J. Beukman, F. K. de Vries, J. van Veen, _et al._ , Edge transport in the trivial phase of inas/gasb, New Journal of Physics 18, 083005 (2016).
  * Du  _et al._ [2015] L. Du, I. Knez, G. Sullivan, and R.-R. Du, Robust helical edge transport in gated $\mathrm{InAs}/\mathrm{GaSb}$ bilayers, Phys. Rev. Lett. 114, 096802 (2015).
  * Mueller  _et al._ [2015] S. Mueller, A. N. Pal, M. Karalic, T. Tschirky, C. Charpentier, W. Wegscheider, K. Ensslin, and T. Ihn, Nonlocal transport via edge states in inas/gasb coupled quantum wells, Phys. Rev. B 92, 081303 (2015).
  * Shojaei  _et al._ [2018] B. Shojaei, A. P. McFadden, M. Pendharkar, J. S. Lee, M. E. Flatté, and C. J. Palmstrøm, Materials considerations for forming the topological insulator phase in inas/gasb heterostructures, Physical Review Materials 2, 064603 (2018).
  * Sengupta  _et al._ [2013] P. Sengupta, T. Kubis, Y. Tan, M. Povolotskyi, and G. Klimeck, Design principles for hgte based topological insulator devices, Journal of Applied Physics 114, 043702 (2013).
  * Livneh  _et al._ [2012] Y. Livneh, P. C. Klipstein, O. Klin, N. Snapi, S. Grossman, A. Glozman, and E. Weiss, $\mathbf{k}\cdot{}\mathbf{p}$ model for the energy dispersions and absorption spectra of inas/gasb type-ii superlattices, Phys. Rev. B 86, 235311 (2012).
  * Wang  _et al._ [1999] L.-W. Wang, S.-H. Wei, T. Mattila, A. Zunger, I. Vurgaftman, and J. R. Meyer, Multiband coupling and electronic structure of $(\mathrm{InAs}{)}_{n}/(\mathrm{GaSb}{)}_{n}$ superlattices, Phys. Rev. B 60, 5590 (1999).
  * Li  _et al._ [2018] C.-A. Li, S.-B. Zhang, and S.-Q. Shen, Hidden edge dirac point and robust quantum edge transport in inas/gasb quantum wells, Physical Review B 97, 045420 (2018).
  * Skolasinski  _et al._ [2018] R. Skolasinski, D. I. Pikulin, J. Alicea, and M. Wimmer, Robust helical edge transport in quantum spin hall quantum wells, Physical Review B 98, 201404 (2018).
  * Luo and Zunger [2010] J.-W. Luo and A. Zunger, Design principles and coupling mechanisms in the 2d quantum well topological insulator $\mathrm{HgTe}/\mathrm{CdTe}$, Phys. Rev. Lett. 105, 176805 (2010).
  * Magri and Zunger [2003] R. Magri and A. Zunger, Predicting interband transition energies for inas/gasb superlattices using the empirical pseudopotential method, Phys. Rev. B 68, 155329 (2003).
  * Fu and Kane [2007] L. Fu and C. L. Kane, Topological insulators with inversion symmetry, Phys. Rev. B 76, 045302 (2007).
  * Wei and Razeghi [2004] Y. Wei and M. Razeghi, Modeling of type-ii inas/gasb superlattices using an empirical tight-binding method and interface engineering, Phys. Rev. B 69, 085316 (2004).
  * Küfner and Bechstedt [2014] S. Küfner and F. Bechstedt, Topological transition and edge states in hgte quantum wells from first principles, Phys. Rev. B 89, 195312 (2014).
  * Anversa  _et al._ [2014] J. Anversa, P. Piquini, A. Fazzio, and T. M. Schmidt, First-principles study of hgte/cdte heterostructures under perturbations preserving time-reversal symmetry, Phys. Rev. B 90, 195311 (2014).
  * Krukau  _et al._ [2006] A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, Influence of the exchange screening parameter on the performance of screened hybrid functionals, The Journal of Chemical Physics 125, 224106 (2006).
  * Wang and Zhang [2014] J. Wang and Y. Zhang, Band-gap corrected density functional theory calculations for inas/gasb type ii superlattices, Journal of Applied Physics 116, 214301 (2014).
  * Sun  _et al._ [2011] W.-F. Sun, M.-C. Li, and L.-C. Zhao, First-principles study of interface relaxation effect on interface and electronic structures of inas/gasb superlattices with different interface types, Superlattices and Microstructures 49, 81 (2011).
  * Garwood  _et al._ [2017] T. Garwood, N. A. Modine, and S. Krishna, Electronic structure modeling of inas/gasb superlattices with hybrid density functional theory, Infrared Physics & Technology 81, 27 (2017).
  * Taghipour  _et al._ [2018] Z. Taghipour, E. Shojaee, and S. Krishna, Many-body perturbation theory study of type-ii inas/gasb superlattices within the gw approximation, Journal of Physics: Condensed Matter 30, 325701 (2018).
  * Yu  _et al._ [2020] M. Yu, S. Yang, C. Wu, and N. Marom, Machine learning the hubbard u parameter in dft+ u using bayesian optimization, npj Computational Materials 6, 1 (2020).
  * Dudarev  _et al._ [1998] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: An lsda+u study, Phys. Rev. B 57, 1505 (1998).
  * Micnas  _et al._ [1990] R. Micnas, J. Ranninger, and S. Robaszkiewicz, Superconductivity in narrow-band systems with local nonretarded attractive interactions, Rev. Mod. Phys. 62, 113 (1990).
  * Hase and Yanagisawa [2007] I. Hase and T. Yanagisawa, Madelung energy of the valence-skipping compound ba bi o 3, Physical Review B 76, 174103 (2007).
  * Nakamura  _et al._ [2009] H. Nakamura, N. Hayashi, N. Nakai, M. Okumura, and M. Machida, First-principle electronic structure calculations for magnetic moment in iron-based superconductors: An LSDA+ negative U study, Physica C: Superconductivity 469, 908 (2009).
  * Persson and Mirbt [2006] C. Persson and S. Mirbt, Improved electronic structure and optical properties of sp-hybridized semiconductors using LDA+U SIC, Brazilian journal of physics 36, 286 (2006).
  * Cococcioni [2012] M. Cococcioni, The LDA+U approach: a simple Hubbard correction for correlated ground states, Correlated Electrons: From Models to Materials Modeling and Simulation; Verlag des Forschungszentrum Jülich: Jülich, Germany (2012).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B 47, 558 (1993).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Perdew  _et al._ [1997] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple [phys. rev. lett. 77, 3865 (1996)], Phys. Rev. Lett. 78, 1396 (1997).
  * Steiner  _et al._ [2016] S. Steiner, S. Khmelevskyi, M. Marsmann, and G. Kresse, Calculation of the magnetic anisotropy with projected-augmented-wave methodology and the case study of disordered ${\mathrm{fe}}_{1-x}{\mathrm{co}}_{x}$ alloys, Phys. Rev. B 93, 224425 (2016).
  * Neugebauer and Scheffler [1992] J. Neugebauer and M. Scheffler, Adsorbate-substrate and adsorbate-adsorbate interactions of na and k adlayers on al(111), Phys. Rev. B 46, 16067 (1992).
  * Yang  _et al._ [2020a] S. Yang, N. Schröter, S. Schuwalow, M. Rajpalk, K. Ohtani, P. KrogstrupGeorg, G. W. Winkler, J. Gukelberger, D. Gresch, G. Aeppli, R. M. Lutchyn, V. N. Strocov, and N. Marom, Electronic structure of inas and insb surfaces: density functional theory and angle-resolved photoemission spectroscopy, arXiv , 2012.14935 (2020a).
  * Broser  _et al._ [1982] I. Broser, R. Broser, H. Finkenrath, R. Galazka, H. Gumlich, A. Hoffmann, J. Kossut, E. Mollwo, H. Nelkowski, G. Nimtz, _et al._ ,  _Physics of II-VI and I-VII Compounds, Semimagnetic Semiconductors/Physik Der II-VI und I-VII-Verbindungen, Semimagnetische Halbleiter Elemente_ , Vol. 17 (Springer, 1982).
  * Wu  _et al._ [2015] Y. Wu, G. Chen, Y. Zhu, W.-J. Yin, Y. Yan, M. Al-Jassim, and S. J. Pennycook, Lda+ u/gga+ u calculations of structural and electronic properties of cdte: Dependence on the effective u parameter, Computational Materials Science 98, 18 (2015).
  * Heyd  _et al._ [2006] J. Heyd, G. E. Scuseria, and M. Ernzerhof, Erratum: “hybrid functionals based on a screened coulomb potential” [j. chem. phys. 118, 8207 (2003)], The Journal of Chemical Physics 124, 219906 (2006).
  * Hey [2003] Hybrid functionals based on a screened Coulomb potential, J. Chem. Phys. 118, 8207 (2003).
  * Vurgaftman  _et al._ [2001] I. Vurgaftman, J. á. Meyer, and L. á. Ram-Mohan, Band parameters for iii–v compound semiconductors and their alloys, Journal of applied physics 89, 5815 (2001).
  * Heyd  _et al._ [2005] J. Heyd, J. E. Peralta, G. E. Scuseria, and R. L. Martin, Energy band gaps and lattice parameters evaluated with the heyd-scuseria-ernzerhof screened hybrid functional, The Journal of chemical physics 123, 174101 (2005).
  * Wei  _et al._ [2000] S.-H. Wei, S. Zhang, and A. Zunger, First-principles calculation of band offsets, optical bowings, and defects in cds, cdse, cdte, and their alloys, Journal of applied Physics 87, 1304 (2000).
  * Kick  _et al._ [2019] M. Kick, K. Reuter, and H. Oberhofer, Intricacies of DFT+U, Not Only in a Numeric Atom Centered Orbital Framework, Journal of Chemical Theory and Computation 15, 1705 (2019).
  * Jiang  _et al._ [2010] H. Jiang, R. I. Gomez-Abal, P. Rinke, and M. Scheffler, First-principles modeling of localized d states with the g w@ lda+ u approach, Physical Review B 82, 045108 (2010).
  * West [1999] A. R. West,  _Basic solid state chemistry_ (John Wiley & Sons Incorporated, 1999).
  * Yang  _et al._ [2020b] S. Yang, C. Wu, and N. Marom, Topological properties of snse/eus and snte/cate interfaces, Phys. Rev. Materials 4, 034203 (2020b).
  * Klein  _et al._ [2012] B. Klein, N. Gautam, S. Myers, and S. Krishna, Temperature-dependent absorption derivative on inas/gasb type ii superlattices, in  _Infrared Technology and Applications XXXVIII_ , Vol. 8353 (International Society for Optics and Photonics, 2012) p. 83530X.
  * Akiho  _et al._ [2016] T. Akiho, F. Couëdo, H. Irie, K. Suzuki, K. Onomitsu, and K. Muraki, Engineering quantum spin hall insulators by strained-layer heterostructures, Applied Physics Letters 109, 192105 (2016).
  * Du  _et al._ [2017] L. Du, T. Li, W. Lou, X. Wu, X. Liu, Z. Han, C. Zhang, G. Sullivan, A. Ikhlassi, K. Chang, _et al._ , Tuning edge states in strained-layer inas/gainsb quantum spin hall insulators, Physical Review Letters 119, 056803 (2017).
  * Tiemann  _et al._ [2017] L. Tiemann, S. Mueller, Q.-S. Wu, T. Tschirky, K. Ensslin, W. Wegscheider, M. Troyer, A. A. Soluyanov, and T. Ihn, Impact of strain on the electronic properties of inas/gasb quantum well systems, Physical Review B 95, 115108 (2017).
  * Qu  _et al._ [2015] F. Qu, A. J. Beukman, S. Nadj-Perge, M. Wimmer, B.-M. Nguyen, W. Yi, J. Thorp, M. Sokolich, A. A. Kiselev, M. J. Manfra, _et al._ , Electric and magnetic tuning between the trivial and topological phases in inas/gasb double quantum wells, Physical Review Letters 115, 036803 (2015).
  * Wang  _et al._ [2018] X.-P. Wang, X.-B. Li, N.-K. Chen, J.-H. Zhao, Q.-D. Chen, and H.-B. Sun, Electric field analyses on monolayer semiconductors: the example of inse, Physical Chemistry Chemical Physics 20, 6945 (2018).
  * Naveh and Laikhtman [1995] Y. Naveh and B. Laikhtman, Band-structure tailoring by electric field in a weakly coupled electron-hole system, Applied physics letters 66, 1980 (1995).