Ferroelectrics at the Nanoscale: A First Principle Approach

The field of ferroelectric materials is driven by its possible use in various micro-electronic devices that take advantage of their multifunctional properties. The existence of a switchable spontaneous polarization (see Figure 1) is at the basis of the design of non volatile ferroelectric random access memories (FERAMs), where one bit of information can be stored by assigning one value of the Boolean algebra ("1" or "0") to each of the polarization states. Also, the high dielectric permitivity of ferroelectrics makes them possible candidates to replace silica as the gate dielectric in metal-oxide-semiconductor field effect transistors (MOSFETs). Their piezoelectric behavior enables them to convert mechanical energy in electrical energy and vice versa. Their pyroelectric properties are the basis for highly sensitive infrared room temperature detectors.


Introduction
The field of ferroelectric materials is driven by its possible use in various micro-electronic devices that take advantage of their multifunctional properties.The existence of a switchable spontaneous polarization (see Figure 1) is at the basis of the design of non volatile ferroelectric random access memories (FERAMs), where one bit of information can be stored by assigning one value of the Boolean algebra ("1" or "0") to each of the polarization states.Also, the high dielectric permitivity of ferroelectrics makes them possible candidates to replace silica as the gate dielectric in metal-oxide-semiconductor field effect transistors (MOSFETs).Their piezoelectric behavior enables them to convert mechanical energy in electrical energy and vice versa.Their pyroelectric properties are the basis for highly sensitive infrared room temperature detectors.
The current miniaturization of microelectronic technology, imposed by the semiconductor industry, raises the question of possible size effects on the properties of the components.Except for the case of the gate dielectric problem in MOSFET's, the thickness of the films used in contemporary applications is still far away from the thickness range where size effects become a concern; therefore the question at the moment is merely academic.Nevertheless, it is very possible that the fundamental limits of materials might be reached in the future.
Despite the efforts and advancement in the field, both theoretically and experimentally (see a recent review in Ref. [1]), many questions still remain open.The main reason for the poor understanding of some of the size effects on ferroelectricity is the vast amount of different effects that compete and might modify the delicate balance between long range dipole-dipole electrostatic interactions and the short range forces, whose subtle equilibrium is known to be at the origin of the ferroelectric instability.The arrows indicate the atomic displacements, exaggerated for clarity.The origin has been kept at the Ba site.Two structures, with polarization along [001], are shown.Application of a sufficiently large electric field causes the system to switch between the two states, reversing the polarization.in principle, we can assign one value of the Boolean algebra ("1" or "0") to each of the polarization states.
The next section will give an overview of the main concepts involved in the modern theory of polarization.A definition of a local polarization will be given in terms of the centers of the wannier functions associated with the band structure of the system.Next, some basic electrostatic notions related with ferroelectric films will be given, in particular the concept of depolarization field and screening by metal contacts.Finally, the results of this research work: applying the layer polarization concept, where we show the hidden structure of the polarization at the nanoscale.

Modern theory of polarization
Classically, the macroscopic polarization in dielectric media is defined to be an intensive quantity that quantifies the electric dipole moment per unit volume [2][3][4].Definitions along these lines work well for finite systems but have important conceptual problems when applied to periodic crystalline systems because there is no unique choice of cell boundaries [5].
In the early 1990s a new viewpoint emerged and lead to the development of a microscopic theory [6][7][8].It starts by recognizing that the bulk macroscopic polarization cannot be determined, not even in principle, from a knowledge of the periodic charge distribution of the polarized crystalline material.This establishes a fundamental difference between finite systems (e.g., molecules, clusters, etc.) and infinite periodic ones.For the first case, the dipole moment can be easily expressed in terms of the charge distribution.While for periodic systems one focuses on differences in polarization between two states of the crystal that can be connected by an adiabatic switching process [6].The polarization difference is then equal to the integrated transient macroscopic current that flows through the sample during the switching process.
Therefore, the macroscopic polarization of an extended system is, according to the modern viewpoint, a dynamical property of the current in the adiabatic limit.The charge density is a property connected with the square modulus of the wavefunction, while the current also has a dependence on the phase.Indeed, it turns out that in the modern theory of polarization [7][8][9], the polarization difference is related to the Berry phase [10,11], defined over the manifold of Bloch orbitals.The theory not only defines what polarization really is, but also proposes a powerful algorithm for computing macroscopic polarizations from first principles.
The modern theory can be equivalently reformulated using localized Wannier functions [12][13][14][15] instead of extended Bloch orbitals.The electronic contribution to the macroscopic polarization P is then expressed in terms of the dipole of the Wannier charge distribution associated with one unit cell.In this way, P is reformulated as a property of a localized charge distribution which goes back to the classic definition of polarization.However, one has to bear in mind that the phases of the Bloch orbitals are essential for building the Wannier functions.They are needed to specify how the periodic charge distribution should be decomposed into localized ones.The knowledge of the periodic charge distribution of the polarized dielectric is not enough to determine the wannier functions.

The main concepts
Here I will review the central topics uncovered in the early 1990s, often known as the modern theory of polarization.Understanding this background is necessary in order follow this chapter.The general idea is to consider the change in polarization of a crystal as it undergoes a slow change and relate this to the current that flows during this adiabatic evolution of the system.These ideas will lead to an expression for the polarization that does not take the form of an expectation value of an operator, but takes the form of a Berry phase, which is a geometrical phase property of a closed manifold (the Brillouin zone) on which a set of vectors, the occupied Bloch states, are defined.
As a classical example of a Berry phase, lett's think of the paralel transport of a vector along a loop on a shpere ( for instance a compass needle carried in a car traveling on the surface of the Earth).After completing a closed path, or loop, the vector will be back at the original starting point but it will be rotated with respect to the direction it was pointing at when it started the trip.The reason for this rotation is purely geometrical topological and intrinsicaly connected to the curvature of the sphere and would not exist if the the vector would be parallel transported along a flat manyfold, like a plane, or cylinder.This rotation angle is related to the integral of the curvature on the surface bounded by the loop.Now lets see the quantum counterpart and its connection to the definition of Polarization in a solid.Let's assume that the crystal Hamiltonian H λ depends smoothly on parameter λ and has Bloch eigenvectors obeying H λ |ψ λ,nk = E λ,nk |ψ λ,nk and that λ changes slowly with time, so that it is correct to use the adiabatic approximation.
Since the spatially averaged current density is just j = dP dt = dP dλ dλ dt we can write the change in polarization during some time interval as △P = j(t)dt (1) and is phrased in terms of the current density that is physically flowing through the crystal as the systems traverses some adiabatic path.It can also be written in terms of the parameter λ as: Then, from Ref. [7] and the rate of change of polarization with λ is a property of the occupied bands only.This expression can be integrated with respect to λ to obtain It can be verified by taking the λ derivative of both sides of Eq.4 and comparing with Eq.3.
The result is independent of the particular path of λ(t) in time, and depends only on the final value of λ as long the change is adiabatically slow.We can associate the physical polarization of state λ with P(λ) and drop the λ label.Eq. 4 can be written then as and this is the electronic contribution to the polarization.To obtain the total polarization it must be added the nuclear (or ionic) contribution where the sum is over atoms s having core charge Z s and spatial position r s in the unit cell of volume Ω.
The equation 5 is the main result of the modern theory of polarization It states that the electronic contribution to the polarization of crystalline insulator may be expressed as a Brillouin zone integral of an operator i∇ k .However, it is not a common quantum mechanical operator, the result of its action on the wavefunctions (i∇ k |u nk ) depends on the relative phases of the Bloch functions at different k.
To understand the nature of the integrand in Eq.5 I want to recall the fundamental paper by Zak [16] in which he postulates the existence of a Berry phase associated with each one of the bands of a one dimensional crystalline solid.He shows how the variation of k over the entire Brillouin zone produces the appearance of a Berry phase.Moreover, he associates this phase, that is a characteristic feature of the whole band, with the band center operator [17].Its eigenvalues turn out to give the average positions of an electron in different bands which coincide with the symmetry centers of the space group of the solid [17] These previous ideas set the stage for the crucial meaning of the expression in Eq.5.
In three dimensions, the Brillouin zone can be regarded as a closed 3-torus obtained by the identifying points u nk = u n,k+G j where G j are the three primitive reciprocal lattice vectors.Then, the electronic polarization contribution of each band can be written as where R j is the real space primitive translational corresponding to G j , and the Berry phase for band n in direction j is More details and discussion, including the equivalent of Eq.7 for the case of connected multiple bands, may be found in Refs.[7] and [8].
Reformulation in terms of Wannier Functions.
We can rewrite Eq.5 in terms of the Wannier functions (WF's), which brings an alternative and more intuitive way of thinking about it.The WF's and Bloch functions can be regarded as two different orthonormal representations of the same occupied Hilbert space.The WF's are localized functions w nR (r) that are labeled by band n and unit cell R and constructed by carrying out a unitary transformation of the Bloch states ψ nk .They are constructed via a Fourier transform of the form where the Bloch states are normalized in one unit cell.There is some freedom in the choice of the these WF's, as the transformation is not unique.In particular, a set of Bloch functions results in WF's w nR ′ which are different from the w nR .Usually, this "gauge" is set by some criterion that keeps the WF's well localized in real space, such as the minimum quadratic spread criterion introduced by Marzari and Vanderbilt [19].However, it is expected that any physical quantity, such as the electronic polarization arising from band n, should be invariant with respect to the phase change β n (k).
Once we have obtained the WF's, we can locate the "wannier centers" r nR = w nR |r|w nR .It can be shown that where φ nj is given by Eq.7.In simple words, the location of the n'th Wannier center in the unit cell is just given by the three Berry phases φ nj of band n in the primitive lattice vector directions R j .The key result, is that the polarization is just related to the Wannier centers by The Berry Phase theory can then be regarded as providing a mapping of the distributed quantum mechanical electronic charge density onto a lattice of negative point charges −e (see Fig. 2).

Layer polarization definition
One issue that has received much attention theoretically is how to quantify the concept of local polarization.This can be very useful in understanding the enhancement or suppression of spontaneous polarization; or even to find new spatial patterns of polarization that would remain hidden otherwise, as is the case that I will show in a following section.It can also be essential for characterizing and understanding interface contributions to such properties.
But first I will give a review of the concept introduced in 2006 by Wu.et al. [20] which is the one used here.As explained in the previous section, the modern theory of polarization establishes that the polarization of a crystal is expressed in the Wannier representation as the contribution of the ionic and electronic charge: where s and m run over ion cores (of charge Q s located at R s ) and Wannier centers (of charge -2e located at r m ), respectively, in the unit cell of volume Ω.
They defined the Layer polarization (LP) along z for superlattices built from II-VI ABO 3 perovskites such as BaTiO 3 , SrTiO 3 , and PbTiO 3 .For these cases, they were able to decompose the system into neutral layers (that is AO and BO 2 subunits) and define a layer polarization in which the sums are restricted to entities belonging to layer j. S is the basal cell area and we are now focusing only on z components.The LP p j thus defined has units of dipole moment per unit area.The total polarization, with units of dipole moment per volume, is exactly related to the sum of LP's via where c = Ω S is the supercell lattice constant along z.They propose two conditions to be able to associate a physical meaning to the LP definition: (i) resolve the arbitrariness associated with the positions of the Wannier centers, and (ii) associate without ambiguity the right number of Wannier centers to each layer.As it will be shown in the next section, we applied this concept to thin ferroelectric films sandwiched between metal contacts.These kind of systems, in which metals and oxides coexist, bring new challenges.Condition (i) was satisfied using maximally localized wannier functions plus a disentanglement procedure (Ref.[21]) and (ii) was satisfied, even close to the metal oxide interfaces.

Depolarization field
When dealing with ferroelectric thin films, a finite polarization normal to the surface will give rise to a depolarizing field and a huge amount of electrostatic energy, enough to be able to suppress the ferroelectric instability.In order to preserve ferroelectricity and minimize the total energy, the depolarization field must be screened, either by free charges coming from metallic electrodes or by breaking into a domain structure.I will go into more details about this, but first a basic electrostatic background will be given.

Influence of the electrical boundary conditions
In the case of a free standing slab of a ferroelectric material, a uniform polarization P with an out of plane direction will appear and originate a surface polarization charge density where n is a unit vector normal to the surface pointing outward.The charge density is proportional to the magnitude of the normal component of the polarization inside the material [22] and is positive one one side of the slab and negative on the other one.In order the determine the electric fields generated by the surface charge density, we need to know the electrical boundary conditions of the problem and solve the Poisson equation.For the general case, let's suppose that an external electric field E ext , perpendicular to the surface, is applied in the vacuum region.There, the electric displacement D = E + 4πP equals the applied electrical field D = E ext as the polarization is zero.The normal component of D is conserved across the ferroelectric/vacuum interface, Now we can see two extreme cases.One of them is that of zero external electric field, E ext = 0, in which case the displacement vector is null while the internal field inside the slab is E = −4πP (see Fig. 3a).The absence of an external electric field and the presence of the surface polarization charges produce an internal field that points in the direction opposite to that of the polarization.Therefore, it tends to restore the paraelectric configuration and that is the reason of its name, Depolarization field.The coupling between the polarization and this internal field is so big, that any atomic relaxation under these circumstances will end up with the atoms back in the centrosymmetric non polar configuration.Another case, a vanishing internal electric field E = 0 corresponds to a case where no field opposes the polarization and a spontaneous polarization might arise (Fig. 3b).In this case, it is the continuity of the normal component electric displacement vector that establishes the value of it, These two cases illustrate how a polarization perpendicular to the surfaces can exist or not depending on the electrical boundary conditions.For simulating a real system, the mechanical boundary conditions should be included via atomic relaxations.These cases have been extensively analyzed [23][24][25][26] and have opened the ground for more realistic systems..

Metal contacts
The next level of approximation to a real case,is the case of having the ferroelectric thin film sandwiched between metal contacts ( short circuit boundary conditions).In the case of a perfect metal, the surface charges σ pol will be exactly compensated at the interfaces by the compensation charges σ com and we will have a null depolarization field (see Fig. 4b).In this case, both σ ′ s will lie in a sheet right at the interface between the metal and the ferroelectric electrode.However, in the case of real electrodes, the compensation of the polarization charges will be incomplete.In this case,the compensation charge σ com will be spread over a finite distance inside the electrode and will create an interface dipole density at each of the ferroelectric/metal interfaces (see Fig. 4c).Due to the short circuit boundary conditions, an amount of charge must flow from one interface to the other, creating a residual depolarization field inside the thin film.There are various models for expressing this field [24,[27][28][29], shown below is one described in Ref. [30]: where d is the film length, and λ is the effective screening length, proportional to the spatial spread of the interface dipole.

Electrostatic energy of the thin film
In the presence of a residual depolarizing field the energy of the thin film U can be approximated by where U is the internal energy under zero field.If we assume that the interface effects are hidden in the screening length parameter, then they only appear through the depolarization field.This is a rough approximation but we will introduce later a model that will relax this condition, but let's keep it for now.Following this, then where m stands for the number of unit cells of the film.U FE can be calculated from bulk DFT calculations.Its dependence on P 0 has typically the shape of a double well.
The electrostatic energy E elec (P 0 ) can be approximated by [1, 31]   E elec The depolarization field can be calculated from first principles by taking the electrostatic average [32,33] or can be expressed using Eq.18.In this case, Eq.19 can be rewritten as The electrostatic energy, the second term, is positive, meaning that the effect of the depolarization field is to suppress the ferroelectric instability by rescaling the quadratic term of the double well.This simple model has been implemented [34,35] in order to extend the  Values of the quantities at the bulk level are represented by dashed lines for the unstrained configuration, and with dotted lines for a geometry under the strain imposed by the substrate.The evolution of the system, from a monodomain configuration at large thickness, (where the depolarization field E d is small), to a 180• stripe domain structure in order to minimize the energy associated with E d , is represented in the inset.From [37] first principle model Hamiltonian to thin films.It also has been used with the result of first principle calculations in Ref. [30].In this work, we will use it as a starting point for a toy model that would mimic the result of our first principle calculations.

Effects of the depolarization field
There are two ways to minimize the energy due to the depolarization field , i) a reduction of the polarization, ii) a reduction of the depolarization field itself.This can be seen by inspecting the first term in Eq.21.
As an example of case i) it is shown in Fig. 5 the experimental results for the tetragonality in function of the film thickness for monodomain PbTiO 3 epitaxially grown on Nb-doped SrTiO 3 [35].The polarization is strongly coupled to strain in perovskites oxides [36], therefore a reduction of the polarization must be accompanied by a reduction of the tetragonality of the system.The incomplete screening of the depolarization field is the driving force for the reduction of the polarization.The structure remains in a monodomain at all thickness.
Case ii) might arise when there are no compensation charges provided by the electrodes or when the compensation charges do not provide an efficient enough screening of the polarization charges.In this case the system might break up into 180• stripe domains to reduce the magnitude of the surface dipole density (see Fig. 6).The reason why some systems remain in a monodomain configuration while other similar heterostructure break up into domains remains an open question that requires further clarification.
In a seminal theoretical paper, Junquera et al. [30] has suggested that the appearance of ferroelectricity in thin films of BaTiO 3 , with SrRuO 3 contacts, should be limited to a thicknesses of more than six unit cells.Below that limit, ferroelectricity would effectively be canceled by the strength of the depolarization field (DF) inside the oxide produced by unscreened charges at the interfaces, which, in that work, were calculated for the first time in a first principles framework.More recently this picture has been broadened, and a variety of studies have shown both experimentally and with more realistic ab initio simulations that thin oxide films can indeed maintain some ferroelectricity at thicknesses even smaller than six unit cells.[37][38][39][40] 5.The hidden nature of ferroelectricity at the nanoscale The size limit of ferroelectricity in ultrathin films can be understood using a simple electrostatics argument: as the thickness of the oxide film is reduced, the intensity of the DF is made stronger, [30] thus increasing the electrostatic energy of the system.In a ferroelectric crystal, this energy contribution can be minimized in two ways: either by breaking the polarization pattern in 180o stripe domains [41][42][43][44] to reduce the magnitude of the surface dipole density, or by a reduction of the ionic polarization while the system remains in a monodomain state.[35,45] The final polarization pattern depends on the individual material and/or structure, and understanding which system modification will occur under what conditions is still a matter of debate.
At small dimensions, the spatial confinement and the interface details become extremely relevant since they determine the spatial localization of the electrons and thus influence directly the ferroelectric characteristics and, more in general, all the electronic properties of the system.[46] Therefore, to clarify the behavior of ferroelectric thin films in the nanometer regime it is crucial to obtain a precise description not only of the geometry and electronic properties of the film but first and foremost of the polarization profile inside the oxide.
In order to obtain a complete description of ferroelectricity at the nanoscale we have exploited the notion of layer polarization (LP).This idea, based on modern theory of polarization and the concept of maximally localized Wannier functions (see Sec.2), was recently introduced to describe the layer-by-layer modulation of polarization in ferroelectric superlattices (Sec.3) .
We have applied this method to a model metal-ferroelectric system (BaTiO 3 /Pt) (see Fig. 7) and obtained a detailed spatial profile of the polarization in the oxide region in the direction of growth that accounts not only for the ionic displacements but also for the spatial rearrangement of the electron density.It is important to note that using more standard methods to evaluate macroscopic polarization, this information would have remained hidden by the averaging procedure, or would have produced misleading results.In fact, by considering exclusively the ionic displacements (rumpling) as a measure of polarization, we do not account for the full electronic contributions to the local polarization and, depending on the system size, we might obtain results that would not describe accurately the physical behavior of the thin film.
Here, we show that the ferroelectric structures associated with small length-scales are more complex than previously thought, mainly due to the redistribution of the electrons under the constraints imposed by the interfaces.In particular, the oxide develops a ferrielectric [47] pattern of polarization that exhibit antiferroelectric properties along the growth axis and with a net total polarization (see Fig. 8).
This unbalanced dipole structure is particularly evident in Fig. 10A, where we show how positive rumpling on all the atomic layers of a thin BaTiO 3 film correspond actually to layer dipoles of opposite signs.This is a consequence of the interplay between the orientation and magnitude of the dipoles with the DF, their mutual interaction and the nature of the interfaces.A simple analytic model of a ferroelectric thin film, where these effects are explicitly taken into account, can capture all these features and it will be discussed in detail at the end.
We have simulated thin films of BaTiO 3 between Pt metal contacts in a (001) stack, where the oxide is terminated with a BaO plane at both interfaces.In this configuration, the metal atoms are directly bonded to the O atoms of the oxide plane.[30,48] The different supercell constructed in this way can be labeled as Pt/(BaO-TiO2)m-BaO/Pt with m=1,2,4,6.Nine atomic planes of Pt have been found to be sufficient to simulate the contacts under short-circuit boundary conditions.The in plane lattice constant of the supercell is set equal to that of BaTiO 3 at 0 K (3.991) [38], and kept it fixed in all calculations while all other geometrical parameters (atomic positions and intra-layer distances) were fully relaxed.All simulations have been performed using Density Functional Theory (DFT) within the Local Density Approximation, with ultrasoft pseudopotentials and a plane waves basis set.We used the expression of Eq.12 to calculate the LP for each plane of the thin film.A is the in plane area of the cell, and we are only interested on the z components (parallel to the growth direction).It is important to stress that in the definition of the LP we have explicitly included the full ionic and electronic information necessary to define unambiguously the local value of the polarization.The essential step for the evaluation of the layer polarizations p j is the determination of the Wannier centers z m .We used the maximally localized Wannier functions (MLWFs) algorithm originally proposed by Marzari and Vanderbilt [19] as implemented in the WanT code 1 .Since the valence bands of the oxide and the bands of the metal are mixed in the full supercell calculation, a disentanglement procedure (Ref.[21]) was applied before starting the localization algorithm [19,49].
The set of MLWFs so obtained cluster around each atomic layer in sets containing the correct number of Wannier function needed to keep the layer charge neutral (Fig. 9).This is a crucial condition for the validity of the definition of LP's in Eq.12 and might not be attainable in systems with more covalent nature.[20] In fact, the meaningful use of this concept is based on the condition of being able to decompose the system into neutral layers (TiO and BaO2 units in the case of perovskite BaTiO 3 ) and resolve the arbitrariness associated with the positions of the Wannier centers.This allows us to associate each sets of Wannier centers to the corresponding crystal plane without ambiguity.By doing so we go back to the classic idea of defining dipoles in spatial neutral units, along the lines of the Clausius-Mossotti limit.vMoreover, we can associate each LP value p j with its correspondent (dimensionally correct) dipole p j A, which gives precise information about the orientation and magnitude of the atomic layer dipoles inside the crystal.
The LP profile calculated for the metal-oxide systems provides a detailed local description of the polarization, which is superior and even contradictory to the information obtained from empirical evaluations based solely on ionic displacements.As an illustration, in Fig. 10A we display the rumpling profile (cation-oxygen normal distance within each atomic layer.)for a thin film comprised of six BaTiO 3 unit cells (m6).From the ionic displacement profile one would be tempted to conclude that the film is in a ferroelectric domain configuration, since the rumpling values (and hence the classical local dipoles) associated to each plane are positive.Instead, if the correspondent LP profile is analyzed, a completely different picture emerges, where the main feature is an alternation of positive and negative values of the LP in the individual atomic layers of the oxide film.This is clearly shown in Fig. 10B, where we plot the LPs values computed for m6.Here the central BaO planes (the BaO planes at the interface behave differently and will be discussed later) have positive dipoles while the TiO2 planes have negative ones.The first are aligned with the direction of the DF, while latter oppose it.
A schematic illustration of the associated dipoles is drawn above the plot.This result implies that what could have been interpreted as a ferroelectric domain using standard geometrical information, in fact displays a much more complex spatial pattern of polarization.We believe that this effect is comparable to the onset of the formation of two-dimensional domain islands commonly observed both theoretically and experimentally in ferroelectric films.[37] In that case, the domains are composed of adjacent 180• domains in order to minimize the electrostatic energy associated with the DF.Here, our results show that the system can reduce its electrostatic energy in an alternative way: the energetic interplay between the dipoles, the DF and the interface bonding drives the individual atomic layer dipoles to arrange themselves in an uncompensated dipole pattern along the growth direction (see inset of Fig. 10A), or, in other words, the system undergoes a ferroelectric-to-ferrielectric phase transition.Similarly to the formation of two-dimensional islands, we also predict the existence of a critical thickness (CT) above which the system exists in a true ferroelectric domain with all the dipoles pointing in the same direction (see Fig8B, where the LP's values for bulk BaTiO 3 are shown with triangles).
As the thickness is reduced, the magnitude of all the dipoles is diminished by the increase of the DF that opposes them.Due to the different physical properties of consecutive layers (BaO and TiO2), their associated dipoles reduce their magnitude at different rates.At the CT, the BaO layers will have null dipoles and eventually they will flip layer polarization direction, aligning with the DF in order to minimize the associated internal energy and creating a ferrielectric dipole pattern (FDP).
The formation and characteristics of the FDP are determined by the bonding properties of the interfaces that pin the LP values of the outer layers of the oxide film.At small film thicknesses the interface effects become dominant and induce the overall spatial variation of the LP's that determines the orientation and magnitude of the layer dipoles along the structure.We can see how the FDP observed in the m6 geometry is increasingly influenced by the locality of the interfaces as the thickness is reduced to m4, m2 and m1 (see Fig. 10C,D,E respectively).The middle crystal planes become more and more influenced by the interface environment, while the FDP becomes irregular, and finally disappear when the film is comprised by just 3 layers (i.e. one unit cell, m1), leaving instead a centro-symmetric paraelectric structure (Fig. 10E).It is important to note that exclusively electronic effects drive this series of regime transitions.
In fact, the ionic polarization, directly proportional to the ionic rumpling, remains almost constant for all the thicknesses considered, and it is only the electronic polarization that changes its values in the different systems.
These results show that the ferroelectric response of a thin film is indeed critically influenced by the interface properties of the system and that the analysis of the local variation of layer polarization captures completely the physical characteristics of the ferroelectric.As the size of the film is increased, the DF decreases and the dipoles that are oriented along the DF direction get smaller, while the interface effects lose relevance.At the critical size, the dipoles flip direction and a ferroelectric domain structure is established again.Eventually, the dipoles will relax to the bulk structure, showed with triangles in Fig. 10 Model We further illustrate these physical principles by developing a simple classical model of a ferroelectric thin film that reproduces the appearance of a FDP.We extended the model introduced in Ref. ([30]) by discretizing the thin film along the growth direction in each of its component unit cells in order to introduce a local spatial description of the polarization with the introduction of the local polarization, P i .Each individual crystal unit cell is comprised of two consecutive layers of BaO and TiO2, therefore P i could be obtained by adding each individual LP and dividing by the lattice constant in the stacking direction, c.The full polarization profile of the film is completely specified by the full set of LPs (p 1 , p 2 , ..., p 2M−1 , p 2M ) associated with each crystal plane of the oxide film (see Fig. 11).We can write the total energy of the thin film composed by M unit cells as: The LP at the interfaces (p L and p R ), are kept fixed.Two distinct interaction parameters θ L and θ R , establish the locality of the interface environment.[50] The interaction parameter θ ij for the internal dipoles is arbitrarily chosen to be unity.
where each term of the first sum accounts for the internal energy of the individual unit cells under zero electric field.The parameters a and b are obtained from an ab initio calculation of the total energy of a BaTiO 3 bulk for different ionic displacements along the soft mode.[51] The contribution from the depolarization field E is included in the second sum [22,30] where Ω is the volume of the unit cell.Each one of its terms favors energetically the local dipoles that follow the same direction than the DF.The third sum represents the classical electrostatic energy for a series of dipoles separated by a distance d ij , The above sum can be separated in three contributions: two that account for the interaction between the interface dipoles with the internal ones, and another one that accounts for the mutual electrostatic interaction between the internal dipoles.We have chosen to assign a fixed magnitude and direction to the interface layer dipoles (p L and p R ), with two distinct interaction parameters θ L and θ R , thus establishing the locality of the interface environment.[50] The interaction parameter θ ij for the internal dipoles is arbitrarily chosen to be unity.
With the constraint imposed by the simplicity of the model, we look for the energetically more favorable spatial polarization profile as a function of the film thickness.The spatial profile of polarization that characterize the film can be expressed as the vector (p L , p TiO 2 , p BaO , p TiO 2 , ..., p BaO , p TiO 2 , p R ) where a ferroelectric domain and a FDP will have all dipoles parallel or antiparallel respectively.Note that the total energy per unit cell (E ′ = E M ) an be expressed as function of only two single variables p TiO 2 and p BaO in the form: where we used Eqs.2324 and we define an overall interface parameter α = (θ R p R + θ L p L ) .Note that the thickness of the film is given by the number of unit cells M and directly affects the magnitude of the depolarization field E. The DF is calculated from an ab initio calculation for a supercell with the same size M and using a well-established averaging technique for the electrostatic potential inside the film.[50] As M → ∞ , the terms in Eq.25 that describe the mutual interaction between the dipoles converge to constant values, the DF vanishes and only the first two terms remain in E?, recovering the bulk properties.Only for small M these terms contribute appreciably to the total energy of the system.This is true in particular for the interface term that contains the parameter (fourth term in Eq. 25) that carries the information that defines the interfaces.We can estimate the interface parameters θ R , θ L , p L , p R from our ab initio simulations 2 and we assume that their numerical values will not depend upon the film thickness (as observed in the actual calculations, at least to first order).
As we vary the thickness of the ferroelectric film, we observe a change of the total energy landscape consequence of the interplay between the different energetic contributions.This is clearly seen in Fig. 12, where we show contour plots of E'(p TiO 2 , p BaO ) for different values of M. On the left, a diagram with the correspondent spatial distribution of dipoles is indicated.For thick films, the spatial distribution of dipoles that minimizes the total energy forms a domain (both LP's in the BaO and TiO2 planes have the same sign), as can be seen in the position of the minimum in the contour plot in Fig. 12A.As the number of unit cells is lowered, the minimum shifts to lower values of p BaO until at a critical thickness (Fig. 12B) it becomes zero.Further reduction of the number of layers M flips only the dipoles belonging to BaO planes thus establishing a FDP (see Fig. 12D).In conclusion, this simple toy model captures the general physics of the thin film, and illustrates the intimate relation between the interface characteristics, the thickness of the film, and the existence of an FDP.
Thus, by exploiting the concept of layer polarization in the description of ferroelectric thin films between metal contacts, we have been able to obtain detailed information on the modulation of polarization at the nanoscale and to understand the constraining effects of the interfaces in the determination of the ferroelectric response of the system.Our results shows that when film thicknesses reach a critical value, the ferroelectric system responds via a transition from the bulk ferroelectric structure to a ferrielectric antidipole pattern, where individual atomic layers acquire uncompensated opposing dipoles.This state arises as consequence of the complex energetic competition between the interface effects, the DF, and orientation and mutual interaction of the layer dipoles.The appearance of this particular ferrielectric state can be understood using a simple phenomenological model where the interface effects are explicitly taken into account.
These results suggest the possibility that such FDPs could be the normal state of a ferroelectric thin film at the nanoscale, even combined with the formation of two-dimensional island domains.It would be tempting to link the formation of such FDPs to the appearance or not of two-dimensional island below the critical thickness, and to understand the ultimate dependence upon the detailed interface structure and the nature of the metal contact.The next section will explore in detail the former point.

Tuning of polarization in metal-ferroelectric junctions
In this section we study the paradigmatic case of a BaTiO 3 film between Pt contacts, already introduced in the previous section, where we modify the ferroelectric properties by a selective control of the chemical species present at the interface.In particular, by inserting a single layer of a different metal (Au,Cu) we demonstrate that we are able to tune the polarization of the ferroelectric film via the modification of the screening properties of the composite metal contacts and through the change in the geometrical constraints at the interface.

Methods and discussion
We have simulated thin films of BaTiO 3 between Pt metal contacts in a (001) supercell.The oxide is terminated with a BaO plane at both interfaces, with the metal atoms directly bonded to the O atoms of the oxide planei,ii.The supercells constructed in this way can be labeled as Pt/(M)/(BaO − TiO2) m − BaO/(M)/Pt with m=1,2,4,6 and M=(Pt, Au or Cu).Nine atomic metal planes (including the intralayer) have been found to be sufficient to simulate the contacts under short-circuit boundary conditions.We used the experimental in-plane lattice parameter of BaTiO 3 (3.99Å),and kept it fixed in all calculations.All simulations have been performed using Density Functional Theory (DFT) within the Local Density Approximation, using ultrasoft pseudopotentials and a plane waves basis set.[52] A detailed analysis of the polarization at the nanoscale is critical for a complete description of the physical properties of the ferroelectric thin film.We have used Modern Theory of Polarization [6,8,18,53,54] and in particular the concept of layer polarization (LP) [20] to evaluate the polarization of the different structures and extract the information on the local profile of polarization at the nanoscale.
Turning to the results, in Fig. 13 we show the total polarization (normalized to the bulk value) in BaTiO 3 films of different thickness between composite metal contacts.These results clearly elucidate our claim: the polarization of the film is critically affected by the contact geometry at all the thicknesses we have considered and a residual polarization can be observed in films as thin as one unit cell.
In fact, a one unit cell thick BaTiO 3 film between plain Pt contacts is still in a ferroelectric state with a polarization 10% of the ideal bulk, and the introduction of a Au intralayer enhances this value up to 70% of that!In contrast, thicker (6 unit cells) oxide films display similar ferroelectric characteristics with 50% of the bulk polarization, a clear indication of the rapid decay of the interface effects with the film thickness.It is worth to note that at variance with the previous cases, a Cu intralayer induces a paraelectric (or almost paraelectric) behavior for all thicknesses.We will come back to this point later in the discussion.
In order to understand better the behavior of ferroelectricity in such ultrathin films, we computed the layer-resolved spatial profile of the polarization, which is quantified by values of the layer polarization along 3 the structure.This is shown in Fig. 14 for m4 films.When only Pt is present at the contact, the oxide develops a pattern of polarization composed of 3 The local polarization is defined as P j = p j c j for each layer j where c j is half the distance between neighboring cations and p i is the calculated layer polarization [20].For the outermost layers, we had to make a somewhat arbitrary choice for c j and used the distance between the outermost cation and the opposite metal plus half the distance between the cation and the one belonging to the second oxide layer.consecutive alternated signs along the direction perpendicular to the interface (see Section 5) From the behavior of the local polarization it is clear that while the pristine Pt (curve with black squares) interface display a distinct ferroelectric behavior, a single layer of Cu (red triangular dots) at the interface between BaTiO 3 and Pt is sufficient to stabilize the system in a paraelectric (non polar) structure (the character of the structure correlates with the symmetry of the polarization profile: an asymmetric pattern correspond to a polar geometry and hence to a ferroelectric behavior; a centro-symmetric pattern on the contrary gives rise to a paraelectric behavior).If we exchange Cu for Au (curve with green triangles), the polarization is restored, although smaller that in the Pt case, and the whole polarization profile is substantially changed.
These effects vary as the size of the film changes.In thinner films the local interface details have strong effects, while they are smoothed out in thicker systems.This can be seen in Fig. 15 A-D where we have plotted the LP profile for different sizes of the oxide and different metal intralayers at the interfaces.In general we notice that the Pt contact maintain the system in a ferroelectric state for all the sizes although the polarization profiles vary greatly with the system size.The same effect is observed with the addition of an Au intralayer.These changes are more noticeable in the thinner systems (m1 and m2) while for the thicker samples both systems share a similar polarization profile with a higher overall polarization in the Au case.This behavior can be directly associated with the screening of the depolarization field 4 at the metallic contacts.We have estimated the depolarization field from the macroscopic average of the electrostatic potential.[30] Indeed we observed a decrease of the DF with the introduction of the Au layer.This clearly demonstrate that Au has better screening properties compared to Pt and in consequence a Au intralayer enhances the overall polarization of the system.A similar reasoning is more difficult to do in the smaller systems due to the strong asymmetry introduced by the proximity of the interfaces and the uncertainty in the evaluation of the macroscopic averages.However, the qualitative behavior remains the same.
Contrary to Pt and Au, Cu intralayer do not stabilize any ferroelectric distortion at any thickness .This is a strong indication that Cu screening is weaker and not sufficient to reduce the depolarization field inside the oxide.
The above behavior is obviously correlated with the redistribution of charge across the interfaces in the different cases.The latter is quantified by the modification of the band alignment induced by the metal intralayer, and its influence on the screening properties of the metal contact.In fact, the knowledge of the band alignment, or Schottky barrier Height (SBH), allows us to define the properties of the interface phase between the metal and oxide.In a previous work, [51] we have demonstrated the correlation between the local structure of the interface composed by a crystalline oxide (BaO) and a d metal, and how we can tune the SBH by controlling the relative overlap of the local density of states of the different atoms close to the interface.Indeed we have found an almost complete similarity between that BaO/metal interface and the BaTiO 3 /metal interface of this work.In both cases we considered a BaO terminated oxide slabs, that show very similar characteristics (as also indicated by the fact that the interfaces between BaO and BaTiO 3 have almost zero valence band offset [52]).
Indeed, we find similar behaviors for the band offsets of the interfaces BaO/M/Pd (see Ref. [52]) and BaTiO 3 /M/Pt when Au, Cu and Pt are used as interlayer M. Furthermore, the SBH for the BaTiO 3 /metal interfaces follow the same ascending order for each metal intralayer (Au (0.8 eV) < Pt (1.1 eV) < Cu (1.4 eV)) in both systems (values given correspond to the case of m4) 5 This trend suggests a relation between the band offsets of the interfaces and the screening properties of the metals, that is to be expected given the strong correlation between the charge transfer at the interface and the SBH.[51] In particular, if a strong dipole is established at the interface, due to a high SBH, as in the case of Cu, fewer charges will be available for screening the depolarization field and the system will not support a ferroelectric distortion.Following Ref. [51] we can state the following phenomenological rule: metal intralayers that reduce the SBH at the interfaces will enhance the screening properties of the contact and stabilize the ferroelectric distortion.

Figure 1 .
Figure 1.Tetragonal ferroelectric structure of BaTiO 3 .Solid, shaded and empty circles represent Ba, Ti, O atoms, respectively.The arrows indicate the atomic displacements, exaggerated for clarity.The origin has been kept at the Ba site.Two structures, with polarization along [001], are shown.Application of a sufficiently large electric field causes the system to switch between the two states, reversing the polarization.in principle, we can assign one value of the Boolean algebra ("1" or "0") to each of the polarization states.

Figure 2 .
Figure 2. From Ref.[18] Mapping of the distributed charge density onto the centers of charge of the Wannier functions.

Figure 3 .
Figure 3. Schematic representation of the planar averaged electrostatic potential (fill line) of a slab with polarization perpendicular to the surface under (a)D = 0 boundary conditions (vanishing external field) and b) a vanishing internal electric field E = 0. Plannar averages are taken on the(x, y) planes parallel to the surface.Dashed lines represent an average over unit cell of the planar averages.An estimate of the macroscopic internal field inside the slab can be obtained from their slope.m stands for the dipole moment parallel to the surface normal and e for the electron charge.From Ref.[23]

Figure 4 .
Figure 4. Depolarization field for an isolated free standing slab in the absence of an external electric field (a), and for a ferroelectric capacitor with perfect screening (ideal metal) (b) and a real metallic electrode (c) under short circuit boundary conditions.From Ref.[1]

Figure 5 .
Figure 5. Evolution of the c/a ratio with the film thickness for monodomain PbTiO 3 films grown epitaxially on top of Nb-doped SrTiO 3 .With circles and squares, the experimental results.The dashed line is the phenomenological theory prediction.Solid lines correspond to the first principles model Hamiltonian results.From Ref. [35].

Figure 6 .
Figure 6.Theoretical predictions of the thickness dependence of the normal average polarization P, the tetragonality c/a, and the out of plane piezoelectric constant d 33 at room temperature for PbTiO 3 thin films grown on a SrRuO 3 /SrTiO 3 substrate.Values of the quantities at the bulk level are represented by dashed lines for the unstrained configuration, and with dotted lines for a geometry under the strain imposed by the substrate.The evolution of the system, from a monodomain configuration at large thickness, (where the depolarization field E d is small), to a 180• stripe domain structure in order to minimize the energy associated with E d , is represented in the inset.From[37]

Figure 7 .
Figure 7.We simulated ferroelectric films at different thicknesses, sandwiched between metal contacts.

Figure 8 .
Figure8.The thin film oxide develops a ferrielectric[47] pattern of polarization that exhibits antiferroelectric properties along the growth axis and with a net total polarization.This pattern would remain hidden by using more standard methods for calculating the local polarization.

Figure 9 .
Figure 9.The set of MLWFs so obtained cluster around each atomic layer in sets containing the correct number of Wannier function needed to keep the layer charge neutral.Here we can see the z position (horizontal lines lines ) of the centers of the Wannier functions and the position of the ions (indicated with circles).To calculate the LP of each layer, we use Eq.12.

Figure 10 .
Figure 10. A. Rumpling (cation-oxygen perpendicular distance (in Å) within each atomic layer) profile for six unit cells (m6) of BaTiO 3 sandwiched between Pt contacts.(on the horizontal axis Ti and Ba indicate the position of the individual TiO2 and BaO planes).B. Solid circles (blue on line): layer polarization (LP) profile for the same system m6 in units of 10 − 10C/m.The values for the equivalent planes in bulk BaTiO 3 are shown with black triangles.The orientation of the individual layer dipoles is displayed by arrows in the bottom of the panel.C.D.E.LP profiles for m4, m2 and m1.As the thickness is reduced the FDP is modified by the increasingly dominant interface effects.

Figure 11 .
Figure 11.Each individual crystal unit cell is comprised of two consecutive layers of BaO and TiO2.A local polarization P i can be associated with it.It is calculated by adding both individual LP's and dividing by the lattice constant in the stacking direction,c.We can associate a dipole p i A to each plane i, separated from its next neighbors by distance d ij .The full polarization profile of the film is completely specified by the full set of dipoles (p L , p 2 , ..., p 2M−1 , p R ).The LP at the interfaces (p L and p R ), are kept fixed.Two distinct interaction parameters θ L and θ R , establish the locality of the interface environment.[50]The interaction parameter θ ij for the internal dipoles is arbitrarily chosen to be unity.

Figure 12 .
Figure 12.Energy landscape E'(p TiO 2 , p BaO ) for different values of θ/M (film thickness).The interaction parameters for both interfaces were considered equal, theta R = θ L = θ and were kept fixed as the film size M is varied.The remaining parameters of the model have been derived from our ab initio calculations.A. θ/M = 0, bulk limit where the usual double well potential produces two equivalent minima with all the layer dipoles parallel.B.θ/M = 3, as the thickness is reduced the interface effects gain relevance and start to modify the energy landscape.The layer dipoles associated with the BaO planes reduce their magnitude in order to minimize the total energy.C. At θ/M = 5 the system reaches a critical thickness where these dipoles become zero and eventually flip orientation.D. For θ/M = 10 the system display a ferrielectric dipole pattern.

Figure 13 .
Figure 13.Normalized polarization for different number of oxide layers and interface phases.The local polarization is calculated for each oxide layer, added up and divided by the BaTiO 3 bulk polarization calculated in the same way with equivalent numbers of layers.

Figure 14 .
Figure 14.Local Polarization for a unit cell of size m4.Only one layer of Cu at the interface betwen the BaTiO 3 and Pt is enough to drive the system from a ferroelectric structure (squared dots) into a paraelectric one (red triangular dots).

Figure 15 .
Figure 15.Layer polarization for different metal interfaces and films thicknesses A)m1, B)m2, C)m4, D)m6.As reference, the LPs for bulk BaTiO 3 are indicated with blue circles and dotted lines for each case.Metals intralayers used are Cu, Au, Pt.Case of m6 with Cu layer: this result is in fact an artificial effect due to the known LDA underestimation of the band gap.The Fermi level of the system overlaps the conduction band and starts filling states that should be empty.We expect the same issue for m > m6, or with the use of metals with a higher work function.The calculations for m < m6 are be correct under the DFT?LDA framework, as the Fermi level is safely below the conduction band.