text:
# The influence of density stratification and multiple nonlinearities on solar torsional oscillations

## 1 Introduction

Recent analyses of the helioseismic data have produced a wealth of information concerning the dynamical modes of behaviour associated with the differential rotation in the solar convection zone (CZ). Previous analyses of data from the Michelson Doppler Imager (MDI) instrument on board the SOHO spacecraft (Howe et al. 2000a) and the Global Oscillation Network Group (GONG) project (Antia & Basu 2000) have found that torsional oscillations penetrate into the convection zone, to depths of at least 10 percent in radius, but did not agree on the behaviour in the poorly resolved region near the bottom of the convection zone. In order to understand these observed dynamical phenomena theoretically, and partly motivated by these apparent discrepancies, we have previously made detailed studies of the observed variations in the differential rotation in the CZ in the framework of nonlinear dynamo models which included the nonlinear action of the azimuthal component of the Lorentz force of the dynamo generated magnetic field on the solar angular velocity (Covas et al. 2000a,b; Covas et al. 2001a,b; Tavakol et al. 2002). These studies produced a number of important robust features, that may be treated as their predictions. These include: (i) the robust existence of torsional oscillations that penetrated all the way down to the bottom of the convection zone for near-critical and moderately supercritical dynamo regimes; (ii) the existence of parameter ranges for which the supercritical models show spatiotemporal fragmentation, with different oscillatory modes of behaviour at the top and the bottom of the CZ; (iii) the existence of equatorial as well as polar branches in the angular velocity residuals; and (iv) the presence of secondary modulations on these equatorial branches. Predictions (i) and (ii) are in principle able to account for different aspects of earlier observations. Interestingly, the most recent inversions by Vorontsov et al. (2002a) demonstrate further important agreement with our predictions. In particular they provide strong evidence that these oscillations do in fact penetrate down to the bottom of the CZ. Furthermore they find clear evidence for the presence of polar branches as well as of modulations on the equatorial branch (see Vorontsov et al. 2002b for preliminary results, and also Vorontsov et al. 2003).

Despite these successes, there are two important issues concerning these models that require further study. Firstly, they ignore an important feature of the real Sun, namely the presence of a strong density stratification in the solar convection zone and, secondly, they involve only one form of nonlinearity – that due to the Lorentz force acting on the large-scale velocity field (here just the rotational flow). A more realistic model would be expected to include other forms of nonlinearity, such as the effects of the Lorentz force on the small-scale motions, for example as naively modelled by a nonlinear $`\alpha `$–quenching. Investigating the consequences of these shortcomings is particularly important, given the limitations on the length of the helioseismic data sets and the resulting uncertainties in the inversions, especially regarding the behaviour in the lower parts of the CZ (see e.g. Vorontsov et al 2002).

The aim of this paper is to improve our previous model by (i) including a density stratification in our previous model, and (ii) by adding a second form of nonlinearity in the form of an explicit $`\alpha `$–quenching. In contrast, in our previous work we imposed a radial variation of alpha that we felt might be consistent with a strong quenching in the overshoot region where the magnetic field was expected to be strong (by putting $`\alpha =0`$ in this region). We can now make a detailed study of both the effects of density stratification and the competing forms of nonlinearity on various observationally important features of our model. Throughout the underlying zero order angular velocity is taken to be consistent with the recent inversions of the helioseismic (MDI) data.

The structure of the paper is as follows. In the next section we outline our model, and Sects. 3 and 4 contain our detailed results. In Sect. 5 we attempt to draw some conclusions.

## 2 The model

### 2.1 Equations

We write the large-scale velocity field in the form $`u=v\widehat{\varphi }\frac{1}{2}\eta `$, and put $`v=v_0+v^{}`$, $`v_0=\mathrm{\Omega }_0r\mathrm{sin}\theta `$ being the underlying (assumed known) solar rotation law. Thus we solve the standard mean-field dynamo equation

$$\frac{B}{t}=\times (u\times B+\alpha B\eta \times B),$$
(1)
for the mean magnetic field and the corresponding equation

$$\frac{v^{}}{t}=\frac{(\times B)\times B}{4\pi \rho r\mathrm{sin}\theta }.\widehat{\varphi }+\rho ^1.\tau ,$$
(2)
($`\tau `$ the stress tensor) for the deviations from the underlying ($``$ mean) rotational velocity $`v_0`$; as usual $`\eta `$ and $`\nu `$ are turbulent transport coefficients. We proceed as described in Moss & Brooke (2000) and in previous papers in this series, Nondimensionalizing in terms of the solar radius $`R`$ and a time $`R^2/\eta _0`$, where $`\eta _0`$ is the maximum value of $`\eta `$, and putting $`\mathrm{\Omega }=\mathrm{\Omega }^{}\stackrel{~}{\mathrm{\Omega }}`$, $`\alpha =\alpha _0\stackrel{~}{\alpha }`$, $`\eta =\eta _0\stackrel{~}{\eta }`$, $`B=B_0\stackrel{~}{B}`$ and $`v^{}=\mathrm{\Omega }^{}R\stackrel{~}{v}^{}`$. This gives dynamo parameters $`R_\alpha =\alpha _0R/\eta _0`$, $`R_\omega =\mathrm{\Omega }^{}R^2/\eta _0`$, $`P_\mathrm{r}=\nu _0/\eta _0`$, where $`\mathrm{\Omega }^{}`$ is the solar surface equatorial angular velocity and $`B_0=(\mathrm{\Lambda }4\pi \rho _0\eta _0\mathrm{\Omega }^{})^{1/2}`$, where $`\mathrm{\Lambda }`$ is an arbitrary scaling factor and $`\rho _0`$ is the density at the lower boundary of the model. As discussed in Moss & Brooke (2000), the physical field $`B_0\stackrel{~}{B}`$ is independent of the value of $`\mathrm{\Lambda }`$ – for ‘historical’ reasons we set $`\mathrm{\Lambda }=10^4`$. $`P_\mathrm{r}`$ is the turbulent Prandtl number.

The model includes a density stratification discussed in Section 2.2 below, as well as a second source of nonlinearity in the form of an $`\alpha `$–quenching discussed in Section 2.3 below.

In this investigation, $`\mathrm{\Omega }_0`$ is given in $`0.64r1`$ by an interpolation on the MDI data obtained from 1996 to 1999 (Howe et al. 2000a). For $`\alpha `$ we put $`\alpha =\alpha _\mathrm{r}(r)f(\theta )`$, where $`f(\theta )`$ was chosen to be $`\mathrm{sin}^2\theta \mathrm{cos}\theta `$. This form has been used previously (see e.g. Rüdiger & Brandenburg 1995) and its choice here is simply to make the butterfly diagrams more realistic. For $`\alpha _\mathrm{r}`$ we took the following choices. In the absence of $`\alpha `$–quenching, $`\alpha _\mathrm{r}`$ was set to be $`\alpha _\mathrm{r}=1`$ for $`0.7r0.8`$ with cubic interpolation to zero at $`r=r_0`$ and $`r=1`$. In presence of $`\alpha `$–quenching, we took $`\alpha _\mathrm{r}=1`$ throughout the convection zone. We use the convention that $`\alpha _\mathrm{r}>0`$ and $`R_\alpha <0`$, to get the correct sense of migration of field patterns. Also, in order to take into token account the likely decrease in the turbulent diffusion coefficient $`\eta `$ in the overshoot region, we allowed a simple linear decrease from $`\stackrel{~}{\eta }=1`$ at $`r=0.8`$ to $`\stackrel{~}{\eta }=0.5`$ in $`r<0.7`$.

At the surface $`r=R`$ we took vacuum boundary conditions, whereby the poloidal field within $`r=R`$ is smoothly joined, by a matrix multiplication, to an external vacuum solution; the azimuthal field there satisfies $`B=0`$. At the inner boundary $`r=r_0`$ we used the same conditions as Tavakol et al. (2002). We note that ‘reasonable’ changes in these boundary conditions do not qualitatively change our conclusions (cf. Tavakol et al. 2002).

We set $`r_0=0.64`$; with the solar CZ proper being thought to occupy the region $`r>\text{ }\mathrm{\hspace{0.33em}0.7}`$, the region $`r_0r<\text{ }\mathrm{\hspace{0.33em}0.7}`$ can be thought of as an overshoot region/tachocline, where the major radial gradients of $`\mathrm{\Omega }`$ occur. Equations (1) and (2) were solved using a modified version of the code described in Moss & Brooke (2000), using the the above boundary conditions, over the range $`r_0r1`$, $`0\theta \pi `$. We used a mesh resolution of $`61\times 101`$ points, uniformly distributed in radius and latitude respectively, i.e. $`\mathrm{\Delta }r=0.06,\mathrm{\Delta }\theta =0.01\pi `$.

### 2.2 Density stratification

The density profile in the solar model computed by Christensen-Dalsgaard et al. (1996; henceforth referred to as CD96) is shown in Fig. 1. We approximated the density stratification in the solar convection zone by a polytropic density profile, writing in general

$$\rho (r)/\rho _0=\left(\frac{r_sr}{r_sr_0}\right)^n,$$
(3)
where

$$r_s=\left(\frac{\rho _{\mathrm{ratio}}^{1/n}r_0}{\rho _{\mathrm{ratio}}^{1/n}1}\right);$$
(4)
we take $`n=1.5`$ below. This density profile corresponds to a polytrope truncated at the appropriate radius to give the required value of $`\rho _{\mathrm{ratio}}`$, with this radius subsequently renormalized to unity. $`\rho _{\mathrm{ratio}}`$ is thus the ratio of the densities at the top ($`r=1`$) and the bottom ($`r=r_0`$) of the computational domain; $`\rho _0=\rho (r_0)`$. Several examples of such profiles, with different values of $`\rho _{\mathrm{ratio}}`$ are also shown in Fig. 1. For $`\rho _{\mathrm{ratio}}100`$, these give a reasonable approximation to the density stratification in the bulk of the solar CZ, apart from the outermost regions near the solar surface

### 2.3 Explicit $`\alpha `$–quenching

The nonlinearity in our original model is through the nonlinearity in the Lorentz force and the subsequent presence of the $`v^{}`$ term in the induction equation. In all of our previous studies the form of $`\alpha _\mathrm{r}`$ was chosen to be independent of the magnetic field $`B`$. It did, however, include an implicit $`\alpha `$-quenching in the overshoot layer in that $`\alpha `$ was prescribed to vanish in $`r0.7`$ and to be relatively small in the layers immediately above where $`|B|`$ was anticipated to be relatively large. Here we study the effects of having an additional nonlinearity in the form of an explicit nonlinear $`\alpha `$-quenching given (reverting temporarily to dimensional units) by

$$\alpha =\alpha _\mathrm{r}(r)\frac{\mathrm{sin}^2\theta \mathrm{cos}\theta }{1+g^{}|B^2|/B_{\mathrm{eq}}^2}$$
(5)
where $`B_{\mathrm{eq}}=(4\pi \rho v_t^2)^{1/2}`$ is the equipartition field strength, $`v_t`$ is the turbulent velocity and $`g^{}`$ is of order unity. (The trigonometric factor is chosen to agree with that used in our earlier work – see above.) In order to estimate how $`B_{\mathrm{eq}}^2`$, appearing in the denominator, varies with depth we recall that to a very good approximation the luminosity $`L=4\pi r^2F_rr^2\rho v_{t}^{}{}_{}{}^{3}`$ is a constant in the CZ, which gives $`v_t\rho ^{1/3}r^{2/3}`$ and hence $`\rho v_{t}^{}{}_{}{}^{2}\rho ^{1/3}r^{4/3}`$. This suggests that the denominator in (5) will be a slowly varying function of the depth, especially away from the solar surface, and this expectation is supported by detailed mixing length models of the solar CZ. Thus we take $`B_{\mathrm{eq}}^2\rho v_t^2=\mathrm{const}`$, and so the $`\alpha `$-quenching expression (5) can be well approximated (now temporarily using dimensionless variables again) by

$$\stackrel{~}{\alpha }=\alpha _\mathrm{r}(r)\frac{\mathrm{sin}^2\theta \mathrm{cos}\theta }{1+g|B^2|},$$
(6)
where $`g=g^{}\mathrm{\Lambda }4\pi \rho _0\eta _0\mathrm{\Omega }^{}/B_{\mathrm{eq}}^2`$. $`B_{\mathrm{eq}}^2`$ is constant, and thus so is $`g`$. (In this connection it may be noted that there is an ongoing controversy regarding the nature and strength of $`\alpha `$-quenching. This is not the place to rehearse the arguments for and against ‘strong’ alpha-quenching (where $`g^{}1`$); the issue is unresolved, and we use (6) as a commonly adopted nonlinearity.)

We can make an estimate of what might be plausible values for $`g`$ by taking $`\eta _0=3\times 10^{11}`$ cm2sec-1n (see Sect. 3 below), $`\mathrm{\Lambda }=10^4`$, $`\rho _0=1`$ gm cm-3, $`\mathrm{\Omega }^{}=2\times 10^6`$ sec-1, $`B_{\mathrm{eq}}=3500`$ G (the latter estimated from an envelope model by Baker & Temesvary (1966) (more recent models do not give significantly different values). Then $`g3\times 10^3g^{}`$. (As noted above, $`B_0^2,g\mathrm{\Lambda }`$, but $`\stackrel{~}{B}\mathrm{\Lambda }^{1/2}`$, so the results in physical variables do not depend on the value adopted for $`\mathrm{\Lambda }`$.)

The inclusion of both types of nonlinearity in the model allows us to study the relative importance of each term and in particular what happens to the amplitudes of the torsional oscillations as $`g`$ increases.

## 3 Results: effects of density stratification

Using the above model, we studied the dynamics in the convection zone, and in particular the changes caused by the introduction of the density stratification and the additional nonlinearity. We concentrated on the effects of these changes on the torsional oscillations and in particular on their amplitudes (i.e. the angular velocity residuals, defined as the current angular velocity minus the background angular velocity $`\mathrm{\Omega }_0`$), since these are the most important indicators to be compared with observations. The model predictions are likely to be of particular importance in deeper regions where the uncertainties in the inversions are largest.

We first study the effects of density stratification in isolation, with $`g=0`$. We fixed $`R_\omega =6\times 10^4`$ (corresponding to $`\eta _03\times 10^{11}`$ cm2sec-1) by requiring that the dynamo period is near to 22 years close to marginal stability. The marginal value of $`R_\alpha `$ is then about $`3.2`$.

### 3.1 Comparison with previous results

To begin with we made a detailed comparison of the results obtained here in the presence of density stratification, with our previous studies which ignored this effect. In particular, we concentrated on one of the findings of these earlier studies, namely the fact that torsional oscillations were found to extend all the way down to the bottom of the CZ.

Our results here show that, for both near-critical, and moderately supercritical dynamo regimes, the oscillations extend all the way down to the bottom of the convection zone, for all $`1\rho _{\mathrm{ratio}}10^6`$. This is in line with our previous results which shows that such behaviour is an extremely robust feature of the models considered to density stratification. An example of such generic penetration is depicted in Fig. 2. This is important observationally since the inversions have so far not been able to resolve unambiguously regimes near the bottom of the convection zone (Howe et al. 2000a; Antia & Basu 2000).

For completeness we also show in Fig. 3 the butterfly diagram (i.e. latitude-time plot) for the toroidal component of the magnetic field $`B`$ near the surface, as well as the corresponding plot for the angular velocity residuals (i.e. the near-surface torsional oscillations) in Fig. 4. As can be seen, the presence of significant density stratification does not alter the qualitative behaviour of the magnetic field and torsional oscillations.

### 3.2 Amplitudes of oscillations as a function of $`\rho _{\mathrm{ratio}}`$

From an observational point of view, an important feature of the dynamics of the CZ is the way the amplitudes of the torsional oscillations vary as a function of model ingredients and parameters, as well as depth in the CZ. A theoretical understanding is crucially important here, especially given the uncertainties in inversions of the helioseismic data in the lower parts of the CZ.

With this in mind, we studied how the amplitudes of the torsional oscillations change as a function of $`\rho _{\mathrm{ratio}}`$, for two values of $`R_\alpha `$; one corresponding to a near-critical and the other to a supercritical dynamo regime. The results are depicted in Figs. 5 and 6. For orientation, we recall that the corresponding amplitudes in the sun, deduced from helioseismological observations, are of the order of one nHz.

The important trend revealed by these figures is that, in both the near–critical and supercritical cases, the amplitudes of the torsional oscillations in the upper part of the CZ increase as $`\rho _{\mathrm{ratio}}`$ increases, while those at the bottom decrease. The decrease in amplitude of the torsional oscillations in the bottom of the convection zone, as $`\rho _{\mathrm{ratio}}`$ increases (i.e. as the density profile becomes more realistic) is of practical importance, as it suggests that the detection of such oscillations may be even more difficult than might appear at present just from the limits to the resolution of the helioseismological inversions at the bottom of the convection zone. The increase of the amplitudes at the top is also important since it improves the agreement between simulated and observed values.

Also of potential interest is the way that the distribution of the perturbation kinetic energy, defined by $`_V\rho (2v_0v^{}+v^2)𝑑V`$, is modified as the density stratification changes. In Fig. 7 we show the variation of this integral with radius, where we take the domain $`V`$ of integration to be the local spherical shell with thickness $`\mathrm{\Delta }r`$. Introducing a stratification markedly reduces the perturbation kinetic energy in the near-surface regions, but the differences between the two stratified models illustrated are relatively small. However the low values of the near-surface density in the stratified models means that the amplitude of the torsional oscillations in this region are not correspondingly decreased – see Figs. 5 and 6. This relative insensitivity of the observed quantity to the stratification is reassuring.

### 3.3 Amplitudes of oscillations as a function of $`R_\alpha `$

To begin with, we verified that for a given density profile (in presence of density stratification), the amplitudes of the oscillations increase as the dynamo number $`R_\alpha `$ is increased (see also Covas et al 2001a). The results are shown in Fig. 8. We then studied how the amplitudes change as we change the density profile by increasing the value of $`\rho _{\mathrm{ratio}}`$. For comparison we also calculated the amplitudes for the density profile of CD96. The results are shown in Fig. 9. As can be seen, the effect of increasing $`\rho _{\mathrm{ratio}}`$ is to increase the amplitudes of oscillations near the surface. The largest amplitudes are attained for the CD96 density distribution.

### 3.4 Amplitudes of oscillations as a function of depth and $`\rho _{\mathrm{ratio}}`$

We also studied the variations of the amplitudes of oscillations as a function of depth for a number of values of $`\rho _{\mathrm{ratio}}`$, for both a near-onset regime (with $`R_\alpha =3.2`$) as well as for more supercritical parameters ($`R_\alpha =4.4`$), and the results are shown in Figs. 10 and 11 respectively. As can be seen, despite the latter having larger amplitudes, the two cases have important similarities. In each the amplitudes near the surface layers increase with increasing $`\rho _{\mathrm{ratio}}`$, whereas they decrease near the bottom.

## 4 Results: comparative study of the nonlinearities

In order to make a comparative study of the effects of the two competing forms of nonlinearities considered here, we begin (see also below) by recalling that, in the absence of any alpha-quenching, increasing the degree of nonlinearity due to the Lorentz force (i.e. increasing $`R_\alpha `$) results in an increase in the amplitudes of the torsional oscillations. The question is what is the effect of increasing the strength of the additional nonlinearity due to $`\alpha `$–quenching on these amplitudes?

We studied this question by making a detailed study of how the amplitudes of torsional oscillations change as a function of $`g`$, the coefficient of the $`\alpha `$–quenching term (Eq. (5)). The results are presented in Fig. 12. A number of trends are apparent. The amplitudes of the torsional oscillations decrease at all depths with increasing $`\alpha `$-quenching coefficent $`g`$. However, the oscillations at the top of the CZ have the highest amplitudes throughout.

To study the influence of increasing the $`\alpha `$–quenching coefficient, we recall that in the case of pure alpha-quenching,

$$\frac{1}{g},$$
(7)
where the angled parentheses denote some sort of global average. We checked that this was exactly satisfied for our code in the quasi-kinematic approximation when the Navier-Stokes equation was switched off. Interestingly, we also found similar behaviour when the additional nonlinear feedback via the Navier-Stokes equation was present. This is shown in Fig. 13, which displays the behaviour of the global magnetic energy as a function of $`g`$.

The amplitudes of the torsional oscillations increase somewhat with $`|R_\alpha |`$ but in the range accessible to our code ($`|R_\alpha |15`$) this increase is insufficient to compensate for the decrease caused by the increase in the $`\alpha `$-quenching coefficient $`g`$ (given by (7)).

We can add here that in our previous, very preliminary, study of the effects of alpha-quenching on the torsional oscillations (Covas et al 2001a), we were only able to find oscillatory solutions for very small values of the alpha-quenching coefficient ($`g<0.1`$). What has enabled us to find the existence of oscillatory solutions at such large values of $`g`$ here is our recent discovery of the existence of multiple attractors in the dynamo system, including both oscillatory and steady solutions, that coexist to $`g1000`$, i.e. in the physically plausible range as estimated in Sect. 2.3.

## 5 Discussion

We have made a detailed study of the effects of including density stratification in the model, as well as an additional form of nonlinearity due to $`\alpha `$–quenching on the behaviour of solar torsional oscillations. We feel that the previous omission of any discussion of these effects reduced the detailed applicability of our results to the observational data.

To study the effects of density stratification, first of all we made a comparison with our previous results and found both important similarities and differences. Concerning similarities, we found that consistent with our previous findings for near-critical and moderately supercritical dynamo regimes, the torsional oscillations extend all the way down to the bottom of the convection zone. This is true in presence of density stratification with all values of $`\rho _{\mathrm{ratio}}`$ in the range $`[1,10^6]`$. Concerning differences, we found that in both the near–critical and supercritical cases, the amplitudes of the torsional oscillations in the upper part of the CZ increase with $`\rho _{\mathrm{ratio}}`$, while those at the bottom decrease. Given that the uncertainties in the inversions are expected to be greatest at the bottom of the CZ, the lowering of amplitudes in these regions with increasing density stratification has potential importance for the analysis of helioseismological data, and for the understanding of the dynamo regimes near the bottom of the convection zone.

When two forms of nonlinearity are simultaneously present, we find that torsional oscillations persist in presence of both nonlinearities, for a large range of values of the $`\alpha `$–quenching coefficient $`g`$ up to $`g>1000`$. The amplitudes of these oscillations decrease however, as they shadow the square of the magnetic field strength, which in turn follows the scaling with $`g`$ given by (7).

We also note that for small $`g`$, the combined effect of increasing $`\rho _{\mathrm{ratio}}`$ and $`g`$ is to leave the amplitude of the near-surface torsional oscillations little changed. However, for larger $`g`$ ($`10`$), the effect of $`\alpha `$-quenching dominates, as the amplitudes of the torsional oscillations scale as $`1/g`$.

Equation (2) suggests, and our models confirm, that variations in $`v^{}`$ closely shadow those in the $`\varphi `$-component of the Lorentz force. It seems a priori difficult to imagine a driver for torsional oscillations that does not involve the Lorentz force, either at a macro- or micro-level (e.g. the ‘nonlinear $`\mathrm{\Lambda }`$–effect’, Kitchatinov & Pipin 1998; Kitchatinov et al. 1999), whatever the prime nonlinearity in the dynamo. Thus the requirement that the Lorentz force shadows the torsional oscillations is plausibly an additional constraint even on dynamo models that do not explicitly include dynamical effects (e.g. conventional $`\alpha `$-quenched models. Brandenburg & Tuominen 1988 appear to have been the first to investigate semi-quantitatively such an idea.) However it is perhaps a little puzzling that values of $`g`$ as large as suggested by the estimates of Sect. 2.3 reduce the amplitudes of the torsional oscillations so substantially, removing the previous agreement with observations.


format: text

text:
# Challenges in Moving the LEP Higgs Statistics to the LHC

## I Introduction

Higgs searches at LEP were based on marginal signal expectations and small background uncertainties. In contrast, Higgs searches at the LHC are based on strong signal expectations and relatively large background uncertainties. Based on our experience with the LEP Higgs search, our group tried to move the tools we had developed at LEP to the LHC environment. In particular, our calculation of confidence levels was based on an analytic computation with the Fast Fourier Transform and the log-likelihood ratio as a test statistic (and systematic errors based on the Cousins-Highland approach). We encountered three types of problems when calculating ATLAS’ combined sensitivity to the Standard Model Higgs Boson: problems associated with large numbers of expected events, problems arising from very high significance levels, and problems related to the incorporation of systematic errors.

Previously, it was shown that the migration of the statistical techniques that were used in the LEP Higgs Working Group to the LHC environment is not as straightforward as one might naïvely expect UW-confidence-VBF . After a brief overview in Section II, those difficulties and their ultimate solution are discussed in Section III. Our group has developed two independent software solutions (both in C++; both with FORTRAN bindings; one ROOT based and the other standalone) which can be found at:

http://wisconsin.cern.ch/software

In Section IV we discuss the incorporation of systematic errors and compare a few different strategies. In Section V we present and discuss the discovery luminosity (the luminosity expected to be required for discovery). Lastly, in Section VI we discuss the statistical notion of power (which is related to the probability of Type II error (the probability we do reject the “signal-plus-background hypothesis” when it is true).

## II The Formalism

Our starting point for this note is a brief review of the techniques that were used at LEP. We refer the interested reader to Kendall  for an introduction to the fundamentals, to Read:1997  for why the likelihood ratio has been chosen as a test statistic, to Junk:1999kv  for a Monte Carlo approach to the calculation and to clfft  for the analytic calculation using Fast Fourier Transform (FFT) techniques. For completeness, we introduce the basic approach below using the notation found in UW-confidence-VBF . For a counting experiment where we expect, on average, $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ background events and $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}`$ signal events, we consider two hypotheses: the null (or background-only) hypothesis in which the number of expected events, $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}`$, is described by a Poisson distribution $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{P}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{;}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}`$ and the alternate (or signal-plus-background) hypothesis in which the number of expected events is described by a Poisson distribution $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{P}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{;}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{+}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}`$. Here the number of events serves the purpose of a test statistic: a real number which quantities an experiment.

It is possible to include a discriminating variable $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{x}$}`$ which has some probability density function (pdf) for the background, $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{f}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{x}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}`$, and some pdf for the signal, $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{f}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{x}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}`$, both normalized to unity. Given an observation at $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{x}$}`$ we can construct the Likelihood Ratio $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{Q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{=}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{f}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{x}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{+}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{f}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{x}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{/}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{f}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{x}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}`$. With several independent observations $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\{}$}\widehat{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{x}$}}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\}}$}`$ we can consider the combined likelihood ratio $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{Q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{=}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{Q}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}`$ . It is possible, and in some sense optimal, to use $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{Q}$}`$ (or in practice $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{=}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\mathrm{ln}}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{Q}$}`$) as a test statistic.

The computational challenge of using the log-likelihood ratio in conjunction with a discriminating variable $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{x}$}`$ is the construction of the log-likelihood ratio distribution for the background-only hypothesis, $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}`$, and for the signal-plus-background hypothesis $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{+}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}`$. In this case, there are not only the Poisson fluctuations of the number of events, but also the continuously varying discriminating variable $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{x}$}`$. In particular, for a single background event the log-likelihood ratio distribution, $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{,}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}`$, must incorporate all possible values of $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{x}$}`$. From these single event distributions we can build up the expected log-likelihood ratio distribution by repeated convolution. This is most effectively done by using a Fast Fourier Transform (FFT) where convolution can be expressed as multiplication in the frequency domain (denoted with a bar). In particular we arrive at:

$`\overline{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}}`$ $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{=}$}`$ $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{e}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{[}$}\overline{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{,}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{]}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\mathrm{and}}$}`$ (1)
$`\overline{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{+}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}}`$ $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{=}$}`$ $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{e}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{+}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{[}$}\overline{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{,}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{+}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{]}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{.}$}`$
From the log-likelihood distribution of the two hypotheses we can calculate a number of useful quantities. Given some experiment with an observed log-likelihood ratio, $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}}`$, we can calculate the background-only confidence level, $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{C}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{L}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ :

$$\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{C}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{L}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{=}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\mathrm{}}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{𝑑}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}}$$
(2)
In the absence of an observation we can calculate the expected $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{C}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{L}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ given the signal-plus-background hypothesis is true. To do this we first must find the median of the signal-plus-background distribution $`\overline{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{+}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}}`$. From these we can calculate the expected $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{C}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{L}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ by using Eq. 2 evaluated at $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{=}$}\overline{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{+}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}}`$.

Finally, we can convert the expected background confidence level into an expected Gaussian significance, $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{N}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$, by finding the value of $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{N}$}`$ which satisfies

$$\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{C}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{L}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\overline{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{+}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{=}$}\frac{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\mathrm{erf}}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{N}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{/}$}\sqrt{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{2}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}}{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{2}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{.}$}$$
(3)
where $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\mathrm{erf}}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{N}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{=}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{2}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{/}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\pi }$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{0}$}^\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{N}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\mathrm{exp}}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{y}$}^\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{2}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{𝑑}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{y}$}`$ is a function readily available in most numerical libraries.

## III  Numerical Difficulties

The methods described in the previous section have been applied to the combined ATLAS Higgs effort with some caveats related to numerical difficulties UW-confidence-VBF . In particular, in the extreme tails of $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}`$, the probability density is dominated by numerical noise. This numerical noise is an artifact of round-off error in the double precision numbers used in the Fast Fourier Transform111We use the FFTW library: http://www.fftw.org. The noise is on the order of $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{10}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{17}$}}`$ (for double precision floating point numbers), which translates into a limit on the significance of about $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{8}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$. For particular values of the Higgs mass, ATLAS has an expected significance well above $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{8}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ with only $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{10}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\mathrm{fb}}$}^\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}`$ of data. In order to produce significance values above the $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{8}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ limit, various extrapolation methods were used in UW-confidence-VBF . We now introduce a definitive solution to this problem based on arbitrary precision floating point numbers.

It should be made clear that the numerical precision problem is not due to the fact that the $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{C}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{L}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ is so small that the evaluation of the integral in Eq. 2 cannot be treated with double precision floating point numbers. Instead, the numerical precision problem is due to the many (approximately $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{2}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{20}$}}`$) Fourier modes which must in total produce a number very close to $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{0}$}`$. In order to rectify this problem we have implemented the Fast Fourier Transform with the arbitrary-precision floating point numbers provided in the CLN library222CLN is available at http://www.ginac.de GiNaC . One might protest that above $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{5}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ we are not interested in the precise value of the significance and that this exercise is purely academic. We refer the interested reader to Sections V & VI for different summaries of an experiments discovery potential.

### III.1 Extrapolation

While the arbitrary precision FFT approach is the definitive solution to the problem of calculating very high expected significance, it is also incredibly time consuming. A much faster, approximate solution is to approximate the $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{C}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{L}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ by fitting the $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ distribution to a functional form. The first method of extrapolation studied was a simple Gaussian fit to the $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ distribution. This method works fairly well, but tends to overestimate the significance. The second method we studied was based on a Poisson fit to the $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ distribution. The Poisson distribution has the desirable properties that it will have no probability below the hard limit $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}`$ and that its shape is more appropriate UW-confidence-VBF . Figure 2 compares these different extrapolation methods.

## IV Incorporating Systematic Uncertainty

One encounters both philosophical and technical difficulties when one tries to incorporate uncertainty on the predicted values $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}`$ and $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ found in Eq. 1. In a Frequentist formalism the unknown $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}`$ and $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ become nuisance parameters. In a Bayesian formalism, $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}`$ and $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ can be marginalized by integration over their respective priors. At LEP the practice was to smear $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ and $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{+}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}}`$ by integrating $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}`$ and $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}`$ with a multivariate normal distribution as a prior. This smearing technique is commonly referred to as the Cousins-Highland Technique, and it is has some Bayesian aspects.

### IV.1 A Purely Frequentist Technique

At the PhysStat2003 conference a purely frequentist approach to hypothesis testing with background uncertainty was presented Cranmer2003freq . This method relies on the full Neyman construction and uses a likelihood ratio similar to the profile method as an ordering rule. In this formalism, a systematic uncertainty at the level of 10% has a much larger effect than when treated with the Cousins-Highland technique.

### IV.2 The Cousins Highland Technique

The Cousins-Highland formalism for including systematic errors on the normalization of the signal and background is provided in Cousins:1992qz  and generalized in Junk:1999kv ; clfft . In particular, for a multivariate normal distribution333In principle, any distribution could be used within this framework. as a prior for the $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}`$ the distribution of the log-likelihood ratio is given by:

$`\overline{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{y}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{=}$}{\displaystyle \colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\mathrm{}}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{e}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}^\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{K}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{[}$}\overline{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{,}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{q}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{]}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\left(}$}\frac{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}}{\sqrt{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{2}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\pi }$}}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\right)}$}^\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{K}$}\frac{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}}{\sqrt{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{|}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{S}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{|}$}}}}`$ (4)
$`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{e}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}^\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{K}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{j}$}^\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{K}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\frac{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}}{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{2}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{S}$}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{j}$}}^\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{j}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{j}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}}{\displaystyle \underset{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}}{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{d}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}`$
where $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{S}$}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{j}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{=}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{i}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{j}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{n}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{j}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}`$. Reference clfft  provides an analytic expression for the resulting log-likelihood ratio distribution including a correlated error matrix; however, this equation was obtained with an integration over negative numbers of expected events and does not hold. Attempts to provide a closed form solution for the positive semi-definite region require analytical continuation of the error function over a wide range of the complex plane. Instead, a numerical integration over the positive semi-definite region has been adopted for our software packages.

## V Discovery Luminosity

Because the calculation of expected significance is technically very difficult at the LHC, other summaries of the discovery potential have been explored. While these techniques are not new, it is important to consider their pros and cons. One such alternate summary of the discovery potential is based on the discovery luminosity”. Define the discovery luminosity, $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{L}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{m}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{H}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}`$, to be the integrated luminosity necessary for the expected significance to reach $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{5}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$. The discovery luminosity is an informative quantity; however, it must be interpreted with some care:

* Collecting an integrated luminosity equal to the nominal discovery luminosity does not guarantee that a discovery will be made. Instead, with $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{L}$}^{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{(}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{m}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{H}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{)}$}`$ of data the median of $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\rho }$}_{\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{s}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{+}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{b}$}}`$ will be at the $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{5}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ level – which corresponds to a 50% chance of discovery. See Section VI for more details.
* In practice an analysis’ cuts, systematic error, and signal and background efficiencies are luminosity-dependent quantities. When we calculate the discovery luminosity, we treat the analysis as constant.

## VI The Power of a $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\mathrm{𝟓}}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ Test

The traditional quantity which is used to summarize an experiment’s discovery potential is the combined significance; however, as was noted in Section III this plot becomes very dificult to make when the significance goes beyond about $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{8}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$. Furthermore, the plot itself starts to loose relevance when the significance is far above $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{5}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$. The discovery luminosity is another possible way of illustrating an experiment’s discovery potential, but it must be interpreted with some care. A third summary of an experiment’s discovery potential which is related to the probability of Type II error: the power. First, it should be noted that the expected significance is a measure of separation between the medians of the background-only and signal-plus-background hypotheses. Thus, when we see the significance curve cross the $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{5}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ line in Fig. 2 there is only a 50% chance that we would observe a $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{5}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ effect if the Higgs does indeed exist at that mass. In practice, we claim a discovery if the observed data exceeds the $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{5}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ critical region, and do not claim a discovery if it doesn’t. The meaning of the $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{5}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ discovery threshold is a convention which sets the probability of Type I error to be $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{2.85}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{10}$}^\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{7}$}`$ . With that in mind, the idea that the significance is $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{20}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ at $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{m}$}_\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{H}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{=}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{160}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\mathrm{GeV}}$}`$ is irrelevant. What is relevant is the probability that we will claim discovery of the Higgs if it is indeed there: that quantity is called the power. The power is defined as $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{1}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\beta }$}`$ where $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\beta }$}`$ is the probability of Type II error: the probability that we reject the signal-plus-background hypothesis when it is true Kendall .

Consider Figure 3 with a background expectation of 100 events. The black vertical arrow denotes the $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{5}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ discovery threshold. The (red) dashed curve shows the distribution of the number of expected events for a signal-plus-background hypothesis with 150 events. Normally, we would say the expected significance is $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{5}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ for this hypothesis; however, we can see that only 50% of the time we would actually claim discovery. The rightmost (blue) curve shows the distribution of the number of expected events for a signal-plus-background hypothesis with 180 events. Normally, we would say the expected significance is $`\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{8}$}\colorbox[rgb]{1,1,1}{$\textcolor[rgb]{0,0,0}{\sigma }$}`$ for this hypothesis; however, a more meaningful quantity – the power – is associated with the probability we would claim discovery which is about 98%. In addition to the power being a germane quantity, it is much easier to calculate.

## VII Conclusion

In conclusion, the migration of the statistical tool-set developed at LEP to the LHC environment is not as straightforward as one might expect. The first difficulties are computational and arise from the combination of channels with many events and channels with few events (these are easily solved). The next difficulties are numerical and arise from the extremely high expected significance of the high-energy frontier. These problems can be solved by brute force; or they can be reinterpreted as conceptual problems, and solved by asking different questions (i.e. power). Lastly, there is a philosophical split related to the Bayesian and Frequentist approach to uncertainty. At the LHC, the choice of the formalism is no longer a second-order effect, and this problem is not so easy to solve.

###### Acknowledgements.

This work was supported by a graduate research fellowship from the National Science Foundation and US Department of Energy Grant DE-FG0295-ER40896.


text:
# Schottky barriers at metal-finite semiconducting carbon nanotube interfaces

## Abstract

Electronic properties of metal-finite semiconducting carbon nanotube interfaces are studied as a function of the nanotube length using a self-consistent tight-binding theory. We find that the shape of the potential barrier depends on the long-range tail of the charge transfer, leading to an injection barrier thickness comparable to half of the nanotube length until the nanotube reaches the bulk limit. The conductance of the nanotube junction shows a transition from tunneling to thermally-activated transport with increasing nanotube length.

Single-wall carbon nanotubes (SWNT) are ideal systems for studying transport in the length scale ranging from the molecular limit as all-carbon cylindrical molecules to the bulk limit as quasi-one-dimensional conductors, Dekker99  and many device concepts have been sucessfully demonstrated on a single-tube basis.Dekker ; Lieber00 ; Avouris02  Among the device physics problems raised, the nature of the electronic transport through a metal-semiconducting carbon nanotube interface stands out Xue99 ; TersoffNT ; Odin00 ; De02  as one of the basic device building blocks. Rhod  Although transport through a metal-long carbon nanotube interface has been studied using the bulk band structure and one-dimensional electrostatics, TersoffNT ; Odin00 ; De02  it is important to investigate device functionality of finite carbon nanotubes with lengths of nanometer range and three-dimensional electrostatics, which will also shed light on the scaling limit of carbon nanotube devices.

Here, we analyze the evolution of the electronic properties of the metal-finite SWNT interface as the length is varied from the molecular to the bulk limit. The model system ia illustrated schematically in Fig. 1, where the open-ended SWNT molecule is attached to the surfaces of the metallic electrodes through the dangling bonds at the end. The device structure chosen represents an atomic-scale analogue of the metal-semiconductor interface since both the interface area and the active device region are atomic-scale. We find that the shape of the potential barrier depends on the charge transfer throughout the junction, which leads to an injection barrier thickness comparable to half of the nanotube length until the nanotube reaches the bulk limit. As a consequence, the SWNT junction conductance shows a transition from tunneling to thermally-activated transport as the nanotube length increases.

We take $`(10,0)`$ SWNT as the protype semiconducting SWNT, whose work function is taken as that of the graphite ($`4.5`$ eV). Dekker  We consider gold (Au) and titanium (Ti) electrodes as examples of high- and low- work function metals ($`5.1`$ and $`4.33`$ eV respectively for polycrystalline materials). We describe the electronic structure of the isolated SWNT using the Extended Huckel Theory (EHT) with non-orthogonal basis sets $`\varphi _m(\stackrel{}{r})`$, Hoffmann88  which gives a bulk (infinitely long) band gap of $`0.9(eV)`$ for $`(10,0)`$ SWNT. Since the SWNT Fermi level is located approximately at mid-gap, before the contact formation, the gold Fermi-level ($`5.1`$ eV) lies below the valence band edge while the titanium Fermi-level ($`4.33`$ eV) lies between the mid-gap and conduction band edge of the bulk SWNT.

Since the screening of Coulomb interaction is relatively ineffective within the SWNT due to the reduced dimensionality, TersoffNT ; Odin00  an atomistic study of the electronic processes throughout the metal-SWNT-metal junction is needed. We use a self-consistent tight-binding theory based on the semi-empirical implementation of the self-consistent Matrix Green’s function (SCMGF) method for first-principles modeling of molecular-scale devices. Xue02  Given the EHT Hamiltonian $`H_0`$ of the isolated SWNT, we calculate the density matrix $`\rho _{ij}`$ and therefore the electron density of the equilibrium SWNT junction from

$`G^R`$ $`=`$ $`\{(E+i0^+)SH\mathrm{\Sigma }_L(E)\mathrm{\Sigma }_R(E)\}^1,`$ (1)
$`\rho `$ $`=`$ $`{\displaystyle \frac{dE}{2\pi }Imag[G^R](E)f(EE_F)}.`$ (2)
where the effect of coupling to the electrodes are included as matrix self-energy operators $`\mathrm{\Sigma }_{L(R)}`$ and calculated using tight-binding parameters (we use a nanotube end-surface distance of $`2.5\AA `$ here). Xue02 ; Papa86  Here $`S`$ is overlap matrix and $`f(EE_F)`$ is the Fermi distribution describing the electrodes (see Ref. Xue02,  for details). The SWNT Hamiltonian is now $`H=H_0+\delta V[\delta \rho ]`$ where $`\delta \rho `$ is the density of transferred charge and $`\delta V`$ is the induced change in the electrostatic potential.

To proceed with self-consistent calculation, we approximate the charge distribution as superposition of atom-centered charge distributions Frau98  $`\delta \rho (\stackrel{}{r})=_i\delta N_i\rho _i(\stackrel{}{r}\stackrel{}{r}_i)`$, where $`\delta N_i=(\rho S)_{ii}N_i^0`$ and $`N_i^0`$ is the number of valence electrons on atomic-site $`i`$ of the bare SWNT. $`\rho _i(\stackrel{}{r})=\frac{1}{N_{\zeta _i}}e^{\zeta _ir}`$ is a normalized Slater-type function ($`𝑑\stackrel{}{r}\rho _i(\stackrel{}{r})=1`$). Frau98  The exponent $`\zeta _i`$ is chosen such that $`𝑑\stackrel{}{r}𝑑\stackrel{}{r}^{}\rho _i(\stackrel{}{r})\rho _i(\stackrel{}{r}^{})/|\stackrel{}{r}\stackrel{}{r}^{}|=I_iA_i`$ Frau98 , where $`I_i(A_i)`$ are the atomic electron affinity (ionization potential). In this way, we obtain $`\delta V(\stackrel{}{r})=_i\delta N_iV_i(\stackrel{}{r}\stackrel{}{r}_i)`$, where $`V_i=𝑑\stackrel{}{r}^{}\rho _i(\stackrel{}{r}^{}\stackrel{}{r}_i)/|\stackrel{}{r}\stackrel{}{r}^{}|`$ can be evaluated analytically Frau98 . We take into account the image-potential effect by including within $`\delta \rho `$ both the tranfered-charge on the carbon atoms and their image charges, rather than imposing an image-type potential correction. The matrix elements of the potential $`\delta V_{mn}=𝑑\stackrel{}{r}\varphi _m^{}(\stackrel{}{r})\delta V(\stackrel{}{r})\varphi _n(\stackrel{}{r})`$ are calculated using two types of scheme: (1) If $`m,n`$ belong to the same atomic site $`i`$, we calculate it by direct numerical integration; (2) if $`m,n`$ belong to different atomic sites, we use the approximation $`\delta V_{mn}=S_{mn}(\delta V_{mm}+\delta V_{nn})/2`$.

The finite SWNT lengths investigated are $`2.0,4.1,8.4,12.6,16.9,21.2`$ and $`25.4`$ (nm), which corresponds to number of unitcells of $`5,10,20,30,40,50`$ and $`60`$ respectively and spans the entire range from the molecular limit to the bulk limit. Due to the dangling bonds at the end, charge transfer occurs between the end and the interior carbon atoms of the SWNT, which should be corrected self-consistently first and give the initial charge configuration $`N_i^0`$. Xue02  The calculated charge transfer and electrostatic potential change for the Au-SWNT-Au and Ti-SWNT-Ti junctions are shown in Fig. 2. The magnified view of the transferred charge and potential shift at the interface and in the middle of the SWNT appear in Fig. 3.

The charge transfer at the metal-SWNT interface reflects the bonding configuration change upon contact to the metallic surfaces, involving mainly the end carbon atoms and decaying rapidly into the interior of the SWNT molecule. Xue02  The electrostatic potential change is instead determined by the transferred charge throughout the metal-SWNT-junction due to the long-range Coulomb interaction, so its magnitude in the middle of the SWNT increases with increasing length until the finite SWNT reaches the bulk limit despite the small magnitude of the transferred charge in the middle of the SWNT. Note  This gives a maximal barrier thickness for electron injection of roughly $`10(nm)`$, i.e., half of the nanotube length where it reaches the bulk limit ($`50`$ unitcells). The magnitude of both the charge transfer and the potential shift at the metal-SWNT interface are approximately the same for all the finite SWNTs studied.

For the SWNT molecules which have approached the bulk limit, the interface perturbation of the electron states in the middle can be essentially neglected. Charging and “band” shift is determined by the shift of the local density of states (LDOS) in the middle of the SWNT relative to the metal Fermi-level. Note that due to the three-dimensional electrostatics of the metal-SWNT-metal junction, the change in the electrostatic potential energy induced by the transferred charge decays moving away from the cylindrical surface of the SWNT. Therefore *the shift of the LDOS along the SWNT axis doesn’t follow the shift in the potential energy*, different from the bulk metal-semiconductor interface. For the 60-unitcell SWNT, we find that for the Au-SWNT-Au junction the Fermi-level is located slightly below (by $`0.05`$eV) the mid-gap, while for the Ti-SWNT-Ti junction it is located above (by $`0.15`$eV) the mid-gap. This is consistent with the fact that in the interior of the SWNT the Fermi-level must lie within the band gap to ensure the small perturbation of the electron states there.

These charge transfer processes are often characterized as “charge-transfer doping”. For the Au-SWNT-Au (Ti-SWNT-Ti) junction, there is a small positive (negative) charge transfer of $`5.7\times 10^4`$ ($`6.4\times 10^5`$) per atom in the middle of the $`60`$-unitcell SWNT. The oscillation of transferred-charge in the middle of the SWNT is due to the two-sublattice structure of zigzag tube, leading to charge transfer among atoms within the unitcell. The SWNT is therefore “hole-doped” by contacting to high work function (Au) and “electron-doped” by contacting to low work function (Ti) electrode. However, within the coherent transport regime, it is clear that both the charge transfer at the interface and the charge transfer inside the SWNT contribute only indirectly to electron transport by modulating the potential landscape acrosss the metal-SWNT-metal junction.

Given the potential shift across the metal-SWNT interface, we can evaluate the length and temperature dependence of the SWNT junction conductance using the Landauer formula

$$G=\frac{2e^2}{h}𝑑ET(E)[\frac{df}{dE}(EE_F)]=G_{Tu}+G_{Th}$$
(3)
where the transmission is calculated from $`T(E)=Tr[\mathrm{\Gamma }_L(E)G^R(E)\mathrm{\Gamma }_R(E)G^A(E)]`$Xue02  Here we have separated the conductance into tunneling contribution $`G_{Tu}=\frac{2e^2}{h}T(E_F)`$ and thermal-activation contribution $`G_{Th}=GG_{Tu}`$. The result for room temperature conductance is shown in Fig. 4. The tunneling conductance (also the zero temperature conductance) for both junctions decreases exponentially with the SWNT length for SWNT longer than $`4.1(nm)`$ appropriate for tunneling transport though potential barriers with identical barrier height. The room-temperature conductance instead saturates with increasing SWNT length. This is because the tunneling is exponentially suppressed for the longer SWNT molecules, while the transport becomes dominated by thermal-acitivation over the top of the potential barrier, whose height is approximately the same for all the finite SWNTs investigated. For Ti-SWNT-Ti junction, this leads to a transition from tunneling to thermally-activated transport at roughly $`5(nm)`$. For Au-SWNT-Au junction, the thermal contribution is already larger than the tunneling contribution for the smallest SWNT studied at room temperature.

The long-range nature of the charge transfer and potential barrier suggests the gate-modulation of junction conductance in SWNT-based transistors may be achieved through the modulation of the injection barrier at the metal-natotube interface, in agreement with recent experimental works. Avouris02  The SWNT junction transport characteristics are sensitive to the shape of the potential change induced by both the gate voltage and source-drain voltage, whose effect increases with increasing SWNT length. Further analysis are needed to achieve a through understanding of nanotube-based devices.

This work was supported by the DARPA Molectronics program, the DoD-DURINT program and the NSF Nanotechnology Initiative.


text:
# Twisted accretion curtains in the intermediate polar FO Aquarii

## 1 Introduction

FO Aqr is an intermediate polar (IP), a magnetic variant of the cataclysmic-variable class of binary stars. These systems consist of a Roche-lobe-filling, late-type secondary and a white-dwarf primary, the latter having a strong magnetic field ($``$1–10 MG), inclined to the rotational axis of the star. The in-falling matter forms an accretion disc which is truncated at the magnetospheric boundary, where it threads onto the white dwarf’s magnetic field. The material is then channeled into ‘accretion curtains’ which guide the material towards the poles, at which stand-off shocks form, heating material to $``$108 K (see Patterson 1994 for an in-depth review).

The accretion geometry of FO Aqr has been the subject of much debate (e.g. Hellier 1991; Norton et al. 1992; Hellier 1993; Beardmore et al. 1998; the last three hereafter referred to as N92, H93 and B98 respectively), but it is now generally accepted that the system shows disc-overflow accretion, with most of the material accreting via disc-fed accretion curtains (causing a pulse at the 20.9-min spin cycle), and a component of the accretion stream overflowing the disc and coupling to the magnetic field directly (causing a pulsation at the beat period between the spin and 4.85-hr orbital cycles).

Whereas many IPs have quasi-sinusoidal spin pulses, FO Aqr shows a more complex pulse, including a broader dip and a narrow ‘notch’ feature, both of which are variable over time (N92, B98). Our aim here is to use the superior spectroscopic capability of the *XMM-Newton* satellite to investigate the spin pulse and thus the way in which material leaves the accretion disc and threads onto field lines.

### 1.1 The spin-cycle ephemeris

The white dwarf in FO Aqr was spinning down during the 1980s (Shafter & Macry 1987), but Steiman-Cameron, Imamura & Steiman-Cameron (1989) reported that the period-lengthening had almost stopped by the time of their 1987 observations. Osborne & Mukai (1989) and Patterson et al. (1998) found that FO Aqr then began spinning up, which is confirmed by the most recent ephemeris by Williams (2003). However, these ephemerides suffer from significant $`OC`$ jitter and develop a cycle-count ambiguity for times later than 1998 (Williams 2003). Thus we cannot compare the phasing of our data (from 2001) with previous reports with any confidence.

## 2 Observations

FO Aqr was observed by the *XMM-Newton* satellite (Jansen et al. 2001) for 37 ks on 2001 May 12. The X-ray MOS-1, MOS-2 (Turner et al. 2001) and PN (Strüder et al. 2001) cameras were operating in small window mode, using the medium filter. The optical/UV OM camera (Mason et al. 2001) was operating in fast mode, observing through the UVM2 filter (2050–2450 Å). The RGS cameras were also in operation; however the count-rate was too low to allow phase-resolved analysis, which is the goal of this paper, thus they are only mentioned in passing here. Some aspects of these data were previously reported by Cropper et al. (2002).

We analysed the data using the xmm-sas software v5.4.1. The source data were extracted from a circular region of radius 3.5 arcsec enclosing the source, with the background being taken from an annulus around this area. Only single or double pixel events with a zero quality flag were selected.

## 3 Lightcurves and Power Spectra

The X-ray and UV lightcurves are presented in Fig. 1, binned at 64 s and 32 s respectively. We show the sum of the MOS lightcurves; the PN lightcurve is very similar to the MOS curves and so is not shown. In any case, the flickering level is higher than the photon noise.

Power spectra of these lightcurves are given in Fig. 2. The dominant pulses are the 4.85-hr and 20.9-min orbital and spin pulses (we label the frequencies $`\mathrm{\Omega }`$ and $`\omega `$ respectively), and there is also power at the $`\omega `$ – $`\mathrm{\Omega }`$ and $`\omega +\mathrm{\Omega }`$ sidebands.

Power at the beat period ($`\omega `$ – $`\mathrm{\Omega }`$) is indicative of spin–orbit interaction in the accretion process. X-ray beat periods are generally interpreted as resulting from accreting material which couples directly to the white dwarf magnetosphere, without passing through an accretion disc (e.g. stream-fed accretion, or disc-overflow accretion; Hellier et al. 1989a; Hellier 1991; Wynn & King 1992).

However, we see approximately equal power at the $`\omega `$ – $`\mathrm{\Omega }`$ and $`\omega +\mathrm{\Omega }`$ sideband frequencies, which may result simply from an amplitude modulation of the spin pulse at the orbital timescale (see Warner 1986), rather than from a stream–magnetosphere interaction. To test this we used the light-curve reconstruction program of H93 to see whether we could reproduce the light curve using a spin pulse multiplied by an orbital modulation, without any intrinsic modulation at the $`\omega `$ – $`\mathrm{\Omega }`$ period.

The resulting model lightcurve is given in the bottom panel of Fig. 1, and the power spectrum of the residuals to this fitted model is shown in the lower panel Fig. 2. The absence of peaks significantly above the flickering implies that the model is adequate in reproducing the periodic components of the lightcurve. Thus we conclude that there was little or no disc-overflow accretion occurring at the time of the *XMM-Newton* observations.

Note that Hellier (1991) and B98 reported that the relative amplitude of the beat and spin pulses change with time, and suggested that the accretion mode of FO Aqr changes from periods with significant disc-overflow accretion to periods where accretion occurs exclusively via the disc. This phenomenon has also been seen in TX Col (Wheatley 1999).

We also note from Fig. 2 that our data show very little power at the 2$`\omega `$ (marked) and 4$`\omega `$ (not shown) frequencies, in contrast to previous data (B98); this indicates substantial changes in the pulse profile over time.

## 4 Spectroscopy

We extracted phase-averaged source and background spectra from the MOS-1, MOS-2 and PN cameras, and generated response and ancillary response files using the sas rmfgen and arfgen tasks. We fitted the data from all three cameras simultaneously in xspec, making no allowance for calibration systematics between the instruments.

The X-ray emission in an IP arises from plasma which is heated to X-ray temperatures at a stand-off accretion shock, and cools as it approaches the white dwarf surface (e.g. Aizu 1973; Cropper et al. 1999). This can be represented by several mekal emitters at different temperatures or by using a cemekl model, which assumes a continuous power-law distribution of mekal emitters.

A single-temperature mekal model (plus absorption) gave a fit with $`\chi _\nu ^2`$ = 1.64 and a temperature of 9 keV; adding a second mekal improved $`\chi ^2`$ by 18 with temperatures 0.2 keV and 9.3 keV; a third mekal improved the $`\chi ^2`$ value by a further 10 ($`\chi _\nu ^2`$ = 1.63), and gave temperatures of 0.2, 3.0 and 14 keV. A fourth mekal gave a negligible improvement in $`\chi ^2`$.

Alternatively, using the cemekl model to give a continuous distribution of temperatures resulted in a comparable $`\chi _\nu ^2`$ of 1.64. We use the 3-mekal model for phase-resolved analysis below; note, however, that we do not interpret the highest temperature as a shock temperature, since this would require the fitting of a full accretion-column model (e.g. Cropper et al. 1999) and is also affected by the likely presence of a reflection component (e.g. Beardmore, Osborne & Hellier 2000). Instead, the 3-mekal model gives an easily-used and flexible representation of the multi-temperature spectrum for studying spectral changes over spin phase.

Blackbody emission may be anticipated in an IP, owing to the heating of the white dwarf surface by the post-shock material, however a blackbody component in the model made negligible difference to the $`\chi ^2`$ value, so was not included.

The spectra show lines of 6.4-keV iron and 0.57-keV oxygen vii that are not reproduced by the 3-mekal model (see Fig. 3). The RGS spectrum confirms that the 0.57-keV line is a real feature. In what follows we added narrow Gaussian components to fit these lines.

As is usual in IPs, simple photoelectric absorption was not adequate to fit the spectrum (giving $`\chi _\nu ^2`$ = 3.44 for a simple photoelectric absorber acting on the 3-mekal emitter described above). Adding a partial covering absorber gave $`\chi _\nu ^2`$ = 1.87 and a second such component reduced the $`\chi _\nu ^2`$ to 1.63. The fitted parameters of the 3-mekal, 3-absorber model are given in Table 1. Note that the amount of absorption changes over spin phase (Section 6) so the parameters fitted to the phase-averaged spectrum will be a weighted average over spin phase.

The above $`\chi ^2`$ values make no allowance for normalisation or calibration systematics between the three instruments. As a test, we have also allowed the parameters of the the 3-mekal, 3-absorber model to optimise for each instrument separately. This reduces the overall $`\chi _\nu ^2`$ from 1.63 to 1.03. When we report the phase-resolved analysis we quote $`\mathrm{\Delta }`$$`\chi ^2`$ values between competing models. Although the absolute $`\chi ^2`$ is lower when allowing the parameters to optimise for each instrument, the $`\mathrm{\Delta }`$$`\chi ^2`$ values are similar.

One result of the above modelling is the finding of a best-fitting abundance of only 0.12 solar. This results from a relative lack of lines in the data compared to that predicted by the mekal code, which assumes that the plasma is optically thin and collisionally ionized. Thus, the different line emission could equally be due to optical depth effects or photo-ionization (see, e.g., Mukai et al. 2003, for such effects in IP spectra). To investigate whether the abundance is genuinely low, or whether the line emission is suppressed, we added a bremsstrahlung emitter for each of the mekal emitters with their temperatures fixed to those of the mekals; and fixed the abundance to 1 (i.e. solar). This allowed a fit to be obtained which was as good as that with a low abundance, but the difference in $`\chi ^2`$ was too low ($``$5) to indicate which of these models is most physical. Applying this model to the higher resolution RGS data showed that either method reproduced the data, provided the abundance was above $``$0.15. We use the variable-abundance model in what follows.

## 5 Orbital modulation

X-ray observations of FO Aqr often show a prominent modulation on the orbital cycle (e.g. N92), as is obvious in our UV and X-ray data (Fig. 1). The standard interpretation is that FO Aqr is a high-inclination system in which disc material, particularly at the stream–disc impact region, periodically obscures the line of sight to the white dwarf (e.g. H93). Note that Hellier, Mason & Cropper (1989b) reported a grazing optical eclipse in FO Aqr, implying an inclination of $``$65, although there is no sign of the white dwarf itself being eclipsed.

Other possible causes of an orbital modulation include stream-fed accretion in which the accretion footprints are linked to orbital phase (N92). This, though, is less able to explain our dataset given the absence of the X-ray beat modulation that would also result (Section 3).

Fig. 4 shows a softness-ratio plot of the orbital modulation, constructed using 1254-s bins to smooth out the spin pulse (we show the data unfolded since we have only $``$2.1 orbital cycles of data). The softness ratio closely follows the lightcurve, implying that the spectrum gets harder in the minima, supporting the absorption idea.

To investigate further we extracted spectra from the three X-ray cameras in the orbital dip ($`\varphi `$ = 0.85–1.2) and out of the dip ($`\varphi `$ = 1.2–0.85). Fitting the spectra with the model discussed in Section 4, we found that the spectra were adequately reproduced by a model in which only the absorption changed during the dip. This required increased partial-covering absorbers with a density near 3$`\times 10^{23}`$ cm-2, and gave a best-fitting $`\chi _\nu ^2`$ of 1.37 (with much of the excess $`\chi ^2`$ arising from systematics between the three X-ray cameras). Allowing other parameters to vary in addition to the absorption improved the $`\chi ^2`$ by only 18 (for $`\nu =2164`$), which is not significant.

From the softness ratio and spectral modelling, and the absence of beat-cycle modulations, we conclude that the orbital modulation is a ‘dipping’ event resulting from absorption by stream or disc material. The material must extend to $``$$`>`$25 above the plane, given the inclination quoted above, and must extend for at least 0.4 in orbital phase. This suggests that the stream/disc interaction region is quite extensive.

## 6 Spin Modulation

Fig. 5 shows the UV and summed X-ray lightcurves folded on the 1254-s spin period, along with the softness ratio of the X-ray data. We have phased our data such that spin-phase zero occurs at the maximum of the UV pulse, located by fitting a sine curve to those data. This phase zero corresponds to 0.23 on the Williams (2003) ephemeris.

The X-ray pulse profile shows a main dip (at phase 0.75 on our plots) followed by a ‘notch’ feature (phases 0.0–0.1), which concurs with previous reports by N92, H93, and B98. The softness ratio of the X-ray spin pulse follows the morphology of the total flux level.

### 6.1 Spectral analysis

To investigate the spin-pulse features noted above we extracted spectra from all three X-ray cameras from spin ‘maximum’ ($`\varphi `$$``$0.12–0.36), ‘minimum’ ($`\varphi `$$``$0.67–0.87), and the ‘notch’ ($`\varphi `$$``$0.99–1.11). We applied the three-mekal model described in Section 4 simultaneously to ‘maximum’ and ‘minimum’, allowing the parameters to vary between regions, and then fixing the emission or absorption, in order to investigate the probable cause of the changes. We then repeated this process between the ‘maximum’ and ‘notch’ phase regions. xspec could not always constrain the mekal temperatures, so we froze these at 0.1, 2 and 30 keV as this had negligible effect on the fit quality. We also tied the abundance across phase regions, since it had a minimal effect on the fit quality ($`\mathrm{\Delta }\chi ^2\mathrm{\hspace{0.17em}5}`$), and fixed the 0.57-keV line at its rest energy with zero width.

Fitting the maximum and minimum phase regions simultaneously without requiring any model parameters to be the same (with the exceptions noted above) yielded a best-fitting model of $`\chi _\nu ^2`$ = 1.2 ($`\nu =1161`$). Note, again, that much of the $`\chi ^2`$ comes from calibration uncertainties; optimising for each instrument separately gave $`\chi _\nu ^2`$= 0.87. Since these systematics dominate $`\chi ^2`$, true uncertainties are much greater than the probabilities drawn from formal statistics; thus we present our results using $`\mathrm{\Delta }\chi ^2`$ values rather than probabilities.

Constraining the model so that the absorption could not vary between maximum and minimum made the fit worse by $`\mathrm{\Delta }\chi ^2=53`$; if instead the absorption is allowed to vary and the emission parameters are fixed, the fit is worse by $`\mathrm{\Delta }\chi ^2=29`$. This suggests that the minimum is mostly the result of absorption.

We then fitted the maximum and notch simultaneously: allowing only the absorption to vary worsened the fit by $`\mathrm{\Delta }\chi ^2=39`$ compared to an unconstrained model, whereas allowing only the emission parameters to vary worsened the fit by only $`\mathrm{\Delta }\chi ^2=4`$. From this we conclude that the notch is probably not an absorption event, but is caused by a change in the underlying emission, such as the occultation of parts of the accretion region.

## 7 Discussion

FO Aqr has a complex X-ray spin-pulse profile. The main feature is a dip, centered at $`\varphi `$$``$0.75, which is caused by increased absorption (section 6.1). Strong absorption dips are characteristic of IPs, and in the standard model occur when the accretion curtains of infalling material obscure the line of sight to the accretion regions (e.g. Rosen, Mason & Córdova 1988; Hellier, Cropper & Mason 1991; Kim & Beuermann 1995).

Our data also show a smaller flux reduction, the ‘notch’ at $`\varphi `$$``$0.0. This notch is persistent in FO Aqr’s X-ray lightcurve (e.g. N92, B98) and is often more prominent than in our data. Our spectral analysis suggests that this is the disappearance of part of the emitting region, rather than an absorption dip.

N92 attributed the notch to the disappearance over the white-dwarf limb of a ‘hot-spot’ within the accretion region at the upper pole. H93 suggested that it is probably the upper pole itself that disappears: the width of the notch implies a region covering 0.001 of the white-dwarf area (N92), which is in line with the estimate of $`<`$ 0.002 for the size of the accreting poles in the eclipsing system XY Ari (Hellier 1997a).

Our simultaneous far-UV lightcurve shows a nearly-sinusoidal spin pulse with a maximum coincident with the notch ($`\varphi `$$``$0.0). In the standard accretion-curtain model, the optical and UV maximum occurs when the upper pole is on the far side of the white dwarf, furthest from the observer, and the accretion curtains are best displayed (Hellier et al. 1987). This is consistent with our interpretation of the notch. We do, though, require an asymmetry in the visibility of the upper and lower poles, otherwise the appearance of the lower pole would compensate for the disappearance of the upper pole and fill in the notch. This point was previously noted by B98, who showed that the notch lightcurve could be reproduced if the magnetic dipole were offset by 0.15 white-dwarf radii from the centre.

Given the above, a simple model would predict that the absorption dip caused by the upper accretion curtain passing through the line of sight would occur half a cycle later, at $`\varphi `$$``$0.5. In fact, it occurs at $`\varphi `$$``$0.77, late by 0.27 cycles. A simple shift in phase, such as accreting only along field lines trailing the pole, would not by itself explain the quarter-cycle discrepancy between phases deduced from the UV and the X-ray. Thus we conclude that the field lines are not an undistorted dipole, but are twisted, being swept backwards by a quarter of a cycle.

We can use B98’s compilation of FO Aqr lightcurves to investigate whether the amount of distortion is variable over time. The notch is prominent in *Ginga* data taken in 1988 and 1990, and in both cases implies an accretion-curtain lag comparable to the amount we see in our data. The notch is either less obvious or absent in *EXOSAT* data from 1983 and 1985, and in *ASCA* data from 1993, but the statistical quality of those lightcurves is much lower. Thus, where both a notch and a dip are seen, the phase-difference between them is compatible with a constant lag, and there is no evidence for episodes when the accretion-curtain twist is markedly different.

### 7.1 Comparison with PQ Gem

While many intermediate polars have quasi-sinusoidal spin pulses, PQ Gem, like FO Aqr, has a more complex pulse (e.g. Mason 1997). However, in that star it is proposed that the accreting field lines *precede* the magnetic pole, the opposite to our conclusion for FO Aqr, based on analysis of the X-ray pulse (Mason 1997), optical polarimetry (Potter et al. 1997), and optical spectroscopy (Hellier 1997b).

A possible line of enquiry is whether the twists in the accretion curtains – opposite in the two stars – are related to the accretion torques on the white dwarfs. In PQ Gem the white dwarf is spinning down on a $`2.4\times 10^4`$-yr timescale (Mason 1997). FO Aqr’s white dwarf has shown a more complex history, having been spinning down between 1981 and 1988, after which it changed to spinning up (Williams 2003).

The spin-pulse profiles of FO Aqr have varied over that time period (e.g. B98), however there is no simple change in the relative locations of dip and notch that could be attributed straightforwardly to the change from spin down to spin up. Thus, although it is tempting to note that in the spinning-down PQ Gem the field lines precede the magnetic pole, whereas in the (currently) spinning-up FO Aqr they trail the magnetic pole, any deductions from this are not straightforward.

Ghosh & Lamb (1979) argue that the net accretion torque on the white dwarf combines the spin-up torque from the angular momentum of the accreted material with the torques of field lines dragging through the accretion disc. The expected result is a spin period which oscillates around an equilibrium period, undergoing phases of spin up and spin down. The field–disc torques could be the cause of the twists in the accretion curtains.

Mason (1997) and Potter et al. (1997), on the basis of optical polarimetry, proposed that in PQ Gem the accreting material threads field lines outside the corotation radius (where the field lines are moving faster than the disc material) and that this leads to the spin-down seen in that star. However, Mason (1997) notes that accretion from outside the corotation radius contradicts the standard theory of how an accretion disc interacts with a magnetic field (e.g. Ghosh & Lamb 1979). Our results from FO Aqr mean that the theory must also be adaptable to explain accretion along trailing field lines, and also the puzzling fact that the lines appear to trail both in episodes of spin-up and of spin-down. Thus any interpretations of the field twists seen in PQ Gem and FO Aqr are currently uncertain.

## 8 Conclusions

Analysis of *XMM-Newton* data on FO Aqr has shown that at the time of our observations:

* There was negligible disc-overflow accretion.
* An orbital dip was probably caused by absorption of the X-ray emission by vertical structure in the accretion disc.
* A narrow ‘notch’ in the X-ray spin pulse is likely to be caused by the upper accretion region disappearing over the limb of the white dwarf.
* The UV spin pulse is quasi-sinusoidal with a maximum coincident with the notch. This will occur when the base of the upper pole is furthest from us, and accords with the accretion curtain model for UV/optical spin pulses.
* A broad dip in the X-ray pulse is caused by absorption, most likely occurring when the outer regions of the upper accretion curtain cross our line of sight to the white dwarf.
* The phase difference between the broad dip and the notch is 0.23, whereas a simple accretion-curtain model would predict 0.5. This implies that the curtains are twisted, with the outer parts trailing the magnetic pole.
* The finding of a trailing curtain contrasts with reports that the accretion curtains lead the magnetic poles in PQ Gem.
* Although PQ Gem’s white dwarf is spinning down, whereas that in FO Aqr is (currently) spinning up, there appears to be no simple relation between the period change and whether the curtain trails or leads.