text:
# Space-VLBI observations of OH maser OH34.26+0.15: low interstellar scattering

## 1 Introduction

OH masers were the first to be discovered as a new class of astronomical phenomenona exhibiting intense, narrow band, polarized, variable emission in interstellar molecular transitions \[Weaver et al. 1965\]. The source of this emission was shown to be small, compact clumps of neutral gas, at the periphery of compact H ii regions created by newly born massive O-stars. Another class of OH masers is related to evoled low- and medium-mass stars, and is not considered here. High angular resolution study of OH masers showed that they consist of a number of bright compact maser spots separated by several arcseconds. The spots themselves are from 2 to 70 milliarcseconds in extent, and in some masers are barely resolved, even with the highest angular resolution available from ground based VLBI arrays \[Moran et al. 1968\].

The goal of the high angular resolution study of OH masers is to determine relative position of the spots, to measure the intensity, and to map the spots. Determination of the size and the shape of the maser spots would be very helpful for understanding the maser emission mechanism and its limiting brightness temperature. The shape of the maser spots may be indicative of the type of physical phenomena responsible for the origin of the maser emission.

Several possible sites for the origin of maser emission; such as the shocks at the border of H ii regions, or the interaction region between a molecular outflow and ambient molecular clouds, as well as protoplanetary accretion disks around young star, have been suggested. The maser spot shape and size are important properties of the maser emission, provided they are intrinsic to the source, and not caused by a propagation effect such as interstellar scattering.

There is a widespread opinion that OH masers are heavily scattered by interstellar turbulence, and that observed images of maser spots are scattering broadened images of essentially point-like sources. This view is based on data obtained from the study of scattering of pulsars and continuum sources. The importance of scattering for OH masers is due to their location within the thin galactic disk where the scattering is most severe. The first observational evidence for the scattering of OH masers obtained by direct VLBI measurements came from Kemball, Diamond and Mantovani \[Kemball et al. 1988\]. They measured the size of the strongest spectral features in 16 galactic masers using a single-baseline interferometer with a fringe spacing of 5 milliarcsec. Only three masers out of 16 observed have shown compact details, less than 2.5 milliarcsec in size. The rest of the masers were completely resolved implying angular size larger than 5 milliarcsec. Kemball et al. \[Kemball et al. 1988\] concluded that their results were consistent with the interstellar scattering angle of about 20 milliarcsec derived from pulsar scattering data, and with the clumpy distribution of electron density fluctuations which are responsible for the large spread of measured angular sizes. Hansen et al. \[Hansen et al. 1992\] made a VLBI survey of interstellar broadening of 20 galactic OH masers, using a 7-station VLBI array. They found a large range of angular sizes from 0.8$`\pm `$0.5 milliarcsec to 74$`\pm `$7 milliarcsec in their sample of OH masers. They fitted the angular size $`\theta `$ versus distance $`D`$ dependence by a relation $`\theta D^{0.65}`$ which is close to the expected scaling with distance for homogeneous turbulence with a Kolmogorov power law spectrum. They concluded that the measured sizes result from interstellar scattering which is strongly non-uniform.

Among the OH masers measured in these VLBI experiments, several are unresolved to the limit of angular resolution. Slysh et al. \[Slysh et al. 1996\] used a 3-station VLBI network to measure the angular size of three such compact OH masers. With a fringe separation of 4.2 milliarcsec and high signal to noise ratio data the angular size of several spectral features in the two main line OH transitions were measured simultaneously. The measured angular size was in the range from 1.4$`\pm `$0.4 milliarcsec to 4.3$`\pm `$0.1 milliarcsec. It was found that the angular size of the spectral features with the same radial velocity, presumably originating in the same region, was larger in the 1667 MHz line than in the 1665 MHz line. This cannot be the result of scattering, since the scattering size should be virtually the same at these two close frequencies (the scaling between 1665 MHz and 1667 MHz is negligible!). It was concluded that the measured angular size at 1667 MHz 4.3$`\pm `$0.1 milliarcsec for OH34.26$`+`$0.15 and 3.5$`\pm `$0.2 milliarcsec for W48 is intrinsic to the sources. These two masers together with OH45.47+0.13 exhibit an angular size for spectral features which is much less than the average broadening estimated from pulsar measurements. The conclusion of that study was that the distribution of the scattering material in the Galaxy is patchy, and that these three masers are located in a direction of low interstellar scattering. For this reason these masers were selected for space-VLBI measurements using the Japanese satellite HALCA and a ground telescope array. In this paper we report the first results of the space-VLBI observations of OH maser OH34.26$`+`$0.15.

## 2 Observations and data reduction

The maser source OH34.26$`+`$0.15 was observed as a part of the Key Science Program of the Japanese Satellite HALCA on March 24-25, 1998, by the space-ground very long baseline interferometer. The satellite radio telescope has an 8-m diameter deployable parabolic mirror and uncooled L-band receiver, orbiting the Earth with a 6 hour period, apogee 21000 km and perigee 560 km \[Hirabayashi et al. 1998\]. The ground radio telescope array in this particular session included the 70-m DSN telescope Tidbinbilla in Australia (Tid), the phased 6$`\times `$22m Australia Telescope (AT), the 22-m Mopra Telescope in Australia (Mop), the Shanghai 25-m Telescope in China (Sh), the 64-m Telescope in Usuda Japan (Usu), and the 64-m Telescope in Bear Lakes (BL) near Moscow (Russia). It was the first time that the Bear Lakes Telescope participated in space-VLBI observations. Table 1 is a list of telescopes with their parameters: system noise temperature Tsys, gain and rms noise in HALCA to ground telescope baselines after 5 minutes of integration. The uv-coverage of the ground telescope array is shown on Fig. 1a, and the uv-coverage of the full array including the space radio telescope is shown on Fig. 1b. It is evident that the addition of the space radio telescope not only increases the uv-coverage by a factor of two due to the high orbit of the space radio telescope, but also significantly fills the uv-plane coverage at lower resolution due to the orbital motion of satellite. For the ground array the synthesized beam was 5.8$`\times `$2.1 milliarcsec, and for the space-ground array the beam was 1.8$`\times `$1.0 milliarcsec.

Left circular polarization data over a 16 MHz bandwidth centered on the 1665 MHz and 1667 MHz OH main lines were recorded using S2 recorders at each station. The data were correlated using the Canadian S2 Space-VLBI correlator111The Canadian S2 Space-VLBI correlator at the National Research Centre of Canada, Dominion Radio Astrophysical Observatory is operated by the Canadian Space Agency. \[Carlson et al. 1999\]. The 16 MHz data were filtered using a FIR bandpass filter and then re-sampled before correlation in order to attain a spectral resolution of 488 Hz channel-1 (0.088 km s-1) over a bandpass of 500 kHz. Amplitude calibration was achieved using the “gain” and Tsys data (Table 1) obtained from ground observatories and from the VSOP operations group. The post-correlation data reduction was performed using the aips package of NRAO.

Phase jitter was discovered at the Bear Lakes telescope which caused splitting of the fringe spectrum into several components. It was found that the phase of the Bear Lakes telescope receiver experienced quasi-periodical excursions with a period of about 200 seconds. The phase jitter was removed by means of phase self-calibration of the ground array with a short solution interval of 30 seconds, which is less than the instrumental, tropospheric and ionospheric phase variation period. The multiple spectral components disappeared, and the central peak increased by a factor of two.

The strongest and most compact spectral feature at the radial velocity 56 km s-1 in the 1665 MHz spectrum was taken as a reference in the self-calibration process. Its position relative to the correlating position was determined by the fringe rate method, using the ground array. The measured offset is 0$`\stackrel{}{.}`$23 in right ascension and 1$`\stackrel{}{.}`$4 in declination, with an accuracy of 0$`\stackrel{}{.}`$1, and our best position for this feature is $`\alpha =18^h53^m18\stackrel{s}{.}685`$, $`\delta =01\mathrm{°}14\mathrm{}59\stackrel{}{.}9`$ (J2000.0). The position of the reference feature A determined here differs from the measured with the VLA by Gaume & Mutel \[Gaume & Mutel 1987\] by about 0$`\stackrel{}{.}`$15 in both coordinates, which is within the combined error limits of the two measurements. After applying this correction the residual fringe rate an all base lines including those to HALCA was found to be close to zero, and the phase variations became slow enough to allow integration time up to 5 minutes.

The delay errors for the ground array telescopes were measured to be less than 300 nanoseconds from observation of the continuum calibration source 1730$``$130, which was observed during the time interval when HALCA was not observing. The maximum residual delay associated with the baselines to HALCA obtained during correlation of the experiment using the Penticton Correlator, which assumed an accurate HALCA orbit and other known apriori information, was less than 300 nanoseconds at both 1665 MHz and 1667 MHz. Thus, any additional delay errors after correlation are unlikely to be larger than 300 nanoseconds which will produce a negligable phase difference over the 5.3 km s-1 spread across maser lines.

The relative positions of the 1665 MHz line maser spots with respect to feature A were determined in two steps. First, an approximate position relative to the reference feature was measured by the fringe rate method with the ground array. These measurements provided positions with an accuracy of about 50 milliarcseconds, with relative displacements of the features reaching 2000 milliarcseconds. Then maps of the separate features were constructed centered on these approximate positions, with the ground only or space-ground array when the signal to noise ratio was large enough. The features mapped by the space-ground array are marked by stars in Table 2. The imaging and cleaning were done with the aips task “imagr”, and the features were analyzed with the aips task “sad” from which the flux density, position and deconvolved size (or size limit) were obtained for each feature. The relative positional accuracy is about 1 milliarcsecond. Result of this analysis is given in Table 2, where for each Gaussian component are given: radial velocity, relative position, flux density, angular size of the major and minor axis of deconvolved Gaussian components as well as position angle of the major axis.

The 1667 MHz OH line data were reduced in the following manner. The 58.70 km s-1 feature in F is the strongest feature and that was used as the phase reference. Images of the other features in F and feature E (after using fringe-rate mapping to determine its approximate position), were then made by imaging and cleaning as for the 1665 MHz line. In order to tie the 1665 MHz and 1667 MHz images to one coordinate system, we assumed that feature E at 1667 MHz was at the position ($``$1350,$``$3530) milliarcseconds displaced from feature A at 1665 MHz, following Gaume & Mutel \[Gaume & Mutel 1987\]. This registration should be accurate to about 50 milliarcseconds. The near coincidence of 1665 MHz features B and G with that of F at 1667 MHz is independent evidence that this registration is accurate.

## 3 Results

The correlated LCP spectra of OH34.26+0.15 on the short Australia Telescope – Mopra baseline on which the source is unresolved, in both main line OH transitions is shown in Fig. 3. Five spectral features A–D and G are present in the 1665 MHz spectrum, and two features E and F in the 1667 MHz spectrum; their relative position is shown in Fig. 3. Fig. 4a shows a map of feature A obtained with the ground array. The feature is unresolved with a synthesized beam of 5.8$`\times `$2.1 milliarcsec (lower left corner). Fig. 4b shows a map of the same feature A obtained with the space-ground array having a synthesized beam of 1.8$`\times `$1.0 milliarcsec (lower left corner). The feature is partly resolved, and this is evident from the different position angles of the beam and the image. A deconvolved fitted Gaussian to this image is an ellipse at a position angle of 70°, with the major axis 2.0 milliarcsec and the minor axis equal to or less than 0.5 milliarcsec. Inspection of Table 2 shows that all components are partly resolved with the beam of the space-ground or ground only array, but most of the components may have structure and consist of several smaller components.

Other 1665 MHz spectral features are offset from the reference feature by up to 1$`\stackrel{}{.}`$9, their position is shown in Fig. 3. Feature B consists of four separate closely spaced spots as can be seen on the map (Fig. 5a) and in Table 2, where radial velocity, relative position, flux, and angular size of all spectral features are given. Feature F (1667 MHz) coincides with feature B within the measurement errors, and its shape is similar to the shape of feature B (Fig. 5). These two features may come from the same maser condensation, a supposition that can be tested with higher accuracy relative 1665/1667 position measurements. If this is the case the radial velocities of features B and F must coincide, and the observed velocity difference of 0.5 km s-1 can be attributed to a difference between Zeeman splittings of 1665 MHz and 1667 MHz lines because of different g-factors of two transitions. If the observed features B and F are $`\sigma `$-components of Zeeman pairs (both are left circular polarized), one can calculate that the velocity difference of 0.5 km s-1 corresponds to the magnetic field of B=4.2 milligauss, which is a value commonly found in OH masers \[Garcia-Barreto et al. 1988\].

## 4 Discussion

### 4.1 Maser spot distribution

OH34.26$`+`$0.15 is located at a distance of 3.8 kpc \[Reid & Ho 1985\] near a cometary H ii region and near two ultra-compact H ii regions. Maser spots A, C and D coincide with the northern ultra-compact H ii region G34.3$`+`$0.2B, while spots B, E, F and G coincide with the cometary H ii region G34.3+0.2C (see Fig. 7 in Gaume, Fey and Claussen 1994). Our map of OH components (Fig. 3) is consistent with the VLA map of the same source obtained by Gaume & Mutel \[Gaume & Mutel 1987\] in 1985. All components A–E have their counterparts on the Gaume & Mutel \[Gaume & Mutel 1987\] map, except component D. On the other hand, some weaker components from the map of Gaume & Mutel \[Gaume & Mutel 1987\] are not present in our map. The splitting of components B and F into four subcomponents is not evident on the Gaume & Mutel \[Gaume & Mutel 1987\] map because of the much lower resolution of the VLA. The relative position of the components on both maps is the same within the VLA position error of about 50 milliarcsec, while the relative positional errors in our map are about 1 milliarcsec. The mean difference in position is 5.5$`\pm `$25 milliarcsec.

### 4.2 Resolution of the feature A

The most compact component is feature A, which is known to be left hand circularly polarized \[Caswell & Haynes 1983\]. It was resolved in one direction, and the angular size in this direction is 2.0 milliarcsec and less than 0.5 milliarcsec in the perpendicular direction. Feature A was fitted with an ellipse in position angle 70° with the ratio of major to minor axis greater than 4. This size is a result of the fitting and deconvolution with the beam of the image shown in Fig. 4b. As an almost independent check of the size of feature A we made a direct fit of a Gaussian model to uv-data using only fringe amplitude information. This procedure is independent in the sense that it does not use the self-calibration and global fringe fitting. The best fit to fringe amplitudes from three telescopes AT, Tid and HALCA was achieved with an elliptical Gaussian, having major and minor axis 2.6$`\pm `$0.5 and 0.3$`\pm `$0.3 milliarcsec and a position angle 61$`\stackrel{}{.}`$5$`\pm `$3$`\stackrel{}{.}`$2 which, within the errors, is the same as obtained from the image. From our data it is not possible to say whether maser spot A has a stripe-like shape, or is composed of a chain of point-like sources. To determine this, one needs higher angular resolution and higher signal to noise data in order to obtain maps with a higher dynamic range.

In this experiment we could not use the full angular resolution of the space-ground array because of the limitations set by the signal to noise ratio of the data. While the orbit of HALCA provided uv-distances up to 135 M$`\lambda `$, the fringes were lost in the noise every time the uv-distance from a ground telescope to the space telescope exceeded about 90 M$`\lambda `$. This can be seen in Fig. 6, where we show fringe amplitude and phase versus time for the baseline AT–HALCA. During the session the uv-distance first increased from 40 M$`\lambda `$ and reached a maximum of 113 M$`\lambda `$ at 11000 seconds, then decreased to 78 M$`\lambda `$ at 16000 seconds. It is evident from Fig. 6 that the fringe amplitude has decreased from about 12 Jy to less than 4 Jy which is below the noise level. This can be seen from the phase behaviour in the upper panel of Fig. 6: the phase varies smoothly from 0 to about 6000 seconds (about 90 M$`\lambda `$), and then becomes chaotic. The smooth phase variation resumes during the second portion beginning at about 23000 seconds, up to the end of this session. During the second portion the uv-distance was between 51 and 64 M$`\lambda `$ and the fringe amplitude was above the noise. This example shows that the low signal to noise ratio prevents us from using the full angular resolution available. If the signal to noise ratio were increased by a factor of 3 it would be possible to implement full angular resolution imaging. The main cause of the low signal to noise ratio is the small size of the space radio telescope. An increase in the signal to noise ratio by a factor of three would require an increase in the space telescope diameter by the same factor, or a combination of larger diameter and lower system noise temperature.

The upper limit of the angular size of feature A cited above corresponds to the lower limit of its brightness temperature of 6$`\times 10^{12}`$ K. Other features have somewhat larger sizes and lower brightness temperature limits, although it cannot be excluded that they consist of several unresolved subcomponents. This value exeeds the brightness temperature predicted by some OH maser models. For example a model by Pavlakis & Kylafis \[Pavlakis & Kylafis 1996\] predicts brightness temperatures about $`10^{11}`$ K, or at least an order of magnitude lower than observed. On the other hand there is an absolute upper bound to the brightness temperature imposed by the stimulated transition line broadening (or resonant Stark effect \[Slysh 1973\]). The stimulated emission rate $`R`$ is

$$R=\frac{\mathrm{\Omega }}{4\pi }\frac{kT_B}{h\nu }A=889\frac{\mathrm{\Omega }}{4\pi }\times T_{B12}\mathrm{s}^1$$
(1)
where $`\frac{\mathrm{\Omega }}{4\pi }`$ is maser emission solid angle, $`A`$ is the spontaneous emission rate ($`7.11\times 10^{11}`$ s-1 for 1665 MHz OH transition), and $`T_{B12}`$ is the brightness temperature in units of 1012 K. For feature A the stimulated emission rate $`R=\frac{\mathrm{\Omega }}{4\pi }5.4\times 10^3`$ s-1, and is comparable to the observed line width which is 0.33 km s-1, or 1.6$`\times 10^3`$ Hz. For an isotropic maser $`\mathrm{\Omega }=4\pi `$, and the contribution of the stimulated emission to the line width is $`\frac{R}{2\pi }=8.5\times 10^2`$ Hz, or about a half of the linewidth. This means that the brightness temperature can not be much higher than the observed upper limit $`6\times 10^{12}`$ K, unless the maser emission is highly directed. For example, if $`\frac{\mathrm{\Omega }}{4\pi }=10^2`$, the brightness temperature can be as high as $`6\times 10^{14}`$ K. Since the maser lines do not show any measurable broadening the maser is probably anisotropic.

Elongated shapes of OH masers features are also seen in other sources. Baudry & Diamond \[Baudry & Diamond 1998\] found a filamentary and arc-like structure in the OH maser W3(OH) at 13.44 GHz. They relate this structure to a shocked environment of the ultra-compact H ii region. One can also speculate that the elongated filaments were formed as a result of an anisotropic compression of the molecular gas into dense condensations in the presence of the magnetic field which is observed in OH masers through Zeeman splitting and polarization of emission lines.

### 4.3 Interstellar scattering

OH masers connected with regions of formation of massive stars are located in the plane of the Galaxy, and are subject to strong interstellar scattering. As a result of scattering, the image of small sources is broadened to an extent that depends on the level of turbulence that produces scattering electron density inhomogeneities. The observed small angular size of maser spots in OH34.26$`+`$0.15 puts a strong upper limit on the turbulence level which is usually characterized by a constant $`C_N^2`$ in the turbulent electron density power spectrum. The scatter-broadened angular size is related to the so-called scattering measure

$$SM=\underset{0}{\overset{L}{}}C_N^2𝑑s\text{,}$$
(2)
where the integral is over distance $`s`$ up to the source at distance $`L`$. For galactic sources the scattered angular diameter $`\theta _s`$ and $`SM`$ are related by the expression \[Taylor & Cordes 1993\]

$$SM=\left(\frac{\theta _s}{71mas}\right)^{\frac{5}{3}}\nu _{GHz}^{\frac{11}{3}}kpcmeters^{\frac{20}{3}}\text{.}$$
(3)
If the measured upper limit of the size of the minor axis of component A of 0.5 milliarcsec is due to scattering, and $`\nu =1.6654`$ GHz, then from equation (3) one has

$$SM=1.68\times 10^3kpcmeters^{\frac{20}{3}}\text{.}$$
(4)
At a distance of 3.8 kpc this gives $`C_N^2=4.4\times 10^4meters^{\frac{20}{3}}`$. This is almost an order of magnitude less than the average found from pulsar scattering measurements, where $`C_N^2=10^{2.5}m^{\frac{20}{3}}`$ (fig. 13 in Cordes, Weisberg & Boriakoff 1985). This means that scattering towards OH34.26$`+`$0.15 is much less than the average in the Galaxy.

There exist measurements for interstellar scattering from a nearby extragalactic source 1849+005, which is only 0$`\stackrel{}{.}`$76 from OH34.26$`+`$0.15 \[Fey, Spangler & Cordes 1991\]. For this source the angular size measured at the frequency 1465 MHz is 560 milliarcsec. Scattering was found to scale with frequency as $`\nu ^{1.65}`$ in the range from 333 MHz to 4860 MHz. From pulse broadening, Frail & Clifton \[Frail & Clifton 1989\] estimated the angular size of a highly scattered pulsar PSR 1849$`+`$00 (which is within 10′ of 1849$`+`$005) as 500 milliarcsec at 1 GHz. If both values are scaled to the OH-line frequency of 1665 MHz, one has 450 milliarcsec from 1849+005 and 216 milliarcsec from the pulsar. These are two to three orders of magnitude larger than the size of the maser spots in OH34.26$`+`$0.15. We note that 1849$`+`$005 is an extragalactic source, and PSR 1849$`+`$00 is a very distant pulsar at a distance of 14.5 kpc \[Frail & Clifton 1989\]. Therefore, the ray path from the sources to the Sun traverses large distances in the Galaxy, with the closest approach to the Galactic Centre of about 4 kpc. The maser OH34.26$`+`$0.15 is on the same path, but is much closer to the Sun, at 3.8 kpc. If most of the scattering material lies beyond that distance, presumably near the point where the path is closest to the Galactic Centre (at 6.6 kpc from the Sun, assuming the distance to the Galactic Centre from the Sun 8 kpc), then the OH maser emission will not be scattered so much. One could argue that the large discrepancy of scattering parameters between OH maser and continuum sources (like pulsars) is due to small scale variations in the distribution of the scattering material. The angular separation between the OH maser and 1849$`+`$005 is 0$`\stackrel{}{.}`$76 and corresponds to 50 pc. This small scale is inconsistent with an almost equally small scattering of two other OH masers in this region OH35.20$``$1.76 and OH45.47$`+`$0.13 \[Slysh et al. 1996\] which are separated by about 10°, or 650 pc. Based on observations of three low scattered OH masers in the same region, we suggest that there is a large scale deviation in the distribution of the scattering material in the Galaxy.

## 5 Summary and conclusions

We report on the first space-VLBI observations of OH masers. It was also the first time that the 64-m telescope at Bear Lakes near Moscow participated in a space-VLBI experiment. We constructed a map of the maser OH34.26$`+`$0.15 with angular resolution better than 1 milliarcsec. The maser consists of small bright spots distributed over several arcseconds, the size of the spots in some cases is less than 0.5 milliarcsec, and the brightness temperature is larger than 6$`\times 10^{12}`$ K. This OH maser source experiences an exceptionally small interstellar scattering, and we propose that there is a strong large scale variation in the distribution of the scattering material in the Galaxy.

## 6 Acknowledgements

We gratefully acknowledge the VSOP project, which is led by the Japanese Institute of Space and Astronomical Science in cooperation with many organizations and radio telescopes around the world. The National Radio Astronomy Observatory is a facility of the National Science Foundation, operated under a cooperative agreement by Associated Universities, Inc.


format: text

text:
# Precision Cosmology and the Density of Baryons in the Universe

\[

## Abstract

Big-bang Nucleosynthesis (BBN) and Cosmic Microwave Background (CMB) anisotropy measurements give independent, accurate measurements of the baryon density and can test the framework of the standard cosmology. Early CMB data are consistent with the longstanding conclusion from BBN that baryons constitute a small fraction of matter in the Universe, but may indicate a slightly higher value for the baryon density. We clarify precisely what the two methods determine, and point out that differing values for the baryon density can indicate either an inconsistency or physics beyond the standard models of cosmology and particle physics. We discuss other signatures of the new physics in CMB anisotropy.

\]

Introduction. Just a decade ago the phrase “precision cosmology” would have been an oxymoron. The COBE FIRAS determination of the temperature of the Cosmic Microwave Background (CMB) to four significant figures, $`T_0=2.725\pm 0.001\mathrm{K}`$  (we quote all errors at $`1\sigma `$) should dispel such thoughts. Cosmologists now foresee a precision era where a flood of high-quality data, from measurements of CMB anisotropy to large-scale structure, pin down cosmological parameters to percent-level precision, decisively testing theories of the early Universe and probing physics at energy scales beyond those accessible in accelerator experiments .

Some of the data can also test the consistency of the standard cosmology. In particular, the longstanding big-bang nucleosynthesis determination of the baryon density, will be checked by CMB anisotropy data at one percent or better . The physics underlying the two measurements could hardly be more different: the baryon density at a time of one second determines how complete the conversion of neutrons and protons to tightly-bound 4He nuclei is, while the baryon density 400,000 years later (at the time of last scattering) determines the amplitude of gravity-driven acoustic oscillations in the baryon-photon fluid.

The determination of the primeval deuterium abundance in a number of high-redshift ($`z24`$) hydrogen clouds  coupled with refined predictions of the standard theory of BBN has led to a determination of the baryon density to an accuracy of about 5%, $`\mathrm{\Omega }_Bh^2=0.019\pm 0.00095`$ . Very recently, the BOOMERanG  and MAXIMA  balloon-borne CMB anisotropy experiments have mapped CMB anisotropy with sufficient angular resolution to make the first CMB determinations of the baryon density, $`\mathrm{\Omega }_Bh^2=0.032_{0.004}^{+0.005}`$ . While the CMB result has a much larger uncertainty (and depends upon the imposed priors and parameters that are allowed to vary ), additional data and new CMB measurements should soon narrow the gap between the two methods.

The early results are encouraging for the consistency of the standard cosmology: the CMB value for the baryon density agrees with the BBN determination within about $`2\sigma `$ and lies far from dynamical determinations of the total mass density, $`\mathrm{\Omega }_Mh^2=0.2\pm 0.04`$ (inferred from $`\mathrm{\Omega }_M=0.35\pm 0.07`$  and $`h=0.7\pm 0.07`$ ). Thus, the CMB measurement strongly supports the case for nonbaryonic dark matter, whose linchpin for twenty years has been the discrepancy between the BBN value for $`\mathrm{\Omega }_B`$ and dynamical measurements of $`\mathrm{\Omega }_M`$ .

The purpose of our Letter is to clarify exactly what is determined, with regard to the baryon density, by BBN and CMB anisotropy and to point out that the two determinations need not agree to be consistent. A discrepancy could very well be the signal of new physics. Here and throughout $`\mathrm{\Omega }_B`$ denotes the fraction of critical density in baryons today and $`h=H_0/100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. The physical baryon density today, $`\rho _B=1.88\times 10^{29}(\mathrm{\Omega }_Bh^2)\mathrm{g}/\mathrm{cm}^3`$.

Comparing apples to apples. In the standard theory of BBN (i.e., isotropic and homogeneous Universe with only the known particle species and the assumption that all three neutrino species are light with mass $`1\mathrm{MeV}`$, and negligible chemical potentials), the yields of BBN depend only upon the baryon-to-photon ratio $`\eta `$, the neutron mean lifetime and eleven key nuclear cross sections . The uncertainties due to the neutron mean lifetime and nuclear data have recently been re-evaluated and significantly reduced . Based upon these predictions and uncertainties, the Burles-Tytler primeval deuterium measurement, $`(\mathrm{D}/\mathrm{H})_P=(3.4\pm 0.3)\times 10^5`$, implies a baryon-to-photon ratio $`\eta _{\mathrm{BBN}}=(5.1\pm 0.2)\times 10^{10}`$ .

In order to infer the present density of baryons one must convert $`\eta _{\mathrm{BBN}}`$ to a baryon density at the time of BBN by multiplying by the photon number density, $`n_\gamma (\mathrm{BBN})=2\zeta (3)T_{\mathrm{BBN}}^3/\pi ^2`$, and the mean mass per baryon ($`\overline{m}`$) and then reduce that density by the volume increase of the Universe since,

$$\rho _B(\mathrm{today})=2\zeta (3)\overline{m}\eta _{\mathrm{BBN}}R_{\mathrm{BBN}}^3T_{\mathrm{BBN}}^3/\pi ^2R_0^3,$$
(1)
where $`R(t)`$ is the cosmic scale factor and baryon-number conservation since BBN has been explicitly assumed.

Because we do not know the value of scale factor at the time of BBN a priori, we cannot proceed without further assumptions. The standard assumption is adiabaticity: the constancy of the electromagnetic entropy per unit comoving volume, which is proportional to $`(RT)^3`$ , since BBN. Allowing for the possibility that the entropy per unit comoving volume has changed it follows that

$$\rho _B(\mathrm{today})=\frac{S_{\mathrm{EM}}(\mathrm{BBN})}{S_{\mathrm{EM}}(\mathrm{today})}\eta _{\mathrm{BBN}}\overline{m}n_\gamma (\mathrm{today}).$$
(2)
The number density of photons today is related to the present temperature of the CMB, $`n_\gamma =(410.5\pm 0.5)\mathrm{cm}^3`$; dividing by the critical density we arrive at the key equation,

$$\mathrm{\Omega }_Bh^2=(0.019\pm 0.00095)\frac{S_{\mathrm{EM}}(\mathrm{BBN})}{S_{\mathrm{EM}}(\mathrm{today})}.$$
(3)
Entropy production after BBN (e.g., by the out-of-equilibrium decay of a massive particle) can increase $`S_{\mathrm{EM}}`$, which would diminish the BBN prediction for the baryon density. It is also possible to reduce $`S_{\mathrm{EM}}`$ (at the expense of more exotic physics ), increasing the BBN prediction for the present baryon density.

The signature of the baryon density in the CMB involves the heights of the “acoustic peaks” in the angular power spectrum. Around the time of last scattering the photon – baryon fluid is undergoing gravity driven acoustic oscillations; Fourier modes caught at maximum compression (odd peaks) or maximum rarefaction (even peaks) produce the highest amplitude temperature fluctuations on the sky, leading to a series of acoustic peaks in the angular power spectrum (see Fig. 1) .

The ratio of the heights of the odd and even peaks increases with baryon density; all other cosmological parameters tend to move the heights of the peaks in unison. Thus, determining the baryon density does not suffer from the cosmic degeneracies that affect other parameters, and an ultimate precision of better than one percent can be expected .

Converting the baryon density at last scattering (redshift $`z_{\mathrm{LS}}1100`$) to the present involves only a factor of $`(1+z_{\mathrm{LS}})^3`$ . While a number of analyses of the BOOMERanG and MAXIMA data have been carried out , the current state of affairs is probably fairly represented by the joint analysis of the two teams ,

$$\mathrm{\Omega }_Bh^2=0.032_{0.004}^{+0.005}.$$
(4)
Changing entropy and its CMB signature. The entropy increase due to the out-of-equilibrium decay (i.e., when $`Tm`$) of a massive particle relic is given in terms of the particle’s lifetime and mass :

$$rS_f/S_i3m_{\mathrm{Pl}}^{1/2}Ym\tau ^{1/2},$$
(5)
where $`Yn/s`$ is the pre-decay abundance (number density per unit comoving entropy density – for a neutrino species, $`Y=0.04`$). We have assumed that the decay products thermalize into photons and increase the EM entropy; to avoid distorting the nearly perfect black-body spectrum of the CMB this must occur before about $`10^6`$sec. In this case the baryon density today is smaller (than $`\mathrm{\Omega }_Bh^2=0.019`$) by a factor of $`r`$.

The increase in entropy has a CMB signature which involves the fact that the energy density in relativistic particles at the time of last scattering is smaller than in the standard scenario. This is because EM entropy production increases the photon-to-neutrino temperature ratio by a factor of $`r^{1/3}`$ over the standard value, $`T_\gamma /T_\nu =r^{1/3}(11/4)^{1/3}`$, and thereby decreases the energy in neutrinos at last scattering (which occurs at a fixed photon temperature, $`T_{LS}0.3\mathrm{eV}`$). This decrease in relativistic energy density can be expressed in terms of a (lower) equivalent number of standard neutrino species (see Fig. 2):

$$N_\nu ^{\mathrm{equiv}}=3/r^{4/3}.$$
(6)
More or less energy in relativistic particles at the time of last scattering affects the spectrum of CMB anisotropy by changing the expansion rate and the rate at which the gravitational potentials associated with density perturbations decay. Less relativistic energy increases the sound horizon at last scattering, thereby shifting the acoustic peaks to larger angular scales (smaller $`\mathrm{}`$). Less relativistic energy also depresses the power around the first acoustic peak by diminishing the integrated Sachs-Wolfe effect , and decreases the damping of the higher-order acoustic peaks (because the peaks have shifted to larger scales relative to the damping length) as can be seen in Fig. 1.

Two points deserve mention here: (1) smaller $`N_\nu ^{\mathrm{equiv}}`$ also changes the power spectrum of matter inhomogeneity by shifting the epoch of matter – radiation equality to an earlier time. (2) Only a change in the EM entropy affects the predicted BBN baryon density; decays without EM decay products (e.g., decay to three neutrinos) do not affect the BBN prediction. If the massive-particle decays produce both EM and non-EM entropy, the EM decay products determine $`r`$, and the non-EM decay products increase the energy density in relativistic particles and lead to an additional term in Eq. (6).

The value of $`N_\nu ^{\mathrm{equiv}}`$ will ultimately be determined from CMB anisotropy to a few percent . This provides a cross check on this explanation: if the present baryon density inferred from BBN exceeds that inferred from CMB by a factor of $`\xi `$, then $`N_\nu ^{\mathrm{equiv}}`$ should be $`3/\xi ^{4/3}`$. Of course, the current data is inconsistent with any post-BBN entropy production, and seems to require $`r<1`$, i.e., entropy reduction.

Entropy reduction. By invoking more exotic physics it is also possible to reduce the EM entropy. In order to do so without violating the 2nd law of thermodynamics, photons must come into thermal contact with a cooler reservoir (denoted as the X sector) after BBN and before last scattering. Such an idea was previous discussed in a different context , and we will briefly discuss its possible relevance here.

The basic idea is simple; owing to the energy dependence of the cross section for X-sector particles to interact with photons (and other familiar particles), the X sector is decoupled (interaction rate per particle $`\mathrm{\Gamma }`$ less than the expansion rate $`H`$) at high temperatures. (This can be achieved provided that $`\sigma v1/T^n`$ with $`n>1`$.) At a temperature $`T_c`$, between the epoch of BBN and the present, the interactions become rapid enough to quickly establish thermal contact and $`T_X=T`$, thereby draining entropy from the photons (assuming that the X sector was cooler).

The authors of Ref.  have shown that a self-consistent model for the X sector can be constructed, and further, that it is consistent with astrophysical considerations that constrain the interaction of unseen particles with photons (e.g., emission of X sector particles through plasmon processes in red giants and supernovae). With that in mind, let us proceed.

For simplicity, assume that the transfer of energy to X particles proceeds quickly, by thermally populating massless degrees of freedom in the X sector with total statistical weight $`g_X`$. The decrease in the EM entropy, $`r=S_{\mathrm{EM}}(\mathrm{today})/S_{\mathrm{EM}}(\mathrm{BBN})`$, follows from energy conservation:

$$r=[2/(2+g_X)]^{3/4}.$$
(7)
The present apparent discrepancy between the BBN and CMB determinations could be explained with $`r2/3`$ ($`g_X1.5`$).

Provided $`T_cT_{LS}`$, this scenario also leads to an increase in energy density in relativistic particles at the time of last scattering and a signature in CMB anisotropy as discussed above. The energy increase arises from: (1) the energy density in X-sector particles; and (2) the higher neutrino-to-photon temperature ratio, $`T_\nu /T=r^{1/3}(4/11)^{1/3}`$. Again, parameterizing this by the equivalent number of standard neutrino species, it follows that

$$N_\nu ^{\mathrm{equiv}}=\frac{4}{7}\left(\frac{4}{11}\right)^{4/3}g_X+3/r^{4/3},$$
(8)
where the two terms correspond to the two effects just mentioned.

We note that the second term is mandatory and robust – depressing the photon temperature necessarily raises the ratio of the neutrino to photon temperature – and is identical in form to the term that arises in the previous case where the entropy is increased, cf. Eq. (6). The first term is model dependent and would be different if the X sector did not reach the same temperature as the photons or if only the massive degrees of freedom were excited. Finally, for $`r=2/3`$, $`N_\nu ^{\mathrm{equiv}}=8.3`$, which as Fig. 1 shows has a dramatic effect on the spectrum of CMB anisotropy.

Concluding remarks. By means of very different physics fine-scale CMB anisotropy and BBN each have the potential to determine the mean baryon density to percent accuracy or better. Currently, the BBN determination has a precision of 5% and the CMB measurement 15%, and the two disagree at about the $`2\sigma `$ level.

As we have emphasized, a disagreement between the two determinations of the baryon density need not indicate inconsistency: if the BBN baryon density is larger (smaller) by some factor, this could be explained by an increase (decrease) in the EM entropy since the time of BBN by the same factor . If the two baryon densities agree, then one can limit any post BBN electromagnetic entropy change. As discussed, any entropy change also has a distinctive testable signature in CMB anisotropy (see Fig. 2).

We have assumed the standard theory of BBN; relaxing one of its assumptions (e.g., a decaying tau neutrino with a mass of $`𝒪(\mathrm{MeV})`$  or large neutrino chemical potentials ) can also change the baryon density inferred from the primeval deuterium abundance. Both possibilities have been discussed . These solutions work by speeding up the expansion around the end of BBN through additional energy density in relativistic particles so that more deuterium remains unburnt. Thus, both solutions also lead to additional relativistic energy at last scattering. It should be noted that the larger expansion rate also increases the helium abundance and hence these solutions require another variable to offset the effect on helium.

A positive chemical-potential for the electron neutrino can achieve this offset. In both the decaying tau neutrino and the large chemical potential scenarios, consistency with all light element abundances can be achieved. We, however, consider the decaying tau neutrino scenario less likely in light of recent results  from the SuperKamiokande collaboration which hint at a mass much less than an MeV for $`\nu _\tau `$. To make a comparison between post-BBN entropy change and non-standard BBN we fix our attention on the large chemical-potential scenario. Using $`\xi `$ to denote the ratio between the actual and the BBN baryon density, we have used the standard BBN code to derive the relation between $`\xi `$ and the minimum $`N_\nu ^{\mathrm{equiv}}`$ required, shown in Fig. 2. The $`\xi `$$`N_\nu ^{\mathrm{equiv}}`$ relation in the non-standard BBN case is distinct from that due to post-BBN entropy decrease.

For more than two decades big-bang nucleosynthesis has been both a critical test of the standard cosmology and a probe of particle physics (e.g., the limit to the number of light neutrino species) and cosmology (e.g., the baryon density). The cross comparison with the baryon density inferred from CMB anisotropy measurements has opened a new window for testing the cosmological framework and exploring physics beyond the standard model. Should the current discrepancy of about $`2\sigma `$ persist, the resolution may well involve non-standard BBN (neutrino chemical potentials) or new physics (entropy change due to new particles). If so, this would be an even more impressive achievement for cosmology than the BBN limit to the number of neutrino species which itself was ultimately confirmed by laboratory experiment.

###### Acknowledgements.

We acknowledge support by the DOE at Chicago, and by the DOE and NASA grant NAG 5-7092 at Fermilab.


text:
# Multiple BeppoSAX Observations of IC 4329A to Probe the Origin of the Compton Reflection Component in Seyfert 1 Galaxies

## 1. Introduction

Ginga, ASCA and recent BeppoSAX X-ray observations of Seyfert 1 galaxies strongly support a general observational framework that includes: ionized absorption by a warm absorber (WA), a steep intrinsic power-law continuum, a Compton reflection hump and associated neutral FeK$`\alpha `$ line (narrow and/or broad), and a high-energy cutoff ($`E_{\mathrm{cutoff}}`$). These observational results have raised the following questions:

* What is the location of the X-ray reprocessor? Is it the accretion disk, the broad-line region (BLR) and/or the molecular torus? Is the answer the same for all Seyfert 1s?
* Where is the location of the WA? Is it within the BLR and/or the molecular torus or is it external to them? Is there only one WA?
* Does $`E_{\mathrm{cutoff}}`$ differ from source-to-source (as currently shown by BeppoSAX, Piro 1999, Matt 2000) or is it variable with time and/or flux in any single object?

In order to address these questions, we have performed multiple (4) observations of the prototypical Seyfert 1 galaxy IC4329A to obtain time-resolved spectral constrains on the WA, reflection continuum, Fe Kα line and, possibly, $`E_{\mathrm{cutoff}}`$. This is potentially the best way (see also Weaver 2000) to establish whether these observational features are produced from material near the source (e.g. the accretion disk, in which case short lags, $`<`$1000 s, are expected between the reflection and continuum variations) or farther away (e.g. in the BLR or torus, in which case longer lags of at least $``$ 106-7 s are expected).

## 2. The Multiple BeppoSAX Observations

IC 4329A $`(z=0.014)`$ is a Seyfert 1 galaxy well studied in X-rays (Miyoshi et al. 1988; Piro, Yamauchi & Matsuoka 1990; Fiore et al. 1992; Cappi et al. 1996) and hard X-rays (Fabian et al. 1993; Madejski et al. 1995; Zdziarski et al. 1995a). It offers the best opportunity to address the above mentioned variability issues because it is among the brightest known Seyferts in X-rays (F(2-10keV)$``$1.6$`\times `$10-10erg cm-2s-1) and shows only moderate variability on short ($`<`$ day) timescales, but large variability on longer timescales (a factor 2-3 in a few days), both below and above 10 keV (Done, Madejski & Zycki 2000; Fabian et al. 1993).

BeppoSAX observed IC 4329A in 1996 for 100 ks (Perola et al. 1999), and then subsequently 3 times in 1998, with the 3 observations spaced approximately 5 days apart. As expected, the source did not vary much within each observation. Unfortunately, it did not vary much between the 4 observations either (see Fig. 1), despite the large and systematic flux variations typically observed in IC4329A by RXTE (Done, Madejski & Zycki 2000).

Spectral ratios were produced between all observations, in search of any model-independent spectral variations (Fig. 2). We find marginal evidence for variability of the spectral features at $`E<2`$ keV and at $`E=8`$–10 keV, but we find no clear variations in the FeKα and reflection component intensities, as well as in the high-$`E_{\mathrm{cutoff}}`$. The variations at low energies are most likely due to variability of the ionization state of the WA responding to continuum variations, but complications due to the presence of a possible scattered component and/or intrinsic soft (excess) component cannot be ruled out.

Spectral fitting of each observation with a complex model including 2 absorption edges, a steep power-law continuum, a reflection component with associated FeK$`\alpha `$, and a high-$`E_{\mathrm{cutoff}}`$ gave best-fit parameters similar to those reported by Perola et al. (1999) with somewhat larger errors (e.g. $`\mathrm{\Gamma }1.85\pm `$ 0.1, R $`0.6\pm 0.2`$, $`E_{\mathrm{cutoff}}`$ $``$ 250 $`\pm `$ 100 keV). Fig. 3 shows the marginal variations of the FeK line intensity and equivalent width. Fig. 4 shows the (similar) best-fit broad-band spectra of two different observations.

The most interesting result of the present analysis is that despite a 30% flux increase between obs. 1 and obs. 2, the line intensity became surprisingly weak (compare Fig. 1 and Fig. 3). During the higher flux obs. 3 the line re-established its typical value of $``$ 100-150 eV, though. In other words, there is some evidence (at a $``$ 2 $`\sigma `$ significance level) that the Fe K$`\alpha `$ line does not follow “instantaneously” (at least within 40 ks) the continuum variations, but it does on a timescale between 4–10 days.

## 3. Conclusions

* Moderate variations of the absorption structure, most likely related to the warm absorber, are apparently detected but require further investigation for more quantitative conclusions.
* Despite the unprecedented statistics of these observations, we find no significant spectral variations of the Compton reflection continuum, the power law spectral index and the high-enegy cutoff. This is probably because BeppoSAX unfortunately caught IC 4329A at a not too different flux level (30% variation) in all four observations hampering, thereby, any spectral variability.
* There is marginal evidence that the FeK$`\alpha `$ line is being produced by the outer regions of the accretion disk in IC 4329A, thereby suggesting a possible explanation for the lack of a clear broad FeK$`\alpha `$ line in this object.

## ACKNOWLEDGEMENTS

M.C. wishes to thank Martin Elvis and Fabrizio Nicastro for usefull discussions. This analysis has made use of data prepared by the BeppoSAX Scientific Data Center.

## REFERENCES

Cappi M., Mihara T., Matsuoka M., et al., 1996, ApJ, 458, 149

Done, C., Madejski, G. M. and Życki, P. T. 2000, ApJ, 536, 213

Fabian, A.C. et al., 1993, ApJ, 416, L57

Fiore, F., Perola, G.C., Matsuoka, et al., 1992, A&A, 262, 37

Miyoshi, S., et al., 1988, PASJ, 40, 127

Madejski, G.M., et al. 1995, ApJ, 438, 672

Matt, G. 2000, these proceedings, astro-ph/0007105

Perola, G.C., Matt, G., Cappi, M., et al., A&A, in press

Piro, L., Yamauchi, M., & Matsuoka, M., 1990, ApJ, 360, L35

Piro, L. 1999, in “Heating and Acceleration in the Universe”, Astron. Nachr., in press

Weaver, K. 2000, these proceedings, astro-ph/0007327


text:
# First principles study of strain/electronic interplay in ZnO; Stress and temperature dependence of the piezoelectric constants.

## I Introduction

Zinc oxide (ZnO) is a tetrahedrally coordinated wide band gap semiconductor that crystallizes in the wurtzite structure (Figure 1). The lack of center of symmetry, combined with a large electromechanical coupling, result in strong piezoelectric properties, and the consequent use of ZnO in mechanical actuators and piezoelectric sensors. In addition, ZnO is transparent to visible light and can be made highly conductive by doping. This leads to applications in surface acoustic wave devices and transparent conducting electrodes. However the piezoelectric properties can change the characteristics of potential energy barriers to mobile charges at interfaces, and hence affect the carrier transport properties. The resulting piezoresistance is at times desirable, for example in ZnO-based metal-oxide varistors which can dissipate large amounts of power in short response times and are commonly found as electrical surge protectors. However the detailed effects of piezoelectrically induced changes on the electrical behavior of ZnO have not yet been well characterized, and as ZnO finds increased application in electronic devices these effects will have large technological impact.

In this paper we present a first-principles study of the strain dependence of the electrical properties of ZnO, using a plane wave pseudopotential (PWPP) implementation of density functional theory (DFT) within the local density approximation (LDA). We calculate and analyze the structural dependence of total energies, band structures and piezoelectric and dielectric constants. The principal result of our analysis is that the piezoelectric constants of ZnO are strongly dependent on conditions of stress and temperature, whereas the dielectric coefficients vary less strongly. Of particular interest is a comparison of the properties of ZnO, in which the Zn $`3d`$ bands are filled with those of PbTiO3, which has empty Ti $`3d`$ states. Both materials have anomalously large piezoelectric coefficients, but in PbTiO3 the dielectric coefficients are also anomalously large, whereas ZnO behaves as an ordinary dielectric. Our hypothesis is that the strong O $`2p`$ \- Zn $`3d`$ hybridization which has been previously noted in PbTiO3 also occurs in ZnO, leading in both materials to strong strain-phonon coupling and consequently large piezoelectric coefficients. Thus a filled $`3d`$ band does not preclude a large peizoelectric response. The large dielectric response in PbTiO3 has a similar origin. In ZnO however, the filled O $`2p`$ and Zn $`3d`$ bands do not allow a large orbital response to an applied electric field, resulting in a normal dielectric response.

The remainder of this paper is organized as follows. In Section II we summarize the results of earlier theoretical studies of ZnO. In Section III we describe the theoretical and computational methods used in this work. In Section IV, we describe our results for the static and response properties of bulk ZnO at its equilibrium lattice constant. In Section V we investigate the dependence of the electronic structure and response functions on external stress or changes in temperature. Finally, in Section VI we present our conclusions and discuss implications for materials growth and device design.

## II Previous theoretical work

First-principles studies of ZnO are computationally challenging. First, the wurtzite structure contains twice as many atoms per unit cell as the zincblende semiconductor structure. In addition, both oxygen and zinc are problematic atoms for the construction of pseudopotentials. In both cases the relevant valence electrons (O $`2p`$ and Zn $`3d`$) have no lower-lying electrons of the same angular momentum to provide an effective repulsive potential from the orthogonalization requirement. As a result they are tightly bound and require a large number of plane waves in their expansion.

The first published band structure of ZnO used the Green’s function KKR method. This was followed by empirical pseudopotential calculations in which the Zn $`3d`$ electrons were placed in the core. In addition to preventing assessment of the Zn $`3d`$ contribution to bonding properties, this approach was later shown to give unsatisfactory results.

Schroer and co-workers circumvented the problem of a large plane wave basis set by combining the use of pseudopotentials with localized Gaussian basis sets containing orbitals of $`s`$, $`p`$, $`d`$ and $`s^{}`$ symmetry. Using this basis they compared the results of LDA calculations using Zn2+ (in which the $`3d`$ electrons are contained in the core) and Zn12+ (in which the $`3d`$ electrons are valence electrons) pseudopotentials, and found that the Zn12+ pseudopotential gave results in good agreement with experiment. For example, their LDA energy minimum volume was 0.6 % below the experimental value and their calculated band structure was in reasonable agreement with angle resolved photoemission measurements. There was a slight discrepancy in the position of the $`d`$ bands which they attributed to inadequacy of the local density approximation in describing these strongly correlated bands.

Dal Corso et al. avoided the use of pseudopotentials entirely by using the all electron full potential linear augmented plane wave (FLAPW) method. They calculated the piezoelectric and polarization properties of ZnO within the LDA, with the purpose of determining the origin of the unusually strong piezoelectric response in ZnO. The principal result of their work was that the contribution to the macroscopic polarization tensor from the relative displacement of the sublattices was large, and only partly canceled by the electronic “clamped-ion” contribution, leading to a large net piezoelectric polarization. (In contrast, in zincblende semiconductors these terms are of similar magnitude and opposite sign, resulting in a small piezoelectric polarization.) As in Ref. , their LDA volume slightly underestimated the experimental value, and their Zn $`d`$ bands were around 4 eV higher in energy than observed in photoemission.

Hartree-Fock calculations have been used successfully to determine the stability of, and transitions between, different phases of ZnO and to calculate the ($`101\overline{0}`$) surface reconstruction for the wurtzite phase. The Hartree-Fock approximation also gives an incorrect energy position for the Zn $`3d`$ bands - this time around 2 eV too low.

## III Computational techniques

### A Pseudopotential construction

The calculations described in this work were performed using a plane wave pseudopotential implementation of density functional theory within the local density approximation. Plane wave basis sets offer many advantages in total energy calculations for solids, including completeness, an unbiased representation, and arbitrarily good convergence accuracy. They also allow for straightforward mathematical formulation and implementation, which is invaluable in the calculation of Hellmann-Feynman forces and in the density functional theory linear response calculations employed here.

However plane wave basis sets necessitate the use of pseudopotentials to model the electron-ion interaction, in order to avoid rapid oscillations of the valence wavefunctions in the region around the ion cores. The difficulties associated with applying the pseudopotential method to tightly bound $`d`$-electrons, which might be expected to require a prohibitively large number of plane waves to expand their pseudopotentials, were mentioned above. An earlier study of cubic ZnS using the smooth Trouiller-Martins pseudopotentials, required plane waves up to 121 Ry in energy to achieve an energy convergence of 0.05 eV. Although feasible for a bulk calculation for the zincblende structure (with only two atoms per unit cell) such a large energy cutoff is undesirable for larger unit cells, such as that of the wurtzite structure, or those required for calculation of surface properties. In this work, we use the optimized pseudopotentials developed by Rappe et al., which allow us to reduce the required energy cutoff to 64 Ry without compromising accuracy or transferability. Optimized pseudopotentials minimize the kinetic energy in the high Fourier components of the pseudo wavefunction, leading to a corresponding reduction in the contribution of high Fourier components in the solid.

For both Zn and O we constructed non-relativistic optimized pseudopotentials. The oxygen pseudopotentials were generated from a $`2s^22p^4`$ reference configuration with core radii, $`r_c`$, of 1.5 a.u. for both $`s`$ and $`p`$ orbitals. They were then optimized using 4 and 3 basis functions with cutoff wave vectors, $`q_c`$, of 7.0 and 6.5 a.u. for $`s`$ and $`p`$ orbitals respectively. $`q_c`$ determines the convergence of the kinetic energy with respect to the plane wave cutoff energy in reciprocal space calculations. These oxygen pseudopotentials were used in earlier calculations for perovskite oxides and gave accurate results. The zinc pseudopotentials were constructed for a neutral Zn atom with reference configuration $`3d^{10}4s^{1.75}4p^{0.25}`$. $`r_c`$ values of 2.0, 1.4 and 1.4 a.u.s were used for $`d`$, $`s`$ and $`p`$ orbitals respectively, with $`q_c`$ values of 8.0, 7.0 and 8.0 Ry (giving a cutoff energy of 64 Ry.) The transferability of the pseudopotential was tested for a variety of +1 and +2 free Zn ions. The pseudo total energies and eigenvalues were in agreement with the all electron values to within 0.001 a.u.s. There was no improvement in transferability, or in agreement with all electron calculations for bulk systems, on inclusion of non linear core corrections. All pseudopotentials were put into separable form using one projector for each angular momentum. For both Zn and O the $`l=1`$ component was chosen as the local potential. The absence of ghost states was confirmed using the ghost theorem of Gonze, Käckell and Scheffler.

### B Density functional theory linear response

We use density functional theory linear response (DFT-LR) to obtain the quadratic couplings between homogeneous strain, internal displacements of atoms and macroscopic electric field. We use a variational formulation of DFT-LR in which the second derivative of total energy is minimized with respect to the first derivatives of Kohn-Sham wave functions with appropriate orthogonality constraints. This method avoids using any finite-difference formulae and yields dielectric or piezoelectric constants with a minimal number of calculations.

To obtain the piezoelectric and dielectric constants, we calculate the DFT-LR of our system to two types of perturbations: (a) phonon (or atomic displacements) and (b) electric field. Using the first-order response wavefunction resulting from (a) in the Hellman-Feynman force formula and in the stress formula, we obtain the dynamical matrix and the coupling between phonons and strain respectively. Similarly, the response wavefunctions resulting from (b) are used to obtain the Born effective charges and the clamped-ion piezoelectric constants respectively.

The symmetry of the wurtzite structure allows three independent piezoelectric constants ($`\gamma _{33},\gamma _{13},\gamma _{14}`$) and two dielectric constants ($`ϵ_{33},ϵ_{13}`$). In the present work, we focus on $`\gamma _{33},\gamma _{13}`$ and $`ϵ_{33}`$. The piezoelectric constants $`\gamma _{33}`$ and $`\gamma _{13}`$ give the polarization along the $`c`$-axis induced by strains $`e_{33}`$ and $`e_{11}`$ respectively. Equivalently, $`\gamma _{33},\gamma _{13}`$ give the stresses $`\sigma _{33},\sigma _{11}`$ induced on the unit cell by an electric field along the $`c`$-axis. The former relationship underlies the earlier work of Dal Corso et al based on finite-difference formulae and geometric phase, and the latter is used in the present work based on DFT linear response.

To obtain the piezoelectric and dielectric constants, we perform two DFT-LR calculations, one with phonon perturbations corresponding to Zn and O displacements along the $`c`$-axis, and another with a perturbing electric field along the $`c`$-axis. The phonon perturbation allows us to obtain both the frequencies of the $`\mathrm{\Gamma }`$-point phonons with $`z`$-polarization (which are proportional to the square roots of the respective force constants, $`\sqrt{K_\alpha }`$), and the strain-phonon couplings, $`L_{\alpha ,zz}`$ and $`L_{\alpha ,xx}`$. The LR calculations with the electric field perturbation yield the optical dielectric constant $`ϵ_{33}^{\mathrm{}}`$, the Born effective charges $`Z_{\alpha ,zz}`$, and the clamped-ion piezoelectric couplings $`\gamma _{33}^0`$ and $`\gamma _{13}^0`$. The dielectric and piezoelectric constants are then given by:

$$ϵ_{33}=ϵ_{33}^{\mathrm{}}+4\pi \underset{\alpha }{}\frac{Z_{\alpha ,zz}Z_{\alpha ,zz}}{K_\alpha }$$

$$\gamma _{33}=\gamma _{33}^0+\underset{\alpha }{}\frac{Z_{\alpha ,zz}L_{\alpha ,zz}}{K_\alpha }$$

$$\gamma _{13}=\gamma _{13}^0+\underset{\alpha }{}\frac{Z_{\alpha ,zz}L_{\alpha ,xx}}{K_\alpha }$$

.

### C Technical Details

All calculations were performed on Silicon Graphics O2 and Origin 200 systems using the conjugate gradient program CASTEP 2.1 and our own related density functional linear response program. We used a plane wave cut off of 64 Ry, which corresponds to around 3000 plane waves per wavefunction in a a single ZnO wurtzite unit cell. A 3x3x2 Monkhorst-Pack grid was used, leading to six k-points in the irreducible Brillouin Zone for the high symmetry structures, and a correspondingly higher number for distorted structures with lower symmetry. The exchange correlation was parameterized using the Perdew-Zunger parameterization of the Ceperley-Alder potential .

## IV Unstrained ZnO; New results and comparison with LAPW calculations

Before presenting our results for strained ZnO, we first discuss bulk, unstrained ZnO, and compare the results of our plane wave pseudopotential calculations with published theoretical and experimental data. A good test of the transferability of the pseudopotentials is that they predict the same minimum energy structure as all-electron calculations which use the same approximation for the exchange-correlation functional. The minimum energy volume obtained by the earlier all-electron LDA calculation (Ref. ) is 45.89 Å3. Ref.  also showed that the LDA $`\frac{c}{a}`$ ratio and $`u`$ value are the same as the corresponding experimental quantities to within 0.5 mRy, and that the energy surface is rather flat around the minimum in the ($`u,\frac{c}{a}`$) plane. The minimum energy volume obtained using our pseudopotentials is 47.6 Å3 (at the experimental $`\frac{c}{a}`$ ratio and $`u`$ value), and that obtained by an earlier pseudopotential calculation using a Gaussian basis is 46.80 Å3. The experimental volume is 47.90 Å3 and we believe the results obtained using our method are within LDA errors.

### A Band structure analysis

In Figure 2 we show the calculated band structure for ZnO at the experimental unit cell volume, ($`a=3.2595`$ Å, $`c=5.2070`$ Å and $`u=0.3820`$), which is approximately equal to the LDA minimum volume for our pseudopotentials. The top of the valence band is set to 0 eV. The band structure is indistinguishable from that presented in Ref.  which was calculated using the FLAPW technique.

The narrow bands between around -6 eV and -5 eV derive largely from the Zn $`3d`$ orbitals and are completely filled. The broad bands between around -5 eV and 0 eV are from the O $`2p`$ orbitals and are again completely filled. The Zn $`4s`$ band is broad and unoccupied, ranging between 1 and 7 eV above the band gap. Figure 3 shows the same band structure along $`\mathrm{\Gamma }`$ to A with the symmetry labels within the $`C_{6v}`$ group added. 1 is the totally symmetric representation, 2 is antisymmetric with respect to $`C_6`$ and $`C_2`$ rotation and $`\sigma _v`$ reflection, 5 and 6 are the doubly degenerate representations, with 5 being antisymmetric with respect to $`C_6`$ and $`C_3`$ rotation, and 6 being antisymmetric with respect to $`C_3`$ and $`C_2`$ rotation. We see that interactions between O $`2p`$ and Zn $`3d`$ orbitals are allowed by symmetry at $`\mathrm{\Gamma }`$ and A, and along the adjoining $`\mathrm{\Delta }`$ line. The Zn $`s`$ orbitals can interact with O $`2p`$ and Zn $`3d`$ bands of 1 and 3 symmetry.

In order to quantify the interactions between the various orbitals we perform a tight-binding analysis along the $`\mathrm{\Gamma }`$ to A direction of the Brillouin zone. Tight-binding parameters are obtained by non-linear-least-squares fitting to the calculated ab initio energies at the high symmetry $`\mathrm{\Gamma }`$ and A points, and at 19 points along the $`\mathrm{\Delta }`$ axis. A good tight-binding fit (rms deviation = 0.11) is obtained when only nearest neighbor interactions between O $`2p`$ and Zn $`3d`$ bands, and O $`2p`$ and Zn $`4s`$ bands are included in the fit. Additional small Zn $`3d`$ \- Zn $`3d`$ interactions are needed to produce dispersion in the upper $`e_g`$ Zn $`3d`$ band along this symmetry axis. The tight-binding parameters which we obtain are given in Table II (in the column labeled ‘structure 1’), and the tight-binding band structure is compared with the ab initio band structure in Figure 4.

The two sets of Zn-O parameters correspond to the two Zn-O distances, the shorter ( $`r_1=1.974`$ Å ) being the separation of the Zn and O atoms lying directly above each other along the $`c`$ direction, and the longer ( $`r_2=1.989`$ Å ) joining Zn and O atoms in adjacent $`c`$-oriented ‘chains’. The Zn $`3d`$ \- Zn $`3d`$ interactions are small. The Zn $`4s`$ \- O $`2p`$ interactions are large, and show the expected increase when the Zn - O spacing is decreased. The Zn $`3d`$ \- O $`2p`$ interactions are also large, and are quite different (even changing sign) for the two different Zn - O atom pair types. In fact both the $`\sigma `$ and $`\pi `$ components are larger for the in-plane pairs which have the larger Zn-O spacing. This distance dependence is unusual for tight-binding parameters, and suggests that the nature of the Zn $`3d`$ \- O $`2p`$ hybridization is different for the Zn - O pairs lying in the $`c`$ axis chains, than for the basal plane Zn - O pairs.

### B Piezoelectric and dielectric properties

In Table I, we present our results for piezoelectric and dielectric compliances and a comparison with those obtained in previous FLAPW calculations and experiments. The agreement between the computational results is good, but the scatter in experimental results indicates that the piezoelectric constants may be sensitive to experimental conditions. We will investigate this issue in the next section. We also report in Table I the frequencies of TO phonons with $`c`$ polarization at the $`\mathrm{\Gamma }`$ point, and note that only one of the TO phonons (395 cm-1) is IR-active; this is the only phonon which contributes to piezoelectric and dielectric constants $`\gamma _{33},\gamma _{13}`$ and $`ϵ_{33}`$, studied in this work. Our results also confirm the conclusion of Ref.  that the piezoelectric constants of ZnO are dominantly contributed by phonons (internal strain). The dielectric constant, on the other hand, has roughly equal contributions from electrons ($`ϵ^{\mathrm{}}`$) and phonons.

## V Effects of strain and temperature on the electronic properties of bulk ZnO

### A Band structures

Next we investigate the effects of strain on the static electronic properties of bulk ZnO. We compare two different strained structures with the equilibrium structure at the experimental lattice constant. First we simulate application of a homogeneous in-plane compression by reducing the $`a`$ lattice constant by 2%, while increasing the $`c`$ lattice constant correspondingly to maintain the same total volume. The value of $`u`$ is held at 0.3820 of the $`c`$ axis. Second we investigate the effect of changing the Zn-O separation along the $`c`$ axis, $`u`$, by increasing it by 5% to 0.4011, while retaining the equilibrium $`a`$ and $`c`$ values.

#### 1 In-plane compression

Application of the in-plane compression to reduce the $`a`$ lattice constant by 2% increases the total energy by 0.1 eV, and creates Hellman-Feynman forces on the atoms in the $`z`$ direction of $`\pm 0.59\frac{\text{eV}}{\text{Å}}`$. The band structure of the compressed structure is shown in Figure 5. We observe a broadening of the bands, as expected from the increased overlap between the orbitals in the basal plane. Note in particular the broadening in the Zn $`3d`$ bands. In spite of the fact that the $`d`$ bands are narrow, they are very sensitive to strain and bond length as a result of $`pd`$ hybridization processes.

In addition to an overall band broadening, there are a number of significant details in the band structure for this compressed structure. First, the O $`2p`$ bands with symmetries 1 and 3 shift down relative to the uppermost O $`2p`$ bands at the $`\mathrm{\Gamma }`$ point. This leads to a reduction in the density of states at the Fermi level, and to an overlapping of the low energy part of these bands with the Zn $`3d`$ bands. The ordering of the lower Zn $`3d`$ bands is reversed compared with those in the unstrained structure, and the Zn $`4s`$ bands shift slightly down in energy, resulting in a smaller band gap.

These observations are consistent with a tight-binding fit along the $`\mathrm{\Gamma }`$ to $`A`$ line, the results of which are given in Table II (structure 2). Again, a good (root mean square deviation = 0.11) tight-binding fit is obtained using the limited interaction set described above. Again the O $`2p`$ \- Zn $`4s`$ hybridization is large and shows the expected variation with bond length. The dependence of the (also large) $`pd`$ parameters on distance does not follow a straightforward pattern. Note that the anomalously low value for the Zn $`4s`$ energy is the result of an attempt by the fitting package to reproduce the down shift in the O $`2p`$ band of 1 symmetry (which also contains a significant Zn $`4s`$ component) within our limited basis set. This energy value shifts up to a more physical positive number if additional interactions are included in the basis.

#### 2 Change in $`u`$ value

Change in the $`u`$ value by 5 % increases the total energy by only 0.07 eV compared with the equilibrium structure, and introduces Hellman Feynman forces of $`\pm 0.59`$ eV in the $`z`$ direction on the ions. The new band structure is shown in Figure 6, and is qualitatively very similar to that of the undistorted structure.

A tight-binding fit using the interaction set described above produces an rms deviation of 0.13, and the parameters listed in Table II. Again the O $`2p`$ \- Zn $`4s`$ interaction is largest for the smallest bond length as expected, but the $`pd`$ hybridizations have a more complicated distance dependence. In this case, the O $`2p`$ bands with 1 and 3 symmetries are shifted down a small amount at $`\mathrm{\Gamma }`$, but not as much as in the previous structure. As a result they do not overlap with the Zn $`3d`$ bands. The Zn $`4s`$ bands are not shifted down relative to their position in the unstrained case.

### B Piezoelectric properties

#### 1 Stress dependence

Finally we explore the piezoelectric and dielectric properties of strained ZnO in detail. As a result of anharmonic couplings between different phonons and strain, the properties of ZnO are strongly structure sensitive. There is a large region in the phase space of structural configurations which is low in energy and therefore contributes to the properties of ZnO. To obtain energies, polarization and other properties for each of these configurations from first-principles is inefficient and impractical. We choose instead to use a model energy functional that captures the physics of the low-energy structural excitations of ZnO.

To keep the model simple, we restrict our analysis to the subspace of degrees of freedom that preserve the symmetry of the wurtzite structure { $`x=e_{xx}=e_{yy}`$, $`z=e_{zz}`$, $`u`$ }. Justification for this simplification rests on the assumption (verifiable through ab initio calculations) that most of the low energy configurations preserve wurtzite symmetry and most of the symmetry breaking distortions are described well at harmonic order and can be integrated out. This is in contrast with earlier work on structural phase transitions where the symmetry-breaking distortions were retained in the subspace.

We write the model energy functional as a symmetry-invariant Taylor expansion in atomic displacements (or normal mode degrees of freedom $`u`$) and strains ($`x`$, $`z`$) defined above. Including the lowest order coupling with external stress $`\sigma _\alpha `$ and electric field along c-axis $`E`$,

$$E_{tot}(x,z,u,\sigma _\alpha ,E)=\frac{1}{2}Ku^2+Au^3+Bu^4+\frac{1}{2}[C_xx^2+C_zz^2+C_{xz}xz]$$

$$+D_xx^3+E_xx^4+D_zz^3+E_zz^4+\underset{\alpha }{}F_\alpha r_\alpha u^2$$

$$+\underset{\alpha ,\beta }{}G_{\alpha \beta }r_\alpha r_\beta u+\underset{\alpha ,\beta }{}H_{\alpha \beta }r_\alpha r_\beta u^2$$

$$\mathrm{\Omega }(2\sigma _xx+\sigma _zz+2\gamma _{13}^0xE+\gamma _{33}^0zE)ZuE\frac{\mathrm{\Omega }}{4\pi }ϵ_{\mathrm{}}E^2,$$
(1)
where the $`\gamma ^0`$’s and $`ϵ_{\mathrm{}}`$ are clamped (or electronic) piezoelectric and dielectric constants, and $`r_\alpha `$ is strain $`x`$ or $`z`$. $`K`$, $`A`$$`H`$, $`\gamma `$ and $`ϵ`$ are the harmonic and anharmonic coupling parameters including force, elastic and mixed coupling constants. These parameters have been determined from DFT total energy and linear response calculations.

The equilibrium state of ZnO under the application of external stress or electric field at zero kelvin is obtained by minimizing the total energy $`E_{tot}`$ with respect to the structural parameters $`x`$, $`z`$ and $`u`$:

$$\frac{E_{tot}}{x}=0,\frac{E_{tot}}{z}=0,\frac{E_{tot}}{u}=0.$$
(2)
For the equilibrium structure, the static dielectric and piezoelectric constants are calculated using the expressions in Section III (also in Ref. ), with various coupling parameters being dependent on the structure. The spontaneous polarization is calculated using the expression

$$P=\frac{E_{tot}}{E}.$$

ZnO is grown in the form of thin solid films on sapphire and there are always interfacial stresses in the films. We use our model to study properties of ZnO as a function of stress in the basal $`ab`$ plane. In Fig. 8, we show our results for the structural parameter $`u`$, dielectric and piezoelectric response of ZnO as the applied stress $`\sigma =\sigma _{xx}=\sigma _{yy}`$ is varied from -1 to 1 GPa.

We find that parameter $`u`$ changes by about 1 percent in this range of basal stress. The dielectric constant $`ϵ_{33}`$ monotonically increases by about 2 to 3 percent. In contrast, the piezoelectric constants are rather sensitive to stress changing by about 15 to 30 %. With access to different contributions to these constants in our model, we discover that most of the dependence of these compliances on stress is due to the phonon contribution.

#### 2 Temperature Dependence

ZnO is also a pyroelectric material, characterized by a temperature dependence of the polarization. This has technological relevance because it leads to the widespread use of ZnO in infra-red detectors. To explore various stress and electric field-dependent properties of ZnO at finite temperature, we obtain free energy functional using a local harmonic model for entropy. In this model, entropy is calculated treating phonons harmonically for given structural parameters. In the present work, we include the optical phonons that are polarized along $`z`$-axis. With these approximations, the free energy is:

$$G(x,z,u,\sigma _x,\sigma _z,E,T)=E_{tot}(x,z,u,\sigma _x,\sigma _z,E)+3k_BTLog(\frac{(K_1K_2K_3)^{\frac{1}{3}}}{k_BT})$$
(3)
where the $`K_i`$’s are the harmonic force constants of the three $`\mathrm{\Gamma }`$phonons with $`z`$-polarization for given values of structural parameters. Due to the anharmonic terms included in the energy expansion Eqn (1), $`K_i`$’s (hence the T-dependent part of the free energy) depend on the structure. The equilibrium state of ZnO under the application of external stress or electric field at finite temperature is then obtained by minimizing the free energy $`G`$ with respect to the structural parameters $`x`$, $`z`$ and $`u`$:

$$\frac{G}{x}=0,\frac{G}{z}=0,\frac{G}{u}=0.$$
(4)
Again, the dielectric and piezoelectric constants are calculated using the expressions in Section III, and the spontaneous polarization is now calculated using the expression

$$P=\frac{G}{E}.$$

We investigate the dependence of various properties of ZnO on temperature. In Fig. 7, we display results for structural parameters, dielectric and piezoelectric constants, and polarization for a range of temperatures from 0 to 450 K. We find that the structural parameters such as the bond-length $`u`$ and volume change only by about 0.3 % from zero kelvin to room temperature.

The pyroelectric constant, which we obtain from our calculated results for polarization, is 20 $`\mu C/m^2/K`$, compared with the experimental value of 9.4 $`\mu C/m^2/K`$. Considering the simplicity in our treatment of temperature through the model entropy function, we find this agreement quite encouraging. In particular, we point out that, while the nonpolar TO phonons with $`z`$-polarization have been omitted from the expansion of energy, they have been included in the entropy term, where we use the expression of entropy treating phonons harmonically, consistent with our energy expansion.

To estimate the effect of these non-polar phonons on the pyroelectric constant, we omitted their contribution to entropy and found a pyroelectric constant of about 40 $`\mu C/m^2/K`$ (almost doubled). If we assume that phonons should generally suppress the pyroelectric constant, the discrepancy between theory and experiment is likely due to our omission of phonons with $`x`$ and $`y`$-polarization from the entropy expression.

The linear dielectric response $`ϵ_{33}`$ of ZnO is quite sensitive to temperature, changing by about 4% in the temperature range considered. The piezoelectric response, on the other hand, is very sensitive to temperature changing by about 20%. This should be an important consideration in designing piezoelectric devices for operation at room temperature.

### C Discussion

It is clear from our results that the piezoelectric response of ZnO is strongly sensitive to both temperature and stress, changing by up to 30 % over the range of parameters considered. This dependence arises from the changes in structural parameters (manifested through the phonon contribution). We saw in Section III that the phonon contribution to the piezoelectric constants arises from the coupling of phonons with strain, $`L`$, and Born effective charge, $`Z`$,

$$\gamma _{phonon}=\frac{LZ}{K},$$

$`K`$ being the force constant. In Fig. 9, we show how $`K`$ and $`L`$ for the polar TO phonon change with temperature. While $`K`$ changes by only 10 %, the coupling with strain, $`L`$, changes by about 25 %. The Born effective charge $`Z`$ (which describes the coupling of this phonon with electric field) is not found to vary much with temperature. The large temperature dependence of the piezoelectric response arises predominantly from that of the coupling of phonon with strain and its force constant. Since only the latter contributes to the dielectric constants, dielectric properties are less sensitive to structure, stress or temperature.

In the Section V A, we found the hybridization between Zn $`d`$ orbitals and O $`p`$ orbitals to be sensitive to structural parameters. The same hybridization was found to be the cause of anomalous Born effective charges in ferroelectric materials such as PbTiO3. While the $`d`$ orbitals of the transition metals in perovskite ferroelectrics are formally unoccupied, those in Zn are fully occupied, leading to normal effective charges. The coupling of phonons with strain, however, is large and structure dependent irrespective of the occupancy of $`d`$ orbitals.

## VI Summary

In summary, we have calculated the electronic and atomic structure of ZnO from first-principles, and analyzed the nature of the bonding using the tight-binding method. We find that hybridization between the Zn $`d`$ orbitals and O $`p`$ orbitals is strongly structure dependent. Using DFT linear response, we have obtained the phonon frequencies, dielectric and piezoelectric constants of ZnO at zero temperature, and have shown that phonons (internal strain) have the dominant contribution to piezoelectricity in ZnO. ¿From DFT linear response and total energy calculations, and a simple model for vibrational entropy, we have constructed an ab initio free energy functional for ZnO to study its properties at finite temperature and under applied stress. Our results show that the piezoelectric properties of ZnO are strongly dependent on both temperature and stress. This clearly has implications for the design of devices intended to operate at room temperature, or under stressed conditions. By analyzing various physical contributions, we have found that this is primarily due to the coupling between phonons and strain. The O $`2p`$ \- Zn $`3d`$ hybridization is the cause of the large magnitude and sensitivity of this coupling.

## VII Acknowledgments

The authors thank Andrew Rappe and Nicholas Ramer for their assistance in generating the Zn pseudopotential, and Paul Verghese and V. Srikant for many useful discussions. Many of the ideas presented in this paper were formulated during the 1998 “Physics of Insulators” workshop at the Aspen Center for Physics. N.A. Hill’s funding was provided through the NSF-POWRE program, grant # DMR-9973859. U. V. Waghmare’s funding was provided by the Office of Naval Research, grant \# N00014-97-1-1068.