text:
# Entropy Analysis of Stochastic Processes at Finite Resolution

## Abstract

The time evolution of complex systems usually can be described through stochastic processes. These processes are measured at finite resolution, what necessarily reduces them to finite sequences of real numbers. In order to relate these data sets to realizations of the original stochastic processes (to any functions, indeed) it is obligatory to choose an interpolation space (for example, the space of band-limited functions). Clearly, this choice is crucial if the intent is to approximate optimally the original processes inside the interval of measurement. Here we argue that discrete wavelets are suitable to this end. The wavelet approximations of stochastic processes allow us to define an entropy measure for the order-disorder balance of evolution regimes of complex systems, where order is understood as confinement of energy in simple local modes. We calculate exact results for the fractional Brownian motion (fBm), with application to Kolmogorov K41 theory for fully developed turbulence.

Many physical systems investigated at present have a complex evolution in time. Frequently, the major information we can obtain on their dynamics comes from time series of noisy appearance consolini ; politzer ; santi ; sreenivasan ; gafarov . These series are samples at finite resolution of subjacent stochastic processes, whose properties are better investigated through a multi-resolution approach, because the realizations of the mentioned stochastic processes are singular everywhere functions. This is, for example, the case of $`1/f^\alpha `$ noises from self-organized systems bak ; durin ; urbach ; cote ; carreras . Thus, singularities are common, and should be interpreted as details that influence function variation at all scales muzy ; mallat . In this paper we propose a method for entropy analysis of arbitrary complex time series that takes advantage of this theoretical standpoint. Moreover, we account for the fact that measurements are made at finite resolution, considering the consequences of sampling, what is not fully accomplished in previous approaches. In general, other formulations are strongly influenced by information-theoretical arguments, applied to the analysis of chaotic behavior bandt . We see difficulties in two main aspects. First, the notion of complexity that is employed in these formulations is based on entropy rates, like the Kolmogorov-Sinai entropy, that measure the degree in which information on initial conditions is lost when the systems evolve. This gives a scale of complexity ranging from zero (non-chaotic deterministic case) to infinite (stochastic case). In this scale, finite values quantify the complexity of deterministic chaos. As a consequence, the problem of defining a proper complexity estimator for stochastic processes is substituted by the statistical investigation of chaotic deterministic dynamical systems eckmann . Second, these formulations do not consider the relationship that scaling has with disorder madalena . Our approach is tailored to face these difficulties. Here, we begin by assuming that the stochastic processes are supported on the real axis, and measured at a discrete set of points with a sampling interval $`\tau `$, during a time $`T`$. The resolution of the measurement is $`N=T/\tau `$, and for convenience we make $`T=1`$, and $`\tau =1/N=2^J`$ (we will work in these units). Such measurement results in a loss of information, which depends on two factors: The resolution, and the interpolation space in which the stochastic processes will be projected unser . The idea is that the mere sampled values say almost nothing about a function. Much more information is conveyed through the hypothesis on how the function varies between the sampled values. This regularity hypothesis, usually implicit when we draw smooth curves between data points, is an essential element of the theory. Without this hypothesis, no information found on the discrete and finite data sets can be attributed to the subjacent model, that is assumed to hold on the continuous support. Now, in the time-scale (time-frequency) domain, it makes sense to search a representation *i*) that is minimally affected by the created end singularities, *ii*) that provides an interpolation based on multi-scale approximation of the actual singularities, and *iii*) that deals with the resolution $`N`$ as a direct experimental parameter. The first and second requirements are the most crucial, and establish the way in which the energy ($`L^2`$-norm squared) is assumed to be distributed inside the interval of measurement, so that this distribution corresponds to the singularity structure of the process, seen at a finite resolution. As a supplementary condition, the best is that the algorithmic complexity grows only linearly with the length of the time series. These requirements are met if we project the stochastic processes in a discrete wavelet space.

Wavelets are associated with time-scale representations. Let $`|F`$ be a vector in a Hilbert space. There are bases $`\left\{|\psi _{jn}\right\}`$ such that $`|F`$ can be expanded as

$$|F=\underset{j=\mathrm{}}{\overset{j=\mathrm{}}{}}\underset{n=\mathrm{}}{\overset{n=\mathrm{}}{}}\psi _{jn}|F|\psi _{jn},$$
(1)
the form of the basis functions being $`\psi _{jn}\left(t\right)=\frac{1}{\sqrt{2^j}}\psi \left(\frac{t2^jn}{2^j}\right).`$ These functions are called discrete wavelets, and are square integrable with zero mean, and may have $`p`$ vanishing moments. The discrete wavelet transform of $`F`$ to this basis is written $`D_{jn}F\psi _{jn}|F,`$ and gives information on the behavior of the function $`F`$ at scale $`2^j`$ and time $`2^jn`$ mallat .

Confining the function in expression (1) to a finite observation window, and measuring at finite resolution, we obtain the approximation

$$|f_J=\varphi |f|\varphi +\underset{j=J+1}{\overset{j=0}{}}\underset{n=0}{\overset{n=2^j1}{}}\psi _{jn}|f|\psi _{jn},$$
(2)
where $`f`$ is the restriction of $`F`$ to $`[0,1]`$. The function $`\varphi `$ is called the scaling function, having the properties of a smoothing kernel mallat . We are interested in the fluctuations, so, without any loss, we employ in our analysis the transform $`fff`$, getting from $`f`$ a version with null time-average. At this point we are ready to express an energy and an entropy for the time series, that are also resolution dependent approximations of the energy density and entropy density for the ideal stochastic process. The energy of $`f_J`$ is defined by the ensemble average of the $`L^2`$-norm squared,

$$E\left(J\right)\overline{f_J|f_J},$$
(3)
that is directly determined with the simple formula $`E\left(J\right)=\overline{_{i=1}^{i=N}f_i^2},`$ where $`f_i`$ are the sampled values of $`f`$. The entropy has the definition

$$S\left(J\right)\underset{j=J+1}{\overset{j=0}{}}\underset{n=0}{\overset{2^j1}{}}P_{jn}\mathrm{log}_2P_{jn}P_\varphi \mathrm{log}_2P_\varphi ,$$
(4)
with $`P_{jn}\overline{|D_{jn}F|^2}/E\left(J\right)`$, and $`P_\varphi \overline{|\varphi |f|^2}/E\left(J\right)`$. This expression has similarities with the basis entropy cost, using the language of reference mallat . In the limit of a single pulse localized around $`t`$ and with duration $`\delta tk\tau `$, with $`k`$ a natural number, only very few wavelet coefficients around the scale $`k\tau `$ and localized near $`t`$ will be appreciably different from zero. This characterizes low entropy processes. On the other hand, the white noises (WN) in which the energy is, in average, equally distributed on all wavelets (the same power dissipation at all times and scales) have the maximum entropy $`S_{\mathrm{max}}=J`$. This entropy functional is useful to quantify the balance between order and disorder (complexity) for different regimes of evolution of complex systems. Actually, we expect that the level of disorder in a process mainly depends on the energy distributed into the details. Consider, for example, the interaction between two surfaces that slip one over the other. When there is no friction we expect no sound. But if we introduce roughness, the energy will be dissipated in microscopic collisions, originating a sound that is classified as a noise, due to the excess of details. Indeed, the intervals between such collisions are so small that no measurement will capture them. The process will continue a noise for ‘high’ resolution instruments, like it is for our ‘low’ resolution ears. In FIG. 1 we see the wavelet decomposition of the sound caused by a piece of plastic pressed and pushed on a rough table. The calculated entropy was $`S=12.52`$ for a maximum of $`S_{\mathrm{max}}=15`$. The random small scale details have amplitudes comparable to that of the function during the whole observation, constituting important components. The great amount of energy distributed on the multitude of details causes this high entropy. This is why we suggest that the proposed entropy is physically appropriate.

A paradigmatic model for the self-affine processes is the fractional Brownian motion (fBm), introduced by Mandelbrot and Van Ness in the realm of stochastic processes with $`1/f^\alpha `$ spectra mandelbrot . The fBm is a mono-parametric stochastic process $`B_H`$ with null ensemble average, Gaussian increments, and whose correlation function has the general form:

$$\overline{B_H\left(t\right)B_H\left(s\right)}=\frac{\sigma _H^2}{2}\left(\left|t\right|^{2H}+\left|s\right|^{2H}\left|ts\right|^{2H}\right),$$
(5)
with $`\sigma _H^2\overline{B_H^2\left(1\right)}`$. The parameter of the process is the Hurst exponent $`H(0,1)`$, which controls the nature of the correlations. The variance of $`B_H`$ depends on time according to $`\overline{B_H^2}=\sigma _H^2\left|t\right|^{2H},`$ that reduces to the Brownian motion form when $`H=1/2`$, representing the uncorrelated case. In this process, for $`H<1/2`$ the increments present infinite range anti-correlations, while for $`H>1/2`$ there are infinite range correlations. The Fourier spectrum of fBm follows the power law

$$S\left(f\right)\frac{1}{f^{2H+1}}.$$
(6)
In the Haar wavelet basis one finds exact and simple expressions for the resolution dependent energy and entropy densities of fBm. The Haar wavelets are periods of square waves written as

$$h_{jn}\left(t\right)=\{\begin{array}{cc}\frac{1}{\sqrt{2^j}},& 2^jnt<2^j\left(n+1/2\right)\\ \frac{1}{\sqrt{2^j}},& 2^j\left(n+1/2\right)t<2^j\left(n+1\right)\end{array}.$$
(7)
For Haar wavelets one proves, using the expression (5), that

$$\overline{\left|\mathrm{\Delta }_{jn}b_{H,J}\right|^2}=\frac{\sigma _H^2\left(12^{2H}\right)}{\left(2H+1\right)\left(2H+2\right)}2^{\left(2H+1\right)j},$$
(8)
where $`\mathrm{\Delta }_{jn}fh_{jn}|f`$, and $`b_H`$ indicates restriction of $`B_H`$ to $`[0,1]`$. This expression gives the scale dependent power spectrum in the Haar basis, from which one can calculate the energy density of $`B_H`$

$$E(H,J)=\overline{b_{H,J}|b_{H,J}}=\frac{\sigma _H^2\left(12^{2HJ}\right)}{\left(2H+1\right)\left(2H+2\right)},$$
(9)
depending on $`\sigma _H`$, that must be determined experimentally . The entropy density has the final form

$$S(H,J)=2^{2H}\left(2H+1\right)\frac{12^{2H(J1)}J+2^{2HJ}(J1)}{\left(12^{2H}\right)\left[12^{2HJ}\right]}+\mathrm{log}_2C(H,J),$$
(10)
where $`C(H,J)=\left[12^{2HJ}\right]/\left(12^{2H}\right)`$.

This entropy density is universal and exact for the fBm processes sampled without rounding error to a given finite resolution. It is very useful since for a great number of practical applications there are just low resolution data, what makes mandatory that the theoretical predictions take resolution into account explicitly. In the limit $`J\mathrm{}`$ the entropy diverges at $`H=0`$ as $`S(H,\mathrm{})1/2H\mathrm{ln}2`$. For other values of $`H`$ the entropy is finite at infinite resolution. At finite resolution $`S(0,J)=(J1)/2+\mathrm{log}_2(J)`$.

The FIG. 2 shows a plot of fBm theoretical entropy for the resolution $`N=2048`$. The theoretical result is compared with $`S(H,J)`$ calculated from samples simulated with the Successive Random Additions (SRA) algorithm peitgen . It shows also how Gaussian white noise (WN) affects the entropy, clarifying what is to be expected from noisy data. For comparison, the dots at the upper part of the figure show the entropy for the simulated pure Gaussian white noise, whose average value is $`S=10.9811`$. In the dotted line we see the effect of a superposed noise with 10% of the amplitude of the fBm. The dash-dotted line shows this effect when the percentage grows to 50%. The entropies for all simulated processes were calculated from 100 samples. The drastic diminution of entropy when $`H`$ increases happens because the details (roughness) become weaker. There is very good agreement between theoretical prediction and simulation.

In conclusion, we propose that the entropy measure (4) is suitable to estimate the complexity of time series by revealing the degree in which *detail modes* are excited in a process. These detail modes are easily represented as localized oscillations or wavelets, which in turn are very effective mathematical instruments for time-scale analysis. Considering such phenomena as fluid or scalar turbulence zsshe ; shraiman , that are characterized by fluctuating cascades from high inertial scales to low dissipative scales, or plasma turbulence krommes , in which some scales are most important for the energy transfer, the proposed entropy seems a natural measure of the intrinsic order-disorder balance. In such phenomena, many structures (vortices, convective cells) appear at intermediate scales, that should be detected as excitations of intermediate scale modes, lowering the entropy, as compared with noises that are completely determined by microscopic dissipation (whose characteristics are indicated in FIG. 1). For example, if a theory with the same correlations of Kolmogorov K41 theory kolmogorov41  is valid, implying $`H=1/3`$, the time series of velocity increments would lead to the entropy $`S(1/3,\mathrm{})=4.27\mathrm{}`$, a number significantly small, compared to the case in FIG. 1. What we learn from the fBm analysis presented here suggests that the measurement of the entropy is a strong method of characterization. Furthermore, it is interesting to obtain experimental curves of entropy for stochastic processes depending on parameters, or in situations where the entropy may change with time. This method has given results (which will be published soon) in three situations: In the analysis of the sound of lungs, in the analysis of heart beat series, and in the analysis of global positioning system errors caused by equatorial spread F in ionospheric plasma. We observe that such entropy measurements could give important information on the non-equilibrium approach to self-organized regimes corral , as well as, on self-organized regimes made unstable, for example, by diffusion newman .

We acknowledge the CAPES, CNPq, FINEP and CTPETRO for financial support.


format: text

text:
# Physical Nucleon Properties from Lattice QCD

## Abstract

We demonstrate that the extremely accurate lattice QCD data for the mass of the nucleon recently obtained by CP-PACS, combined with modern chiral extrapolation techniques, leads to a value for the mass of the physical nucleon which has a systematic error of less than one percent.

preprint: ADP-03-112/T550

Hadronic physics presents fascinating theoretical challenges to the understanding of strongly interacting systems in terms of their fundamental degrees of freedom in QCD, quarks and gluons. Traditional perturbative methods in quantum field theory cannot account for the highly non-trivial structure of the QCD vacuum – which is dominated by large scale quark and gluon correlations. Lattice gauge theory Wilson:1974sk  has so far provided the only rigorous method of solving non-perturbative QCD. We will show that recent progress within the field AliKhan:2001tx , together with advances in effectively field theory (EFT) Young:2002ib , now permit the accurate, model-independent extraction of the physical properties of hadrons from lattice QCD simulations, even though it is not yet feasible to make calculations at the physical quark mass.

It is now possible to make extremely accurate calculations of the masses of hadrons, such as the nucleon, with dynamical fermions. Indeed the CP-PACS group has just reported data with a precision of order 1% AliKhan:2001tx . However, such precise evaluations are limited to quark masses an order of magnitude larger than those found in nature. In order to compare with experiment, which is after all one of the main aims in the field, it is therefore necessary to extrapolate in quark mass. It is well known that such an extrapolation is complicated by the unavoidable non-analytic behaviour in quark mass, which arises from Goldstone boson loops in QCD with dynamically broken chiral symmetry Pagels:1975se .

Early work motivated by the important role of Goldstone bosons led to the construction of chiral quark models Thomas:1984kv , which incorporated this non-analytic behavior. An alternative, systematic approach, designed to avoid reference to a model, involved the construction of an effective field theory (EFT) to describe QCD at low energy Weinberg:1979kz . The application to baryons has developed to the point where chiral perturbation theory ($`\chi `$PT) is now understood as a rigorous approach near the chiral limit Gasser:1988rb ; Bernard:1995dp .

Because it is defined as an expansion in momenta and masses about the chiral limit, $`\chi `$PT provides an attractive approach to the problem of quark mass extrapolation for lattice QCD. The advantages of formulating $`\chi `$PT with a physical regularisation scheme, as opposed to the commonly implemented dimensional regularisation, have been demonstrated by Donoghue et al. Donoghue:1998bs . Early implementations of a finite-range regulator (FRR) to evaluate chiral loop integrals in $`\chi `$PT suggested that, in the context of the extrapolation of lattice data from relatively large quark masses, FRR provides a more reliable procedure Leinweber:1999ig . There has been considerable debate on whether current lattice data is within the scope of dimensionally-regularised $`\chi `$PT or whether the form of the FRR chosen introduces significant model-dependence Bernard:2002yk . However, this issue has now been resolved Young:2002ib  through a recent detailed study of numerous regularisation schemes in $`\chi `$PT, both dimensional and FRR. This study quantified the applicable range of the EFT Young:2002ib  and established that all of the FRR considered provided equivalent results over the range $`m_\pi ^20.8\mathrm{GeV}^2`$.

Here we demonstrate that by adopting the FRR formulation of $`\chi `$PT one can improve its convergence properties to the point where the chiral extrapolation to the physical quark mass can be carried out in a model-independent way with a systematic uncertainty of less than one percent. Given this remarkable result, the current errors on extrapolated quantities are dominated by the statistical errors arising from the large extrapolation distance – the lightest simulated pion mass being typically $`m_\pi ^20.27\mathrm{GeV}^2`$, in comparison with the physical value, $`0.02\mathrm{GeV}^2`$.

In the usual formulation of effective field theory the nucleon mass as a function of the pion mass (given that $`m_\pi ^2m_q`$ Gell-Mann:1968rz ) has the formal expansion:

$`M_N`$ $`=`$ $`a_0+a_2m_\pi ^2+a_4m_\pi ^4+a_6m_\pi ^6+\mathrm{}`$ (1)
$`+\sigma _{\mathrm{N}\pi }+\sigma _{\mathrm{\Delta }\pi }.`$
In principle the coefficients, $`a_n`$, can be expressed in terms of the parameters of the underlying effective Lagrangian to a given order of chiral perturbation theory. In practice, for current applications to lattice QCD, the parameters must be determined by fitting to the lattice results themselves. The additional terms, $`\sigma _{\mathrm{N}\pi }`$ and $`\sigma _{\mathrm{\Delta }\pi }`$, are the two self-energy terms involving the (Goldstone) pion which yield the leading and next-to-leading non-analytic behaviour of $`M_N`$. As these terms involve the coupling constants in the chiral limit, which are essentially model-independent Li:1971vr , the only additional complication they add is that the ultra-violet behaviour of the loop integrals must be regulated in some way.

Traditionally, one uses dimensional regularisation which (after infinite renormalisation of all even powers of $`m_\pi `$) leaves only the non-analytic terms – $`c_{\mathrm{LNA}}m_\pi ^3`$ and $`c_{\mathrm{NLNA}}m_\pi ^4\mathrm{ln}(m_\pi /\mu )`$, respectively. Within dimensional regularisation one then arrives at a truncated power series for the chiral expansion,

$`M_N`$ $`=`$ $`c_0+c_2m_\pi ^2+c_{\mathrm{LNA}}m_\pi ^3+c_4m_\pi ^4`$ (2)
$`+c_{\mathrm{NLNA}}m_\pi ^4\mathrm{ln}{\displaystyle \frac{m_\pi }{\mu }}+c_6m_\pi ^6+\mathrm{},`$
where the bare parameters, $`a_i`$, have been replaced by the finite, renormalised coefficients, $`c_i`$. Through the chiral logarithm one has an additional mass scale, $`\mu `$, but the dependence on this is eliminated by matching $`c_4`$ to “data” (in this case lattice QCD). Provided the series expansion in Eq. (2) is convergent over the range of values of $`m_\pi `$ where the lattice data exists, one can then use Eq. (2) to evaluate $`M_N`$ at the physical pion mass. Unfortunately there is considerable evidence that this series is not sufficiently convergent Young:2002ib ; Bernard:2002yk ; Thomas:2002sj ; SVW ; Hatsuda:1990tt .

The physical reason for the lack of convergence seems to be that, apart from the usual symmetry breaking scale in $`\chi `$PT, $`\mathrm{\Lambda }_{\chi \mathrm{SB}}4\pi f_\pi `$ 1 GeV, there is another mass scale in the problem which does not appear explicitly in a dimensionally regularised theory. This additional scale, $`\mathrm{\Lambda }`$, is the inverse of the size of the nucleon, the source of the pion cloud Detmold:2001hq . Of course, this scale must appear implicitly in the correct effective field theory as a correlation between the coefficients $`c_n`$. However, since we want to access pion masses above 0.5 GeV and the nucleon size of 0.5–1.0 fm suggests $`\mathrm{\Lambda }0.4`$ GeV, it is clear that $`m_\pi /\mathrm{\Lambda }>1`$ and the naive series should not be expected to converge over the range of pion masses of interest. This problem can be overcome through the use of a FRR Young:2002ib ; Donoghue:1998bs .

Clearly, in extracting information from lattice data for physical hadrons one must avoid model dependence. In line with the implicit $`\mu `$-dependence of the coefficients in the familiar dimensionally regulated $`\chi `$PT, in the case of a FRR the systematic expansion of the nucleon mass is:

$`M_N`$ $`=`$ $`a_0^\mathrm{\Lambda }+a_2^\mathrm{\Lambda }m_\pi ^2+a_4^\mathrm{\Lambda }m_\pi ^4`$ (3)
$`+\sigma _{\mathrm{N}\pi }(m_\pi ,\mathrm{\Lambda })+\sigma _{\mathrm{\Delta }\pi }(m_\pi ,\mathrm{\Lambda }),`$
where the dependence on the shape of the regulator is implicit. As for the mass scale, $`\mu `$, entering the usual chiral logarithms, the dependence on the value of $`\mathrm{\Lambda }`$ and the choice of regulator is once again eliminated, to the order of the series expansion, by fitting the coefficients, $`a_n^\mathrm{\Lambda }`$, to lattice data. The clear indication of success in eliminating model dependence and hence having found a suitable regularisation method, is that the higher order coefficients (for the best fit value of $`\mathrm{\Lambda }`$) should be small and that the result of the extrapolation should be insensitive to the shape of the ultraviolet regulator.

In order to investigate the model dependence associated with the self-energies, several regulators are considered. We evaluate the loop integrals in the heavy baryon limit

$$\sigma _{BB^{}}^\pi =\frac{3}{16\pi ^2f_\pi ^2}G_{BB^{}}_0^{\mathrm{}}𝑑k\frac{k^4u^2(k)}{\omega (k)(\omega _{BB^{}}+\omega (k))},$$
(4)
taking $`u(k)`$ to be either a sharp cut-off, $`\theta (\mathrm{\Lambda }k)`$, a dipole, $`(1+k^2/\mathrm{\Lambda }^2)^2`$, a monopole, $`(1+k^2/\mathrm{\Lambda }^2)^1`$ or finally a Gaussian, $`\mathrm{exp}(k^2/\mathrm{\Lambda }^2)`$. We note that the dipole is the most physical, in that the observed axial form factor of the nucleon has this shape Thomas:2001kw . In Eq. (4) for $`B=B^{}=N`$ we have $`G_{NN}=g_A^2`$ (with $`g_A=1.26`$), while $`G_{N\mathrm{\Delta }}=32g_A^2/25`$, the SU(6) value. In addition, $`\omega (k)=\sqrt{k^2+m_\pi ^2}`$, $`\omega _{NN}=0`$ and $`\omega _{N\mathrm{\Delta }}=292`$ MeV, the physical $`\mathrm{\Delta }`$-$`N`$ mass splitting. Clearly the four choices of regulator have very different shapes, with the only common feature being that they suppress the integrand for momenta greater than $`\mathrm{\Lambda }`$.

In addition to the FRR expansions of Eq. (3) and the standard dimensionally regulated expansion of Eq. (2), we also consider the case where Eq. (2) is modified to maintain the square-root branch-cut at $`m_\pi =\omega _{N\mathrm{\Delta }}`$ Young:2002ib . One can now compare the expansion about the chiral limit for these six different regularisation schemes in order to assess their rate of convergence. It turns out that all the FRR expansions precisely describe the dimensional regularization expansion over the range $`m_\pi ^2(0,0.7)`$ GeV2. Furthermore, the smooth, FRR formulations are consistent with each other, for the renormalized chiral coefficients, $`c_{0,2,4}`$, to an extraordinarily precise level Young:2002ib . This ensures a reliable extrapolation to the regime of physical quark masses.

We now exploit the accuracy of the most recent CP-PACS data AliKhan:2001tx  by using it to determine either the first four unknown parameters, $`c_0c_6`$, appearing in Eq. (2) (i.e. using dimensional regularisation, with and without the $`\mathrm{\Delta }N\pi `$ branch point) or the four parameters $`a_0,a_2,a_4`$ and $`\mathrm{\Lambda }`$ in Eq. (3), for each choice of ultra-violet regulator r0scale . It is only the extreme accuracy of the data which makes the determination of as many as four parameters possible.

Figure 1 shows the resulting fits to the lattice data over the range $`m_\pi ^2(0,1.0)\mathrm{GeV}^2`$, with the corresponding parameters given in Table 1. It is remarkable that three of the four curves based on Eq. (3), the dipole, monopole and Gaussian cases, are indistinguishable on this plot. Furthermore, we see from Table 1 that the coefficient of $`m_\pi ^4`$ in all of those cases is quite small – an order of magnitude smaller than the higher order coefficients obtained using Eq. (2), which also oscillate in sign. This indicates that the residual series, involving $`a_{0,2,4}`$, is converging when the chiral non-analytic behaviour is evaluated with a FRR.

As explained by Donoghue et al. Donoghue:1998bs , one can combine the order $`m_\pi ^{0,2,4}`$ terms from the self-energies with the “bare” expansion parameters, $`a_{0,2,4}`$, to obtain physically meaningful renormalized coefficients. These are shown in Table 2, in comparison with the corresponding coefficients found using Eq. (2). Details of this renormalization procedure are given in Ref. Young:2002ib .

The degree of consistency between the best-fit values found using all three smooth choices of FRR is remarkably good, whereas the large numerical values and the change of sign in the coefficients found from the fit using Eq. (2) suggest that the latter is nowhere near convergence. We can understand the problem very simply; it is not possible to accurately reproduce the necessary $`1/m_\pi ^2`$ behaviour of the chiral loops (for $`m_\pi >\mathrm{\Lambda }`$) with a 3rd order polynomial in $`m_\pi ^2`$. The final point to notice about the finite range regulators is that all but the sharp cutoff naturally yield a correction to the Goldberger-Treiman relation of approximately the right size. As shown by Birse and McGovern McGovern:1998tm , this correction dominates the non-analytic term of order $`m_\pi ^5`$ which only arises at two loops in dimensionally regulated $`\chi `$PT. The absence of the correct $`m_\pi ^5`$ term in the case of the sharp cutoff explains its slight deviation from the results of the other FRR.

It is clear that the use of an EFT with FRR enables one to make a model-independent extrapolation of the nucleon mass as a function of the quark mass. Although minimal deviation is seen between the best fit curves, we need to determine how well these curves are in fact constrained by the statistical uncertainties of the lattice data. As all data points are statistically independent, the one-sigma deviation from the best-fit curve is defined by the region for $`(\chi ^2\chi _{\mathrm{min}}^2)/dof<1`$. We use a standard $`\chi ^2`$ measure, weighted by the squared error of the simulated data point and $`\chi _{\mathrm{min}}^2`$ corresponds to the optimum fit to the data. To preserve the infrared physics of EFT, we put a lower bound of $`300\mathrm{MeV}`$ (just above the $`N\mathrm{\Delta }`$ mass splitting) on the regulator parameter, $`\mathrm{\Lambda }`$.

We show the one-sigma variation from the best-fit dipole curve by the shaded region in Figure 2. The extrapolated values for the nucleon mass are shown in Table 3. It is clear that the large lower-deviation on the extrapolated values are a result of the regulator mass not being constrained from above. Alternatively, for the FRR cases, one could fix $`\mathrm{\Lambda }`$ and vary an additional parameter, $`a_6`$, in the spirit of dimensional regularisation. In that case the error bounds are symmetric and similar in magnitude to those found for the dimensionally regulated cases reported in Table 3. The primary source of the large error band is the large extrapolation distance.

We also quote the $`\sigma `$-commutator in Table 3 – calculated as $`m_\pi ^2M_N/m_\pi ^2`$ at the physical pion mass Leinweber:2000sa . Again the statistical errors are large. In this case, as well as for the nucleon mass itself, it is clear that a simulated point at a pion mass $`m_\pi ^20.1\mathrm{GeV}^2`$ would greatly reduce the statistical error in the extrapolation. We note that our study has focused on the systematic errors associated with chiral extrapolation, while those arising from finite lattice spacing and volume Young:2002cj ; SVW  remain to be quantified.

To summarize, we have shown that the extremely precise dynamical simulation data from CP-PACS permits us to determine four parameters in the chiral extrapolation formulas, Eqs. (2) and (3). Whereas the former (involving dimensional-regularization) does not appear to be convergent over the required mass range, the improved convergence properties of a finite range regulator yield an excellent description of the data over the full mass range, regardless of the functional form chosen for the vertex regulator. Table 3 summarises the resulting values of the nucleon mass and the $`\sigma `$-commutator evaluated at the physical pion mass. The nucleon mass is constrained by the chiral extrapolation technique based upon a FRR procedure to a systematic error less than $`1\%`$. The $`\sigma `$-commutator is similarly constrained to within a systematic error of $`3\%`$. The values obtained are both consistent with earlier extractions from lattice data using a similar technique Leinweber:2000sa  and, within the relatively large statistical errors, in agreement with the experimental values Gasser:1991ce . We note also that the systematic uncertainty in the extraction of the low energy constant, $`c_0`$, is less than 1%, while for $`c_2`$ it is at the level of a few percent. With the issue of model independence resolved, there is an urgent need for high precision lattice QCD simulations at $`m_\pi ^20.1`$ GeV2 in order to reduce the present statistical error on the extrapolation. Such simulations should be feasible with the new generation of lattice QCD computers currently under construction.

###### Acknowledgements.

We thank M. Birse, J. McGovern and S.V. Wright for helpful conversations. This work was supported by the Australian Research Council.


text:
# First-principles study of (BiScO3)1-x-(PbTiO3)x piezoelectric alloys

## I Introduction

Perovskite alloys based on PbTiO3 (PT) are of considerable interest for applications as piezoelectric actuator materials. The phase diagrams of the most technologically important alloys are characterized by a morphotropic phase boundary (MPB) separating the PT-rich tetragonal phase, in which the polarization lies along a $`001`$ direction, from the PT-poor rhombohedral phase, in which the polarization is along a $`111`$ direction. Because the structural transition between these two phases brings about a large electromechanical response, materials with a composition lying close to the MPB are the preferred ones for applications. Examples of such materials are PZT or PZN-PT, in which PT is alloyed with PbZrO3 (PZ) or Pb(Zn1/3Nb2/3)O3, respectively.

Recently, attention has been drawn to a new class of materials in which PT is alloyed with Bi-based perovskites, the best-studied example being (BiScO3)1-x-(PbTiO3)x (BS-PT). The appeal of BS-PT is twofold: (i) it appears to have piezoelectric properties comparable to those of PZT and PZN-PT in quality, and (ii) its dielectric and piezoelectric properties should be more robust, over a wider temperature range, than those of PZT and PZN-PT. The physical reason behind (ii) is that the Curie temperature within the MPB composition range ($`T_C^{\mathrm{MPB}}`$) is higher in BS-PT ($``$ 450C) than in PZT ($``$ 400C) or PZN-PT ($``$ 200C), suggesting that BS-PT at room temperature will age more slowly, its properties will be more temperature-independent, etc. A possible drawback is that BS-PT is evidently not thermodynamically stable in the perovskite crystal structure over the entire range of composition. However, for applications purposes it is only important that the perovskite structure can be obtained for compositions near the MPB, and this is indeed the case.

The authors of Ref.  synthesized BS-PT guided by an empirical crystal-chemistry rule suggesting that since Bi3+ is too small an ion to form a cubic perovskite with Sc3+, it will have a strong tendency to move off-center from its high-symmetry position, yielding a high structural transition temperature. While perhaps overly simplistic, this reasoning is undoubtedly partially correct and is likely to be useful for materials engineering. However, it is unclear how to explain the large piezoelectric responses found in BS-PT based on ion-size considerations alone. The ferroelectric properties of materials like PZT are known to be related to the partial covalency of some bonds, and the situation in BS-PT should be similar. On the other hand, most of the relevant piezoelectric alloys have only Pb on the A site of their perovskite structure, and the chemistry of Pb is known to be very important to their properties. The substitution of Pb by Bi in BS-PT is thus of particular interest in view of the good piezoelectric properties of this material.

We present here a first-principles study of BS-PT. We find that Bi plays a crucial role, and in particular that hybridization between Bi-6$`p`$ and O-2$`p`$ orbitals is the driving force for the ferroelectric instabilities of BS-PT, allowing for very large polarizations and responses. Our results are generally consistent with experiment, but they also suggest that the intrinsic ferroelectric and piezoelectric properties of BS-PT alloys might be even better than those measured experimentally to date.

We present our results in the form of a systematic comparison of the properties of BS-PT with those of PZT. Section II describes the first-principles methods employed. In Sec. III we study the ferroelectric instabilities of pure BS, PT, and PZ, introducing some corresponding calculations on BiYO3 (BY) as an aid to the discussion. In Sec. IV we discuss the ferroelectric instabilities of the alloys and make contact with the experimental results, while Sec. V is devoted to the piezoelectric properties of the alloys near the MPB. We summarize and conclude in Sec. VI.

## II Methods

The calculations were performed within a plane-wave implementation of the local-density approximation (LDA) to density-functional theory, using ultrasoft pseudopotentials  to represent the ionic cores. The following electronic states were included in the calculation: the 2$`s`$ and 2$`p`$ states of O, the 3$`s`$, 3$`p`$, 3$`d`$, and 4$`s`$ states of Sc and Ti, the 4$`s`$, 4$`p`$, 4$`d`$, and 5$`s`$ states of Y and Zr, and the 5$`d`$, 6$`s`$, and 6$`p`$ states of Pb and Bi. For the 5-atom unit cell calculations, we used a 35 Ry plane-wave cutoff and a $`6\times 6\times 6`$ k-point grid for Brillouin zone integrations. We checked that these choices yield converged results. We also made some calculations with a rocksalt-ordered (fcc) 10-atom unit cell, for which we used a 30 Ry cutoff and a $`4\times 4\times 4`$ k-point grid. Atomic relaxations were considered to be converged when the residual forces were smaller than $`5\times 10^5`$ a.u. Polarizations were computed using the Berry-phase expression of King-Smith and Vanderbilt  and always using a dense-enough k-point mesh along the direction of the strings. The character of the electron energy bands was estimated by computing the partial density associated with the pseudopotential projectors of specific atoms (see Eq. (19) of Ref.  for guidance).

Calculations using the virtual-crystal approximation (VCA) were implemented as in Ref. . In the VCA one works with virtual atoms defined by a pseudopotential that is a weighted average of those of real atoms. Such a procedure may lead to meaningful estimates of the average properties of an alloy provided that the valence orbitals of the real atoms forming the virtual atom are similar in character. This can be expected to be true for the pairs Bi-Pb and Sc-Ti in BS-PT. Also, both groups BiSc and PbTi have a total nominal ionic charge of $`+6`$, allowing them to be mixed in arbitrary proportions. These requirements are also satisfied by PZT, which has previously been studied within the VCA. It is important to note that the VCA is not expected to render an accurate description of the electronic structure of a truly disordered alloy, but rather to provide us with qualitative energetic and structural trends. Also, the VCA is expected to work better for isovalent pairs of atoms (e.g., Zr-Ti) than for heterovalent pairs (e.g., Bi-Pb or Sc-Ti), which makes BS-PT a particularly challenging system. For this reason, we assess the reliability of our VCA results by comparing them with some supercell calculations and available experimental information.

## III Pure compounds

As a preliminary step towards understanding BS-PT, we first consider pure BS (in the perovskite structure) and compare with corresponding calculations for pure PT and PZ. The first three lines of Table I show the structural and polarization results obtained by relaxing these systems within ferroelectric phases of tetragonal ($`P4mm`$) or rhombohedral ($`R3m`$) symmetry, always maintaining a 5-atom unit cell. The atomic relaxations corresponding to the tetragonal phases of BS and PT are given in Table II.

For BS, we find the rhombohedral phase to be lower in energy than the tetragonal phase, which agrees with the experimental fact that BS-rich BS-PT is rhombohedral, and we find that the relaxations are very large. This is reflected both in the structural data (we obtain $`c/a=1.29`$ for the tetragonal phase of BS, compared with 1.05 for PT; see also the atomic displacements in Table II) and in the energetics ($`\mathrm{\Delta }E_{\mathrm{tet}}`$ and $`\mathrm{\Delta }E_{\mathrm{rho}}`$ of BS are one order of magnitude larger than for PT and PZ). These results suggest that the cubic perovskite structure is not a very natural one for BS, a conclusion which is reinforced by experimental indications that pure BS at zero pressure is not even thermodynamically stable.

The polarization associated with the tetragonal phase of BS is large, while that of the rhombohedral phase is relatively small. In both cases, the largest contribution to $`P`$ comes from the displacement of Bi relative to its neighboring O ions (Ox and Oy in Table II). More specifically, in the cubic phase of BS all effective charges are normal, i.e., with values close to the nominal ionic charges, except for those of Bi (6.4) and Ox and Oy ($``$3.2 along the polarization direction). Moreover, these are the ions with the largest relative displacements associated with the ferroelectric instability (see Table II). The differnce in $`P`$ between the two phases reflects a significantly smaller charge transfer when the Bi ion approaches three O ions (rhombohedral phase) instead of four (tetragonal phase).

The fundamental role of Bi in the ferroelectric instabilities of BS is best appreciated in the electronic band-structure data shown in Fig. 1. As in PT (panel c), the top valence states of the cubic phase of BS (panel a) are O-2$`p`$ in character. However, the lowest unoccupied states of BS are Bi-6$`p`$ like, at variance with the Ti-3$`d`$ bands we find in PT. This fact, together with the relatively small gap of cubic BS, results in ferroelectric instabilities mainly characterized by a hybridization between Bi-6$`p`$ and O-2$`p`$ orbitals, not the one involving $`d`$ bands of the B cation that is typical of most other ferroelectric perovskites. In the tetragonal phase of BS (panel b), the gap opens considerably and the top valence bands increase their Bi-6$`p`$ character. The strong Bi-6$`p`$–O-2$`p`$ hybridization is also reflected in the fact that the effective charge of Bi in the tetragonal phase recovers a normal value (2.9 along the polarization direction).

As indicated earlier, Eitel and co-workers chose to study BS-PT because of the large size mismatch between Bi and Sc. Such a mismatch is conveniently quantified in terms of the tolerance factor

$$t=\frac{R_\mathrm{A}+R_\mathrm{O}}{\sqrt{2}(R_\mathrm{B}+R_\mathrm{O})},$$
(1)
where $`R_\mathrm{O}`$, $`R_\mathrm{A}`$ and $`R_\mathrm{B}`$ are the ionic radii of oxygen and of the A and B cations respectively. That is, $`t`$ is the ratio between the ideal cubic lattice parameters based on A–O vs. B–O bonding alone, with $`t`$=1 corresponding to high perovskite stability. As can be seen from Table I, the relatively low value of $`t`$=0.907 for BS implies that the Bi ion should be very unstable in the high-symmetry position of the cubic perovskite structure, consistent with our findings. We also studied BiYO3 (BY), for which $`t`$=0.845, and found even stronger ferroelectric instabilities that are very similar in character to those of BS (see Tables I and II). However, $`P_{\mathrm{tet}}`$ and $`P_{\mathrm{rho}}`$ of BY are almost identical to those of BS, not larger. Also, the distance between Bi and its bonded O neighbors is roughly the same in the ferroelectric phases of BS and BY (4.33 and 4.35 a.u., respectively, in the tetragonal phase). These two facts seem to indicate that the onset of the Bi-$`6p`$–O-$`2p`$ hybridization (which also dominates the development of polarization in BY) is not very $`t`$-dependent, suggesting that the large polarizations of BS and BY cannot be attributed mainly to a size-mismatch effect. Instead, they are strongly linked to the particular chemistry of Bi, i.e., its strong tendency to bond covalently with the surrounding oxygen anions.

## IV Alloys: Ferroelectric instabilities

We now turn to a study of the BS-PT alloy within the virtual crystal approximation (VCA). We also carry out parallel VCA calculations on the PZT alloy as a relevant reference.

Figure 2 shows the energies of the tetragonal and rhombohedral ferroelectric phases of BS-PT and PZT. (The cubic phase, which is material- and composition-dependent, is chosen as the zero of energy.) Structural and polarization results corresponding to those obtained for the pure compounds are shown in Tables I and II.

Associating the MPB composition with the crossing of the energy curves for the two ferroelectric phases in Fig. 2, we estimate $`x_{\mathrm{MPB}}0.65`$ for BS-PT, while the experimental value is 0.64. Similarly, for PZT we get $`x_{\mathrm{MPB}}0.55`$, to be compared with experimental values from about 0.45 to about 0.52 (the range of stability of the intermediate monoclinic phase of PZT). Taking into account all the approximations involved in our estimation, we must regard the agreement in the case of BS-PT as partly fortuitous. Nevertheless, these results suggest that our approach is adequate.

Figure 2 also shows that the energy reductions associated with the ferroelectric instabilities of BS-PT are considerably larger than those of PZT. In particular, at the MPB we have $`\mathrm{\Delta }E^{\mathrm{MPB}}=`$0.234 eV for BS-PT and $``$0.056 eV for PZT. These results are reminiscent of the highly unstable cubic phase of BS previously discussed.

It is tempting to try to connect these numbers with the Curie temperatures of BS-PT and PZT at the MPB. We know that, in principle, the thermal energy required for the system to remain in the high-symmetry paraelectric phase is of the order of $`|\mathrm{\Delta }E^{\mathrm{MPB}}|`$, which suggests that $`T_C^{\mathrm{MPB}}`$ of BS-PT is about 5 times larger than that of PZT, i.e., higher than 1000C. However, we also know that in systems in which there is competition between different structural instabilities, such an argument can grossly overestimate the transition temperature. In PZT, for example, it is known that there is a competition between the ferroelectric rhombohedral (FE-rho) phase and an antiferrodistortive (AFD) phase involving rotations of the oxygen octahedra. Such a competition also exists in BS-PT, since we find that this material presents a strong AFD instability. Taking PZT as a reference, one could try to estimate the effect of this competition in BS-PT in the following way. The AFD phase supposes a strain of the high-symmetry structure that is opposite to that of the FE-rho phase. At $`x`$=0.5 for BS-PT the volume per 5-atom cell of the FE-rho phase at the MPB is 415 a.u., while it is 396 a.u. for the AFD phase, resulting in a volume ratio of 1.05. This ratio is 1.01 in PZT, which suggests that the FE-rho–AFD competition will be stronger in BS-PT than in PZT. Energetically, though, we find the opposite trend. The FE-rho phase is more favorable by 0.15 eV in the case of BS-PT, while in PZT FE-rho is preferred to AFD by only 0.02 eV. Hence, we are unable to determine in which of the two systems the reduction of $`T_C^{\mathrm{MPB}}`$ caused by the FE-rho–AFD competition should be larger. On the other hand, there is another mechanism missed by the VCA and which could significantly lower the transition temperatures of BS-PT. Because BS-PT is a disordered heterovalent alloy, internal electric fields will be present in it. Such fields will cause local frustrations and enhance competition between different polarization directions, resulting in a decrease of $`T_C`$. In isovalent PZT the internal fields and their effects are bound to be much smaller. Quantifying these effects would require a full statistical-mechanics study of the disordered alloys, but that is beyond the scope of this work. In conclusion, all we can say about the issue of the transition temperatures is that our results indicate that $`T_C^{\mathrm{MPB}}`$ of BS-PT could well be significantly larger than that of PZT.

Let us now discuss the structural and polarization data. The ionic and strain relaxations of BS-PT are much larger than those of PZT (see Tables I and II). Importantly, at the MPB composition we find $`c/a=`$1.079 and $`P_{\mathrm{tet}}=0.875`$ C/m2 for BS-PT, to be compared, respectively, with 1.027 and 0.687 C/m2 obtained for PZT. The very large values calculated for BS-PT in the MPB region suggest that its piezoelectric response could be very large too, assuming the usual picture in which this arises from easy polarization rotation associated with the near-degeneracy between tetragonal and rhombohedral phases. In fact, according to our results, it would not be surprising to find that BS-PT has even better piezoelectric properties than PZT.

Here we have to note that, while our results for PZT are in reasonable agreement with experiment (at the MPB, $`c/a`$ and $`P`$ of PZT are measured to be about 1.025 and 0.75 C/m2, respectively), our results for BS-PT significantly deviate from experimental values. Ref.  reports $`c/a1.023`$ and $`P_{\mathrm{tet}}0.32`$ C/m2 at the MPB of BS-PT, values that are much smaller than those that result from our VCA calculations. In view of this discrepancy, we decided to asses the reliability of the VCA by studying a 10-atom supercell of BS-PT at 50/50 composition. We have chosen an fcc supercell in which there is rocksalt ordering both on the A sublattice (Bi/Pb) and on the B sublattice (Sc/Ti). The resulting supercell has $`T_d`$ symmetry, so that all tetragonal ferroelectric domain orientations remain degenerate, while the the rhombohedral domains split in two groups of four degenerate orientations. (Here we report results for the lower-energy rhombohedral domains only.) The results, summarized in Table I, display a significant reduction in the $`c/a`$ of the tetragonal phase, but not large enough to reconcile theory with experiment. The VCA error is presumably maximum at $`x=0.5`$, so at $`x_{\mathrm{MPB}}`$ we would expect that $`c/a`$ should be at least 1.045, still significantly larger than the experimental value. Regarding the polarization, the 10-atom supercell values are actually larger than the VCA results.

In Table I we also included results corresponding to the analogous 10-atom supercell calculations for PZT. As expected, the VCA is more accurate for PZT (isovalent alloy) than for BS-PT (heterovalent). This is further confirmed by the atomic relaxation data (not shown here).

Hence, for BS-PT a relatively large discrepancy between theory and experiment remains. Its origins are most probably related to the differences between the idealized crystalline samples we consider and the ceramic samples actually grown (the internal fields present in real disordered samples are likely to play a particularly relevant role). In any case, we hope that the very large $`c/a`$ and polarizations we obtain at the MPB will stimulate further experimental study of this promising system.

## V Alloys: Piezoelectric response near the MPB

Finally, we compute the piezoelectric response of BS-PT, and compare with that of PZT. First we investigate the coefficients $`e_{ij}=dP_i/d\eta _j`$ that describe the so-called inverse piezoelectric effect ($`\eta _j`$ is the $`j`$-th component of the strain tensor in Voigt notation). We report here our results for $`e_{15}`$ of the tetragonal phase, the coefficient known to be responsible for the large piezoelectric response of PZT.  We focus on VCA alloys close to the MPB on the tetragonal side; specifically, we consider BS-PT at $`x=0.7`$ (BS-PTx=0.7) and PZT at $`x=0.6`$ (PZTx=0.6).

Following Ref. , we split the piezoelectric coefficient into two parts

$$e_{15}=e_{15,\mathrm{i}}+e_{15,\mathrm{c}},$$
(2)
where $`e_{15,\mathrm{c}}`$ is computed at fixed internal atomic coordinates (“clamped-ion” term) and

$$e_{15,\mathrm{i}}=\underset{s}{}\frac{ea}{v}Z_{11}^{}(s)\frac{du_1(s)}{d\eta _5},$$
(3)
(“internal-strain” term). Here $`s`$ runs over the atoms in the unit cell, $`e`$ is the magnitude of the electron charge, $`Z_{11}^{}`$ is the effective charge associated with displacements $`u_1`$ along the $`x`$ direction, $`a`$ is the lattice constant in the $`x`$ direction, and $`v`$ is the unit cell volume.

For BS-PTx=0.7 we obtain $`e_{15,\mathrm{c}}=`$0.93 C/m2 and $`e_{15,\mathrm{i}}=8.18`$ C/m2, which add up to $`e_{15}=7.25`$ C/m2. For PZTx=0.6 we obtain $`e_{15,\mathrm{c}}=`$0.66 C/m2, $`e_{15,\mathrm{i}}=`$8.24 C/m2, and $`e_{15}=7.58`$ C/m2. Hence, both systems have a similar $`e_{15}`$, and in both the internal-strain contribution strongly dominates the response. To understand this, we note that a large value of $`e_{15,\mathrm{i}}`$ can arise in Eq. (3) either from a “dielectric effect” (large $`Z^{}`$) or an “elastic effect” (large response of the internal coordinates to strain, i.e., large $`du_1/d\eta _5`$). Table III shows the calculated values of $`Z_{11}^{}`$ and $`du_1/d\eta _5`$ for all atoms in the unit cell of VCA BS-PTx=0.7 and PZTx=0.6. In both systems the response is mainly associated with an enhanced covalent bonding between the A ion and oxygens Oy and Oz, as pictorially shown in Fig. 3. Our results further suggest that the substitution of Pb by Bi enhances such a mechanism, consistent with the fact that Bi is more electronegative than Pb. The B ion contributes by shifting closer to Ox, but for the B–Ox pair the elastic effect is substantially smaller than it is in the case of A–Oy and A–Oz. We also note that when these results are compared with the analogous ones for the smaller $`e_{33}`$ coefficient (not shown here – see Refs.  and  for related PZT results), one finds that it is the elastic effect that makes $`e_{15}`$ particularly large.

The above results on $`e_{15}`$ have revealed much about the chemical origin of the piezoelectric response of these systems. However, the energetics of the piezoelectric response, as well as the actual technological interest of a material, is best characterized by the direct piezoelectric coefficients $`d_{ij}=dP_i/d\sigma _j`$ reflecting the polarization that develops in direction $`i`$ by application of a stress with Voigt component $`j`$ (or equivalently, the strain $`j`$ induced by an electric field along $`i`$). In particular, $`d_{15}`$ measures the ease of polarization rotation associated with the large responses at the MPB. Technical difficulties have traditionally hampered the direct calculation of $`d_{ij}`$ coefficients from first principles. In the present case, however, it is not hard to show that $`d_{15}`$ is just given by $`e_{15}/\stackrel{~}{C}_5`$, where $`\stackrel{~}{C}_5`$ is the shear elastic constant associated with $`\eta _5`$ when the atomic positions are allowed to relax in response to the strain (i.e., the dressed elastic constant) at fixed electric field. We calculated $`d_{15}=`$168 pC/N (with $`\stackrel{~}{C}_5`$=0.62 a.u.) for BS-PTx=0.7 and $`d_{15}=`$196 pC/N (with $`\stackrel{~}{C}_5`$=0.56 a.u.) for PZTx=0.6. These results agree in magnitude with the piezoelectric coefficients of ceramic samples of BS-PT and PZT measured at room temperature and at the MPB, as well as with the semiempirical calculation of Wu and Krakauer for an ordered PZT supercell.

Thus, our calculations clearly provide evidence that BS-PT can be expected to be a piezoelectric of comparable quality to PZT. This conclusion, as well as the leading role of Bi and/or Pb in the response, is further supported by calculations of $`e_{33}`$ in the tetragonal phase of the above-mentioned 10-atom supercells (i.e., by calculations in which the VCA was not used).

## VI Summary and conclusions

Our first-principles study of BS-PT indicates that this material displays very large structural distortions and polarizations at the MPB. In particular, we obtain $`c/a`$ between 1.05 and 1.08 and $`P_{\mathrm{tet}}0.9`$ C/m2. We also find that the piezoelectric response of BS-PT near the MPB is comparable to that of PZT.

Our calculations show that the large polarization and piezoelectric responses in BS-PT are mainly related with the onset of Bi/Pb-6$`p`$–O-2$`p`$ hybridization, a mechanism that is enhanced upon substitution of Pb by Bi, since Bi is a more covalent atom. Our results also provide evidence of the dominant role of Pb in the piezoelectric response of PZT, in agreement with recent experimental studies.  In both BS-PT and PZT, this predominance of the A-site ions is related to a large elastic effect, i.e., a large displacement of the ions in response to strain.

Our results are in reasonably good agreement with experiment. The quantitative discrepancies between theory and experiment (in particular, for the values of $`c/a`$ and the polarizations at the MPB) suggest that the intrinsic piezoelectric properties of BS-PT alloys may be even better than those measured experimentally to date. We thus hope that our results will stimulate further experimental work on this and related materials systems.

###### Acknowledgements.

J. I. thanks J. B. Neaton for many stimulating discussions. This work was supported by ONR Grant No. N0014-97-1-0048 and by the Center for Piezoelectrics by Design (CPD) under ONR Grant N00014-01-1-0365. Computational facilities for the work were also provided by the CPD.