First-principles study of (Ba,Ca)TiO$_3$ and Ba(Ti,Zr)O$_3$ solid solutions
Danila Amoroso Andres Cano Philippe Ghosez
cond-mat.mtrl-sci





Physique Théorique des Matériaux, Q-MAT, CESAM, Université de Liège (B5), B-4000 Liège, Belgium
ICMCB, CNRS, Université de Bordeaux, UPR 9048, F-33600 Pessac, France

ICMCB, CNRS, Université de Bordeaux, UPR 9048, F-33600 Pessac, France

Physique Théorique des Matériaux, Q-MAT, CESAM, Université de Liège (B5), B-4000 Liège, Belgium





(Ba,Ca)TiO_3 and Ba(Ti,Zr)O_3 solid solutions are the building blocks of lead-free piezoelectric materials that attract a renewed interest.
We investigate the properties of these systems by means of first-principles calculations, with a focus on the lattice dynamics and the competition between different ferroelectric phases. We first analyze the four parent compounds in order to compare their properties and their different tendency towards ferroelectricity. 
The core of our study is systematic characterization of the binary systems (Ba,Ca)TiO_3 and Ba(Ti,Zr)O_3 within both the Virtual Crystal Approximation and direct supercell calculations. 
In the case of Ca doping, we find a gradual transformation from A-site to B-site ferroelectricity due to steric effects that largely determines the behavior of the system. In the case of Zr doping, in contrast, the behavior is eventually dominated by cooperative Zr-Ti motions and the local electrostatics.  
In addition, our comparative study reveals that the specific microscopic physics of these solid sets severe limits to the applicability of the Virtual Crystal Approximation for these systems.



First-principles study of (Ba,Ca)TiO_3 and Ba(Ti,Zr)O_3 solid solutions
    Philippe Ghosez
    December 30, 2023
=======================================================================





§ INTRODUCTION


The potential interest of BaTiO_3-based solid solutions for piezoelectric applications has been known since the sixties,<cit.> although actual piezoelectric devices mainly rely on lead zirconate titanate. However, 
the interest in alternative lead-free solid solutions has been renewed since the 2003 european directive on the restriction of the use of certain hazardous substances in electrical and electronic 
equipment.<cit.> 
Among all the possible combinations, 
the partial homovalent substitution of Ba by Ca and Ti by Zr in BaTiO_3 
(Ba_1-xCa_xTi_1-yZr_y hereafter BCTZ) is of special interest.
In their seminal paper, Liu and Ren<cit.> reported a high piezoelectric coefficient of 620 pC/N for the ceramic system 
(1-x)Ba(Ti_0.80Zr_0.20)O_3-x(Ba_0.70Zr_0.30). This has been ascribed to the presence of an special point in the phase diagram in which the tetragonal and rhombohedral ferroelectric phases meet the paraelectric cubic one. Subsequently, 
Keeble et al.<cit.> revisited this phase diagram 
and observed an intermediate ferroelectric orthorhombic phase that also meets the three previous ones. 

As in the case of Pb-based system, the achievement of high piezoelectric response in BCTZ is believed to be linked to the existence of a so-called phase convergence region, 
in which the lack of an energy barrier between different ferroelectric states makes the landscape isotropic and polarization free to rotate<cit.>. 

Despite a great number of experimental studies have been reported,<cit.> the theoretical investigation of BCTZ compounds remain comparatively very limited and, to the best of our knowledge, no comprhensive study 
based on direct Density Functional Theory (DFT) has been carried out. In order to fill this gap and to clarify the microscopic mechanisms involved 
in the appearence and competition of different polar phases in 
(Ba,Ca)(Ti,Zr)O_3 solid solutions, a systematic ab-initio investigation is reported here.
Specifically, we report a first-principles study of the structural and dynamical properties of different compositions of 
Ba_1-xCa_xTiO_3 and BaTi_1-yZr_yO_3 solid solutions and 
related parent compounds. 
We perform calculations using the virtual crystal approximation (VCA) 
and also directly on supercells. The comparison of the results allows us to determine the applicability of VCA for these particular solids solutions. In addition, we also study role of the local atomic arrangements by considering various BaZrO_3/mBaTiO_3 supercells.

The paper is organized as follows. In Sec. <ref> we introduce the technical details of our DFT calculations. Then, 
in Sec. <ref>, we carefully characterize the origin of the main instabilities via the analysis of the interatomic force constants 
and the energetics of several possible stable and metastable phases in the four parent compounds. 
Next, in Sec. <ref> and Sec <ref>, we address the two A- and B- site substitutions separately in (Ba,Ca)TiO_3 and Ba(Ti,Zr)O_3 binary-systems, respectively. 
Finally, global discussion and conclusions are reported in Sec. <ref> and Sec. <ref>. 



§ TECHNICAL DETAILS 





 §.§ Computational Details


We use the general framework of the density functional theory (DFT) to compute the structural properties and the electric polarization. 
Dynamical and piezoelectric properties have been calculated using 
Density Functional Perturbation Theory (DFPT)<cit.>, as implemented in the ABINIT package<cit.>. 
The exchange-correlation potential was evaluated within the generalized gradient approximation (GGA) using the Wu-Cohen (WC) functional<cit.> for all simulations.
Optimized Norm-Conserving Pseudopotentials<cit.> have been employed with the following orbitals considered as the valence states:
5s, 5p and 6s for Ba, 2s, 2p and 4s for Ca, 3s, 3p, 4s and 3d for Ti, 4s, 4p, 4d and 5s for Zr and 2s and 2p for O.
The energy cutoff for the expansion of the electronic wavefunctions has been fixed to 45 Ha. 
Before performing computationally demanding calculations on ordered-alloy supercells by using a standard DFT approach, a first investigation has been addressed
by using the Virtual Crystal Approximation (VCA) within the straightforward combination of the reference atomic pseudopotentials<cit.>.
For the optimization of the cubic pervoskite structures, phonons and polarization calculations within VCA, we used a 6x6x6 k-points mesh for the Brillouin zone sampling 
for which energy is converged to 0.5 meV, 
whereas for the optimization of the polar P4mm, R3m and Amm2 and piezoelectric repsonse calculations different sampling from 8x8x8 to 8x6x6 were used.
For the optimization of the supercells and for the associated phonons calculations we used the 8x8x8 k-mesh. For all the tetragonal superlattices we used 8x8x6 k-mesh sampling. 
In order to allow comparison between different structures, energy will be reported in meV/f.u. (i.e. per 5 atoms).
The q-points for the phonon dispersion curves and IFCs of the four pure compounds included Γ, X, M, R and the Λ point halfway from Γ to R of the 
simple cubic Brillouin zone.



 §.§ Structure of the parent compounds


The optimized the lattice parameters of the cubic structure are reported in Table <ref>. These values are used in Sec. <ref> for the calculation of the phonon dispersion curves and interatomic force constants. The obtained values are in excellent agreement with experimental data for BaTiO_3 and BaZrO_3 <cit.>
(within 1%), whereas for CaTiO_3 the understimation is of about 2% <cit.>. The larger error for CaTiO_3 can be assigned to the fact that 
the 0 K DFT results are compared to measurements at 1720 K.
For CaZrO_3 the optimized lattice parameter overstimates the experimental value <cit.> of about 2%. 

Atoms in the ABO_3 perovskite structure are labeled according to Fig <ref>.



 §.§ Supercell structures







We have considered different supercells describing the Ca and Zr doping of BaTiO_3 and performed DFPT calculations on the corresponding relaxed high-symmetry
paraelectric structures.

For Ba_7/8Ca_1/8TiO_3 and BaTi_7/8Zr_1/8O_3, we have used the smallest reference structure which is the
cubic Pm3̅m (O^h_1) with 40 atoms. The GGA-WC optimized lattice constants are 7.917Å  and 8.003Å  respectively (Fig. <ref>).

For Ba_0.50Ca_0.50TiO_3 and BaTi_0.50Zr_0.50O_3
three distinct geometric arrangements
of the 40-atoms cell were considered:
columns of same cations along the [001] direction, layers of same cations parallel to the [110] planes and rocksalt configuration, Fig. <ref>.
The paraelectric reference of the two first “ordered" structures are characterized by a tetragonal symmetry within the P4/mmm (D^4h_1) space-group and they can be reduced by symmetry to
10-atoms cell, whereas the rocksalt configuraton has the cubic symmetry within the Fm3̅m (O^h_5) space-group.
Structural relaxation have been done on the 2x2x2 supercells (40 atoms) for each geometrical ordering, whereas DFPT calculations 
and structural relaxation of lower symmetry structures 
have been performed on the reduced 10-atoms cell for the P4/mmm structures and on the 40-atoms cell for the Fm3̅m one.
The optimized lattice parameters of the high-symmetry references are reported in Table <ref>. 
Figures have been produced by using the VESTA package<cit.>. 



§ PARENT COMPOUNDS



We start our study by considering individually the BaTiO_3, CaTiO_3, CaZrO_3 and BaZrO_3 parent compounds. 
The dynamics and the energetics of these systems is fundamental to understand the key features of corresponding solid solutions.

BaTiO_3 is one of the most studied perovskite, both from the theoretical <cit.> and 
experimental <cit.> points of view. This perovskite is characterized by a tolerance factor t greater than 1, t≃ 1.06<cit.>, 
that 
allows to predict a polar distorted ground-state <cit.>. In fact, while stable at high temperature in the centrosymmetric cubic (Pm3̅m) phase, it 
undergoes <cit.> ferroelectric structural phase transitions to a tetragonal (P4mm) structure at ≃ 393 K, to a orthorhombic (Amm2) phase at ≃ 278 K 
and to a rhombohedral (R3m) ground-state at 183 K. 

CaTiO_3, on the contrary, has a tolerance factor t smaller than 1, t≃ 0.97<cit.>. Accordingly, this compound displays a non-polar orthorhombic (Pnma) ground-state.
 It  exhibits at least two observed phase transitions at high temperatures: from a cubic (Pm3̅m) to a tetragonal phase (I4/mcm) at ≃ 1634 K and from a 
tetragonal (I4/mcm) to the orthorhombic ground-state (Pnma) 
at ∼ 1486 K <cit.>.

CaZrO_3, that has t≃ 0.91, was observed only in two different structures: 
the high temperature cubic (Pm3̅m) form and its orthorhombic (Pnma) ground-state with the transition temperature at ≃ 2170 K.<cit.>

BaZrO_3, with 
t≃ 1.00<cit.>, is 
experimentally not known to undergo any structural phase transition 
and remains in the paraelectric cubic phase down to 2 K. <cit.>

In Section <ref>, we provide a comparison of the dynamical properties (phonon dispersion curves and interatomic force constants) 
of the parent compounds in their cubic phase. 
Then, in Section <ref>, we compare the energetics of various metastable phases of lower symmetry arising from the condensation
of individual and combined distortions related to the unstable phonon modes identified in Sec. <ref>.



 §.§ Phonon dispersion curves and Interatomic Force Constants





The identification of imaginary phonon frequencies and the corresponding displacement patterns allow us to pinpoint the main instabilities behind the structural phase transitions in our systems.
Thus, we considered the cubic reference structure and we computed the phonon dispersion curves along selected high-symmetry lines in the simple Brillouin zone. 
The results are shown in Fig. <ref>. 
The presence of imaginary phonon frequencies (shown as negative values in Fig. <ref>) reveals the structural instabilities of the cubic phase.
The nature of the corresponding transition is usually determined by the character of the main unstable modes and the related eigendisplacement vectors indicate the atomic displacements
that spontaneously appear to reach the most stable configuration.

We complement this analysis with the calculation of the interatomic force constanst in real space. These constans are defined as 
C_kα,k'β(l,l')=∂^2E_tot/∂τ^l_kα∂τ^l'_k'β, where E_tot is the total energy of a period crystal 
and τ_kα^l is the displacement of the atom k in the cell l along direction α from its equilibrium position.
For a pair of distinct atoms, the IFC can be interpreted as minus the harmonic spring constant between them,
so that negative values correspond to stable interactions. For the “on-site” forces, on the contrary, positive values correspond to stable interactions.
For a more detailed description of the physical meaning see Refs. <cit.>. 
In our analysis, we further separate the contribution of the dipole-dipole interaction (DD) from that of
the short-range forces (SR) in order to identify the key mechanisms that lead the system to display or not specific phonon instabilities.

In the following we address a systematic description of the phonon dispersion curves reproduced in Fig. <ref>. 






  §.§.§ BaTiO_3



Ferroelectricity in BaTiO_3 is known to be related to the Ti off-centering driven by an unstable polar mode at Γ <cit.>. 
This ferroelectric (FE) instability is such that it expand over the entire Γ-X-M planes of the simple cubic Brillouin zone, as can be seen in Fig. <ref>(a). 
While each atom is at a position stable against individual displacements (Table <ref>),
the origin of this distortion has to be primarily ascribed to the destabilizing Ti-O interaction, reflected in the positive value of the interatomic force constant along the 
bond's direction (Table <ref>). An additional contribution comes from the strong interaction of subsequent Ti-atoms along the B-B' chain direction compared to the small 
value in the transverse component. 
These anisotropic couplings give rise to dispersion of the unstable mode almost flat in the Γ-X-M plane and highly dispersive along Γ-R. This
reflects the chain-like nature of the instability in real space. The polar distortion requires cooperative atomic displacement along Ti-O chains<cit.>.
Furthermore, the negligible contribution of Ba-atom to the ferroelectric distortion, has to be ascribed to its sizable on-site force constant and very weak Ba-O_1 coupling.

In Table <ref>, we also compare the eigendiplacements associated to the unstable F_1u(TO_1) mode to those of the 
stable F_1u(TO_2) and F_1u(TO_3) ones and we report the related mode effective charges, as described in Ref.<cit.>. 
The mode effective charge is giant for the unstable mode in comparison to the others<cit.>.
 In fact, the very anomalous Born effective charges on Ti and O_∥ 
combined through the specific pattern of eigendiplacements associated to the TO_1 mode in order to produce a large spontaneous polarization, while for the TO_2 and TO_3 modes, 
the motions are coupled so that Ti and O generate polarizations that partly compensate.



  §.§.§ CaTiO_3
 


CaTiO_3 displays two main non-polar instabilities at the R and M points of the cubic Brillouin zone
related to antiferrodistortive (AFD) motions (ω≃ 209i cm^-1 and ω≃ 198i cm^-1, respectively). 
These correspond to cooperative rotations of oxygen octahedron, around the B-atoms, with consecutive octahedron along the rotation axis being in-phase at M 
(a^0a^0c^+ in Glazer's notation <cit.>) and anti-phase at R (a^0a^0c^- ).
As such, the M-instability appearing in the phonon spectrum is a continuation of the instability at R, 
while in BaTiO_3 it is a continuation of the polar instability at Γ. A detailed description of such AFD-instabilities is reported in Ref.<cit.>.

In addition, there is an unstable mode at the Γ-point that is also polar. 
This FE instability is now restricted to a region around the Γ-point highlighting a larger and more isotropic ferroelectric correlation volume.
This mode enables the condensation of a polar distortion. 
However, the character of the corresponding ferroelectricity is rather different 
compared to that in BaTiO_3 where it is dominated by the Ti displacements (B-site ferroelectricity). 
In fact, in CaTiO_3, this instability turns out to be dominated by the Ca displacements (A-type ferroelectricity) as can be seen from the eigendisplacements displayed in  
Table <ref> and the red color of the corresponding phonon in Fig. <ref>(b).
Despite this important involvement of the A-cation in the polar distortion, 
its Born effective charge is not strongly anomalous. 
This has to be ascribed to fact that the polar distortion of Ca is driven by a steric effect. 
In fact, we already mentioned that CaTiO_3 has a tolerance factor less than one. Therefore, the small size of the Ca ionic radius allows its distortion in the cubic perovskite structure.
Consequently, the involvement of the oxygens in the distortion is inverted with respect to the BaTiO_3 case, as 
the O_2/3 lying on the plane perpendicular to the direction of the distortion are now more involved than the apical oxygen (see Table <ref>).

These results can be rationalized by looking at the effect of the substitution of Ba by Ca on the on-site and interatomic force constants.
In fact, the “on-site" force constant of Ca as well as the A-A' interaction are 
significantly smaller than the corresponding ones in BaTiO_3, while the “on-site" force constant of Ti is increased.
Additionally, the destabilizing A-O_1 interaction
 (xx', Table <ref>) become significantly positive in agreement with the opposite direction
of the respective atomic eigendisplacements (Table <ref>), whereas along the Ti-O_1 chain a stronger repulsive interaction prevents an important participation of
the B-cation to the polar distortion. 
 
These previous observations can be also related to the concurrent appearance of AFD distortions in the system.
Although the IFCs between the reported pairs of oxygen are remarkably similar in BaTiO_3 and CaTiO_3 (and also in BaZrO_3 and CaZrO_3 as we will see below), 
the oxygens tilting is favored by the fast decrease of the oxygen “on-site” forces (linked to the destabilizing A-O_1 interaction) 
in the directions perpendicular to the B-O chains and to the increasing stiffness
in the parallel directioni (Table <ref>).
Moreover, the phonon dispersion curves appear
substantially flat along R-M suggesting the absence of coupling between the oxygens in different planes, but, as can be seen in Tables <ref>
and <ref>, the transverse interactions are far to be negligible. As proposed in Ref.<cit.>, this could be due to the joint action
of the A-O coupling with a compensation between different interplane interaction of the oxygens (see also Ref.<cit.>). 



  §.§.§ CaZrO_3
 


CaZrO_3 exhibits much more intricate phonon branches than CaTiO_3.
However, the dynamical properties of both compounds show some similarities. 
As in CaTiO_3, the strongest instabilities of CaZrO_3 are at the R and M points of the Brillouin zone 
and associated to AFD oxygens rotations (ω≃ 238i cm^-1 and ω≃ 233i cm^-1, respectively). Also, the dispersion curve along 
 R-M appears completely flat. 

Nevertheless, the lowest polar instability is no more confined at Γ, but the polar-antipolar instability extends all over the cubic Brillouin zone. 
This reflects a FE instability very localized in real space.

In CaZrO_3, the distortion is fully driven by the A-cation (Ca) and the O_2/3 anions. 
In fact, Zr on the B-site does not participate to the polar distortion with also subsequent reduced participation of the apical O_1 (see Table <ref>-(TO_1 mode)). 
Finally, the dynamical effective charges of both Zr and the related O appear less anomalous than in the titanates resulting in 
a lower value of the mode effective charge. 

As before, these dynamical properties can be understood in terms of IFCs. In fact, the “on-site” force associated to Ca and the in-plane “on-site” force of O are drastically reduced 
while for Zr it is strengthened by making the A-cation and the O_2/3 anions highly unstable and the B-cation almost fixed with respect to their high-symmetry positions. 
Accordingly, the Ca-O_1 interatomic force constants are largely dominated by the dipole-dipole interaction, 
while the Zr-O_1 interaction is strongly repulsive. 
The interplay between all these features allows the understanding of the strong A-site driven character of the instabilities.
In fact individual Ca displacements are nearly unstable (Table <ref>).

A proper comparison of CaZrO_3 an CaTiO_3 with respect to PbZrO_3 and PbTiO_3 respectively, can be done via the lattice dynamics analysis presented in Ref.<cit.>.



  §.§.§ BaZrO_3



Starting from the analyzed dynamics in BaTiO_3 and CaZrO_3, the dynamical properties of BaZrO_3 are fairly predictable. 
The substitution of the Ti atoms with Zr on the B-site joined to
the presence of Ba on the A-site leads to the stability of the cubic phase, also confirmed by steric arguments, since the tolerance factor is close to 1.
In fact, unlike the previous perovskite systems, the phonon spectrum of BaZrO_3 shows only a very weak instability at the R-point 
with associated phonon frequency ω≃ 60.2i cm^-1 and no unstable polar modes. 

By looking at the specific quantities reported in Table <ref>, 
the softest polar mode (TO_1) displays the smallest polarity (Z̅^*) with respect to the corresponding modes in the other compounds. 
Moreover, unlike Ti in BaTiO_3, a tendency of Zr to decrease the spontaneous polarization results from the specific combination of the associated pattern of distortion with the dynamical charges. 
Conversely, the second stable polar mode (TO_2), even if stiffer, displays bigger polarity because of the additive contribution of Zr and O. 

This behaviour can be justified by means of the force constants. In fact, the B-cation experiences an increased “on-site” term as well as O in its transversal direction
with respect to the case of BaTiO_3.
Accordingly, the Zr-O_1 is strongly dominated by the short-range forces.
On the other hand, the A-cation shows an “on-site” force constant in between BaTiO_3 and CaTiO_3, 
but the A-O interaction is still too weak to destabilize Ba (see Tables <ref> and <ref>). 


For all the four compounds, the interatomic force constants associated to the A-B cations and to the O-O interactions have the same sign and
 order of magnitude. This trend points out the key role played by the on-site force constants in generating different kinds of instabilities according
to the environment that each atom experiences in the different perovskites. In fact, changes in the
A-O and B-O interactions (Tables <ref> and <ref>) are primarily associated to changes of that values.








 §.§ Energetics of metastable phases






In the previous section we have seen that the parent compounds can develop two types of instability: a FE instability associated to a polar mode at Γ and 
AFD instabilities associated to unstable modes at M and R.
Now we quantify the energy difference associated to the condensation of the corresponding (polar and non-polar) distortions with respect to the cubic phase.
The results are shown in Fig. <ref>.
A careful crystallographic explanation of the symmetry changes associated to these distortions can be found in Ref.<cit.>.



  §.§.§ BaTiO_3

In BaTiO_3 the energy sequence of the polar phases arising from the condensation of the polar mode is in line with the experimental sequence of phase transitions when lowering 
the temperature (Pm3̅m → P4mm → Amm2 → R3m). At the same time, we do not find any other distorted structure as a possible metastable phase. 
To be noticed is that the energy landscape is relatively flat. The R3m ground-state is only 20 meV/f.u. below the paraelectric Pm3̅m phase and 
the energetic difference between the polar phases is ≃ 5 meV/f.u.



  §.§.§ CaTiO_3

Interestingly, in CaTiO_3 the energy sequence of the polar phases, arising from the condensation of the unstable polar mode, 
is reversed. In addition, the relative energy difference is larger than in BaTiO_3.
In fact, the tetragonal P4mm structure is the lowest energy configuration among the polar ones with a relative energy gain of about 70 meV/f.u.  
Moreover, additional non-polar structures arising from the condensation of the unstable AFD modes at M and R 
appear as metastable or eventually stable. 
As usual in perovskites with t<1 <cit.>, in fact, the ground-state corresponds to the Pnma (a^-b^+a^-) structure which is ≃ -350 meV per f.u. below the  reference cubic phase. 
One-tilt structures, such as P4/mbm (a^0a^0c^+) or I4/mcm (a^0a^0c^-), 
and the two-tilts Imma (a^0b^-b^-) are intermediate metastable structures that appear ≃ 100 meV and ≃ 50 meV per f.u. above the ground-state, respectively.



  §.§.§ CaZrO_3

In CaZrO_3 the energy sequence of both the polar and AFD phases is the same as CaTiO_3, but they are shifted down in energy
becoming very spread. The AFD-Pnma ground-state is lower by about -1 eV/f.u. than the Pm3̅m phase. 
This is in tune with the very high transition temperature experimentally observed for the sequence Pm3̅m → Pnma, even if the phonon frequencies related to the unstable modes are 
close to the ones in CaTiO_3. 



  §.§.§ BaZrO_3

For BaZrO_3 
we found an interesting energy competion between the antiferroditortive R3̅c (a^-a^-a^-), Imma (a^0b^-b^-) and I4/mcm (a^0a^0c^-) structures. 
According to the very tiny value of the instability (Fig. <ref>(c)), 
the condensation of the oxygens rotations provides an energy gain relative to the cubic phase of about ≃ 2.5 meV per f.u., 
while the three distorted phases have an energy that differs by less than 0.1 meV/f.u. 
In spite of the negligible energy gain, the amplitude of the AFD distortion is significant. 
The biggest distortion is provided by the tetragonal phase with an angle of rotation of ∼ 4^∘ about the [001] direction. 
 



  §.§.§ Ferroelectric phases & polar modes 

In order to interlock the optimized polar structures with the lattice dynamics, we evaluated the contribution 
of each polar mode to the condensed distortion. 
The overlap matrix is reported in Table <ref>. It is interesting to notice, that the three ferroelectric states are mostly due to the condensation of the unstable optical mode
 for all the three perovskites BaTiO_3, CaTiO_3 and CaZrO_3, so that it is possible to establish a nearly one-to-one correspondence with the 
pattern of distortion associated to the unstable TO_1 mode and the displacements as obtained from the structural optimization,
 while the contribution of the TO_2 and TO_3 modes remains very small.

An important remark is the huge difference in the total distortion τ between BaTiO_3 and the Ca-based perovskites, that allow for possible bigger spontaneous polarization in the latter 
compounds even if the A-cation at play is not either stereochemically active or involved in the hybridization with the oxygens, like Ti in BaTiO_3. 





The previous analysis of the dynamics and energetics associated to the four parent compounds 
has emphasized some similarities and differences. This is helpful for a better understanding 
of properties arising while mixing cations at the A- and B-sites in (Ba,Ca)TiO_3 and Ba(Ti,Zr)O_3 solid solutions, respectively.
We will now present results coming from a systematic characterization of the latter systems by also testing and comparing two different approaches:
 the “virtual crystal approximation” (VCA) and supercell-based calculations.
 
We note that the parent BaTiO_3 and CaTiO_3 compounds display inverted sequence of polar phases. Then, for the (Ba,Ca)TiO_3 solid solution, 
it can be expected the emergence of a region with strong competition between these phases and a crossing point in energetics. This will be confirmed in the next Section.



§ SOLID SOLUTIONS: (BA,CA)TIO_3



In the last few years, Ba_1-xCa_xTiO_3 (BCT) has started to arouse curiosity in the experimental community as the Ca off-centering seems to play an important role in stabilizing 
ferroelectricity against chemical pressure effects <cit.>. 
In particular, in Ref.<cit.>, Fu and Itoh have characterized the single crystals of BCT in a temperature range from 2 K to 400 K and for compositions ranging from 
x=0.00 up to x=0.34. They found that the Curie point 
is nearly independent of the Ca-concentration for the Pm3̅m → P4mm transition, whereas there is a shift of the P4mm → Amm2 and Amm2 → R3m phase transitions 
toward lower temperatures. Accordingly, the consequence of Ca-substitution is the stabilization of the tetragonal ferroelectric phase. 

Let us now analyze the dynamics and energetics as predicted by means of first-principles calculations. We first report results from the VCA approach in Sec. <ref>, 
then from the use of supercells in Sec. <ref>.





 §.§ Virtual Crystal Approximation (VCA)





  §.§.§ Lattice parameter


First, we report the evolution of the lattice parameter of the Pm3̅m-cubic phase as obtained from the structural optimization within VCA. 
Because of the reduced volume of CaTiO_3 with respect to BaTiO_3 (Table <ref>), the lattice parameter decreases monotonically with the Ca-concentration, 
but the trend deviates from the linearity of the Vegard's law, as shown in Fig. <ref>(a). 

As for the pure compounds, for solid solutions we can define the tolerance factor as t= (r_A+r_O)/√(2)(r_B+r_O), with r_A=(1-x)r_Ba+(x)r_Ca. 
Values of the ionic radii, r_i, for the pure atoms have been taken from Ref.<cit.>. For increasing x, t decreases from 1.06 to 0.97 
reaching 1 at x=0.6, as reported on top of Fig. <ref>.



  §.§.§ Dynamical properties




Then, we analyze the evolution of the lowest phonon frequencies at the high symmetry points of the Brillouin zone of the cubic reference. 
DFPT calculations have been performed on the optimized structures.
We find that, for increasing Ca-concentration, the original polar instabilities of BaTiO_3 
changes character by becoming confined to the center of the Brillouin zone (AFE instabilities at X and M disappear), while the AFD modes
become unstable, as shown Fig. <ref>(b). Changes in the phonon instabilities are linked to the evolution of t. 
In fact, the polar instability dominates until t becomes smaller than one and the AFD modes become largely unstable. 
Specifically, 
the lattice dynamics is BaTiO_3-like for 0 ≤ x ≤ 0.2. The unstable polar mode at Γ is mainly sustained by the destabilizing 
Ti-O interaction, while the A-O interaction stays mostly repulsive by resulting in the inactivity of the A-site, as reproduced in Fig. <ref>(c,d,e). 
Accordingly, the instabilities at the X- and M-points are also due to Ti-O polar motion. For 0.3 ≤ x ≤ 0.6, the scenario starts to change. 
The polar instability becomes weaker as it results from the important reduction in the amplitude of the associated frequencies. 
Particularly, the polar distortion remains unstable at Γ, while it becomes progressively stable at X and M. 
These smooth changes can be related to smooth changes in the type of interaction between the cations and the oxygens. A change of sign is observed along 
both the Ti-O and A-O interaction that corresponds to a strong competition between attractive and repulsive forces in determing the nature and character of the polar instability: 
a reduction of the amplitude of the Ti-motion corresponds to increasing negative value of the Ti-O interaction, while an increasing contribution to the polar distortion of the A-site 
arises from positive values of the A-O interaction, as shown in Fig. <ref>(c,d,e). 
Finally, for 0.7 ≤ x ≤ 1.0, the lattice dynamics becomes CaTiO_3-like. Accordingly, instabilities related to oxygens rotations appear at the M- and R-points of 
the cubic Brillouin zone with larger associated frequencies than the polar one, that remains unstable only at Γ. This scenario results from strong destabilizing A-O interactions and 
largely repulsive Ti-O interactions. Accordingly, the character of the polar distortion also changes by becoming largely driven by the motion of the A-site with respect to the B-one.

Moreover, changes in the phonon behaviour also affect the dynamics of the oxygens. 
The increasing contribution to the distortion from the A-site, i.e. increasing long-range forces between the A-cation and oxygens,
 favors the motion of the planar oxygens with respect to the apical one (labelled 
O_2/3 and O_1 in Fig. <ref>(e), respectively). Specifically, the contribution is reversed when going from t>1 to t<1. 


The latter analysis of the evolution of the dynamical properties when going from the BaTiO_3-rich region to the CaTiO_3-rich one in the Ba_(1-x)Ca_(x)TiO_3 “virtual-system" 
reveals the presence of two parallel mechanisms: the progressive weakening of the long-range forces between Ti and O atoms in favor of their strengthing between the A-site and oxygens. 
In terms of character of the phonons instabilities, this change corresponds to a smooth evolution from the B-driven into the A-driven distortions
associated to the parents BaTiO_3 and CaTiO_3, respectively. 
The reasons behind these two phenomena are also different and complementary. In fact, Ca-doping of the virtual A-site produces: 
(i) varying interatomic force constants between the A-O atomic pair in favors of destabilizing long-range forces sustained by the progressive lowering of the A-cation stiffness 
when going from BaTiO_3 to CaTiO_3 (see Table <ref>); 
(ii) reduction of the volume, that can be considered as increasing isotropic pressure on the ATiO_3 system and, therefore, shortnening of the Ti-O bond lengths. 
The latter effect produces increasing stiffness of Ti-atoms, that results in the change of sign of the Ti-O interatomic force constants and associated weakening of the B-type ferroelectricity.
In fact, as investigated at the first-principle level in Refs. <cit.> and as it turns out when comparing graphs (a) and (c) reported in Fig. <ref>, 
the balance between short-range and long-range forces between the Ti and O atoms is strongly sensitive to pressure. 
At variant, the varying composition does not affect significantly the A-B interaction, that in fact remains almost the same as in the pure parent compounds 
(see Tables <ref> and <ref>).



  §.§.§ Energetics landscape




As we previously did for the pure parent compounds, now that the main instabilities of BCT “virtual-system" have been identified,
 we look at the energy competition between different phases arising from the condensation of the corresponding modes in order to obtain an overview of the energetics 
trend as a function of varying concentration. 
As shown in Fig. <ref>, the main effect of the Ca-substitution 
is to reduce the energy gain both between the three ferroelectric phases and with respect to the cubic phase by making them strongly competitive, while 
additional antiferrodistortive structures appear as metastable in the CaTiO_3-rich region by providing a non-polar ground-state.

Such behaviour proceeds from the evolution of the dynamical properties previously analyzed. In fact, an important energy competition between polar phases starts at x=0.3, that corresponds 
to the concentration where we observe the change of sign of the Ti-O IFCs (see Fig. <ref>(c)). 
Then, the energies of the polar phases smoothly converge to that of the cubic one in line with the progressive weakening of the polar instability. 
Noteworthy is that the inversion in the energy sequence of the three polar phases happens at x=0.9 where we 
observe the inversion of the dominant character of the polar instability between virtual A-cation and Ti (see Fig. <ref>(e)).
 
Additional structures with oxygens rotations appear with an higher energy gain from x=0.7 proceeding from the appearence of the AFD phonons instabilities 
(see Figs. <ref>(b) and <ref>). 
It is also noteworthy that the contribution to the polar distortion arising from the apical and planar oxygens is inverted from this point, 
with the dominant motions of the latter ones (see Fig. <ref>(e)).



  §.§.§ Polarization and piezoelectric response




Proceeding from the previous observation that for 0.0 ≤ x ≤ 0.6 the BCT “virtual-system" experiences 
a progressive weakening of the ferroelectric instability 
and that no other phonons are unstable, we now evaluate the evolution of the spontaneous polariation, P_s, and piezoelectric coefficients, d_ij. 
According to the VCA results,
we perfomed calculations on the lowest energy rhombohedral-R3m structure (Fig <ref>). Values are normalized to the calculated ones of pure R3m-BaTiO_3, 
that are P_s^BTO≃ 38 μ C/cm^2 and d_11^BTO≃ 76  pC/N, d_15^BTO≃ 270  pC/N and d_33^BTO≃ 15  pC/N. 

Specifically, in Fig. <ref>, we report P_s as obtained both via first-principles calculations by means 
of the Berry Phase Theory <cit.> (blue closed circles) and via the Born effective charges by means of the approximation 
P_s,α=1/Ω∑_k,βZ^*_k,αβΔτ_k,β <cit.> (green open squares).
It results that P_s decreases as function of increasing Ca-composition in tune with the decreasing energy gain associated to the polar phases in this range of composition, 
i.e. the weakening of the polar instability detected within VCA (Figs. <ref> and <ref>). 
In order to clarify if changes in the polarization are mostly due to varying effective charges or atomic distortions, we disentangled the two contributions.
In the first case, we kept constant the atomic displacements of pure R3m-BaTiO_3 and took into account the different effective charges associated to each BCT compositions (orange open rhombus). 
In the second case, we kept constant the Born effective charges of pure R3m-BaTiO_3 and considered the distortion arising from different compositions (red open circles). 
It results that the evolution of P_s has to be widely ascribed to decreasing atomic distortion in BCT “virtual-system" in such range of composition, as shown in Fig. <ref> (top). 
In fact, this behaviour is in line with the evolution of the eigendisplacements associated to the unstable polar mode at Γ 
with a pronounced reduction of the titanium and oxygens motion Fig. <ref>(e).

Despite the decrease of polarization, the piezoelectric response increases with the Ca-concentration (Fig. <ref> bottom). 
This trend is due to the fact that the free-stress dielectric response, ε_ij, remarkably increases as well, as shown in the inset of Fig. <ref>. 
By writing down the piezoelectric coefficient as d_ik=de_s,i/dE_k∝(χ_j,k P_s,l+P_s,jχ_l,k) <cit.>, it is easier to understand 
such behaviour: near a phase transition, where the lowest-frequency polar mode, i.e. the soft mode, goes to zero, the dielectric response diverges <cit.>. 
Accordingly, the calculated values of the lowest-frequency polar modes in the R3m-phase evolve like 169, 160, 124, 73 (in cm^-1)
for 0.0, 0.2, 0.4 and 0.6 x-composition, respectively.
Therefore, the softening of the polar mode overcomes the reduction of P_s in the piezoelectric response within the VCA approach.
 



 §.§ Supercell Calculation



In order to check the validity of VCA and to better characterize the impact of the (Ba,Ca) substitution on the dynamics of the system
we performed direct DFT calculations on the Ba_0.875Ca_0.125TiO_3 and Ba_0.50Ca_0.50TiO_3 compositions by means of supercells.
Details about the atomic arrangements have been provided in Sec. <ref>.



  §.§.§ Ba_0.875Ca_0.125TiO_3




Phonons calculations performed at the Γ-point of the cubic Pm3̅m-supercell of Ba_0.875Ca_0.125TiO_3 (Fig. <ref>(a)) have revealed several instabilities 
related to polar modes. The most unstable one, with ω≃ 169i cm^-1, is associated to a polar distortion driven by all the Ti-atoms and the single Ca-atom 
against the oxygens. Ba-atoms are almost fixed as in the parent BaTiO_3. 
Therefore, Ca is polar active on the A-site already at low concentration.
The different dynamics can be analyzed in terms of interatomic force constants between the A-cations and oxygens in the investigated solid solution.
In fact, the A-O interaction results to be opposite if we focus on Ba or Ca atoms in line with the dynamical properties of the respective pure compounds (Table <ref>). 
As reported in Table <ref>, the Ba-O interaction remains dominated by repulsive forces, while the Ca-O interaction is largely dominated by the destabilizing long-range interaction. 
Moreover, by looking at the atomic pair properties reported in Table <ref>, a further effect of the Ca-presence is to make oxygens 
inequivalent even in the high centro-symmetric reference. This effect is due to an internal 
structural relaxation around the Ca-atom because of its smaller size than Ba. This mechanism strongly affects the Ti-O interaction. In fact, the Ti atom experiences 
shorter and longer bonds with the corresponding apical oxygen, that produce strong competition between long- and short-range forces: along the longer bond the dipole-dipole interaction 
overcomes the repulsion, while it is the opposite along the shorter bond (see values for the supercell case in Table <ref>). 
For this specific BaTiO_3-rich concentration, such delicate balance 
between the attractive and repulsive interactions along the Ti-O bond results in the weakening of the B-site driven ferroelectricity.

The intriguing manifestation of this phenomenon is the achievement of a quasi-degeneracy between different polar states.
In fact, we analyzed five different polar states corresponding to the tetragonal, orthorhombic, rhombohedral, monoclinic and triclinic phases with P4mm, Amm2, R3m, Pm and P1 
space-group symmetry, respectively. 
Specifically, the everage gain in energy with respect to the Pm3̅m-cubic phase is about -20 meV/f.u. as for the parent BaTiO_3. Nevertheless, the maximum energy difference between 
the different polar phases is of about 0.5 meV/f.u.

According to the latter isotropic energy landscape, the spontaneous polarization results to be essentially the same for all these ferroelectric phases, that is 
P_s≃ 39.5 μ C/cm^2 as obtained via Berry phase calculations. This value is slightly higher than P_s of pure BaTiO_3 in the R3m-phase showed in Fig. <ref> of 
the previous section. 
Moreover, we also calculated the d_ij coefficients of the piezoelectric tensor in the stable R3m-phase. The d_33 component, parallel to the polarization direction, remains 
unchanged with respect to BaTiO_3, while the d_11 and d_15 components, tranversal to the polar axis, are considerably enhanced. 
Specifically they are: d_33≃ 15, d_11≃ 344 and d_15≃ 1458 in pC/N. 





  §.§.§ Ba_0.5Ca_0.5TiO_3


First result in Ba_0.5Ca_0.5TiO_3 ((Fig. <ref>(b-d))) is the preference of an ordered configuration composed by alternating layers of the same type of A-cations 
with respect to a columnar and rocksalt ordering. 
The relative energy gain of the high-symmetry columnar and layered configurations with respect to the rocksalt one is about -35.5 meV/f.u. and -81.4 meV/f.u., respectively.



Then, as previously, we performed DFPT calculations on the three high-symmetry reference supercells.
All three configurations display several instabilities.
In all of them the most unstable mode corresponds to a polar distortion driven by both the Ti and Ca atoms like in Ba_0.875Ca_0.125TiO_3, but 
with an average major contribution from calcium. The phonon frequency associated to that mode 
is ω≃ 141i cm^-1 in all three cases. 
Additionally, AFD modes related to oxygens rotation become unstable. However, the energy gain carried to the system by the condensation of the latter distorions alone is lower 
than the one coming from the polar distortions, that determine the ground-state. 

Accordingly, the most stable polar phase results to be the P2mm-structure with the polarization along the [100] direction of the layered atomic arrangement. 
The associated P_s is of about 51 μ C/cm^2 and it is considerably greater than the one calculated in P4mm-BaTiO_3, that is P_s≃ 34 μ C/cm^2.
Moreover, in order to see eventual effects of the atomic configuration on the energy competition between different polar states, 
we condensed different patterns of polar distortions in the three arrangements. The associated energy landscape is quite spread as reported in Fig. <ref>(c). 
It results that, for x=0.50, the system prefers a state with only “one-component” polarization independently of the atomic configuration. 
The R3m-like state with polarization along the [111] pseudo-cubic direction is largely penalized.

Also in this case, changes in the dynamical properties (i.e. substantial contribution from Ca-atoms to the polar distortion and the appearance of unstable AFD modes) and in the 
energetics of the polar phases can be analyzed in term of interatomic force constants. As made clear in Tables <ref> and <ref>, 
the main features are the same as the ones described for the x=0.125 Ca-concentration. 
However, the higher concentration of calcium produces here: (i) a global weakening of the long-range interaction between Ti and O atoms in favor of repulsive forces via 
further shortened Ti-O bonds lenghts, i.e. further reduction of the [TiO_6]-octahedral volume; this can be the reason behind the tendency of the system to prefer the 
“tetragonal-like” polar state; (ii) the presence of unstable AFD modes because of the reduced number compressed Ba-O distances otherwise preventing the octahedral rotations.



 §.§ VCA vs SUPERCELL calculations


By comparing results from supercells-based calculations with the ones from VCA,
it seems that the latter approach is successful to provide an average and qualitative identification of the main changes occuring in the (Ba,Ca)TiO_3 solid solution. 
However, it cannot provide a proper description of the microscopic mechanism behind properties. 
The reason of such a limitation is that, by construction, VCA considers all A atoms as equivalent and forces them to behave similarly while, in fact, 
they want to adopt different behaviors.
The importance of considering the actual local cations is clear by comparing the quantities reported in Tables <ref> and <ref>, where we clearly see different results are obtained depending on the method. 
In fact, as deeply discussed above, evolutions in the dynamics as well as in the energy landscape are strongly related to the presence of calcium that 
directly interacts with the oxygens and, indirectly, affects the Ti-O interaction via local changes of the structure.
Therefore, the progressive weakening of ferroelectricity, within VCA, is not representative of the real behaviour. 
In fact, the actual weakening of the Ti-driven ferroelectricity due to decreasing volume is progressively compensated by the Ca-driven ferroelectricity. 
Accordingly, the inversion in the sequence of polar phases when going from BaTiO_3 to CaTiO_3 appears already at x=0.50 within the supercell, while it is expected for x≥0.9 within VCA.
Additionally, an opposite trend of the spontaneous polarization is also obtained within the two approaches: decreasing values for larger x within VCA, while increasing ones within the supercells.

As a global result, the success of VCA in detecting the weakening of the B-driven ferrolectricity is due to the fact that the Ca-doping mainly induces a steric effect 
because of the volume contraction when the Ca-concentration increases, while the electronic properties are not primary. 
In fact, the top of the valence is mostly characterized 
by the O p and Ti d states entering the hybridization mechanism, while the s and p states of Ba and Ca stay well below the Fermi level. 
Nevertheless, 
the trends and microscopic mechanisms involved as predicted by means of supercells calculations are in line with experimetal observations<cit.>. 
This remark confirms that, in order to provide a proper first-principles characterization of (Ba,Ca)TiO_3, VCA-based approaches 
can provide some trend, but are not appropriate to explain the underlying physics; for the latter, 
supercells-based calculations taking explicitly into account the different nature of the cations are more suitable. 



§ SOLID SOLUTION: BA(TI,ZR)O_3


Ba(Ti,Zr)O_3 (BTZ) solid solutions involve the homovalent substitution between Ti^4+ and Zr^4+, that are nevertheless quite different atoms both for the ionic radii 
(0.605 Å and 0.72 Å respectively <cit.>) and the electronic configuration (Ti 4s^2 3d^2 and Zr 4d^2 5s^2) entering the hybridization mechanism, that for instance
leads the ferroelectricity in BaTiO_3. 

The first experimental investigation of the phase diagram of the BaTiO_3-BaZrO_3 binary system dates back to the 1956 with the work of Kell and Hellicar<cit.>, 
reporting the abrupt effect of Zr-concentration on the decreasing of the Curie point for the three ferroelectric phases of BaTiO_3. Subsequently, 
other solid experimental investigations came out in the late '90s, 
when J. Ravez et al. in Refs.<cit.> provided a
clear distinction of phases in the BaTi_1-yZr_yO_3 ceramics diagram: classical ferroelectric BaTiO_3-like for 0.00≤ y≤0.10,
only one ferroelectric-paralectric transition observed in the range 0.10≤ y≤0.27 and relaxor ferroelectric behavior for 0.27< y≤0.42.
Then, in 2004, A. Simon et al. also
investigated the crossover from a ferroelectric to a relaxor state in lead-free solid solutions in Ref.<cit.>, confirming that beyond a definite concentration y,
BTZ ceramics show relaxor properties. For BaTi_0.80Zr_0.20O_3 only one resonance of the permittivity at the ferroelectric-paraelectric T_C of about 315 K
 is observed, whereas for y=0.35 the dielectric anomaly is broad and frequency-dependent as a function of temperature. Moreover in Ref.<cit.>, they provide
an EXAFS study of BTZ systems and conclude that BTZ is relaxor and the avarage crystal structure is cubic (Pm3̅m space group)
in the range 0.25≤ y≤0.50. In addition, in a complementary work on BCTZ<cit.>, 
the X-ray diffraction has revealed (110) peaks in BaTi_0.80Zr_0.20O_3 ceramic but a weak tetragonality (i.e. a/c∼ 1), that means closeness to the cubic phase. 
However, no experimental data on single crystal samples are available for direct comparison with ab-initio results.


As for the previous case of (Ba,Ca)TiO_3, we first investigate the Ba(Ti,Zr)O_3 system first by means of the virtual crystal approximation. 
Then, we go beyond by using supercells-based calculations.



 §.§ Virtual Crystal Approximation (VCA)





  §.§.§ Lattice parameters


We first report the trend of the lattice parameter of the Pm3̅m-cubic phase as a function of Zr-composition. Because of the 
bigger volume of BaZrO_3 with respect to BaTiO_3 (Table <ref>), 
the trend is monotonically increasing with the amount of zirconium and the agreement with the linearity of the Vegard's law is quite 
satisfactory (Fig. <ref>(a)).
Contrariwise, because of the bigger ionic radius of Zr than Ti, the tolerance factor[also in this case, 
we defined an average radius for the B-cation as r_B=(1-y)r_Ti+(y)r_Zr] decreases without going below 1, that is the value for pure BaZrO_3, 
as reported on top of Fig. <ref>. 



  §.§.§ Dynamical properties




Then, we focus on the dynamics by 
looking at the evolution of the lowest phonon frequencies at the high-symmetry points of the cubic Brillouin zone.
Calculations have been performed on the VCA-optimized Pm3̅m structure. 
It results that the main effect of the Zr-doping is to abruptly stabilize the cubic phase, as clear from Fig. <ref>(b). 
In fact, already in the Ti-rich region, i.e. y>0.10, the original polar instability vanishes and an AFD unstable mode BaZrO_3-like appears at the R-point only for y≥0.90.
The reason behind such abrupt weakening of ferroelectricity in the BTZ “virtual-system" at low y-concentration 
can be traced back to abrupt changes in the type of interaction between the virtual cation at the B-site and 
the oxygens. In fact, the total IFC along the B-O direction is of about -0.001 Ha/Bohr^2 at y≥0.05. Therefore, the 
 change of sign with respect to the corresponding interaction in pure BaTiO_3 (Tables <ref>) means that the Zr-substitution strongly favors 
short-range repulsive forces.

Differently from BCT systems, it is not possible to relate such behaviour to the evolution of the tolerance factor or to simple volumetric reasons. 
In fact, in line with the properties of the two parent compounds, t remains larger than one and the volume increases for each intermediate composition. 
The main mechanism involved here is the weakening of the Ti-O interaction when introducing Zr on the same B-site, as this corresponds to either 
direct changes in the O 2p and B d hydridization mechanism from the electronic perspective or
breaking correlated Ti-O chains necessary to sustain ferroelectricity in BaTiO_3 from the lattice dynamics. In fact, we have already seen in Sec.<ref>, that the Zr-O interaction in 
BaZrO_3 is strongly repulsive. Moreover, differently than the role played by CaTiO_3 in BCT, no contribution from the A-site is observed, 
as the Ba-O interaction is too weak to sustain alone a polar instability in BTZ.

As a consequence, it is finally not surprising that the “virtual”-BTZ system does not display any instabilities for a very large range of composition.



  §.§.§ Energetics landscape
 



In line with the drastic changes of the dynamical properties in BTZ predicted by VCA, 
the energy difference both between the three polar phases and with respect to the cubic one is abruptly reduced, as shown in Fig. <ref>.
In fact, the energy competition is within 0.1 meV/f.u. at y=0.10 and beyond this concentration no stable or 
metastable polar phases are allowed within VCA. 



  §.§.§ Polarization and piezoelectric response




After identifying the abrupt weakening of ferroelectricity in the BTZ “virtual-system”, we now evaluate the associated effect on the piezoelectric response. 
Therefore, we calculate the spontaneous polarization P_s and the piezoelectric coefficients d_ij for the 
stable phase, that is the R3m-phase for y=0.05 and 0.10 (Fig.<ref>). Also in this case, as we are mostly interested in the qualitative trend of these 
properties as a function of Zr-doping, we report values normalized to those of R3m-BaTiO_3. 

In Fig. <ref> we report values of P_s calculated both via Berry phase (blue closed circles) and 
Born effective charges (green open squares). This is to be compared with the same calculations for BCT in Fig <ref>. 
In line with the strong weakening of the polar instability, P_s is drastically reduced by the Zr-doping. 
Similarly to BCT, from the futher analysis of the distinguished contribution from effective charge variation and atomic displacements, it results that the 
evolution of P_s has to be widely ascribed to largely reduced atomic motions. However, 
by looking at y=0.10, it is important to note that there is a disagreement between the value of polarization 
by fixing the Born effective charges to that of pure R3m-BaTiO_3 (red open circles) 
or to the ones of R3m-BTZ “virtual-system” (green open squares) as 
the Zr-doping also affects the Born effective charges. In fact, we have already shown that 
the Z^* associated to Zr is less anomalous than the one of Ti (Table <ref>). 
Therefore, the average effective charge associated to the “virtual” B-cation is also reduced with increasing Zr-concentration. 
This is an important warning about the ability of VCA in describing the properties arising from the mixing of Zr and Ti on the B-site, as they have quite 
different electronic properties (see Sec.<ref> below).

Such P_s, actually, favors piezoelectricity. In fact, the huge increase of dielectric response, ε_ij, 
associated to the drastic softening of the polar mode,
largely compensates the decrease of P_s. This results in the enhancement of the piezoelectric response with respect to pure BaTiO_3,  
as reported in Fig. <ref>. Accordingly, the values of the lowest-frequency polar mode in the stable R3m-phase evolve like 169, 146, 69 (in cm^-1) 
for y=0.00, 0.05 and 0.10, respectively.


Despite the overall predictions provided by VCA seem to be reliable when comparing the behaviour of the two parent compounds and 
with respect to experiments, supercells-based calculations are better required to clarify the local effect of Zr-doping on the dynamics of BaTiO_3.



 §.§ Supercell Calculations
 


The dynamical properties of BaTi_0.875Zr_0.125O_3 and BaTi_0.50Zr_0.50O_3 supercells are here investigated. Structural details are provided
 in Sec. <ref>





  §.§.§ BaTi_0.875Zr_0.125O_3




BaTi_0.875Zr_0.125O_3 supercell (Fig. <ref>(e))
hosts several instabilities in the  cubic Pm3̅m phase, all related to ferroelectric and antiferroelectric distortions.
The most unstable, with ω≃ 214i cm^-1, is a polar mode mostly related to the opposite motion of Ti and O, 
while the Zr atom slightly moves in the same direction as the oxygens. 
The difference in these relative displacements affects the neighboring atom lying on the polar direction, that 
slightly moves with the oxygens as well. 
That is, in terms of the labels indicated in Fig. <ref>(a), if the polar distorion is along the x direction, for instance, Ti_1 slightly off-centers with Zr 
in the same direction as the oxygens, but opposite to the other Ti-atoms. Ba-atoms are almost fixed as in the parent BaTiO_3. 
Such different dynamics is made clear when looking at the interatomic force constants reported in Table <ref>. 
In fact, the original destabilizing Ti-O interaction of BaTiO_3 is not preserved for each direction in the system. 
The Ti-atoms having Zr as first-neighbors experience a strong repulsive Ti-O interaction along the O-Zr-O-Ti-O chain direction, 
while the long-range forces remain the major ones 
along the other directions as well as for the other Ti-atoms. 
This is strongly ascribed to the “chain-like" character of the polar instability of BaTiO_3, as discussed in Sec. <ref>:
 the persistence of the polar instability requires a minimum number of correlated Ti-O displacements <cit.> that, however, are prevented along the chain containing Zr. 
In fact, according to the properties of BaZrO_3 described in Sec. <ref>, the Zr-O interaction remains strongly repulsive also in this Ti-rich solid solution. 
However, along the preserved O-Ti-O-Ti-O chains, the long-range interaction is even indirectly strenghted by the Zr-presence. In fact, the bigger size of Zr 
(i.e. of the [ZrO_6]-octahedral volume) makes the Ti-O bond's length perpendicular to the O-Zr-O-Ti-O chain direction longer than in pure BaTiO_3. This mechanism 
locally favors the long-range interaction. 
Accordingly, the “on-site" force constants tensor associated to the Ti-atoms adjacent to Zr are anisotropic: by referring to Ti1 of Fig. <ref>(a), 
the “on-site" terms are (+0.247, +0.132, +0.132) Ha/Bohr^2. Therefore, they are much more stiffer along the Zr-direction, while softer along the tranversal directions than in 
BaTiO_3 (see Table <ref>). The other Ti-atoms display a slightly anisotropic tensor, but similar to the one in BaTiO_3: by referring to Ti2 of Fig. <ref>(a), 
values are (+0.162, +0.151, +0.162) Ha/Bohr^2.
In contrast, Zr displays an isotropic tensor BaZrO_3-like with diagonal values of about +0.229 Ha/Bohr^2,
as it experiences the same interactions in each direction.
Concerning Ba at the A-site, the dynamics is the same as in pure BaTiO_3. In fact, the A-O interatomic force constant is close to zero along parallel direction of the coupling meaning 
that Ba distortion does not lead ferroelectricity in BaTi_0.875Zr_0.125O_3. The A-B interaction is also similar to the one in the two parent compounds 
reported in Table <ref>.

Such complex dynamics results in slightly reduced total polarization and more spread energy landscape with respect to the case of the parent BaTiO_3.
On one hand, the calculated spontaneous polarization is of about 28 μ C/cm^2 and 35 μ C/cm^2 for the P4mm and R3m phases, respectively, 
while it is of about 34 μ C/cm^2 and 38 μ C/cm^2 for the corresponding phases in pure BaTiO_3. This has to be assigned to the 
local depolarizing contribution arising from the O-Zr-O-Ti-O chain.
On the other hand, the energy gain of the three polar phases relative to the cubic phase is higher than in BaTiO_3 and it appears much more pronounced for the 
Amm2 and R3m phases as reported in Fig. <ref>(c). This effect traces back to the overall volume increase, 
that further favors the rhombohedral phase. 



  §.§.§ BaTi_0.50Zr_0.50O_3



