From Slater to Mott physics by epitaxially engineering electronic correlations in oxide interfaces

Using spin-assisted ab initio random structure searches, we explore an exhaustive quantum phase diagram of archetypal interfaced Mott insulators, i.e. lanthanum-iron and lanthanum-titanium oxides. In particular, we report that the charge transfer induced by the interfacial electronic reconstruction stabilises a high-spin ferrous Fe2+ state. We provide a pathway to control the strength of correlation in this electronic state by tuning the epitaxial strain, yielding a manifold of quantum electronic phases, i.e. Mott-Hubbard, charge transfer and Slater insulating states. Furthermore, we report that the electronic correlations are closely related to the structural oxygen octahedral rotations, whose control is able to stabilise the low-spin state of Fe2+ at low pressure previously observed only under the extreme high pressure conditions in the Earth’s lower mantle. Thus, we provide avenues for magnetic switching via THz radiations which have crucial implications for next generation of spintronics technologies.


INTRODUCTION
Iron compounds and minerals are more often encountered in the ferric (Fe 3+ ) than ferrous (Fe 2+ ) oxidation state. Under normal conditions the latter is generally observed in its high-spin (HS) (t 4 2g e 2 g ) state, indeed the distribution of electrons among 3d orbitals for six-fold coordinated Fe 2+ results in a high-spin configuration due to Hund's rule and high-spin-pairing energy. However, established geological models [1][2][3] propose that the extreme pressures deep inside Earth's mantle lead to the collapse of the atomic orbitals of iron from the high-spin to the low-spin state. In particular, for the case of ferropericlase under increasing pressure, a transition towards diamagnetic Fe 2+ is expected in the pressure range of 40-55 GPa [4][5][6] . Although the latter transition has been widely studied in minerals, both through measurements and ab initio calculations 7,8 , a realisation of such a transition in synthetic compounds remain unclear. Recent progress in the field of pulsed laser deposition 9 for oxides nano-engineering has opened possibilities to study the magnetic transitions of Fe in a controlled fashion. In particular, straining interfaces in oxide heterostructures provides a fertile ground for emergent properties such as ferroelectricity 10,11 , high-temperature superconductivity 12 , piezoelectricity 13 , magnetoresistance 14 , structural reconstructions 15 , multiferroicity 16 and charge transfer (CT) 17,18 . CT is one important consequence of electronic reconstructions in oxide interfaces [19][20][21][22] . Indeed, in these systems, the structural and electronic continuity defines a set of band alignment rules, which yield electronic transfer from the high to the low-energy bands. In correlated systems, where the single-electron picture does not hold, CT is closely connected to the strength of electronic correlations. A classification of oxide materials has been provided early on by Zaanen et al. 23 , which distinguishes materials in classes of Mott-Hubbard, CT, and Slater-insulators. The three phases are characterised by different electronic and magnetic states and transitions across them are of great technological interest. Of particular significance in these materials are the low-lying metastable states, manifesting exotic properties attributed to their many-body interactions. Exploring this phase space is paramount to fully classifying the practical importance of these materials. In this study, we report an exhaustive set of exceptional and unprecedented electronic and structural properties in an archetypal oxide superlattice comprising single interfaces of lanthanum-iron and lanthanum-titanium oxides. Using a spinassisted ab initio random structure algorithm 24,25 , we predict a multitude of metastable states with their relative likelihood, and categorise them according to their contrasting magnetic features. We identify a robust antiferromagnetic ferrous ground state that can withstand significant amounts of biaxial strain that is influenced by the strength of many-body interactions, and a Slater to Mott transition under strain is revealed. Furthermore, we unveil spin state transitions of the ferrous Fe that are induced at low pressure by octahedral rotations in the THz regime, analogous to what is found at high pressure in the centre of the Earth.

RESULTS
Epitaxial engineering of robust high-spin ferrous oxides In Fig. 1a, we show a schematic band diagram of the electronic reconstruction at the LaTiO 3 /LaFeO 3 interface. Considering only the bulk counterparts, we note that Ti and Fe are both 3+ and HS configuration, with d 1 and d 5 filling, respectively. When the LaTiO 3 /LaFeO 3 interface is formed, band alignment of the oxygen states leads to an effective CT of one electron from Ti (3d) to Fe (3d) states. Thus, the CT fixes Ti 4+ and Fe 2+ oxidation states in the superlattice, where the ferrous iron has the freedom to occupy either a high-or low-spin configuration.
Using spin-assisted random structure searches we can explore the phase space of the accessible (ground and metastable) magnetic states. We use the AIRSS package 24,25 interfaced with Quantum Espresso to execute the search whose primary constraint is the starting spin configuration. As shown in Fig. 1b, we obtained a classification of the possible spin-states into four categories. The most likely candidates found have an antiferromagnetic ordering, occurring in 57% of cases and are followed by the ferromagnetic structures at 0.07 eV higher in energy, which 1 Theory where a manifold of M ππ and M 00 like spin-density waves (SDW) emerge at 0.2 eV. We emphasise that such SDW metastable configurations, as revealed by our statistical analysis, are not often mentioned in the literature, where the usually studied phases are AFM G-, C-, A-types. Finally, at 0.7 eV the phase space is dominated by states with only a PCT occurring, wherein residual moments on one of the Ti ions give rise to a more intricate magnetic ordering. Thus, the spin-assisted random structure search outlines the ground state of the heterostructure, with CT and Fe2+ HS, and highlights the importance of the low-lying ferromagnetic configurations, being the most likely metastable state for LaTiO 3 /LaFeO 3 with a different magnetic order, and it ultimately provides a rich phase diagram (see Supplementary Discussion for further analysis).
In Fig. 1c, we characterise the orbital resolved density of states of the Fe ion. The electron transferred from Ti (3d) to Fe (3d) sketched in Fig. 1a forms a localised Fe t 2g singlet at the top of the valence band. Near the Fermi level (minority, −0.5 eV) there is a strongly localised magnetic moment, induced by a spin blocking effect, where the electronic correlations localise the singlet state by a process where minority spin electron (which has been transferred from Ti to Fe) cannot tunnel to any neighbouring Fe sites, due to their large anti-parallel magnetic moment. Indeed, the minority Fe (3d) spin is localised by the saturated magnetic moments of its neighbours, a specific feature of the antiferromagnetic ground state in contrast to additional magnetic metastable states (see Supplementary Fig. S1). The optimised HS configuration can be determined through an energetic competition between the crystal field, super-exchange and Hund's coupling contributions. For an octahedrally coordinated ferrous iron atom coupled to its nearest neighbour in a crystalline lattice, we can consider all density-density interactions U = V + 2J, crystal field splitting Δ and nearest neighbour super-exchange coupling. We obtain the energy difference between the low and HS configuration: Using typical values of, e.g. Δ~2 eV, J~0.45 eV, U~4.8 eV and t~2 eV, 19 the balance between magnetic interactions and crystal field splitting stabilises the high-spin configuration. Thus, the superlattice inherits the magnetic ground state of the parent LaFeO 3 constituent, as the in-plane super-exchange coupling is protected upon the formation of the (1/1) interface. We emphasise that these features persist also in a (2/2) interface (see Supplementary  Fig. S11).
We demonstrate in Fig. 1d that through the process of epitaxial engineering the Fe 2+ HS configuration is particularly robust to inplane strain. Throughout a range of realistic epitaxial strain amounts (from −5% to +5%), there is a consistent ≈1-3 eV energy difference between magnetic and nonmagnetic configurations. Using a spin-assisted random structure search (see Supplementary  Fig. S6) several new phases are discovered, the majority being ferromagnetic that lie between 1 and 10 meV above the antiferromagnetic phase, as shown in Fig. 1d. We notice that, as the in-plane states are increasingly localised through tensile strain, then the relative energy difference between the ferromagnetic and antiferromagnetic configurations becomes smaller, which is typical of the onset of a Mott transition. In the opposite limit, where the in-plane states are compressed, the states are far in energy, indicative of itinerant magnetism (see Supplementary Fig.  S6c). Additionally, two new phases are predicted and illustrated in Fig. 1c, the first with no charge transfer (NCT) and the other with a partial charge transfer (PCT), both lying about 1 eV above the ground state. We note that the spin-assisted random structure search always finds structures with a final configuration that is HS thus verifying the robustness of the antiferromagnetic ground state, finding that it comprises~74% of all predicted structures.

Ionic radii and correlation strength
We further extend the analysis of the phenomena stemming from the epitaxial strain, focusing on the strong coupling between electronic and lattice degrees of freedom. We study as a key structural parameter the ionic radii. We estimate the atomic radii by averaging over the equatorial and apical bonds of the TiO 6 and FeO 6 octahedrals of all the compounds shown in Fig. 1d. Note that we investigate both the ground state and the metastable configurations. Interestingly we notice in Fig. 2a 26 . Similarly all the structures fall close to their nominal atomic configuration.
Outlier cases emerge in the limits of epitaxial compression (YAlO 3 , LuAlO 3 ) and expansion (LaLuO 3 ) in the CT HS configuration, where the ionic radius for Fe has a large effective value for YAlO 3 and LuAlO 3 , whereas the ionic radius of Ti increases for LaLuO 3 . Since the obtained ionic radii are related to the volumes, in Fig. 2b we analyse its trend as a function of the selected substrates. We notice that a decrease of the lattice constant is coupled to a volume collapse of LaTiO 3 /LaFeO 3 with CT HS. However, if a < 3.79 Å, the ground state volume does not follow the expected trend (dashed line) but instead expands abruptly, albeit within the stability range of the perovskites (see Supplementary Fig. S8a). A close comparison between the energies of the relaxed configurations and the unrelaxed structure corresponding to the extrapolated volume reveals that volume expansion together with increasing of octahedral rotations (see Supplementary Fig. S8c) favour the lowest energy configurations for LaTiO 3 / LaFeO 3 constrained on LuAlO 3 and YAlO 3 . In the opposite limit of strain, when the structure is clamped on LaLuO 3 the large tensile strain triggers an asymmetric in-plane extension of the Ti-O bonds, which allows the positively charged Ti to effectively screen its neighbouring oxygen anions. Further analysis of the O-Ti-O angles (see Supplementary Fig. S8b) reveals that for substrates with a > 3.905 Å, Ti prefers a tetrahedral coordination.
We extend the analysis by considering their density of states. The fractional density of states in Fig. 2c shows that if the heterostructure is clamped on LaLuO 3 , the insulating gap is between d(Ti) and d(Fe), hence the configuration is Mott-Hubbard type. As the in-plane lattice constant is reduced to the LaAlO 3 one, the band distribution at the Fermi level changes. The Fe (3d) bands, albeit keeping the localised character as in LaLuO 3 , are more hybridised with the O (2p), resulting in an overall CT Mott insulator. As the in-plane parameter is further reduced (up to LuAlO 3 substrate), the Fe (3d) states show a broader distribution and a larger hybridisation with the oxygen states. A further clarification on the nature of the heterostructure clamped on LuAlO 3 is provided by looking at the fractional density of states of the ferromagnetic configuration in Fig. 2d. We observe metastable metallic behaviour, separated by 1-10 meV from the antiferromagnetic states which are instead insulators. Furthermore, we unveil a metal-to-insulator transition depending on the amount of epitaxial strain. We have observed that the superlattice uses a cooperative mechanism between magnetism and structural deformation (volume expansion, increased octahedral rotation and enhanced La polarisation) to lower the energy and stabilise an insulator with itinerant magnetism, bearing similarities with the characteristic features of a Slater insulator. Therefore, the CT HS region in the ionic radii phase diagram is characterised by the transition between three different magnetic configurations: from the Slater insulator-like behaviour stemming from the itinerant magnetism to a more localised Mott CT insulator with an O(2p)-Ti (3d) band gap and finally up to a Mott-Hubbard insulating d − d gap. Our results confirm the expected trend that the effect of correlation increases as we go from compressive to tensile epitaxial strain, due to the reduction of the bandwidth. Notably, there is also a simultaneous structural deformation associated with the local coordination of both La and Ti (see Supplementary Fig.  S8). Moreover, we find that in the Mott phase at room temperature, via paramagnetic dynamical mean-field theory (DMFT) calculations, the CT, local magnetic moment, and spectral weight remain consistent with the zero temperature DFT magnetic calculations (see Supplementary Figs. S9 and S10).

Pressure induced spin-state crossover
We have highlighted that the ground state has a Fe 2+ HS configuration, while the LS one represents a metastable state for all of the structures analysed. An unexplored question that naturally follows is whether a volume-induced spin crossover can be observed in the heterostructure at hand. The spin crossover in iron has been widely studied in the context of geochemical modelling of the Earth's deep interior. In particular, crystal field theory and band theory predicted for ferropericlase and ironbearing magnesium silicate perovskites that a HS to LS transition occurs as a result of compression [1][2][3] . Although the latter transition has been widely studied in minerals, a realisation of such a transition in deposited compounds remains unclear. Based on similar studies of the structural-induced switching of spin-state ordering [27][28][29] , we analyse first the stability of the HS state in absence of substrate under the uniaxial compressive strain along the c-axis.
For the current analysis, there are significant changes in the volume and pressure of the structures studied and energy no longer becomes the thermodynamical quantity of choice, but instead the enthalpy. Thus, in Fig. 3a we show the enthalpies associated with the configurations obtained under the uniaxial strain. We predict for increased compression and consequent volume reduction that the enthalpy difference between the LS and the HS configuration decreases, which suggests the progressive destabilization of the latter. In particular, we report that the HS state remains robust until very large pressure. At a pressure of 60 GPa, the system undergoes a transition towards the LS configuration, as the LS state accommodates smaller volumes. Furthermore, as a consequence of the uniaxial compression, the octahedral distortions change and we notice that the high to LS transition is assisted by a quench of the in-plane octahedral rotations. This observation reveals that octahedral distortions play a pivotal role in stabilising a LS configuration in the heterostructure.
We proceed with the comparative study of the enthalpy of structures with and without octahedral rotations shown in Fig. 3b and we consider uniaxial compression (along the c-axis) of the LaTiO 3 /LaFeO 3 constrained on LaAlO 3 substrate. Firstly, we focus on the low symmetry structure. Despite the introduction of further in-plane pressure, due to the reduction of a from 4.01 to 3.79 Å, the configuration with Fe 2+ HS remains the most stable one in a pressure range 0-40 GPa. A more interesting picture emerges if we isolate the effect of inducing octahedral rotations. Hence, we suppress all octahedral rotations, promoting the starting ground state to a structure in higher symmetry, with non-zero forces only along the z-axis. The configuration with LS is the most stable structure even at small pressure values.

Terahertz-induced spin-state transition mediated by octahedral distortions
Our results so far indicate that the spin transition is weakly affected by the c-axis compression, albeit particularly sensitive to octahedral tilts. We now focus on probing the spin-state dependence on octahedral rotations only, in a manner that is independent of strain or pressure and accessible with THz light pumps. Our motivation is linked to advances in state-of-the-art epitaxial growth techniques which have allowed for the atomic control of octahedral rotations, particularly in their capabilities to suppress them 30,31 . Supplementary to these advances, ultrafast optical spectroscopy can now probe selective low-energy phonon modes by a THz pump 32 . We demonstrate that through the manipulation of octahedral rotations in the LaTiO 3 /LaFeO 3 superlattice, a spin-state transition can be achieved on the LaAlO 3 substrate. We analyse the effect of octahedral rotations by constructing a linearly interpolated isovolumetric pathway between the high symmetry undistorted superlattice and low symmetry distorted structures, with in-plane lattice parameters constrained to a = 3.79 Å. At each point along the path, the lowest energy spin-state configuration is chosen such that the energy is given relative to the high-spin distorted structure, as shown in Fig. 4a. α is the Fe-O-Fe angle and is a measure of octahedral rotation in the superlattice where for an induced tilt of ∼6°the spin-state transition occurs. A high energy THz pump can induce lattice vibrations on scales of fractions of an Ångstrom and this corresponds to the appropriate range defined by α. With this in mind, we explore the phonon spectrum, achieved by doing Γ point DFPT calculations, to identify vibrational modes within the THz regime that are compatible with octahedral rotations. The phonon density of states, shown in Fig. 4b, is distinctly split between the strongly coupled heavier anions and cations and lighter oxygen ions with a phase boundary at 6 THz. We identify a representative A g phonon candidate present at 7.42 THz, illustrated in Fig. 4c, that is similar to the tilts induced in Fig. 4a  and of primarily oxygen character only. Moreover, there are additional modes beyond 7.42 THz that extend into the Raman spectrum and beyond the range of THz spectroscopy. In summary, while the spin-state crossover is hard to achieve through uniaxial compression, since the large THz radiation needed would influence structural, dynamic and magnetic stability of the system, a more accessible pathway is provided by manipulating the octahedral rotations which involves transitions at a lower energy scale.

DISCUSSION
We have shown that a remarkably robust HS antiferromagnetic ferrous oxide emerges in heterostructure due to interfacial electronic reconstruction. In particular, we illustrate that through epitaxial engineering the electronic correlation strength is tunable and causes a rich phase diagram to arise, with a Mott-to Slater-like transition, predicted to remain stable at room temperature. The realisation of this effect has direct technology relevance for switchable oxide devices at the nanoscale. Moreover, we show that a magnetic phase transition in the superlattice can be triggered by octahedral rotations for which geological scales of pressure are not required. As such, we identify a series of compatible phonon modes at the THz scale in the vibrational spectra that can excite these rotations that are capable of activating the magnetic phase transition. From our results, it is clear that a plethora of fascinating properties emerge beyond what can be attributed solely to electronic interfacial reconstruction alone when additional tunable degrees of freedom are manipulated at the interface of two Mott insulators.

Calculation details
All DFT calculations were performed using the plane wave code Quantum Espresso 33 , version 6.4.1, together with the GGA-PBE exchange-correlation functional 34 . Atomic cores were treated within the ultrasoft, nonlinear core correction approach 35 with valence configurations La(4f5s 5p5d6s6p), Ti (3s3p3d4s), Fe(4p4s3d3p3s) and O(2s2p). The plane wave basis representation is used for the wavefunctions, with a cutoff of 800 eV. Results are converged by sampling a 6 × 6 × 4Γ-centred k-point mesh in the Brillouin zone. We apply U = 4.8 eV on Fe and U = 3.0 eV on Ti within the DFT + U approximation, where we specify the effective value of the Coulomb interaction U eff = U − J 36 as implemented in the QUANTUM Espresso package 37 . Our values reflect those commonly used in the literature for the study of a similar system 19 . Both LaTiO 3 and LaFeO 3 belong to the Pbmn spacegroup, with a − a − c + octahedral distortions in Glazer notation. We therefore constructed the (1/1) heterostructure starting from unit cell of LaTiO 3 with 20 atoms which has ffiffi ffi 2 p a pc ffiffi ffi 2 p a pc 2c pc and then substituting two in-plane Ti ions with two Fe. The heterostructure is then fully relaxed into its (1/1) ground state configuration. We further investigated the stability of the ground state obtained for the unstrained heterostructure, extending the calculation where we apply both uniaxial and epitaxial strain. The former is achieved by varying the c-axis of the heterostructure unit cell and allowing for internal relaxation. Regarding the latter, we focus on epitaxial strain constraining the orthorhombic in-plane lattice a and b to the pseudocubic parameter a sub of a set of different substrates, while the out-of-plane lattice parameter c is free to relax. Therefore, the initial structure for all the simulations with different substrates, characterised by in-plane lattice constant a sub , is represented by a 20 atoms unit cell (with 12 oxygen, 2 Fe, 2 Ti and 4 La ions) with the Gd-orthoferrite distortions inherited from the bulk counterparts (LaTiO 3 and LaFeO 3 ) and in-plane lattice constants constrained to ffiffi ffi 2 p a sub ffiffi ffi 2 p a sub . The allowed relaxations involve internal coordinates of the atoms and the lattice vector along z-axis. We clarify that we did not directly include the substrates in the calculations and as a result, we are dealing with only a single interface between two perovskites. Structural degrees of freedom are relaxed until all forces are smaller than 1 mRyd/a.u. In our calculation we consider different initial magnetic conditions for the superlattice in analysis as shown in Supplementary Fig. S1. We calculated the phonon frequencies in a THz range, using a finite-difference method 38 . We use the AIRSS package 24,25 interfaced with Quantum Espresso to execute the spinassisted random structure search. In doing so, we constrain the starting spin configurations and in-plane lattice parameters to a range of epitaxial substrates. We perform one-shot DFT + DMFT calculations at T = 290 K = 1/40 eV. The low-energy tight-binding Hamiltonian is constructed in the basis of maximally localised Wannier functions using the Wannier90 code 39 . The DMFT self-consistency cycle is implemented using the TRIQS/ DFTTools libraries 40 . For different symmetry inequivalent correlated sites the effective impurity problem is solved using the continuous-time Quantum Monte Carlo hybridisation-expansion solver as implemented in TRIQS/CTHYB 41 . The full frequency spectral function A(ω) is obtained via analytically continuing the Matsubara Green's function using The TRIQS/ MaxEnt library 42 . All results are performed within the "one-shot" DFT + DMFT approximation and neglect the effect of charge self-consistency. Further details can be found in the Supplementary Discussion.

DATA AVAILABILITY
The datasets generated and/or analysed during the current study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The code used during the current study are available from the corresponding author upon reasonable request. a b c A g A g Fig. 4 Terahertz-induced spin-state transition mediated by octahedral distortions. a Magnetic moment per Fe (Relative energỹ E ¼ E distortion À E LaAlO3 ) of distorted structures shown on the left (right) axis, as a function of octahedral rotations α. b Atomically resolved phonon density of states for the structure with π − α = 27 ∘ . c A g vibrational mode at 7.42 THz which is compatible with octahedral rotations shown in a.