First-principles study of two-dimensional electron and hole gases at the head-to-head and tail-to-tail 180$^\circ$ domain walls in PbTiO$_{3}$ ferroelectric thin films

We study from first-principles the structural and electronic properties of head-to-head (HH) and tail-to-tail (TT) 180$^\circ$ domain walls in isolated free-standing PbTiO$_{3}$ slabs. For sufficiently thick domains ($n$ = 16 unit cells of PbTiO$_{3}$), a transfer of charge from the free surfaces to the domain walls to form localized electron (in the HH) and hole (in the TT) gases in order to screen the bound polarization charges is observed. The electrostatic driving force behind this electronic reconstruction is clearly visible from the perfect match between the smoothed free charge densities and the bound charge distribution, computed from a finite difference of the polarization profile obtained after the relaxation of the lattice degrees of freedom. The domain wall widths, of around six unit cells, are larger than in the conventional 180$^\circ$ neutral configurations. Since no oxygen vacancies, defects or dopant atoms are introduced in our simulations, all the previous physical quantities are the intrinsic limits of the system. Our results support the existence of an extra source of charge at the domains walls to explain the enhancement of the conductivity observed in some domains walls of prototypical, insulating in bulk, perovskite oxides.


I. INTRODUCTION
Domain walls (DWs) in oxide ferroelectric materials, defined as nanometer-scale interfaces that separate regions with different orientations of the spontaneous electric polarization, have attracted a lot of attention during the last few years due to the novel functional properties they might display, radically departing from those observed in the parent materials 1 . One of these surprises comes from the observation of a room-temperature electronic conductivity at ferroelectric DWs in a thin-film insulating multiferroic BiFeO 3 2 . Following this work, Guyonnet and co-workers 3 demostrated that 180 • DWs in much simpler, purely ferroelectric tetragonal perovskite Pb(Zr 0.2 Ti 0.8 )O 3 thin films were also conducting, suggesting that this phenomenon is more general than initially envisaged. Indeed, the presence of a conductive domain wall is not restricted to thin films, but equally applies to millimeterthick widebandgap single crystals, such as LiNbO 3 4 or WO 3 5 . The microscopic origin of the conduction is controversial. It was first ascribed to structural changes at the DW that produced a polarization discontinuity, leading to steps in the electrostatic potential, and a concomitant lowering of the band gaps 2 . More recently, it has been suggested that the conductivity could be linked with extrinsic factors, such as oxygen vacancies 3,6 .
But other sources for a metallic type-conductivity at DWs have been proposed. Sluka et al. 7,8 observed an enhancement of conductivity up to 10 9 times that of the parent matrix in charged DWs of the prototypical ferro-electric BaTiO 3 . This was attributed to the presence of a stable degenerate electron gas that would be formed to screen polarization divergences in head-to-head (HH) and tail-to-tail (TT) DWs 9 . In conventional 180 • domains in tetragonal ferroelectrics, the walls run parallel to the tetragonal polarization axis. In such a configuration, the polarization charge (defined as the dot product of the polarization times the unitary vector normal to the DW) vanishes, yielding to a neutral configuration. However, tilting of the DWs would lead to polarization-induced bound charges that would result in a large and unfavorable electrostatic energy 10 . The extreme case, where the polarization points in the same direction as the normal vector to the wall, leads to the formation of head-tohead DWs (if the polarization is directed towards the DW; positive bound charges), or tail-to-tail DWs (if the polarization is directed outwards from the DW; negative bound charges). The properties of charged DWs might be very different from those of the neutral configurations. Different mechanisms have already been postulated in the literature to explain the screening of the polarization bound charges and stabilization of HH and TT 180 • DWs in ferroelectrics. The theoretical background for the formation of charged domain walls in proper ferroelectrics based on a phenomenological model was established in Ref. 8. Phase field simulations support the existence of a degenerate quasi-two-dimensional electron gas, localized due to the potential well formed by the polarization charges 7 . The application of a Landau model to address the problem of charged DWs in an isolated ferroelectric, including the polarization profile, width and formation energy of the domains was undertaken in Ref. 9. From a first-principles perspective, Wu and Vanderbilt 11 proposed a periodically repeated chemical "delta" doping in PbTiO 3 superlattices, where the Ti atoms at certain layers are replaced by donor (Nb) or acceptor (Sc) atoms drawn from neighboring columns of the Periodic Table. Rahmanizadeh et al. 12 suggested for the same material (i) the formation of a conducting layer at the domain wall in a superlattice geometry; (ii) the appearance of polarons due to the localization of one electron on a Ti ion (requiring for this the addition of an on-site Coulomb repulsion term) also in a superlattice geometry; and (iii) the presence of defects with various amounts of concentrations in thin film geometries. But in all the former first-principles simulations, the presence of substitutional atoms at one or the two HH and TT domain walls in the superlattices, or a varying defect concentration at the tails in the thin film geometry were assumed, and this fully determines the electric displacement, and therefore the polarization, within the domain as was shown in Ref. 13. Questions like what is the intrinsic critical thickness for the stabilization of the two-dimensional conducting layers for the screening of the polarization charges, the spatial extension of the hole and electron gases, or how large the polarization would be, if the HH and TT domain walls are induced in an ideal slab are not accessible from the previous computations. Similar issues regarding the screening of polarization discontinuities with electronic reconstructions have been the subject of intense studies in the celebrated LaAlO 3 /SrTiO 3 interfaces (Ref. 13 and references therein), and theoretically in KNbO 3 /ATiO 3 (A=Sr,Ba,Pb) superlattices 14,15 . They have been also faced by Aguado-Puente et al. 16,17 in a similar system to ours: an interface between a ferroelectric material like PbTiO 3 and a dielectric like SrTiO 3 . Above a critical thickness of 16 unit cells, the ferroelectric layers could be stabilized in the out-of-plane monodomain configuration due to the electrostatic screening provided by the free carriers. We wonder whether the same results hold when the dielectric is replaced by an opposite domain, doubling the amount of polarization charge at the DW.
The rest of the paper is organized as follows: In Sec. II we describe the method on which the simulations are based. In Sec. III we present the results on the polarization profiles (Sec. III A), the electronic properties including the layer-by-layer projected density of states (Sec. III B), and the transfer of charge that yields to the formation of the two-dimensional hole and electron gases (Sec. III C), together with the coupling between the electronic properties (density of free charge) and the structural properties (divergence of the polarization). Future perspectives are summarized in Sec. IV.

II. METHODS
We carried out simulations of HH and TT domain walls in a free standing slab of PbTiO 3 , within a numerical  atomic orbital method as implemented in the Siesta code 18 . Exchange and correlation were treated within the local density approximation (LDA) 19 to the density functional theory (DFT) 20,21 . Since we were interested in the intrinsic properties of the domains, no Hubbard-term to deal with the on-site Coulomb repulsion on the Ti d states was considered. Such a term would reduce the Ti 3d-O 2p hybridization, at the origin of the ferroelectric behaviour of PbTiO 3 , and would favour a localization of the electrons at the interface, that would strongly depend on the value of the U -term 12 . Core electrons were replaced by ab initio norm conserving pseudopotentials, generated using the Troullier-Martins scheme 22 in the Kleinman-Bylander fully nonlocal separable representation 23 . Due to the large overlap between the semicore and valence states, the semicore 3s and 3p electrons of Ti, and the 5d electrons of Pb were considered as valence electrons and explicitly included in the simulations. Ti and Pb pseudopotentials were generated scalar relativistically. The reference configuration and cutoff radii for each angular momentum shell for the pseudopotentials used in this work can be found in Ref. 24 for Ti and O, and in Table I for Pb.
The one-electron Kohn-Sham eigenstates were expanded in a basis of strictly localized numerical atomic orbitals 25,26 . We used a single-ζ basis set for the semicore states of Ti and Pb, and double-ζ plus polarization for the valence states of all the atoms. For Pb, an extra shell of 5f orbitals was added. All the parameters that define the shape and range of the basis functions were obtained by a variational optimization of the energy 27 in bulk cubic BaTiO 3 (for Ti, and O). For Pb, another optimization was performed in bulk PbTiO 3 , frozen in the atomic orbitals of Ti and O to those previously optimized in BaTiO 3 28 . The electronic density, Hartree, and exchangecorrelation potentials, as well as the corresponding matrix elements between orbitals, were calculated in a uniform real-space grid, with an equivalent plane-wave cutoff of 1200 Ry in the representation of charge density. For the Brillouin zone integrations we use a Monkhorst-Pack sampling 29 equivalent to 6 × 6 × 6 in a five atom perovskite unit cell. A Fermi-Dirac distribution was chosen for the occupation of the one-particle Kohn-Sham electronic eigenstates, with a smearing temperature of 250 K for the HH and 500 K for the TT domain structures.
To simulate both domain structures we used a tetragonal (1×1) supercell, periodically repeated in space, of the type vacuum/(PbO-TiO 2 ) 2n -PbO/vacuum with n = 16, as shown in Fig. 1 (simulations with n=12 showed that the system went back to the paraelectric structure without stabilizing the metallic gases). Both surfaces were terminated with a PbO layer, since this is the only stable surface termination found by first-principles 30 . The vacuum thickness was equivalent to roughly eight unit cells of the perovskite. A dipole correction was introduced to avoid spurious interaction between periodic images of the slab in the out-of-plane direction. In order to simulate the effect of the mechanical boundary conditions due to the strain imposed by an hypothetical substrate, the in-plane lattice constant was fixed to the theoretical equilibrium lattice constant of bulk SrTiO 3 (a 0 = 3.874 A), one of the most common substrates used to grow oxide heterostructures. In constrained bulk PbTiO 3 under this compressive strain, the spontaneous polarization amounts to P S = 77.1 µC/cm 2 . As the initial point, an ideal structure was defined stacking along the [001] direction 2n unit cells of PbTiO 3 in a tetragonal paraelectric configuration (i. e. zero rampling in all the atomic layers) with a tetragonality of c/a = 1.0325. Then, on top of this paraelectric structure, the PbTiO 3 atoms were moved in opposite directions in the two halves of the slab geometry following the displacement pattern of the bulk tetragonal soft mode scaled to 62%. That is the value of the polarization of a meta-stable structure of a monodomain PbTiO 3 slab where the polarization charge can be screened by an electronic reconstruction and the formation of two-dimensional metallic gases at the surfaces 16 . By construction, the central PbO layer is a mirror symmetry plane. Starting from this geometry, a conjugate gradient minimization was performed till the maximum component of the force on any atom was smaller than 0.04 eV/Å.
Once the atomic structure is relaxed and the oneparticle density matrix converged, in order to compute the density of states, a non-self-consistent calculation was carried out with a much denser sampling of 150 × 150 × 5 Monkhorst-Pack mesh.
To establish the notation, we shall call the plane parallel to the interface the (x, y) plane, whereas the perpendicular direction will be referred to as the z axis.   After relaxation of the slabs with n = 16 unit cells in each domain, we observe how the two polar configurations are meta-stable, while the most stable phase is the centrosymmetric unpolarized structure. A similar result was found in the PbTiO 3 /SrTiO 3 interfaces 17 . The TT configuration is more stable than the HH, as can be seen from the differences in energy with respect the centrosymmetric paraelectric phase shown in Table II.

A. Polarization profiles
In Fig. 2 we plot the layer-by-layer effective polarization profile in the HH [ Fig. 2(a)] and TT [ Fig. 2(b)] domain structures for PbTiO 3 . The layer polarizations have been calculated using the approximate expression based on the re-normalized bulk dynamical charges, following the recipe given in Ref. 31. The polarization pro- where δ is the typical half-width of the domain wall (i.e. to make the transition from −P 0 to P 0 we need a space of four times this length δ. The results of the fitting of the profiles shown in Fig. 2 to the model summarized in Eq. (1) can be found in Table II. Clearly, the domain wall widths are larger than in the neutral configuration, where they were found to be very narrow, of the order of the lattice constant 32 , also supporting the conclusions of Ref. 9.
These 180 • domain walls, where the wall is perpendicular to the polarization, give rise to polarization-induced bound charges, that lead to the local concentration of a large electrostatic energy that destabilizes this configuration. The stabilization of the domains observed in Fig. 2 implies that the positive (for the HH) and negative (for the TT) polarization charges at the domain wall have been screened. Since in our calculations we do not consider the presence of defects or dopants 11 , the most plausible scenario is the metallization of the domain wall and the free surfaces due to the neutralization by free carriers, yielding to the formation of thin quasi-two-dimensional metallic layers. The fingerprint to characterize such scenario would be the presence of a substantial amount of free charge populating the band edges (top of the valence band and bottom of the conduction band) both at the domain wall and at the free surfaces of the slab. Fig. 3 is adjacent to the free surface and the top one lies close to the domain wall. Only the PDOS for one of the domains are shown, since there is a mirror symmetry plane at the domain wall. All the PDOS curves were calculated following the recipe given in Sec. III.A.1 of Ref. 31 with the Dirac delta functions for the eigenvalues in the PDOS computations replaced by smearing normalized Gaussians with a finite width that is twice as large, as suggested in Appendix B of Ref. 31. On top of the heterostructure PDOS we superimpose the bulk PDOS (dashed red lines), calculated with an equivalent k-point sampling. The bulk reference calculation was computed from a periodic bulk calculation with a five atom per unit cell PbTiO 3 structure, with the atomic distortions and out-of-plane strain extracted from a region of the HH or TT domains where the relaxed atomic structure (Fig. 2) has converged into a regular pattern. Then we extract the PDOS for the same TiO 2 atomic layer. Finally,we note that the superposition of the bulk PDOS and supercell PDOS at each layer was done by matching the sharp peaks of the Ti(3s) semicore band. First of all, it is remarkable to check how for the per-ovskite material studied in this work, PbTiO 3 , and even upon the electronic reconstruction and metallization, a well-defined energy gap persists in all layers between the bands that are mostly O-p in character (top of the valence band at the bulk level), and the bands that are mostly Ti-t 2g in character (bottom of the conduction band at the bulk level). The local metallization of the free-surface layers and the domain wall is evident from the location of the valence-band top and the conduction-band bottom of these layers with respect to the Fermi level of the whole structure, taken as zero in Fig. 3. For the HH domain configuration [ Fig. 3(a)], the non-vanishing PDOS at the Fermi level for the three topmost TiO 2 layers points to the electronic population of the conduction-band bottom of the PbTiO 3 unit cells close to the domain wall. This free electronic charge has been transferred from the free surface, as can be clearly seen from the generation of holes at the bottom-most layers (see the crossing of the Fermi level with the PDOS at the valence-band top for the three bottom-most TiO 2 layers). For the TT domain configuration [ Fig. 3(b)], the same behaviour is found, but with the role of the electrons and holes interchanged (i.e. generation of holes at the top of the valence band in the layers adjacent to the domain wall, and the population of the bottom of the conduction band with electrons close to the free surfaces).

Fig. 3 shows the layer (spatially) resolved PDOS on all the atomic orbitals of the Ti and O atoms at the TiO 2 layers of the different PbTiO 3 unit cells. The bottom layer plotted in
The PDOS of the conduction and valence bands converges fairly quickly to the bulk curve when moving away from the interface, and the PDOS vanish at the Fermi level, which implies that the system is locally insulating. Furthermore, the PDOS in each layer appears rigidly shifted with respect to the neighboring two layers, consistent with the presence of an internal remnant depolarization field resulting from the incomplete screening within the domains (see black dashed lines in Fig. 3). Since the bulk calculations were performed at zero external electric field, while in the slab calculations there is a residual depolarizing field, the coincidence of the two curves far enough from the free surfaces and the domain wall means that the main effect on the layer by layer PDOS comes from the lattice distortions, and not from macroscopic depolarizing fields. This residual depolarizing field will play a major role in the determination of the atomic orbitals involved in the electronic reconstruction, as will be discussed below.

C. Coupling between density of free carriers and the polarization profiles
To estimate the amount of charge that has been transferred from the free-surfaces to the domain wall in order to screen the depolarizing fields, we plot in Fig. 4 the planar average of the free charge, ρ free (z), defined as  where S is the area of the interface unit cell, and the corresponding nanosmoothed version 33 , ρ free (z). For the electronic charge density populating the bottom of the conduction band ρ free (r) has been computed as while for the hole charge density populating the top of the valence band, it amounts to where w k is the weight of a special k point in the discrete Monkhorst-Pack mesh, f nk are the occupation numbers of the eigenstate ψ nk and the sum is restricted to the states with eigenvalue E nk higher than E 0 , a value of the energy corresponding to the center of the gap between valence and conduction band. The results are shown in Fig. 4. Again, we clearly see that the system is not locally charge neutral at the domain wall and at the freesurfaces. If the profiles of the nanosmoothed average free charge densities are integrated along z, we can extract the free-carrier per unit area, labeled as σ e (for the electrons) and σ h (for the holes). For the HH domains, there is an additional electron density at the domain wall populating the bottom of the conduction bands that amounts to σ HH e = 1.033 electrons [ Fig. 4(a)]. Those electrons are transferred from the free-surfaces, where the integration of the free hole charge in Fig. 4 A further insight on the origin of these electron and hole free charges can be drawn after its decomposition into different orbital components, as it is done in Fig. 5. The current wisdom is that the d xy orbitals are the first to be populated in the related LaAlO 3 /SrTiO 3 34 or PbTiO 3 /SrTiO 3 17 interfaces, where two-dimensional electron gases are formed to screen the polar catastrophe. This situation is similar to that found in bulk polarized PbTiO 3 and also observed in the free surface of our system [ Fig. 5(c)]. However, when discussing the metallic states formed around the domain walls we can observe, on the one hand, that the HH electron gas (Fig. 5b) is, for the most part, formed by d xz and d yz orbitals, while d xy ones host very little charge. This can also be observed in the doubled-peaks of Fig. 4b. On the other hand, the TT hole gas has a primary contribution from the π-orbitals in the PbO plane (i.e. p x and p y orbitals in the oxygen at the PbO layer) and a smaller, but also significant, contribution from both σ (i.e. the p x , respectively p y , orbital of the equatorial O located along x, respectively y, with respect the Ti) and π-bonded orbitals (i.e. the p y and p z , respectively p x and p z , orbitals of the equatorial O located along x, respectively y, with respect the Ti) in TiO 2 planes. Both facts are consistent with the effects of localization due to the remnant electric potential energy within the domains, whose shape resembles a ∨ for the HH and ∧ for the TT, with the peak located right at the domain wall in both cases. Those are confinement potentials for the respective carriers. For the bands with majority weights on orbitals directed parallel to the domain wall, there is no effect on the dispersion; only the center of the subbands might be displaced following exactly the shape of the potential. This is the case of the d xy bands of Ti, and on the O-σ orbitals. However, the bands with majority weights on orbitals directed perpendicular to the domain wall, both the dispersion and the local center of mass of the bands are strongly modified due to the confinement potential: the range of the possible interactions is strongly reduced. As a consequence the hybridization is decreased, and the corresponding bonding orbitals (O-π orbitals in PbO layers) increase the energy and the antibonding orbitals (d xz ,d yz ) suffer a concomitant reduction. Those are the first orbitals to be filled with electrons [ Fig. 5(b)] or emptied to create holes [ Fig. 5(d)] around the domain wall.
In order to elucidate how this electronic reconstruction is a result of the screening of the polarization-induced charges at the 180 • domain walls, we also plot in Fig. 4 a numerical differentiation of the polarization profile shown in Fig. 2. There is an essentially perfect matching between the nanosmoothed free charge densities and the bound charge profile, ρ b (blue lines in Fig. 4), defined as ρ b (z) = −dP/dz. Therefore, the polarization discontinuity at the domain wall is the driving force for the electronic reconstruction, and this transfer of charge between the domain wall and the free surface provides an effective screening mechanism that stabilizes the monodomain configurations in both domains.

IV. CONCLUSIONS
Our density-functional simulations have shown how 180 • HH and TT charged domains can be stabilized in the out-of-plane monodomain configuration in an ideal PbTiO 3 slab (no dopants nor vacancies), thanks to the electrostatic screening provided by free carriers. Above a critical size of the domains, in our case n = 16 unit cells, an electronic reconstruction takes place, where charge is transferred from the free-surface of the slab to the domain wall, and generates quasi-two-dimensional electron and hole gases. The width of the domain wall of around 6 (for the HH) and 7 (for the TT) unit cells, is much larger than the typical width of neutral domains (of the order of one lattice constant). A clear electrostatic coupling between the density of free-charge and the polarization profile, whose negative gradient equals the bound charge, is found.
Our results support previous phenomenological models based on Landau and semiconductor theory, regarding the existence of quasi-two-dimensional charge densities at the 180 • HH-TT domain walls in isolated ferroelectric 9 , and postulated them to be the driving force for the metallic like conductivity found at the walls between insulating polar materials 7 .
These first-principles simulations might serve as benchmarks to perform second-principles computations that include all the relevant electronic and lattice degrees of freedom 35 in very large systems under operating conditions of field and temperature. Once the atomistic and electronic models for the second-principles calculations will be validated, then more challenging calculations including tilting domains like the ones that appear during domain switching, or the conductivity at domain walls can be tackled.