text:
# Local structure study of 𝐈𝐧_𝐱𝐆𝐚_{𝟏-𝐱}𝐀𝐬 semiconductor alloys using High Energy Synchrotron X-ray Diffraction

## I Introduction

Semiconductor alloys are known as technologically important materials for their wide applications in optoelectronic devices such as lasers and detectors.  The local structure information is of fundamental importance in understanding the alloy systems because their physical properties are strongly influenced by the local atomic displacements present in the alloys. For example, it is known that chemical and compositional disorder strongly affect the electronic structure of zinc-blende type alloys  and their enthalpies of formation.

In this paper we present a detailed study of the local and average structure of the $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$ alloy series. The average structure of $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$ was studied by Woolley.  The structure is of the zinc blende type ($`F4\overline{3}m`$ over the entire alloy range. The lattice parameters, and therefore the average In-As and Ga-As bond lengths, interpolate linearly between the values of the end-members according to Vegard’s law.  However, consideration of the local structure reveals a very different situation. The local structure of $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$ was first studied by Mikkelson and Boyce using extended x-ray absorption fine structure (XAFS).  According to this experiment, the individual nearest neighbor (NN) Ga-As and In-As distances in the alloys are rather closer to the pure Ga-As and In-As distances. Further XAFS experiments showed that this is quite general behavior for many zinc-blende type alloy systems.  Since then a number of theoretical and model studies have been carried out on the semiconductor alloys to understand how the alloys accommodate the local displacements.

Until now these models and theoretical predictions are tested mainly by the comparison with XAFS data. The XAFS results give information about the nearest neighbor and next nearest neighbor distances in the alloys but imprecise information about bond-length distributions and no information about higher-neighbor shells. This limited structural data makes it difficult to differentiate between competing models for the local structure. For example, even a simple radial force model  rather accurately predicts the nearest neighbor distances of $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$ alloys in the dilute limit. Therefore, one needs more complete structural information including nearest neighbor, far-neighbor distances, and bond length distributions to prove the adequacy of model structures for these alloys.

The atomic pair distribution function (PDF), $`G(r)`$, measures the probability of finding an atom at a distance $`r`$ from another atom.  One of the advantages of the PDF method over other local probes such as XAFS is that it gives both local and intermediate range information because both Bragg peaks and diffuse scattering are used in the analysis. It is also possible to obtain information about the static bond length distribution from the PDF peak width and about correlations of atom displacements.

In this paper, we present a detailed X-ray diffraction study of $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$, $`(0x1)`$. A preliminary analysis of the data has been published elsewhere.  Using high energy synchrotron x-rays, we measured the total scattering structure function, $`S(Q)`$, of the $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$ alloy system extended to high $`Q`$ ($`Q_{max}=45`$ Å-1) where $`Q`$ is the magnitude of the momentum transfer of the scattered x-rays ($`Q=4\pi \mathrm{sin}\theta /\lambda `$ for elastic scattering). From these structure functions we obtained the corresponding high real-space resolution PDFs through a Fourier transform according to

$$G(r)=\frac{2}{\pi }_0^{\mathrm{}}Q[S(Q)1]\mathrm{sin}QrdQ.$$
(1)
In these PDFs, the first peak is clearly resolved into two sub-peaks corresponding to the Ga-As and In-As bond lengths.  The evolution of the bond-length with doping gives good agreement with XAFS. For the far-neighbor peaks, the peak widths are much broader in the alloy samples than those of the pure end-members reflecting the increased disorder. We model the local structure of $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$ alloys using a supercell model  based on the Kirkwood potential  which gives good agreement with the alloy data with no adjustable parameters. The results of the modeling have been analyzed to reveal the average atomic static distribution on the As and (In,Ga) sublattices. Finally, we have calculated the diffuse scattering that one would get from the Kirkwood model. This compares qualitatively well with published diffuse scattering results from $`\mathrm{In}_{0.53}\mathrm{Ga}_{0.47}\mathrm{As}`$

## II Experimental details

### A Data collection

The alloy samples, with compositions $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$, (x = 0, 0.17, 0.33, 0.5, 0.83, 1) were prepared by a melt and quench method. An appropriate fraction of InAs and GaAs crystals were powdered, mixed and sealed under vacuum in quartz ampoules. The samples were heated beyond the liquidus curve of the respective alloy  to melt them and held in the molten state for 3 hours before quenching them in cold water. The alloys were powdered, resealed in vacuum, and annealed just below the solidus temperature for 72-96 hours to increase the homogeneity of the samples. This was repeated until the homogeneity of the samples, as tested by x-ray diffraction, was satisfactory. X-ray diffraction patterns from all the samples showed single, sharp diffraction peaks at the positions expected for the nominal alloy similar to the results obtained by Mikkelson and Boyce.

High energy x-ray powder diffraction measurements were conducted at the A2 wiggler beamline at Cornell High Energy Synchrotron Source (CHESS) using intense x-rays of 60 KeV ($`\lambda =0.206`$ Å). The incident x-ray energy was selected using a Si(111) double-bounce monochromator. All measurements were carried out in flat plate symmetric transmission geometry. In order to minimize thermal atomic motion in the samples, and hence increase the sensitivity to static displacements of atoms, the samples were cooled down to 10 K using a closed cycle helium refrigerator mounted on the Huber 6 circle diffractormeter. The samples were uniform flat plates of loosely packed fine powder suspended between thin foils of kapton tape. The sample thicknesses were adjusted to achieve sample absorption $`\mu t1`$ for the 60 KeV x-rays, where $`\mu `$ is the linear absorption coefficient of the sample and $`t`$ is the sample thickness.

The experimental data were collected up to $`Q_{max}=45`$ Å-1 with constant $`\mathrm{\Delta }Q`$ steps of 0.02 Å-1. This is a very high momentum transfer for x-ray diffraction measurements. For comparison, $`Q_{max}`$ from Cu $`K_\alpha `$ x-ray tube is less than 8 Å-1. This high $`Q_{max}`$ is crucial to resolve the small difference ($`0.14`$ Å) in the In-As and Ga-As bond lengths.

To minimize the measuring time, the data were collected in two parts, one in the low $`Q`$ region from 1 to 13 Å-1 and the other in mid-high $`Q`$ region from 12 to 50 Å-1. Because of the intense scattering from the Bragg peaks, in the low $`Q`$ region the incident beam had to be attenuated using lead tape to avoid detector saturation. The maximum intensity was scaled so that the count rate across the whole detector energy range in the Ge detector did not exceed $`2\times 10^4s^1`$. At these count-rates detector dead-time effects are significant but can be reliably corrected as we describe below. To reduce the random noise level below 1%, we repeated runs until the total *elastic scattering* counts become larger than 10,000 counts at each value of $`Q`$. Also, to obtain a better powder average the sample was rocked with an amplitude of $`\pm 0.5^{}`$ at each $`Q`$-position. The scattered x-rays were detected using an intrinsic Ge solid state detector. The signal from the Ge detector was processed in two ways. The signal was fed to a multi-channel analyzer (MCA) so that a complete energy spectrum was recorded at each data-point. The signals from the elastic and Compton scattered radiation could then be separated using software after the measurement. In parallel, the data were also fed through single-channel pulse-height analyzers (SCA) which were preset to collect the elastic scattering, Compton scattering, and a wider energy window to collect both the elastic and Compton signals. For normalization, the incident x-ray intensity was monitored using an ion chamber detector containing flowing Ar gas.

For the SCAs, the proper energy channel setting for the elastic scattering is crucial. Any error in the channel setting could cause an unknown contamination by Compton scattering and make data corrections very difficult. There’s no such problem in the MCA method since the entire energy spectrum of the scattered radiation is measured at each value of $`Q`$. The main disadvantage of the MCA method is that it has a larger dead-time, although this can be reliably corrected as we show below. Fig. 1 shows a representative MCA spectrum

taken from the InAs sample at $`Q=45`$ Å. It is clear that the Compton and elastic scattering are well resolved at this high momentum transfer. The elastically scattered signal, which contains the structural information, is obtained by integrating the area under the elastic scattering peak.

### B Data Analysis

The measured x-ray diffraction intensity may be expressed  by

$$I^{mea}(Q)=PA[N(I_{eu}^{coh}+I_{eu}^{inc}+I_{eu}^{mul})],$$
(2)
where P is the polarization factor, A the absorption factor, N the normalization constant, and $`I_{eu}^{coh},I_{eu}^{inc},I_{eu}^{mul}`$ are the coherent single scattering, incoherent (Compton), and multiple scattering intensities, respectively, per atom, in electron units. The total scattering structure function, $`S(Q)`$, is then defined as

$$S(Q)=[I_{eu}^{coh}(f^2f^2)]/f^2,$$
(3)
where $`f=f(Q)`$ is the sample average atomic form factor and $`f^2`$ is the sample average of the square of the atomic form factor. Therefore, to obtain $`S(Q)`$ from the measured diffraction data, we have to apply corrections such as multiple scattering, polarization, absorption, Compton scattering and Laue diffuse corrections on the raw data.

The corrections were carried out using a home-written computer program, PDFgetX  that is able to utilize the MCA data. The results obtained using the MCA approach are very similar to those obtained using the SCA approach.  It appears that both approaches work well for quantitative high energy x-ray powder diffraction. One possible advantage of the MCA method is that energy windows of interest can be set after the experiment is over which is precluded if data are only collected using SCA’s.

We briefly describe some of the features of the data correction using PDFgetX. Data are first corrected for detector dead-time. In this experiment, we used the pulser method.  A pulse-train from an electronic pulser of known frequency is fed into the detector preamp. The voltage of the pulser pulses are set so that the signal appears in a quiet region of the MCA spectrum. The measured counts in the pulser signal in the MCA (or, indeed, in an SCA window set on the pulser signal) is then recorded for each data point. The data dead-time correction is then obtained by scaling the raw data by the ratio of the known pulser frequency and the measured pulser counts. This method accounts for dead-time in the preamp, amplifier and MCA/SCA electronics but not in the detector itself. However, in general the dead-time is dominated by the pulse-shaping time in the amplifier or the analogue-digital conversion in the MCA or SCA and so this method gives rather accurate dynamic measurement of the detector dead-time. An alternative dead-time correction protocol for correcting MCA data is to use the MCA real-time/live-time ratio. This works reasonably well if the MCA conversion time is the dominant contribution to the detection dead-time. This approach gave similar results to the pulser correction in this case.

Multiple (mainly double) scattering can be a problem if samples are relatively thick and the radiation is highly penetrating as in the present case. The multiple scattering contribution contains no usable structural information and must be removed from the measured intensity. It depends on sample thickness and many other sample dependent factors such as attenuation coefficient, atomic number and weight of sample constituent.  It increases as the sample becomes thicker in both transmission and reflection geometry. The multiple scattering correction was calculated using the approach suggested by Warren  in the isotropic approximation. Calculation of the multiple scattering intensity is considerably simplified when the elastic and Compton signals are separated as is done here since only completely elastic multiple scattering events need to be considered. In $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$ samples, the multiple scattering ratio was around 10$`\%`$ maximum at high $`Q`$ in transmission geometry. This result suggests that the proper multiple scattering correction becomes important in the high $`Q`$ region. Fig. 2 shows the double scattering ratio calculated for the $`\mathrm{In}_{0.5}\mathrm{Ga}_{0.5}\mathrm{As}`$ sample.

The x-ray polarization correction is almost negligible for synchrotron x-ray radiation because the incident beam is almost completely plane polarized perpendicular to the scattering plane. As a result there is virtually no angle dependence to the measured intensity due to polarization effects.

The Compton scattering correction is very important in high-energy x-ray diffraction data analysis. It can become larger than the coherent scattering intensity at high $`Q`$, as is evident in Fig. 1. Even a small error in determining the Compton correction can lead to a significant error in the coherent scattering intensity in the high $`Q`$ region. However, in this region of the diffraction pattern the elastic and Compton-shifted scattering are well separated in energy and can be reliably separated using the energy resolved detection scheme we used here. At low-$`Q`$ the Compton-shift is small and the Compton and elastic signals cannot be explicitly separated unless a higher energy resolution measurement is made, for example, using an analyzer crystal. However, the Compton intensity is much lower and the coherent scattering intensity is much larger. In this region a theoretically calculated Compton signal can be subtracted from the data containing both elastic and Compton scattering. Uncertainties in this process have a very small effect on the resulting $`S(Q)`$. Fig. 3 shows the signals from the Compton and elastic

scattering in the $`\mathrm{In}_{0.5}\mathrm{Ga}_{0.5}\mathrm{As}`$ sample. At low $`Q`$, some contamination from the elastic scattering is apparent in the Compton channel. For the Compton scattering correction, we followed two steps. In the high $`Q`$ region the Compton scattered signal was directly removed by integrating a narrow region of interest in the MCA spectrum which only contained the elastic peak. In the low $`Q`$ region we calculated the theoretical Compton scattering  and subtracted this from the combined (unresolved) Compton plus elastic scattering signal. These two regions were smoothly interpolated using a window function, following the method of Ruland  in which the theoretical Compton intensity is smoothly attenuated with increasing $`Q`$.

At very high $`Q`$ values, due to the Debye-Waller factor, the Bragg peaks in the elastic scattering signal disappear and the normalized intensity asymptotes to $`f(Q)^2`$. This fact allows us to obtain an absolute data normalization by scaling the data to line up with $`f(Q)^2`$ in the high-$`Q`$ region of the diffraction pattern. Finally, the total scattering structure function, $`S(Q)`$, is then obtained using Eq. 3 and the corresponding PDFs, $`G(r)`$ are obtained according to Eq. 1.

## III Results

Fig. 4 shows the experimental reduced total scattering structure functions, $`F(Q)=Q[S(Q)1]`$, for the $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$ alloys measured at 10 K.

It is clear that the Bragg peaks are persistent up to $`Q35`$ Å-1 in the end-members, GaAs and InAs. This reflects both the long range order of the crystalline samples and the small amount of positional disorder (dynamic or static) on the atomic scale. In the alloy samples, however, the Bragg peaks disappear at much lower $`Q`$-values but still many sharp Bragg peaks are present in the mid-low $`Q`$ region. Instead, oscillating diffuse scattering which contains local structural information is evident in high $`Q`$ region. The observation of Bragg peaks reflects the presence of long-range crystalline order in these alloys. The fact that the Bragg peak intensity disappears at lower $`Q`$-values in the alloys than the end-members reflects that there is significant atomic scale disorder in the alloys as expected. The oscillating diffuse scattering in the high-$`Q`$ region originates from the stiff nearest-neighbor In-As and Ga-As covalent bonds.

Fig. 5 shows the corresponding reduced PDFs, $`G(r)`$, obtained using Eq. 1.

In the alloys, it is clear that the first peak is split into a doublet corresponding to shorter Ga-As and longer In-As bonds.  The position in $`r`$ of the left and right peaks does not disperse significantly on traversing the alloy series. This shows that the local bond lengths stay close to their end-member values and do not follow Vegard’s law, in agreement with the earlier XAFS and PDF reports. However, already by 10 Å the structure is behaving much more like the average structure. For example, the doublet of PDF peaks around 11 Å in GaAs (Fig. 5) remains a doublet (it doesn’t become a quadruplet in the alloys) and disperses smoothly across the alloy series to its position at around 12 Å in the pure InAs. This shows that already by 10 Å the structure is exhibiting Vegard’s law type behavior.

It is also notable that for the nearest neighbor PDF peak, the peak widths are almost the same in both alloys and end-members but for the higher neighbors, the peaks are much broader in the alloys than in the end-members.

## IV Atomic displacements in the alloys

### A Modeling

A simplified view of the structural disorder in the $`\mathrm{A}_\mathrm{x}\mathrm{B}_{1\mathrm{x}}\mathrm{C}`$ type tetrahedral alloys can be intuitively visualized by considering simple tetrahedral clusters centered about C sites (the unalloyed site). In the random alloy this site can have 4 A-neighbors (type-I), 3 A- and 1 B-neighbors (II), 2 A- and 2 B-neighbors (III), 1 A and 3 B-neighbors (IV) or 4 B neighbors (V). We assume that the mixed site (A,B) atoms stay on their ideal crystallographic positions. By considering each cluster type in turn we can predict the qualitative nature of the atomic displacements present in the alloy. Let the A atoms be larger than the B atoms.

In clusters of type I and V the C atom will not be displaced away from the center of the tetrahedron. As is shown in Fig. 6, in type II clusters the C atom will displace away from the center directly towards the B atom. This is a displacement in a $`111`$ crystallographic direction. In type III clusters it will displace in a direction between the two B atoms along a $`100`$ crystallographic direction. Finally, in type IV clusters it will again be a $`111`$ type displacement but this time in a direction directly away from the neighboring A atom. Such a cluster model was used to make quantitative comparisons with the nearest neighbor bond distances observed in XAFS measurements over the whole alloy series. However, it was later shown that the agreement was largely accidental due to the boundary conditions chosen that all the clusters should have the average size determined from Vegard’s law. Nonetheless, it is interesting to compare the prediction of this simple cluster model with the nearest-neighbor PDF peaks measured here since, for the first time, we have an accurate measurement of the bond length distributions as well as the bond lengths themselves.

Each cluster is independently relaxed according to the prescription of Balzarotti et al.  to get the bond-lengths within each cluster type. Assuming a random alloy the number of each type of cluster that is present can be estimated using a binomial distribution. This gives the static distribution of bond lengths predicted by the model. These are then convoluted with the broadening expected due to thermal motion. This was determined by measuring the width of the nearest neighbor peaks in the end-member compounds, InAs and GaAs. The result is shown in Fig. 7(a). It is clear that, although the

cluster model gets the peak positions reasonably correct as exemplified by the agreement it gets with XAFS data,  it does rather a poor job of explaining the shape of the measured pair distribution. The major discrepancy is that too much intensity resides at, or close to, the undisplaced position leading to an unresolved broad first PDF peak in sharp contrast to the measurement. In contrast, we show in Fig. 7(b) the nearest neighbor atomic pair distribution in the Pauling limit,  again broadened by thermal motion. The peak positions were obtained by using the bond-lengths of the end-member compounds. It is clear that this actually does a better job than the cluster model, though it slightly, and not surprisingly, overemphasizes the splitting. The dashed line in this figure shows the peak profile that we obtain if we make the assumption that the nearest neighbor bond length changes in the alloy as seen in the Z-plot,  but there is no increase in the bond length distribution. Again, this gives rather good agreement emphasizing the fact that there is very little inhomogeneous strain to the covalent bond length due to the alloying.

A better model for the structure of these alloys  is obtained from a relaxed supercell of the alloy system using a Kirkwood potential.  The potential contains nearest neighbor bond stretching force constants $`\alpha `$ and force constants $`\beta `$ that couple to the change in the angle between adjacent nearest neighbor bonds. In this relaxed supercell model, the force constants were adjusted to fit the end-members  with $`\alpha _{\mathrm{Ga}\mathrm{As}}`$ = 96N/m, $`\alpha _{\mathrm{In}\mathrm{As}}`$ = 97N/m, $`\beta _{\mathrm{Ga}\mathrm{As}\mathrm{Ga}}`$ = $`\beta _{\mathrm{As}\mathrm{Ga}\mathrm{As}}`$ = 10N/m and $`\beta _{\mathrm{In}\mathrm{As}\mathrm{In}}`$ = $`\beta _{\mathrm{As}\mathrm{In}\mathrm{As}}`$ = 6N/m. The additional angular force constants required in the alloy are taken to be the geometrical mean, so that $`\beta _{\mathrm{Ga}\mathrm{As}\mathrm{In}}`$ = $`\sqrt{(\beta _{\mathrm{Ga}\mathrm{As}\mathrm{Ga}}.\beta _{\mathrm{In}\mathrm{As}\mathrm{In}})}`$. The PDFs for the alloys could then be calculated in a self-consistent way for all the alloys with no adjustable parameters.  In this model, the lattice dynamics are also included in a completely self-consistent way. Starting with the force constants and the Kirkwood potential, the thermal broadening of the PDF peaks at any temperature can be determined directly from the dynamical matrix and this is how the PDFs were calculated in the present case. The model-PDF is plotted with the data in Fig. 8 with the nearest-neighbor peak shown on an expanded scale in Fig. 7(c). The excellent agreement with the data over the entire alloy range suggests that the simple Kirkwood potential provides an adequate starting point for calculating distorted alloy structures in these III-V alloys.

Note that in comparing with experiment, the theoretical PDF has been convoluted with a Sinc function to incorporate the truncation of the experimental data at $`Q_{max}=45`$ Å.

### B 3-D atomic probability distribution

Now, we analyze the relaxed supercell of alloy system obtained using a Kirkwood potential to get the average three dimensional atomic probability distribution of As and (In,Ga) atoms. Fig. 9 shows iso-probability surfaces for the As site in the $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$ alloy.

The probability distributions were created by translating atomic positions of the displaced arsenic atoms in the supercell ($`20\times 20\times 20`$ cubic cell) into a single unit cell. To improve statistics, this was done 70 times. The surfaces shown enclose a volume where the As atom will be found with 68 % probability. The probability distribution is viewed down the  axis. It is clear that the As atom displacements, though highly symmetric, are far from being isotropic. The same procedure has been carried out to elucidate the atomic probability distribution on the (In,Ga) sublattice. The results are shown in Fig. 10, plotted on the same scale as in Fig. 9.

In contrast to the As atom static distribution, the (In,Ga) probability distribution is much more isotropic and sharply peaked in space around the virtual crystal lattice site.

In all compositions, the As atom distribution is highly anisotropic as evident in Fig. 9 with large displacements along $`100`$ and $`111`$ directions. This can be understood easily within cluster model as we discussed in Section IV A. The $`100`$ displacements occur in type III clusters and the $`111`$ displacements occur in type II and IV clusters. This also explains why, in the gallium rich alloy in which the three and four Ga cluster is dominant, the major As atom displacements are along , \[1$`\overline{1}\overline{1}`$\], \[$`\overline{1}`$1$`\overline{1}`$\], and \[$`\overline{1}\overline{1}`$1\] as we observed in Fig. 9(a). On the contrary, in the indium rich alloy, the major displacements are along \[$`\overline{1}\overline{1}\overline{1}`$\], \[$`\overline{1}`$11\], \[1$`\overline{1}`$1\], and \[11$`\overline{1}`$\], as can be clearly seen in Fig. 9(d).

The atomic probability distribution obtained from the Kirkwood model for the (In,Ga) sublattice is shown in Fig. 10. As we discussed, this is much more isotropic (though not perfectly so), and more sharply peaked than the As atom distribution. However, contrary to earlier predictions,  and borne out quantitatively by the supercell modeling, there is significant static disorder associated with the (In,Ga) sublattice. In order to compare the magnitude of the static distortion of the (In,Ga) sublattice with that of the As sublattice, we calculated the standard deviation, $`\sigma `$, of the As and (In,Ga) atomic probability distributions. This was calculated using $`\sigma _i=\sqrt{\frac{1}{N1}_{k=1}^N(d_i(k))^2},(i=x,y,z)`$, where $`d_i`$ refers to the displacement from the undistorted sublattice of atoms in the model supercell in $`x`$, $`y`$, and $`z`$ directions, and $`N`$ is the total number of atoms in the supercell. Table I summarizes the values of $`\sigma `$ for the As and (In,Ga) atomic probability distributions in the alloys.

It shows that for all compositions the static disorder on the (In,Ga) sublattice is around 60% of the disorder on the As sublattice. These static distortions give rise to a broadening of PDF peaks as described in Ref. 27 and evident in Fig. 5 of this paper. To evaluate the static contribution to the PDF peak broadening, $`\sigma _D`$, from the $`\sigma ^{}s`$ reported in Table I we used the following expression:

$$\sigma _D^2=\sigma _a^2+\sigma _b^2,$$
(4)
where $`a`$, $`b`$ can be As, or (In,Ga). For example, for $`x`$ = 0.5 alloy, we get $`\sigma _D^2`$ = 0.0188(4) Å2 for As-As peaks in the PDF, 0.0130(4) Å2 for As-(In,Ga) peaks and 0.0072(3) Å2 for (In,Ga)-(In,Ga) peaks. These values are in good agreement with the mean square static PDF peak broadening of As-As, As-(In,Ga) and (In,Ga)-(In,Ga) peaks, shown in Fig. 4 of Ref. 27, of 0.0187(1) Å2, 0.0128(1) Å2, and 0.0053(1) Å2 respectively.

## V Correlated atomic displacements

We have shown that on the average, atomic displacements of As atoms in $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$ alloy are highly directional. In this section, we would like to address the question whether these atomic displacements are correlated from site to site. To investigate this we have calculated theoretically the diffuse scattering intensity which would be obtained from the relaxed Kirkwood supercell model and compare it with the known experimental diffuse scattering.

The Fig. 11 shows diffuse scattering of

$`\mathrm{In}_{0.5}\mathrm{Ga}_{0.5}\mathrm{As}`$ alloy calculated using the DISCUS program. In this calculation the Bragg-peak intensities have been removed. Strong diffuse scattering is evident at the Bragg points in the characteristic butterfly shape pointing towards the origin of reciprocal space. This is the Huang scattering which is peaked close to Bragg-peak positions and has already been worked out in detail.

In addition to this, clear streaks are apparent running perpendicular to the  direction. The diffuse scattering calculations on $`(hkl)`$ planes where $`l0,\mathrm{integer}`$ (Fig. 12),

show that these diffuse streaks are extended along the $`[00l]`$ direction consisting of sheets of diffuse scattering perpendicular to the  direction of reciprocal space. Diffuse scattering with exactly this (110) symmetry was observed in the TEM study of $`\mathrm{In}_{0.53}\mathrm{Ga}_{0.47}\mathrm{As}`$. Careful observation of our calculated diffuse scattering indicates that the diffuse scattering has a maximum on the low-$`Q`$ side of the ($`hk0`$) planes passing through the Bragg points, with an intensity minimum on the high-$`Q`$ side of these planes. This is characteristic size-effect scattering obtained from correlated atomic displacements due to a mismatch between chemically distinct species as recently observed in a single-crystal diffuse scattering study on Si1-xGex, for example. This asymmetric scattering was clearly observed in the earlier diffuse scattering study on $`\mathrm{In}_{0.53}\mathrm{Ga}_{0.47}\mathrm{As}`$.

The single-crystal diffuse scattering intensity which is piled up far from the Bragg-points is giving information about intermediate range ordering of the atomic displacements. It is interesting that it is piled up in planes perpendicular to  whereas the local atomic displacements are predominantly along $`100`$ and $`111`$ directions. This observation underscores the complementarity of single-crystal diffuse scattering and real-space measurements such as the PDF. The real-space measurements are mostly sensitive to the direction and magnitude of local atomic displacements and less sensitive to how the displacements are correlated over longer-range (though this information is in the data). On the other hand, the single crystal diffuse scattering immediately yields the intermediate range correlations of the displacements but one has to work harder to extract information about the size and nature of the local atomic displacements. Used together these two approaches, together with XAFS, can reveal a great deal of complementary information about the local structure of disordered materials.

The single crystal diffuse scattering suggests that atomic displacements are most strongly correlated (i.e., correlated over the longest range) along  directions although the displacements themselves occur along $`100`$ and $`111`$ directions. The reason may be that the zinc-blende crystal is stiffest along  directions because of the elastic anisotropy in the cubic crystal. This was shown for the case of InAs and was used to explain why the 5th peak in the PDF (coming from In-As next neighbor correlations along  direction) was anomalously sharp in both experiments and calculations. If the material is stiffer in this direction, one would expect that strain fields from displacements will propagate further in these directions than other directions in the crystal correlating the displacements over longer range. These are consistent with the displacement pair correlation function calculation by Glas  which shows that the correlation along $`110`$ directions is larger than correlations along $`100`$ and $`111`$ and extends further.

## VI conclusions

In conclusion, we have obtained high real-space resolution PDFs of $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$ ($`0x1`$) alloys using high energy synchrotron x-ray diffraction. For this purpose, we developed a data analysis technique adequate for high energy synchrotron x-ray diffraction. The PDFs show a clearly resolved doublet corresponding to the Ga-As and In-As bond lengths in the first peak of the alloys. Far-neighbors peaks are much broader in the alloys than that of the pure end members.

We show that As atom displacements are highly directional and can be represented as a combination of $`100`$ and $`111`$ displacements. On the contrary, the (In,Ga) atomic distribution is much more isotropic. The magnitude of (In,Ga) sublattice disorder is less than, but rather comparable ($`\sigma _{(In,Ga)}0.6\sigma _{As}`$) to, the As sublattice disorder. Also, the single crystal diffuse scattering shows that atomic displacements are correlated over the longest range in  directions although the displacements themselves occur along $`100`$ and $`111`$ directions.

All of the available data, including previous XAFS studies, the present data, differential PDF data, and diffuse scattering on a closely related system are well explained by a relaxed supercell model based on the Kirkwood potential. This study also underscores the importance of having data from complementary techniques when studying the detailed structure of crystals with significant disorder.

###### Acknowledgements.

We gratefully acknowledge M. F. Thorpe and J. S. Chung for making their supercell calculation program available and giving valuable help. We would like to acknowledge Th. Proffen for discussions about diffuse scattering in $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}\mathrm{As}`$ alloys. This work was supported by DOE through grant DE FG02 97ER45651. CHESS is supported by the National Science Foundation through grant DMR97-13424


format: text

text:
# Decay properties of the one-particle Green function in real space and imaginary time

## I Introduction

The GW approximation for the electronic self-energy is known to produce quasiparticle band structures in very good agreement with experimental data for a wide range of crystalline materials that include semiconductors, simple metals, and even transition metals. It corresponds to a summation of ring diagrams to infinite order and may be thought of as an extended Hartree-Fock scheme with dynamically screened exchange. Despite the long-recognized success for bulk systems there are still comparably few applications to more complicated geometries, however, because of the rapid increase in computational cost. Unless additional simplifications like a plasmon-pole model for the screened Coulomb interaction are employed, conventional reciprocal-space algorithms scale like $`N^4/b^2`$ for a calculation of the complete spectral function, where $`N`$ is a measure of the system size, such as the number of atoms in a bulk material, and $`b`$ is inversely proportional to the number of frequency mesh points. The scaling problem is further compounded if diagrammatic vertex corrections beyond the GW approximation are included. The design of more efficient implementations is therefore a pressing task.

In order to achieve a more favorable scaling, Rojas et al. recently proposed an algorithm based on an alternative representation of the Green function in real space and imaginary time. This approach allows the calculation of the irreducible polarizability and the self-energy as simple products rather than numerically expensive convolutions. Furthermore, it directly exploits locality, i.e., the fact that the Green function decays to zero as its spatial or temporal arguments move apart. As a result, with a fixed cutoff energy this so-called space-time method only scales like $`N^2/\gamma ^3b`$, where $`\gamma `$ denotes the exponential decay constant in real space. Logarithmic terms are neglected. Under the assumption that the frequency mesh maps one-to-one onto the relevant region of the imaginary time axis, $`b`$ can there be identified with the inverse period $`\frac{1}{2}k_BT`$ in a finite-temperature formalism and with the exponential decay constant otherwise. The locality of the Green function is a manifestation of the more general “nearsightedness” principle and due to destructive quantum interference in many-electron systems.

While the behavior of the Green function in reciprocal space and the frequency domain has been analyzed in great detail, there has so far been no comprehensive study of the decay properties in real space and, for zero temperature, on the imaginary time axis. In fact, even qualitative features are not always recognized correctly: Ref. , for instance, describes the asymptotic tail of the Green function as being proportional to $`|𝐫𝐫^{}|^2`$. While this assertion is correct for the homogeneous electron gas at zero temperature, I argue here that real materials in general exhibit exponential rather than algebraic decay. Attempts to model the self-energy using results obtained for the homogeneous electron gas, which are motivated by seemingly universal features in the short-range part of the nonlocality and the success of the local-density approximation in density-functional theory, should appreciate such fundamental differences.

The first objective of this paper is to systematically study the decay properties of the one-particle Green function in solids, separately for metals and insulators, and to obtain analytic expressions for the asymptotic behavior in real space and imaginary time. Besides intrinsic interest, these can be directly exploited in electronic-structure calculations. In the space-time method, for instance, an analytic treatment of the asymptotic tail on the imaginary time axis, in combination with a customized Gauss-Legendre grid, reduces the number of mesh points where the Green function must be evaluated numerically by one order of magnitude. I show here that the exponential fit proposed in Ref.  is only valid for systems with a finite band gap, however, while metals at zero temperature instead require an algebraic fitting function. The second objective is to relate the decay constants $`\gamma `$ and $`b`$ to known system parameters. This is then used for a more detailed scaling comparison between the space-time approach and conventional implementations.

Due to the central role of locality in $`O(N)`$ methods within density-functional theory, the spatial decay rate of the one-particle density matrix, which is a special element of the Green function, has recently received some attention. Physically, the nonlocality of the density matrix stems entirely from the delocalized quantum character of the wave functions, while the Green function also describes the actual transport of particles. The results obtained for the density matrix can therefore not be generalized in a straightforward way, although they are recovered here if the imaginary-time argument approaches zero from below. Furthermore, the one-particle density matrix contains no information about the behavior of the Green function on the imaginary time axis.

This paper is organized as follows. Sec. II introduces the Green function in real space and imaginary time. In Sec. III an analytic solution for the homogeneous electron gas is presented. I derive asymptotic formulas for real metals and insulators in Secs. IV and V. The results are summarized in Sec. VI together with a discussion. Atomic units are used throughout.

## II The Green function

The imaginary-time or thermal Green function is defined as

$`𝒢(𝐫,𝐫^{};\tau )`$ $`=`$ $`\mathrm{\Theta }(\tau )e^{(\widehat{H}\mu \widehat{N})\tau }\widehat{\psi }^{}(𝐫^{})e^{(\widehat{H}\mu \widehat{N})\tau }\widehat{\psi }(𝐫)`$ (2)
$`\mathrm{\Theta }(\tau )e^{(\widehat{H}\mu \widehat{N})\tau }\widehat{\psi }(𝐫)e^{(\widehat{H}\mu \widehat{N})\tau }\widehat{\psi }^{}(𝐫^{}),`$
where $`\widehat{\psi }^{}(𝐫)`$ and $`\widehat{\psi }(𝐫)`$ denote the creation and annihilation operator for an electron at $`𝐫`$. The Hamiltonian $`\widehat{H}`$ is modified by the chemical potential $`\mu `$ and the electron number operator $`\widehat{N}`$, the angular brackets indicate the thermodynamic average, and $`\mathrm{\Theta }(\tau )`$ denotes Heaviside’s step function. Spin degrees of freedom are suppressed to simplify the notation. At finite temperature the real variable $`\tau `$ is restricted to $`\beta <\tau <0`$ with $`\beta =1/k_BT`$ for the hole part and $`0<\tau <\beta `$ for the electron part, and the Green function is periodically repeated. Due to the additional antiperiodicity $`𝒢(\tau )=𝒢(\tau +\beta )`$ it suffices to examine the hole part. At zero temperature the period becomes infinite, however, and the properties of both parts must be studied separately. Note that Refs.  employ an alternative prefactor convention $`G(i\tau )=i𝒢(\tau )`$ for the imaginary-time Green function at zero temperature, which stems from analytic continuation in the time domain rather than the frequency domain.

For real materials it is convenient to rewrite the hole part of the Green function in the form

$$𝒢(𝐫,𝐫^{};\tau )=\underset{n}{}𝒢_n(𝐫,𝐫^{};\tau ),$$
(3)
$$𝒢_n(𝐫,𝐫^{};\tau )=\underset{𝐑,𝐑^{}}{}w_n(𝐫𝐑)F_n(𝐑𝐑^{};\tau )w_n^{}(𝐫^{}𝐑^{}),$$
(4)
$$w_n(𝐫𝐑)=\frac{1}{\mathrm{\Omega }_B}e^{i𝐤𝐑}\psi _{n𝐤}(𝐫)d^3k,$$
(5)
$$F_n(𝐑𝐑^{};\tau )=\frac{1}{\mathrm{\Omega }_B}e^{i𝐤(𝐑𝐑^{})}\frac{e^{(ϵ_{n𝐤}\mu )\tau }}{1+e^{(ϵ_{n𝐤}\mu )\beta }}d^3k$$
(6)
and analogous for the electron part. Here $`𝐑`$ labels the lattice vectors, $`w_n(𝐫𝐑)`$ indicates the Wannier orbitals corresponding to a set of Bloch wave functions $`\psi _{n𝐤}(𝐫)`$, and the integrals extend over the first Brillouin zone with volume $`\mathrm{\Omega }_B`$. Wannier orbitals are widely accepted to exhibit exponential localization. In addition to a rigorous proof for isolated bands, there is compelling numerical evidence to this effect for the general case of composite bands. Crucially, the localization persists even if defects or surfaces break the periodicity of the solid. As the decay length of Wannier orbitals is of the order of the interatomic spacing, only a small number of terms contributes significantly to the sum (4).

## III The homogeneous electron gas

For the homogeneous electron gas, the hole part of the Green function is given by

$$𝒢(r,\tau )=\frac{1}{8\pi ^3}e^{i𝐤𝐫}\frac{e^{(ϵ_k\mu )\tau }}{1+e^{(ϵ_k\mu )\beta }}d^3k$$
(7)
with $`ϵ_k=\frac{1}{2}k^2`$. The range of integration can be extended to infinity, because the spherical Fermi surface is fully contained inside the Brillouin zone. The angular integrals are straightforward and lead to

$$𝒢(r,\tau )=\frac{i}{4\pi ^2r}_{\mathrm{}}^{\mathrm{}}ke^{ikr}\frac{e^{(ϵ_k\mu )\tau }}{1+e^{(ϵ_k\mu )\beta }}𝑑k.$$
(8)
The remaining integral is most simply solved by closing the contour across the upper complex half-plane. The integrand has relevant first-order poles at $`\pm k_l+i\gamma _l`$ with

$$k_l=\sqrt{\mu +\sqrt{\mu ^2+\omega _l^2}},$$
(9)
$$\gamma _l=\sqrt{\mu +\sqrt{\mu ^2+\omega _l^2}},$$
(10)
and the fermion Matsubara frequencies $`\omega _l=(2l+1)\pi /\beta `$. Evaluation of the residues yields the final result

$$𝒢(r,\tau )=\frac{1}{\pi \beta r}\underset{l=0}{\overset{\mathrm{}}{}}e^{\gamma _lr}\mathrm{cos}(k_lr\omega _l\tau ),$$
(11)
which explicitly shows the periodicity along the imaginary time axis. At finite temperature the Green function is essentially a superposition of exponentially damped oscillations, whose long-range behavior is dominated by the $`l=0`$ contribution. For sufficiently small $`T`$ the chemical potential approximately equals the Fermi energy $`ϵ_F=\frac{1}{2}k_F^2`$. The damping constants are given by $`\gamma _l\omega _l/k_F`$ in this regime, and the wave vectors of the corresponding oscillations become $`k_lk_F`$. The overall decay constant is hence proportional to the temperature and given by $`\gamma =\gamma _0=\pi k_BT/k_F`$.

As the limit $`T0`$ is nontrivial, the Green function at zero temperature is most simply calculated from Eq. (7) with a step function instead of the Fermi distribution. The integral can be solved analytically and yields

$`𝒢(r,\tau )`$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}(k_Fr)}{2\pi ^2r\tau }}+{\displaystyle \frac{1}{2(2\pi |\tau |)^{3/2}}}\mathrm{exp}\left({\displaystyle \frac{k_F^2\tau }{2}}{\displaystyle \frac{r^2}{2\tau }}\right)`$ (13)
$`\times \left[\text{erfi}\left({\displaystyle \frac{k_F\tau +ir}{\sqrt{2|\tau |}}}\right)+\text{erfi}\left({\displaystyle \frac{k_F\tau ir}{\sqrt{2|\tau |}}}\right)\right].`$
For the electron part $`\text{erfc}(z)=1\text{erf}(z)`$ replaces $`\text{erfi}(z)=\text{erf}(iz)/i`$. The limiting formulas are identical in either case and follow from the asymptotic representation of the error function $`\text{erf}(z)`$ in the complex plane. For $`r\mathrm{}`$ the Green function becomes

$$𝒢(r,\tau )\frac{k_F\mathrm{cos}(k_Fr)}{2\pi ^2r^2}+\frac{1k_F^2\tau }{2\pi ^2r^3}\mathrm{sin}(k_Fr),$$
(14)
which decays algebraically rather than exponentially. On the imaginary time axis the Green function decays like

$$𝒢(r,\tau )\frac{\mathrm{sin}(k_Fr)}{2\pi ^2r\tau }\frac{\mathrm{cos}(k_Fr)}{2\pi ^2k_F\tau ^2}$$
(15)
for large positive or negative $`\tau `$.

## IV Metals

The explicit solution for the homogeneous electron gas illustrates more general principles that also apply to real metals. Bands above or below the Fermi level yield contributions to the Green function analogous to the case of insulators, which is discussed below. They are short-ranged, so that for $`|𝐫𝐫^{}|\mathrm{}`$ the sum (3) is dominated by terms arising from partially occupied bands that cross the Fermi level. Furthermore, the strong localization of the Wannier orbitals implies that the spatial decay properties of the relevant $`𝒢_n`$ are in turn determined by $`F_n`$. If the homogeneous electron gas is used for guidance, then the decay length $`1/\gamma `$ for sodium with an average density $`r_s=3.99`$ (Ref. ) at room temperature $`T=298\mathrm{K}`$ may be estimated as 85.8 Å. In contrast, the decay length of the temperature-independent Wannier orbitals is comparable to the interatomic spacing 3.72 Å in the bcc cell and more than an order of magnitude smaller. The high-temperature region where $`F_n`$ becomes sufficiently short-ranged such as not to dominate the asymptotic behavior is irrelevant for solid-state physics, because it lies above the melting point of the crystal lattice.

$`F_n`$ has the mathematical form of a three-dimensional Fourier transform and is long-ranged because of rapid oscillations of the integrand in Eq. (6). Quantitatively, at finite temperature the integrand falls from near unity to zero in a narrow region around the Fermi surface. Its width is given by $`\delta kk_BT/|ϵ_{n𝐤}|`$, where the gradient is taken at the Fermi level. Directional effects due to the anisotropic crystal structure are ignored here, and it is assumed that the energy dispersion is analytic and smooth on the scale of $`\delta k`$. For sufficiently small $`T`$ the latter condition is always fulfilled. According to basic Fourier analysis, the extent of this rapid variation is directly proportional to the exponential decay constant of the transform, so that $`\gamma \delta k`$ for $`F_n`$ and hence for $`𝒢`$.

At $`T=0`$ the integrand changes qualitatively, because the Fermi distribution has a discontinuity at the Fermi surface. As the Fourier transform of a discontinuous function decays algebraically rather than exponentially, both the electron and the hole part of the Green function now assume the asymptotic form $`𝒢(𝐫,𝐫^{};\tau )|𝐫𝐫^{}|^\eta `$ with $`\eta >0`$. Their decay on the imaginary time axis is likewise algebraic $`𝒢(𝐫,𝐫^{};\tau )|\tau |^\kappa `$ with $`\kappa >0`$, because the exponential $`\mathrm{exp}(ϵ_{n𝐤}\tau )`$ in the integrand reaches the Fermi level and thus exhausts the prefactor $`\mathrm{exp}(\mu \tau )`$.

## V Insulators

Temperature effects are negligible for insulators if the thermal energy $`k_BT`$ is much smaller than the band gap $`\mathrm{\Delta }`$. As this condition is almost always fulfilled below the melting point, it suffices to examine the case $`T=0`$. If no band crosses the chemical potential, then all $`F_n`$ are spatially short-ranged. For $`\tau 0`$ the integral (6) has the straightforward solution $`F_n(𝐑𝐑^{};0)=\delta _{𝐑,𝐑^{}}`$. At finite imaginary times the only significant contribution comes from the vicinity of the band edges. The energy dispersion can hence be replaced by its harmonic approximation $`|ϵ_{n𝐤}\mu |\frac{1}{2}\mathrm{\Delta }+\frac{1}{2}_{i,j}k_i(m^{})_{ij}^1k_j`$, where $`m^{}`$ denotes the effective mass tensor at the top of the valence band for the hole part and at the bottom of the conduction band for the electron part. The range of integration may now be extended to infinity, and $`F_n`$ is readily seen to be of Gaussian type. The overexponential falloff implies that for insulators the asymptotic properties of the Green function in real space are determined by the Wannier orbitals, which are exponentially localized. In the weak-binding case, where the band structure can be obtained by perturbation from that of the homogeneous electron gas, their decay rate is given by $`\gamma a\mathrm{\Delta }`$, where $`a`$ is the lattice constant. Many common semiconductors have gaps substantially smaller than their band width and are thus expected to fall into the weak-binding regime, e.g., for silicon $`\mathrm{\Delta }=1.15\mathrm{eV}`$ (Ref. ) as opposed to $`\frac{1}{2}(2\pi /a)^2=5.1\mathrm{eV}`$. The proportionality factor is direction-dependent but about unity, so that the decay length $`1/\gamma `$ may be estimated as 1.22 Å, which is of the order of the interatomic spacing 2.35 Å in silicon, as expected, and compatible with numerical calculations. In the tight-binding regime the Wannier orbitals also decay exponentially, but the dependence of $`\gamma `$ on the gap is indeterminate and subject to details of the atomic potential. This case may apply for large-gap Mott-Hubbard insulators like NiO but is less relevant in practice, because the GW approximation anyway shows serious deficiencies for such strongly correlated systems. The falloff on the imaginary time axis is exponential with a decay constant $`b=\frac{1}{2}\mathrm{\Delta }`$ both for the electron and the hole part, because the finite difference between the chemical potential and the band edges implies an incomplete cancellation of the prefactor $`\mathrm{exp}(\mu \tau )`$.

## VI Summary

In this paper I have investigated the asymptotic properties of the one-particle Green function in real space and imaginary time. An analytic solution for the homogeneous electron gas shows the Green function as a superposition of damped harmonic oscillations with a wave vector approaching $`k_F`$ as $`T0`$. The damping $`𝒢\mathrm{exp}(\gamma |𝐫𝐫^{}|)`$ is exponential at finite temperature but changes qualitatively to algebraic falloff $`𝒢|𝐫𝐫^{}|^\eta `$ at $`T=0`$. While the thermal Green function is by definition periodic on the imaginary time axis with the period $`2\beta `$, it decays like $`𝒢|\tau |^\kappa `$ when the periodicity is lifted at $`T=0`$. I have shown that these results also apply to real metals. In contrast, for insulators at zero temperature the Green function decays exponentially like $`𝒢\mathrm{exp}(\gamma |𝐫𝐫^{}|)`$ in real space and $`𝒢\mathrm{exp}(b|\tau |)`$ in imaginary time.

The limiting formulas can be exploited in electronic-structure calculations based on a space-time representation and allow an analytic treatment of the long-range tails. The asymptotic behavior of derived propagators within a perturbation scheme follows directly from the Green function. In particular, the self-energy in the GW approximation is the product of the Green function $`𝒢`$ and the screened Coulomb interaction $`W`$. In insulators screening is incomplete with $`W1/ϵ|𝐫𝐫^{}|`$ at large distances, where $`ϵ`$ is the dielectric constant, but the self-energy still decays exponentially in real space because the Green function does. On the other hand, for the homogeneous electron gas at zero temperature $`W\mathrm{cos}(2k_F|𝐫𝐫^{}|)/|𝐫𝐫^{}|^3`$, so that the self-energy only decays algebraically. This distinct asymptotic behavior should be appreciated if the electron-gas self-energy is applied to real materials. The problem is naturally solved by the ad hoc introduction of a gap in the otherwise metallic spectrum.

Finally, the generic dependence $`\gamma k_BT/|ϵ_{n𝐤}|`$ and $`b=\frac{1}{2}k_BT`$ for metals as well as $`\gamma a\mathrm{\Delta }`$ and $`b=\frac{1}{2}\mathrm{\Delta }`$ for insulators can be used to assess the scaling of different algorithms for self-energy calculations in the GW approximation. If all other parameters remain fixed, then conventional reciprocal-space implementations scale like $`N^4/T^2`$ for metals at finite temperature and $`N^4/\mathrm{\Delta }^2`$ for insulators at zero temperature, whereas the space-time method scales like $`N^2/T^4`$ and $`N^2/\mathrm{\Delta }^4`$, respectively. This difference in efficiency gain as the temperature or the band gap increase is an important factor to be taken into account for benchmarking purposes.


text:
# Looking for a varying 𝛼 in the Cosmic Microwave Background

## I Introduction

There has been a recent growth of interest in theories where some of the usual constants of nature are actually time- and/or space-varying quantities. Most notably, the possibility of a time-varying fine-structure $`\alpha `$, has been the subject of a considerable amount of work, both at the theoretical and experimental/observational level.

From the theoretical point of view, the motivation comes from the recent work on the higher-dimensional theories , which are thought to be required to provide a consistent unification of the know fundamental interactions. In such theories the ‘effective’ three-dimensional constants are typically related to the ‘true’ higher-dimensional constants via the radii of the (compact) extra dimensions . On the other hand, these radii often have a non-trivial evolution, naturally leading to the expectation of time (or even space) variations of the ‘effective’ coupling constants we can measure .

There are a number of different ways in which a variation of $`\alpha `$ can be modelled. From a ‘theoretical’ point of view, the more convenient one appears to be to interpret it as a variation in the speed of light $`c`$ , but other alternatives have been explored . It is also possible to analyse the consequences of the variation of $`\alpha `$ in a more phenomenological context, as was done in .

On the observational level, the situation is at present somewhat confusing—see  for a brief summary. The best limit from laboratory experiments (using atomic clocks) is

$$|\dot{\alpha }/\alpha |<3.7\times 10^{14}yr^1.$$
(1)
Measurements of isotope ratios in the Oklo natural reactor provide the strongest geophysical constraints ,

$$|\dot{\alpha }/\alpha |<0.7\times 10^{16}yr^1,$$
(2)
although there are suggestions  that due to a number of nuclear physics uncertainties and model dependencies a more realistic bound is $`|\dot{\alpha }/\alpha |<5\times 10^{15}yr^1`$. Note that these measurements effectively probe timescales corresponding to a cosmological redshift of about $`z0.1`$ (compare with astrophysical measurements below).

Three kinds of astrophysical tests have been used. Firstly, big bang nucleosynthesis  can in principle provide rather strong constraints at very high redshifts, but it has a strong drawback in that one is always forced to make an assumption on how the neutron to proton mass difference depends on $`\alpha `$. This is needed to estimate the effect of a varying $`\alpha `$ on the $`{}_{}{}^{4}He`$ abundance. The abundances of the other light elements depend much less strongly on this assumption, but on the other hand these abundances are much less well known observationally. Hence one can only find the relatively weak bound

$$|\mathrm{\Delta }\alpha /\alpha |<2\times 10^2,z10^910^{10}.$$
(3)
Secondly, observations of the fine splitting of quasar doublet absorption lines probe smaller redshifts, but should be much more reliable. Unfortunately, the two groups which have been actively studying this topic report different results. Webb and collaborators  were the first to report a positive result,

$$\mathrm{\Delta }\alpha /\alpha =(1.9\pm 0.5)\times 10^5,z1.01.6$$
(4)
Note that this means that $`\alpha `$ was smaller in the past. Recently the same group reports two more (as yet unpublished) positive results , $`\mathrm{\Delta }\alpha /\alpha =(0.75\pm 0.23)\times 10^5`$ for redshifts $`z0.61.6`$ and $`\mathrm{\Delta }\alpha /\alpha =(0.74\pm 0.28)\times 10^5`$ for redshifts $`z1.62.6`$. On the other hand, Varshalovich and collaborators  report only a null result,

$$\mathrm{\Delta }\alpha /\alpha =(4.6\pm 4.3\pm 1.4)\times 10^5,z24;$$
(5)
the first error bar corresponds to the statistical error while the second is the systematic one. This corresponds to the bound

$$|\dot{\alpha }/\alpha |<1.4\times 10^{14}yr^1$$
(6)
over a timescale of about $`10^{10}`$ years. It should be emphasised that the observational techniques used by both groups have significant differences, and it is presently not clear how the two compare when it comes to eliminating possible sources of systematic error.

Finally, a third option is the cosmic microwave background (CMB) . This probes intermediate redshifts, but has the significant advantage that one has (or will soon have) highly accurate data.

The reason why the Cosmic Microwave Background is a good probe of variations of the fine-structure constant is that these alter the ionisation history of the universe . The dominant effect is a change in the redshift of recombination, due to a shift in the energy levels (and, in particular, the binding energy) of Hydrogen. The Thomson scattering cross-section is also changed for all particles, being proportional to $`\alpha ^2`$. A smaller effect (which has so far been neglected) is expected to come from a change in the Helium abundance.

As is well known, CMB fluctuations are typically described in terms of spherical harmonics,

$$T(\theta ,\varphi )=\underset{\mathrm{}m}{}a_\mathrm{}mY_\mathrm{}m(\theta ,\varphi ),$$
(7)
from whose coefficients one defines

$$C_{\mathrm{}}=<|a_\mathrm{}m|^2>.$$
(8)
Increasing $`\alpha `$ increases the redshift of last-scattering, which corresponds to a smaller sound horizon. Since the position of the first Doppler peak (which we shall denote as $`\mathrm{}_{peak}`$) is inversely proportional to the sound horizon at last scattering, we see that increasing $`\alpha `$ will produce a larger $`\mathrm{}_{peak}`$ . This larger redshift of last scattering also has the additional effect of producing a larger early ISW effect, and hence a larger amplitude of the first Doppler peak . Finally, an increase in $`\alpha `$ decreases the high-$`\mathrm{}`$ diffusion damping (which is essentially due to the finite thickness of the last-scattering surface), and thus increases the power on very small scales.

The authors of  provide an analysis of these effects and conclude that future CMB experiments should be able to provide constraints on a varying $`\alpha `$ at the recombination epoch (that is, at redshifts $`z1000`$) at the level of

$$|\dot{\alpha }/\alpha |<7\times 10^{13}yr^1,$$
(9)
or equivalently

$$|\alpha ^1d\alpha /dz|<9\times 10^5,$$
(10)
which seems to indicate that these constraints can only become competitive in the near future.

Here we analyse these effects for the BOOMERanG  and MAXIMA  data. We briefly review the method and then discuss the results in the next section. We find that this data tends to prefer a value of $`\alpha `$ that was lower in the past. However, we strongly emphasise that there are interesting and so far unnoticed degeneracies in the physics of the problem which imply that this method of determining the fine-structure constant can only produce strong constraints if other cosmological parameters are independently known. We will comment on this point in section III.

While this paper was being finalised, another preprint appeared , containing an independent analysis of the same data. It should be noticed that there are some significant differences in the two analysis procedures, as well as in the results, which we will point out along the way. In the cases where a direct comparison is possible, our work confirms their results, while in the other cases we provide some physical motivation for the differences.

## II Data analysis

We perform a likelihood analysis of the recently released BOOMERanG  and MAXIMA  data, allowing for the possibility of a time-varying fine-structure constant. The method used follows the procedure described in . The angular power spectrum $`C_l`$ was obtained using a modified CMBFAST algorithm which allows a varying $`\alpha `$ parameter. We have changed the subroutine RECFAST  according to the extensive description given in .

We vary the power spectrum normalisation $`C_2`$ within the $`95\%`$ limits for the COBE 4-year data  The space of model parameters spans

$$\mathrm{\Omega }_0=(0.11.0),$$
(11)
$$\mathrm{\Omega }_bh^2=(0.010.028),$$
(12)
$$H_0=(5080),$$
(13)
$$\alpha /\alpha _0=(0.91.1),$$
(14)
and the $`C_l`$ normalisation

$$bias=C_2/C_{2,COBE}=(0.831.16).$$
(15)
Note that $`\alpha _0`$ is the value of the fine structure constant today. The basic grid of models was obtained considering parameter step sizes of 0.1 for $`\mathrm{\Omega }_0`$; 0.003 for $`\mathrm{\Omega }_bh^2`$; 5 for $`H_0`$; 0.01 for $`\alpha /\alpha _0`$ and finally 0.01 for the $`bias`$. In order to compute the maxima and the confidence intervals for the 1-dim marginalised distributions we have increased the grid resolution of each of the model parameters using interpolation procedures.

All our models have $`\mathrm{\Omega }_{total}=1`$ and no tilt. We point out that this is in agreement  with the best-fit model for the $`\mathrm{\Omega }_{total}=1`$ case for the combined analysis of the BOOMERanG and MAXIMA data. Somewhat surprisingly, the authors of  seem to find that the same data prefer tilted models even in the ‘standard’ case, that is without considering a varying fine-structure constant. This will obviously affect their results, as we will discuss below.

We also emphasise that $`\alpha `$ may or may not have had a different value at the nucleosynthesis epoch, which in particular will affect the value of $`\mathrm{\Omega }_bh^2`$. Since in the present paper we do not treat this effect, the correct approach  is to accept the observationally determined value. This is the reason why we ignore excessively high values of $`\mathrm{\Omega }_bh^2`$ which seem to be required by some analyses of the CMB data .

For each of these $`C_l`$ spectra we compute the flat band power estimates of the CMB anisotropies obtaining a simulated observation for each data point. These estimates cover the multi-pole range sampled by both BOOMERanG and MAXIMA experiments. The values of $`\mathrm{\Omega }_b`$ are obtained as function of $`H_0`$ in order to satisfy the range for $`\mathrm{\Omega }_bh^2`$ defined above.

### A The full dataset

Interesting conclusions can be drawn by comparing the nominal calibration case for both experiments with the Likelihood obtained after marginalising over the calibration uncertainties of BOOMERanG ($`20\%`$) and MAXIMA ($`8\%`$). As mentioned in  the fitting to a lower second Doppler peak can be achieved either by increasing the value of $`\mathrm{\Omega }_bh^2`$ or decreasing the value of $`\alpha /\alpha _0`$. Note also that there is an additional constraint on $`\alpha /\alpha _0`$ coming from the position of the main acoustic peak . We shall comment on the relative importance of the two constraints below.

For the nominal calibration case we obtain a best fit model with

$$H_0=55,\mathrm{\Omega }_bh^2=0.025,\mathrm{\Omega }_0=0.5,$$
(16)
$$\alpha /\alpha _0=0.94,bias=0.87,\chi ^2=20.38.$$
(17)
In Fig 1 we plot the marginalised distributions for all the parameters.

For this case the maxima and confidence intervals for the marginalised distributions are as follows (all $`1\sigma `$): $`\alpha /\alpha _0=0.96_{0.03}^{+0.02}`$, $`H_0=50_{0.00}^{+16.02}`$, $`\mathrm{\Omega }_bh^2=0.026_{0.001}^{+0.002}`$, $`\mathrm{\Omega }_0=0.3_{0.1}^{+0.3}`$; $`bias=0.9\pm 0.04`$.

We then marginalised the 5-dim Likelihood function over the calibration uncertainties assuming a Gaussian prior. In this case we obtain a best fit model with

$$H_0=55,\mathrm{\Omega }_bh^2=0.025,\mathrm{\Omega }_0=0.4;$$
(18)
$$\alpha /\alpha _0=0.93,bias=0.91,\chi ^2=15.55.$$
(19)
In Fig 2 we plot the marginalised distributions for all the parameters for the distribution marginalised over the calibration errors. Comparing with Fig 1 we notice that in the case of the marginalised distribution the distributions for $`\mathrm{\Omega }_0`$ and $`\alpha /\alpha _0`$ are slightly shifted towards lower values, while those of $`\mathrm{\Omega }_bh^2`$ and $`H_0`$ are not significantly affected.

The maxima and confidence intervals for the marginalised distributions are as follows (again, all are $`1\sigma `$): $`\alpha /\alpha _0=0.94_{0.02}^{+0.03}`$, $`H_0=50_{0.0}^{+16.62}`$; $`\mathrm{\Omega }_bh^2=0.026\pm 0.002`$, $`\mathrm{\Omega }_0=0.3_{0.14}^{+0.17}`$, $`bias=0.94\pm 0.07`$.

If we consider the best calibration case assuming a uniform prior we get the same best fit set of parameter values as for the calibration marginalised likelihood, apart from the bias (which now has the value $`bias=0.84`$). The best calibration factor is 1.0 (nominal) for BOOMERanG and 0.92 (ratio with respect to the nominal case) for MAXIMA (this corresponds to $`\chi ^2=15.73`$). This is just telling us that if we keep BOOMERanG at the nominal calibration case and lower the height of the MAXIMA data points we force the normalisation of the models to decrease.

If instead we consider the best calibration case assuming a Gaussian prior we get a best fit with a calibration factor of 1.1 for BOOMERanG and 1.0 (nominal) for MAXIMA ($`\chi ^2=15.75`$ before weighting). As should be expected we get a higher best fit value for the models normalisation of $`bias=0.92`$. Meanwhile for both cases we observe a decrease on the value of $`\mathrm{\Omega }_0`$ from 0.5 to 0.4 and on the value of $`\alpha /\alpha _0`$ from 0.94 to 0.93, when compared with the nominal case.

Note that both pushing the BOOMERanG data up or pushing the MAXIMA data down provide for a better overlap of the two data sets. Then the different overall normalisations (in particular the height of the first Doppler peak) account for the different values of the cosmological parameters.

In fig 3 we plot the 2-dim Likelihood functions obtained after marginalising over the remaining three parameters. In Fig 4 we plot the likelihood surface for $`H_0`$ and $`\alpha /\alpha _0`$ as well as for $`\mathrm{\Omega }_bh^2`$ and $`H_0`$. Similarly, Fig 5 contains the corresponding likelihoods for the nominal calibration case. This highlights the fact that there are some non-trivial degeneracies in the problem . We shall return to this point below.

It is of interest to investigate the case where no variation of the fine structure constant is allowed. For that purpose we considered the conditional distribution for $`\alpha /\alpha _0=1.0`$ to obtain a best fit model with

$$H_0=75,\mathrm{\Omega }_bh^2=0.028,\mathrm{\Omega }_0=0.3;$$
(20)
$$bias=1.0,\chi ^2=17.49.$$
(21)
In Fig 6 we plot this distribution marginalised over $`\mathrm{\Omega }_0`$ and the bias. Increasing the value of $`\alpha `$ seems to force a higher best fit value of $`H_0`$ and of $`\mathrm{\Omega }_bh^2`$ and a lower value of $`\mathrm{\Omega }_0`$ with a best fit COBE normalised model. Again, this is consistent with  (which also find a tilt $`n_s=0.99`$, while  find $`n_s=0.92`$).

If we instead condition our distribution to a value of $`\mathrm{\Omega }_bh^2=0.019`$ we get a best fit model with

$$H_0=50,\mathrm{\Omega }_0=0.3,\alpha /\alpha _0=0.9;$$
(22)
$$bias=0.89,\chi ^2=18.23.$$
(23)
This result emphasises the rather obvious point that reducing the value of $`\mathrm{\Omega }_bh^2`$ requires a lower value of $`\alpha /\alpha _0`$ to account for a low second acoustic peak.

Once more, we emphasise that the values of $`\mathrm{\Omega }_b`$ are obtained as function of $`H_0`$ in order to satisfy the range for $`\mathrm{\Omega }_bh^2`$ defined above. This means that we should expect to observe some correlation when plotting $`H_0`$ against $`\mathrm{\Omega }_bh^2`$. This is indeed confirmed by Fig 3.

Therefore we have so far confirmed the fact that to fit a low second Doppler peak we need a high baryonic content  and a lower fine structure constant in the past. However, an important question still remains: what’s the weight of the second acoustic peak relative to the main peak in drawing the above conclusions? We recall that in  it was shown that the position of the first Doppler peak can by itself provide a constraint on $`\alpha `$. Can it happen that the main acoustic peak is still a heavy factor in determining the above best fit parameters?

### B The first Doppler peak

In order to answer this important question, we considered the set of data points sampling the multi-pole region to up $`l400`$. Hence this data set consists now on the first 8 Boomerang and the first 5 Maxima data points with the nominal calibration ($`BM_{dp}`$). We then repeat the likelihood analysis for this new data set.

We obtain a best fit model with

$$H_0=50,\mathrm{\Omega }_bh^2=0.019,\mathrm{\Omega }_0=0.4,$$
(24)
$$\alpha /\alpha _0=0.9,bias=0.87,\chi ^2=14.91.$$
(25)
This is rather encouraging, particularly because the best-fit value for $`\mathrm{\Omega }_bh^2`$ is precisely the one found by observations. However, the maxima and confidence intervals for the marginalised distributions are as follows (again, all are $`1\sigma `$): $`\alpha /\alpha _0=0.97\pm 0.04`$, $`H_0=71.4_{13.5}^{+6.1}`$; $`\mathrm{\Omega }_bh^2=0.024\pm 0.003`$, $`\mathrm{\Omega }_0=0.4_{0.21}^{+0.37}`$, $`bias=1.0_{0.08}^{+0.05}`$.

We conclude that most of the best fit model parameters do not lay within the $`1\sigma `$ range around the maximum of the marginalised distribution for the corresponding parameter. Therefore the 5-dim likelihood must have a narrow peak around this best model with an enlarged surface around the remaining values which height is not significantly smaller then the absolute peak.

In Fig 7 we plot the marginalised distributions for all the parameters, while in fig 8 we plot the 2-dim Likelihood functions obtained after marginalising over the remaining three parameters.

Comparing Fig 7 with Fig 1 we immediately notice a number of extremely interesting points. Firstly, even though the dataset for the first Doppler peak favours a smaller $`\alpha `$ in the past, there is a non-negligible likelihood for larger values as well. The inclusion of information from the second Doppler peak all but eliminates this possibility (compare this with Fig. 4 of ). Secondly, the full dataset increases the preferred values of $`\mathrm{\Omega }_bh^2`$ (or more accurately, decreases the probability for low values). And thirdly, data from the first Doppler peak alone is basically insensitive to $`H_0`$ and $`\mathrm{\Omega }_0`$ (recall that all our models have $`\mathrm{\Omega }_{total}=1`$), while the full dataset tends to favour low values of $`H_0`$ and also narrows the distribution for $`\mathrm{\Omega }_0`$ around a value of $`\mathrm{\Omega }_0=0.3`$ (and most notably reduces the probability of lower values such as $`\mathrm{\Omega }_0=0.1`$ which would be allowed by the reduced dataset).

This might explain the differences in the contour plots of the 2-dim distribution of ($`\alpha /\alpha _0`$, $`\mathrm{\Omega }_0`$) in Fig 8 and Fig 5; the plot in Fig 8 shows a correlation between ($`\alpha /\alpha _0`$, $`\mathrm{\Omega }_0`$) which disappears when including the other Doppler peaks. These 2-dim plots do also indicate correlations between ($`\alpha /\alpha _0`$, $`H_0`$); ($`\alpha /\alpha _0`$,$`\mathrm{\Omega }_bh^2`$) and ($`\alpha /\alpha _0`$,bias) which do exist in both situations.

Finally, we also consider the case with no variation of the fine structure constant allowed for the reduced dataset. We obtain a best fit model with

$$H_0=75,\mathrm{\Omega }_bh^2=0.025,\mathrm{\Omega }_0=0.5;$$
(26)
$$bias=1.02,\chi ^2=15.24.$$
(27)
If now we condition our distribution to a value of $`\mathrm{\Omega }_bh^2=0.019`$ we get a best fit model with

$$H_0=50,\mathrm{\Omega }_0=0.4,\alpha /\alpha _0=0.9;$$
(28)
$$bias=0.87,\chi ^2=14.91.$$
(29)
In Fig 9 we plot these conditional distributions marginalised over $`\mathrm{\Omega }_0`$ and the bias.

## III Discussion and conclusions

In this paper we have performed a likelihood analysis of the combined BOOMERanG and MAXIMA datasets, allowing for the possibility of a time-varying fine-structure constant, for which there is further observational evidence elsewhere . We have confirmed the intuitively obvious expectation that this data prefers a value of $`\alpha `$ that was smaller in the past by a few per cent.

However, we wish to emphasise that this is not the same as saying that the CMB can readily provide and unambiguous measurement of the fine-structure constant. As we hopefully made clear above, there are some interesting degeneracies in the problem which imply that other cosmological parameters could still fairly easily mimic a varying $`\alpha `$. Hence this method of measurement of $`\alpha `$ is still far from being ‘competitive’, in the sense that statements about $`\alpha `$ will not be possible until independent accurate determinations of $`\mathrm{\Omega }_bh^2`$ and $`H_0`$ (and possibly other parameters) are available.

We have also shown that the main reason behind the preferred lower value of $`\alpha `$ in the present dataset still comes mainly from the data points around the first Doppler peak. The main effect of the high-$`\mathrm{}`$ data points is to increase the preferred value for $`\mathrm{\Omega }_bh^2`$ and eliminate the possibility of a larger fine-structure constant in the past. A secondary (from this perspective) effect of the small angular scale data is to tighten the constraints on other parameters. Furthermore, we believe that this relative dominance of the low-$`\mathrm{}`$ measurements will remain even in the post-MAP era.

###### Acknowledgements.

We thank Bruce Bassett, Pedro Ferreira, João Magueijo, Anupam Mazumdar and Paul Shellard for useful discussions and comments. C.M. and G.R. are funded by FCT (Portugal) under ‘Programa PRAXIS XXI’, grant no. PRAXIS XXI/BPD/11769/97 and PRAXIS XXI/BPD/9990/96, respectively. We thank Centro de Astrofísica da Universidade do Porto (CAUP) for the facilities provided. GR also thanks the Dept. of Physics of the University of Oxford for support and hospitality during the progression of this work.