Electrical polarization induced by atomically engineered compositional gradient in complex oxide solid solution

Control of inhomogeneity in materials in order to avoid unexpected effects to the system remains a challenge. In this study, we seek to engineer inhomogeneity in materials and anticipate new properties. Through precise control of composition at the atomic scale, an electrical polarization is induced in the composition-graded LaAlO3–SrTiO3 solid solution epitaxially deposited on NdGaO3 substrates. By tailoring the direction of compositional gradient, the relationship between structure and electrical polarization is simulated via phase-field modeling and revealed by a combination of scanning transmission electron microscopy and scanning probe microscopy. The analysis of the results indicates that the induced polarization is due to the flexoelectric effect in the compositional gradient system. The results of this study provide a new pathway for obtaining a new material genome. Moreover, by a suitable design of the new genome, that is, by using different combinations of compositional gradient geometries, local conduction can be obtained and manipulated, providing a new approach to obtain the desired properties. Controlling the crystal quality of oxide thin films offers a way to engineer unexpected and novel material properties. Altering the distance between atoms in a crystalline film, by growing it on a different substrate for example, induces a physical strain and alters its properties. The characteristics of the film can also be altered by varying its composition. Researchers led by Chun-Gang Duan from the East China Normal University, Shanghai, China, and Ying-Hao Chu from National Chiao Tung University, Hsinchu, Taiwan, combined these two approaches by controlling the atomic strain in a combination of strontium titanate and lanthanum aluminate. The team’s computer model indicated that this allowed the creation of a material with so-called flexoelectric properties, where the electric and elastic properties are interlinked. This was verified using electron and piezoresponse force microscopy. In this work we report the electrical polarization can be induced by composition-graded LaAlO3-SrTiO3 solid solution. The analysis indicates that the induced polarization is attributed to the flexoelectric effect in the compositional gradient system. By a suitable combination of compositional gradient in geometry, a local conduction can be created and manipulated, delivering a new approach to obtain unveiled properties. delivering a new pathway to acquire material genome based on material inhomogeneity.


Introduction
The use of new materials opens a new era for human beings. For example, silicon-based semiconductor technology introduced the era of electronic devices, revolutionizing our daily lives. Thus, exploration of new materials is very important for developing nextgeneration electronics. However, it usually takes several decades to explore and use new materials for practical applications. To discover, manufacture, and deploy advanced materials twice as fast at a fraction of the cost, the US government has started the "Materials Genome Initiative" program 1 . Phase diagrams play a very important role in the materials development process, typically serving as a cooking recipe. In these phase diagrams, the structures and phases with different orderings can be determined once the temperature and composition are identified, assisting scientists and engineers in the design of fabrication processes 2,3 . Typically, these phase diagrams are constructed based on three-dimensional homogenous compositions 4 . To tailor these phase diagrams, researchers have started to use dimensionality control to explore the evolution of phase diagrams as the dimensionality decreases from three dimensions (bulk), to two dimensions (thin films) 5 , one dimension (nanowires) 6 , and zero dimensions (nanocrystals) 7 , leading to an increasingly dominant role of the surface effects. However, due to their thermodynamic stability, the existence of inhomogeneities in materials cannot be avoided. Material inhomogeneities, such as vacancies, dislocations, domain walls, and grain boundaries, are typically viewed as defects 8 . Furthermore, the physical properties of these inhomogeneities can be very different from those of the homogenous bulk, so the properties of materials are determined by the combination of its homogenous and inhomogeneous components 9 . Therefore, inhomogeneity of materials can be viewed as a hidden degree of freedom in material design. In this study, atomic engineering of inhomogeneity in materials together with the exploration of physical properties of these inhomogeneities is demonstrated.
One key approach to manipulate the ground state of materials is to impose strain on material systems. The magnitude of the strain is defined as the difference in the lattice parameters between the strained and ground states. In an epitaxial system, the substrate can be used to generate epitaxial strain because the lattice parameters of the substrate and film are different Strain ¼ Δa a bulk ¼ a film Àa bulk a bulk [10][11][12] .
The reference point is the bulk lattice parameter, which is typically fixed when the composition is determined. However, if one can vary the composition across the film thickness, the bulk lattice parameter can be modulated according to Vegard's law 13 . Moreover, when the thickness of such a system is on the order of a few nanometers, the inplane lattice of the whole film is coherent with the substrate. Since the bulk value is changing, a large-strain gradient can be generated. In the presence of a strain gradient, it is expected that new functionalities, such as flexoelectric, flexomagnetic, and flexomagnetoelectric effects, can be generated [14][15][16] . In this study, a solid solution composed by two parent compounds, LaAlO 3 (LAO) and SrTiO 3 (STO), is adopted as a model system. To engineer the inhomogeneity on the atomic scale, a compositional gradient at the atomic level was modeled using phase-field simulations. Epitaxial growth of the STO-LAO solid solution was carried out to acquire the desired heterostructure. A detailed study based on transmission electron microscopy and scanning probe microscopy was carried out to verify the prediction of the phase-field simulations. To tailor the functionalities of the system, the local conduction can be engineered by a proper design of the system, opening a new avenue for the manipulation of the dielectric and transport properties of materials. Thus, material inhomogeneity can form a new material genome for the design and development of new functional materials for next-generation technology.

Results and discussion
We first developed a phase-field model to examine the possibility of inducing polarization by chemical and strain gradients in the compositionally graded STO-LAO solid solution. STO is modeled as an incipient ferroelectric material that remains paraelectric in bulk at room temperature. Notably, ferroelectricity may be induced by misfit strain due to substrate clamping in STO thin films.
However, based on the previous phase-field simulations, the misfit strain in the present study (−1.5 to 1.5%) cannot induce an appreciable ferroelectric polarization at room temperature 17 . LAO is treated as a regular dielectric material, without a ferroelectric instability. A quasi-threedimensional simulation with spatial inhomogeneity was performed on a 128Δx × 2Δx × 128Δz grid, where 36 grid points along the z-direction are used to represent thẽ 3.6 nm (9 u.c.) STO-LAO solid solutions and 80 grid points are used for the substrate. Each grid point within the film represents either pure STO or LAO to which the corresponding Landau theory parameters and elastic parameters are assigned [18][19][20] . To construct the graded composition as model inputs, a probability is assigned to each layer of the system to specify the STO fraction in that layer. Then, a random distribution of STO and LAO is generated within each layer. The constructed compositions of the upgraded and downgraded solid solution systems are mapped out in Fig. 1a, b, respectively, where the grids along the out-of-plane direction have been coarsened into nine layers corresponding to the samples with nine unit cells in the experiments, as shown below.
Given the inhomogeneous structures generated above, the equilibrium polarization distribution P(r) was obtained by evolving the time-dependent Ginzburg Landau (TDGL) equation, ∂PðrÞ ∂t ¼ ÀL δF δPðrÞ , where r is the position vector, L is the kinetic coefficient, and F is the total free energy. The elastic, electric, flexoelectric, and gradient contributions are considered in the total free energy F. The mechanical equilibrium equation with film boundary conditions was solved 21 with the eigenstrain originating from electrostriction and the experimentally measured graded strain across the solid solution. Shortcircuit boundary conditions are assumed in the electrostatic equilibrium equation 22 . A more detailed description of the modeling is provided in the Methods section. Starting from a tiny random polarization representing the initial paraelectric states, the two compositional-graded systems are evolved at room temperature to minimize the total free energy using the TDGL equation. After reaching equilibrium, a distribution of electric polarization appears in both scenarios, as plotted in Fig. 1a, b. The upgraded system in Fig. 1a exhibits a generally upward polarization that tends to concentrate at the STO-rich surface. By contrast, the downgraded counterpart shown in Fig. 1b is downward polarized with larger polarization locating at the STO-rich bottom. The direction of the spontaneous polarization agrees with our analysis in terms of strain gradients of the two solid solution systems and flexoelectric-generated symmetry breaking.
To reveal the mechanism of the spontaneous polarization, we averaged the in-plane and out-of-plane strain components, ε 11 and ε 33 , and the out-of-plane polarization, P 3 , within each layer and plotted their profiles along the thickness direction in Fig. 1c for both graded solid solutions. In both scenarios, P 3 reaches highest values (~1.0 μC/cm 2 ) at the STO-rich side and shrinks tõ 0.5 μC/cm 2 at the LAO-rich side. Since the in-plane and out-of-plane strain gradients are nearly constant across the film, their combination should give rise to a uniform polarization across the film assuming a homogeneous flexoelectric property. Therefore, the nonuniformity in the polarization magnitude may arise from local chemical variations. Moreover, it is hypothesized that the strain gradient breaks the symmetry of the systems, producing the overall polarization via flexoelectricity. We can test this conjecture by deliberately turning off the flexoelectric effect in our phase-field simulations while keeping all the other conditions unchanged. As shown by the dashed lines in Fig. 1c, the spontaneous polarization is strongly decreased, suggesting the governing role of flexoelectricity in creating polarity in the mixture of the otherwise nonpolar STO and LAO due to chemically induced strain gradients. A one-dimensional simulation can be found in Fig. S1 in the Supplementary information to confirm that the polarization can be induced purely by the strain gradient.
To fabricate the STO-LAO solid solution in which the mixing ratio is designed on the scale of a monolayer, it is crucial to control the growth mode and the amount of deposited materials. The samples were prepared by pulsed laser deposition with a dual-target system. STO and LAO c Layer-averaged out-of-plane strain ε 33 (solid lines) and in-plane strain ε 11 (dotted lines) in the downgraded (red) and upgraded (blue) STO-LAO solid solutions. d Out-of-plane polarization in the downgraded (red) and upgraded (blue) STO-LAO solid solutions with (solid lines) and without (dashed lines) considering the flexoelectric effect targets were mounted on a computer-controlled exchanging stage to allow the laser pulses to hit the respective target surfaces alternatively to fabricate the STO-LAO solid solution. The growth mode of the materials can be tuned as the layer-by-layer mode by adjusting the growth temperature and oxygen pressure 23 . The deposition rate of an oxide monolayer can be precisely determined by the reflection high-energy electron diffraction (RHEED), and therefore, the composition ratio of xSTO-(1-x) LAO in the monolayer can be controlled by the number of laser pulse once the deposition rates of STO and LAO are determined through the analysis of RHEED data 24 . Two different compositional gradient samples were prepared, one with more LAO in the bottom of the film, the upgraded compositional gradient system, and the other with more STO in the bottom of the film, the downgraded compositional gradient system, as schematically illustrated in Fig. 2a. Additionally, a reference sample with a homogenous mixture of LAO and STO was also prepared. The total thickness of the STO-LAO solid solutions is 3.5 nm, corresponding to nine unit cells. As the composition is varied across each monolayer in the solid solutions, the bulk lattice parameter with different compositional ratios can be calculated according to Vegard's law 25 . The lattice constants of the STO-LAO solid solutions with different composition ratios are given in Fig. 2b. The lattice constants increased from 3.80 Å of the 10% STO-90% LAO layer to 3.89 Å of the 90% STO-10% LAO, corresponding to 1.55% tensile strain and 0.79% compressive strain, respectively. Moreover, in such solid solutions, due to the ultrasmall thickness, the in-plane lattice parameters of these layers are very close to those of the NGO substrates due to the coherent heteroepitaxial growth as evidenced by X-ray diffraction (Supplementary information Fig. S2). The misfit strain of each layer relative to NGO shows a linear relationship with its composition, as shown in Fig. 2c. Therefore, the upgraded and downgraded compositional gradients in these ultrathin films present opposite strain gradients, and an electrical polarization is induced due to the coupling between the strain gradient and the ionic displacement 14,26 .
To confirm the composition distribution in both graded samples and in the reference sample, the atomicresolution X-ray energy dispersive spectroscopy (EDS)mapping technique with a spherical aberration (C s ) corrected scanning transmission electron microscope (STEM) was employed. Figure 3a, b show the high-angle annular dark-field (HAADF) images and the corresponding element distributions of Al, Ti, La, Sr, Ga, Nd, and O at atomic scale for the upgraded and downgraded compositional gradient systems, respectively. For the upgraded compositional gradient system shown in Fig. 3a, the contrasts of Ti and Sr increase from bottom to up while that of Al and La decrease. In contrast, for the downgraded compositional gradient system in Fig. 3b, the contrasts of Ti and Sr decrease from bottom to top, while that of Al and La increase. Additionally, a discernible contrast of Nd appeared in the first La layer in Fig. 3a, indicating that a slight interdiffusion occurred between the Nd and La layers, possibly due to the similarity of the valence and ionic radius of Nd and La. In contrast, there was no interdiffusion between the Nd and Sr layer, as shown by the sharp interface in Fig. 3b. To quantify the compositional gradient, an EDS spectrum was acquired and quantified from each Ti-Al-O layer in the downgraded solid solution. The atomic percentage of Ti from each layer was given in Fig. 3c, and matches the designed composition very well. The plot of the composition versus the number of layers shows perfect linearity, as shown in Fig. 3d, indicating that a uniform compositional gradient was obtained. Moreover, the HADDF and EDS images show that Sr and La are at the A site and Ti and Al are at  Fig. S3 (a), confirming that the desired solid solution can be grown. Figure 4a shows a typical HAADF image of the upgraded compositional gradient film along the NGO [100] zone axis. The interface between the thin film and the NGO substrate is very sharp, as indicated by the overlaid atomic structure. In this perovskite solid solution system, the A sites are occupied by Sr and La and the B sites are occupied by Ti and Al. Both the A sites and B sites are clearly observed in the simultaneously recorded annular bright-field (ABF) image together with the oxygen columns, as shown in Fig. 4b. Taking advantage of the lightelement sensitivity and superior spatial resolution of the ABF imaging technique, the cation displacement and electrical polarization can be obtained from the ABF image 27,28 . Here, to obtain more information, the ABF image was quantified by an improved peaks finding method reported by Yin et al. 29 . The displacement vectors of the B sites were mapped unit-cell by unit-cell, as shown in Fig. 4c. The yellow arrows demonstrate an obvious variation of directions and amplitudes, indicating that the polarizations fluctuate drastically in both lateral and vertical directions. When the vertical components of the displacements were considered, the downward displacements were dominant, revealing an upward net out-ofplane polarization, as schematically indicated by a white arrow in Fig. 4b. It should be noted that obvious in-plane displacements of the B site are present in the bottom layers (LAO-rich), corresponding to a strong in-plane tensile strain. For the downgraded compositional gradient film, Fig. 4d, e show the typical HAADF and ABF images along the NGO [100] zone axis, respectively. Figure 4f is the unit-cell by unit-cell mapping of the displacement vectors of the B sites quantified from Fig. 4e. In this sample, the amplitudes of the displacements are smaller than those in Fig. 4c. However, the upward out-of-plane components of the B-site displacements are dominant, revealing a downward net out-of-plane polarization, as schematically indicated by the white arrow in Fig. 4e. Because the B sites are occupied by both Ti and Al elements with the ratios varying from the surface to the bottom, the amplitudes of the B-site displacements cannot be directly related to the amplitudes of the electrical polarization. However, careful STEM observations convincingly show that the upgraded compositional gradient system has the upward net polarization, whereas the downgraded compositional gradient system has the downward net polarization. Finally, B-site displacements analyses were also performed for the homogeneous solid solution, demonstrating that no clear polarization can be observed, as shown in Fig. S3 (b) and (c).
Since the average vertical displacement of B-site atoms can be acquired from the STEM results, the polarization magnitude can be simple estimated as P = q × s/V, where q is the valence of the B-site atoms, s is the displacement of B-site atoms, and V is the volume of unit cell. We have calculated the Born effective charge of Ti in SrTiO 3 (~7) and Al in LaAlO 3 (~3) as a function of lattice parameter as shown in Supplementary information Fig. S4. The value of q is equal to the weighted average of the Ti and Al valence states. Thus, it varies across the whole solid solution. The volume of the unit cell (V) is equal to c × a 2 f . a f is the inplane lattice constant, 3.863Å, because the solid solution system is coherent to the NGO substrate, and the variations of c axis for the two solid solution systems are shown in Fig. 4g, h for the upgraded and downgraded compositional gradient solid solution, respectively, as obtained from the ABF images shown in Fig. 4c, f. The average vertical displacement of the B-site atoms in both compositional gradient solid solutions were also quantified from the corresponding ABF images and are plotted in Fig. 4i, j. Therefore, the polarization for the upgraded compositional gradient system is~16.5 μC/cm 2 and that for the downgraded compositional gradient system is 6.8 μC/cm 2 . In BaTiO 3 -based perovskite structure, the increase in the displacement of the Ti 4+ cation leads to the enhancement of the polarization 30 . In our case, the upgraded compositional gradient system shows larger displacements of the B-site cation compared with those of the downgraded system, and therefore, the upgraded compositional gradient system possesses the larger magnitude of polarization than the downgraded system. Such an asymmetric polarization may be attributed to the different strain relaxation behaviors of the two samples, as revealed in Fig. 4g, h.
Because the STO-LAO solid solutions are only nine unit cells thick, the in-plane lattice of the whole film is coherent to that of the NdGaO 3 substrate, suggesting a single and uniform in-plane lattice parameter for the solid solution 31 . Since the lattice parameter of the bulk varies with the corresponding compositional ratio, the strain for each monolayer is different and a large-strain gradient can be generated. Therefore, the theoretical magnitude of the spontaneous polarization based on the flexoelectric effect can be calculated. The strain gradient of the system is 0.00675 nm −1 , and the magnitude of the spontaneous polarization based on the flexoelectric effect is~4 μC/cm 2 based on the flexoelectric coefficient of 6.1*10 −9 μC/m for STO 32 . Different directions of the strain gradients (upgraded and downgraded) can lead to different net outof-plane polarizations. When the STO-LAO solid solutions are examined by piezoresponse force microscopy (PFM), it is found that both types of compositional gradient films possess electric-field switchable PFM responses, as shown in Fig. 5a. The observable differences between the results of upgraded and downgraded samples are the asymmetric voltages required to switch the PFM phases (Supplementary information Fig. S5). However, the hysteresis behaviors of the PFM responses can be attributed to various mechanisms, including the switching of spontaneous polarizations and surface-charging effects. In Fig. 5a, the area poled with a positive voltage has a relatively weak PFM amplitude, which is a conventional The displacement of B-site atom of (i) upgraded and (j) downgraded compositional gradient STO-LAO solid solutions. In these figures, the layers on the right-hand side correspond to the interface indication that the phase inversion is due to the local surface charging. By further examining the sample with contact Kelvin probe force microscopy (cKPFM) 33,34 , where the PFM signal is read at various dc bias V read values after the sample is written by different V write , we found the PFM hysteresis is due to the charge injection. As presented in Fig. 5b, c, the linear PFM responses with V read indicate nonferroelectric behaviors, and the nearly constant slopes are due to the similar mechanism of the electrostatic effects in the traditional surface potential measurement. The shift of each line is the source of PFM hysteresis, and the V read for minimizing the PFM signal under different V write corresponds to the level of the surface charge injection. In the upgraded STO-LAO solid solution, the zero cKPFM signals correspond to the sample bias V read ranging from −1.3 V to 0.8 V under V write from −20 V to + 20 V, while in the downgraded STO-LAO solid solution, the zero cKPFM signals correspond to V read ranging from −2.0 V to −0.2 V under the same conditions. These results reveal that the downgraded sample shows easier positive surface charge injection than the upgraded sample, satisfying the boundary conditions based on predicted spontaneous polarization directions of the phase-filed simulations and in agreement with the atomically resolved STEM observations. As shown in Supplementary information Fig. S6, for the reference sample, the homogeneous solid solution has easier positive surface charge injection than the upgraded solid solution, but has harder positive surface charge injection than the downgraded solid solution, which supports the presence of the additional surface charge effect induced by the polarization arising from compositional gradients.
These new inhomogeneities in materials can be viewed as a new material genome. To further demonstrate the implications of our results, we design new heterostructures with the combination of the up-and downgraded solid solutions. First, an upgraded structure was grown on the NGO substrate. Then, a downgraded structure was fabricated on top of the upgraded structure to form a head-to-head polarization configuration, as shown in the left panel of Fig. 6a. On the other hand, by first growing a downgraded structure and then an upgraded structure on top, we obtain a tail-to-tail polarization configuration, as shown in the right panel of Fig. 6a. Local conduction can be expected, since the polarization accumulates charge carriers due to the screening effect 35 . These two structures (head-to-head and tail-to-tail) were measured and revealed show local conduction with the resistance of the head-to-head configuration lower than that of the tail-to-tail configuration, as shown in Fig. 6b. The resistance increased with decreasing temperature, indicating insulating behavior. It should be noted that the effect of oxygen vacancy should not vary much since each layer always maintains charge neutrality. We also analyzed the Ti-L 2,3 ELNES from the STO-rich layers and LAO-rich layers to search for the oxygen vacancies as shown in Supplementary information Fig. S7. Both spectra show the feature of Ti 4+ without obvious differences, indicating that there were no detectable oxygen vacancies in both thin films. Furthermore, no conduction was detected in the heterostructure with the homogenous composition. The results demonstrate that oxygen vacancies should not be the dominating factor for the conductance in the samples. We also performed temperature-dependent Hall measurements to characterize the charge-carrier type and carrier concentration of the structures. The carrier type is n-type in both cases. The calculated mobility is 0.5 cm 2 V −1 s −1 for the head-to-head structure and 1.3 cm 2 V −1 s −1 for the tail-to-tail structure at 300 K, respectively, as shown in Fig. 6c. Furthermore, the carrier concentration given by the inverse relation of the measured Hall coefficient can be extracted as shown in Fig. 6d. The carrier concentration, n, is~10 14 cm 2 . In addition, the carrier concentration shows almost no temperature dependence.

Conclusions
The compositional gradient of the STO-LAO solid solution is successfully fabricated by pulsed laser deposition with the assistance of RHEED. Atomically resolved EDS maps show the distribution of each element, confirming the existence of compositional gradients. Furthermore, the results of the average displacements of B-site atoms obtained from C s -corrected ABF-STEM images suggest that the upgraded compositional gradient system shows an upward polarization and the downgraded compositional gradient system displays a downward polarization, which is supported by the phase-field simulations and further confirmed by PFM observations. The polarization direction is controlled by the strain gradient in the solid solutions. Furthermore, these upgraded and downgraded structures were combined to form the different head-to-head and tail-to-tail polarization configurations. Although the charge carriers in both of these structures are n-type, the mobility in these two structures is different, giving rise to the difference in the conductivity, which is absent in the homogenous solid solution. These results suggest that electrical polarization can be induced through the compositional gradient, and can be used to affect the local conduction behavior, providing a new pathway for acquiring a material genome based on material inhomogeneity.

Experimental section
Film growth STO-LAO solid solutions were grown on an NGO substrate by pulsed laser deposition. NGO substrates were preannealed at 1000°C for 2 h to eliminate the contamination and obtain a flat surface. The deposition process was carried out at the substrate temperature of 750°C with the oxygen pressure maintained at 2 × 10 −5 Torr. A KrF excimer laser (λ = 248 nm) at a repetition rate of 3 Hz and a fluence of ∼1.75 J cm −2 was used to ablate the stoichiometric LAO and STO targets. RHEED was adopted to monitor the growth condition and growth rate in real time. At these conditions, LAO and STO layers were found to follow the layer-by-layer growth mode. Estimating from the RHEED oscillation data, the deposition rate was determined as 0.028 u.c./pulse for LAO and 0.042 u.c./pulse for STO. Then, a dual targets process was adopted to fabricate the compositional gradient samples by controlling the laser pulses of each target. For a detailed description of the dual targets process, please refer refs. 36,37 .

STEM observations
The thin film sample was first coated with a protective layer of amorphous carbon (≈50 nm). A cross-sectional TEM specimen was prepared with a dual-beam focused ion beam (FIB) system (NB5000, Hitachi, Japan) using Gaion accelerating voltage ranging from 2 to 40 kV, and then ion-milling (Gatan 691, Gatan, USA) operated from 1.5 to 0.5 kV was carried out in order to remove the residual amorphous layer.
STEM observations were performed using a Cscorrected STEM (JEM-2100F, JEOL, Co., Tokyo, Japan) operated at 200 kV and equipped with a spherical aberration corrector (CEOS Gmbh, Heidelberg, Germany), providing a minimum probe of~1 Å in diameter. During STEM imaging, a probe convergence angle of 25 mrad was used. The detection angles of 92-228 mrad and 8-25 mrad were used for HAADF and ABF imaging, respectively.
Atomic-resolution EDS analyses were performed using an aberration-corrected scanning transmission electron microscope (JEM-ARM200CF, JEOL Co. Ltd) operated at 200 keV. STEM-EDS maps were acquired by scanning the beam over different regions, and using the NSS3 software developed by Thermo Fisher Scientific Inc. The STEM-EDS system was equipped with a silicon drift detector (SDD) and the solid angle for the whole collection system was~1.7 sr. The probe size was 1.2 Å with a probe current of~60 pA. The spectral peaks used for EDS mapping were K α = 4.510 keV for Ti, K α = 14.163 keV for Sr, K α = 1.486 keV for Al, L α = 4.650 keV for La, K α = 9.250 keV for Ga, L α = 5.229 keV for Nd, and K α = 0.532 keV for O. The individual element maps were also summed to produce a combined map. To obtain sufficiently strong signals for individual spectra, the total acquisition time for each map was~100 min with a dwell time of 20 ms. The measurements were repeated to ensure reproducibility and were found to be little affected by the variation of the experimental conditions reported here.

Phase-field modeling
The STO-LAO solid solution system was modeled using the phase-field method for ferroelectrics based on the Ginzburg-Landau-Devonshire theory 38 . The free energy density including the bulk, electric, elastic, flexoelectric, and gradient energy densities of STO and LAO, assuming a linear modification by the local concentration of LAO, x(r), can be written as The detailed expressions of each term in the total free energy density can be found in our previous works [17][18][19] , and the same parameters for STO were adopted here. The flexoelectric coefficients of STO calculated by firstprinciples theory were used 39 . The LAO was treated as a regular dielectric, where the bulk energy is described by a second-order term of polarization f LAO bulk ¼ a LAO 1 P 2 , where a LAO 1 is the dielectric stiffness of LAO. For simplicity, the parameters in the other energy terms of LAO are assumed to be identical to those of STO.
The heterogeneous structure of STO-LAO was generated by randomly assigning x(r) to each grid point of the simulated system with x = 1 representing LAO and x = 0 representing STO. The concentration gradient along the film thickness was controlled by specifying the statistical possibility of x = 0 at each layer. After generating the heterogeneous systems, we started from a distribution of random-noise polarization to mimic the initial paraelectric state. Then, the polarization evolves at room temperature as governed by the TDGL equation to minimize the total free energy until an equilibrium is reached.

Piezoresponse force microscopy
The piezoresponse force microscope (PFM) images and contact Kevin probe force microscopy (cKPFM) results were obtained using Bruker Multimode using nanosensor probes (spring constant~7 N/m), the signal was read with AC 1 V. The phase-voltage hysteresis loops of the compositional gradient samples without the bottom electrode in the Supplementary information were measured by superimposing a 0.2 Hz triangular square-stepping wave with a bias window up to 70 V with AC 1 V to read the signal.

Electrical transport measurement
The electrical contacts were fabricated in the van der Pauw configuration by photolithography. The contacts were made by 5 nm Ti and 25 nm Au with the contact dimensions of 0.5 mm × 0.5 mm. The temperaturedependent resistance and Hall measurements were performed using a Physical Property Measurement System (PPMS) from Quantum Design.