text:
# CDF/PUB/STATISTICS/PUBLIC/6031 30 June 2002 Signal Significance in Particle Physics

## 1 What Particle Physicists Mean by Significance

When one of your colleagues approaches you and declares that she has made a “significant” observation, intuitively that means that she has observed some phenomena whose interpretation allows her to eliminate or falsify one or more hypotheses, and usually support one or a small number of alternative hypotheses. We furthermore expect that the observation has sufficient statistical power that we expect that additional observations are unlikely to change these conclusions. Scientists have therefore attempted to identify a consistent statistical framework in which we can quantify this concept of “significance.”

In particle physics, this concept of statistical significance has not been employed consistently in the most important discoveries made over the last quarter century. Examples of the major discoveries made over an approximately 10 year period between the late 1970’s and the late 1980’s illustrate this point.

Let us consider first the discovery of the $`\mathrm{{\rm Y}}`$ meson (and the $`b`$ quark) in 1977 by L. Lederman and colleagues . This was made through the observation of $`\mu ^+\mu ^{}`$ final states in high-energy proton-nucleus collisions at Fermilab, where a large resonant signal was observed on top of a steeply falling background of dimuon candidates. The experimenters estimated that they observed a signal of approximately 770 events on top of a non-resonant background of 350 candidates. They characterized the signal as “significant” but made no attempt to quantify or explain exactly what they meant.

The discovery of the $`W^{}`$ boson at CERN in 1983 by the UA1 collaboration  was made by observing 6 events produced in proton-antiproton collisions where a high energy electron or antielectron was observed in coincidence with a signature for a recoiling energetic neutrino. The collaboration estimated the background to these 6 events as being “negligible” and claimed discovery of the expected charged weak intermediate vector boson. This observation was subsequently confirmed by the UA2 collaboration.

The discovery of $`B`$ mesons in 1983 by the CLEO collaboration  was performed by carefully reconstructing a variety of different decay modes and illustrating an invariant mass peak at 5.4 $`\mathrm{GeV}/\mathrm{c}^2`$. The collaboration observed a total event rate of 17 events on a background of between 4 and 7 events. They claimed definitive observation of a new particle, but made no statement that quantified the statistical power of the observation.

As a final example, I note the discovery of $`B^0`$ meson flavour mixing in 1987 by the ARGUS collaboration . The experimenters observed $`24.8\pm 7.6\pm 3.8`$ unexpected same-sign dilepton events versus a total of $`25.2\pm 5.0\pm 3.8`$ opposite-sign dilepton candidates. They characterized this as a “3 $`\sigma `$” observation, namely, that the probability that the observed number of same-sign dilepton events could have been as great or greater than the observed value was equivalent to the probability of a Gaussian statistic being observed at least 3 standard deviations from its expected mean (a probability of $`1.35\times 10^3`$).

This brief review illustrates that quantifying the statistical significance or power in seminal particle physics measurements is not uniformly done. It also illustrates that in at least the one case in which it was done, the significance was defined as the probability222Unless otherwise noted, “probability” in this article refers to the frequentist definition of this concept. of the “null hypothesis” having been responsible for the observation.

In this paper, I will first review briefly the formal concept of “statistical significance.” I will then discuss several examples that illustrate the use of this concept in particle physics. I do not have the opportunity to review all of the techniques that have been in recent use, but refer the reader to other articles in these proceedings (for example, the review of the $`CL_S`$ method by A. Read).

## 2 Formal Definitions of Significance

### 2.1 The Frequentists Perspective

The concept of statistical signficance is formally introduced in the context of hypothesis testing . Suppose that we have two hypotheses, $`H_0`$ and $`H_1`$, and a measurement whose value is a test statistic $`X`$ that, as a random variable, provides some discrimination between these two hypotheses. Let $`f_0(X)`$ and $`f_1(X)`$ represent the probability distribution functions for $`X`$ associated with the two hypotheses.

Prior to making a measurement of $`X`$, we would identify a “critical region,” $`w`$, such that we would select hypothesis $`H_1`$ if $`Xw`$ and $`H_0`$ otherwise. We now have four possible outcomes when we make a measurement of $`X`$. If $`Xw`$ and the hypothesis $`H_1`$ is true, then we have selected the correct hypothesis. If $`Xw`$ and $`H_0`$ is true, then we have incorrectly concluded that $`H_1`$ is true. This is known as a mistake of the first kind, and the probability for this decision is

$`{\displaystyle _{Xw}}f_0(X)𝑑X`$ $`=`$ $`\alpha .`$ (1)
The probability $`\alpha `$ is known as the “significance” of the test.

We have two other possibilities. The first is if we measure $`Xw`$ when $`H_0`$ is true. In this case, we would have made the correct inference. Finally, we have the case where $`Xw`$ and $`H_1`$ is true. This is known as a mistake of the second kind, and the probability for that decision is

$`{\displaystyle _{Xw}}f_1(X)𝑑X`$ $`=`$ $`\beta .`$ (2)
The probability $`1\beta `$ is known as the “power” of the test. The situation is illustrated in Fig. 1a). The significance $`\alpha `$ is therefore a measure of the ability of a test to avoid mistakes of the first kind, whereas the power $`1\beta `$ measures the ability of a test to avoid mistakes of the second kind. In defining an “optimimum” test, one would like to choose $`X`$ and the region $`w`$ such that $`\alpha `$ and $`\beta `$ are as small as possible.

### 2.2 Significance in Particle Physicists– The P-Value

The statistical definition of significance is made in the context of choosing between two hypotheses. However, the use of significance in particle physics discoveries is in a different context. The typical case is that an experiment makes a measurement of the test statistic $`X`$, say $`X_0`$. Furthermore, using the same notation as before, we assume $`X`$ has a probability density $`f_0(X)`$ if the hypothesis $`H_0`$ is true. We further assume that we can categorize observations of $`X`$ into those that are more and less consistent with $`H_0`$ (for the sake of discussion, I will assume that values of $`X`$ greater than $`X_0`$ are less likely given $`H_0`$). A measure of the inconsistency of the observed value $`X_0`$ with the hypothesis $`H_0`$ is then the probability

$`{\displaystyle _{X>X_0}}f_0(X)𝑑X.`$ (3)
This is formally identical to the definition of the significance in Eq. 1 if we now define the critical region to be $`w=\{X|X>X_0\}`$, i.e. the probability of observing a value of $`X`$ equal to or greater than our observation. This probability is formally known as the “p-value” of the observation, a convention that the Particle Data Group now has adopted . The advantage of using the formal term for this quantity is that it avoids confusion with the concept of significance defined in hypothesis tests, where the region $`w`$ is defined a priori, i.e. before the measurement is made.

The p-value for a given measurement and a specific hypothesis has a number of features. First, it only depends on the measurement and the probability density for the hypothesis. It is not a hypothesis test. It only provides a measure of the consistency of the hypothesis and the measurement. In that sense, it is most often quoted when one has made a measurement that appears to be inconsistent with a single hypothesis. A very small p-value is then used to support the inference that the specific hypothesis should be rejected.

Referring then to the example of the discovery of $`B^0`$ mixing given in Section 1, we can now say that the p-value for the observation for the non-mixing hypothesis was $`1.35\times 10^3`$. From a frequentist perspective, if one rejected the non-mixing hypothesis at this p-value and it was always true, then one would expect to be wrong (i.e., reject the correct hypothesis) on average 1 out of every 740 times.

### 2.3 A Few More Comments on Hypothesis Testing

Although the literature uses the p-value of an observation as a measure of its statistical significance, the concept of hypothesis testing is an important one in particle physics. One sees it most often used in the context where one is designing or proposing an experiment and wishes to characterize the experiment’s ability to distinguish existing and known physics phenomena (such as that predicted by the Standard Model) from possible new physics .

In those cases, a crucial aspect of the experiment design is the selection of the optimal statistic $`X`$ and the optimal critical region. For a specific measurement, such as the observation of a process above some expected background rate, the choice of $`X`$ will depend on the measurement and the ingenuity of the experimenter. One would like to identify test statistics that have quite different probability density functions for the hypotheses you wish to distinguish in order to be able to define a critical region with the smallest possible $`\alpha `$ and $`\beta `$. At the same time, the decision should be informed by the effect of any systematic uncertainties that may degrade the separation between two hypotheses. To that extent, one often attempts to identify statistics that are not expected to be affected by systematic uncertainties.

### 2.4 Neyman-Pearson Theorem

Hypothesis testing has an extensive literature, but relatively few general results have been identified that can guide our judgement. One result, known as the Neyman-Pearson Theorem, is surprisingly useful, and it is worth reviewing here for the insight it provides.

Suppose we have two hypotheses, $`H_0`$ and $`H_1`$, and we have defined a test statistic $`X`$. Then for a given significance $`\alpha `$, we can define the region $`w`$ which gives us optimal $`\beta `$ (i.e., the smallest value of $`\beta `$ and therefore the greatest statistical power) by choosing $`w`$ as follows. We first form the ratio of probability density functions for the two hypotheses

$`I_N(X)`$ $``$ $`{\displaystyle \frac{f_0(X)}{f_1(X)}}.`$ (4)
The Neyman-Pearson Theorem then concludes that the optimal region $`w`$ is the one over which $`I_N(X)`$ is maximal, namely that we find the value $`c_\alpha `$ such that when

$`w=\{X|I_N(X)<c_\alpha \},`$ (5)
the probability of observing $`Xw`$ is

$`{\displaystyle _{Xw}}f_0(X)𝑑X`$ $`=`$ $`\alpha .`$ (6)
This construction is illustrated in Fig. 1b).

The Neyman-Pearson test has one signficant limitation–it is only valid for what are known as “simple hypotheses,” or hypotheses where there are no unknown parameters that would be estimated from the data. In addition, since it is only applicable when comparing two hypotheses, it cannot be employed in cases where you have multiple alternative hypotheses to consider. However, despite these limitations, this theorem gives us considerable insight into the definition of the critical region. For example, we can relate the ratio of probabilities to the ratio of likelihoods of the two hypotheses:

$`I_N(X)`$ $``$ $`{\displaystyle \frac{f_0(X)}{f_1(X)}}{\displaystyle \frac{L_0(X)}{L_1(X)}},`$ (7)
where $`L_i(X)`$ are the likelihood functions defined for the two hypotheses $`i=0`$ and $`i=1`$. This suggests that the likelihood ratio is one source of guidance for defining critical regions that have significant (though perhaps not optimal) power.

## 3 The Bayesian Perspective

Our considerations up to this point have been from a frequentist perspective, using the standard definition of a frequentist probability. In calculating a p-value for a measurement, one has to assume a hypothesis and then determine the probability (or probability density) for all possible outcomes of the measurement.

A Bayesian statistician does not consider data other than the single measurement. However, for each hypothesis, the Bayesian could define a credibility interval that reflects his or her degree-of-belief in each hypothesis, and the ratio of these credibility intervals–what is called the “Bayes discriminant factor”– becomes a measure of the relative confidence one has in the two hypothesis. Formally, this ratio is

$`{\displaystyle \frac{P(H_0|X)}{P(H_1|X)}}={\displaystyle \frac{L_0(X)}{L_1(X)}}\dot{{\displaystyle \frac{\pi _0(X)}{\pi _1(X)}}},`$ (8)
where $`\pi _0(X)`$ and $`\pi _1(X)`$ are the prior probabilities associated with the two hypotheses.

This ratio can be used to reject one of the two hypotheses. The Bayesians would argue that there is no benefit in attempting to make anything but a relative statement about the degree-of-belief of the two hypotheses. Thus, there is no direct analogy to the p-value in this framework. The advantage of this perspective is that it avoids the need to understand the probability density of all possible outcomes for a given hypothesis. It also has the advantage that any inferences you draw are less sensitive to an outcome that has a low probability regardless of the hypothesis. In such cases, the Bayes discriminant factor still provides information, whereas the p-value is no longer very informative and could in fact be misleading.

The disadvantages with this Bayesian approach are, however, that one has to assume prior distributions for each hypothesis, and one is only allowed to make relative confidence statements about two hypotheses. For these reasons, one finds very limited use of the Bayes discriminant factor in particle physics.

## 4 P-Values and Experimental Design

The definition of significance in terms of a p-value for an observation immediately makes clear the importance of a priori decisions on the random variables one will measure and how one will define those observations that prefer one hypothesis over another. A carefully designed experiment will identify these and optimize their choice before any data is analyzed.

However, many particle physics experiments make unique measurements using general-purpose apparatus designed to study a large class of processes. Thus, it is difficult, and often impossible, to anticipate what one will observe and how. In fact, early studies of the data will often guide the experimenters to focus in specific features that appear unusual or unexpected. In this context, the evaluation of a p-value may prove very difficult.

A simple example illustrates this problem. Suppose one measures an invariant mass spectrum in a specific region, say $`[m_1,m_2]`$, and one observes a narrow enhancement in a small mass interval, say $`\mathrm{\Delta }m`$ wide, of $`N_o`$ events above an expected background of $`N_b`$ events. In this case, it would be natural to assume that the hypothesis we wish to test is the “null” hypothesis where we expect $`N_b`$ events in this mass interval $`\mathrm{\Delta }m`$ and then determine the probability of observing at least $`N_o`$ events. Assuming that the background rate is well known (and so we can ignore its uncertainty), the p-value for this observation would be given by the Poisson probability for observing at least $`N_o`$ events when the mean rate is $`N_b`$, or

$`\alpha ={\displaystyle \underset{n=N_o}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{exp}\left(N_b\right)\left(N_b\right)^n}{n!}}.`$ (9)
However, this probability does not take into account the fact that we are considering all possible choices of mass interval $`\mathrm{\Delta }m`$ in the region $`[m_1,m_2]`$.

A proper estimate of this p-value would then have to include the likelihood of observing at least $`N_o`$ events in any possible interval $`\mathrm{\Delta }m`$. This increases the p-value of the observation, and changes the possible inferences one can make. For example, a Monte Carlo calculation where $`\delta m`$ is 1% of the interval, $`N_o=8`$ events and $`N_b=100`$ (i.e., the average number of events in any $`\delta m`$ interval is one) gives a p-value that is 500 times larger than the result in Eq. 9.

### 4.1 Blind Analyses

The prevalence of the p-value in making inferences rests on the assumption that it is possible to estimate the frequency of all observations of the test statistic, and that it is possible to identify the class of observations that are less consistent with a given hypothesis (the critical region in the language of hypothesis testing). This is inherently difficult in cases where one allows the definition of test statistic and critical region to depend on the actual experimental outcome itself. A tactic to eliminate such bias is the “blind analysis,” where one defines the critical region and the statistic without knowledge of the relevant data .

The ideal experiment is one in which the measurement and any calculation of its p-value does not have to be informed by the data itself. No choices with regard to selection of data, modifications in the test statistic or choice of critical region would then be allowed once data collection has started. This approach avoids the possibility of selecting, consciously or unconsciously, a critical region or test statistic that tend to favour or disfavour a given hypothesis based on the data observed.

A number of celebrated failures of inference in particle physics over the last half-century illustrate what happens when the experimenter allows the data to guide his or her choices in making inferences about data . In all these cases, the quoted p-value has been assessed incorrectly because it has failed to take into account how the frequency of a given observation would be affected by making choices about the test statistic and critical region based on the actual distribution of the data itself.

### 4.2 Use and Limitations of Blind Experiments

The simple example of “bump hunting” illustrates the fundamental problem in particle physics where one is searching for evidence of new phenomena; it is inherently difficult to identify a priori what class of observations one would expect to use in such a search. Besides the difficulty of defining in advance all possible means of separating “signal” from “background,” it is also difficult to limit access to data when one also has to verify that the instrumentation is working correctly and that any artefacts created by effects such as miscalibration and errors in bookkeepping are identified and mitigated. The experiment design also has to allow the experimenter access to the data to measure the rate of background events in the signal sample.

Despite these challenges, the elimination of certain biases that are otherwise difficut to control make a blind analysis an attractive approach given the benefits of being able to make straightforward estimates of p-values for the possible outcomes. This technique is reviewed in another contribution to these proceedings .

## 5 P-Values in a Counting Experiment

### 5.1 General Considerations

A common particle physics experiment involves the search for new phenomena by observing a unique class of events in particle interactions that cannot be described by background hypotheses. One usually can reduce this problem to that of a “counting experiment,” where one identifies a class of events using well-defined criteria, counts up the total number of observed events, $`N_o`$, and estimates the average rate of events, $`N_b`$, that come from the various background processes. One can then perform a straightforward estimate of the p-value of a given observation of $`N_o`$ events, assuming that the probability density for the random variable $`N_o`$ follows a Poisson distribution, i.e. the formula in Eq. 9.

There are several issues that even this simple problem has to address. First, one has to be sure that the criteria used to select the class of events was not in itself biased by how $`N_o`$ varied as the criteria were modified. Here is where a blind analysis has its greatest benefit, since this bias is explicitly guarded against. Second, one has to take into account possible uncertainties in the estimate of the background rate $`N_b`$. It is beyond the scope of this article to discuss this issue, and the interested reader is referred to the growing literature on this topic  (a typical frequentist approach is to extend the ensemble of possible measurements to include those experiments with different values of $`N_b`$ consistent with the knowledge of $`N_b`$). Third, the careful experimenter has to make sure that all information relevant to the search is used in the measurement. It is at best inefficient and at worst misleading to ignore relevant data (for example, a possible channel in which the number of observed events can provide additional information on the process being studied).

As a concrete example of the calculation of a p-value for a typical counting experiment, I will summarize the techniques used by the CDF and D0/ collaborations in their search for top quark production.

### 5.2 The Top Quark Search

The top quark was discovered by pair-production in proton-antiproton collisions at an energy of 1.8 TeV . The top quark decays predominantly via the process $`tWb`$, with the $`W`$ boson subsequently decaying either leptonically via $`Wl\nu _l`$ (where “$`l`$” can be either an electron, muon or tau lepton) or hadronically via $`Wq\overline{q}^{}`$ (the quark final states are either $`u\overline{d}`$ or $`c\overline{s}`$). This results in three categories of possible final states with different topologies, efficiencies and background rates:

1. the lepton+jets channel, involving one high energy lepton, a neutrino and three or more jets from the hadronic decay of the $`W`$ and the $`b`$ quarks,
2. the dilepton channel, involving two high energy leptons, evidence for two neutrinos, and two or more jets from the $`b`$ quarks, and
3. the hadronic channel, involving at least six jets.

In both experiments, one had to use additional criteria to improve the signal-to-noise ratios in the final candidate event samples. For CDF, the most effective way to do this was to require evidence that at least one of the jets arose from a $`b`$ quark using two different “b-tagging” techniques. Thus, one could characterize the final states by the number of $`b`$ tags, with the events with one or two $`b`$ tags having increasing purity. For D0/, the most effective way to reduce backgrounds was by imposing more stringent kinematic criteria (a topological selection) and using soft muon $`b`$-tagging.

The searches used data samples of increasing sensitivity. The first reported data came after the CDF and D0/ collaborations had recorded 19.6 and 15.0 $`\mathrm{pb}^1`$, respectively . At that time, the experiments had not completed analysis of the hadronic channels, which were expected to be dominated by background. The results of these analyses are summarized in Table 1, where we list the number of observed events, the estimated background rates, and the branching ratio times efficiency of observing a $`t\overline{t}`$ decay in each mode.

The collaborations evaluated the statistical significance of their data by using a Monte Carlo calculation to estimate the frequency that the expected background processes would create a combined signal that was at least as large as that observed. The Monte Carlo calculation created an ensemble of experiments that modelled the possible observations in all channels assuming the Standard Model background hypothesis. For a given channel, the estimated background rate was used as the mean of a Poisson distribution of observed events. In order to account for uncertainties in the background rate, the mean value used to generate a new member of the ensemble was obtained by sampling a Gaussian distribution with the mean and width of the estimated background rate . The results of these p-value calculations are summarized in Table 2. The collaborations concluded that the individual observations did not provide sufficient evidence to exclude the background hypothesis.

The collaborations proceeded to determine how likely their set of observations were assuming the background hypothesis by identifying a statistic that combined the observations in the individual channels. In the case of a counting experiment involving several channels, the maximum likelihood estimate of the rate of the process is simply the sum of the event rates in each channel. Thus, the natural statistic to evaluate the combined significance of the observations was the observed sum of events in all channels. However, the CDF collaboration noted that the most sensitive measure of the cross section was not the total number of observed events in their sample, but the total number of observed $`b`$-tags (since there was a much larger probability of observing two $`b`$-tags in a signal event than in an event from a background process). Thus, CDF chose as its statistic the sum of the number of $`b`$-tags in the lepton+jet events combined and the number of dilepton events. Since the D0/ data relied less on $`b`$-tagging, the collaboration chose to use the total number of observed events.

The calculation of the p-value of the observation assuming the background hypothesis was performed by a Monte Carlo procedure that effectively created a set of “pseudo-experiments.” In each pseudo-experiment, the number of $`b`$-tags and dilepton events from the different background sources was drawn from a Poisson distribution that had as its mean value the estimated background rate for the process. The uncertainty in the various background components was taken into account as described above, as was the correlation in the different background sources. This correlation arose from the fact that a number of background sources contributed both types of $`b`$-tags, whereas others did not. In effect, this increased the frequency of observing a larger number of $`b`$-tags (since now the fluctuations in the two components were correlated).

The resulting p-values are summarized in Table 2. One sees that the single most significant p-value was $`2.6\times 10^3`$. If one had not taken into account the correlations between the background sources, the combined p-value would have been $`1.6\times 10^3`$, or a factor of almost two smaller. Alternatively, the combined p-value determined by just counting events would have been approximately $`10^2`$. This demonstrates the sensitivity of a p-value calculation to the approximations used to determine it. Given all this information, both experiments concluded that the observations were not sufficiently compelling statistically to exclude the background hypothesis.

### 5.3 Significance Required for Discovery

In the search for the top quark, the CDF and D0/ collaborations argued that observations with p-values of order $`10^3`$ were not sufficiently significant to be used to claim discovery of a new phenomenon. Although this is clearly a matter of opinion, it is roughly consistent with the practice in the field, where typically the “$`5\sigma `$” standard is used as rough rule of thumb to define the sensitivity necessary for discovery. This corresponds to a p-value equivalent to between $`5.7\times 10^7`$ and $`2.8\times 10^7`$, depending on whether you are searching for a deviation from a mean or a one-sided fluctuation from the mean.

As a concrete example, the two Tevatron collaborations used an identical analysis procedure when approximately a factor of two more data had been recorded by both experiments. The resulting p-values of the CDF and D0/ observations assuming the background hypothesis were $`1\times 10^6`$ and $`2\times 10^6`$, respectively . Both experiments concluded that the background hypothesis could be excluded and claimed observation of top quark pair production.

## 6 P-Values for Continuous Test Statistics

High-energy physics measurements often examine statistical variables that are continuous in nature. In fact, to identify a sample of events enriched in the signal process, one often imposes selection requirements on such continuous variables. Often, it is important to take into account the entire distribution of a given variable for a set of events, and not just whether the events lie in a given range of values.

The general problem can be posed in the following way. Suppose we have a set of event data each characterized by a set of statistics $`\stackrel{}{X}_i`$, where $`i=1`$ to $`N`$. In addition, one has a hypothesis to test that predicts the distribution of $`\stackrel{}{X}`$, say $`f(\stackrel{}{X};\stackrel{}{\alpha })`$, where we assume this function to be normalized to unity between $`X_{min}`$ and $`X_{max}`$, the minimum and maximum values of $`X`$, and $`\stackrel{}{\alpha }`$ is a set of parameters that are either known or estimated directly from the data. Then the general problem is to define a statistic that gives a measure of the consistency of the distribution of data with the distribution given by the hypothesis.

### 6.1 Possible Tools

The most widely used such statistic in the 1-dimensional case is a form of a “runs test,” which compares the predicted cumulative distribution

$`g(X)={\displaystyle _{X_{min}}^X}f(X^{})𝑑X^{}`$ (10)
with the observed cumulative distribution $`h(X)`$. The most common test is the Kolomogorov-Smirnov (K-S) test , which makes this comparison by first finding the K-S distance

$`\delta =\mathrm{max}\{|g(X)h(X)|,X(X_{min},X_{max})\},`$ (11)
namely the largest difference between the two cumulative distributions. This test statistic has a characteristic distribution that can be calculated analytically to provide one with a p-value, specifically the probability that one would observe a value of this test statistic as large as or larger than the observed value.

The K-S test gives a distribution-free measure of the consistency of a 1-dimensional continuous variable and is often used in the particle physics literature. Although there are a number of other tests that could be used in this case, all with similar properties , the K-S test has become a reference standard to employ.

### 6.2 Extension to Higher Dimensions

The K-S test (and other runs tests) are in principle limited to 1-dimensional distributions, but there are extensions to the case of several dimensions, though with a number of restrictions. The extension requires one to assume that the probability distribution predicted by the hypothesis can be factorized, so that

$`f(\stackrel{}{X})=f_1(X_1)f_2(X_2)\mathrm{}f_n(X_n),`$ (12)
where $`n`$ is the number of continuous variables being compared. This in effect requires each of the variables to be uncorrelated, a strong assumption and one that has to be verified in practice. With this assumption, however, one can then define a set of independent statistics $`\delta _i`$, $`i=1`$ to $`n`$, and the associated p-value for each observed K-S statistic $`p_i`$. Then one can combine these independent p-values into a single measure of significance.

### 6.3 Example: CDF “Superjets”

A concrete example of this technique is a recent analysis of hadron collider data performed by the CDF collaboration. A study was performed of events that were consistent with the production of one or more hadronic jets and a $`W`$ boson decaying to a lepton-neutrino pair. The collaboration defined a subsample of these events where at least one jet was identified as a “superjet”, namely a $`b`$-quark candidate jet with both the presence of a secondary vertex in the jet displaced from the interaction vertex and the presence of a second lepton associated with the jet consistent with coming from the semileptonic decay of a $`b`$ hadron .

The collaboration found 13 such events in the 1992-96 Tevatron Collider data, where they estimated that they would have expected $`4.4\pm 0.6`$ events from Standard Model background sources. This observation has a p-value of 0.001, treating it as a counting experiment and using the techniques introduced above. The authors then proceeded to examine nine separate kinematic variables that had distributions that were predicted to be largely uncorrelated, but that might distinguish between the Standard Model backgrounds and a variety of exotic sources of events. A typical example of such a comparison is given in Figure 2, where the observed distribution of the lepton pseudorapidity ($`\eta \mathrm{ln}\mathrm{tan}(\theta /2)`$, where $`\theta `$ is the angle of the lepton relative to the incoming proton beam axis) is compared with the predicted $`\eta `$ distribution.333The authors chose background distributions for these figures obtained using Monte Carlo calculations, but used background distributions for their p-value calculations obtained by “bootstrapping,” using a complementary data sample that had no signal events and that was argued to provide a good characterization of the expected Standard Model backgrounds. The Standard Model Monte Carlo calculations resulted in similar p-value estimates. The plots on the right-hand side are the distributions of the K-S distance as determined from a Monte Carlo calculation.

The p-values from each of the distributions were determined and range from 0.001 to 0.15. The authors comment that “given the a posteriori selection of the 9 kinematical variables, the combined statistical significance cannot be unequivocally quantified.” However, we can determine a combined p-value by calculating the product of the 9 p-values, $`p_{tot}`$, and determining how likely it would be to obtain this product value assuming the background hypothesis. This is given by

$`P_{tot}={\displaystyle \underset{m=1}{\overset{9}{}}}\left[{\displaystyle \underset{k=0}{\overset{m1}{}}}{\displaystyle \frac{(\mathrm{ln}p_{tot})^k}{k!}}\right],`$ (13)
and equals $`1.6\times 10^6`$ assuming you set aside the reservations of the authors.

This estimate of the overall p-value raises a number of comments. First, are the variables sufficiently uncorrelated that any residual correlations can be ignored? Various tests were made of this assumption by the authors, but no rigorous argument was presented. Second, uncertainties in the Standard Model predictions have not been incorporated into the p-value calculation. These may have some effect on the overall result, but it is unclear how large this might be. Third, the effect on the p-value estimate of the a posteriori choice of variables is virtually impossible to assess. A study of a series of alternate variables were made by the authors, but no firm conclusion could be drawn.

Of these, perhaps the third is the most vexing. It is true that the choice of the 9 variables for this analysis was made after the 13 event data sample had been identified as being unusual. In that sense, it is no longer possible to argue that the quoted p-value is an unbiased measure of the significance of the observation.

In this case, the best strategy is to repeat the measurement with an independent data sample to determine if the same effect is observed. However, this analysis serves as a good example of the issues one must face in making such a multi-variate estimate of significance.

## 7 Observations on Current Practice and Summary

Particle physicists have increasingly relied on numerical estimates of statistical significance. The literature is replete with the use of the p-value, and this appears to have developed into one common measure, as illustrated by the examples provided above. Other measures of significance are often quoted, such as the equivalent number of standard deviations a measurement lies from the value predicted by a hypothesis. This is, of course, just a p-value under a different name.

More significantly, there are consistent attempts in the literature to include in p-value estimates more complete information about a given measurement, such as the sensitivity of the estimate to systematic uncertainties and information from several statistics. The more difficult problem of avoiding unconscious bias in the selection of statistics is addressed through the use of “blind analyses,” but the effective application of such techniques to truly serendipitous discoveries is problematic. Here, the time-honoured technique of testing specific hypotheses developed through the study of one data set by creating and analyzing an independent data set with at least comparable statistical power remains the most effective tool for separating what we would call the “statistical fluctations” from first evidence for truly new phenomena.

Finally, what is an appropriate criteria for claiming a discovery on the basis of the p-value of the null hypothesis? The recent literature would suggest a p-value in the range of $`10^6`$, comparable to a “$`5\sigma `$” observation, provides convincing evidence. However, the credibility of such a claim relies on the care taken to avoid unconscious bias in the selection of the data and the techniques chosen to calculate the p-value.


format: text

text:
# Elementary excitations of trapped Bose gas in the large-gas-parameter regime

## I Introduction

The mean-field Gross-Pitaevskii (GP) theory gp  has been quite successful in explaining both static and dynamic properties dalfovo ; fetterrv ; castinrv  of Bose-Einstein condensates (BEC) produced in alkali atoms anderson ; bradley ; davis  . The reason for the success of the GP theory is mainly the satisfaction of dilute gas condition $`na^3<<1`$ ( where $`n`$ is atomic density, $`a`$ is the s-wave scattering length of interatomic potential and the parameter $`na^3`$ is called the gas parameter) in the above mentioned experiments. Typically the values of gas parameter in these experiments were in the range of $`10^510^4`$. The gas parameter can have a large value when the number of atoms $`N`$ in the condensate or the scattering length is large. Thus in the later experiments reported in Refs. ketterle ; stamperkurn  $`N`$ was of the order of $`10^8`$ and the corresponding value of gas parameter was $`10^3`$. On the other hand, in a very recent experiment cornish  the maximum value of the peak gas parameter was $`n(0)a^310^2`$ ($`n(0)`$ is the peak density of the condensate) although the condensate contained only $`10^4`$ 85 Rb atoms. In this experiment the scattering length $`a`$ could be varied continuously from a negative to a very large positive value by exploiting the Fesbach resonance. These kind of condensates, then, naturally raise the question about the validity of the mean-field GP approach. Theoretical investigations of the effects of going beyond the GP (BGP) theory have already been reported in the literature timmermans ; braaten ; nunes ; polls ; banerjee . Recently fabrocini  ground-state properties of BEC with large values of gas parameter (as achieved in the experiments cornish ) have been theoretically studied to assess the accuracy of the GP theory. It has been shown that for large values of the gas parameter ($``$ $`10^2`$) the total energy and the chemical potential are modified significantly by the BGP theory. It is well known that these changes are also reflected in the frequencies of collective excitations which can be measured with a greater accuracy. Therefore, by measuring the frequencies of collective excitations it should be possible to observe the effect of the BGP theory. This has motivated us to calculate the corrections to the frequencies of collective oscillations due to the BGP theory. We focus our attention on the condensates with values of gas parameters, spanning a range similar to that achieved in the experiment cornish  and make an assessment of the accuracy of the GP theory in predicting the frequencies of collective oscillations of such condensates.

In the past corrections to the frequencies of collective oscillations of trapped bosons arising due to the BGP theory have been estimated pitaevskii ; braaten2  within the Thomas-Fermi (TF) approximation which is valid in the large $`N`$ limit (more precisely when $`Na/a_{h0}>>1`$ is satisfied, where $`a_{h0}`$ is the characteristic length of trapping potential and it is defined in the next section). In the aforesaid references corrections to the frequencies were obtained by solving the time-dependent second-order differential equation perturbatively for the density. Thus the analytical expressions for the fractional change in the frequencies are correct only up to the first order in $`\sqrt{n(0)a^3}`$. It is, then, expected that these expressions for the corrections will not be accurate in the large-gas-parameter regime, as achieved in the recent experiment cornish . In this paper we also derive an analytical expression for the correction within the TF approximation by making use of the sum-rule approach. We find that the expression for the correction obtained by us is valid for somewhat larger range of gas parameters than those of Refs. pitaevskii ; braaten2 . As expected, for very small gas parameters the expression derived by us correctly reduces to those of Refs. pitaevskii ; braaten2 . We will return to the question of validity of these analytical results later.

It is well known that the TF approximation does not reproduce the surface region of density profile of trapped condensates accurately dalfovo . On the other hand, the collective oscillation frequencies characterizing the response properties of the trapped gas crucially depend on the tail or the surface region of the density profile bertsch . Therefore, it becomes necessary to go beyond the TF approximation especially for the calculation of frequencies of collective oscillations braaten2 . To accomplish this task we make use of well established variational method banerjee ; byam ; fetter1 ; mp  to determine the ground-state densities. The collective oscillation frequencies are then obtained by using the sum-rule approach bohigas ; lipparini  of many-body response theory. In the sum-rule approach the calculation of frequencies requires the knowledge of the ground-state wave function (or the ground-state density) of the many-body system only. Therefore the accuracy of frequencies determined by the sum-rule approach depends on the accuracy of the ground-state wave function (ground-state density) employed for the calculations. Keeping this in mind our aim is also to obtain reasonably accurate ground-state densities of the trapped condensates. To this end we employ variational approach by using an energy functional of $`N`$ interacting trapped bosons within the local density approximation (LDA) along with a judiciously chosen variational ansatz for the ground-state density banerjee ; mp .

The paper is organized in the following manner. In Section II we describe the theoretical method for the calculation of the frequencies of collective oscillations. Section III contains the analytical estimate of the corrections to the frequencies within the TF approximation, and results of our variational calculation. The paper is concluded in Section IV.

## II Variational Sum-rule Approach

The central quantity in the variational approach is the energy functional for a condensate of $`N`$ bosons, each with mass $`m`$, confined in a trap potential $`V_t(𝐫)`$:

$$E[n]=𝑑𝐫\left[\frac{\mathrm{}^2}{2m}\sqrt{n(𝐫)}^2\sqrt{n(𝐫)}+V_t(𝐫)n(𝐫)+ϵ(n)n(𝐫)\right],$$
(1)
where, the first, the second, and the third terms represent the kinetic energy of bosons, the energy due to the trapping potential, and the energy due to interatomic correlation within the LDA, respectively. Following the usual practice we use the perturbative expansion for $`ϵ(n)`$ in terms of the gas parameter $`na^3`$,

$$ϵ(n)=\frac{2\pi \mathrm{}^2an}{m}\left[1+\frac{128}{15\sqrt{\pi }}\left(na^3\right)^{\frac{1}{2}}+8\left(\frac{4\pi }{3}\sqrt{3}\right)(na^3)\mathrm{ln}\left(na^3\right)+𝒪\left(na^3\right)\right].$$
(2)
The first term in the above expansion, which corresponds to the energy of the homogenous Bose gas within the mean-field theory as considered in the GP theory, was calculated by Bogoliubov bogoliubov . The second term was obtained by Lee, Huang and Yang (LHY) lee , while the third term was first calculated by Wu wu  using the hard-sphere model for the interatomic potential. Although, it has been emphasized in the literature that the above expansion is valid only for $`na^3<<1`$, it is only recently that the range of validity of above expansion has been systematically investigated by Giorgini et al. giorgini . They have used diffusion Monte Carlo (DMC) method to calculate the ground-state of uniform gas of bosons interacting through different model potentials with $`a>0`$. It has been found that the expansion (2) is valid as long as $`na^3<10^3`$. However, for the values of the gas parameter beyond $`10^3`$ the inclusion of the logarithmic term leads to severe mismatch with the DMC simulation results. On the other hand, expansion (2) up to the LHY term gives an accurate representation of the DMC calculations even for the gas parameter of the order of $`10^2`$. Consequently, we do not consider the logarithmic term in the expansion (2) for all the calculations in this study.

The trapping potential $`V_t(𝐫)`$ is taken to be axially symmetric characterized by two angular frequencies $`\omega _{}^0`$ and $`\omega _z^0`$ ($`\omega _x^0=\omega _y^0=\omega _{}^0\omega _z^0`$). It is given by

$$V_t(𝐫)=\frac{m\omega _{}^{0}{}_{}{}^{2}}{2}\left(x^2+y^2+\lambda _0^2z^2\right),$$
(3)
where $`\lambda _0=\omega _z^0/\omega _{}^0`$ is the anisotropy parameter of the trapping potential ($`\lambda _0=1`$ corresponds to a spherically symmetric trap).

The ground-state properties which are needed to calculate the collective frequencies are determined by minimizing the energy functional (1) with respect to the density $`n(𝐫)`$hk  satisfying the normalization condition

$$n(𝐫)𝑑𝐫=N.$$
(4)
To implement the variational scheme we take the ansatz for the ground-state density as banerjee ; mp :

$$n(𝐫_1)=\frac{p}{2\pi \mathrm{\Gamma }(\frac{3}{2p})}\lambda ^{\frac{1}{2}}\left(\frac{\omega _{}}{\omega _{}^0}\right)^{\frac{3}{2}}e^{\left(\frac{\omega _{}}{\omega _{}^0}\right)^p\left(r_1^2+\lambda z_1^2\right)^p},$$
(5)
where $`\lambda `$, $`\omega _{}`$ and $`p`$ are the variational parameters which are obtained by minimizing $`E[n]`$ with respect to these parameters. In the above equation we have used the scaled co-ordinates $`𝐫_1=𝐫/a_{h0}`$ with $`a_{ho}=\left(\mathrm{}/m\omega _{}^o\right)^{\frac{1}{2}}`$. We note that $`n(𝐫_\mathrm{𝟏})`$ is normalized to unity and from Eq. (4) it can be seen that

$$n(𝐫)=\frac{N}{a_{h0}^3}n(𝐫_\mathrm{𝟏}).$$
(6)
We make this choice of variational form for the density as it has been shown convincingly that this ansatz is capable of describing the ground state of Bose gas in a trap quite accurately both within mp  and beyond banerjee  the GP theory for a wide range of particle numbers. Using this variational form the physical observables can be expressed analytically in terms of three variational parameters banerjee . The total energy $`E_1=E/\mathrm{}\omega _{}^o`$ can be written as banerjee

$$\frac{E_1}{N}=T+U+E_{int}^1+E_{int}^2.$$
(7)
Here $`T`$ and $`U`$ denote the average kinetic and the trapping potential energies per particle, respectively. The term E$`{}_{int}{}^{}{}_{}{}^{1}`$ gives the interaction energy per particle in the mean-field approximation as considered in the GP theory while E$`{}_{int}{}^{}{}_{}{}^{2}`$ gives the correction in the interaction energy due to the LHY term in the expansion (2) arising from the BGP theory.

The ground-state energy components obtained from the variational calculation are employed to compute the frequencies of collective oscillations by using the sum-rule approach of many-body response theory. In the following we briefly outline the sum-rule approach and present some of the results which are relevant for the present work.

According to the basic result of the sum-rule approach bohigas ; lipparini  the upper bound of the lowest excitation energy is given by

$$\mathrm{}\mathrm{\Omega }_{ex}=\sqrt{\frac{m_3}{m_1}},$$
(8)
where

$$m_p=\underset{n}{}|0|F|n|^2\left(\mathrm{}\omega _{n0}\right)^p,$$
(9)
is the $`p`$-th order moment of the excitation energy associated with the excitation operator $`F`$ and $`\mathrm{\Omega }_{ex}`$ is the frequency of excitation. Here $`\mathrm{}\omega _{n0}=E_nE_0`$ is the excitation energy of eigenstate $`|n`$ of the Hamiltonian $`H`$ of the system. The upper bound given by Eq. (8) is close to the exact lowest excited state when this state is highly collective, that is, when the oscillator strength is almost exhausted by a single mode. This condition is satisfied by the trapped bosons in most of the cases. Moreover, Eq. (8) can be used for computation of the excitation energies by exploiting the fact that the moments $`m_1`$ and $`m_3`$ can be expressed as expectation values of the commutators between $`F`$ and $`H`$ in the ground state $`|0`$ bohigas ; lipparini :

$`m_1`$ $`=`$ $`{\displaystyle \frac{1}{2}}0|[F^{},[H,F]]|0,`$
$`m_3`$ $`=`$ $`{\displaystyle \frac{1}{2}}0|[[F^{},H],\left[[H,[H,F]]\right]|0.`$ (10)
The main advantage of the sum-rule approach is that it allows us to calculate the dynamic properties like excitation frequencies of many-body systems with the knowledge of ground state $`|0`$ (or the ground-state density) only. Thus, with the knowledge of a reasonably accurate ground-state density we calculate the expectation values given in Eq.(10) to determine the frequencies of collective oscillations. It is worth mentioning that sum-rule approach has been extensively used to calculate collective frequencies of the trapped condensate stringari ; kimura .

For the purpose of calculation one needs to choose an appropriate excitation operator $`F`$. In this paper we concentrate only on the compressional collective mode characterized by the $`z`$-component of angular momentum index, $`m^{}=0`$. For the axially symmetric trap the mode $`m^{}=0`$ involves coupling of the monopole ($`l=0`$) and the quadrupole ($`l=2`$) modes. Following Kimura et al. kimura  the operator $`F`$ for the excitation of $`m^{}=0`$ collective mode is written as

$$F=\underset{i}{}\left(x_i^2+y_i^2\alpha z_i^2\right),$$
(11)
where $`\alpha `$ is a parameter characterizing the coupling of the two modes due to the axial symmetry of the trap. This parameter is obtained by making the excitation energy given by Eq.(8) extremal. Here, we note that $`F`$ given above can also be used to describe excitation of the monopole ($`\alpha =1`$) and the quadrupole ($`\alpha =2`$) modes separately. For spherically symmetric trap ($`\lambda _0=1`$) the two modes get decoupled and $`\alpha =1`$ and $`\alpha =2`$ correspond to the minimum and the maximum values of the excitation energies, respectively.

Using the energy functional given by Eq. (1) along with the expansion (2) we find after a tedious although straightforward algebra the following expressions for the moments $`m_1`$ and $`m_3`$:

$$m_1=\frac{4\mathrm{}^2}{m^2}\left(\frac{U_{}}{\left(\omega _{}^0\right)^2}+\alpha ^2\frac{U_z}{\left(\omega _z^0\right)^2}\right),$$
(12)
$$m_3=\frac{8\mathrm{}^4}{m^2}\left[T_{}+U_{}+\alpha ^2\left(T_z+U_z\right)+\left(1\alpha /2\right)^2\left(E_{int}^1+\frac{9}{4}E_{int}^2\right)\right].$$
(13)
In the above equations $`T_{}`$ and $`U_{}`$ ( $`T_z`$ and $`U_z`$ ) are the transverse components ($`z`$-component) of the kinetic and potential energies, respectively. By substituting Eqs. (12) and (13) in Eq. (8) we get following expression for the frequency of $`m^{}=0`$ mode

$$\mathrm{\Omega }_{m^{}=0}^2=2\frac{\left[T_{}+U_{}+\alpha ^2\left(T_z+U_z\right)+\left(1\alpha /2\right)^2\left(E_{int}^1+\frac{9}{4}E_{int}^2\right)\right]}{\left(\frac{U_{}}{\left(\omega _{}^0\right)^2}+\alpha ^2\frac{U_z}{\left(\omega _z^0\right)^2}\right)}.$$
(14)
The frequency of collective oscillations is determined by substituting the values of energy components from the variational calculation and making $`\mathrm{\Omega }_{m^{}=0}^2`$ extremal with respect to $`\alpha `$. The variation with respect to $`\alpha `$ leads to a quadratic equation in $`\alpha `$ and the two roots, corresponding to the maximum ($`\alpha _+`$) and the minimum ($`\alpha _{}`$) values of $`\mathrm{\Omega }_{m^{}=0}^2`$, are

$$\alpha _\pm =\frac{B}{2}\pm \frac{1}{2}\sqrt{B^2+4C},$$
(15)
with

$$B=\frac{\left(2f_2f_4+f_3f_4/22f_1f_52f_3f_5\right)}{f_3f_5},$$

and

$$C=\frac{f_4}{f_5},$$
(16)
where $`f_1=T_{}+U_{}`$, $`f_2=T_z+U_z`$, $`f_3=E_{int}^1+\frac{9}{4}E_{int}^2`$, $`f_4=U_{}/(\omega _{}^0)^2`$ and $`f_5=U_z/(\omega _z^0)^2`$. In the next Section we first apply these results to estimate analytically the correction to the collective frequency due to the LHY term within the TF approximation followed by the results of our variational calculation.

## III Results and discussion

### III.1 Analytical estimate within TF approximation

Before discussing the results of our variational calculation we present the analytical estimate of collective frequencies within the TF approximation. We generalize the perturbative results of Ref. pitaevskii ; braaten2  and also discuss the validity regime of these results. For convenience we consider the case of spherically symmetric trap ($`\lambda _0=1`$). Generalization for the axially symmetric trap is straightforward. It is easy to verify that for spherically symmetric trap $`B=1`$ and $`C=2`$ and the two roots are $`\alpha =2`$ and $`\alpha =1`$ corresponding to the quadrupole and the monopole modes, respectively. It is well known that the frequency of quadrupole mode is not affected by the LHY term pitaevskii . Therefore we focus our attention on the monopole mode. By substituting $`\alpha =1`$ in Eq. (14) the frequency of the monopole mode of bosons in a spherical trap ($`U=U_{}+U_z`$ ) within the TF approximation ($`T_{}=T_z=0`$) can be written as

$$\frac{\mathrm{\Omega }^2}{\left(\omega _{}^0\right)^2}=\left[\frac{U+\frac{9}{4}E_{int}^1+\frac{81}{16}E_{int}^2}{U}\right].$$
(17)
Hereafter we drop for convenience the subscript $`m^{}=0`$ in $`\mathrm{\Omega }`$. To eliminate $`U`$ from the above equation we make use of the virial relation banerjee

$$2T2U+3E_{int}^1+\frac{9}{2}E_{int}^2=0,$$
(18)
and obtain within the TF approximation

$$\frac{\mathrm{\Omega }^2}{\left(\omega _{}^0\right)^2}=\left[\frac{\frac{15}{4}E_{int}^1+\frac{117}{16}E_{int}^2}{\frac{3}{2}E_{int}^1+\frac{9}{4}E_{int}^2}\right].$$
(19)
We note that the expression given by Eq. (19) is exact within the TF limit. To derive an analytical expression for the frequency we express $`E_{int}^1`$ and $`E_{int}^2`$ in terms of parameter $`x=\sqrt{n(0)a^3}`$. This is done by using the analytical expressions for the total energy and chemical potential dalfovo ; banerjee  obtained within the TF approximation with the LHY term included in the expansion $`ϵ(n)`$. These expressions are:

$`E_{int}^1`$ $`=`$ $`{\displaystyle \frac{2}{7}}\mu _{TF}{\displaystyle \frac{9}{16}}\mu _{TF}x,`$
$`E_{int}^2`$ $`=`$ $`{\displaystyle \frac{5}{8}}\mu _{TF}x,`$ (20)
where $`\mu _{TF}`$ is the chemical potential of trapped bosons obtained within the TF approximation dalfovo . It is important to note that the expressions for $`E_{int}^1`$ and $`E_{int}^2`$ given above (Eq. (20)) are approximate in nature and neglect terms of order $`𝒪(na^3)`$ and, thus, expected to be valid in the small gas parameter regime. It has been recently shown fabrocini  that the expressions in Eq. (20) are not valid for very large gas parameters and lead to results which differ considerably from the exact results . In the large-gas-parameter regime one needs to solve the nonlinear equation of the density more accurately with the help of numerical methods. Consequently, to compute the frequencies of collective oscillation (for that matter any observable) in the large-gas-parameter regime within the TF approximation one needs a non-iterative or non-perturbative method to solve the nonlinear equation. In this paper we take a different route to achieve this goal. Instead of solving the time-dependent differential equation for the density (or the time-dependent GP equation) we exploit the sum-rule approach of the many-body response theory along with the accurate ground-state density determined variationally to compute the frequencies of collective oscillations. We describe the results of such a calculation in the next subsection.

Now by substituting Eq. (20) in Eq. (19) we get following expression for the frequency of the monopole collective mode within the BGP theory:

$$\frac{\mathrm{\Omega }}{\omega _0}=\left[\frac{5\left(1+\frac{147}{64}x\right)}{1+\frac{21}{16}x}\right]^{1/2},$$
(21)
where $`\omega _0`$ is the angular frequency of the spherically symmetric trap. In the limit of $`x0`$ above expression (Eq.(21)) correctly reduces to the mean-field result $`\mathrm{\Omega }_{mf}=\sqrt{5}\omega _0`$ within the TF approximation stringari .

It is customary to define the fractional change in the frequency as pitaevskii ; braaten2 ,

$$\delta =\frac{1}{2}\left(\frac{\mathrm{\Omega }^2\mathrm{\Omega }_{mf}^2}{\mathrm{\Omega }_{mf}^2}\right).$$
(22)
By using Eq. (21) we get

$$\delta =\frac{\frac{63}{128}x}{1+\frac{21}{16}x}.$$
(23)
In the limit of $`x<<1`$ we recover the results of Refs. pitaevskii ; braaten2 . Fig. 1 shows the fractional change $`\delta `$ as a function of dimensionless parameter $`a/a_0`$ (where $`a_0`$ is the Bohr radius of hydrogen atom) for $`N=10^4`$ 85Rb atoms trapped in a spherical trap with $`\omega _0/2\pi =\left(\left(\omega _{}^0\right)^2\omega _z^0\right)^{1/3}/2\pi =12.83Hz`$, where $`\omega _{}^0=2\pi \times 17.5`$Hz and $`\omega _z^0=2\pi \times 6.9`$Hz cornish  . The values of the parameter $`a/a_0`$ are chosen such that the condition for the TF approximation, that is , $`Na/a_{h0}>>1`$ is satisfied. Moreover, as discussed above Eq.(20) (and hence Eq. (21)) is valid only for small gas parameter values and accordingly we have restricted ourselves to $`a/a_0=5000`$. It is evident from Fig.1 that the correction to the monopole frequency due to inclusion of LHY term in the interatomic correlation energy obtained by us is slightly lower than the corresponding results of Refs. pitaevskii ; braaten2 . In particular for $`a/a_0=5000`$ we obtain $`7\%`$ change in the monopole frequency as compared to the correction of $`9\%`$ predicted in Ref. pitaevskii ; braaten2 . It is only for very small values of $`x`$ (less than $`10^3`$) that fractional corrections obtained in Ref. pitaevskii ; braaten2  match with Eq. (23).

### III.2 Variational results

Before discussing the results for the excitation frequencies obtained by variational method, we first demonstrate that this method is capable of yielding sufficiently accurate results for the ground-state properties of BEC even for large gas parameter values, as achieved in the experiment of Cornish et al. cornish . The parameters used in this calculation are similar to the ones used in the previous section. We present these results in Table I and II for the spherically symmetric and the axially symmetric traps respectively. To check the accuracy, we compare our variational results with those of Ref. fabrocini  (presented in parenthesis in Table I and II) obtained by numerically solving the modified GP equation. The results in Table I and II clearly show that the variational method gives quite accurate results for the ground-state properties as they compare very well with the numbers of Ref. fabrocini . Therefore the variational method is indeed suitable for describing the ground-state properties even in the regime of large gas parameters. The accuracy of our variational results was further checked by verifying that the generalized virial relation as given in Eq. (18) is satisfied. In all the variational calculations the virial relation was satisfied up to the fifth decimal place or better. This indicates that not only the total energy but its different components are also obtained accurately by the variational method, which is quite crucial for the calculation of the frequencies of collective oscillations.

Now we focus our attention on the results for the frequencies of collective oscillations of trapped bosons. Figs. 2 and 3, show the upper ($`\mathrm{\Omega }_u`$) and the lower ($`\mathrm{\Omega }_l`$) frequencies respectively, associated with the breathing mode as a function of the dimensionless parameter $`a/a_0`$. The results presented in Figs. 2 and 3 correspond to the axially symmetric trap which is more relevant from the experimental point of view. Notice that the maximum value of the parameter $`a/a_0`$ is chosen to be $`10^4`$ in accordance with the experimental number cornish  and the corresponding value of the peak gas parameter is $`10^2`$. For comparison we also show in these figures the corresponding results obtained within the GP theory. It is important to note here that the GP results shown in Figs. 2 and 3 are also obtained variationally by dropping the LHY term from the expansion (2) and thus, retaining only the first term. The results corresponding to the GP case obtained by us match quite well with those already reported in the literature stringari ; kimura ; edwards ; singh ; esry ; sinha ; dalfovoex ; you ; hutchinson . It can be clearly seen from Fig. 2 that for small values of the gas parameter the numbers obtained within the GP and the BGP theory are quite close. The difference between these two results, however, increases with $`a`$. For example, for $`a/a_0=8000`$ the BGP number for $`\mathrm{\Omega }_u`$ is $`6\%`$ higher than the corresponding mean-field GP result. Thus we conclude that the introduction of LHY term in the interatomic correlation energy results in a significant correction to the upper frequency $`\mathrm{\Omega }_u`$ of $`m^{}=0`$ mode of the collective oscillations. Considering the accuracy of current experiments stamperkurn  in measuring the frequencies of the collective excitations, it should be possible to observe these corrections experimentally. It is also important to note that the frequency $`\mathrm{\Omega }_u`$ obtained within the GP theory remains almost constant with the increase in the parameter $`a/a_0`$. In contrast, the BGP theory predicts $`\mathrm{\Omega }_u`$ growing monotonically with $`a/a_0`$. Therefore a signature of the predictions of the BGP theory can be observed by measuring $`a/a_0`$ dependence of the excitation frequency $`\mathrm{\Omega }_u`$.

Next we focus our attention on the lower frequency $`\mathrm{\Omega }_l`$ of the breathing mode which is shown in Fig. 3. For small values of the gas parameter the frequency of the lower branch obtained via GP theory is very close to the corresponding BGP theory result. In the same regime of the gas parameter the frequency $`\mathrm{\Omega }_l`$ decreases with the gas parameter both for the GP and the BGP cases. The decrease of the lower frequency with an increase in the gas parameter value has already been theoretically predicted stringari ; edwards ; esry ; you  and experimentally verified as well jin . However, in contrast to the GP case, frequency of the lower branch for the BGP case starts increasing, albeit slowly, after a critical value of $`a/a_02200`$. It is also noteworthy that the change in the frequency $`\mathrm{\Omega }_l`$ obtained from the BGP theory is considerably lower than that in the upper frequency $`\mathrm{\Omega }_u`$. For example, for $`N=10^4`$ and $`a/a_0=8000`$ the change in the lower frequency is just $`0.8\%`$ while that in the upper frequency is $`6\%`$.

We conclude this section by making a comparison of the frequencies obtained by using the analytical expression as give in Eq. (21) with the results obtained by the variational approach. The comparison is shown in Fig. 4 where we plot the frequency of the monopole mode of a condensate in a spherically symmetric trap as a function of $`a/a_0`$ obtained using Eq. (21) (dashed line) and using the variational method (solid line). The range of values of $`a/a_0`$ is chosen such that the TF approximation (i.e. $`Na/a_{h0}>>1`$ ) is satisfied. Accordingly the minimum value of $`a/a_0`$ in Fig. 4 is chosen to be $`10^3`$ which corresponds to $`Na/a_{h0}=174`$. It can be seen that the results obtained using Eq.(21) are systematically higher than the corresponding results from variational calculation and the difference between the two grows as the value the gas parameter is increased. Moreover, we also note that the two results are reasonably close for values of $`a/a_0`$ of the order of 3500 ($`n(0)a^33.77\times 10^3`$). Therefore, we emphasize that the correction obtained within the TF approximation as given by Eq. (21) is valid as long as the gas parameter is of the order of $`10^3`$. As discussed earlier, for condensates with the values of the gas parameter beyond $`10^3`$ one needs to evaluate $`E_{int}^1`$ and $`E_{int}^2`$ in a non-perturbative manner. The variational approach employed in this paper accomplishes exactly this task.

## IV Conclusion

In this paper we have calculated the frequencies of collective excitation of a trapped condensates in the large-gas-parameter regime using the BGP theory . To go beyond the GP theory we have included a higher-order term (LHY term) in the interatomic correlation energy obtained from the ground-state energy of a uniform Bose gas. The calculations of the frequencies were performed by employing the sum-rule approach of many-body response theory along with the variational method to obtain the ground-state energy components. We have also derived an analytical expression for the corrections to the frequencies of the compressional mode valid for small gas parameters within the TF approximation. The results of our variational calculations clearly show that for large gas parameter values, the frequencies of the compressional mode characterized by the index $`m^{}=0`$ are significantly altered. Although the change in the lower branch of $`m^{}=0`$ mode, $`\mathrm{\Omega }_l`$, is small, we observe a new qualitative feature. The frequency $`\mathrm{\Omega }_l`$ increases slowly as the gas parameter exceeds a critical value in contrast to a monotonically decrease predicted within the GP theory. It would be very interesting to verify the predictions of the BGP theory by measuring collective excitation frequencies for condensates with large gas parameter as reported in Ref. cornish .

## Acknowledgement

It is a pleasure to thank Dr. S. C. Mehendale for helpful discussions and critical reading of the manuscript.


text:
# 1 Iron deficiency in PG1159 stars

## 1 Iron deficiency in PG1159 stars

The origin of hot hydrogen-deficient post-AGB stars (spectral types \[WC\] and PG1159) is supposed to be a late He-shell flash. Detailed summaries on their spectroscopic characteristics and quantitative analyses, and relevant evolutionary calculations can be found e.g. in Werner (2001) and Herwig (2001), respectively.

In the last white dwarf workshop, held in Delaware two years ago, we presented for the first time clear evidence for iron deficiency in a PG1159 star (Miksa et al. 2001). Because PG1159 central stars are very hot ($`T_{\mathrm{eff}}`$ $`>`$75 000 K), the metals are highly ionized. The dominant ionization stage of iron in the line formation region is Fe VII and most of its lines are located in the FUV spectral region. Contrary to our expectation we were not able to detect any Fe VII line in the FUSE FUV spectrum of the pulsating central star of K1-16, which means that iron is deficient by at least one dex. This has initiated an effort to utilize archival IUE and HST spectra as well as new FUSE spectra to look for iron in other PG1159 stars (Miksa et al. 2002). It turned out that IUE and HST spectra do show a hint of iron deficiency in some other objects, but the data quality is not sufficient for reliable quantitative analyses. The same holds for a number of FUSE spectra from PG1159 stars. However, in two more cases, where FUSE spectra with sufficiently high S/N were obtained, we were successful in proving that iron is deficient. The first one is the central star of NGC 7094, which is classified as a hybrid-PG1159 star, because it exhibits H Balmer lines. No iron lines were detected and we concluded that iron is depleted by at least 1.5 dex (Miksa et al. 2002). The next case is the central star of Abell 78, which is classified as one of the rare \[WC\]–PG1159 transition objects. For Abell 78 we find an iron deficiency of 1.5 dex from the lack of iron lines as well (Werner et al. 2002a).

What is the origin of the iron deficiency? We think that iron has been transformed to heavier elements. The high C and O abundances in PG1159 stars result from envelope mixing caused by a late He-shell flash. This event also modifies the near-solar abundance ratios of iron-peak elements in the envelope by dredging up matter in which s-process elements were built-up by n-capture on 56Fe seeds during the AGB phase. This scenario can and will be tested by analyzing the resulting Fe/Ni abundance ratio, because it is significantly changed in the intershell region in favor of Ni by the conversion of 56Fe into 60Ni. More quantitative results from nucleosynthesis calculations in appropriate stellar models have been presented recently (Herwig et al. 2002) and inclusion of nuclear networks in evolutionary model sequences will become available in the near future.

Interestingly, Asplund et al. (1999) have found that in Sakurai’s object, which is thought to undergo a late He-shell flash, iron is reduced to 0.1 solar and Fe/Ni$``$3. This and other s-process signatures should also be exhibited by Wolf-Rayet type central stars and PG1159 stars. And in fact, latest results confirm that iron deficiency among H-deficient post-AGB stars is not restricted to PG1159 stars. It appears that three Wolf-Rayet type central stars are iron deficient, too. Gräfener et al. (2002) report a low iron abundance in SMP 61, an early type \[WC5\] central star in the LMC. Its abundance is at least 0.7 dex below the LMC metallicity. Crowther et al. (2002) find evidence for an iron underabundance of 0.3–0.7 dex in the Galactic \[WC\] type central stars stars NGC 40 (\[WC8\]) and BD+303639 (\[WC9\]). In this context it is also interesting to remark that the hot DO white dwarf RE 0503-289 is iron deficient and rich in nickel (Barstow et al. 2000). If this is related to the iron deficiency in PG1159 stars is an open question.

## 2 Chandra spectroscopy of H1504+65

H1504+65 is the most extreme PG1159 star, not only because of its extraordinarily high effective temperature, which is close to 200 000 K. We have shown that H1504+65 is not only hydrogen-deficient but also helium-deficient. From optical spectra it was concluded that the atmosphere is primarily composed of carbon and oxygen, by equal amounts. Neon lines were detected in soft X-ray spectra taken with the EUVE satellite and in an optical-UV Keck spectrum, and a neon abundance of 2–5% (mass fraction) was derived (Werner & Wolff 1999).

H1504+65 was observed with Chandra on Sep. 27, 2000, with an integration time of 7 hours. Flux was detected in the range 60Å–160Å and the spectral resolution is about 0.1Å. The spectrum displays astonishing details (Fig. 1). At first sight the spectrum appears to be rather noisy, but comparison with model spectra suggests that most absorption features in the observation are probably real. Mg lines are discovered for the first time in H1504+65 (see Werner et al. (2002b) for more details). Exploratory models including opacities of heavy metals (Ca–Ni) were computed. This introduces a very large number of additional lines, however, identifying individual heavy metal lines will be problematic or even impossible, because the vast majority of lines in the available lists have uncertain wavelength positions.

The Mg abundance is of the order 2%, similar to the Ne abundance determined earlier. This allows two alternative interpretations. Either we see a naked C-O white dwarf with 3$`\alpha `$ processed matter, which is, according to evolutionary models, dominated by C and O with an admixture of 22Ne and 25Mg of the order 1%. Or we see a naked O-Ne-Mg white dwarf which has undergone carbon burning. A comparison with the respective evolutionary models (Ritossa et al. 1999) suggests that we could look onto such a naked remnant where 20Ne and a mixture of several Mg isotopes (24Mg,25Mg,26Mg) should be present. We have already speculated earlier, that the onset of C burning and subsequent evolutionary processes might explain the uniqueness of H1504+65. In order to confirm this idea, we will also look for Na absorption lines, because the stellar models predict a relatively high amount of 23Na in this case.

To conclude, H1504+65 might have been one of the “heavyweight” intermediate-mass stars (9 M$`{}_{}{}^{}M`$ 11 M) which form white dwarfs with electron-degenerate O-Ne-Mg cores resulting from carbon burning. This speculation is corroborated by the fact that H1504+65 is one of the most massive PG1159 stars known (0.86$`\pm `$0.15 M).