# Coupling between tilts and charge carriers at polar-nonpolar perovskite
interfaces

Daniel Bennett<EMAIL_ADDRESS>Theory of Condensed Matter Group, Cavendish
Laboratory, University of Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE,
United Kingdom  Physique Théorique des Matériaux, QMAT, CESAM, University of
Liège, B-4000 Sart-Tilman, Belgium  Pablo Aguado-Puente  CIC Nanogune BRTA,
Tolosa Hiribidea 76, 20018 San Sebastian, Spain  Emilio Artacho  Theory of
Condensed Matter Group, Cavendish Laboratory, University of Cambridge, J. J.
Thomson Avenue, Cambridge CB3 0HE, United Kingdom  CIC Nanogune BRTA, Tolosa
Hiribidea 76, 20018 San Sebastian, Spain  Donostia International Physics
Center, Paseo Manuel de Lardizabal 4, 20018 San Sebastian, Spain  Ikerbasque,
Basque Foundation for Science, 48011 Bilbao, Spain  Nicholas Bristowe  School
of Physical Sciences, University of Kent, Canterbury, Kent CT2 7NH, United
Kingdom  Centre for Materials Physics, Durham University, South Road, Durham
DH1 3LE, United Kingdom

###### Abstract

The phenomenological theory for the polar instability giving rise to a two-
dimensional electron gas at perovskite interfaces is hereby extended to
include the coupling to perovskite tilts. A Landau theory for homogeneous
tilts is first explored, setting the scene for the further, more realistic
Landau-Ginzburg theory describing varying tilt amplitudes across a thin film.
The theory is also generalized to account for the response to an applied
electric field normal to the interface, which allows a finer control on phase
transitions. The conventionally described physics of a single metal-insulator
transition is substantially enriched by the coupling, the model describing
various scenarios with one or two transitions, possibly continuous or
discontinuous. First-principles calculations permit the estimation of the
parameters defining the model, which have been calculated for the interface
between lanthanum aluminate and strontium titanate.

##  I Introduction

The discovery of an insulator–metal transition at the interface between two
insulating perovskites sparked considerable interest in the scientific
community Ohtomo  _et al._ (2002); Ohtomo and Hwang (2004). When thin films of
lanthanum aluminate ($\mathrm{LaAlO_{3}}$, LAO) were grown on substrates of
strontium titanate ($\mathrm{SrTiO_{3}}$, STO), a two-dimensional electron gas
(2DEG) was found to appear at the interface in order to screen the polar
discontinuity there. This 2DEG has been found to be associated with
interesting phenomena such as enhanced capacitance Li  _et al._ (2011a),
superconductivity Reyren  _et al._ (2007) and magnetism Brinkman  _et al._
(2007), even at the same time Li  _et al._ (2011b); Bert  _et al._ (2011). It
also has potential for applications in field effect transistors (FET) Thiel
_et al._ (2006); Cen  _et al._ (2008); Förg  _et al._ (2012), sensors Xie  _et
al._ (2011), photodetection Irvin  _et al._ (2010), thermoelectrics Pallecchi
_et al._ (2010); Filippetti  _et al._ (2012) and solar cells Assmann  _et al._
(2013); Liang  _et al._ (2013).

The physical origin and character of this insulator–metal transition have been
debated for many years, and discussed in several reviews Hwang (2006); Pauli
and Willmott (2008); Huijben  _et al._ (2009); Mannhart and Schlom (2010);
Pentcheva and Pickett (2010); Chen  _et al._ (2010); Warusawithana  _et al._
(2013); Chambers (2011); Pentcheva  _et al._ (2012); Gabay  _et al._ (2013);
Stemmer and James Allen (2014); Bristowe  _et al._ (2014). In response to a
polar discontinuity at the interface, first proposed in the original papers on
LAO/STO Ohtomo  _et al._ (2002); Ohtomo and Hwang (2004), two of the most
popular theories for the appearance of the 2DEG are indirectly supported by
experimental evidence: the first is electronic reconstruction, where the 2DEG
forms via a transfer of electrons from the valence band at the surface of the
thin film to the conduction band at the interface in order to screen the polar
discontinuity Nakagawa  _et al._ (2006); Popović  _et al._ (2008); Pentcheva
and Pickett (2009). The second achieves the same screening by means of redox
defects, where oxygen vacancies or hydrogen adatoms, for example, form at the
surface, creating charge carriers that can move to the interface Kalabukhov
_et al._ (2007); Siemons  _et al._ (2007); Herranz  _et al._ (2007); Eckstein
(2007); Shibuya  _et al._ (2007); Chen  _et al._ (2011); Bristowe  _et al._
(2011); Lee  _et al._ (2012). In Ref. Liu  _et al._ , 2013 it was predicted
that both mechanisms are possible, depending on the experimental conditions
such as the oxygen pressure and the growth temperature Lee  _et al._ (2012);
Li  _et al._ (2015).

The character of the transition with film thickness, i.e. whether the carriers
appear continuously or discontinuously after a critical thickness, is still
debated. It has been observed in experimental studies Ohtomo and Hwang (2004)
and suggested in theoretical studies Bristowe  _et al._ (2014); Yu and Zunger
(2014) based on redox effects that a discontinuous transition occurs at a
critical thickness $d_{\text{c}}$ between 3 and 4 unit cells of LAO. However,
film thickness is not a good parameter with which to make any conclusions
about the order of the transition, since it is fixed for each sample and is
changed discretely by a number of unit cells. Thus, it is impossible to
conclude whether the transition is continuous or discontinuous, and the
problem is only of theoretical interest. A more realistic approach to
investigate this would be to apply an electric field, using top and back gates
Chen  _et al._ (2016), to a sample which is close to the critical thickness.
An electric field can be used to switch the 2DEG on Kim  _et al._ (2015) and
off Hosoda  _et al._ (2013) in a single sample, which is a desirable feature
in practical applications. It can also enhance the properties of the 2DEG such
as the superconductivity Bert  _et al._ (2012); Eerkes  _et al._ (2013); Davis
_et al._ (2018) magneto-transport Liu  _et al._ (2015) and optical behavior
Lei  _et al._ (2014). An electric field could be tuned with more precision and
thus would be more suitable for studying the character of the insulator–metal
transition experimentally.

A phenomenological theory at the mean-field level, which treats the carriers
as a homogeneous charge distribution $\sigma$, predicts that the onset of
carriers with thickness is continuous, with a critical exponent for $\sigma$
of 1 Bristowe  _et al._ (2014). When thinking about the redox defects proposed
in Ref. Bristowe  _et al._ , 2014, an assumption of non-interacting defects
suggests a drastically discontinuous transition, switching on directly to
$\sigma\sim P_{s}$ Yu and Zunger (2014), $P_{s}$ being the polar
discontinuity, half an electron per primitive unit cell surface area for
LAO/STO; when it becomes favorable for one vacancy-carrier pair to form, it is
favorable for all of them to form, giving the mentioned discontinuity at
$d_{\text{c}}$. One can go beyond the mean-field level by considering the
interactions between the traps Bristowe  _et al._ (2011). The vacancy at the
surface and carrier at the interface act like a dipole, and thus there would
be dipole-dipole interactions between the vacancies. This predicts a
transition which is still continuous, but with a critical exponent of
$\frac{2}{3}$.

Something that has to our knowledge not yet been considered is the influence
of other phase transitions, such as the antiferrodistortive (AFD) rotations of
oxygen octahedra (tilts), which both LAO and STO can exhibit. Tilts compete
with the polar mode, so it is reasonable to expect that they may indirectly
interact with the carriers, which appear to screen the polar discontinuity.
The competition between tilts and the polar mode in LAO/STO has been
considered in a previous study Gazquez  _et al._ (2017), although the
influence of the tilts on the appearance of the carriers at the
phenomenological level has not been considered.

In this paper we generalize the phenomenological model of carrier formation at
polar-nonpolar perovskite interfaces to account for coupling to tilts in the
thin film. We show that, upon coupling to homogeneous tilts, four new distinct
types of transitions are possible, depending only on the energetics of the
tilts, polar discontinuity and the coupling between them. These include
continuous and discontinuous transitions of the carriers, facilitated by the
tilts, and both simultaneous and distinct transitions of tilts and carriers.
Using first-principles calculations, we can make predictions about the type of
transition which occurs at the LAO/STO interface. We then generalize the model
to allow for inhomogeneous tilts in the polar thin film, using Ginzburg-Landau
theory. The inhomogeneity of the tilts is determined by the correlation length
in the film and the extrapolation lengths, which describe the relative energy
differences between the tilts in the interior and at the surfaces. We show
that two of the transitions predicted at the homogeneous level are unaffected
by the inhomogeneity of the tilts, other than the values at which the
transitions occur being renormalized. For the other two transitions, their
character can be changed by decreasing the correlation length, leading to two
entirely new types of transitions which are not possible for homogeneous
tilts.

##  II Coupling to homogeneous tilts

###  II.1 Phenomenological theory of carrier formation

Figure 1: Sketch of the electrostatics of a polar-nonpolar interface. A thin
film of thickness $d$ is placed on top of a semi-infinite substrate. An
electric field is applied using a back gate under the substrate and a biased
tip. For simplicity, we assume that the tip is much wider than the thickness
of the film, and that the applied field is homogeneous. A 2DEG forms at the
interface and a corresponding 2DHG or redox defects form at the surface of the
thin film as the combined result of the polar discontinuity and applied field.
$D_{\text{film}}$ is the displacement field of the film and $P_{s}$ is the
polar discontinuity. $D_{\text{air}}$ is the displacement field of the region
of air between surface of the film and the tip, and $D_{\text{sub}}$ is the
displacement field of the substrate.

First, we review the phenomenological model of carrier formation at polar-
nonpolar perovskite interfaces Bristowe  _et al._ (2014), introducing the
applied electric field in a similar manner to Ref. Aguado-Puente  _et al._
(2015). A sketch of the system is shown in Fig. 1. We have a polar thin film
of thickness $d$ on top of a non-polar substrate, with e.g. air above, similar
to the geometry considered for ferroelectric thin films in Ref. Bennett  _et
al._ (2020). Experimentally, a potential bias is usually applied between the
tip of an atomic force microscope (AFM) at the surface of the film and a back
gate at the bottom of the substrate Caviglia  _et al._ (2008); Stengel  _et
al._ (2009); Förg  _et al._ (2012); Hosoda  _et al._ (2013). Realistically,
this would result in inhomogeneous electric fields which decay away from the
tip. However, we note for example that, in the case of fields applied with an
AFM tip, the lateral extent over which the field can be considered homogeneous
($\sim 100$ of nm) is much larger than the usual thickness of the polar film
(few nm). This is at least in part due to the meniscus of water and other
surface adsorbates that form around the AFM tip Paruch and Guyonnet (2013). In
the interest of simplicity we consider an homogeneous field, assuming the tip
is much wider than the thickness of the film. With the usual LaO/TiO2
termination, a 2DEG of carrier concentration $\sigma$ forms between the film
and substrate – and a corresponding two-dimensional hole gas (2DHG) or redox
defects form at the surface – in order to screen the polar discontinuity.

The displacement fields in the film, substrate and air region between the film
and tip are
$D_{\text{film}}=\epsilon_{\text{film}}\mathcal{E}_{\text{film}}+P_{s}$,
${D_{\text{sub}}=\epsilon_{\text{sub}}}\mathcal{E}_{\text{sub}}$, and
$D_{\text{air}}=\epsilon_{\text{air}}\mathcal{E}_{\text{air}}$, respectively,
where $P_{s}$ is the polar discontinuity. We consider the displacement field
in the substrate as a tunable parameter with which the field (or potential
difference) applied to the system is measured. The displacement field in the
film could also be used Gazquez  _et al._ (2017). All vectors are normal to
the interface, so vector notation is omitted. From Maxwell’s equations, the
boundary condition at interface relates the displacement fields in the film
and the substrate to the free charge, $\sigma$:

$\begin{gathered}\left(D_{\text{sub}}-D_{\text{film}}\right)\cdot\hat{n}=-\sigma\\\
\epsilon_{\text{sub}}\mathcal{E}_{\text{sub}}-\epsilon_{\text{film}}\mathcal{E}_{\text{film}}-P_{s}=-\sigma\\\
\implies\mathcal{E}_{\text{film}}=\frac{\epsilon_{\text{sub}}}{\epsilon_{\text{film}}}\mathcal{E}_{\text{sub}}-\frac{1}{\epsilon_{\text{film}}}(P_{s}-\sigma)\end{gathered}\;,$
(1)

where $\hat{n}$ points outwards from the film. Assuming that the boundaries of
the film have equal and opposite carrier concentration, the boundary condition
at the surface between the film and the air region is also satisfied by
enforcing finite-$D$ boundary conditions: $D_{\text{sub}}=D_{\text{air}}$ (not
the electric fields).

The electrostatic energy of the film is Jackson (1999)

$\mathcal{F}_{\text{elec}}=\int_{0}^{D_{\text{film}}}\mathcal{E}_{\text{film}}\cdot\delta
D_{\text{film}}=\frac{1}{2}\epsilon_{\text{film}}\mathcal{E}_{\text{film}}^{2}\;.$
(2)

Using Eq. (1) we can write the electrostatic energy in terms of
$\mathcal{E}_{\text{sub}}$, $P_{s}$ and $\sigma$:

$\mathcal{F}_{\text{elec}}=\frac{1}{2}\frac{\epsilon_{\text{sub}}^{2}}{\epsilon_{\text{film}}}\mathcal{E}_{\text{sub}}^{2}-\frac{\epsilon_{\text{sub}}}{\epsilon_{\text{film}}}\mathcal{E}_{\text{sub}}\cdot(P_{s}-\sigma)+\frac{1}{2\epsilon_{\text{film}}}(P_{s}-\sigma)^{2}$
(3)

The first term corresponds to the energy of the applied field in the
substrate, which does not affect the carrier concentration and can be
neglected (the energy of the substrate is also neglected for the same reason).
The second term represents the coupling between the applied field and the
polar mode and the third term is the electrostatic energy of the polar
discontinuity. Writing Eq. (3) in terms of reduced variables
$\sigma^{\prime}=\frac{\sigma}{P_{s}}$ and
$\mathcal{E}^{\prime}=\frac{\mathcal{E}_{\text{sub}}}{(P_{s}/\epsilon_{\text{sub})}}$,
which is the displacement field in the substrate, in units of the polar
discontinuity, we get

$\mathcal{F}_{\text{elec}}(\sigma^{\prime},\mathcal{E}^{\prime})=\frac{P_{s}^{2}}{2\epsilon}(1-\sigma^{\prime})^{2}-\frac{P_{s}^{2}}{\epsilon}\mathcal{E}^{\prime}\cdot(1-\sigma^{\prime})\;,$
(4)

where we drop the subscript $\epsilon_{\text{film}}\to\epsilon$, since the
permittivity of the substrate has been absorbed into the scale for the applied
field and no longer appears in the free energy. For $\mathcal{E}^{\prime}=0$,
the energy is minimized when $\sigma^{\prime}=1$, i.e. the polar discontinuity
is fully screened by the carriers. For LAO/STO, this corresponds to a carrier
concentration of exactly half an electron per unit cell surface area. However,
we have neglected the formation energy of the carriers which appear the
interface. Including the cost of generating an electron-hole pair across a gap
$\Delta$, the free energy is

$\mathcal{F}_{\sigma}(\sigma^{\prime},\mathcal{E}^{\prime})=\frac{P_{s}^{2}}{\epsilon}\left[\frac{d_{\text{c}}}{d}\sigma^{\prime}+\frac{1}{2}(1-\sigma^{\prime})^{2}-\mathcal{E}^{\prime}\cdot(1-\sigma^{\prime})\right]\;,$
(5)

where $d_{c}=\frac{{\epsilon}_{0}\Delta}{P_{s}}$ is the critical thickness,
above which it is favorable for carriers to appear. Minimizing with respect to
$\sigma^{\prime}$ gives:

$\sigma^{\prime}(d,\mathcal{E}^{\prime})=1-\left(\frac{d_{\text{c}}}{d}+\mathcal{E}^{\prime}\right)\;,$
(6)

which implies that the carrier concentration increases with the thickness of
the film, and can be tuned, and even switched on/off with an electric field.
When an electric field is applied, the critical thickness for carriers to
appear is

$d_{\text{c}}(\mathcal{E}^{\prime})=\frac{d_{\text{c}}}{1-\mathcal{E}^{\prime}}\;,$
(7)

which can be reduced or increased, depending on the sign of
$\mathcal{E}^{\prime}$; the sign convention used for electric field is so that
it is aligned with the polar discontinuity when positive. In Fig. 2 (a) we
show the carrier concentration as a function of thickness at zero field, and
at positive and negative values of $\mathcal{E}^{\prime}$. When
$\mathcal{E}^{\prime}=0$, the appearance of carriers is as described in Ref.
Bristowe  _et al._ , 2014: $\sigma^{\prime}$ switches on at $d_{\text{c}}$ and
approaches 1 asymptotically from below. When a negative field is applied, the
polar discontinuity enhanced, and the carrier transition occurs at a reduced
thickness. In this case there is a second critical thickness where
$\sigma^{\prime}$ reaches 1 and saturates:

$d_{\text{c},2}(\mathcal{E}^{\prime})=-\frac{d_{\text{c}}}{\mathcal{E}^{\prime}}\;.$
(8)

Figure 2: (a): Carrier concentration as a function of thickness for zero
(black), positive and negative (red) values of reduced field
$\mathcal{E}^{\prime}$. When $\mathcal{E}^{\prime}$ is positive,
$d_{\text{c}}$ increases to $d_{\text{c}+}$. When $\mathcal{E}^{\prime}$ is
negative, $d_{\text{c}}$ is reduced to $d_{\text{c}-}$, and $\sigma^{\prime}$
saturates to 1 at a second critical thickness. (b): Reduced carrier
concentration with applied field at different thicknesses. The critical
thickness $d=d_{\text{c}}$ is shown in black, and the sub- and super-critical
thicknesses, $d=\frac{1}{2}d_{\text{c}}$ and $d=2d_{\text{c}}$, respectively,
are shown in red.

When $\mathcal{E}^{\prime}=-1$, the critical thicknesses are
$d_{\text{c}}(-1)=\frac{d_{\text{c}}}{2}$ and
$d_{\text{c},2}(-1)=d_{\text{c}}$, i.e. the carriers appear at half the
original critical thickness and saturate at the original critical thickness.
When a positive field is applied, the polar discontinuity is partially
screened, and the critical thickness increases.

The other possibility is to fix the thickness and induce the carriers with an
applied field. From Eq. (7), at a fixed thickness there is a critical field
value for which the carriers appear:

$\mathcal{E}_{\text{c}}^{\prime}(d)=1-\frac{d_{\text{c}}}{d}\;.$ (9)

If $d<d_{\text{c}}$, the film is sub-critical and
$\mathcal{E}_{\text{c}}^{\prime}(d)<0$. If $d\geq d_{\text{c}}$, the carriers
will already be present at zero field. Applying a positive field partially
screens the polar discontinuity, reducing the need for carriers. There is a
second critical field value at which the carrier concentration saturates:

$\mathcal{E}_{\text{c},2}^{\prime}(d)=-\frac{d_{\text{c}}}{d}\;,$ (10)

which is always negative for finite $d$. We have the relation between the two
critical field values:

$\mathcal{E}_{\text{c}}^{\prime}(d)-\mathcal{E}_{\text{c},2}^{\prime}(d)=1\;.$
(11)

The different scenarios for inducing or suppressing carriers with an applied
field are summarized in Fig. 2 (b).

###  II.2 Coupling to tilts

In this section we consider the effect of coupling to homogeneous tilts on the
formation of carriers. The simplest way to consider this is to add the
independent free energies of the tilts and carriers plus a coupling term:
$\mathcal{F}=\mathcal{F}_{\sigma}+\mathcal{F}_{\phi}+\mathcal{F}_{\sigma\phi}$.

Bulk LAO undergoes a transition from cubic to rhombohedral at $T_{C}\sim
541^{\circ}$\mathrm{C}$$ Lehnert  _et al._ (2000). The tilt pattern observed
is, in Glazer notation Glazer (1972), $a^{-}a^{-}a^{-}$. For a thin film of
LAO, the tilt pattern changes to $a^{-}a^{-}c^{0}$ due to the tensile biaxial
strain imposed by clamping to the STO substrate Hatt and Spaldin (2010);
Fister  _et al._ (2014). If a compressive biaxial strain were applied by using
a substrate with a smaller lattice constant than LAO, the observed tilt
pattern changes to $a^{0}a^{0}c^{-}$ Hatt and Spaldin (2010). For simplicity
we assume that the thin film undergoes a simple displacive transition from
untilted to tilted below some temperature $T_{C}$. Thus, the free energy of
tilts can be described by a double well:

$\mathcal{F}_{\phi}(\phi^{\prime},T)=\frac{1}{2}X_{\phi}^{-1}\left[\frac{1}{4}\phi^{\prime
4}-\frac{1}{2}\left(1-\frac{T}{T_{C}}\right)\phi^{\prime 2}\right]\;,$ (12)

where $\phi^{\prime}=\frac{\phi}{\phi_{0}}$ is the tilt angle in units of
$\phi_{0}$, the bulk equilibrium tilt angle, and $X_{\phi}^{-1}$ is the
curvature about the minima of the double well, in units of energy per unit
volume, playing the role of an inverse susceptibility.

The simplest coupling we can introduce is a biquadratic coupling between the
tilts and the polar mode Benedek and Fennie (2013); Salje (1992); Gazquez  _et
al._ (2017):

$\mathcal{F}_{\sigma\phi}=\frac{1}{2}A\phi^{\prime
2}\left(\frac{\mathcal{E}_{\text{film}}(\sigma^{\prime})}{P_{s}/\epsilon}\right)^{2}\;,$
(13)

where $A$ is the biquadratic coupling coefficient, in units of energy per unit
volume. Combining Eqs. (5), (12) and (13), the total free energy is

$\begin{gathered}\mathcal{F}_{\text{tot}}=\frac{P_{s}^{2}}{\epsilon}\left[\frac{d_{\text{c}}}{d}\sigma^{\prime}+\frac{1}{2}\left(1+\frac{A}{P_{s}^{2}/\epsilon}\phi^{\prime
2}\right)\left(1-\sigma^{\prime}-\mathcal{E}^{\prime}\right)^{2}\right]\\\
+\frac{1}{2}X_{\phi}^{-1}\left[\frac{1}{4}\phi^{\prime
4}-\frac{1}{2}\left(1-\frac{T}{T_{C}}\right)\phi^{\prime
2}\right]\end{gathered}\;.$ (14)

Note that there are three independent energy scales $X_{\phi}^{-1}$,
$\frac{P_{s}^{2}}{\epsilon}$ and $A$, corresponding to the tilts, the polar
discontinuity, and the biquadratic coupling, respectively.

Setting $T=0$ and minimizing Eq. (14) with respect to $\sigma^{\prime}$ and
$\phi^{\prime}$, we get

$\begin{split}\phi^{\prime}&=\sqrt{1-2AX_{\phi}(1-\sigma^{\prime}-\mathcal{E}^{\prime})^{2}}\\\
\sigma^{\prime}&=1-\mathcal{E}^{\prime}-\frac{d_{\text{c}}/d}{1+\frac{A}{(P_{s}^{2}/\epsilon)}\phi^{\prime
2}}\end{split}\;.$ (15)

Note that when we set $A=0$, we recover the solutions of the uncoupled order
parameters: Eq. (6) for $\sigma^{\prime}$ and $\phi^{\prime}=1$. Eq. (15) is a
pair of coupled equations which does not have analytic solutions. However, we
can use physical constraints to understand the behavior of $\sigma^{\prime}$
and $\phi^{\prime}$ upon coupling. Firstly, the square of the tilts cannot be
negative:

$\begin{gathered}\phi^{\prime
2}=1-2AX_{\phi}(1-\sigma^{\prime}-\mathcal{E}^{\prime})^{2}\geq 0\\\ \implies
1-\mathcal{E}^{\prime}-\frac{1}{\sqrt{2AX_{\phi}}}\leq\sigma^{\prime}\leq
1-\mathcal{E}^{\prime}+\frac{1}{\sqrt{2AX_{\phi}}}\end{gathered}\;.$ (16)

Note that Eq. (16) provides bounds on $\mathcal{E}^{\prime}$, since the upper
bound must be greater than or equal to 0, and the lower bound must be less
than or equal to 1:

$-\frac{1}{\sqrt{2AX_{\phi}}}\leq\mathcal{E}^{\prime}\leq
1+\frac{1}{\sqrt{2AX_{\phi}}}\;,$ (17)

although the upper bound is redundant because
${\sigma^{\prime}(\mathcal{E}^{\prime}>1)<0}$, which is unphysical. The lower
bound in Eq. (16) can also lead to situations where $\sigma^{\prime}<0$. This
leads to a change in behaviour when the following condition relating the
coefficients $X_{\phi}^{-1}$ and $A$ is held:

$A=\frac{1}{2(1-\mathcal{E}^{\prime})^{2}}X_{\phi}^{-1}\;.$ (18)

Above this line, for $\sigma^{\prime}=0$ we would have $\phi^{2}<0$, so the
tilts must be zero in the absence of carriers. Thus, the carriers will appear
at $d_{\text{c}}$ as in the uncoupled model, and the tilts will appear at some
thickness, which we call $d_{\phi}$, which is greater than $d_{\text{c}}$.
Below this line, there are finite tilts in the absence of carriers. In this
case, carriers appear earlier at $d_{\phi}<d_{\text{c}}$, facilitated by a
change in the tilts.

To summarize, above the straight line in Eq. (18), there are two
possibilities: an uncoupled appearance of carriers at $d_{\text{c}}$ and a
kink in the carrier concentration facilitated by the appearance of tilts at
$d_{\phi}>d_{\text{c}}$. Below the line, both order parameters appear
simultaneously at $d_{\phi}<d_{\text{c}}$, i.e. there is one critical
thickness which is reduced by tilting.

While the uncoupled carrier transition at $d_{\text{c}}$ is always continuous,
we can investigate whether the transition at $d_{\phi}$ is continuous or
discontinuous. If we insert the solution for $\phi^{\prime}$ into
$\sigma^{\prime}$ in Eq. (15), we obtain a cubic equation,

$\begin{split}f(\sigma^{\prime})=2AX_{\phi}(1-\sigma^{\prime}-\mathcal{E}^{\prime})^{3}&-\left(1+\frac{P_{s}^{2}/\epsilon}{A}\right)(1-\sigma^{\prime}-\mathcal{E}^{\prime})\\\
&+\left(\frac{P_{s}^{2}/\epsilon}{A}\right)\frac{d_{\text{c}}}{d}=0\end{split}\;,$
(19)

the roots of which determine the value of $\sigma^{\prime}$. We cannot obtain
the roots analytically, but we can plot $f(\sigma^{\prime})$ at different
thicknesses and gain some insight about the transition, see Fig. 3.

Figure 3: Sketch of Eq. (19) for different values of thickness. The black dots
mark the maxima/minima of the curve, i.e. Eq. (20). A non-zero carrier
concentration is obtained when the polynomial has a root between 0 and 1,
indicated by the red dots. For the highest curve, there is no root between 0
and 1, so the $\sigma=0$. When $\sigma^{\prime}_{\text{min}}$ crosses the zero
axis, we will have $\sigma>0$. Since $\sigma^{\prime}_{\text{min}}>0$,
$\sigma$ jumps from zero to a positive value, and the transition is
discontinuous. As the curve is shifted further downwards, the root moves to
the right until it reaches 1, where it saturates.

Since $d$ only appears in the constant term, the curve is shifted up or down
by changing thicknes, and the positions of the extrema are unaffected. The
maximum and minimum are uniquely determined by the energy densities
$X_{\phi}^{-1}$, $A$ $\frac{P_{s}^{2}}{\epsilon}$, and the field strength
$\mathcal{E}^{\prime}$:

$\sigma^{\prime}_{\text{min/max}}=1-\mathcal{E}^{\prime}\pm\sqrt{\frac{1}{6AX_{\phi}}\left(1+\frac{P_{s}^{2}/\epsilon}{A}\right)}\;.$
(20)

Again, Eq. (20) leads to bounds on the applied field:

$-\sqrt{\frac{1}{6AX_{\phi}}\left(1+\frac{P_{s}^{2}/\epsilon}{A}\right)}\leq\mathcal{E}^{\prime}\leq
1+\sqrt{\frac{1}{6AX_{\phi}}\left(1+\frac{P_{s}^{2}/\epsilon}{A}\right)}\;,$
(21)

and again the upper bound is redundant. Eq. (20) provides a clear picture of
how the carrier transition occurs: the curve in Eq. (19) is shifted by
changing $d$, and carriers appear when there is a root between 0 and 1. If
$\sigma^{\prime}_{\text{min}}<0$, the carriers appear continuously and if
$0<\sigma^{\prime}_{\text{min}}<1$ they appear discontinuously. Thus,
$\sigma^{\prime}_{\text{min}}=0$ defines a boundary between first and second
order transitions:

$6(1-\mathcal{E}^{\prime})^{2}A^{2}-X_{\phi}^{-1}\left(A+\frac{P_{s}^{2}}{\epsilon}\right)=0\;,$
(22)

which is quadratic in $A$ and has one positive solution:

$A=\frac{1}{12(1-\mathcal{E}^{\prime})^{2}}X_{\phi}^{-1}\left(1+\sqrt{1+24(1-\mathcal{E}^{\prime})^{2}\frac{(P_{s}^{2}/\epsilon)}{X_{\phi}^{-1}}}\right)\;.$
(23)

Above this curve, the transition at $d_{\phi}$ is discontinuous and below it
is continuous. Combining Eqs. (18) and (23) results in four different possible
sequences of transitions, depending only on the energy densities
$X_{\phi}^{-1}$, $\frac{P_{s}^{2}}{\epsilon}$, $A$ and the field
$\mathcal{E}^{\prime}$.

###  II.3 Transitions with thickness: zero field

At zero field, Eqs. (18) and (23) become

$\begin{gathered}A=\frac{1}{2}X_{\phi}^{-1}\\\
A=\frac{1}{12}X_{\phi}^{-1}\left(1+\sqrt{1+24\frac{(P_{s}^{2}/\epsilon)}{X_{\phi}^{-1}}}\right)\end{gathered}\;,$
(24)

and the appearance of carriers and tilts depends only on the energy densities
$X_{\phi}^{-1}$, $\frac{P_{s}^{2}}{\epsilon}$, $A$. These transitions are
summarized in Fig. 4. The two lines in Eq. (24) intersect at a tetra-critical
point:

$\begin{split}A^{*}&=\frac{P_{s}^{2}}{2\epsilon}\\\
{X_{\phi}^{-1}}^{*}&=\frac{P_{s}^{2}}{\epsilon}\\\ \end{split}\;.$ (25)

In Fig. 5 we plot the order parameters as a function of thickness in each
region of Fig. 4, as well as contours of the total free energy as a function
of $\sigma^{\prime}$ for several thicknesses.

Figure 4: Diagram summarizing the coupled transitions of carriers and tilts.
Region I, the line $A=0$, describes the uncoupled order parameters. The tilts
and carriers appear simulataneously in regions II (continuous) and IV
(discontinuous). The tilts appear after the carriers in regions III
(continuous) and V (discontinuous). First-principles calculations, summarized
in Table 1, were used to place LAO/STO on the diagram. The smaller dots show
the effect of applying a small amount of biaxial strain. Figure 5: Summary of
the carrier transitions with film thickness for the four coupled regions in
Fig. 4. The left panel for each region is an order parameter plot of
$\phi^{\prime}$ and $\sigma^{\prime}$ with thickness. The critical thickness
from the uncoupled theory is $d_{\text{c}}$, and the additional critical
thickness at which the tilts switch on/jump is $d_{\phi}$. The right panel
shows an effective theory, $\mathcal{F}_{\text{eff}}(\sigma^{\prime})$,
obtained by using Eq. (15) to write Eq. (14) solely in terms of
$\sigma^{\prime}$. Several energy curves are shown for different thicknesses
in each panel, and the dashed line shows a continuation of the uncoupled
theory. The value of $\sigma^{\prime}$ at each value of $d$ is given by the
horizontal position of the minimum on each curve.

###  II.4 Transitions with thickness: nonzero field

At a finite applied field, the different possible transitions with thickness
are determined by the lines

$\begin{gathered}A=\frac{1}{2(1-\mathcal{E}^{\prime})^{2}}X_{\phi}^{-1}\\\
A=\frac{1}{12(1-\mathcal{E}^{\prime})^{2}}X_{\phi}^{-1}\left(1+\sqrt{1+24(1-\mathcal{E}^{\prime})^{2}\frac{(P_{s}^{2}/\epsilon)}{X_{\phi}^{-1}}}\right)\end{gathered}\;,$
(26)

with a tetra-critical point

$\begin{gathered}A^{*}=\frac{1}{2}\frac{P_{s}^{2}}{\epsilon}\\\
{X_{\phi}^{-1}}^{*}=(1-\mathcal{E}^{\prime})^{2}\frac{P_{s}^{2}}{\epsilon}\end{gathered}\;.$
(27)

Interestingly, Eq. (27) reduces to Eq. (25) when $\mathcal{E}^{\prime}=0,+2$,
although $\mathcal{E}^{\prime}=+2$ is not physical. A positive field moves the
tetra-critical point to the left, making region II the largest region as
$\mathcal{E}^{\prime}\to 1$. A negative field moves the tetra-critical point
to the right. Thus, using a negative or positive field, the boundary lines can
be moved, deforming the regions, and it may be possible to tune the order of
the carrier and tilt transitions with thickness at a polar-nonpolar perovskite
interface. The phase transition diagram is sketched in Fig. 6 for
$\mathcal{E}^{\prime}=-0.2$. We can see that the tetra-critical point moves to
the right, changing the positions of the regions.

Figure 6: Phase transition lines at $\mathcal{E}^{\prime}=-0.2$. A positive
field makes them move up and a negative field makes them move down. So it
might be possible to tune the transitions with thickness by applying a fixed
field. The dashed lines show the boundaries between regions at zero field.

###  II.5 Transitions with temperature

In this section we investigate the possibility of inducing a carrier
transition with temperature. When the order parameters are uncoupled, at a
nonzero temperature we have

$\begin{split}\sigma^{\prime}(d)&=1-\frac{d_{\text{c}}}{d}\\\
\phi^{\prime}(T)&=\sqrt{1-\frac{T}{T_{C}}}\end{split}\;,$ (28)

where we set $\mathcal{E}^{\prime}=0$ for simplicity. For nonzero $A$,
minimizing Eq. (14) gives

$\begin{split}\phi^{\prime}&=\sqrt{\left(1-\frac{T}{T_{C}}\right)-2AX_{\phi}(1-\sigma^{\prime})^{2}}\\\
\sigma^{\prime}&=1-\frac{1}{1+\frac{A}{(P_{s}^{2}/\epsilon)}\phi^{\prime
2}}\frac{d_{\text{c}}}{d}\end{split}\;.$ (29)

If we require that $\phi^{2}\geq 0$ as before, we get

$1-\frac{1}{\sqrt{2AX_{\phi}(T)}}\leq\sigma^{\prime}\leq
1+\frac{1}{\sqrt{2AX_{\phi}(T)}}\;,$ (30)

where $X_{\phi}(T)=\left(1-\frac{T}{T_{C}}\right)^{-1}X_{\phi}$. As in Eq.
(16), the inequality on the left hand side can lead to situations where
$\phi^{2}<0$. This implies that there is an additional critical temperature,
$T_{C}^{\prime}$:

$T_{C}^{\prime}\equiv T_{C}\left(1-2AX_{\phi}\right)\;.$ (31)

This is the reduced critical temperature of the thin film, a phenomenon which
is typically seen in phase transitions when going from bulk to thin films
Lubensky and Rubin (1975); Kretschmer and Binder (1979); De Gennes (2018);
Chandra and Littlewood (2007). Interestingly, we have $T_{C}^{\prime}=0$ when

$A>\frac{1}{2}X_{\phi}^{-1}\;,$ (32)

i.e. there is no tilt transition at all. This is the same straight line
obtained for transitions with thickness and applied field at zero temperature.
However, it has a different meaning for transitions with temperature: below
the line there is a tilt transition at a reduced critical temperature
$T_{C}^{\prime}$, and above the line there is no tilt transition. There is a
second transition temperature, $T_{C}^{\prime\prime}$, below which the
carriers appear. Inserting $\phi^{\prime}$ into $\sigma^{\prime}$ in Eq. (29)
and letting $\sigma^{\prime}(T_{C}^{\prime\prime})=0$, we get

$T_{C}^{\prime\prime}\equiv
T_{C}^{\prime}-T_{C}\left(\frac{P_{s}^{2}/\epsilon}{A}\right)\left(\frac{d_{\text{c}}}{d}-1\right)$
(33)

We can see that $T_{C}^{\prime\prime}<T_{C}^{\prime}$ for sub-critical films,
and $T_{C}^{\prime\prime}=T_{C}^{\prime}$ when $d=d_{c}$. Thus, two
transitions are possible by decreasing temperature in sub-critical films and
only one transition is possible otherwise. Additionally, it is possible for
the transition at $T_{C}^{\prime\prime}$ to be first or second order.
Inserting $\phi^{\prime}$ into $\sigma^{\prime}$ as before, we obtain another
cubic equation in $\sigma^{\prime}$:

$\begin{split}f_{T}(\sigma^{\prime})=\underbrace{2AX_{\phi}}_{\Lambda}(1-\sigma^{\prime})^{3}&-\underbrace{\left(1-\frac{T}{T_{C}}+\frac{P_{s}^{2}/\epsilon}{A}\right)}_{\Gamma(T)}(1-\sigma^{\prime})\\\
&+\left(\frac{P_{s}^{2}/\epsilon}{A}\right)\frac{d_{\text{c}}}{d}=0\end{split}\;,$
(34)

with extrema

$\sigma^{\prime}_{\text{max/min}}(T)=1\pm\sqrt{\frac{\Gamma(T)}{3\Lambda}}\;.$
(35)

Figure 7: Phase transition diagram for transitions driven by temperature. The
line $A=0$ describes the uncoupled order parameters (region I). The straight
line is the same as in Fig. 4. Above this line, no tilt transitions with
temperature are possible: $T_{C}^{\prime}=0$. The curved line is given by Eq.
(38) and determines whether the carrier transition is continuous (below) or
discontinuous (above). The number of transitions, i.e. whether or not
$T_{C}^{\prime}$ and $T_{C}^{\prime\prime}$ are different is determined by the
ratio $\frac{d}{d_{\text{c}}}$ only: if $\frac{d}{d_{\text{c}}}<1$ there are
two transitions and if $\frac{d}{d_{\text{c}}}\geq 1$ there is only one
transition. Figure 8: Summary of the carrier transitions with temperature for
the two coupled regions in Fig. 7. The left hand side shows transitions for
thin films at the uncoupled critical thickness ($d=d_{\text{c}}$) and the
right hand side shows transitions for sub-critical thin films
($d=\frac{1}{2}d_{\text{c}}$). The right panels shows energy contours for
different temperatures and are as described in Fig. 5, but an effective theory
in $\phi^{\prime}$ instead of $\sigma^{\prime}$ is used.

In the previous section, the positions of the extrema were not dependent on
the quantity that induced the transitions. Assuming the transition occurs at
some temperature $T^{*}$, we can determine the character of the transition by
requiring that $f_{T}(\sigma_{\text{min}}(T^{*}))=0$,

$\begin{split}-\Lambda\left(\sqrt{\frac{\Gamma(T^{*})}{3\Lambda}}\right)^{3}&+\Gamma(T^{*})\left(\sqrt{\frac{\Gamma(T^{*})}{3\Lambda}}\right)\left(\frac{P_{s}^{2}/\epsilon}{A}\right)\frac{d_{\text{c}}}{d}=0\\\
\implies&\Gamma(T^{*})=\sqrt[3]{\frac{27}{4}\left(\left(\frac{P_{s}^{2}/\epsilon}{A}\right)\frac{d_{c}}{d}\right)^{2}\Lambda}\end{split}\;,$
(36)

and also requiring that $\sigma_{\text{min}}(T^{*})>0$,

$\frac{1}{4}\left(\frac{P_{s}^{2}/\epsilon}{A}\right)\frac{d_{\text{c}}}{d}<\Lambda\;.$
(37)

Thus, we get

$A>\sqrt{\frac{1}{4}\frac{d_{\text{c}}}{d}\left(\frac{P_{s}^{2}}{\epsilon}\right)X_{\phi}^{-1}}$
(38)

The lines in Eqs. (32) and (38) form a diagram for transitions with
temperature, shown in Fig. 7. There are only two distinct regions (II and
III), but the number of transitions in each is determined by the ratio
$\frac{d}{d_{\text{c}}}$, so four different types of carrier transitions are
possible. Plots of the order parameters for the four different scenarios are
shown in Fig. 8.

###  II.6 Transitions with applied field

For transitions with applied field, note from the left-hand side of Eq. (16)
that if $\sigma=0$, there is a critical field value at which the tilts switch
on or off:

$\mathcal{E}_{\phi}=1-\frac{1}{\sqrt{2AX_{\phi}}}\geq 0$ (39)

We get $\mathcal{E}_{\phi}=0$ when $A=\frac{1}{2}X_{\phi}^{-1}$, which is the
same condition for transitions with thickness at zero field. In order to
determine whether transitions with applied field are continuous or
discontinuous, we must solve Eq. (19), assuming there is some critical field
at which carriers appear, i.e. we must solve
$f(\sigma^{\prime}_{\text{min}}(\mathcal{E}^{\prime*}))=0$ for
$\mathcal{E}^{\prime*}$ in a similar manner to the previous section. The
boundary between continuous and discontinuous transitions of the tilts with
applied field is:

$X_{\phi}^{-1}=\frac{27}{2}\frac{\left(d_{\text{c}}/d\right)^{2}}{\left(1+\frac{P_{s}^{2}/\epsilon}{A}\right)^{3}}\frac{1}{A}\;,$
(40)

which in general depends on the ratio $\frac{d}{d_{\text{c}}}$. For
$\frac{d}{d_{\text{c}}}=1$, Eq. (40) touches the straight line at the point
$(A^{*},X_{\phi}^{-1*})$. When $\frac{d}{d_{\text{c}}}<1$, the two lines
intersect before $(A^{*},X_{\phi}^{-1*})$. When $\frac{d}{d_{\text{c}}}>1$,
the two lines never intersect. The diagram describing the different possible
transitions with applied field is sketched in Fig. 9 for
$\frac{d}{d_{\text{c}}}=1$, and plots of the order parameters for each region
are shown in Fig. 10. In all except region II, the tilts are zero at zero
field. When a negative field is applied, the carriers appear as in the
uncoupled model. In region II the carriers grow linearly with a negative
applied field, and can be switched off with a positive applied field. In all
regions, tilts appear when a positive field is applied, because a positive
field suppresses the polar mode. Although first-order transitions of the tilts
should be observed in regions III and IV, we did not find any first-order
transitions for physically meaningful values of the order parameters or
applied field in region III. A first-order transition of the tilts and
carriers for a positive applied field was observed in region IV, however.

Figure 9: Diagram for transitions with applied field. The curved line depends
on the ratio $\frac{d}{d_{\text{c}}}$. The dashed black line is
$\frac{d}{d_{\text{c}}}=1$, which touches the straight line at
$(A^{*},X_{\phi}^{-1*})$. Curves for $\frac{d}{d_{\text{c}}}=0.8,1.2$ are also
sketched. First-order transitions should be observed below the curved line but
above the straight line, i.e. regions III and IV. Figure 10: Summary of the
carrier transitions with applied applied field for the four coupled regions in
Fig. 9. The dashed blue line shows the value of $\sigma^{\prime}$ without
coupling to tilts. All plots are for $\frac{d}{d_{\text{c}}}=1$, and the
following coorindates in Fig. 9: (a): $X_{\phi}^{-1}$ = 1, $A$ = 0.25, (b):
$X_{\phi}^{-1}$ = 0.25, $A$ = 0.15, (c): $X_{\phi}^{-1}$ = 1.75, $A$ = 1, (d):
$X_{\phi}^{-1}$ = 1, $A$ = 1.

##  III Coupling to inhomogeneous tilts

Figure 11: Sketch of inhomogeneous tilts in a thin film grown on a substrate,
e.g. LAO/STO, and the effect of (a): $\lambda$, the correlation length of the
tilts in the LAO film (b): $\delta_{\text{surface}}$, the relative difference
in the energy of the tilts at the surface with respect to the bulk and (c):
$\delta_{\text{interface}}$, the relative difference in the energy of the
tilts at the interface with respect to the bulk. The effect of the energetics
of the tilts in the substrate may be absorbed into
$\delta_{\text{interface}}$.

We can allow for inhomogeneous tilts by generalizing from a Landau theory to a
Ginzburg-Landau theory, i.e. the free energy is expanded in powers of both the
order parameter and its gradient. We do this following the methodology which
has been used to describe semi-infinite ferromagnetic Lubensky and Rubin
(1975), ferroelectric Kretschmer and Binder (1979); Chandra and Littlewood
(2007); Luk’yanchuk  _et al._ (2009) and superconducting De Gennes (2018)
systems in which the order parameter is inhomogeneous due to the presence of a
free surface. Considering only the tilts in the thin film, Eq. (12) becomes an
integral over the volume of the film:

${\mathcal{F}^{\prime}}_{\phi}^{\text{GL}}=\frac{1}{V}\int_{V}\left[\lambda^{2}\left(\nabla\phi^{\prime}(\bm{r})\right)^{2}+\frac{1}{4}\phi^{\prime}(\bm{r})^{4}-\frac{1}{2}\phi^{\prime}(\bm{r})^{2}\right]\mathop{}\\!\mathrm{d}{V}\;,$
(41)

where ${\mathcal{F}^{\prime}}\equiv\mathcal{F}/\frac{1}{2}X_{\phi}^{-1}$ and
$\lambda>0$ is the correlation length of the tilts, which sets the length
scale for the variations of the tilts, see Fig. 11 (a). The correlation length
depends in general on the temperature of the system, but here we assume a
fixed value for simplicity. Eq. (41) is written so that each term in the
integrand is dimensionless, and the free energy is in units of
$\frac{1}{2}X_{\phi}^{-1}$, the energy scale of the homogeneous tilts.

Next, we reintroduce the carriers. We take the order parameter to be the
carrier concentration averaged over an in-plane unit cell:
$\sigma^{\prime}\equiv\left<\sigma^{\prime}\right>=\frac{1}{A}\int\sigma^{\prime}(x,y)\mathop{}\\!\mathrm{d}{A}$.
If we assume that the in-plane variance in $\sigma^{\prime}$ is small, so that
$\left<(1-\sigma^{\prime})\right>^{2}\approx(1-\left<\sigma^{\prime}\right>)^{2}$,
then
${\mathcal{F}^{\prime}}_{\sigma}^{\text{GL}}={\mathcal{F}^{\prime}}_{\sigma}$.
It should be noted that this is a significant assumption, especially in the
low carrier regime.

For inhomogeneous tilts, the biquadratic coupling term becomes

${\mathcal{F}^{\prime}}_{\sigma\phi}^{\text{GL}}=\frac{1}{V}\int\frac{1}{2}(2AX_{\phi})\phi^{\prime}(\textbf{r})^{2}(1-\sigma^{\prime})^{2}\mathop{}\\!\mathrm{d}{V}\;,$
(42)

which, since we assume the in-plane variance in $\sigma^{\prime}$ is small,
reduces to

${\mathcal{F}^{\prime}}^{\text{GL}}_{\sigma\phi}=\frac{1}{2}(2AX_{\phi})(1-\sigma^{\prime})^{2}\left<\phi^{\prime
2}\right>\;,$ (43)

where

$\left<\phi^{\prime
2}\right>=\frac{1}{V}\int\phi^{\prime}(\textbf{r})^{2}\mathop{}\\!\mathrm{d}{V}$
(44)

is the mean of the square of the tilt throughout the film.

Under these assumptions, the free energy can be minimized with respect to
$\sigma^{\prime}$ as in the homogeneous case, yielding

$\sigma^{\prime}=1-\frac{1}{1+\frac{A}{P_{0}^{2}/\epsilon}\left<\phi^{\prime
2}\right>}\left(\frac{d_{\text{c}}}{d}+\mathcal{E}^{\prime}\right)\;,$ (45)

which is the same as the second line of Eq. (15), but with ${\phi^{\prime
2}\to\left<\phi^{\prime 2}\right>}$. We only consider transitions with
thickness for illustrative purposes, although it would be simple enough to
include an applied field as we did in the previous section.

For the tilts, we have

$\begin{split}{\mathcal{F}^{\prime}}_{\phi}^{\text{GL}}+{\mathcal{F}^{\prime}}_{\sigma\phi}^{\text{GL}}=&\frac{1}{V}\int_{V}\left[\lambda^{2}\left(\nabla\phi^{\prime}\right)^{2}+\frac{1}{4}\phi^{\prime
4}\right.\\\
&\left.\qquad\quad-\frac{1}{2}\left(1-2AX_{\phi}(1-\sigma^{\prime})^{2}\right)\phi^{\prime
2}\right]\mathop{}\\!\mathrm{d}{V}\end{split}\;.$ (46)

Following a similar treatment of ferroelectric thin films with inhomogeneous
polarization Chandra and Littlewood (2007), we split the free energy into
interior and surface contributions:
${\mathcal{F}^{\text{GL}}=\mathcal{F}^{\text{GL}}_{\text{interior}}+\mathcal{F}^{\text{GL}}_{\text{boundary}}}$,
where $\mathcal{F}^{\text{GL}}_{\phi,\text{boundary}}$ is given in terms of an
expansion of a local order parameter at the boundaries. First, we use the
reverse product rule to rewrite the gradient term:

$(\nabla\phi^{\prime})^{2}=\nabla\cdot(\phi^{\prime}\nabla\phi^{\prime})-\phi^{\prime}\nabla^{2}\phi^{\prime}\;,$
(47)

then we use Stokes’ theorem on the first term in Eq. (47):

$\begin{split}\int(\nabla\phi^{\prime})^{2}\mathop{}\\!\mathrm{d}V&=\int\left[\nabla\cdot(\phi^{\prime}\nabla\phi^{\prime})-\phi^{\prime}\nabla^{2}\phi^{\prime}\right]\mathop{}\\!\mathrm{d}V\\\
&=\int\phi^{\prime}\nabla\phi^{\prime}\mathop{}\\!\mathrm{d}S-\int\phi^{\prime}\nabla^{2}\phi^{\prime}\mathop{}\\!\mathrm{d}V\end{split}\;,$
(48)

where for a thin film on a substrate the surface integral has two terms: one
at the free surface and one at the interface with the substrate:

$\begin{gathered}{\mathcal{F}^{\prime}}^{\text{GL}}=\frac{1}{V}\int_{V}\bigg{[}\frac{1}{4}\phi^{\prime
4}-\frac{1}{2}\left(1-2AX_{\phi}(1-\sigma^{\prime})^{2}\right)\phi^{\prime
2}\\\
-\lambda^{2}\phi^{\prime}\nabla^{2}\phi^{\prime}\bigg{]}\mathop{}\\!\mathrm{d}{V}\\\
+\frac{1}{V}\int_{\text{s}}\left[\lambda^{2}\left(\hat{n}\cdot\nabla\phi^{\prime}\right)\phi^{\prime}\right]\mathop{}\\!\mathrm{d}{S}\\\
+\frac{1}{V}\int_{\text{i}}\left[\lambda^{2}\left(\hat{n}\cdot\nabla\phi^{\prime}\right)\phi^{\prime}\right]\mathop{}\\!\mathrm{d}{S}\end{gathered}\;,$
(49)

where s and i refer to the surface and interface, respectively. Note that the
surface free energies are also in units of $\frac{1}{2}X_{\phi}^{-1}$, and the
surface integrals are dimensionless.

Next we claim that there must be a difference in the energies at the
boundaries compared to the interior of the film Lubensky and Rubin (1975); De
Gennes (2018); Chandra and Littlewood (2007), which can be described using an
expansion of the local order parameters at the boundaries:

$\begin{split}{\mathcal{F}^{\prime}}^{\text{GL}}_{\phi,\text{boundary}}&=\frac{1}{V}\int_{\text{s}}\left[\lambda^{2}\left(\hat{n}\cdot\nabla\phi^{\prime}\right)\phi^{\prime}+\frac{1}{2}\delta_{\text{s}}\phi^{\prime
2}\right]\mathop{}\\!\mathrm{d}{S}\\\
&+\frac{1}{V}\int_{\text{i}}\left[\lambda^{2}\left(\hat{n}\cdot\nabla\phi^{\prime}\right)\phi^{\prime}+\frac{1}{2}\delta_{\text{i}}\phi^{\prime
2}\right]\mathop{}\\!\mathrm{d}{S}\end{split}\;.$ (50)

The quadratic coefficient of the surface order parameter describes the
relative difference between the energy of the tilts at the surface and the
energy per unit volume of the tilts in the bulk, and therefore must have units
of length. We call this the extrapolation length, $\delta$.
$\delta_{\text{s}}$ and $\delta_{\text{i}}$ are the surface and interface
extrapolation lengths, respectively, which describe the difference in energy
at the surface and interface with respect to the bulk, see Figs. 11 (b) and
(c). In general, higher order terms could be included in this expansion, but
we assume that the leading term is the quadratic one Chandra and Littlewood
(2007). Also, the boundaries may have a temperature dependence which differs
from the bulk, which could lead to a change in sign of the extrapolation
lengths.

Minimizing the total free energy in the bulk and at both boundaries, we get

$\begin{gathered}\phi^{\prime
3}-\left(1-2AX_{\phi}(1-\sigma^{\prime})^{2}\right)\phi^{\prime}-\lambda^{2}\nabla^{2}\phi^{\prime}=0\\\
(\hat{n}\cdot(\nabla\phi^{\prime}))+\frac{\delta_{\text{i}}}{\lambda^{2}}\phi^{\prime}=0,\qquad
z=0\\\
(\hat{n}\cdot(\nabla\phi^{\prime}))+\frac{\delta_{\text{s}}}{\lambda^{2}}\phi^{\prime}=0,\qquad
z=d\\\
\sigma^{\prime}=1-\frac{1}{1+\frac{A}{P_{0}^{2}/\epsilon}\left<\phi^{\prime
2}\right>}\left(\frac{d_{\text{c}}}{d}+\mathcal{E}^{\prime}\right)\end{gathered}\;.$
(51)

This is the generalization of Eq. (15) but with inhomogeneous tilts. The
expression for $\sigma^{\prime}$ is almost identical, but with the square of
the tilt replaced with the mean of the square of the tilt throughout the film.
$\phi^{\prime}$ is described by a second order differential equation plus two
boundary conditions, one for the free surface and one for the interface. These
boundary conditions are referred to as Robin boundary conditions, i.e. a
linear combination of the tilt and its gradient at each boundary. The
equations for $\phi^{\prime}$ and $\sigma^{\prime}$ are again self-consistent,
making them difficult to solve. In the previous section we used physical
constraints to determine the character of the carrier and tilt transitions
under various conditions, summarized in Fig. 4, but this would be difficult or
impossible for Eq. (51). Thus, our only option is obtaining numerical
solutions, solving for $\sigma^{\prime}$ and $\phi^{\prime}$ iteratively until
both are converged below a suitable tolerance. This approach was used in the
previous section to obtain the order parameter plots shown in Figs. 5, 8 and
10.

Eq. (51) is solved as follows: the tilt profile inside the film is obtained by
solving the first three lines using finite difference methods. The mean of the
tilt squared in the film is then used to calculate the carrier concentration,
which is re-inserted into the ODE. The two are solved self-consistently until
their relative changes between successive iterations are below a suitable
tolerance.

For a given set of system parameters, we can obtain order parameter plots as a
function of thickness or applied field, similar to Figs. 5, 8 and 10. The main
difference with the inhomogeneous tilts is that $\lambda$, $\delta_{\text{s}}$
and $\delta_{\text{i}}$ need to be specified, in addition to $X_{\phi}^{-1}$,
$A$ and $d_{\text{c}}$. Thus, for inhomogeneous tilts, Fig. 4 becomes five
dimensional, and systemically navigating such a diagram would be impractical.

It is also difficult to estimate $\lambda$, $\delta_{\text{s}}$ and
$\delta_{\text{i}}$. Previous studies have fixed $\lambda$ to a small number
of unit cells Kretschmer and Binder (1979). In Ref. Gazquez  _et al._ , 2017,
$\lambda$ is estimated to be around 1-2 unit cells in bulk LAO from
calculations of the dispersion of the phonon branch corresponding to the
$a^{-}a^{-}c^{0}$ tilt. To our knowledge, it is not possible to directly
measure the extrapolation lengths. In principle, one could fit measurements of
tilts in a thin film, either from experiment or first-principles calculations,
to Eq. (51), but the tilts at the boundaries are strongly affected by both the
correlation length and the extrapolation length, which would lead to
overfitting.

It is helpful to note that, since we only consider indirect coupling, the
carrier concentration depends only on the mean of the square of the tilts, and
not on the specific shape of the tilts. Thus, we can understand the influence
of the inhomogeneity of the tilts by fixing the extrapolation lengths and
allowing the correlation length to change.

Figure 12: Order parameter plots in the four coupled regions in Fig. 4, for
both homogeneous (dashed lines) and inhomogeneous (solid lines) tilts. (a):
Region II (red): $X_{\phi}^{-1}=1.36454$, $A=0.407$. (b): Region III (green):
$X_{\phi}^{-1}=0.2$, $A=0.5$. $\sigma$ and $\sigma_{\text{GL}}$ are identical
between $d_{\text{c}}$ and $d_{\phi}$. (c): Region IV (yellow):
$X_{\phi}^{-1}=3$, $A=1.4$. (d): Region V (blue): $X_{\phi}^{-1}=0.5$,
$A=1.5$. $\sigma$ and $\sigma_{\text{GL}}$ are identical between
$d_{\text{c}}$ and $d_{\phi}$. In all regions, $d_{\text{c}}=2$ unit cells.

First, we examine the behavior of the order parameters in the different
regions in Fig. 4 and compare with the results obtained for homogeneous
tilting. For each region, we set $\lambda=2$ unit cells, $d_{\text{c}}=2$ unit
cells, the thickness at which a density of trapped Ti 3$d$-like states has
been observed in LAO/STO Sing  _et al._ (2009); Berner  _et al._ (2010);
Takizawa  _et al._ (2011); Slooten  _et al._ (2013), $\delta_{\text{i}}=5$
unit cells and $\delta_{\text{s}}=0$. The order parameter plots are shown in
12 (solid lines), alongside the corresponding result from the homogeneous
theory (dashed lines). We see that, for this choice of parameters, all of the
points except the one in region V shift to region III when the tilts become
inhomogeneous. In all cases, the carriers appear continuously at
$d_{\text{c}}$, and $d_{\phi}$ increases with respect to the corresponding
scenario from the homogeneous theory. This makes sense, since at smaller
thicknesses there would be a large gradient in the tilts and hence a large
energy penalty. Thus, the tilts can only appear at larger thicknesses, where
the magnitude is roughly constant inside the bulk of the film and only changes
in a small region near both boundaries.

Figure 13: Order parameter plots in region II for several values of $\lambda$,
for both homogeneous (dashed lines) and inhomogeneous (solid lines) tilts.
(a): $\lambda=2$ unit cells, (b): $\lambda=1$ unit cell, (c): $\lambda=0.5$
unit cells, and (d): $\lambda=0.5$ unit cells. Figure 14: Order parameter
plots in region IV for several values of $\lambda$, for both homogeneous
(dashed lines) and inhomogeneous (solid lines) tilts. (a): $\lambda=2$ unit
cells, (b): $\lambda=1$ unit cell, (c): $\lambda=0.5$ unit cells, and (d):
$\lambda=0.5$ unit cells.

Next, we examine the effect of changing $\lambda$. We found that for regions
III and V, the transitions remain in regions III and V, but the value of
$d_{\phi}$ is sensitive to $\lambda$. For regions II and IV, the points move
through different regions as $\lambda$ decreases, shown in Figs. 13 and 14,
respectively. As $\lambda$ decreases, and there is less of a penalty for the
inhomogeneity of the tilts, $d_{\phi}$ decreases, and eventually the order of
the transitions reverse, which occurs continuously in region II (Figs. 13 (c)
and (d)) and discontinuously in region IV (Fig. 14 (d)). These types of
transitions do not correspond to any of regions in Fig. 4, so we label them
regions VI and VII. In region VI, tilts first appear continuously at
$d_{\phi}$, and then carriers at a critical thickness between $d_{\phi}$ and
$d_{\text{c}}$. In region VII, the order of the transitions with thickness is
the same, but the carriers appear discontinuously.

In Fig. 15 we summarize the transitions observed in the inhomogeneous theory.
Regions III and V remain unchanged, with the exception that $d_{\phi}$ is
renormalized by $\lambda$. For regions II and IV, the type of transition can
change, depending on the value of $\lambda$ (and the other parameters). The
different transitions that can occur by changing $\lambda$ but fixing the rest
of the parameters are shown. For $\lambda\to\infty$, the homogeneous case is
realized. As $\lambda$ decreases, the character of the transition changes. For
sufficiently small $\lambda$, two entirely new types of transitions occur,
which we label region VI (Figs. 13 (c) and (d)) and region VII (Fig. 14 (d)).

Figure 15: Summary of the additional carrier and tilt transitions which can
occur when the tilts are inhomogeneous. The sequence of transitions which can
occur in regions II and IV by changing $\lambda$ (fixing the rest of the
parameters) are shown below and above the diagram, respectively.

##  IV Discussion and Conclusion

In this paper, we introduced a phenomenological theory which considers the
influence of homogeneous tilts on the appearance of carriers at polar-nonpolar
perovskite interfaces. We show that, upon coupling, four new types of
transitions of tilts and carriers with film thickness or applied field are
possible, depending only on the energetics of the tilts, the polar
discontinuity, and the biquadratic coupling between the tilts and the polar
mode. These include the simultaneous appearance of tilts and carriers at a
reduced thickness compared to the uncoupled theory, which can be either
continuous or discontinuous, and separate transitions, with carriers appearing
first and then tilts, the second transition being either continuous or
discontinuous.

First-principles calculations of bulk LAO were performed to predict the
character of the transitions at the LAO/STO interface. For various exchange-
correlation functionals, all of the calculations predicted that at zero
temperature and for homogeneous tilts, LAO/STO lies in region II of Fig. 4,
i.e. a single continuous transition of both tilts and carriers. We also
investigated the possibility of tuning the character of the transition via
biaxial strain by changing the in-plane lattice constant by up to $\pm$1% and
repeating the calculations. We found that the position of the point in Fig. 4
can be changed by adjusting the strain, but not by enough to change the
character of the transitions. However, in general, if a system were closer to
a boundary between regions, it may be possible to change the regions using
biaxial strain. Additionally, if compressive strain is applied to
$a_{\text{LAO}}$, the tilt pattern changes to $a^{0}a^{0}c^{-}$ Hatt and
Spaldin (2010), which would change the competition with the polar mode. This
could be achieved by using a nonpolar perovskite with a smaller lattice
constant than LAO as the substrate material.

In addition to changing the strain and the tilt pattern as mentioned above,
changing the materials which form the polar-nonpolar interface can also change
the chemistry at the interface, and hence the polar discontinuity. Different
interfaces have been investigated experimentally Li  _et al._ (2015); Kumar
_et al._ (2020); Wang  _et al._ (2016); Moetakef  _et al._ (2011),
computationally Fan  _et al._ (2014); Maznichenko  _et al._ (2019); Kim  _et
al._ (2018), and high throughput searches for new interfaces have been done
Yang  _et al._ (2016). Carrier transitions have also been observed at spinel-
perovskite interfaces Chen  _et al._ (2013); Niu  _et al._ (2017); Yang  _et
al._ (2020). Thus, it may be possible to observe a wide variety of behavior in
the carrier transitions at different polar-nonpolar perovskite interfaces for
different combinations of materials. More specifically, by estimating where
different systems would lie on Fig. 4, it may be possible to establish an
atlas of transitions across different possible polar-nonpolar perovskite
interfaces.

We also showed how the appearance of carriers can be tuned via temperature
through the temperature dependence of the tilts. At a finite temperature, the
curvature of the tilt double well is
${X_{\phi}^{-1}(T)=X_{\phi}^{-1}\left(1-\frac{T}{T_{C}}\right)}$, so
increasing the temperature would move any point on the diagram in Fig. 4 to
the left. We would expect the dependence of the biquadratic coupling term
$A(T)$ on temperature to be similar but weaker. Thus, it may be possible to
move a system from one point in the diagram to another if the initial point is
sufficiently close to a boundary. In Ref. Fister  _et al._ , 2014, when grown
on a substrate of cubic STO, LAO was found to undergo a transition from
untilted to tilted below $695^{\circ}$\mathrm{C}$$ for a thickness of 24 unit
cells and $540^{\circ}$\mathrm{C}$$ for a thickness of 9 unit cells, where
both samples were under the same strain conditions. Because both films were
supercritical there was no carrier transition associated with the tilt
transitions, but their observations are in agreement with our prediction of a
second transition temperature, $T_{C}^{\prime\prime}$, which is sensitive to
the ratio $\frac{d}{d_{\text{c}}}$. Assuming that the films were sufficiently
thick so that the effect of gradients are negligible, this may provide
evidence for the additional transition temperature arising from the coupling
between tilts and carriers.

We also investigated the influence of the inhomogeneity of the tilts on the
coupling between tilts and carriers. Upon coupling to inhomogeneous tilts, the
diagram in Fig. 4 becomes five dimensional, with additional axes for the
correlation length of the tilts in the polar film, and extrapolation lengths
for the interface and the free surface of the film. Practically, it would be
very demanding to explore this five dimensional space numerically, and to our
knowledge, it can’t be done analytically. However, since the biqaudratic
coupling between tilts and carriers only depends on the mean of the tilt
squared throughout the film, we can examine the influence of inhomogeneous
tilts by examining the effect of the correlation length on transitions in each
region in Fig. 4, with the rest of the parameters fixed. We found that, for
the regions corresponding to the separate appearances of tilts and then
carriers, the order of the transitions is unaffected by the correlation
length, except for the value of the second transition shifting slightly. For
the regions corresponding to simultaneous transitions, we found that a finite
correlation length can actually change the character of the transition, and
can lead to two new types of transitions which are not possible for
homogeneous tilts, where the tilts appear before carriers, either continuously
or discontinuously.

Due to the large number of system parameters required, many of which cannot be
directly measured, it remains difficult to conclusively determine the order of
the carrier transition at polar-nonpolar perovskite interfaces. However, we
have shown that, even using a phenomenological description, an extremely rich
variety of behavior is possible for the appearance of carriers at polar-
nonpolar perovskite interfaces. This includes first-order transitions, which
have been observed experimentally Ohtomo and Hwang (2004), and have been
claimed to occur in other theoretical studies Yu and Zunger (2014), but
without a detailed explanation.

## Acknowledgments

DB would like to acknowledge funding from the EPSRC Centre for Doctoral
Training in Computational Methods for Materials Science under grant number
EP/L015552/1, St. John’s College (University of Cambridge), and the University
of Liége under special funds for research (IPD-STEMA fellowship programme).
PAP would like to acknowledge funding from the Diputación Foral de Gipuzkoa
through Grant 2020-FELL-000005-01. Funding from the Spanish MICIN is also
acknowledged, through grant PID2019-107338RB- C61/AEI/10.13039/501100011033,
as well as a María de Maeztu award to Nanogune, Grant CEX2020-001038-M funded
by MCIN/AEI/ 10.13039/501100011033.

## Appendix A: First-Principles Calculations

We can approximate $X_{\phi}^{-1}$, $A$ and $\frac{P_{s}^{2}}{\epsilon}$ for
LAO/STO by performing first-principles calculations using bulk LAO, biaxially
strained to the lattice parameter of cubic STO. We follow the methodology in
Ref. Gazquez  _et al._ , 2017, where similar calculations were performed.

First-principles density functional theory (DFT) calculations were performed
using the abinit code Gonze  _et al._ (2002); Gonze (2005); Gonze  _et al._
(2009, 2016, 2020), version 8. We used both Perdew-Burke-Edwards (PBE) Perdew
_et al._ (1996) and PBEsol Perdew  _et al._ (2008); Marques  _et al._ (2012)
exchange-correlation functionals within the generalized gradient approximation
(GGA) using psml García  _et al._ (2018) norm-conserving Hamann (2013)
pseudopotentials, obtained from pseudo-dojo Van Setten  _et al._ (2018). We
also performed calculations using the Perdew-Wang (PW92) Perdew and Wang
(1992) exchange-correlation functional within the local density approximation
(LDA), using the projector augmented-wave method Blöchl (1994); Torrent  _et
al._ (2008) (PAW) in order to compare to Ref. Gazquez  _et al._ , 2017. For
the electronic configurations, we included 11 valence electrons for La
($\mathrm{5s^{2}5p^{6}5d^{1}6s^{2}}$), 11 for Al
($\mathrm{2s^{2}2p^{6}3s^{2}3p^{1}}$) and 6 for O ($\mathrm{2s^{2}2p^{4}}$),
explicitly including the semicore $\mathrm{2s^{2}2p^{6}}$ states in the
valence configuration of Al. When the semicore electrons were included for Al,
a cutoff of $2500\ $\mathrm{e}\mathrm{V}$$ was required to adequately converge
the total energy of the 5-atom unit cell. A Monkhorst-Pack $k$-point grid
Monkhorst and Pack (1976) of $6\times 6\times 6$ was used for the 5-atom
calculations and a grid of $4\times 4\times 3$ was used for the 20-atom
calculations.

Functional |  $a_{\text{LAO}}$ (Å) |  $a_{\text{STO}}$ (Å) |  $c$ (Å) |  $E_{0}$ (meV/$\Omega_{5}$) |  $X_{\phi}^{-1}$ (meV/$\Omega_{5}$) | $\epsilon_{5}$ | $\epsilon_{20}$ |  $\frac{P_{s}^{2}}{\epsilon_{5}}$ (meV/$\Omega_{5}$) |  $A$ (meV/$\Omega_{5}$)
---|---|---|---|---|---|---|---|---|---
PAW LDA Gazquez  _et al._ (2017) | 3.75 | 3.85 | 3.68 | 27.7 | 221.6 | - | - | - | 96.1
PAW LDA | 3.75 | 3.86 | 3.67 | 35.9 | 287.8 | 26.1 | - | 425.7 | -
PBE | 3.81 | 3.93 | 3.72 | 60.5 | 483.9 | 40.4 | 26.5 | 268.9 | 141.2
PBEsol | 3.77 | 3.89 | 3.69 | 54.2 | 433.8 | 34.6 | 24.6 | 317.9 | 129.4
Table 1: Summary of results from first-principles calculations.

We first optimized the geometry of cubic STO in order to obtain the lattice
parameter, $a_{\text{STO}}$. The in-plane lattice parameters of LAO were then
fixed to $a_{\text{STO}}$, and the out-of-plane lattice parameter $c$ was
allowed to relax. It was found that the lattice parameter was not affected
much more by the tilt or the electric fields applied, so it was fixed to the
value obtained in the untilted case in the more time-consuming calculations.
The untilted calculations were performed using the 5-atom primitive cell. To
allow for tilts, we used a $\sqrt{2}\times\sqrt{2}\times 2$ supercell
containing 20 atoms, which is the smallest cell required to allow for the
$a^{-}a^{-}c^{0}$ tilt pattern which appears in LAO when biaxially strained to
$a_{\text{STO}}$ Hatt and Spaldin (2010). The depth of the tilt double well,
$E_{0}$, is the energy difference per 5 atom unit cell between the tilted and
untilted systems. Using $\mathcal{F}_{\phi}(\phi^{\prime}=\pm 1)=E_{0}$, we
get

$X_{\phi}^{-1}=-8E_{0}\;.$ (52)

$E_{0}$ and $X_{\phi}^{-1}$ are given in Table I for all of the functionals
used, in units of
$\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\Omega_{\mathrm{5}}$, where
$\Omega_{5}$ is the volume of the 5-atom unit cell. The LDA results are close
to those obtained in Ref. Gazquez  _et al._ , 2017, but the results obtained
using PBE and PBEsol differ by a factor of $\sim 2$.

Figure 16: Polarization as a function of applied electric field, from the
PBEsol calculations. The solid lines are results obtained from the untilted
system and the dashed lines are results obtained for the tilted system. The
blue lines indicate the static case, where the ions were fixed and the red
lines indicate the relaxed case, where a geometry relaxation was performed at
each value of $\mathcal{E}$.

$A$ and $\frac{P_{s}^{2}}{\epsilon}$ were obtained by calculating the
dielectric constants in the untilted and tilted systems. This was done by
performing a set of bulk calculations using a finite electric field parallel
to the $c$-axis Nunes and Vanderbilt (1994); Nunes and Gonze (2001); Souza
_et al._ (2002); Umari and Pasquarello (2002); Diéguez and Vanderbilt (2006).
In Ref. Gazquez  _et al._ , 2017 a finite displacement field was used Stengel
_et al._ (2009), but either type of field can be used to calculate $A$. The
polarization, calculated using Berry phases Berry (1984); King-Smith and
Vanderbilt (1993), was measured in the presence of a small applied field
$\mathcal{E}$, with the ions both fixed and allowed to relax. Plots of
polarization as a function of applied field are shown in Fig. 16. In Gaussian
units, the relation between polarization and electric field, and dielectric
constant is

$\begin{split}P_{i}&=\chi_{ij}\mathcal{E}_{j}+\mathcal{O}(\mathcal{E}^{2})\\\
\epsilon_{ij}&=1+4\pi\chi_{ij}\end{split}\;,$ (53)

The slope of $P(\mathcal{E})$ is reduced in the tilted case, because the tilts
compete with the polar mode. For the static calculations, the tilted and
untilted systems gave identical results. When the tilts are allowed to relax,
the relation between $P$ and $\mathcal{E}$ is

$P(\mathcal{E},\phi(\mathcal{E}))=\left(\frac{1}{\epsilon_{5}}+\frac{A}{P_{s}^{2}}\phi^{\prime}(\mathcal{E})^{2}\right)^{-1}\mathcal{E}\;,$
(54)

where $\epsilon_{5}$ is the dielectric constant of the untilted system. The
dielectric constant of the tilted system, $\epsilon_{20}$, is obtained in the
limit $\mathcal{E}\to 0$:

$\frac{1}{\epsilon_{20}}=\frac{1}{\epsilon_{5}}+\frac{A}{P_{s}^{2}}\;,$ (55)

where we used $\phi^{\prime}(\mathcal{E}\to 0)=1$. Rearranging Eq. (55) gives

$\frac{\epsilon_{20}^{-1}-\epsilon_{5}^{-1}}{\epsilon_{5}^{-1}}=\frac{A}{(P_{s}^{2}/\epsilon_{5})}\;,$
(56)

which is exactly the vertical axis of the phase transition diagram in Fig. 4.
Values of $\epsilon_{5}$, $\epsilon_{20}$, $A$ and
$\frac{P_{s}^{2}}{\epsilon_{5}}$ are given in Table 1 for all of the
functionals used. For the LDA (PAW) calculations, the electric field
calculations for the tilted system failed to converge, so we used the value of
$A$ from Ref. Gazquez  _et al._ , 2017 in Fig. 4.

The results in Table 1 were used to place the LAO/STO interface on the phase
transition diagram in Fig. 4. Although there is a relatively large distance
between the three points, each exchange-correlation functional predicts that
the LAO/STO interface is in region II, i.e. a single second order transition
at a reduced critical thickness $d_{\phi}$. The difference in LDA can be
attributed to the shallower tilt double well. The difference between PBE and
PBEsol appears to arise from the difference in the untilted dielectric
constants.

We also performed calculations to investigate the effect of strain on the
carrier transition, using the PBEsol exchange-correlation functional.
Calculations were repeated but with small amounts of compressive and tensile
strain applied to $a_{\text{STO}}$, up to $\pm$1%. The results are included in
Fig. 4, indicated by the smaller dots. We found that a compressive strain
moved the point towards the origin and a tensile strain moved the point up and
to the right, although it appears that a significantly larger amount of strain
than the ones investigated here would be required to change the order of the
carrier transition. It may be possible to change the character of the
transitions using biaxial strain for a system with a point which lies closer
to a boundary between regions, however.

## References

  * Ohtomo  _et al._ (2002) A. Ohtomo, D. Muller, J. Grazul, and H. Y. Hwang, Nature 419, 378 (2002).
  * Ohtomo and Hwang (2004) A. Ohtomo and H. Hwang, Nature 427, 423 (2004).
  * Li  _et al._ (2011a) L. Li, C. Richter, S. Paetel, T. Kopp, J. Mannhart, and R. Ashoori, Science 332, 825 (2011a).
  * Reyren  _et al._ (2007) N. Reyren, S. Thiel, A. Caviglia, L. F. Kourkoutis, G. Hammerl, C. Richter, C. Schneider, T. Kopp, A.-S. Rüetschi, D. Jaccard, _et al._ , Science 317, 1196 (2007).
  * Brinkman  _et al._ (2007) A. Brinkman, M. Huijben, M. Van Zalk, J. Huijben, U. Zeitler, J. Maan, W. G. van der Wiel, G. Rijnders, D. H. Blank, and H. Hilgenkamp, Nature materials 6, 493 (2007).
  * Li  _et al._ (2011b) L. Li, C. Richter, J. Mannhart, and R. C. Ashoori, Nature Physics 7, 762 (2011b).
  * Bert  _et al._ (2011) J. A. Bert, B. Kalisky, C. Bell, M. Kim, Y. Hikita, H. Y. Hwang, and K. A. Moler, Nature Physics 7, 767 (2011).
  * Thiel  _et al._ (2006) S. Thiel, G. Hammerl, A. Schmehl, C. Schneider, and J. Mannhart, Science 313, 1942 (2006).
  * Cen  _et al._ (2008) C. Cen, S. Thiel, G. Hammerl, C. W. Schneider, K. Andersen, C. Hellberg, J. Mannhart, and J. Levy, Nature materials 7, 298 (2008).
  * Förg  _et al._ (2012) B. Förg, C. Richter, and J. Mannhart, Applied Physics Letters 100, 053506 (2012).
  * Xie  _et al._ (2011) Y. Xie, Y. Hikita, C. Bell, and H. Y. Hwang, Nature communications 2, 1 (2011).
  * Irvin  _et al._ (2010) P. Irvin, Y. Ma, D. F. Bogorin, C. Cen, C. W. Bark, C. M. Folkman, C.-B. Eom, and J. Levy, Nature Photonics 4, 849 (2010).
  * Pallecchi  _et al._ (2010) I. Pallecchi, M. Codda, E. G. d’Agliano, D. Marré, A. Caviglia, N. Reyren, S. Gariglio, and J.-M. Triscone, Physical Review B 81, 085414 (2010).
  * Filippetti  _et al._ (2012) A. Filippetti, P. Delugas, M. Verstraete, I. Pallecchi, A. Gadaleta, D. Marré, D. Li, S. Gariglio, and V. Fiorentini, Physical Review B 86, 195301 (2012).
  * Assmann  _et al._ (2013) E. Assmann, P. Blaha, R. Laskowski, K. Held, S. Okamoto, and G. Sangiovanni, Physical review letters 110, 078701 (2013).
  * Liang  _et al._ (2013) H. Liang, L. Cheng, X. Zhai, N. Pan, H. Guo, J. Zhao, H. Zhang, L. Li, X. Zhang, X. Wang, _et al._ , Scientific reports 3, 1 (2013).
  * Hwang (2006) H. Y. Hwang, Science 313, 1895 (2006).
  * Pauli and Willmott (2008) S. Pauli and P. Willmott, Journal of Physics: Condensed Matter 20, 264012 (2008).
  * Huijben  _et al._ (2009) M. Huijben, A. Brinkman, G. Koster, G. Rijnders, H. Hilgenkamp, and D. H. Blank, Advanced Materials 21, 1665 (2009).
  * Mannhart and Schlom (2010) J. Mannhart and D. Schlom, Science 327, 1607 (2010).
  * Pentcheva and Pickett (2010) R. Pentcheva and W. E. Pickett, Journal of Physics: Condensed Matter 22, 043001 (2010).
  * Chen  _et al._ (2010) H. Chen, A. M. Kolpak, and S. Ismail-Beigi, Advanced Materials 22, 2881 (2010).
  * Warusawithana  _et al._ (2013) M. Warusawithana, C. Richter, J. A. Mundy, P. Roy, J. Ludwig, S. Paetel, T. Heeg, A. Pawlicki, L. F. Kourkoutis, M. Zheng, _et al._ , Nature communications 4, 1 (2013).
  * Chambers (2011) S. A. Chambers, Surface Science 605, 1133 (2011).
  * Pentcheva  _et al._ (2012) R. Pentcheva, R. Arras, K. Otte, V. G. Ruiz, and W. E. Pickett, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 4904 (2012).
  * Gabay  _et al._ (2013) M. Gabay, S. Gariglio, J.-M. Triscone, and A. Santander-Syro, The European Physical Journal Special Topics 222, 1177 (2013).
  * Stemmer and James Allen (2014) S. Stemmer and S. James Allen, Annual Review of Materials Research 44, 151 (2014).
  * Bristowe  _et al._ (2014) N. Bristowe, P. Ghosez, P. B. Littlewood, and E. Artacho, Journal of Physics: Condensed Matter 26, 143201 (2014).
  * Nakagawa  _et al._ (2006) N. Nakagawa, H. Y. Hwang, and D. A. Muller, Nature materials 5, 204 (2006).
  * Popović  _et al._ (2008) Z. S. Popović, S. Satpathy, and R. M. Martin, Physical review letters 101, 256801 (2008).
  * Pentcheva and Pickett (2009) R. Pentcheva and W. E. Pickett, Physical review letters 102, 107602 (2009).
  * Kalabukhov  _et al._ (2007) A. Kalabukhov, R. Gunnarsson, J. Börjesson, E. Olsson, T. Claeson, and D. Winkler, Physical Review B 75, 121404 (2007).
  * Siemons  _et al._ (2007) W. Siemons, G. Koster, H. Yamamoto, W. A. Harrison, G. Lucovsky, T. H. Geballe, D. H. Blank, and M. R. Beasley, Physical review letters 98, 196802 (2007).
  * Herranz  _et al._ (2007) G. Herranz, M. Basletić, M. Bibes, C. Carrétéro, E. Tafra, E. Jacquet, K. Bouzehouane, C. Deranlot, A. Hamzić, J.-M. Broto, _et al._ , Physical review letters 98, 216803 (2007).
  * Eckstein (2007) J. N. Eckstein, Nature Materials 6, 473 (2007).
  * Shibuya  _et al._ (2007) K. Shibuya, T. Ohnishi, M. Lippmaa, and M. Oshima, Applied Physics Letters 91, 232106 (2007).
  * Chen  _et al._ (2011) Y. Chen, N. Pryds, J. E. Kleibeuker, G. Koster, J. Sun, E. Stamate, B. Shen, G. Rijnders, and S. Linderoth, Nano letters 11, 3774 (2011).
  * Bristowe  _et al._ (2011) N. Bristowe, P. Littlewood, and E. Artacho, Physical Review B 83, 205405 (2011).
  * Lee  _et al._ (2012) S. W. Lee, Y. Liu, J. Heo, and R. G. Gordon, Nano letters 12, 4775 (2012).
  * Liu  _et al._ (2013) Z. Liu, C. Li, W. Lü, X. Huang, Z. Huang, S. Zeng, X. Qiu, L. Huang, A. Annadi, J. Chen, _et al._ , Physical Review X 3, 021010 (2013).
  * Li  _et al._ (2015) C. Li, Z. Liu, W. Lü, X. R. Wang, A. Annadi, Z. Huang, S. Zeng, T. Venkatesan, _et al._ , Scientific reports 5, 1 (2015).
  * Yu and Zunger (2014) L. Yu and A. Zunger, Nature communications 5, 5118 (2014).
  * Chen  _et al._ (2016) Z. Chen, H. Yuan, Y. Xie, D. Lu, H. Inoue, Y. Hikita, C. Bell, and H. Y. Hwang, Nano letters 16, 6130 (2016).
  * Kim  _et al._ (2015) S. K. Kim, S.-I. Kim, H. Lim, D. S. Jeong, B. Kwon, S.-H. Baek, and J.-S. Kim, Scientific reports 5, 8023 (2015).
  * Hosoda  _et al._ (2013) M. Hosoda, Y. Hikita, H. Y. Hwang, and C. Bell, Applied Physics Letters 103, 103507 (2013).
  * Bert  _et al._ (2012) J. A. Bert, K. C. Nowack, B. Kalisky, H. Noad, J. R. Kirtley, C. Bell, H. K. Sato, M. Hosoda, Y. Hikita, H. Y. Hwang, _et al._ , Physical Review B 86, 060503 (2012).
  * Eerkes  _et al._ (2013) P. Eerkes, W. G. van der Wiel, and H. Hilgenkamp, Applied physics letters 103, 201603 (2013).
  * Davis  _et al._ (2018) S. Davis, Z. Huang, K. Han, T. Venkatesan, V. Chandrasekhar, _et al._ , Physical Review B 98, 024504 (2018).
  * Liu  _et al._ (2015) W. Liu, S. Gariglio, A. Fête, D. Li, M. Boselli, D. Stornaiuolo, and J.-M. Triscone, APL materials 3, 062805 (2015).
  * Lei  _et al._ (2014) Y. Lei, Y. Li, Y. Chen, Y. Xie, Y. Chen, S. Wang, J. Wang, B. Shen, N. Pryds, H. Hwang, _et al._ , Nature communications 5, 1 (2014).
  * Gazquez  _et al._ (2017) J. Gazquez, M. Stengel, R. Mishra, M. Scigaj, M. Varela, M. Roldan, J. Fontcuberta, F. Sánchez, and G. Herranz, Physical review letters 119, 106102 (2017).
  * Aguado-Puente  _et al._ (2015) P. Aguado-Puente, N. Bristowe, B. Yin, R. Shirasawa, P. Ghosez, P. B. Littlewood, and E. Artacho, Physical Review B 92, 035438 (2015).
  * Bennett  _et al._ (2020) D. Bennett, M. Muñoz Basagoiti, and E. Artacho, Royal Society open science 7, 201270 (2020).
  * Caviglia  _et al._ (2008) A. Caviglia, S. Gariglio, N. Reyren, D. Jaccard, T. Schneider, M. Gabay, S. Thiel, G. Hammerl, J. Mannhart, and J.-M. Triscone, Nature 456, 624 (2008).
  * Stengel  _et al._ (2009) M. Stengel, N. A. Spaldin, and D. Vanderbilt, Nature Physics 5, 304 (2009).
  * Paruch and Guyonnet (2013) P. Paruch and J. Guyonnet, Comptes Rendus Physique 14, 667 (2013).
  * Jackson (1999) J. D. Jackson, “Classical electrodynamics,” (1999).
  * Lehnert  _et al._ (2000) H. Lehnert, H. Boysen, J. Schneider, F. Frey, D. Hohlwein, P. Radaelli, and H. Ehrenberg, Zeitschrift für Kristallographie-Crystalline Materials 215, 536 (2000).
  * Glazer (1972) A. Glazer, Acta Crystallographica Section B: Structural Crystallography and Crystal Chemistry 28, 3384 (1972).
  * Hatt and Spaldin (2010) A. J. Hatt and N. A. Spaldin, Physical Review B - Condensed Matter and Materials Physics 82, 1 (2010), 0808.3792 .
  * Fister  _et al._ (2014) T. Fister, H. Zhou, Z. Luo, S. S. A. Seo, S. Hruszkewycz, D. Proffit, J. Eastman, P. Fuoss, P. Baldo, H. N. Lee, _et al._ , APL Materials 2, 021102 (2014).
  * Benedek and Fennie (2013) N. A. Benedek and C. J. Fennie, The Journal of Physical Chemistry C 117, 13339 (2013).
  * Salje (1992) E. K. Salje, Physics reports 215, 49 (1992).
  * Lubensky and Rubin (1975) T. Lubensky and M. H. Rubin, Physical Review B 12, 3885 (1975).
  * Kretschmer and Binder (1979) R. Kretschmer and K. Binder, Physical Review B 20, 1065 (1979).
  * De Gennes (2018) P.-G. De Gennes,  _Superconductivity of metals and alloys_ (CRC Press, 2018).
  * Chandra and Littlewood (2007) P. Chandra and P. B. Littlewood, in  _Physics of ferroelectrics_ (Springer, 2007) pp. 69–116.
  * Luk’yanchuk  _et al._ (2009) I. A. Luk’yanchuk, L. Lahoche, and A. Sené, Phys. Rev. Lett. 102, 147601 (2009).
  * Sing  _et al._ (2009) M. Sing, G. Berner, K. Goß, A. Müller, A. Ruff, A. Wetscherek, S. Thiel, J. Mannhart, S. Pauli, C. W. Schneider, _et al._ , Physical Review Letters 102, 176805 (2009).
  * Berner  _et al._ (2010) G. Berner, S. Glawion, J. Walde, F. Pfaff, H. Hollmark, L.-C. Duda, S. Paetel, C. Richter, J. Mannhart, M. Sing, _et al._ , Physical Review B 82, 241405 (2010).
  * Takizawa  _et al._ (2011) M. Takizawa, S. Tsuda, T. Susaki, H. Hwang, and A. Fujimori, Physical Review B 84, 245124 (2011).
  * Slooten  _et al._ (2013) E. Slooten, Z. Zhong, H. Molegraaf, P. Eerkes, S. De Jong, F. Massee, E. Van Heumen, M. Kruize, S. Wenderich, J. Kleibeuker, _et al._ , Physical Review B 87, 085128 (2013).
  * Kumar  _et al._ (2020) S. Kumar, J. Kaswan, B. Satpati, A. Shukla, B. Gahtori, J. Pulikkotil, and A. Dogra, Applied Physics Letters 116, 051603 (2020).
  * Wang  _et al._ (2016) Y. Wang, W. Tang, J. Cheng, M. Behtash, and K. Yang, ACS applied materials & interfaces 8, 13659 (2016).
  * Moetakef  _et al._ (2011) P. Moetakef, T. A. Cain, D. G. Ouellette, J. Y. Zhang, D. O. Klenov, A. Janotti, C. G. Van de Walle, S. Rajan, S. J. Allen, and S. Stemmer, Applied Physics Letters 99, 232116 (2011).
  * Fan  _et al._ (2014) X. Fan, W. Zheng, X. Chen, and D. J. Singh, PLOS ONE 9, 1 (2014).
  * Maznichenko  _et al._ (2019) I. V. Maznichenko, S. Ostanin, A. Ernst, J. Henk, and I. Mertig, physica status solidi (b) , 1900540 (2019).
  * Kim  _et al._ (2018) Y. Kim, Y. M. Kim, J. Shin, and K. Char, APL Materials 6, 096104 (2018).
  * Yang  _et al._ (2016) K. Yang, S. Nazir, M. Behtash, and J. Cheng, Scientific reports 6, 34667 (2016).
  * Chen  _et al._ (2013) Y. Chen, N. Bovet, F. Trier, D. Christensen, F. Qu, N. H. Andersen, T. Kasama, W. Zhang, R. Giraud, J. Dufouleur, _et al._ , Nature communications 4, 1 (2013).
  * Niu  _et al._ (2017) W. Niu, Y. Zhang, Y. Gan, D. V. Christensen, M. V. Soosten, E. J. Garcia-Suarez, A. Riisager, X. Wang, Y. Xu, R. Zhang, _et al._ , Nano letters 17, 6878 (2017).
  * Yang  _et al._ (2020) X. Yang, H. Su, and G. Wu, Physical Review Materials 4, 016001 (2020).
  * Gonze  _et al._ (2002) X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, _et al._ , Computational Materials Science 25, 478 (2002).
  * Gonze (2005) X. Gonze, Zeitschrift für Kristallographie-Crystalline Materials 220, 558 (2005).
  * Gonze  _et al._ (2009) X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken, F. Bottin, P. Boulanger, F. Bruneval, D. Caliste, R. Caracas, M. Côté, _et al._ , Computer Physics Communications 180, 2582 (2009).
  * Gonze  _et al._ (2016) X. Gonze, F. Jollet, F. A. Araujo, D. Adams, B. Amadon, T. Applencourt, C. Audouze, J.-M. Beuken, J. Bieder, A. Bokhanchuk, _et al._ , Computer Physics Communications 205, 106 (2016).
  * Gonze  _et al._ (2020) X. Gonze, B. Amadon, G. Antonius, F. Arnardi, L. Baguet, J.-M. Beuken, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, _et al._ , Computer Physics Communications 248, 107042 (2020).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Perdew  _et al._ (2008) J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Physical review letters 100, 136406 (2008).
  * Marques  _et al._ (2012) M. A. Marques, M. J. Oliveira, and T. Burnus, Computer physics communications 183, 2272 (2012).
  * García  _et al._ (2018) A. García, M. J. Verstraete, Y. Pouillon, and J. Junquera, Computer Physics Communications 227, 51 (2018).
  * Hamann (2013) D. Hamann, Physical Review B 88, 085117 (2013).
  * Van Setten  _et al._ (2018) M. Van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. Hamann, X. Gonze, and G.-M. Rignanese, Computer Physics Communications 226, 39 (2018).
  * Perdew and Wang (1992) J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
  * Blöchl (1994) P. E. Blöchl, Physical review B 50, 17953 (1994).
  * Torrent  _et al._ (2008) M. Torrent, F. Jollet, F. Bottin, G. Zérah, and X. Gonze, Computational Materials Science 42, 337 (2008).
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Physical Review B 13, 5188 (1976).
  * Nunes and Vanderbilt (1994) R. W. Nunes and D. Vanderbilt, Phys. Rev. Lett. 73, 712 (1994).
  * Nunes and Gonze (2001) R. W. Nunes and X. Gonze, Phys. Rev. B 63, 155107 (2001).
  * Souza  _et al._ (2002) I. Souza, J. Íñiguez, and D. Vanderbilt, Phys. Rev. Lett. 89, 117602 (2002).
  * Umari and Pasquarello (2002) P. Umari and A. Pasquarello, Phys. Rev. Lett. 89, 157602 (2002).
  * Diéguez and Vanderbilt (2006) O. Diéguez and D. Vanderbilt, Phys. Rev. Lett. 96, 056401 (2006).
  * Berry (1984) M. V. Berry, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 392, 45 (1984).
  * King-Smith and Vanderbilt (1993) R. King-Smith and D. Vanderbilt, Physical Review B 47, 1651 (1993).