Polar morphologies from first principles: PbTiO$_3$ films on SrTiO$_3$ substrates and the $p(2 \times \Lambda)$ surface reconstruction

Low dimensional structures comprised of ferroelectric (FE) PbTiO$_3$ (PTO) and quantum paraelectric SrTiO$_3$ (STO) are hosts to complex polarization textures such as polar waves, flux-closure domains and polar skyrmion phases. Density functional theory (DFT) simulations can provide insight into this order, but, are limited by the computational effort needed to simulate the thousands of required atoms. To relieve this issue, we use the novel multi-site support function (MSSF) method within DFT to reduce the solution time for the electronic groundstate whilst preserving high accuracy. Using MSSFs, we simulate thin PTO films on STO substrates with system sizes $>2000$ atoms. In the ultrathin limit, the polar wave texture with cylindrical chiral bubbles emerges as an intermediate phase between full flux closure domains and in-plane polarization. This is driven by an internal bias field born of the compositionally broken inversion symmetry in the [001] direction. Since the exact nature of this bias field depends sensitively on the film boundary conditions, this informs a new principle of design for manipulating chiral order on the nanoscale through the careful choice of substrate, surface termination or use of overlayers. Antiferrodistortive (AFD) order locally interacts with these polar textures giving rise to strong FE/AFD coupling at the PbO terminated surface driving a $p(2 \times \Lambda)$ surface reconstruction. This offers another pathway for the local control of ferroelectricity.


INTRODUCTION
With the advent of advanced deposition techniques [1,2] has come a revolution in the engineering of thin film perovskite oxides and layered heterostructures for a variety of applications in nanoelectronics. These advancements have propelled research into the great variety of physical phenomena such systems can exhibit. These include enhanced colossal magnetoresistance [3], high-temperature superconductivity [4], the formation of interfacial two-dimensional electron and hole gases (2DEG/HG) [5,6] and the emergence of negative capacitance [7]. While these are becoming well-documented [8][9][10], new emergent phenomena resulting from exotic electrical polarization textures, including polar waves, vortices and polar skyrmion phases [11][12][13][14][15] are less well understood. The toroidal moment born from the chirality of these polar morphologies can give rise to strong electrotoroidic, pyrotoroidic and piezotoroidic effects [13,16,17] all of which show promise to be exploited in new low dimensional functional devices. It is the purpose of this work to investigate these polar morphologies in low dimensional ferroelectric (FE) PbTiO 3 (PTO) on SrTiO 3 (STO) substrates using state-of-the-art methods based on density functional theory (DFT) designed for demanding calculations with large numbers of atoms [18,19].
Further complex phenomena can arise by the interacting order parameters of the system. Notably, many perovskites (and heterostructures) are susceptible to both the antiferrodistortive (AFD) and FE distortions. In the bulk, these two modes were thought to suppress one another, although, recent evidence suggests that a cooperative regime may also exist [20,21]. At surfaces and interfaces we see phase coexis-tence. For example, at the PbO terminated [001] surface of PTO, we observe the AFD c(2 × 2) surface reconstruction. This is characterised by strong antiphase rotations of the TiO 6 octahedra about the [001] axis (or a 0 a 0 c − in Glazer's notation [22,23]) and is known to coexist with and mutually enhance in-plane ferroelectricity [24][25][26]. It is now popular to interface PTO with STO in the repeating (STO n /PTO n ) N superlattice for n alternating perovskite unit cells repeated N times in a layered heterostructure. In the case of an ultrashort period (n = 1), hybrid-improper ferroelectricity can arise from the coupling of AFD and FE modes. [27].
The (STO n /PTO n ) N heterostructure continues to be studied from a theoretical perspective. Most frequently considered is full periodic boundary conditions with infinitely repeated layers. This has been studied with first principles DFT [28] where the focus has been on the interactions of FE and AFD instabilities in the monodomain configurations. More recently, however, polydomain configurations have been studied and ferroelectric flux-closure domains have been shown to be stable [29]. Since these simulations are of the bulk superlattice, there can be no account for the c(2 × 2) surface reconstruction. Polydomain simulations of pure TiO 2 terminated PTO films (which do not give rise to enhanced surface AFD modes [24,25]) have been performed but neglect the STO substrate by instead choosing to adopt the fictitious free-standing film geometry [30]. This work then cannot make any account for substrate/interface effects or the intrinsically broken inversion symmetry for epitaxially deposited films on substrates. To the author's knowledge, there has been only a single DFT-based study treating both the surface and the STO substrate. This study probed the nature of the theoretically proposed 2DEG/HG pair [6] where thicker films (≥ 14 unit cells) of the FE monodomain out-of-plane configuration were considered as large depolarising fields are known to suppress this configuration in thinner films [31]. Multiple FE domains are considered a competing phase but no 2DEG/HG pair emerges as the resulting flux closure domains are an alternate mechanism for screening the depolarising field [32].
The polydomain configuration is challenging to simulate from the perspective of DFT. The difficulty arises from the increased number of atoms N in the supercell as a result of treating a fixed domain period Λ. If we are also to treat a finite thickness of film t, the Kittel scaling law requires that the equilibrium Λ too increases since Λ ∝ √ t [33,34]. We must also account for a significant amount of substrate as the large ferroelectric polarization of PTO may be able polarise a few layers of STO. Finally, if we wish to include AFD modes in our simulations (like the c(2 × 2) surface reconstruction), the periodicity in the [010] direction must be doubled as to not frustrate the octahedral rotations. To include all of all of these effects, calculations with a supercell containing a few thousand atoms are needed rather than the few hundred that conventional planewave DFT is able to handle [35,36]. To overcome these limitations, some resort to alternative methods including phase fields, shell models, Monte Carlo and second principles [12,13,[37][38][39][40][41][42]. Whilst these approaches each have their merits, they cannot universally serve as a replacement for full DFT. Most of these methods implicitly accept that full DFT offers superior accuracy since they either are or can be parameterised by the theory [38,41,42]. Further, all of these methods make no account (or little account [38]) for the electronic structure, limiting the ability of simulations to adapt to new chemical environments like at surfaces and interfaces with a substrate. In this work, we utilise a novel variation on full local orbital DFT which allows us to consider systems of a few thousand atoms whilst preserving high chemical accuracy [19].
We provide a full first principles study on ultrathin films of PTO down to a single PbO monolayer (N z = 0) up to 9 unit cells (N z = 9) where N z is the number of perovskite unit cells stacked in the [001] direction. We treat explicitly the STO substrate and the PbO termination invoking the c(2 × 2) surface reconstruction. As well as treating paralelectric films, we treat two monodomain configurations of the polarization (P || [100] and P || [110]) and stripe domain patterns comprised of alternating P || [001] and P || [001] domains. We do not consider monodomain polarization oriented in the out-of-plane directions (P || [001] or P || [001]) since a previous work has shown this is suppressed by the depolarising field within our range of thicknesses (also verified to be true within our simulations). We characterise the competing energetics of the different film geometries as well detailing the resulting polar morphologies, interactions with the surface reconstruction and other structural features. This work is then organised as follows: section II A constitutes an overview of the multi-site support function method as applied in this study whilst Section II B presents finer details associated with the simulation method. Section II C details the treated film geometries for each thickness of film. In section III A we compare the competing energetics between each film geometry as a function of thickness. In section III B, the local polarization fields of the various film geometries are discussed and in Section III C, we detail the amplitudes of local AFD modes, including the c(2 × 2) surface reconstruction, its coupling with surface polarization and the asymmetrical rumpling of surface/interfacial Pb cations. We conclude this work in Section IV providing a summary of the findings of this investigation, commenting on the impact this work has on the field and suggesting new areas for which our novel method is applicable.

II. THEORETICAL METHOD
Simulations based on DFT have provided accurate firstprinciples descriptions of countless condensed matter and molecular systems since the 1970s. The method is centered upon the self-consistent calculation of the charge density n(r) for a system of interacting ions and electrons occupying independent Kohn-Sham [43] orbitals. The ground-state total energy of the system is evaluated as a uniquely defined functional of this n(r) as detailed in the Hohenberg-Kohn theorem [44]. This method would allow for exact calculations if it weren't for the small electronic exchange and correlation terms for which the exact functional is unknown, thus, (increasingly accurate) approximations must be made. In this work, calculations were performed using the implementation of DFT present in the CONQUEST code [18,45]. This code is designed for large-scale, massively parallel simulations on thousands (even millions [46]) of atoms, utilizing the latest in high performance computing (HPC). While we now focus on the multi-site support functions used in this work (Section II A), we refer the reader to [18,35,36] for further details on the CONQUEST methodology and to [45] for a recent (freely available) public release of the code.

A. Multi-site support functions
This section seeks to provide a convenient overview of the multi-site support function method (MSSF) method as implemented in CONQUEST [18,45]. For a more detailed discussion, the reader is referred to [19,47] and references therein.
We begin by considering the Kohn-Sham equations in a non-orthogonal basis in terms of the Hamiltonian and overlap matrix elements for ionic sites i, j orbital indices α, β, eigenvalues nk and wavevector k. Here, the Kohn-Sham eigenstates are written as a linear combination of support functions φ iα (r) for support function coefficients c nk iα . A support function is a local orbital which we define to become zero at some cutoff radius r c (φ iα (r) = 0, |r − R i | ≥ r c ). It is then possible to represent these support functions in a compact, efficient basis of atomic-centered functions whose spatial cutoff leads to the onset of matrix sparsity (important to the O(N ) mode of operation). How we choose to do so is paramount as this determines the accuracy and computational effort required for the calculation of the electronic groundstate. For accurate calculations we generally employ multiple-ζ basis-sets [48][49][50]. The most accurate evaluation of the support function coefficients here comes from the mapping of one ζ to one support function and the subsequent diagonalisation of a large H k . This, however, comes with the burden of the greatest computational cost since this operation scales cubically with the total number of basis functions. Another choice is to make a basis set expansion of φ iα (r) where index χ iζ is ζth atomic basis function on atomic site i. The coefficients d can be found through an optimisation method such as the conjugate gradients scheme. Another choice for equation 2 is to build a support function from all of the basis functions on atomic site i. This is a simple extension of eq 3 where instead of having ζ ∈ α, α spans all ζ on site i. While performing any of the previous two contractions does indeed reduce the number of required support functions, we are still unable to achieve a minimal representation due to atomic symmetry constraints [51].
To achieve this minimal representation, we use the multisite contraction where k is an atomic site enclosed within a sphere of radius r MS about target atom i, inclusive of i. As the coefficients C now span a subspace of a local molecular orbital (MO), we are now relinquished from the atomic symmetry constraints. To find the optimal balance between efficiency and accuracy, we use the local filter diagonalisation method [47,52] to evaluate the coefficients. This requires the solution of a local eigenproblem in a subset s enclosed by a sphere of radius r LD (≥ r MS ) for each target atom.
For subspace eigenvalues s , MO coefficients C s and where H s and S s are the local subsets of the system-wide Hamiltonian and overlap matrices respectively. By projecting C s onto a set of trial vectors t, we obtain the contraction coefficients C Lastly, C is mapped onto the corresponding position in the contraction coefficient vectors, whose elements are extracted for use in eq 4.
The MSSF method therefore introduces two new adjustable parameters to the DFT calculation; r LD and r MS . The former expands the space used for the local filter diagonalisation and the latter modifies the representation of φ iα (r). The accuracy of the contraction can be improved by increasing both of these subject to r MS ≤ r LD . In practice, convergence can be achieved for relatively modest values of both (1-2 lattice constants of the bulk). Using this method, high accuracy calculations with a few thousand atoms can feasibly be performed on most standard HPC systems. Structural relaxations of ≈ 1000 atoms (to a stringent force tolerance, as seen in this work) can be expected to take 1-2 weeks of wall-time on ≈ 200 physical cores. Simulations of ≈ 2000 atoms may need ≈ 500 physical cores to meet the memory requirements, dependant on the memory-per-node on the HPC system.

B. Simulation details
Within the MSSF mode of operation described in section II A, we find that r MS = r LD = 6.35Å sufficiently converges relevant quantities (discussed in the supporting information). All simulations use the local density approximation (LDA) as parameterised by Ceperley, Alder, Perdew and Zunger [53,54] to describe exchange and correlation effects. This functional has recently been shown to perform well against higher rung functionals for describing the properties of the perovskite oxides [55]. Although both hybrid and meta-GGA functionals better predict the experimental structural properties, they still overestimate the bulk polarization of PTO. Since this is the primary order parameter of our system, LDA is a good choice as the magnitude of the overestimate is much less.
Optimised Vanderbilt norm-conserving pseudopotentials are used to replace core electrons. [56,57]. We use the scalar-relativistic variety available in the PseudoDojo library (v0.4) [58]. These are used as an input for the ONCVPSP code (v3.3.1) [56]. We then generate a double-ζ plus double-polarization (DZDP) basis set of pseudo-atomic orbitals (PAOs) describing the Pb 5d 6s 6p 6d p , Sr 4s 4p 5s 4d p , Ti 3s 3p 4s 3d 4p p and the O 2s 2p 3d p orbitals respectively. Those orbitals superscripted with p are polarization orbitals aimed at increasing the angular flexibility of the basis. Ti 3s, Ti 3p and Pb 5d orbitals are treated as semi-core states described only with a single ζ. These PAOs are used to represent the support functions in section II A. The PAO basis constructed using this set of pseudopotentials has been shown to provide high accuracy lattice constants, bulk moduli and electronic properties [48,49].
Momentum space integrals are performed on a non-Γcentered 6/N x × 6/N y × 1 uniform mesh as described by Monkhorst & Pack [59] where N x,y are the number of perovskite unit cells included in the [100] and [010] directions of the supercell respectively. While many real-space inte-grals are performed using intuitive analytic operations [51], some are performed on a fine, regular integration grid with a plane-wave equivalent cutoff of 300 Ha. This yields the cubic P m3m lattice constants of PTO and STO (a PTO and a STO ) as 3.904Å and 3.874Å respectively. Further, the tetragonality (c/a) of FE P 4mm PTO is obtained as 1.04 producing a spontaneous polarization of 79.02 µC/cm 2 as calculated with Resta's method [60,61] within the modern theory of polarization. This method is equivalent to the Berry phase [62] formula of King-Smith and Vanderbilt [63] in the limit of large cells; we use 10 PTO unit cells in the direction of the FE distortion to achieve convergence. Each of the structures in section II C are fully relaxed (subject to constraints) with quenched molecular dynamics until the maximum absolute value of the force on every atom falls below 0.01 eV/Å. Since our slabs feature asymmetric surface terminations, a small dipole moment can develop in the out-of-plane direction. This propagates to a spurious electric field across the supercell which we correct using the scheme of Bengtsson [64].

C. Supercell configurations
The film geometries considered in this work are displayed in Figure 1. Universal to all structures is a fixed amount of STO substrate. We find that 8/7 alternating SrO/TiO 2 monolayers is sufficient to converge the relaxed ionic positions of the PTO films (demonstrated in the supporting information, in agreement with a similar classical study [42]). An interfacial region a I = 1/2(a STO + c PTO ) exists between the first PbO layer and last SrO layer which we assign to belong to neither the film or substrate. The in-plane components of the supercells are held to integer multiples of a STO . The two bottom-most monolayers of the substrate have their atomic positions fixed in structural relaxations. These measures ensure we are applying a realistic mechanical strain to the PTO film as well as simulating the effects of a semiinfinite substrate. To limit unfavorable interactions with images of the film in the [001] direction, we introduce a total of 20Å of vacuum. All supercells have the generalised di- We simulate films of thickness N z = 0, 1, 2, 3, 5, 7 and 9 formula units of PTO where a thickness of 0 corresponds to a single PbO monolayer. Films of N z ≥ 3 unit cells were increased in steps of two unit cells such that the equilibrium domain period predicted by Kittel scaling could increase by an even integer number of unit cells (a requirement for domains to have an equal number of unit cells). This range of film thicknesss could encompass different energetically stable geometries including a transition between FE monodomain, polydomain and possible intermediate phases. It also spans low dimensional films with strong interface coupling with such effects decreasing with increased thickness.
We choose to treat the following supercell configurations with and without the influence of AFD modes. Note also that structures treated without these modes do not show the c(2×2) surface reconstruction. To do so, we must set N y = 2 FIG. 1. The initial supercell configurations for the Nz = 3 films before structural relaxation. Shown here are the supercells not including AFD modes. Each configuration is however also treated with AFD modes following the explanation in section II C. a) The paraelectric supercell constrained such that spontaneous polarization cannot emerge. b) The monodomain in-plane ferroelectric case (P || [100] is shown here, but we also treat P || [110]) constrained such that spontaneous polarization cannot develop in the out-of-plane direction. c) The polydomain ferroelectric case with equally sized up and down domains for the ferroelectric polarization. Shown here is the Λ = 6 case. and have N x set to multiples of 2. The important AFD distortions for both STO and PTO are the zone-boundary R + 4 and M + 3 modes. The former is equivalent to a a 0 a 0 c − Glazer tilt system [22,23] while the latter is a 0 a 0 c + . The former has been found previously to be the most energetically favorable so we do not treat the M + 3 mode [65]. The R + 4 mode can be defined with a single rotation angle θ z as shown in the inset of Figure 2. We find that the angles θ z = 6.21 • and θ z = 4.52 • (Figure 2) define minima in the energy of the cubic supercells of STO and PTO respectively. Since simulations are performed with the in-plane lattice constants of STO, we initialise supercells with the optimal strained PTO rotation angle θ z,σ = 5.20 • as found in the tetragonal supercell. These results indicate that strains of < 1% in PTO are able to both increase the optimal θ z and increase the depth of the energy minima. These angles overestimate what is seen in experiment. For I4/mcm STO, θ z ≈ 2.1 • [66] while AFD modes are not observed in the PTO bulk. Calculations using hybrid functionals have been able improve this angle for STO (θ z = 1.9 • [67]) but are not used in this study in part because of the computational cost but also due to the aforementioned overestimate in the bulk polarization of PTO.

Paraelectric
First we consider films of the high symmetry, paraelectric, non-polar variety ( Figure 1a). These supercells are initialised by PTO formula units with the high symmetry cubic fractional atomic positions and optimal cubic out-of-plane lattice constant c PTO = 3.904Å. Although spatial inversion symmetry is intrinsically broken by the composition of our supercells, we constrain atomic relaxations of these films to the lowest symmetry space group (P 4mm) carefully preventing any cation-anion counter motion. This allows for the most degrees of freedom in relaxations without incurring any in-trinsic ferroelectric polarization. It is the purpose of these films to provide a baseline for the relative stability of polar phases.

Monodomain in-plane ferroelectric
In this supercell, we allow monodomain ferroelectric polarization to develop in the [100] and [110] directions (Figure 1b). PTO formula units retain the same c PTO as the paraelectric case but atomic positions now correspond to the Γ − 4 mode of the primitive cubic PTO unit cell orientated in the [100] direction. This induces a ferroelectric polarization of 56.9 µC/cm 2 in the infinite bulk as calculated with the method described in Section II B. To create [110] polarization, we simply displace the metal cations along the [110] direction and displace Oxygen anions along the [110] direction proportional to the magnitude of the bulk Γ − 4 mode. During structural relaxation, we apply symmetry constraints to prevent the development of out-of-plane polarization [25] although the non-trivial depolarising field that would result naturally suppresses it.

Polydomain ferroelectric
Here we consider films intialised with a striped domain structure consisting of alternating regions of PTO polarised in the [001] and [001] directions which we will from this point onwards refer to as up and down domains respectively ( Figure 1c). Up and down domains are equal in size (N up = N down ) and together form a full domain period Λ. Equivalently, N up + N down = N x ≡ Λ. Domain walls are chosen to be centered on the PbO plane. This choice, however, is arbitrary as the energy difference between this and the TiO 2 centering is found to be ∼ 1 meV per unit cell [30]. This is (slightly) beyond the resolution of our calculations, a fact noted in a comparable study. [29].
It is found that there is a minimum thickness of film for this type of ferroelectricity to occur. Theoretical results of the free standing PTO film have shown that the polydomain ferroelectric film is lower in energy than the paraelectric configuration [30] but make no account of possible monodomain in-plane orientations for the polarization. In experiments conducted with an STO substrate [34] it was found that that the polydomain configuration was only observed above 3 unit cells in thickness (in the temperature range of 311-644K). We also confirm that no out-of-plane component of P remains after structural relaxation of polydomain N z = 1 and 2 films (Λ = 4). We then only consider this configuration for those film thicknesses N z ≥ 3. PTO unit cells are initialised with the strained FE P 4mm unit cell. When the in-plane constants are constrained to a STO , this results in c = 4.049Å and a slightly enhanced polarization of 80.07 µC/cm 2 .
An accurate account for the polydomain structure requires us to work at the equilibrium domain period. To avoid the need to manually find this domain period for each thickness (which would require us to simulate several domain periods for each film thickness), we use knowledge from previous experimental and theoretical data. In particular, we find that both theory and experiment agree upon Λ = 6 when N z = 3 [30,34]. We use X-ray diffraction data [33] for the remaining film thicknesses. We choose the nearest even number of unit cells (to preserve the periodicity of the AFD rotation pattern) corresponding to the experimental domain periods. We then have for the N z = 3, 5, 7, 9 films, domain periods of Λ = 6, 8, 10, 12 unit cells respectively. Films of this configuration are free from constraints during structural relaxation.
In a study of monodomain out-of-plane polarization, it was found that P [001] (towards the STO substrate) becomes more polar than P [001] (towards the vacuum). For polydomain films, this could mean that up domains are less polar than down domains. This would result in a small net dipole moment in the [001] direction (and a spurious electric field) developing during structural relaxation which is corrected for using the scheme described in Section II B.

A. Competing phases
The film geometries considered have energetics which evolve as a function of film thickness. Figure 3 displays this behaviour indicating the favorability of different phases. We measure this favorability as the energy difference ∆E between the geometry of film in question and the paraelectric film. We choose to measure ∆E in meV/atom since the (more common) definition of meV/formula unit would vary with film thickness. This is because as the film thickness increases, the Pb/Sr fraction (upper axis of Figure 3) increases as an artefact of a fixed amount of STO substrate. As a result, we expect to see a component energy lowering from ferroelectric phases as we increase N z since the energy of PTO is lowered with the onset of ferroelectricity (by 9.58 meV/atom in the bulk). We then expect to observe a rise in energy for purely AFD phases since the the fraction of STO, favoring AFD modes, has decreased. The energetics of film thickness N z = 0 are not present on Figure 3 since no ferroelectric phase was stable. Adding AFD rotations does however lower the energy compared to the paraelectric film by an amount similar to N z = 1.
Considering first the monodomain in-plane ferroelectric films, Figure 3 shows that the favored axis for the polarization is always [110] as was shown in a DFT study of the free standing film under compressive strain [68]. This is true with and without the influence of AFD modes. This favorable direction seems to diminish with film thickness becoming almost degenerate with [100] polarization at N z = 7 for the films not influenced by AFD modes. When AFD modes are taken into account, the degree of favorability for [ Figure  1 versus the film thickness in PTO unit cells, Nz. Since we use a fixed amount of STO substrate, the formula unit of the film alters with Nz. This is accounted for by the upper x-axis indicating the Pb/Sr fraction. The area in grey indicates the domain of Nz for which the polydomain configuration is not considered.
supercell (of length √ 2a STO ) when compared with a distortion parallel to one of the pseudocubic axes (of length a STO ) relieves this stunting. We can also deduce whether coupling between AFD and FE modes is cooperative or competitive. The sum of ∆E for the FE [100] curve and the AFD curve is always lower than the combined FE [100] + AFD curve. This indicates that the the coupling is competitive with AFD modes suppressing FE ones and vice versa as is usually true for bulk modes. When making the same comparison for [110] polarization (which is not observed in the bulk), however, it very closely mirrors the combined FE [110] + AFD curve. This suggests that FE and AFD modes are at worst independent of one another, but, for N z = 2 or 3 are mildly cooperative, with the FE [110] + AFD curve being lower in energy than the sum by ≈ 0.2 meV/atom (close to the resolution of the simulation).
Remarkably, the polydomain configuration is not univer-sally the ground state. Monodomain [110] polarization is the lowest in energy until a film thickness of 4-5 unit cells. This is close to the experimental observation of a polydomain structure at a thickness of 3 unit cells [34] with the difference perhaps being an artefact of finite temperature in experiment. We note also that these results agree with the theoretical findings of Shimada [30] in that at a thickness of 3 unit cells, the polydomain configuration is lower in energy the the paraelectric film. In this work, however, while the energy is lowered by this geometry, monodomain in-plane [110] polarization is favored at this depth. It is important to note that this work treats the PbO termination whilst the work of Shimada treats the TiO 2 termination. As we discuss later (in Section III B), the N z = 3 polydomain film does not condense the flux-closure domain morphology but instead shows a polar wave (Figure 5b). This clearly will have an impact on the energetics. Comparing the energy of the combined polydomain FE + AFD curve with the sum of the polydomain FE and AFD curves, we find that the latter is always lower in energy (a gap which widens with increasing film thickness). This suggests that that polydomain FE competes with AFD modes. This effect is not as drastic as the competition between FE modes and monodomain [100] polarization however.

B. Polarization morphologies
In this section we analyse and compare the polar morphologies of the different films using a metric known as the local polarization, P (i) . We define P (i) per 5-atom perovskite unit cell using a linear relationship between the local mode and the Born effective charges Z * in the manner first suggested by Resta [61] (now used in similar works [29,69]). We discuss this method in more detail in the supporting information. A vector field of this quantity has been calculated for all polar structures presented in Figure 4 and 5.
In Figure 4a we show the Ti-centered local polarization along the up and down domain centroids. The domain centroid here is a string of Ti centered unit cells in the vertical direction located at the centre of a domain. It is at the down domain centroid that the maximal local polarization can be found, buried in the upper third of the PTO film. Indeed, there is a discrepancy in polarization between the up and down domains throughout the entire film, leading to a small net dipole moment in the [001] direction. This effect can be explained by the compositionally broken inversion symmetry present in even the highest symmetry films (the relaxed geometry described in Section II C 1). The result is that P [001] is favored by a built-in bias field directed towards the substrate. While we find that this field is local to the first few PTO surface monolayers, we suggest that the enhanced local P [001] modes at the surface spread to the rest of the domain due to a finite correlation length associated with the polar atomic displacements [70]. This leads to an indentation of the substrate (as seen in figure 7) creating extra volume for the down domains and enhanced local tetragonality which mutually couples with the polarisation. The opposite argument is also true for the local P [001] modes within the up domain whereby the internal bias field is now depolarising. We predict that such a disparity between the up and down domains will diminish with increasing film thickness, tending to zero as the film thickness becomes much larger than the correlation length of the enhanced local polar modes. The scenario explained above is analogous to the observation of built-in bias fields present three component or tricolor superlattices [71,72] where our third component is the vacuum region .
We find that the maximal polarization increases with film thickness for the polydomain films. Although we only display N z = 3 and N z = 9 in Figure 4a, the rise is gradual when considering the maximal polarization of the N z = 5 and N z = 7 films also. This reflects the results of synchrotron Xray diffraction [34] whereby satellite peak intensity (indicative of domain polarization) is an increasing function of film thickness. Such a phenomenon is likely a result of decreasing depolarising field strength with increasing film thickness as is observed for thicker P [001] films [73]. It is notable that by N z = 9, this maximal polarization (for the case without AFD modes) exceeds the strained PTO bulk figure by 26%. The same enhancement is not seen for the N z = 3 which suffers a 12% reduction. It can be seen for all film geometries that the influence of AFD modes tends to reduce the polarization in the film. Figure 4a shows that allowing for AFD modes produces a local polarization reduction of ≈ 15 µC/cm 2 for all films. This reduction, however, becomes more severe as we reach the surface where we suggest that the enhanced local rotational modes at the surface reconstruction compete strongly with the polarization.
The majority of the relaxed polydomain structures, as seen in similar works, form the flux closure morphology ( Figure  5a). That is, the local polarization gradually rotates through 180 • across the domain wall at the top of the film and rotates through -180 • at the bottom as a mechanism for screening the depolarising field. The result is counter-tilting vortexlike domain walls (Figure 5a, red areas), with a small vertical area for ≈ 0 polarization at the vortex center. This domain morphology forms fully for the N z = 5, 7, 9 films with and without the influence of AFD modes. We note that whilst these domain structures share similarities with other works, they do have some key differences. The vortex cores are not located centrally (in the z-direction) in the film. They are instead shifted towards the substrate. We also see that at the surface, flux is closed more abruptly than at the interface with the substrate. At the interface, the polarization penetrates (≈ 2 unit cells) into the substrate helping to close the flux. Together with the previously mentioned tendency for a stronger polarization in the [001] direction, this makes for a more asymmetrical morphology than those observed in the superlattice or free standing film arrangement [29,30]. In contrast to the work of Shimada [30], the N z = 3 film does not appear to form complete a flux closure domain structure. Examining Figure 5b we see that the flux does not fully close at the interface with STO. Instead, the polar morphology is S -like or wave-like, in this case orientated in the [100] direction giving the film a net in-plane macroscopic polarization. We deduce that at the surface, the flux does close by analysing the displacements of the terminating PbO layer (whose polarization vectors are not calculable with our method as explained in the supporting information). This gives rise to small cylindrical vortices near the surface which are sometimes known as cylindrical chiral bubbles (Figure 5b, red area) [11,12].
Such an instability has been predicted as an intermediate phase in phase field simulations under applied electric fields in the superlattice arrangement [13] as well as at zero field in high-angle annular dark field (HAADF) images [11]. This polar topology is also observed in thin PZT films as predicted with Monte-Carlo simulations of a first principlesbased Hamiltonian [12]. Although such a polar topology has not yet been observed in thin PTO films, we note that pre-vious XRD studies do not explicitly rule out the coupling to the characteristic in-plane wavevectors for the N z = 3 film [33,34,74]. In this work, the polar-wave morphology appears as an intermediate phase between full flux-closure domains and in-plane polarization like in reference [13] but in the absence of an applied field. This is replaced by the built-in bias field emergent from the compositionally broken inversion symmetry. Since this bias field is directed to suppress up domains but enhance down domains, up domains now have an increased critical thicknesses for the absolute suppression of out-of-plane polarization compared to down domains. We suggest that at N z = 3, we are below this critical thickness for up domains (where in-plane FE modes are now favored) but still above it in down domains. The resulting polar wave texture is then emergent from the closing of flux between P [001] and P [100] (although P a) b) [001] [010] [100] [100] is equally favored) domains as shown in figure 5b. This finding suggests that control over polar morphologies can be achieved in ultrathin films by careful engineering of the boundary conditions. Specifically, we suggest that the builtin bias field can be manipulated through the choice of substrate, film surface termination or use of overlayers. This principle of design offers a promising avenue for the manipulation of chiral order in low-dimensional devices operating through the control of toroidal moments. Figure 4b and 4c show the polarization profiles for the films initialised with uniform polarization oriented in the [110] and [100] directions respectively. For the former, we find that for most cases, P x ≈ P y (hence P remains aligned along [110]) apart from the thinnest N z = 1 and 2 films. N z = 1 in particular, when coupled with AFD modes, (red line) polarization rotates strongly to be mostly polarised in [001]. For both [110] and [100] polarised films, we see that these films become polar as we cross the STO/PTO interface, quickly adopting a bulk like value. Then, for films without AFD modes, we see a polar enhancement near the surface. For those with AFD modes, we see a polar reduction followed by enhancement near the surface contrasting with the strong reduction for the polydomain films. Computing the norm P 2 x + P 2 y of the polarization for the [110] and [100] films in the bulklike region shows us that the [110] films have a polarization slightly enhanced (+2µC/cm 2 ) compared to just [100] polarization. What is notable for the monodomain films is that (past a thickness of two unit cells at least) there is no trend for the behaviour of the local polarization for increasing film thickness. Each film shows the same bulk-like polarization then the same surface reduction or enhancement.

C. The p(2 × Λ) surface reconstruction
An analysis of the competition between local FE and AFD modes can provide valuable insight into the design of new low dimensional devices whereby FE and AFD modes can be tuned to enhance their functional properties. In this section, we analyse the interaction of FE and AFD modes with a focus on strong coupling within the reconstructed surface layers of the PbO terminated films. Figure 6 shows the local evolution of the R + 4 octahedral rotation angle θ z along the [001] direction (Figure 6a) and across a domain period Λ in the [100] direction (Figure 6b) for the N z = 9 (Λ = 12) film. The rotational behaviour in Figure 6a (left) is similar to the behaviour reported by Bungaro [25] whereby the reconstructed area couples cooperatively with in-plane [100] polarization. We show that this also true for P || [110] (polarization is locally enhanced at the surface as shown by Figure 4b and 4c) with the rotation pattern θ z (z) being almost indistinguishable from the [100] polar film. We note that whilst this mutual AFD/FE enhancement is active at the surface, for the rest of the film, AFD/FE modes are mutually reduced when compared to the paraelectric structure and the bulk. The θ z (z) trend for the polydomain films (Figure 6a, right) is more complex. Rotations are almost completely quenched (reduced to ≈ 2 • ) at the up and down domain centroids compared to the domain wall (apart from the reconstructed region). We suggest that this occurs since the maximal polarisation at the centroids outcompetes AFD modes. We see similar behaviour for all N z for the effects listed above, even for the strength of the reconstruction angle. This is remarkably resistant to finite size effects, even at a single PbO monolayer (N z = 0), θ z persists at around 12 • despite the change in chemical environment for the surface TiO 6 octahedra. Figure 6b shows a coupling between local polarization at the surface (upper) and the surface reconstruction angle (lower). Across a domain period Λ, we see that θ z (x) modulates by 1.75 • . |θ z (x)| peaks close to P x = 0 which is remarkable since this is precisely the opposite behaviour to the monodomain in-plane films where FE and AFD mutually support one another. Since P x (x) and P z (x) are π/2 radians out of phase with each other in the surface layer, it can also be said that the peaks in |θ z (x)| coincide with extrema in |P z (x)|. We see also the the magnitude of P z (x) has an impact on the height of the peaks in |θ z (x)|. That is, since the down domains are more polar than the up domains, the larger down polarization in the surface layer reduces the |θ z (x)| peak and vice versa for the for the smaller polarization in the up domain. We note that such AFD/FE couplings are different to those found in a recent shell model study of the free standing film [41]. In this study, the symmetrical [001] boundary conditions for the film modulate θ z (x) over Λ/2 instead. We also note that whilst the strongest octahedral rotations θ z are about the [001] axis, we also see smaller rotations of ≈ 2-3 • about the [100] axis in the surface reconstruction. In addition to this, at the surface of the polar wave film for N z = 3, we observe rotations about all three pseudocubic axes in the surface reconstruction. The strongest rotations are still about the [001] axis however. The broken inversion symmetry leads to further structural effects for the polydomain films. We can measure the distortion of the films in the vertical direction by considering the Pb displacements in the terminating PbO layer with the vaccuum and at the STO interface. 7 shows an exemplar calculation for the N z = 7 (Λ = 10) film. The other film thicknesses all resemble the behaviour of this film. With the exception of a small dip in ∆z at the domain wall for the surface rumpling (a mechanism to reduce the domain wall formation energy), we see that interface and surface rumpling behave spatially like two sinusoids π/2 radians out-of-phase with one-another, reflected about the [100] direction. This is in contrast to the PTO/STO superlattice configuration [29] where the intact inversion symmetry preserves symmetrical rumpling (in-phase sinusoids). Due to the periodicity of octahedral rotations and Pb cation surface rumpling in the [001] direction, it is clear that for the polydomain films we must reconsider the labelling of the surface reconstruction. For monodomain and paraelectric films, the Wood's notation [75] of c(2 × 2) remains correct but an analysis of the symmetry for the the polydomain case reveals this must alter to become p(2 × Λ). We suggest that these fine FE/AFD interactions could now be observable in experiment thanks to recent advances in integrated differential phase contrast (iDPC) imaging [76]. In contrast to HAADF, iDPC images yield the positions of the metal cations and oxygen anions resolved at a subunit cell level [76]. This allows for the direct measurement of local FE and AFD modes offering an exciting new avenue for direct comparison with atomistic results.

IV. CONCLUSIONS
We have used large-scale DFT calculations to simulate the energetics, interaction with AFD modes and other structural features for various polar morphologies present in thin PTO films on STO substrates. Our method is successful in providing accurate first principles results for systems comprised of a few thousand atoms, well beyond what is feasible with traditional plane-wave based methods, venturing towards sys- tems sizes usually simulated with Monte Carlo or phase field techniques. This method has allowed for the explicit simulation of the STO substrate, multiple ferroelectric domains and doubled periodicity in the [010] direction used include AFD modes and the surface reconstruction. This has ensured that our simulations represent realistic experimental conditions when compared to previous works which neglect one or more of these features. Such simulations can be performed on standard HPC systems with high accuracy calculations of ≈ 1000 − 2000 atoms becoming feasible using ≈ 200 − 500 physical cores.
We have demonstrated the stability of the polydomain film geometry compared to monodomain phases. We have found that the polydomain case becomes more energetically favorable between 4-5 unit cells in thickness, close to the experimental observation at 3 unit cells. We find that the general effect of including AFD modes is to lower the energy significantly (in most cases more than the FE distortion) and to suppress the amplitudes of local polar modes (apart from at the surface of monodomain in-plane FE films, where they are mutually cooperative).
We show that polydomain films display the flux-closure domain morphology for N z = 5, 7 and 9 whilst the N z = 3 film shows the similar polar wave morphology with cylindrical chiral bubbles as an intermediate phase between full flux closure domains and in-plane ferroelectricity. We find that local polarization is enhanced at the domain centroids when compared to PTO bulk; a trend which increases as a function of N z . Most notably, down domains feature enhanced local polar modes promoted by the internal bias field born of the compositionally broken inversion symmetry. Equally, this bias field acts to depolarise up domains leading to different critical thicknesses for the total suppression of out-ofplane ferroelectricity for the two domains. Since these critical thicknesses are the outcome of the strength and direction of the bias field, engineering this with a careful choice of substrate, surface termination or overlayers allows for control over polar textures at the nanoscale. This finding is especially important for next-generation functional devices reliant upon the control of toroidal order.
There are consequences also for the periodicity of the AFD surface reconstruction. While we find that the reconstruction for the monodomain in-plane FE and paraelectric films is c(2 × 2), we find that coupling between surface polarization and AFD modes in the polydomain geometry modifies this. In addition to surface rumpling of Pb cations, the surface reconstruction angle modulates up to 1.75 • across a domain period Λ. We then suggest that for the polydomain films, the correct label for the AFD surface reconstruction is p(2 × Λ) in Wood's notation. This provides direct evidence that the strength of AFD modes can be locally controlled by the strength of FE modes and vice versa. Such knowledge could motivate new principles of design in low dimensional functional devices whereby FE and AFD modes are locally tuned by their interactions.
We suggest that the MSSF method implemented within the CONQUEST code can now be used to solve a plethora of problems within the perovskite oxides. This could include the simulating other potentially possible polar morphologies in the PTO/STO system such as skyrmion phases [14] and disclinations [13]. Since our method is general, we can extend to other problems in the perovskite oxides (and beyond) such as realistic defect concentrations and highly disordered configurations the popular solid solution families AB 3 where DFT methods used to circumvent the need for large supercell calculations fail in the reproduction of local structural distortions [77]. This document provides supporting information for the article "Polar morphologies from first principles: PbTiO 3 films on SrTiO 3 substrates and the p(2 × Λ) surface reconstruction". We show here our tests for convergence in the parameters of our simulations. In Section 1, we show the convergence of the parameters of the multi-site support function (MSSF) method [1]. In Section 2, we study the vertical displacement of metal cations as a function of the number of included substrate monolayers. In Section 3, we study the convergence of the total energies and energy differences for the fineness of integration grid and Brillouin zone sampling. Lasty, in Section 1 we detail the method used for the calculation of local polarization. This method is used to produce the polarization vector field plots in the main text. For details about the generation of the basis sets used in this work we refer the reader to reference [2] -particularly towards its supplementary material. For further details about the CONQUEST code we refer the reader to a recently written a review article [3] and to [4] for a public release of the code under an MIT license. If any further information or raw data are required, do not hesitate to contact the authors via the email contacts at the bottom of this page † .

Multi-site support function ranges
This work uses the multi-site support function (MSSF) method [1] to contract the basis set size and thus the dimensions of the Hamiltonian matrix (along with other matrices). This significantly reduces the computational effort required for its diagonalisation. Since the computational times associated with standard matrix diagonalisation routines scale asymptotically as O(N 3 basis ) [5] (where N basis is the total number of basis functions) the MSSF contraction is particularly important. This is since the MSSF contraction typically reduces the total number of basis functions in a calculation † Jack S. Baker: jack.baker.16@ucl.ac.uk, David R. Bowler: david.bowler@ucl.ac.uk Exact diagonalisation Plane-wave (recentered) Figure 1: The equations of state E(V ) for the cubic phases of STO (left) and PTO (right) using different r M S = r LD . We also include E(V ) for the exact diagonalistion of PAOs (where the MSSF contraction has not been made) and the E(V ) (PW) result obtained from a plane-wave calculation using the same pseudopotentials. The PW curve has V 0 shifted to agree with the minimum of the exact diagonalisation calculation so a better comparison can be made of the energetics. The true lattice constants can be found in Table 1.
by a factor of 3 or 4. Such savings allow our calculations to venture into the region of ≈ 1000-3000 atoms with high accuracy.
While we refer to reader to the main text for an explanation of the method (and to references [1] and [6]) we recall that the MSSF contraction introduces two new adjustable parameters to the calculation; r M S and r LD . We must carefully converge these two parameters to ensure we retain the same accuracy as the uncontracted PAO basis (which we refer to as 'exact diagonalisation'). We do so by calculating the total energy of the high temperature cubic P m3m phases of PTO and STO as a function of volume and fitting the Birch-Murnaghan equation of state [7]. We also consider the relaxed geometry of the P 4mm ferroelectric phase of PTO. Figure 1 shows the equations of state for cubic PTO (Figure 1 right) and STO (Figure 1 left) for various r M S = r LD and exact diagonalisation. We also include on each figure the equation of state obtained with a plane-wave (PW) basis set using the same pseudopotential as performed with the ABINIT code (v8.10.3) [8,9]. Plane-wave and PAO calculations use a 6 × 6 × 6 Monkhorst-Pack mesh for Brillouin zone integrals (convergence is displayed in figure 4). Plane-wave calculations use a cut-off energy of 40Ha whilst PAO calculations use a regular integration grid (for wavefunctions and the charge density) with a plane-wave equivalent cutoff of 300Ha (convergence is displayed in figure 3). We see that by r M S = r LD = 6.350Å excellently agrees with the exact diagonalistation and PW calculation. The parameters of the fit are shown in table 1. We see once again that by r M S = r LD = 6.350Å, all the parameters of the fit closely resemble the exact diagonalisation result. We see also that the lattice constants (calculated from the EOS fit -not lattice vector optimisation which produces a slightly different result) beyond r M S = r LD = 6.350Å agree with the plane-wave calculations to an error of +0.51% and +0.52% for PTO and STO respectively. We see also that the Bulk modulus is well described by the PAO basis sets offering errors of -3.55% and -0.05% respectively. Finally, we note that all PAO calculations overestimate (by a small amount) the bulk modulus derivative B 0 = (∂B/∂P ) P =0 . This is of little impact however since the Birch-Murnaghan equation of state is not particularly sensitive to B 0 . Note, for example, that the exact diagonalisation and plane-wave curves of figure 1 are rather indistinguishable but differ by ≈ 13% in B 0 . We consider now the relaxed structure of ferroelectric P 4mm PTO. We firstly note that the structure of P 4mm PTO is unchanged between the exact diagonalisation result and the r M S = r LD = 6.350Å within the resolution of our method. That is, a static calculation using the MSSF contraction on the structure optimised with exact diagonalisation shows that our force tolerance (0.01eV/Å) is still met. Further, the phase transition energies (∆E = E P 4mm − E P m3m ) for the MSSF contracted (-47.38meV) and uncontracted (-47.40meV) basis sets differ only by 0.02meV. The phase transition energy is correct to the plane-wave result (-47.06 meV) to -0.7% showing that our method is structurally and energetically accurate.

The amount of SrTiO 3 substrate
In this section, we present a convergence study for the amount of STO substrate for the simulations in the main text. We perform our tests on the paraelectric phase of the N z = 2 film by including successively more SrO and TiO 2 monolayers in the substrate and relaxing each film (with 20Å of vacuum space) to a force tolerance of 0.01 eV/Å. To two bottom-most monolayers have fixed positions during relaxation to simulate a semi-infinite substrate. We use a uniform integration grid with a plane-wave equivalent cut-off of 300 Ha and a Monhorst-Pack Mesh [10] of dimensions 6 × 6 × 1 for momentum space integrals. We then measure the vertical displacements ∆z on the film metal cations from their initial positions and check for convergence as function of the number of substrate monolayers. As can be seen in Figure 2, Pb cations relax towards the substrate whilst Ti Cations relax away from it. This behaviour is only qualitatively achieved beyond N M L ≈ 4 with the exact positions stabilising beyond N M L ≈ 11. This suggests that N M L ≥ 11 should provide a sufficient amount of substrate for our simulations. We do, however, decide to include an additional 4 monolayers (N M L = 15) to allow further 'breathing room' for ferroelectric distortions that could penetrate into the substrate when simulating polar phases. Although we find that significant polarization (of the order of the PTO bulk spontaneous polarization) penetrates only about 4-5 monolayers into the substrate, a much smaller background level does penetrate deeper into the substrate. Our inclusion of a larger amount of substrate then allows for this effect.

Integration grids and Monkhorst-Pack meshes
In this section we show convergence in the total energy (per atom, to be consistent with the main text) and energy differences for increasingly fine real-space integration grids [5] and momentum-  Figure 3: A convergence study for the the fineness of integration grid (measured in a plane-wave equivalent cut-off energy, E cut ) for the bulk P m3m and P 4mm phases of PTO. The upper and middle panels are for convergence of the total energy whilst the lower panel measures the energy difference ∆E between the phases in meV/atom.   : The convergence of reciprocal space integrals using uniform Monkhorst-Pack meshes for the bulk phases of PTO and the N z = 2 film (detailed in the main article). The upper panels show the convergence of the total energy with increasing grid dimensions. The middle panel shows the same test but as a difference with the highly converged (to ≈ 5 × 10 −7 eV in the total energy) 12 × 12 × 12 mesh. The lower panel also compares with the 12 × 12 × 12 mesh, but measures the energy difference between the considered phases. space Monkhorst-pack [10] meshes. Figure 3 shows energy convergence for the fineness of integration grid for the bulk phases of PTO. We see that our choice of a 300Ha cut-off is correct to the (very fine) 1200Ha cut-off grid to 0.02 meV/atom in the energy difference (lower panel). Figure 4 shows the convergence in the total energy and energy differences with respect to the fineness of Monkhorst-Pack mesh for the bulk phases of PTO (left column) and the N z = 2 film (right column, geometries explained in the main text). It is interesting to point out that energy convergence for the film geometry (a 2D mesh of N × N × 1) is achieved at a coarser mesh than the bulk geometry (a 3D mesh of N × N × N ). We suggest that this is because convergence in the the out-of-plane direction for the film geometry is (close to) perfectly achieved with a single k-point as ensured by the large vacuum region separating periodic images of the film. We note that there is a small effect related to energy differences between Γ-centered and non-Γ-centred meshes. For the film geometry, wee see that the energy difference is lower by 0.19meV/atom for the 6 × 6 × 1 mesh (Figure 4, bottom right) when one phase is Γ-centred and the other is not. This effect reduces to a difference of 0.1meV/atom if we consider the energy difference between a paraelectric film (not centred on Γ) and a film with AFD modes (which is centered on Γ since we use a 6/N x × 6/N y × 1 mesh and N x = N y = 2). We therefore conclude that energy differences between geometries without AFD modes and with AFD modes will all be lower in energy by 0.1meV/atom. This extra degree of energy lowering should be taken into accoutn when examining Figure 3 in the main text, but, makes no difference to the conclusions drawn in the main text (i.e, the favorability hierarchy of different phases).

The local polarization
We define the local polarization P (i) in the manner suggested by Meyer and Vanderbilt [11] for local unit cell i, local unit cell volume Ω c , atom index α, Born effective charge tensor Z * α , local atomic displacements from the idealised cubic bulk positions u (i) α and atomic weight factor w α which is related to how many atoms α belong to the local cell i. We note that whilst the concept of a local polarization itself is rather troublesome (since the electrical polarization is a macroscopic quantity), equation 1 is a useful probe of the local mode and the resulting dipole moment. We must carefully note, however, that there are indeed several choices for the local unit cell i. These choices are discussed in the work of Meyer and Vanderbilt [11]. We chose to work with two definitions, the Ti-centered unit cell (where Pb or Sr is at the origin) and the Pb/Sr-centered unit cell (where Ti is at the origin). This allows for the calculation of local polarization vectors centered at the both the metal cation sites. For the Ti-centered cells, we have the atomic weight factors w Ti = 1, w Pb = 1/8 and w O = 1/2. For the Pb-centered cells we have w Ti = 1/8, w Pb = 1 and w O = 1/2.
The Born effective charge tensors were calculated for cubic bulk PTO and STO using finite differences in the macroscopic polarization. We do so by displacing each of the symmetry inequivalent sites within a 5 × 5 × 5 supercell (625 atoms) of PTO and STO a small amount (0.005Å) and calculating the resulting polarization with Resta's method [12]. For consistency, we use the r MS = r LD = 6.35Å MSSF contraction in the calculation. For PTO, the 3x3 symmetrical diagonal tensor has elements Z * Pb, PTO = 3.89 and Z * Ti, PTO = 7.08 whilst Oxygen is We note also the the elements of the Oxygen tensors reorder based on which site in the local unit cell we are considering. A drawback of this method is that a local unit cell cannot always be found. That is, at the film surface, there will always be an incomplete unit cell. For this reason, our calculations always finish with a vector centered on Ti and not Pb, despite treating the PbO termination. At the PTO/STO interface, there is also a hybrid unit cell with half its A-sites being Pb and half being Sr. For this area, we define no local polarization vectors.