Understanding the metal-to-insulator transition in La1−xSrxCoO3−δ and its applications for neuromorphic computing

Transition metal oxides that exhibit a metal-to-insulator transition (MIT) as a function of oxygen vacancy concentration are promising systems to realize energy-efficient platforms for neuromorphic computing. However, the current lack of understanding of the microscopic mechanism driving the MIT hinders the realization of effective and stable devices. Here we investigate defective cobaltites and we unravel the structural, electronic, and magnetic changes responsible for the MIT when oxygen vacancies are introduced in the material. We show that, contrary to accepted views, cooperative structural distortions instead of local bonding changes are responsible for the MIT, and we describe the subtle interdependence of structural and magnetic transitions. Finally, we present a model, based on first principles, to predict the required electric bias to drive the transition, showing good agreement with available measurements and providing a paradigm to establish design rules for low-energy cost devices.


INTRODUCTION
The search for computing architectures with low-power consumption is an active field of research, and neuromorphic architectures, which aspire to mimic the human brain 1 , have attracted much attention lately. The realization of neuromorphic devices imposes specific requirements on the materials to be used: the response to an input signal should be dependent on the signal's intensity, frequency and the previous status of the material, in a way similar to neurons and synapses in the brain, which combine computation and memory in one unit. In addition, the energy consumption to transmit and process signals should be extremely low to be affordable on a large scale; remarkably there are~10 11 neurons and~10 14 synapses in the brain, with a power consumption of~10 W for daily cognitive tasks 2 .
In the last decades, several transition metal oxides (TMO) have been proposed as promising resistive switching materials 3,4 , i.e. systems showing tunable resistance states, induced by an external electrical bias. These TMO exhibit a metal-to-insulator transition (MIT) as a function of pressure, temperature, or doping, which may be designed to mimic the behavior of neurons and synapses in the presence of stimuli [5][6][7] .
In spite of extensive work, the mechanisms underlying the MIT in cobaltites remain elusive, thus hampering the control of the properties of these materials to design efficient neuromorphic devices. Here we report an investigation of the MIT in LSCO (and the semiconductor-to-insulator transition in LCO) using first principles calculations; we unravel the complex interplay between structural, electronic, and magnetic properties of the material as the transition is driven from a perovskite to a BM phase, as a function of oxygen vacancies. Our calculations show that cooperative, global structural distortions leading to a decrease of the elastic energy of the solid, which in turn are accompanied by changes in the magnetic state of the material, are the main factors driving the MIT. Our results allow for the identification of general descriptors to design materials for neuromorphic computing applications, including multiple resistive states. Finally, we present a model based on first principles to predict the electric bias necessary to drive the MIT, which shows good agreement with experiment, and is general and applicable to broad classes of TMO.

RESULTS
We used a combination of electronic structure techniques based on density functional theory (DFT) with a Hubbard correction (U). In order to choose the most appropriate energy functional and U parameter, we computed the structural, electronic, and magnetic properties of the perovskite LCO and BM SCO 2.5 phases, using U values from 3 to 7 eV and several functionals (LDA 21 , PBE 22 , and SCAN 23 ). We found that the PBE functional with U = 3 eV represents the optimal choice to correctly describe most measured properties of both phases, and we adopted the same protocol to study the topotactic transition as a function of oxygen vacancies. (Supplementary Discussion 1 contains a detailed study  leading to our choice of the DFT parameters, see Supplementary  Fig. 1 and Supplementary Tables 1 and 2).

Structural transitions
We start by analyzing the structural distortions associated with the topotactic transformation between perovskite and BM phases in LCO and LSCO, as a function of oxygen vacancy concentrations (V O , including 4.2, 8.3, and 12.5%, where V O = δ/3 × 100%). We carried out calculations in a supercell with four layers of octahedral units (Fig. 1), and we considered symmetryinequivalent oxygen vacancy positions (namely 2 configurations at 4.2% V O , 12 at 8.3% V O , and 1 at 12.5% V O ), and different magnetic states at each concentration (ferromagnetic, FM, and several AFM states including G-AFM, A-AFM, and C-AFM) (see Supplementary Figs. 2 and 3). We found that the lowest-energy structure corresponds to ordered vacancies belonging to the same layer (see Supplementary Fig. 4) and we observed a transformation from octahedral CoO 6 units to pyramidal CoO 5 and then to tetrahedral CoO 4 ones within a layer; such transformations then occurred in alternating layers, eventually leading to the BM phase. Our results are consistent with previous in situ TEM observations, which indicated that the perovskite to BM transition occurs as a nucleation process of oxygen vacancies from layer to layer, in both LCO 24 and La 2/3 Sr 1/3 MnO 3 (ref. 25 ). We also found that the perovskite lattice expands with increasing vacancy concentration, as shown in Fig. 2a; the same trend was observed in various perovskite materials [26][27][28] . Our calculations show that upon Sr doping, the lattice expansion is less pronounced, consistent with experimental data 29 reported for 50% Sr doping (green line in Fig.  2a). Indeed, Sr doping leads to an increase of the Co oxidation state, counterbalancing, at least in part, the effect of the vacancy's extra electrons. As a result, the Co 3d electrons mainly occupy the t 2g manifold (see Supplementary Fig. 5 for energy diagram) and hence the Jahn-Teller (JT) distortion following the introduction of vacancies is smaller compared to LCO, where the e g band is occupied. Another notable change revealed by our calculations occurs in the Co-O-Co angle, which substantially decreases, as a function of V O (Fig. 2b). We expect these results to be only weakly dependent on the distribution of Sr dopants in the crystal. For example, we computed the structural, magnetic, and electronic structure properties of non-defective LSCO samples as a function of Sr and we found that the standard deviation of the calculated lattice parameters for different Sr positions is less than 0.5% (see Supplementary Figs. 6 and 7 and Supplementary Table 3); we also found that all samples are FM metals, with similar magnetic moments of the Co atoms (see Supplementary Fig. 8).
Following the analysis of ref. 30 we characterized the distortions of the perovskite phase introduced by oxygen vacancies as a combination of JT distortions, a breathing mode of the octahedral units and a change of the Co-O-Co tilt angle, as shown in Fig. 3. The combined JT/breathing distortion of the octahedra is characterized by an increase of in-plane (ac plane in Fig. 1) and out-of-plane (along the b-axis in Fig. 1) Co-O bond lengths, with the latter being more pronounced (see Supplementary Tables 4  and 5 and Supplementary Fig. 9 for lattice parameters change and bond length and bond angle change).
The Co-O-Co angle variation, or rotation distortion, is reflected in the decrease of the angle in the out-of-plane direction, more . The introduction of one oxygen vacancy in the 40 atom perovskite cell (4.2% concentration) considered in our calculations leads to the transformation of an octahedral CoO 6 unit into a pyramidal CoO 5 unit in a given layer (layer 4); when a second oxygen vacancy is introduced (8.3% concentration) the CoO 5 units are transformed into tetrahedral CoO 4 ones. Upon introduction of three (12.5%) and four (16.7%) vacancies, an octahedral-to-pyramidal-totetrahedral transformation is observed in alternating layers (layer 2) and eventually a topotactic transformation to a brownmillerite (BM) phase is achieved. b Top-down view along the b-axis of the perovskite and BM structures. In the BM structure, octahedra were removed to clearly show the rotation pattern of tetrahedra in alternating layers of the defective perovskite that led to the BM formation. The BM has Pmnb crystal symmetry 58 as the tetrahedra in alternating layers (layers 2 and 4) have opposite handed rotation (shown by curved arrows). prominent than in other directions. Such a decrease is caused both by the rotation of octahedra and by the Co ion displacement observed in defective polygons.
While our results on JT/breathing distortions are validated, at least qualitatively, by the volume expansion found in experiments 15,29 , there are no experimental measurements of the Co-O-Co tilt angle. However, there are indirect effects of bond angle changes that can be measured, e.g., the atomic distance change between cations (La/Sr) in the out-of-plane direction. When octahedra rotates, La/Sr moves closer to its neighboring oxygen ion, thus increasing their covalent bonding strength 31 and shortening the intra-layer La(Sr)-La separation. We found that such separation decreases (in even number rows in our cell: see Supplementary Fig. 10), in good agreement with the experimental trend reported for LCO 24 .
In order to further understand the physical picture underlying the distortions found in our calculations, we modeled the atomic deformations of the perovskite using a Hamiltonian written as a sum of three terms: In Eq. (1) the index r denotes the position of each octahedron in the unit cell; a T , a D , and a R are the stiffness of the Jahn-Teller, breathing, and rotation distortions and T r , D r , and R r are the respective deformation magnitudes. We assumed a T , a D , and a R to be site independent; they were computed as averages over all the octahedra/polygons in the unit cell. The Hamiltonian in Eq. (1) is written in the harmonic approximation. Octahedral units in layers adjacent to oxygen-deficient layers (vacancy position marked by green circle) expand as a combination of Jahn-Teller and breathing mode distortions (expansion indicated by arrows within the octahedra). The octahedra rotate to accommodate the expansion with a rotation angle Φ 0 . Oxygen ions connecting the octahedral and oxygen-deficient layers move as a result of rotations, shown by the arrows on these oxygen ions. Co ion shifts to the center of the oxygen-deficient polygons, and a slight JT expansion in the in-plane direction occurs in the oxygendeficient polygons. Red spheres: oxygen ions, blue spheres: Co ions.

Fig. 4
Computed stiffness E (meV Å −2 ) of Jahn-Teller distortions in momentum space (k) with cooperative distortion effects included. a With rotation angle Φ 0 = 0°and b With rotation angle Φ 0 = 10°in perovskite and BM LCO. Φ 0 is the octahedral rotation angle defined in Fig. 3, and it is equal to the deviation of the Co-O-Co tilt angle (Fig. 2) from 180°divided by two. The color scale represents the scaled value (E/a T : a T is the Jahn-Teller stiffness parameter a T computed for the perovskite phase; see text) in the perovskite (left panels) and BM phases (right panels). Note the different scale used for Φ 0 = 0°and Φ 0 = 10°, which shows the notable decrease in energy for Φ 0 = 10°. Such decrease is more significant in the brownmillerite than in the perovskite case. L is the octahedral lattice spacing.
We first estimated a i from the average vibrational frequencies of the corresponding mode i which we obtained by carrying out ab initio phonon calculations (using the Q UANTUM ESPRESSO 32,33 code): Here ν is the average wave number of the ith mode and μ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p is the reduced mass; m 1 and m 2 denote the mass of Co and O atoms, respectively. We found that a T and a D are on the order of 10 4 meV Å −2 , while a R is two orders of magnitude smaller (see Supplementary Table 6), indicating that the JT/breathing distortion requires much higher energy to occur than rotations. The small energy required to trigger these rotations stems from phonon modes of very low energy, as identified in our phonon calculations as well as detected experimentally 34 . Following the protocol of ref. 30 (see Supplementary Equations 1 and 2), we rewrote Eq. (1) by applying specific structural constraints, so as to describe the interaction between the JT/breathing mode and rotations. In particular, we rewrote Eq. (1) by applying two deformation constraints: octahedra are corner-sharing and the parallelogram shape is assumed to be maintained. In doing so, we allowed for the coupling between the deformation of neighboring octahedra, and we could describe the distortions of each single octahedron using a combination of different phonon modes. By expressing the coupled distortions in momentum space, we could easily represent phonon-phonon interactions. Figure 4 shows the stiffness of the JT/breathing mode in momentum space, which includes cooperative effects and thus depends on rotation. By increasing the rotation angle Φ 0 , the stiffness of the JT/breathing mode decreases, showing that rotations can effectively counterbalance the elastic energy increase caused by JT distortions. Coupling between different distortions has been previously observed in other pervoskite systems such as manganites 35 and nickelates 36 . The stiffness shown in Fig. 4 progressively decreases from the perovskite to the BM phase, due to more facile rotations as V O increases. Such cooperative distortion effects occur also in the presence of Sr doping; however, the JT/breathing distortion is less pronounced than in LCO (see Fig.  2a), and so are rotations, as expected (see Fig. 2b).

Magnetic transitions
The structural deformations described above are accompanied by important changes in the magnetic properties of the system. With increasing V O , we observed a transition from a non-magnetic structure (LCO) or FM (LSCO) to an FM-local AFM (for V O = 4.2% and 8.3% respectively), where some of the Co atoms neighboring oxygen vacancies give rise to an AFM state, while the surrounding ones form an FM state. Eventually the system undergoes a transition to a complete AFM phase. Note that in the spin-polarized case, where one has partially filled d orbitals (full occupation of majority spin channel and partially filled minority spin channel), the superexchange interaction between Co-O-Co with almost 180°, favors AFM ordering, according to Goodenough-Kanamori-Anderson rules. In turn this ordering implies lower Co-O orbital hybridization and longer Co-O bond length, consistent with a previous study on SrCoO 3 (ref. 37 ). As shown in Fig. 5a, the absolute and total magnetic moment differences increase, consistent with the decrease in total magnetic moment observed in experiments 15 . Our results are also consistent with those of previous theoretical and experimental studies that reported different magnetic states at different vacancy concentrations 37,38 . For example, in strained LCO epitaxial film with 11.1% ordered vacancy concentration, first-principle calculations 38 predicted a local AFM structure developing next to vacancies embedded in an FM environment.
Interestingly, if we constrain the system to be non-magnetic (a configuration which is less energetically favorable than that of a magnetic state), we do not observe any volume expansion (see Fig. 5b) or any increase in the Co-O bond length, indicating that the combined JT/breathing distortions are suppressed. This result can be understood noting that the Co 3d-O 2p hybridization is higher (see Supplementary Figs. 11 and 12 of PDOS), thus leading to a shorter Co-O bond length 37 . As a consequence, the changes of angles in the layer with tetrahedral units and in the out-of-plane direction are greatly suppressed: the out-of-plane Co-O-Co angle is the same as in the perovskite phase, and the in-plane Co-O-Co angle in tetrahedral units approaches 180°(see Table 1), indicating that rotations are inhibited. If the system is instead constrained in a FM state, we observed only a slight modification of the deformation compared to that in the ground state; the volume expansion is within the error bar of our calculations at V O = 16.7%  S. Zhang and G. Galli (see Fig. 5b), and the Co-O-Co angle is just slightly larger (<3%) than in the ground state (see Table 1).

Metal-insulator transition
The structural and magnetic transitions identified as a function of increasing V 0 are responsible for the MIT (Fig. 6a) found in our calculations and observed in experiments in the case of LSCO 15 , and for the semiconductor-to-insulator transition in the case of LCO. We emphasize that the MIT is intimately connected to the magnetic state transition described above (see Fig. 5a). Note that doping LCO with Sr decreases the band gap, consistent with Cheng et al.'s 14 results obtained by combining experiments and DFT calculations; in LSCO the band gap opening is first observed at a higher V O than in LCO, and in both cases it is clearly associated with the FM-to-AFM transition point. The trends in the density of states (DOS) as a function of V O shown in Fig. 6b reveals in detail how the electronic structure of the system evolves (DOS for LSCO data is shown in Supplementary Fig. 13). Interestingly, a gap opens up first between states associated to oxygen-deficient layers. As oxygen vacancies are added to alternating layers, states originating from oxygen-deficient Co ions (i.e. from Co less than sixfold coordinated) as well as from neighboring octahedral Co are increasingly more localized. This finding is supported by the analysis of Born effective charge (see Supplementary Table 8), which shows that the charge transfer from Co to O decreases with increasing V O , and it is consistent with the observed decrease of the dielectric constant. Furthermore, we observed that the band gap increases as a function of structural deformations leading to pyramidal and then tetrahedral coordination. The main contribution to band edge states (see Supplementary Fig. 11) originates from the octahedral Co's e g manifold (majority spin channel) for V O ≤ 8.3 %; as V O increases, the main contribution comes from the pyramidal Co's e g manifold (majority spin channel) for V O = 12.5%, and from the minority spin channel of Co for V O = 16.7%, with contributions from each d orbital of either octahedral Co or tetrahedral Co. These results highlight the importance of configurational as well as valence state changes for a band gap to develop in the system. A detailed Co spin configuration analysis ( Fig. 6c and see Supplementary Table 7 and Supplementary Fig. 15 for details) reveals how the electronic structure of the system is related to structural deformations and magnetic state change. We found that the crystal field splitting of oxygen-deficient Co depends not only on the oxygen vacancy configuration but also on the spin configuration of neighboring octahedral sites. For example, the spin-up pyramidal Co site present at V O = 4.2% experiences an octahedral crystal field splitting with six valence electrons, resembling that of its neighboring octahedral Co site; the spindown pyramidal Co site exhibits instead a pyramidal crystal field splitting and seven electrons. This indicates that the effect of oxygen vacancies on the electronic structure is not a local effect. In a consistent fashion, the variation in the number of electrons belonging to the 3d orbitals of Co caused by oxygen vacancies (see Fig. 6c) occurs not only in oxygen-deficient layers but also in neighboring octahedral layers. To understand this variation in c The corresponding oxygen-deficient structures and Co spin configurations. The color of the dashed border corresponds to the color legend in b and it is used to distinguish Co with different oxygen coordination. High-spin and low-spin octahedral Co is also distinguished (see Supplementary Fig. 15 for magnetic moment of each Co). Pink arrow represents increased valence electrons in the 3d orbital caused by oxygen vacancies (due to both extra valence electrons from the oxygen vacancy as well as the decreases charge transfer from Co to O, see more explanation in the main text).
number of electrons, it is important to consider the contribution of charge from oxygen vacancies as well as the decrease of electron transfer from Co to O due to an increase in electron localization of the Co orbitals. The charge variation at oxygen-deficient Co sites is mainly due to the extra electrons, while electron localization influences the charge at octahedral sites the most, and more so for low-spin state Co than spin-polarized one. The relative importance of these effects (extra charges and decrease in electron transfer due to localization) depends on V O , and in the BM phase, both effects contribute equally to determine the electronic structure of tetrahedral and octahedral Co sites. The increase in number of electrons leads to transitions from LS to HS or IS spin states at octahedral Co sites, and triggers JT/breathing distortions in layers with octahedral sites, followed by rotations to decrease the elastic energy. The final, stable spin configurations, magnetic structures, and structural distortions are thus the result of a global redistribution of the electrons in the lowest-energy configurations of the system. This redistribution is influenced both by the oxygen vacancy configuration and by their concentration. Hence the materials electronic structure is intimately connected to the magnetic state transition and to structural deformations.
In addition, we note that in the projected DOS, a fairly large gap is present between states associated with oxygen-deficient layers, as soon as oxygen vacancies are introduced in the system, indicating that electron correlation increases at oxygen-deficient Co sites. The gap between states associated to octahedral layers decreases when a small amount of vacancies is introduced, as it is related to extra electrons occupying the e g orbitals. However, as the concentration of V O increases, the effect of correlation between electrons at octahedral Co sites increases as well, eventually leading to the opening of a gap. The presence of Sr partially compensates the electron doping effect from oxygen vacancies; hence, correlation effects are expected to be weaker and likely weakly influenced by an increase in oxygen vacancy concentration, as reflected by less significant structural distortions and smaller band gap changes in LSCO than in LCO.
To further understand how electronic correlation is affected by the magnetic state and associated structural deformations in the material, we carried out calculations by constraining LCO to be in specific magnetic states, corresponding to structural distortions different from those observed in the ground state (see Supplementary Fig. 12). We found that when the system is constrained to be non-magnetic, the band gap vanishes instead of increasing. The projected DOS shows that the band gap closure is mainly due to the energy levels arising from oxygen-deficient Co moving closer to the Fermi level. Note that in the absence of spin polarization (see Supplementary Fig. 12a), tetrahedral and octahedral Co sites exhibit the same crystal field splitting and hence one observes a decrease in electron correlation relative to the spin-polarized case. When constraining the system in an FM state, we observed instead a finite, albeit smaller band gap, and the difference compared to the AFM ground state (Fig. 6b) is mainly in the localization of the states; in particular, the band stemming from tetrahedrally coordinated Co ions is less localized than in the ground state. Overall our results indicate that due to spin polarization, especially in the presence of AFM ordering, the oxygen-deficient Co sites experience a crystal field splitting different from the one of the octahedral sites, leading to an increase in localization of the d orbitals and thus to an increase in correlation. As a result, defect states introduced by oxygen vacancies are pushed well inside the valence band of the material, thus allowing for a band gap to open.
We further analyzed the effect of structural distortions on the band gap by constraining the angle of the G-AFM BM LCO structure to be close to the value it has in the initial perovskite structure (160°of all Co-O-Co angle). Lattice expansion was allowed in order to stabilize the G-AFM state. We found that the hybridized octahedral Co 3d-O 2p energy levels moved upwards, towards the Fermi level. The crystal field splitting of tetrahedral Co sites becomes similar to the octahedral case, showing the tight link between the system's crystal field splitting and rotation distortions. These results suggest that lattice expansion alone cannot lead to an MIT, and that rotation is a key factor for the opening of the band gap. Hence in order to have a complete opening of the gap all elements identified here-volume expansion, rotation, and associated magnetic state change-are essential.
From mechanisms to devices We close by discussing the impact of our results on device applications. Experiments have demonstrated that applying electric fields 17 , e.g. by the ion-gel gating method 29,39 , is a promising way to control V O in LSCO and in other perovskite materials such as nickelate where a synaptic behavior was generated by electrical stimuli 40 ; hence it is interesting to estimate the required applied potential to modulate the vacancy concentration. The field applied to a non-metallic material (E 0 ) can be evaluated using the following equation: where E is the induced field, P polar is the permanent polarization vector (which is zero for the non-defective perovskite phase of LCO, but non zero in the presence of defects); ε 1 is the longwavelength dielectric constant of the system and Eq. (3) is derived by assuming that the ionic and electronic components of the induced polarization are approximately the same (see Supplementary Discussion 6). Indeed, we know from experiments that the zero-frequency dielectric constant of LCO is about twice as ε 141 . We computed ε 1 using density functional perturbation theory (DFPT) and P polar using the Berry phase approach. We obtained the ith component of E as E i = η/d i , where d i is a supercell unit length and η is the electrochemical potential necessary to trigger the perovskite to BM phase transition. We computed η by requiring that vacancies can form spontaneously, i.e. that the oxygen vacancy free energy formation (F vac ) equals to zero; this requirement yields the condition: η j j ¼ 1 2 Fvac ne where F vac was obtained by adding to the total energy of formation a vibrational free energy contribution computed from phonon frequencies (see Supplementary Equation 6). Figure 7a reports F vac /n (free energy per vacancy) as a function of V O for both LCO and LSCO. To estimate the electrochemical potential required for a complete perovskite to BM phase transition, we considered the highest value of F vac /n, yielding η j j = 0.6 V for LCO and 0.4 V for LSCO, respectively. We found that Sr doping generally decreases F vac /n by~0.4-0.5 eV relative to LCO, thus leading to a lower value of η j j. These results indicate that Sr-V O segregation may occur during the topotactic transformation as the V O formation energy is lowered with Sr doping; however, we do not expect that a lower oxygen vacancy formation energy will affect the MIT mechanism described earlier, since also in the case of 100% segregation (i.e., SCO 3−δ ), a similar MIT transition has been shown to occur 16,37 , accompanied by an FM-to-AFM change during the perovskite to BM topotactic transition. We also note that the energy required to create oxygen vacancies is higher than that of the barrier for oxygen diffusions (the migration enthalpy in LCO is estimated to be~0.7 eV from DFT+U calculations 42 and experiments 43 ); thus, the diffusion process is unlikely to be the bottleneck of topotactic phase transitions.
In Fig. 7b, we show the applied electric field E 0 and the corresponding voltage as a function of V O assuming the maximum value of η computed above. Applying an electric field along the out-of-plane axis b leads to the lowest polarization value in the material and thus minimizes the difference between E 0 and E. The upper limit of the voltage corresponding to E 0 averaged along the three axes is estimated to be~1.2 V for LCO (at V O = 8.3%) and 0.8 V for LSCO (at V O = 12.5%). The latter value is approximate, as it was obtained by assuming the same polarization for LSCO and LCO at 12.5% V O , and considering the value of η computed for LSCO.
The voltage estimated in our calculations is a promising value, qualitatively similar to that found experimentally (0.5 V) to trigger a lattice expansion in 50% Sr-doped LCO films by the introduction of oxygen vacancies 29 . It is a value comparable to the operation voltage (0.8 V) in VO 2 -based artificial neurons, which is considered low 44 compared to that of CMOS-based neurons. A typical CMOS component requires a minimum operating voltage around 0.5 V, while a CMOS implementation of a neuron typically requires eight or more active components, thus leading to an increase of the total energy consumption 44 . Moreover, for an onset of the MIT or a transformation between different resistive states, a complete perovskite to BM phase transition is not necessary; therefore, a small variation of V O near the transition point suffices, e.g., near V O = 8.3% in LCO, which would further decrease the required potential value by~0.14 V, due to a smaller value of F vac per vacancy. Varying the Sr concentration may also help modulate the required potential, which we find decreases the upper bound by 0.4 V for Sr = 37.5%. In fact, in the 100% Sr case (SCO x ), the relative stability of the perovskite and BM phases is reversed and a voltage of only 30 mV is found to be necessary to trigger the topotactic phase transition 17 .

DISCUSSION
Using a series of first principles calculations, we characterized, at the microscopic level, the MIT in La 1−x Sr x CoO 3−δ (x = 0 and x = 0.375) as a function of oxygen vacancy concentration (δ from 0 to 0.5). We found that the introduction of oxygen vacancies in the perovskite phase of LCO and LSCO lead to structural distortions, accompanied by a change of the magnetic state of the material and of the Co ion's oxidation state. In particular, in LCO (LSCO) we observed a non-magnetic (FM) to AFM transition during the topotactic perovskite to BM transformation. Importantly, when constraining the system to be non-magnetic, all structural distortions, including the combined JT/breathing distortions and rotations, are suppressed, thus preventing a transition to an insulating state. We identified several structural transition pathways from the perovskite to the BM, with the lowest-energy pathway well separated (~0.8 eV per vacancy lower) from other metastable pathways. We found that the combined JT/ breathing distortions lead to a localization of states near oxygendeficient Co ions and thus an increased electron correlation effect, enabling the opening of a band gap between states associated to oxygen-deficient layers; the rotations triggered to lower the elastic energy of the system are further responsible for the opening of the gap between states connecting octahedral and tetrahedral layers, eventually leading to an MIT. We also showed that the band gap may be tuned by regulating the oxygen vacancy concentration, which is a necessary prerequisite to form different resistive states. These findings are consistent with the experimentally measured resistivity increase in the material, as a function of oxygen vacancy concentration and can be further validated by future optical measurements. We emphasize that in order to have a complete opening of the gap upon introduction of vacancies, all elements identified in our calculationsvolume expansion, lattice distortions, and associated magnetic state change-are essential. The cooperative distortions predicted in our calculations, in particular the rotations of Co-oxygen units may be probed by X-ray and TEM measurements. It is interesting to note that the microscopic mechanism identified here in LSCO differ from that reported for the MIT in defective nickelates 45 . In nickelates, local electronic structure changes arising from excess electron in different crystal field splitting environments of the Ni ions can explain the band gap opening, in the absence of magnetic state changes and cooperative structural distortions.
Finally, we showed that the topotactic transformation in LCO and LSCO may be induced by an application of an electric field of the order of that found experimentally, and we presented a general model applicable to broad classes of TMOs. We note that a similar non-volatile resistive switching process may also be triggered in VO 2 (ref. 46 ), a broadly studied TMO for neuromorphic applications 3,47 . However, in this case the Fermi level is moved to the conduction band 48 as soon as vacancies are formed, and no intermediate resistive states may be created. In addition, in LSCO the formation of oxygen-deficient layers during the topotactic transformation facilitates the migration of oxygen ions 49,50 , which in turn increases the switching speed from a metal to an insulator, relative to vanadium oxide compounds.
The results obtained here, especially the interplay between structural, electronic, and magnetic properties identified in our calculations, fully characterize the MIT for the first time, and may be used to engineer specific resistive states in cobaltite films, e.g., by interfacial heterostructures 51 and by controlling strain effects. It is well known that strain plays a key role in the presence of interfaces 38,[52][53][54] ; the combined structural and magnetic effects identified in our calculations may thus be used to engineer interfaces with desired properties for neuromorphic applications in cobaltites and likely manganites as well 55 , and in general to understand topotactic transformations in several classes of TMOs.

DFT+U calculation parameters
We used DFT, the PBE generalized gradient approximation with a Hubbard U correction and the projected augmented wave PP from PSlibrary 56 , and performed calculations using Q UANTUM ESPRESSO (QE) code (V6.2 and V6.4.1) 32,33 . To obtain the ground state structures for different vacancy concentrations and La/Sr ratios, a ffiffi ffi 2 p 4 ffiffi ffi 2 p perovskite supercell (La S. Zhang and G. Galli (Sr) 8 Co 8 O 24 , 40 atoms) was used, which was then optimized for each configurations within the orthorhombic lattice symmetry, using a planewave cutoff of 90 Ry and a 6 × 3 × 6 Monkhorst-Pack k-point grid. Sr doping positions were determined by sampling all symmetry-inequivalent positions as A-site cation (see Supplementary Fig. 6).

DATA AVAILABILITY
Data that support the findings of this study will be available through the Qresp 57 curator.