text:
# 1 INTRODUCTION

## 1 INTRODUCTION

Many measurements at the Tevatron rely on parton distribution functions (PDFs) for significant portions of their data analysis as well as the interpretation of their results. For example, in cross section measurements the acceptance calculation often relies on Monte Carlo (MC) estimates of the fraction of unobserved events. As another example, the measurement of the mass of the $`W`$ boson depends on PDFs via the modeling of the production of the vector boson in MC. In such cases, uncertainties in the PDFs contribute, by necessity, to uncertainties on the measured quantities. Critical comparisons between experimental data and the underlying theory are often even more dependent upon the uncertainties in PDFs. The uncertainties on the production cross sections for $`W`$ and $`Z`$ bosons, currently limited by the uncertainty on the measured luminosity, are approximately 4%. At this precision, any comparison with the theoretical prediction inevitably raises the question: How “certain” is the prediction itself?

A recent example of the importance of PDF uncertainty is the proper interpretation of the measurement of the high-$`E_T`$ jet cross-section at the Tevatron. When the first CDF measurement was published , there was a great deal of controversy over whether the observed excess, compared to theory, could be explained by deviations of the PDFs, especially the gluon, from the conventionally assumed behavior, or could it be the first signal for some new physics .

With the unprecedented precision and reach of many of the Run I measurements, understanding the implications of uncertainties in the PDFs has become a burning issue. During Run II (and later at LHC) this issue may strongly affect the uncertainty estimates in precision Standard Model studies, such as the all important $`W`$-mass measurement, as well as the signal and background estimates in searches for new physics.

In principle, it is the uncertainties on physical quantities due to parton distributions, rather than on the PDFs themselves, that is of primary concern. The latter are theoretical constructs which depend on the renormalization and factorization schemes; and there are strong correlations between PDFs of different flavors and from different values of $`x`$, which can compensate each other in the convolution integrals that relate them to physical cross-sections. On the other hand, since PDFs are universal, if we can obtain meaningful estimates of their uncertainties based on analysis of existing data, then the results can be applied to all processes that are of interest in the future.

One can attempt to assess directly the uncertainty on a specific physical prediction due to the full range of PDFs allowed by available experimental constraints. This approach will provide a more reliable estimate for the range of possible predictions for the physical variable under study, and may be the best course of action for ultra-precise measurements such as the mass of the $`W`$ boson or the $`W`$ production cross-section. However, such results are process-specific and therefore the analysis must be carried out for each case individually.

Until recently, the attempts to quantify either the uncertainties on the PDFs themselves (via uncertainties on their functional parameters, for instance) or the uncertainty on derived quantities due to variations in the PDFs have been rather unsatisfactory. Two commonly used methods are: (1) Comparing the predictions obtained with different PDF sets, e.g., various CTEQ , MRS  and GRV  sets; (2)

Within a given global analysis effort, varying individual functional parameters ad hoc, within limits considered to be consistent with the existing data, e.g. . Neither method provides a systematic, quantitative measure of the uncertainties of the PDFs or their predictions.

As a case in point, Fig. 1 shows how the calculated value of the cross section for $`W`$ boson production at the Tevatron varies with a set of historical CTEQ PDFs as well as the most recent CTEQ  and MRST  sets. Also shown are the most recent measurements from DØ and CDF111It is interesting to note that much of the difference between the DØ and CDF $`W`$ cross sections is due to the different values of the total $`p\overline{p}`$ cross sections used.

While it is comforting to see that the predictions have remained within a narrow range, the variation observed cannot be characterized as a meaningful estimate of the uncertainty: (i) the variation with time reflects mostly the changes in experimental input to, or analysis procedure of, the global analyses; and (ii) the perfect agreement between the values of the most recent CTEQ5M1222 CTEQ5M1 is an updated version of CTEQ5M differing only in a slight improvement in the QCD evolution (cf. note added in proof of ). The differences are completely insignificant for our purposes. Henceforth, we shall refer to them generically as CTEQ5M. Both sets can be obtained from the web address http://cteq.org/. and MRS99 sets must be fortuitous, since each group has also obtained other satisfactory sets which give rise to much larger variations of the $`W`$ cross section. The MRST group, in particular has examined the range of this variation by setting a variety of parameters to some extreme values . These studies are useful but can not be considered quantitative or definitive. What is needed are methods that explore thoroughly the possible variations of the parton distribution functions.

It is important to recognize all potential sources of uncertainty in the determination of PDFs. Focusing on some of these, while neglecting significant others, may not yield practically useful results. Sources of uncertainty are listed below:

* Statistical uncertainties of the experimental data used to determine the PDFs. These vary over a wide range among the experiments used in a global analysis, but are straightforward to treat.
* Systematic uncertainties  within each data set. There are typically many sources of experimental systematic uncertainty, some of which are highly correlated. These uncertainties can be treated by standard methods of probability theory *provided* they are precisely known, which unfortunately is often not the case – either because they may not be randomly distributed and/or because their estimation in practice involves subjective judgements.
* Theoretical uncertainties arising from higher-order PQCD corrections, resummation corrections near the boundaries of phase space, power-law (higher twist) and nuclear target corrections, etc.
* Uncertainties due to the  parametrization of the non-perturbative PDFs, $`f_a(x,Q_0)`$, at some low momentum scale $`Q_0`$. The specific choice of the functional form used at $`Q_0`$ introduces implicit correlations between the various $`x`$-ranges, which could be as important, if not more so, than the experimental correlations in the determination of $`f_a(x,Q)`$ for all $`Q`$.

Since strict quantitative statistical methods are based on idealized assumptions, such as random measurement uncertainties, an important trade-off must be faced in devising a strategy for the analysis of PDF uncertainties. If emphasis is put on the “rigor” of the statistical method, then many important experiments cannot be included in the analysis, either because the published errors appear to fail strict statistical tests or because data from different experiments appear to be mutually exclusive in the parton distribution parameter space . If priority is placed on using the maximal experimental constraints from available data, then standard statistical methods may not apply, but must be supplemented by physical considerations, taking into account experimental and theoretical limitations. We choose the latter tack, pursuing the determination of the uncertainties in the context of the current CTEQ global analysis. In particular, we include the same body of the world’s data as constraints in our uncertainty study as that used in the CTEQ5 analysis; and adopt the “best fit” – the CTEQ5M1 set – as the base set around which the uncertainty studies are performed. In practice, there are unavoidable choices (and compromises) that must be made in the analysis. (Similar subjective judgements often are also necessary in estimating certain systematic errors in experimental analyses.) The most important consideration is that quantitative results must remain robust with respect to reasonable variations in these choices.

In this Report we describe preliminary results obtained by our group using the two approaches mentioned earlier. In Section 3 we focus on the error matrix, which characterizes the general uncertainties of the non-perturbative PDF parameters. In Sections 4 and 5 we study specifically the production cross section $`\sigma _W`$ for $`W^\pm `$ bosons at the Tevatron, to estimate the uncertainty of the prediction of $`\sigma _W`$ due to PDF uncertainty. We start in Section 2 with a review of some aspects of the CTEQ global analysis on which this study is based.

## 2 Elements of the Base Global Analysis

Since our strategy is based on using the existing framework of the CTEQ global analysis, it is useful to review some of its features pertinent to the current study .

### Data selection:

Table 1 shows the experimental data sets included in the CTEQ5 global analysis, and in the current study. For neutral current DIS data only the most accurate proton and deuteron target measurements are kept, since they are the “cleanest” and they are already extremely extensive. For charged current (neutrino) DIS data, the significant ones all involve a heavy (Fe) target. Since these data are crucial for the determination of the normalization of the gluon distribution (indirectly via the momentum sum rule), and for quark flavor differentiation (in conjunction with the neutral current data), they play an important role in any comprehensive global analysis. For this purpose, a heavy-target correction is applied to the data, based on measured ratios for heavy-to-light targets from NMC and other experiments. Direct photon production data are not included because of serious theoretical uncertainties, as well as possible inconsistencies between existing experiments. Cf.  and . The combination of neutral and charged DIS, lepton-pair production, lepton charge asymmetry, and inclusive large-$`p_T`$ jet production processes provides a fairly tightly constrained system for the global analysis of PDFs. In total, there are $``$1300 data points which meet the minimum momentum scale cuts which must be imposed to ensure that PQCD applies. The fractional uncertainties on these points are distributed roughly like $`dF/F`$ over the range $`F=0.0030.4`$.

### Parametrization:

The non-perturbative parton distribution functions $`f_a(x,Q)`$ at a low momentum scale $`Q=Q_0`$ are parametrized by a set of functions of $`x`$, corresponding to the various flavors $`a`$. For this analysis, $`Q_0`$ is taken to be 1 GeV. The specific functional forms and the choice of $`Q_0`$ are not important, as long as the parametrization is general enough to accommodate the behavior of the true (but unknown) non-perturbative PDFs. The CTEQ analysis adopts the functional form

$$a_0x^{a_1}(1x)^{a_2}(1+a_3x^{a_4}).$$

for most quark flavors as well as for the gluon.333 An exception is that recent data from E866 seem to require the ratio $`\overline{d}/\overline{u}`$ to take a more unconventional functional form. After momentum and quark number sum rules are enforced, there are 18 free parameters left over, hereafter referred to as “shape parameters” $`\{a_i\}`$. The PDFs at $`Q>Q_0`$ are determined from $`f_a(x,Q_0)`$ by evolution equations from the renormalization group.

### Fitting:

The values of $`\{a_i\}`$ are determined by fitting the global experimental data to the theoretical expressions which depend on these parameters. The fitting is done by minimizing a global “chi-square” function, $`\chi _{\mathrm{global}}^2`$. The quotation mark indicates that this function serves as a *figure of merit* of the quality of the global fit; it does not necessarily have the full significance associated with rigorous statistical analysis, for reasons to be discussed extensively throughout the rest of this report. In practice, this function is defined as:

$`\chi _{\mathrm{global}}^2`$ $`=`$ $`{\displaystyle \underset{n}{}}{\displaystyle \underset{i}{}}w_n\left[\left(N_nd_{ni}t_{ni}\right)/\sigma _{ni}^d\right]^2`$ (1)
$`+`$ $`{\displaystyle \underset{n}{}}\left[\left(1N_n\right)/\sigma _n^N\right]^2`$
where $`d_{ni}`$, $`\sigma _{ni}^d`$, and $`t_{ni}`$ denote the data, measurement uncertainty, and theoretical value (dependent on $`\{a_i\}`$) for the $`i^{\mathrm{th}}`$ data point in the $`n^{\mathrm{th}}`$ experiment. The second term allows the absolute normalization ($`N_n`$) for each experiment to vary, constrained by the published normalization uncertainty ($`\sigma _n^N`$). The $`w_n`$ factors are weights applied to some critical experiments with very few data points, which are known (from physics considerations) to provide useful constraints on certain unique features of PDFs not afforded by other experiments. Experience shows that without some judiciously chosen weights, these experimental data points will have no influence in the global fitting process. The use of these weighing factors, to enable the relevant unique constraints, amounts to imposing certain prior probability (based on physics knowledge) to the statistical analysis.

In the above form, $`\chi _{\mathrm{global}}^2`$ includes for each data point the random statistical uncertainties and the combined systematic uncertainties in uncorrelated form, as presented by most experiments in the published papers. These two uncertainties are *combined in quadrature* to form $`\sigma _{ni}^d`$ in Eq. 1. Detailed point to point correlated systematic uncertainties are not available in the literature in general; however, in some cases, they can be obtained from the experimental groups. For global fitting, uniformity in procedure with respect to all experiments favors the usual practice of merging them into the uncorrelated uncertainties. For the study of PDF uncertainties, we shall discuss this issue in more detail in Section 5.

### Goodness-of-fit for CTEQ5M:

Without going into details, Fig. 2 gives an overview of how well CTEQ5m fits the total data set. The graph is a histogram of the variable $`x(dt)/\sigma `$ where $`d`$ is a data value, $`\sigma `$ the uncertainty of that measurement (statistical and systematic combined), and $`t`$ the theoretical value for CTEQ5m. The curve in Fig. 2 has no adjustable parameters; it is the Gaussian with width 1 normalized to the total number of data points (1295). Over the entire data set, the theory fits the data within the assigned uncertainties $`\sigma _{ni}^d`$, indicating that those uncertainties are numerically consistent with the actual measurement fluctuations. Similar histograms for the individual experiments reveal various deviations from the theory, but globally the data have a reasonable Gaussian distribution around CTEQ5M.

## 3 Uncertainties on PDF parameters: The Error Matrix

We now describe results from an investigation of the behavior of the $`\chi _{\mathrm{global}}^2`$ function at its minimum, using the standard error matrix approach . This allows us to determine which combinations of parameters are contributing the most to the uncertainty.

At the minimum of $`\chi _{\mathrm{global}}^2`$, the first derivatives with respect to the $`\{a_i\}`$ are zero; so near the minimum, $`\chi _{\mathrm{global}}^2`$ can be approximated by

$$\chi _{\mathrm{global}}^2=\chi _0^2+\frac{1}{2}\underset{i,j}{}F_{ij}y_iy_j$$
(2)
where $`y_i=a_ia_{0i}`$ is the displacement from the minimum, and $`F_{ij}`$ is the *Hessian*, the matrix of second derivatives. It is natural to define a new set of coordinates using the complete orthonormal set of eigenvectors of the symmetric matrix $`F_{ij}`$ as basis vectors. These vectors can be ordered by their eigenvalues $`e_i`$. Each eigenvalue is a quantitative measure of the uncertainties in the shape parameters $`\{a_i\}`$ for displacements in parameter space in the direction of the corresponding eigenvector. The quantity $`\mathrm{}_i1/\sqrt{e_i}`$ is the distance in the 18 dimensional parameter space, in the direction of eigenvector $`i`$, that makes a unit increase in $`\chi _{\mathrm{global}}^2`$. If the only measurement uncertainty were uncorrelated gaussian uncertainties, then $`\mathrm{}_i`$ would be one standard deviation from the best fit in the direction of the eigenvector. The inverse of the Hessian is the error matrix.

Because the real uncertainties, for the wide variety of experiments included, are far more complicated than assumed in the ideal situation, the quantitative measure of a given increase in $`\chi _{global}^2`$ carries

little true statistical meaning. However, qualitatively, the Hessian gives an analytic picture of $`\chi _{\mathrm{global}}^2`$ near its minimum in $`\{a_i\}`$ space, and hence allows us to identify the particular degrees of freedom that need further experimental input in future global analyses.

From calculations of the Hessian we find that the eigenvalues vary over a wide range. Figure 3 shows a graph of the eigenvalues of $`F_{ij}`$, on a logarithmic scale. The vertical axis is $`\mathrm{}_i=1/\sqrt{e_i}`$, the distance of a “standard deviation” along the $`i^{\mathrm{th}}`$ eigenvector. These distances range over 3 orders of magnitude.Large eigenvalues of $`F_{ij}`$ correspond to “steep directions” of $`\chi _{\mathrm{global}}^2`$. The corresponding eigenvectors are combinations of shape parameters that are well determined by current data. For example, parameters that govern the valence $`u`$ and $`d`$ quarks at moderate $`x`$ are sharply constrained by DIS data. Small eigenvalues of $`F_{ij}`$ correspond to “flat directions” of $`\chi _{\mathrm{global}}^2`$. In the directions of these eigenvectors, $`\chi _{\mathrm{global}}^2`$ changes little over large distances in $`\{a_i\}`$ space. For example, parameters that govern the large-$`x`$ behavior of the gluon distribution, or differences between sea quarks, properties of the nucleon that are not accurately determined by current data, contribute to the flat directions. The existence of flat directions is inevitable in global fitting, because as the data improve it only makes sense to maintain enough flexibility for $`f_a(x,Q_0)`$ to fit the available experimental constraints.

Because the eigenvalues of the Hessian have a large range of values, efficient calculation of $`F_{ij}`$ requires an adaptive algorithm. In principle $`F_{ij}`$ is the matrix of second derivatives at the minimum of $`\chi _{\mathrm{global}}^2`$, which could be calculated from very small finite differences. In practice, small computational errors in the evaluation of $`\chi _{\mathrm{global}}^2`$ preclude the use of a very small step size. Coarse grained finite differences yield a more accurate calculation of the second derivatives. But because the variation of $`\chi _{\mathrm{global}}^2`$ varies markedly in different directions, it is important to use a grid in $`\{a_i\}`$ space with small steps in steep directions and large steps in flat directions. This grid is generated by an iterative procedure, in which $`F_{ij}`$ converges to a good estimate of the second derivatives.

From calculations of $`F_{ij}`$ we find that the minimum of $`\chi _{\mathrm{global}}^2`$ is fairly quadratic over large distances in the parameter space. Figures 4 and 5

show the behavior of $`\chi _{\mathrm{global}}^2`$ near the minimum along each of the 18 eigenvectors. $`\chi _{\mathrm{global}}^2`$ is plotted on the vertical axis, and the variable on the horizontal axis is the distance in $`\{a_i\}`$ space in the direction of the eigenvector, in units of $`\mathrm{}_i=1/\sqrt{e_i}`$. There is some nonlinearity, but it is small enough that the Hessian can be used as an analytic model of the functional dependence of $`\chi _{\mathrm{global}}^2`$ on the shape parameters.

In a future paper we will provide details on the uncertainties of the original shape parameters $`\{a_i\}`$. But it should be remembered that these parameters specify the PDFs at the low $`Q`$ scale, and applications of PDFs to Tevatron experiments use PDFs at a high $`Q`$ scale. The evolution equations determine $`f(x,Q)`$ from $`f(x,Q_0)`$, so the functional form at $`Q`$ depends on the $`\{a_i\}`$ in a complicated way.

## 4 Uncertainty on $`\sigma _W`$: the Lagrange Multiplier Method

In this Section, we determine the variation of $`\chi _{\mathrm{global}}^2`$ as a function of a single measurable quantity. We use the production cross section for $`W`$ bosons ($`\sigma _W`$) as an archetype example. The same method can be applied to any other physical observable of interest, for instance the Higgs production cross section, or to certain measured differential distributions. The aim is to quantify the uncertainty on that physical observable due to uncertainties of the PDFs integrated over the entire PDF parameter space.

Again, we use the standard CTEQ5 analysis tools and results  as the starting point. The “best fit” is the CTEQ5M1 set. A natural way to find the limits of a physical quantity $`X`$, such as $`\sigma _W`$ at $`\sqrt{s}=1.8`$ TeV, is to take $`X`$ as one of the search parameters in the global fit and study the dependence of $`\chi _{\mathrm{global}}^2`$ for the 15 base experimental data sets on $`X`$.

Conceptually, we can think of the function $`\chi _{\mathrm{global}}^2`$ that is minimized in the fit as a function of $`a_1,\mathrm{},a_{17},X`$ instead of $`a_1,\mathrm{},a_{18}`$. This idea could be implemented directly in principle, but a more convenient way to do the same thing in practice is through Lagrange’s method of undetermined multipliers. One minimizes, with respect to the $`\{a_i\}`$, the quantity

$$F(\lambda )=\chi _{\mathrm{global}}^2+\lambda X(a_1,\mathrm{},a_{18})$$
(3)
for a fixed value of $`\lambda `$, the Lagrange multiplier. By minimizing $`F(\lambda )`$ for many values of $`\lambda `$, we map out $`\chi _{\mathrm{global}}^2`$ as a function of $`X`$. The minimum of $`F`$ for a given value of $`\lambda `$ is the best fit to the data for the corresponding value of $`X`$, i.e., evaluated at the minimum.

Figure 6 shows $`\chi _{\mathrm{global}}^2`$ for the 15 base experimental data sets as a function of $`\sigma _W`$ at the Tevatron. The horizontal axis is $`\sigma _W`$ times the branching ratio for $`W`$ leptons, in nb. The CTEQ5m prediction is $`\sigma _WBR_{\mathrm{lep}}=2.374`$ nb.

The vertical dashed lines are $`\pm 3`$% and $`\pm 5`$% deviations from the CTEQ5m prediction.

The two parabolas associated with points in Fig. 6 correspond to different treatments of the normalization factor $`N_n`$ in Eq. 1. The dots are variable norm fits, in which $`N_n`$ is allowed to float, taking into account the experimental normalization uncertainties, and $`F(\lambda )`$ is minimized with respect to $`N_n`$. The justification for this procedure is that overall normalization is a common systematic uncertainty. The square points are fixed norm fits, in which all $`N_n`$ are held fixed at their values for the global minimum (CTEQ5m). These two procedures represent extremes in the treatment of normalization uncertainty. The parabolas are just least-square fits to the points in the two cases.

The other curve in Fig. 6 was calculated using the Hessian method. The Hessian $`F_{ij}`$ is the matrix of second derivatives of $`\chi _{\mathrm{global}}^2`$ with respect to the shape parameters $`\{a_i\}`$. The derivatives (first and second) of $`\sigma _W`$ may also be calculated by finite differences. Using the resultant quadratic approximations for $`\chi _{\mathrm{global}}^2(a)`$ and $`\sigma _W(a)`$, one may minimize $`\chi _{\mathrm{global}}^2`$ with $`\sigma _W`$ fixed. Since this calculation keeps the normalization factors constant, it should be compared with the fixed norm fits from the Lagrange multiplier method. The fact that the Hessian and Lagrange multiplier methods yield similar results lends support to both approaches; the small difference between them indicates that the quadratic functional approximations for $`\chi _{\mathrm{global}}^2`$ and $`\sigma _W`$ are only approximations.

For the quantitative analysis of uncertainties, the important question is: How large an increase in $`\chi _{\mathrm{global}}^2`$ should be taken to define the likely range of uncertainty in $`X`$? There is an elementary statistical theorem that states that $`\mathrm{\Delta }\chi ^2=1`$ in a constrained fit corresponds to 1 standard deviation of the constrained quantity $`X`$. However, the theorem relies on the assumption that the uncertainties are gaussian, uncorrelated, and correctly estimated in magnitude. Because these conditions do not hold for the full data set (of $``$ 1300 points from 15 different experiments), this theorem cannot be naively applied quantitatively.444It has been shown by Giele et.al. , that, taken literally, only one or two selected experiments satisfy the standard statistical tests. Indeed, it can be shown that, if the measurement uncertainties are correlated, and the correlation is not properly taken into account in the definition of $`\chi _{\mathrm{global}}^2`$, then a standard deviation may vary over the entire range from $`\mathrm{\Delta }\chi ^2=1`$ to $`\mathrm{\Delta }\chi ^2=N`$ (the total number of data points – $``$ 1300 in our case).

## 5 STATISTICAL ANALYSIS WITH SYSTEMATIC UNCERTAINTIES

Fig. 6 shows how the fitting function $`\chi _{\mathrm{global}}^2`$ increases from its minimum value, at the best global fit, as the cross-section $`\sigma _W`$ for $`W`$ production is forced away from the prediction of the global fit. The next step in our analysis of PDF uncertainty is to use that information, or some other analysis, to estimate the uncertainty in $`\sigma _W`$. In ideal circumstances we could say that a certain increase of $`\chi _{\mathrm{global}}^2`$ from the minimum value, call it $`\mathrm{\Delta }\chi ^2`$, would correspond to a standard deviation of the global measurement uncertainty. Then a horizontal line on Fig. 6 at $`\chi _{\mathrm{min}}^2+\mathrm{\Delta }\chi ^2`$ would indicate the probable range of $`\sigma _W`$, by the intersection with the parabola of $`\chi _{\mathrm{global}}^2`$ versus $`\sigma _W`$.

However, such a simple estimate of the uncertainty of $`\sigma _W`$ is not possible, because the fitting function $`\chi _{\mathrm{global}}^2`$ does not include the correlations between systematic uncertainties. The uncertainty $`\sigma _{ni}^d`$ in the definition (1) of $`\chi _{\mathrm{global}}^2`$ combines in quadrature the statistical and systematic uncertainties for each data point; that is, it treats the systematic uncertainties as uncorrelated. The standard theorems of statistics for Gaussian probability distributions of random uncertainties do not apply to $`\chi _{\mathrm{global}}^2`$.

Instead of using $`\chi _{\mathrm{global}}^2`$ to estimate confidence levels on $`\sigma _W`$, we believe the best approach is to carry out a thorough statistical analysis, including the correlations of systematic uncertainties, on individual experiments used in the global fit for which detailed information is available. We will describe here such an analysis for the measurements of $`F_2(x,Q)`$ by the H1 experiment  at HERA, as a case study. In a future paper, we will present similar calculations for other experiments.

The H1 experiment has provided a detailed table of measurement uncertainties – statistical and systematic – for their measurements of $`F_2(x,Q)`$.  The CTEQ program uses 172 data points from H1 (requiring the cut $`Q^2>5`$ GeV2). For each measurement $`d_j`$ (where $`j=1\mathrm{}172`$) there is a statistical uncertainty $`\sigma _{0j}`$, an uncorrelated systematic uncertainty $`\sigma _{1j}`$, and a set of 4 correlated systematic uncertainties $`a_{jk}`$ where $`k=1\mathrm{}4`$. (In fact there are 8 correlated uncertainties listed in the H1 table. These correspond to 4 pairs. Each pair consists of one standard deviation in the positive sense, and one standard deviation in the negative sense, of some experimental parameter. For this first analysis, we have approximated each pair of uncertainties by a single, symmetric combination, equal in magnitude to the average magnitude of the pair.)

To judge the uncertainty of $`\sigma _W`$, as constrained by the H1 data, we will compare the H1 data to the global fits in Fig. 6. The comparison is based on the true, statistical $`\chi ^2`$, including the correlated uncertainties, which is given by

$$\chi ^2=\underset{j}{}\frac{\left(d_jt_j\right)^2}{\sigma _j^2}\underset{kk^{}}{}B_k\left(A^1\right)_{kk^{}}B_k^{}.$$
(4)
The index $`j`$ labels the data points and runs from 1 to 172. The indices $`k`$ and $`k^{}`$ label the source of systematic uncertainty and run from 1 to 4. The combined uncorrelated uncertainty $`\sigma _j`$ is $`\sqrt{\sigma _{0j}^2+\sigma _{1j}^2}`$. The second term in (4) comes from the correlated uncertainties. $`B_k`$ is the vector

$$B_k=\underset{j}{}\frac{\left(d_jt_j\right)a_{jk}}{\sigma _j^2},$$
(5)
and $`A_{kk^{}}`$ is the matrix

$$A_{kk^{}}=\delta _{kk^{}}+\underset{j}{}\frac{a_{jk}a_{jk^{}}}{\sigma _j^2}.$$
(6)
Assuming the published uncertainties $`\sigma _{0j}`$, $`\sigma _{1j}`$ and $`a_{jk}`$ accurately reflect the measurement fluctuations, $`\chi ^2`$ would obey a chi-square distribution if the measurements were repeated many times. Therefore the chi-square distribution with 172 degrees of freedom provides a basis for calculating confidence levels for the global fits in Fig. 6.

Table 2 shows $`\chi ^2`$ for the H1 data compared to seven of the PDF fits in Fig. 6. The center row of the Table is the global best fit – CTEQ5m. The other rows are fits obtained by the Lagrange multiplier method for different values of the Lagrange multiplier. The best fit to the H1 data, i.e., the smallest $`\chi ^2`$, is not CTEQ5m (the best global fit) but rather the fit with Lagrange multiplier 1000 for which $`\sigma _W`$ is 0.8% smaller than the prediction of CTEQ5m. Forcing the $`W`$ cross section values away from the prediction of CTEQ5m causes an increase in $`\chi ^2`$ for the DIS data. At $`\sqrt{s}=1.8`$ TeV, $`W`$ production is mainly from $`q\overline{q}W^+W^{}`$ with moderate values of $`x`$ for $`q`$ and $`\overline{q}`$, i.e., values in the range of DIS experiments. Forcing $`\sigma _W`$ higher (or lower) requires a higher (or lower) valence quark density in the proton, in conflict with the DIS data, so $`\chi ^2`$ increases.

The final column in Table 2, labeled “probability”, is computed from the chi-square distribution with 172 degrees of freedom. This quantity is the probability for $`\chi ^2`$ to be greater than the value calculated from the existing data, if the H1 measurements were to be repeated. So, for example, the fit with Lagrange multiplier $`3000`$, which corresponds to $`\sigma _W`$ being 3.2% larger than the CTEQ5m prediction, has probability 0.092. In other words, if the H1 measurements could be repeated many times, in only 9.2% of trials would $`\chi ^2`$ be greater than or equal to the value that has been obtained with the existing data. This probability represents a confidence level for the value of $`\sigma _W`$ that was forced on the PDF by setting the Lagrange multiplier equal to -3000. At the 9.2% confidence level we can say that $`\sigma _WBR_{\mathrm{lep}}`$ is less than 2.450 nb, based on the H1 data. Similarly, at the 21.2% confidence level we can say that $`\sigma _WBR_{\mathrm{lep}}`$ is greater than 2.294 nb.

Fig. 7 is a graph of $`\chi ^2/N`$ for the H1 data compared to the PDF fits in Table 2. This figure may be compared to Fig. 6. The CTEQ5 prediction of the $`W`$ production cross-section is shown as an arrow, and the vertical dashed lines are $`\pm 3`$% away from the CTEQ5m prediction. The horizontal dashed line is the 68% confidence level on $`\chi ^2/N`$ for $`N=172`$ degrees of freedom. The comparison with H1 data alone indicates that the uncertainty on $`\sigma _W`$ is $`3`$%.

There is much more to say about $`\chi ^2`$ and confidence levels. In a future paper we will discuss statistical calculations for other experiments in the global data set. The H1 experiment is a good case, because for H1 we have detailed information about the correlated uncertainties. But it may be somewhat fortuitous that the $`\chi ^2`$ per data point for CTEQ5m is so close to 1 for the H1 data set. In cases where $`\chi ^2/N`$ is not close to 1, which can easily happen if the estimated systematic uncertainties are not textbook-like, we must supply further arguments about confidence levels. For experiments with many data points, like 172 for H1, the chi-square distribution is very narrow, so a small inaccuracy in the estimate of $`\sigma _j`$ may translate to a large uncertainty in the calculation of confidence levels based on the absolute value of $`\chi ^2`$. Because the estimation of experimental uncertainties introduces some uncertainty in the value of $`\chi ^2`$, it is not really the absolute value of $`\chi ^2`$ that is important, but rather the relative value compared to the value at the global minimum. Therefore, we might study ratios of $`\chi ^2`$’s to interpret the variation of $`\chi ^2`$ with $`\sigma _W`$.

## 6 CONCLUSIONS

It has been widely recognized by the HEP community, and it has been emphasized at this workshop, that PDF phenomenology must progress from the past practice of periodic updating of *representative* PDF sets to a systematic effort to map out the uncertainties, both on the PDFs themselves and on physical observables derived from them. For the analysis of PDF uncertainties, we have only addressed the issues related to the treatment of experimental uncertainties. Equally important for the ultimate goal, one must come to grips with uncertainties associated with theoretical approximations and phenomenological parametrizations. Both of these sources of uncertainties induce highly correlated uncertainties, and they can be numerically more important than experimental uncertainties in some cases. Only a balanced approach is likely to produce truly useful results. Thus, great deal of work lies ahead.

This report described first results from two methods for quantifying the uncertainty of parton distribution functions associated with experimental uncertainties. The specific work is carried out as extensions of the CTEQ5 global analysis. The same methods can be applied using other parton distributions as the starting point, or using a different parametrization of the non-perturbative PDFs. We have indeed tried a variety of such alternatives. The results are all similar to those presented above. The robustness of these results lends confidence to the general conclusions.

The Hessian, or error matrix method reveals the uncertainties of the shape parameters used in the functional parametrization. The behavior of $`\chi _{\mathrm{global}}^2`$ in the neighborhood of the minimum is well described by the Hessian if the minimum is quadratic.

The Lagrange multiplier method produces constrained fits, i.e., the best fits to the global data set for specified values of some observable. The increase of $`\chi _{\mathrm{global}}^2`$, as the observable is forced away from the predicted value, indicates how well the current data on PDFs determines the observable.

The constrained fits generated by the Lagrange multiplier method may be compared to data from individual experiments, taking into account the uncertainties in the data, to estimate confidence levels for the constrained variable. For example, we estimate that the uncertainty of $`\sigma _W`$ attributable to PDFs is $`\pm 3`$%.

Further work is needed to apply these methods to other measurements, such as the $`W`$ mass or the forward-backward asymmetry of $`W`$ production in $`p\overline{p}`$ collisions. Such work will be important for analyzing high precision measurements involving the electroweak bosons in the Tevatron Run II era.

## Acknowledgement

This research is supported in part by the National Science Foundation.


format: text

text:
# Anomalous Coherent Backscattering of Light from Opal Photonic Crystals

## I Introduction

Expreimental studies of coherent backscattering (CBS) of light from disordered media have been completed by many groups during the last decade.1-8 CBS measurements are of great interest because they have allowed to test the physics of weak localization, which was first uncovered in electronic transport, when applied to light waves and also to determine the light mean-free-path, $`l^{}`$ in random media. In CBS measurements the normalized reflected light intensity (also called the albedo, $`\alpha `$) is monitored versus the deviation angle, $`\delta \theta `$ from the backscattering direction of the incident laser beam. In the multiple scattering regime interference effects are averaged out except in a small angular range around $`\delta \theta =0`$, where constructive interference originating from reciprocal light paths enhance the diffusive reflection from the sample. As a result a coherent scattering intensity cone is formed in the albedo, $`\alpha (\delta \theta )`$ for small $`\delta \theta `$. From the full width, $`\mathrm{\Delta }\theta `$ of the CBS cone at quarter maximum, $`l^{}`$ is obtained using the relation $`\mathrm{\Delta }\theta =\lambda _L/2\pi l^{}`$, where $`\lambda _L`$ is the impinging laser wavelengh.

The first CBS measurements were performed on suspensions of polystyrene spheres. Since the suspension particles were very small in diameter, their Brownian motion provided a natural ensemble averaging. For solid materials, however CBS usually shows a noisy but repeatable mesoscopic pattern. Under these conditions the CBS cone around $`\delta \theta =0`$ can be obtained by applying an elaborate ensemble-averaging technique in order to suppress the dominant mesoscopic fluctuations. The averaging techniques employed to get the cone include spot variation, tilt angle variation, and sample spin-rotation. All the above experimental studies1-8 were performed on completely disordered, isotropic and homogeneous media.

Photonic crystals have recently attracted a lot of attention because of potential optical device applications. These applications pertain to photonic crystals with complete band gaps i.e. with no propagating solutions of the Maxwell equations within a certain frequency range. However, most of the available three dimensional photonic crystals nowadays have incomplete gaps in their photonic band structure, with light propagation forbidden only along certain directions. Along these directions and for frequencies within the gap the light intensity decays with a certain decrement $`L_B`$-that is the Bragg attenuation length. In the ideal case (when disorder is negligibly small) the length $`L_B`$ can be extracted from spectroscopic reflectivity measurements. Since $`L_B`$ is directly related to the gap width, then its value can be deduced from the angular (or spectral) width of the Bragg peak in reflectivity. However, in the more realistic situation disorder is present in the sample and this broadens the Bragg band. We also note that even if $`L_Bl^{}`$ broadening can be substantial due to the sensitivity of the Bragg diffraction to the properties of the surface boundary. Under these conditions a separate determination of $`L_B`$ and $`l^{}`$ is a challenge.

In general, despite the conceptual elegance of the CBS technique, the main quantitative information it has provided2-8 so far is $`l^{}`$ value in macroscopically homogeneous and isotropic disordered media. We note that $`l^{}`$ value in such media can be also determined from more conventional measurements of the transmission versus thickness that do not rely on the interference of clockwise and anticlockwise scattering paths.

In the present work we show that for optically periodic structures with disorder the CBS measurements provide a unique opportunity to separate the two effects, namely Bragg attenuation and light diffusion lengths. This reveals the tremendous potential of the CBS technique as compared to the more conventional optical spectroscopies.

## II Experimental Results

For our studies of periodic structures we have chosen the synthetic $`SiO_2`$ opals that are actually three-dimensional photonic crystals. The synthetic opal single crystals were fabricated from a colloidal suspension of mono-dispersed silica ($`SiO_2`$) spheres with a mean diameter, $`𝒟`$=$`295nm`$ and dispersion in $`𝒟`$ of about $`4\%`$. In contrast to the natural gem opals that exist in the market the space between $`SiO_2`$ spheres in artificial opals is filled with air. A sintering process leads to the formation of inter-sphere necks that mechanically connect neighboring silica spheres in a face-centered-cubic (fcc) lattice. The opals with air-filled voids are chalk-like diffuse light scatterers that show some iridescence due to Bragg diffraction from their crystal planes. The iridescence is dramatically enhanced, however if the voids are filled with an index matching liquid.

In opals there exist photonic gaps, $`\mathrm{\Delta }\nu _{hkl}`$, in the frequency spectrum around each $`\nu _{hkl}`$, where $`(hkl)`$ are the Miller indices of the crystal planes. The central frequencies, $`\nu _{hkl}`$, are determined by the Bragg condition: $`\lambda _B=c\nu _{hkl}^1=2n_{eff}d_{hkl}`$, where $`d_{hkl}`$ is the inter-planar spacing, $`n_{eff}`$ is the opal effective index of refraction ($`n_{eff}=1.376`$ for air-filled opals) and $`c`$ is the speed of light. In the absence of disorder the stop band width, $`\mathrm{\Delta }\nu _B`$ is related to the Bragg attenuation length, $`L_B`$, by the relation $`\mathrm{\Delta }\nu _B/\nu _B=2d/\pi L_B`$ (with implied $`hkl`$ indices), where $`\nu _B=c/\lambda _B`$. As mentioned above, however disorder in opals dramatically broadens $`\mathrm{\Delta }\nu _B`$ so that this relation is no longer valid.

For our CBS measurements, we cut a centimeter-size opal single crystal into planes orthogonal to several $`\left[hkl\right]`$ directions including the fast growing $`\left[111\right]`$ direction and a prominent $`\left[220\right]`$ direction. We note that the $`\left(200\right)`$ planes in opals are rather disordered.

We used a standard experimental set-up with angular resolution in $`\delta \theta `$ of less than $`1mrad`$. The coherent light beam was directed from continuous wave lasers such as $`Ar^+`$ and $`HeNe`$, as well as various well-collimated semiconductor lasers. The laser beam was linearly polarized with transverse electric field ($`TE`$) polarization. The reflected intensity from the opal $`\left(hkl\right)`$ surface with $`TE`$ polarization was measured versus $`\delta \theta `$ using a silicon detector and phase-sensitive techniques. The incident beam ($`5mm`$ in diameter) was directed at various angles, $`\theta `$ between $`30^{}`$ to $`70^{}`$ relative to the surface normal, and $`\delta \theta `$ was varied within $`200mrad`$. To obtain specle-free CBS cones we performed configuration averaging using three different methods. These include different illuminated areas of the opal surface, different inclination angles $`\phi <2^{}`$ of the illuminated surface with respect to the laser beam, and sample rotation about the normal direction to the crystal surface.

In Fig. 1 we show typical CBS cones measured at various incident angles $`\theta `$ with respect to $`\left[111\right]`$ that we obtained using an $`Ar^+`$ laser at $`\lambda _L=515nm`$. The CBS cones are apparent in all cases and their width, $`\mathrm{\Delta }\theta `$ at quarter maximum versus $`\theta `$ is summarized in Fig. 2. Except for a small $`\theta `$ interval, $`\gamma `$ (about $`5^{}`$) around $`\theta =42^{}`$, $`\mathrm{\Delta }\theta `$ monotonously increases with $`\theta `$. Within the interval $`\gamma `$ around $`\theta =42^{}`$, however there is an anomalous increase in $`\mathrm{\Delta }\theta `$ of about a factor of two (Fig. 2).

Insight into the origin of the observed anomaly in $`\mathrm{\Delta }\theta `$ ($`\theta `$) may be obtained from Fig. 3 in which the angle-dependent reflectivity spectra, $`R(\lambda )`$, from the $`\left(111\right)`$ opal surface are shown. It is seen that for $`\lambda _L=515nm`$ the angle $`\theta =\theta _B=42^{}`$, at which the anomaly in the CBS cones occurs, in fact corresponds to the Bragg conditions.

Another evidence that the CBS anomaly is related to the Bragg conditions can be observed from the obtained CBS cones at $`\theta `$ around $`\theta _B=0`$ (Fig. 4). These CBS cones were measured using a $`HeNe`$ laser at $`\lambda _L=633nm`$ that is close to $`\lambda _B`$ at $`\theta =0`$ (Fig. 3, inset). Fig. 4 inset shows that $`\mathrm{\Delta }\theta `$ for CBS cones measured at normal incidence stays approximately constant at about $`17mrad`$ for laser wavelengths outside the Bragg conditions. However, for CBS measured at $`\theta \theta _B=0`$, $`\mathrm{\Delta }\theta `$ in the CBS cone increases by a factor of $`4`$. Also the incident angle interval ($`\gamma 20^{}`$, Fig. 4) in which the CBS anomalies are observed with laser illumination at $`\lambda _L=633nm`$ is much larger than that obtained at $`\lambda _L=515nm`$ ($`\gamma 5^{}`$, Fig. 2). Both observations, namely the anomaly in $`\mathrm{\Delta }\theta `$ ($`\theta `$) at $`\theta =\theta _B`$ and strengthening of the anomaly at $`\theta _B0`$ require an explanation. A possible model is given in the next Section.

## III Discussion

The Bragg-related anomalies of the CBS cones show that light transport and localization is substantially affected by the $`\left[111\right]`$ photonic band gap in the opal, where light becomes more localized. To understand these results we first recall the most transparent derivation of the CBS cones in the literature and apply it to laser frequencies outside the photonic gap. We then modify the theory to interpret the anomalous CBS results at intra-gap laser frequencies at $`\theta \theta _B`$.

The coherent scattering albedo, $`\alpha _c(\theta ,\theta ^{})`$, at an angle $`\theta ^{}`$ from a semi-infinite slab at $`z=0`$ due to an incident laser beam at an angle $`\theta `$ is determined by the first and last light scattering events occurring at the spatial points $`𝐫`$ and $`𝐫^{}`$, respectively, and by the light diffusion process that occurs between these events. If $`z`$ and $`z^{}`$ are the projections of $`𝐫`$ and $`𝐫^{}`$ on the $`z`$-axis, and $`\rho `$ is the projection of $`(𝐫𝐫^{})`$ on the sample interface plane, then the coherent part of albedo is given by the expression

$$\alpha _c(\theta ,\theta ^{})_0^{\mathrm{}}𝑑z𝑑z^{}d^2\rho 𝒫_\varphi (z)\left\{\mathrm{cos}\left[k_{}\rho +k(\mathrm{cos}\varphi \mathrm{cos}\varphi ^{})(zz^{})\right]\right\}\mathrm{\Pi }(\rho ,z,z^{})𝒫_\varphi ^{}(z^{}).$$
(1)
Here $`k=2\pi n_{eff}/\lambda _L`$ and $`k_{}=k(\mathrm{sin}\varphi \mathrm{sin}\varphi ^{})`$, where the angles $`\varphi `$ and $`\varphi ^{}`$ are, respectively, the propagation directions of the incoming and outgoing light inside the slab. These angles are related to $`\theta `$ and $`\theta ^{}`$ through the Snell law $`n_{eff}\mathrm{sin}\varphi =\mathrm{sin}\theta `$. In Eq. (1) $`𝒫_\varphi (z)`$ is the probability density for the incident light to propagate up to the depth $`z`$ before the occurence of the first scattering event. Analogously, $`𝒫_\varphi ^{}(z^{})`$ is the surviving probability density after the occurence of the last scattering event. $`𝒫_\varphi ^{}(z^{})`$ has the same form as $`𝒫_\varphi (z)`$, where $`\varphi `$, $`z`$ are replaced by $`\varphi ^{}`$, $`z^{}`$ respectively. Since $`𝒫_\varphi (z)`$ satisfies the Poisson statistics, it is given by $`𝒫_\varphi (z)=(l^{}\mathrm{cos}\varphi )^1\mathrm{exp}\left(z/l^{}cos\varphi \right)`$. The propagator $`\mathrm{\Pi }(\rho ,z,z^{})`$ in Eq. (1) describes the diffusive motion of light between the first and last scattering events. To take into account the fact that diffusion occurs in the presence of the boundary $`z=0`$, it is customary to choose the propagator $`\mathrm{\Pi }`$ in the form

$$\mathrm{\Pi }(\rho ,z,z^{})=\mathrm{\Pi }_0(\rho ,z,z^{})\mathrm{\Pi }_0(\rho ,z,z^{}),$$
(2)
where $`\mathrm{\Pi }_0(\rho ,z,z^{})=\left[4\pi D\left(\left(zz^{}\right)^2+\rho ^2\right)^{1/2}\right]^1`$ is the diffusion propagator in an unbounded medium; $`D=l^{}c/3`$ is the light diffusion coefficient, and $`z^{}`$ is the mirror image of $`z^{}`$ with respect to the plane $`z=z_0`$, where $`z_0=0.7l^{}`$. The choice (2) of the propagator $`\mathrm{\Pi }`$ insures that the boundary conditions $`\mathrm{\Pi }(\rho ,z_0,z^{})=\mathrm{\Pi }(\rho ,z,z_0)=0`$ are obeyed. The evaluation of integrals in Eq. (1) for small $`\delta \theta =\theta \theta ^{}`$ yields

$$\alpha _c\left(\delta \theta \right)\left(\frac{\mathrm{cos}\varphi }{1+k_{}l^{}\mathrm{cos}\varphi }\right)^2\left(\frac{1e^{2k_{}z_0}}{k_{}l^{}}+\mathrm{cos}\varphi \right),$$
(3)
where $`\delta \theta `$ enters into the r.h.s. of Eq. (3) through $`k_{}=2\pi \mathrm{cos}\theta \delta \theta /\lambda _L`$ and $`\theta `$ and $`\varphi `$ are related through Snell’s law. Eq. (3) contains only a single free parameter, namely $`l^{}`$.

We used Eq. (3) to fit the CBS cones in $`\alpha _c(\delta \theta )`$ at various $`\theta `$ outside the anomaly angle interval $`\gamma `$ around $`\theta =\theta _B`$. The excellent fits shown in Fig. 1 for $`\theta =10^{},34^{}`$ and $`60^{}`$ were obtained with a single $`l^{}=7.2\mu m`$. From this value for the optical mean-free-path length and Eq. (3), we calculated the width of the albedo for the CBS cone at quarter maximum, $`\mathrm{\Delta }\theta `$, versus $`\theta `$, as shown in Fig. 2. Except near $`\theta _B`$ the overall agreement between the calculated and the experimentally determined $`\mathrm{\Delta }\theta `$ is excellent showing that light transport is quite ordinary for frequencies outside the photonic gap.

For describing the CBS process at $`\theta \theta _B`$ and correspondingly at $`\varphi \varphi _B`$ we note that the diffusive propagation of light between the first and last scattering events that is given by $`\mathrm{\Pi }(\rho ,z,z^{})`$ in Eq. (1) is affected very little by the Bragg diffraction. The reason for this is the small angle interval, $`\beta (\theta _B)`$, around $`\theta _B`$ where Bragg condition is satisfied. An upper bound estimate for $`\beta (\theta _B)`$ can be obtained by relating it to the spectral width of the Bragg gap $`\delta \lambda _B`$, where $`\delta \lambda _B/\lambda _L=\beta (\theta _B)\mathrm{sin}2\theta _B/2\left(n_{eff}^2\mathrm{sin}^2\theta _B\right)`$. If we choose for $`\delta \lambda _B/\lambda _L`$ the measured value of the reflectivity peak $`\delta \lambda _B/\lambda _L=0.09`$, then for $`\lambda _L=515nm`$ and $`\theta _B42^{}`$ we calculate $`\beta `$ to be $`0.25`$ radian. This forms a stereo-angle of about $`0.07`$ stereo-radian around the Bragg direction that is only about $`0.5\%`$ of the complete $`4\pi `$ stereo-angle. Therefore, even if the light beam is incident at $`\theta =\theta _B`$, it is diverted after the first scattering event and is then followed by an almost regular diffusion; this is true since the probability for light to propagate close to $`\theta _B`$ again is very small. This permits the generalization of the CBS theory to photonic crystals with relatively small stop bands, simply by only modifying the probability densities $`𝒫_\varphi (z)`$ and $`𝒫_\varphi ^{}(z^{})`$ in Eq. (1). For $`\theta \theta _B`$, light intensity decays as $`\mathrm{exp}\left(z/L_B\right)`$ due to the Bragg diffraction from the photonic crystal planes. Therefore $`𝒫_\varphi (z)`$ can be written at $`\theta =\theta _B`$ as a product of the surviving probabilities: $`\mathrm{exp}\left(z/l^{}cos\varphi \right)`$ and the Bragg attenuation, $`\mathrm{exp}\left(z/L_B\right)`$. Equivalently, we may represent $`𝒫_\varphi (z)`$ at $`\theta =\theta _B`$ as a single exponent, $`\mathrm{exp}\left(z/l_{eff}^{}\mathrm{cos}\varphi _B\right)`$, where $`\left(l_{eff}^{}\left(\varphi _B\right)\right)^1=\left(l^{}\right)^1+\mathrm{cos}\varphi _B\left(L_B\right)^1`$. This expression shows that Eq. (3) describing the CBS cone is still valid at $`\theta =\theta _B`$, but with a smaller effective mean-free-path, $`l_{eff}^{}`$. A smaller $`l_{eff}^{}`$ leads to a broader CBS cone at $`\theta =\theta _B`$, as indeed seen in Fig. 2 for $`\theta _B42^{}`$.

If $`\theta `$ and $`\theta ^{}`$ are close, but not equal to $`\theta _B`$, then the Bragg attenuation becomes weaker. It can be shown that this effect is taken into account by the following modification of $`l_{eff}^{}(\varphi )`$

$$\frac{1}{l_{eff}^{}\left(\varphi \right)}=\frac{1}{l^{}}+\frac{\mathrm{cos}\varphi _B}{L_B}\left[1\frac{\pi ^2L_B^2\mathrm{sin}^2\left(2\theta _B\right)}{\lambda _L^2\left(n_{eff}^2\mathrm{sin}^2\theta _B\right)}\left(\theta \theta _B\right)^2\right]^{1/2}.$$
(4)
Since the argument under the square root in Eq. (4) must be positive, we use this constraint to estimate the angle interval $`\gamma `$ around $`\theta _B`$ where the CBS anomalies due to the photonic crystal are observed. This leads to the expression: $`\gamma (\theta _B)=\lambda _L^2/\pi L_Bd\mathrm{sin}\left(2\theta _B\right)`$. Fig. 2 shows that an excellent fit is obtained for both $`\mathrm{\Delta }\theta `$ enhancement and $`\gamma `$, using Eqs. (3) and (4) with $`L_B=5.1\mu m`$ and the previously obtained $`l^{}=7.2\mu m`$. Using these values of $`l^{}`$ and $`L_B`$ along $`\left[111\right]`$ we evaluated $`\gamma (42^{})=4^{}`$, which agrees well with the data of Fig. 2. We also calculated the “pure” Bragg band $`\delta \lambda _B`$ to be $`\delta \lambda _B/\lambda _L0.03`$. We note, that this value is 3 times smaller than $`\delta \lambda /\lambda _L0.09`$ extracted from the Bragg reflectivity band in Fig. 3. This illustrates the ambiguity in determining the photonic band structure form the more conventional reflectivity spectroscopy.

For $`\theta 0`$ i.e. close to normal incidence and $`\lambda _L\lambda _B`$, Eq. (4) is no longer valid, since $`|\theta \theta _B|\theta _B`$. In this case it can be shown that the expression for the effective light mean-free-path is given by

$$\frac{1}{l_{eff}^{}\left(\varphi \right)}=\frac{1}{l^{}}+\frac{1}{L_B}\left[1\left(\frac{\pi L_B\theta ^2}{\lambda _Ln_{eff}}\right)^2\right]^{1/2}.$$
(5)
Again the constraint that the argument under the square root is positive now yields $`\gamma (0^{})=26^{}`$. This value is much larger than the obtained $`\gamma (42^{})=4^{}`$ for $`\lambda _L=515nm`$. However this much larger value is consistent with our measurements at $`\lambda _L=633nm`$ (Fig. 4). Moreover, from Eq. (5) we calculated theoretical curve $`\mathrm{\Delta }\theta `$ $`(\theta )`$ shown in Fig. 4. This curve fits very well the width enhancement of the CBS cone at $`\theta _B=0`$, when $`\theta `$ is changed for a fixed $`\lambda _L`$ (Fig. 4 for $`\lambda _L=\lambda _B=633nm`$). Eq. (5) is also consistent with the obtained results when $`\lambda _L`$ is changed for a fixed $`\theta =0`$ (Fig. 4 inset). The agreements obtained validate both our modified CBS model and the extracted optical parameters for the opal photonic crystal. We again emphasize that consistent values of $`L_B`$ were inferred from both the magnitude $`\mathrm{\Delta }\theta `$ $`(\theta _B)`$ and the width $`\gamma (\theta _B)`$ of the CBS anomaly.

The widths of CBS cones from $`\left(220\right)`$ planes measured with $`\lambda _L=365nm`$ are shown in Fig. 5, where $`\mathrm{\Delta }\theta `$ is plotted versus the incident angle $`\theta `$ measured from $`\left[220\right]`$ direction. For $`\lambda _L=365nm`$ the Bragg CBS anomaly in $`\mathrm{\Delta }\theta `$ occurs at $`\theta =\theta _B=\mathrm{cos}^1\left(\lambda _L/2d_{220}n_{eff}\right)20^{}`$, and is clearly seen in Fig. 5. Fitting this CBS anomaly in a similar way as that at $`\theta =42^{}`$ in Fig. 2, we deduced a Bragg length, $`L_B=4.3\mu m`$ for the photonic gap in $`\left[220\right]`$ direction. In our fitting we used the value of mean-free-path, $`l^{}=4.6\mu m`$ that was extracted from the widths of the CBS cones for incident light beams close to the normal direction to the $`\left(220\right)`$ planes. We note that, in addition to $`L_B`$ along $`\left[220\right]`$, $`l^{}`$ is also different from that extracted before from measurements along $`\left[111\right]`$.

## IV Microscopic origin of the scatterers

Apart from the Bragg-related anomalies we note that three qualitative features of the data presented above indicate that the description of CBS from opal photonic crystal within the conventional model of a medium with randomly positioned point-like scatterers is not completely adequate. These features are:

(i) The Rayleigh scattering law $`l^{}\lambda _L^4`$ that is characterisitc for point-like scatterers does not hold in opals. This fact becomes evident from Fig. 4 inset; except for a narrow $`\lambda _L`$ interval near $`\lambda _L\lambda _B`$, the CBS cone width $`\mathrm{\Delta }\theta `$ remains approximately unchanged with $`\lambda _L`$ for $`\lambda _L`$ interval between $`800nm`$ to $`360nm`$.

(ii) The value of $`l^{}=7.2\mu m`$ for $`\left[111\right]`$ direction appreciably exceeds $`l^{}=4.6\mu m`$ for $`\left[220\right]`$ direction.

(iii) The cone width $`\mathrm{\Delta }\theta `$ increases with $`\theta `$ for CBS along $`\left[111\right]`$ (Fig. 2) but shows no dependence on $`\theta `$ for CBS along $`\left[220\right]`$ (Fig. 5).

The question that we would like to address here is whether these CBS features can help in pinpointing the size, shape, and orientation of the generic defect scatterers in opal.

Firstly, independence of $`l^{}`$ on $`\lambda _L`$ suggests that the size of the defects is larger than $`\lambda _L`$. On the other hand, defects with size $`\lambda _L`$ scatter most of the light in the forward direction. More precisely, the scattered light intensity is restricted to a cone with aperture $`\lambda _L/\pi `$. As a direct consequence of this almost forward scattering, we should have $`l_{tr}^{}l^{}`$, where $`l_{tr}^{}`$ is the transport mean-free-path that determines the diffusion coefficient ($`D=cl_{tr}^{}/3`$) in this case.

The conclusion that $`l_{tr}^{}`$ is much bigger than $`l^{}`$ in opal has some dramatic consequencies. As can be seen from the general expression for the coherent albedo \[Eq. (1)\], the dependence $`\alpha _c(\delta \theta )`$ is determined by two factors: the surviving probabilities $`𝒫_\varphi `$, $`𝒫_\varphi ^{}`$ and the propagator $`\mathrm{\Pi }`$. For point-like scatterers both ingredients change within the same characteristic spatial scale, namely $`l^{}`$. However if $`l_{tr}^{}l^{}`$, then the scale of the propagator that describes the diffusion process is $`l_{tr}^{}`$, whereas the scale of $`𝒫_\varphi `$, $`𝒫_\varphi ^{}`$ is still $`l^{}`$. This would make the cone width $`\mathrm{\Delta }\theta \lambda _L/l_{tr}^{}`$ rather than $`\mathrm{\Delta }\theta \lambda _L/l^{}`$ as in the usual case. Consequently $`\mathrm{\Delta }\theta `$ would be independent of the incident angle $`\theta `$, and, correspondigly insensitive to the Bragg anomaly, contrary to the data.

A possible resolution of this apparent contradiction with the experiment that we can suggest here is to assume that the scatterers are highly anisotropic, such as disk-like with radius $`\lambda _L`$ and thickness $`\lambda _L`$. Moreover, in order to insure that CBS is different from the $`\left(111\right)`$ and $`\left(220\right)`$ opal surfaces as found experimentally, we have to assume that the disks are not oriented randomly, but rather their orientation is correlated. Although we did not extend the quantitative theory of CBS to this situation, intuitively it seems likely that the disks normals be parallel to the $`\left[111\right]`$ direction. For such an orientation the light would diffuse easier along the $`\left(111\right)`$ planes. Under these conditions of strongly anisotropic diffusion, longer paths do not necessarily lead to narrower cones as in the case of isotropic scatterers. Since CBS is essentially the “structure factor” of the propagation along the boundary, then the easier diffusion along $`\left(111\right)`$ would lead to narrower cones, i.e. to $`\mathrm{\Delta }\theta _{\left[111\right]}<\mathrm{\Delta }\theta _{\left[220\right]}`$, in agreement with the experimental results (Figs. 2 and 5).

A possible realization of the geometry and arrangement of scatterers as conjectured above can be provided by stalking faults defects. Such defects can be easily created along the $`\left[111\right]`$ opal fastest growing direction. Also, this type of defects has been recently shown to dominate disorder in colloidal crystals, and we note that synthetic opals grow from silica colloids.

## V Conclusions

In conclusion, we have demonstrated that CBS is a powerful technique for investigating photonic crystals with incomplete photonic band gaps. These measurements show that the remarkable universality of CBS characteristics for diverse disordered materials is not valid for photonic crystals at near Bragg conditions. Using the major effects on light transport near Bragg conditions and their influence on CBS cones, we deduce the values of several fundamental parameters of an opal photonic crystal. Remarkably, despite the existence of substantial scattering due to disorder, nevertheless we showed in this work that the Bragg attenuation length and, thus the photonic band gap width of an “ideal” crystal along $`\left[111\right]`$, as well as light mean-free-path can be readily extracted from the CBS measurements, even if they are comparable in magnitude (unlike the situation in Ref. ).

The reason why for $`L_Bl^{}`$ disorder is unable to completely mask the contribution of Bragg scattering to the CBS inside the photonic crystal, is that at Bragg conditions the CBS involves $`\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}`$ light propagation before the first and after the last scattering events, when the existing disorder is not yet relevant. This makes the CBS measurements a much more sensitive tool for characterizing the inherent band gap features for disordered photonic crystals as compared to the more conventional measurements of optical reflectivity.

The only other study of CBS from photonic crystals that we know of was recently reported in Ref. . The CBS measurements in that work were performed on photonic crystals based on polystyrene spheres and also on crystals composed of air spheres in a $`TiO_2`$ matrix. In contrast to the present study that has been performed on a single opal sample where the stop-band was “swept” by changing the angle of incindence, $`\theta `$; in Ref. , on the contrary the incident angle was kept close to normal. Under this condition, in order to explore the CBS within the stop-band frequency range in the polystyrene opals, the authors measured CBS on a set of samples with different sphere diameters using three different wavelenghts for the incident light beam. Also the band edge frequencies were determined from reflectivity measurements, rather than by CBS as in our work. Even under these assumptions, the authors failed to measure CBS cones inside the stop-band.

For air-sphere $`TiO_2`$ crystal, however, the authors used a setup allowing continuous tuning of $`\lambda _L`$. They observed that between the red and blue edges of the stop-band, $`\mathrm{\Delta }\theta `$ decreases by a factor of $`2`$. This observation, which actually contradicts our own study, was explained on the basis of a previous theoretical work. It was claimed that for light frequencies above the lower band-edge the electromagnetic wave, in course of the diffusion process, can occasionally “hit” the boundary interface at the Bragg angle $`\theta =\theta _B`$; thus rather than exiting the sample, it experiences instead the Bragg diffraction to the interior. As a result of such events the diffusive trajectories become longer (on average) leading to the narrowing of the CBS cone.

In light of the CBS model that was outlined in Section III, the above explanation corresponds to the modification of the propagator $`\mathrm{\Pi }(\rho ,z,z^{})`$ in Eq. (1). In contrast, our analysis suggests that $`\mathrm{\Pi }`$ is not appreciably affected by the Bragg conditions, whereas the prime effect of the photonic band structure is the modification of the surviving probability densities $`𝒫_\varphi (z)`$ and $`𝒫_\varphi ^{}(z^{})`$ in Eq. (1). The apparent contradiction between our results and those of Ref.  can be probably accounted for by the difference in the “strenghts”of the Bragg diffraction process, or $`L_B`$ in synthetic $`SiO_2`$ opals and air-sphere crystals in $`TiO_2`$. We note, however that the relative width of the Bragg reflectivity peaks in both samples is comparable $`\delta \lambda _B/\lambda 0.1`$.

###### Acknowledgements.

This work was supported by NSF grant DMR 973280, the Petroleum Research Fund under grant ACS-PRF #34302-AC6 and the Army Research Office.