format: text

text:
# Atomic model of susy Hubbard operators

## Abstract

We apply the recently proposed susy Hubbard operators \[P. Coleman, C. Pépin and J. Hopkinson, Phys. Rev. B 63, 140411(R) (2001)\] to an atomic model. In the limiting case of free spins, we derive exact results for the entropy which are compared with a mean field + gaussian corrections description. We show how these results can be extended to the case of charge fluctuations and calculate exact results for the partition function, free energy and heat capacity of an atomic model for some simple examples. Wavefunctions of possible states are listed. We compare the accuracy of large N expansions of the susy spin operators \[P. Coleman, C. Pépin and A. M. Tsvelik, Phys. Rev. B 62, 3852 (2000); Nucl. Phys. B 586, 641 (2000)\] with those obtained using ‘Schwinger bosons’ and ‘Abrikosov pseudo-fermions’. For the atomic model, we compare results of slave boson, slave fermion and susy Hubbard operator approximations in the physically interesting but uncontrolled limiting case of N$``$2. For a mixed representation of spins we estimate the accuracy of large N expansions of the atomic model. In the single box limit, we find that the lowest energy susy saddle-point reduces to simply either slave bosons or slave fermions, while for higher boxes this is not the case. The highest energy saddle-point solution has the interesting feature that it admits a small region of a mixed representation, which bears a superficial resemblance to that seen experimentally close to an antiferromagnetic quantum critical point.

1. INTRODUCTION

One of the interesting challenges arising in systems exhibiting strong electronic correlations is the ability to adequately develop a microscopic description of metallic physics in regions where Fermi liquid theory breaks down. The observations of $`\frac{B}{T}`$ scaling, linear resistivity over three decades in temperature coinciding with a T$`\mathrm{ln}[\frac{T_0}{T}]`$ divergence in the specific heat capacity in YbRh2Si2 is reminiscent of properties shown by CeCu6-xAux which additionally exhibits $`\frac{E}{T}`$ scaling. In these heavy fermion systems close to an antiferromagnetic quantum critical point only one energy or length scale seems to be important–which generally only occurs in theories which are below their upper critical dimension. Evidence for a qualitatively similar phase diagram for the high-Tc cuprates has been advanced by Tallon and Loram who have assembled a large amount of experimental data on T\* and interpret it in terms of an underlying quantum critical point, presumably responsible for the linear T dependence seen in the normal state of these systems that has given rise to the marginal Fermi liquid phenomenology of Varma and others. Recent $`\mu `$SR measurements by Panagopoulos et al and an as yet unsubstantiated report of high temperature superconductivity in a FET of CaCuO2 by Schon et al have supported this interpretation. Again, the culprit for this unusual metallic behaviour is believed to be an underlying quantum critical point. If the entire Fermi surface fundamentally breaks down in this region, one might imagine that the competition between spins trying to magnetically order and those trying to form heavy quasiparticles could lead to a new kind of quasi-particle excitation.  That is, if the Fermi temperature approaches 0 as we approach a quantum critical point, we may see the formation of novel kinds of states; perhaps states of higher symmetry classes could become important even at these low temperatures. The key question that must be addressed is: how do spin and charge interact at the brink of magnetism?

Theoretical approaches to this problem are hindered by the difficulty of capturing the profound transformation in spin correlations that develops at the boundary between antiferromagnetism and paramagnetism. Usually we model this behavior by representing the spin as a boson in a magnetic phase, or as a fermion in a paramagnetic phase, but by making this choice, the character of spin and charge excitations which appear in an approximate field theory is restricted and lacks the flexibility to describe the co-existence of strong magnetic correlations within a paramagnetic phase.

These considerations have motivated the development of new methods to describe the spin and charge excitations of a strongly correlated material which avoid making the choice between a bosonic or fermionic spin. Recently, we introduced a new representation of the Hubbard operators which attempts to treat magnetism and paramagnetism on an equal footing within an approximate field theory allowing a description of physics in terms of operators which obey canonical commutation relations, yet avoid forcing us to specify the nature of a system’s ground state. In this approach, we introduce vector fields

$$F_a=(f_1,\mathrm{},f_N,\varphi )B_a=(b_1,\mathrm{},b_N,\chi )$$
(1)
where $`f_\sigma `$ and $`b_\sigma `$ are Abrikosov pseudo-fermions and Schwinger bosons respectively, while $`\varphi `$ and $`\chi `$ are slave bosons and slave fermions. For $`|a>|\sigma >,|0>`$, the Hubbard operators can be written

$$|a><b|=X_{ab}=B_a^{}B_b+F_a^{}F_b$$
(2)
To guarantee that the Hubbard operator representation is irreducible, we need to set the values of the linear and quadratic Casimirs of the group (here $`g=Diag[1\mathrm{}1,1]`$),

$$C^{(1)}=\text{Tr}[X]=Q,C^{(2)}=\text{Tr}[XgX]=Q(N1Y).$$
(3)
The first constraint corresponds to summing the total number of spin and hole states at a site (henceforth labeled by Q), while the latter describes the degree of antisymmetry of a given representation (which is fixed by Y = h-w (see Fig.1)). In terms of canonical operators, we can express Q and Y as:

$$Q=n_b+n_\varphi +n_f+n_\chi ,$$
(4)
$$Y=n_\varphi +n_f(n_b+n_\chi )+\frac{1}{Q}[\theta ,\theta ^{}],$$
(5)
where $`\theta =_\sigma b_\sigma ^{}f_\sigma \chi ^{}\varphi `$ is an operator interconverting fermions and bosons for the corner state. While the physically interesting case (Fig.1 (a)) corresponds to Q=1, Y=0 and N=2, one can hope to learn new qualitative features from large N expansions about different symmetry classes as has been done previously for Q = Q0, Y = 1-Q0 ((Fig.1(c))slave fermions–magnetism); and Q = Q0, Y = Q0 \- 1 ((Fig.1(d))slave bosons-paramagnetism). One now can also treat cases like (Fig.1(b) Y=0 –paramagnetism + magnetism?).

In this paper we first consider the limit of free spins with strong correlations in which charge fluctuations are totally suppressed. In this limit, Q becomes the total number of spins ($`n_b+n_f`$) and the energy of the empty state is taken to infinity. We derive exact results for the entropy of a general Q,Y. Using a mean field approach plus gaussian corrections, we compute the entropy of the general free spin case and compare this with $`\frac{1}{N}`$ expansions of the exact results (Stirling’s approximation). We then compare the supersymmetric spin formalism with earlier large N approximations. Representations of spins in terms of N-component ‘Schwinger’ bosons and ‘Abrikosov pseudo-fermions’ become exact in the large N limit, but in the physical limit N $``$ 2, provide only approximate descriptions of the physics. Calculation of the exact entropy of free spin states, allows us to estimate the error inherent in mean field descriptions, mean field plus gaussian corrections, and self-consistent mean field + correction approximations to each of the three methods. Surprisingly, self-consistent approximations seem to not always do as well as mean field approximations at low N.

Re-introducing charge fluctuations, we next treat an atomic model. In such a model there is an excitation energy for a spin to leave its localized site, creating an empty or vacuum state. The Hamiltonian for such a system is particularly simple, being given by simply the energy difference between the unoccupied site and the occupied spinful state. In terms of Hubbard operators, this can be written:

$$H=E_dX_{\sigma \sigma }=E_dQE_dX_{00}=E_dX_{00}$$
(6)
where $`E_d`$ is the energy of the localized level, and we have dropped a constant term in the last term for convenience. Now, we have charge fluctuations in addition to spin fluctuations allowed by the system. For small N we can write exact wavefunctions corresponding to possible states of the system and their associated energies. Then, we can write correspondingly exact free energies, entropies and heat capacities for such a system. These exact values allow us at low N to compare the efficacy of the supersymmetric Hubbard operator formalism with that of previously existing slave boson and slave fermion approximations. We find that even in this uncontrolled regime, bosonic descriptions (slave fermions) provide the best approximation to spin states, while fermionic descriptions of the spin (slave bosons) provide the best approximation to charge states at N=2. The lowest energy saddle-point of the susy solution at N=2 corresponds to this solution, while the highest energy saddle-point generates a phase diagram which bears a superficial resemblance to the phase diagram of the cuprates and includes a small region where the solution is neither fermionic nor bosonic but mixed. For higher N, we show how one can generate the exact partition function and estimate the accuracy in calculations of the entropy and heat capacity of this new large N approach which allows representations of the spin neither bosonic nor fermionic in character.

2. FREE SPIN LIMIT

First, let us consider the free spin. This corresponds to taking the limiting case when we send $`E_d\mathrm{}`$.

A. Entropy

1. Exact Result

We are left with SU(N) spins as the only states of our system. To count the number of states available, we make use of the fact that the dimension of an irreducible representation of SU(N) can be expressed as:

$$\lambda _N=\underset{i,j}{}\frac{(N+ji)}{g_{ij}}$$
(7)
where $`g_{ij}`$ is the hook length of the (i,j) box of the Young tableau (see Fig. 2 a)) defined as the number of neighbors to the right and bottom of the box (see Fig. 2 b)). For the specific case of L-shaped Young tableaux, this yields,

$`\lambda _N`$ $`=`$ $`{\displaystyle \frac{N}{h+w1}}{\displaystyle \underset{i=2}{\overset{h}{}}}({\displaystyle \frac{N+1i}{hi+1}}){\displaystyle \underset{j=2}{\overset{w}{}}}({\displaystyle \frac{N+j1}{wj+1}})`$ (8)
$`=`$ $`{\displaystyle \frac{1}{h+w1}}{\displaystyle \frac{N!}{(Nh)!(h1)!}}{\displaystyle \frac{(N+w1)!}{(w1)!N!}}`$ (9)
$`=`$ $`\left(\begin{array}{c}N\\ h\end{array}\right)\left(\begin{array}{c}N+w1\\ w\end{array}\right){\displaystyle \frac{wh}{N(w+h1)}}`$ (14)
From this formula, one can deduce an exact expression for the entropy depending only on N (SU(N)), w (width of the Young tableau), and h (height of Young tableau). This equation is:

$$S=\mathrm{ln}\left[\left(\begin{array}{c}N\\ h\end{array}\right)\left(\begin{array}{c}N+w1\\ w\end{array}\right)\frac{wh}{N(w+h1)}\right]$$
(15)
Since we work in the large N limit it would be useful to know the leading dependence of this exact result. We can use Stirling’s approximation to expand the factorials, i.e.

$$\mathrm{ln}(x!)=x\mathrm{ln}(\frac{x}{e})+\frac{1}{2}\mathrm{ln}(2\pi x)$$
(16)
which, after some work, implies that

$$S_{Stirling}=N(S_f[\stackrel{~}{h}]+S_b[\stackrel{~}{w}])\frac{1}{2}(\mathrm{ln}(2\pi N\stackrel{~}{w}(1+\stackrel{~}{w}\frac{1}{N}))+\mathrm{ln}(2\pi N\stackrel{~}{h}(1\stackrel{~}{h})))+\mathrm{ln}(\frac{wh}{N(w+h1)})$$
(17)
where $`\stackrel{~}{h}=\frac{h}{N},h=\frac{Q+Y+1}{2},\stackrel{~}{w}=\frac{w}{N},w=\frac{QY+1}{2}`$ and

$`S_f[x]=[x\mathrm{ln}x+(1x)\mathrm{ln}(1x)]`$ (18)
$`S_b[x]=[x\mathrm{ln}x(1+x)\mathrm{ln}(1+x{\displaystyle \frac{1}{N}})]`$ (19)
In order to express h, and w as functions of N, we can remove the leading dependence of Q and Y on these variables and expand in it. We can choose the new variable (other than N) in such a way that its sign agrees with that of Y, and its maximal amplitude is 1 for the special case when $`Q=\frac{N}{2}`$. That is, we define the variable $`\xi `$ such that $`w=\frac{N}{4}(1\xi )+\frac{1+\xi }{2}`$ and $`h=\frac{N}{4}(1+\xi )+\frac{1\xi }{2}`$. With this assignment, for the case we will consider ($`Q=\frac{N}{2}`$), $`1\xi 1`$ and $`Y=(\frac{N}{2}1)\xi `$. The physics at $`N=2`$ will be invariant of the value of $`\xi `$. It does have the problem that $`N\mathrm{}`$ and $`\xi \pm 1`$ do not commute.

Following the above-developed procedure, we express the exact entropy of the description as ($`Q=\frac{N}{2}`$):

$$\frac{2S_{exact}}{N}=\frac{2}{N}\mathrm{ln}\left(\left(\begin{array}{c}N\\ h\end{array}\right)\left(\begin{array}{c}N+w1\\ w\end{array}\right)\frac{wh}{N(\frac{N}{2})}\right)$$
(20)
and expansion of this according to Stirling’s approximation gives us:

$`{\displaystyle \frac{2S_{Stirling}}{N}}`$ $`=`$ $`({\displaystyle \frac{1}{2}}+\xi ({\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{N}}))\mathrm{ln}(({\displaystyle \frac{1}{4}}+{\displaystyle \frac{1}{2N}})+\xi ({\displaystyle \frac{1}{4}}{\displaystyle \frac{1}{2N}}))`$ (21)
$``$ $`({\displaystyle \frac{3}{2}}\xi ({\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{N}}))\mathrm{ln}(({\displaystyle \frac{3}{4}}{\displaystyle \frac{1}{2N}})\xi ({\displaystyle \frac{1}{4}}{\displaystyle \frac{1}{2N}}))({\displaystyle \frac{1}{2}}+\xi ({\displaystyle \frac{1}{N}}{\displaystyle \frac{1}{2}}))\mathrm{ln}(({\displaystyle \frac{1}{4}}+{\displaystyle \frac{1}{2N}})`$ (22)
$`+`$ $`\xi ({\displaystyle \frac{1}{2N}}{\displaystyle \frac{1}{4}}))+({\displaystyle \frac{5}{2}}+\xi ({\displaystyle \frac{1}{N}}{\displaystyle \frac{1}{2}}))\mathrm{ln}(({\displaystyle \frac{5}{4}}{\displaystyle \frac{1}{2N}})+\xi ({\displaystyle \frac{1}{2N}}{\displaystyle \frac{1}{4}})){\displaystyle \frac{2}{N}}\mathrm{ln}(\pi N){\displaystyle \frac{1}{N}}\mathrm{ln}(1{\displaystyle \frac{1}{N}})`$ (23)
extracting terms of order $`\frac{1}{N}`$ from inside the logs, this becomes the leading order approximation:

$`{\displaystyle \frac{2S_{Stirling,lead}}{N}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1+\xi (1{\displaystyle \frac{2}{N}}))\mathrm{ln}(1+\xi ){\displaystyle \frac{1}{2}}(3+\xi (1+{\displaystyle \frac{2}{N}}))\mathrm{ln}(3\xi )`$ (24)
$``$ $`{\displaystyle \frac{1}{2}}(1+\xi (1+{\displaystyle \frac{2}{N}}))\mathrm{ln}(1\xi )+{\displaystyle \frac{1}{2}}(5+\xi (1+{\displaystyle \frac{2}{N}}))\mathrm{ln}(5\xi ){\displaystyle \frac{2}{N}}\mathrm{ln}(\pi N){\displaystyle \frac{2}{N}}`$ (25)
2. Mean field + Gaussian Corrections

(a)Susy. Since we cannot enforce the constraints $`Q=Q_0`$ and $`Y=Y_0`$ rigorously at the mean field level, (Y has a four-fermion term of order $`\frac{1}{N}`$), we would like to see how inclusion of these terms consistently affects the saddle point approximation we make. For free spins, inclusion of the constraints as Lagrange multipliers leads to an action of the form

$$S=d\tau (\underset{\sigma }{}\overline{f}_\sigma (\tau +\lambda +\zeta )f_\sigma +\underset{\sigma }{}\overline{b}_\sigma (\tau +\lambda \zeta )b_\sigma +\frac{\zeta }{Q}\underset{\sigma ,\beta }{}(\overline{b}_\sigma f_\sigma \overline{f}_\beta b_\beta \overline{f}_\sigma b_\sigma \overline{b}_\beta f_\beta )\lambda Q_0\zeta Y_0$$
(26)
The total free energy then is made up of five terms to order $`\frac{1}{N}`$,

$$F_{tot}=𝒟[f,\overline{f},b,\overline{b}]e^S=F_{mft}+F_{\frac{1}{N}g.c.}=NF_f+NF_b+F_\eta \lambda Q\zeta Y+F_{\delta \lambda _f}+F_{\delta \lambda _b}$$
(27)
where the last three terms correspond to calculating contributions due to

$$\mathrm{\Pi }_\eta =\frac{Q_0}{2\zeta }+\text{},\mathrm{\Pi }_{\delta \lambda _f}=\text{},and\mathrm{\Pi }_{\delta \lambda _b}=\text{},$$
(28)
where the two terms of order N would ordinarily define the location of the saddle point, and the three latter terms be evaluated at the location of the mean field saddle point. As this would correspond to setting $`\frac{1}{Q}<[\theta ,\theta ^{}]>=0`$ at the saddle point, we choose to include $`F_\eta `$ in determining the location of the saddle point itself. This has the effect of coupling the gaussian fluctuations in $`\lambda _f`$ and $`\lambda _b`$ such that the sum $`F_{\delta \lambda _f}`$ \+ $`F_{\delta \lambda _b}`$ is replaced by a single term $`F_{\delta _{\lambda _f},\delta _{\lambda _b}}`$. In the above, $`\lambda _f=\lambda +\zeta `$ and $`\lambda _b=\lambda \zeta `$ arise from Lagrange multipliers $`\lambda `$ and $`\zeta `$ which have been introduced to enforce $`Q=h+w1`$ and $`Y=hw`$ respectively at the mean field level. $`\eta `$ is a Hubbard-Stratonovich field introduced to decouple the O($`\frac{1}{N}`$) term from Y. Details of the calculation are given in Appendix A. Grouping all our terms from the Gaussian approximation to the entropy (calculated in Appendix A), we have

$`F`$ $`=`$ $`NT\mathrm{ln}(1+e^{\beta \lambda _f})+NT\mathrm{ln}(1e^{\beta \lambda _b})+T\mathrm{ln}({\displaystyle \frac{\stackrel{~}{n}_b+\stackrel{~}{n}_f}{2\stackrel{~}{n}_b\stackrel{~}{n}_f}})\lambda (h+w)\zeta (hw)`$ (29)
$`+`$ $`{\displaystyle \frac{T}{2}}\mathrm{ln}((2\pi )^2(N\stackrel{~}{n}_b(1+\stackrel{~}{n}_b)N\stackrel{~}{n}_f(1\stackrel{~}{n}_f)n_\alpha (1n_\alpha )(N\stackrel{~}{n}_b(1+\stackrel{~}{n}_b)+N\stackrel{~}{n}_f(1\stackrel{~}{n}_f))))`$ (30)
Since there is no energy term here $`F=TS`$ to calculate the entropy we simply need to divide by the temperature although we can also perform the constrained derivative $`S=\frac{F}{T}|_n=\frac{F}{T}\frac{F}{\lambda }\frac{\lambda }{T}|_n`$.

$`{\displaystyle \frac{2S_{mfc}}{N}}`$ $`=`$ $`2\mathrm{ln}(1+e^{\beta \lambda _f})2\mathrm{ln}(1e^{\beta \lambda _b})+2{\displaystyle \frac{\beta \lambda _b}{N}}w+2{\displaystyle \frac{\beta \lambda _f}{N}}h){\displaystyle \frac{2}{N}}\mathrm{ln}({\displaystyle \frac{\stackrel{~}{n}_b+\stackrel{~}{n}_f}{\stackrel{~}{n}_b\stackrel{~}{n}_f}})`$ (31)
$``$ $`{\displaystyle \frac{1}{N}}\mathrm{ln}((2\pi )^2(N\stackrel{~}{n}_b(1+\stackrel{~}{n}_b)N\stackrel{~}{n}_f(1\stackrel{~}{n}_f)n_\alpha (1n_\alpha )(N\stackrel{~}{n}_b(1+\stackrel{~}{n}_b)+N\stackrel{~}{n}_f(1\stackrel{~}{n}_f))))+{\displaystyle \frac{2}{N}}(\mathrm{ln}(2))`$ (32)
where we have used the extra $`\beta \lambda `$ from the constraint on Q to create the third to last term as shown in Appendix A. Here we have $`n_\alpha =\frac{1}{e^{\beta 2\zeta }+1}`$, $`\stackrel{~}{n}_f=\frac{1}{e^{\beta \lambda _f}+1}`$ and $`\stackrel{~}{n}_b=\frac{1}{e^{\beta \lambda _b}1}`$ where $`\lambda _f=\lambda +\zeta `$ and $`\lambda _b=\lambda \zeta `$ are the Lagrange multipliers as they affect fermions and bosons respectively. The last term we remove as the factor of two from the double counting of the corner box. Since $`2\stackrel{~}{h}=q+y+\frac{1}{N}`$ and $`2\stackrel{~}{w}=qy+\frac{1}{N}`$, the extra $`\lambda `$ contribution allows us to neatly group prefactors in terms of $`\stackrel{~}{h}`$ and $`\stackrel{~}{w}`$. That is:

$`{\displaystyle \frac{2S_{mfc}}{N}}`$ $`=`$ $`2(1\stackrel{~}{h})\mathrm{ln}(1\stackrel{~}{n}_f)2\stackrel{~}{h}\mathrm{ln}\stackrel{~}{n}_f+2(1+\stackrel{~}{w})\mathrm{ln}(1+\stackrel{~}{n}_b)2\stackrel{~}{w}\mathrm{ln}\stackrel{~}{n}_b{\displaystyle \frac{2}{N}}\mathrm{ln}({\displaystyle \frac{\stackrel{~}{n}_b+\stackrel{~}{n}_f}{\stackrel{~}{n}_b\stackrel{~}{n}_f}})`$ (33)
$``$ $`{\displaystyle \frac{1}{N}}\mathrm{ln}((2\pi )^2(N\stackrel{~}{n}_b(1+\stackrel{~}{n}_b)N\stackrel{~}{n}_f(1\stackrel{~}{n}_f)n_\alpha (1n_\alpha )(N\stackrel{~}{n}_b(1+\stackrel{~}{n}_b)+N\stackrel{~}{n}_f(1\stackrel{~}{n}_f))))`$ (34)
Taking $`\frac{F_{mft}}{\stackrel{~}{n}_f}=0=\frac{F_{mft}}{\stackrel{~}{n}_b}`$ where the mean field contribution to the free energy shown in Eq. (18) is given by:

$$F_{mft}=NT((1\stackrel{~}{h})\mathrm{ln}(1\stackrel{~}{n}_f)+\stackrel{~}{h}\mathrm{ln}(\stackrel{~}{n}_f)(1+\stackrel{~}{w})\mathrm{ln}(1+\stackrel{~}{n}_b)+\stackrel{~}{w}\mathrm{ln}(\stackrel{~}{n}_b))+T\mathrm{ln}(\frac{\stackrel{~}{n}_b+\stackrel{~}{n}_f}{\stackrel{~}{n}_b\stackrel{~}{n}_f})$$
(35)
yields the saddle-point equations:

$`\stackrel{~}{n}_f+{\displaystyle \frac{1}{N}}(1n_\alpha )=\stackrel{~}{h}`$ (36)
$`\stackrel{~}{n}_b+{\displaystyle \frac{1}{N}}n_\alpha =\stackrel{~}{w}`$ (37)
where $`n_\alpha =\frac{1}{e^{2\beta \zeta }+1}=\frac{1}{\frac{\frac{1}{\stackrel{~}{n}_f}1}{\frac{1}{\stackrel{~}{n}_b}+1}+1}`$ is the gaussian contribution to the $`\frac{1}{Q}<[\theta ,\theta ^{}]>`$ term of the constraints which we need to include in the saddle-point (ie. at the mean field level). With this assignment, we see that at the mean field level we recover the results given by the earlier Abrikosov pseudofermion and Schwinger boson approaches. Analytic solution of these equations shows a recovery of the constraints for these two cases when the Young tableau is a column or a row respectively, while for Y = 0 we have in general a mixed case (see Fig. 3b).

(b) Abrikosov fermions. The exact entropy (normalized) of a column of height $`\frac{N}{2}`$ and width 1, is given by:

$$\frac{2S}{N}=\frac{2}{N}\mathrm{ln}\left(\begin{array}{c}N\\ \frac{N}{2}\end{array}\right)$$
(38)
Using Stirling’s approximation to expand the factorials, we obtain

$$\frac{2S_{Stirling^{}s}}{N}=2\mathrm{ln}2\frac{1}{N}\mathrm{ln}(\frac{\pi N}{2})$$
(39)
which should be the best that a $`\frac{1}{N}`$ theory can do. With self-consistent corrections, one might hope to approach the exact solution.

For Abrikosov pseudofermions the mean field + gaussian corrections free energy is given by:

$`F`$ $`=NT\mathrm{ln}(1+e^{\beta \lambda })\lambda Q+{\displaystyle \frac{T}{2}}\mathrm{ln}(2\pi N\stackrel{~}{n}_f(1\stackrel{~}{n}_f))`$ (40)
$`=NT\mathrm{ln}(1\stackrel{~}{n}_f)NqT\mathrm{ln}({\displaystyle \frac{1\stackrel{~}{n}_f}{\stackrel{~}{n}_f}})+{\displaystyle \frac{T}{2}}\mathrm{ln}(2\pi N\stackrel{~}{n}_f(1\stackrel{~}{n}_f))`$ (41)
where the last term is the gaussian corrections to the mean field result. Taking the usual mean field saddle-point $`\frac{F_{mft}}{n_f}=0`$ (mfc) we obtain $`\stackrel{~}{n}_f=q`$, such that the entropy (at $`q=\frac{1}{2}`$),

$$\frac{2S_{mfc}}{N}=2\mathrm{ln}2\frac{1}{N}\mathrm{ln}(\frac{\pi N}{2})$$
(42)
reproduces the Stirling’s approximation to the exact result. An attempt at a fully ’self-consistent’ (fsc) description might set $`\stackrel{~}{n}_f=\frac{q+\frac{1}{2N}}{1+\frac{1}{N}}`$ but such attempts do not seem to improve the accuracy.

(c) Schwinger Bosons. For the simple case of a purely bosonic description of spins, the exact entropy (for a row of width $`\frac{N}{2}`$ and height 1) is given by:

$$\frac{2S_{exact}}{N}=\frac{2}{N}\mathrm{ln}\left(\begin{array}{c}N+\frac{N}{2}1\\ \frac{N}{2}\end{array}\right)$$
(43)
which using Stirling’s approximation simplifies to:

$`{\displaystyle \frac{2S_{Stirling}}{N}}=(3{\displaystyle \frac{1}{N}})\mathrm{ln}({\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{N}})+\mathrm{ln}2{\displaystyle \frac{1}{N}}\mathrm{ln}(\pi N)+(2+{\displaystyle \frac{1}{N}})\mathrm{ln}(1{\displaystyle \frac{1}{N}})`$ (44)
which to leading order in N (ie dropping $`\frac{1}{N}`$ inside log’s) yields

$$\frac{2S_{Stirling,lead}}{N}=3\mathrm{ln}(\frac{3}{2})+\mathrm{ln}(2)\frac{1}{N}\mathrm{ln}(\frac{3\pi N}{2})$$
(45)
The mean field + gaussian corrections free energy is given by

$`F`$ $`=`$ $`NT\mathrm{ln}(1e^{\beta \lambda })\lambda Q+{\displaystyle \frac{T}{2}}\mathrm{ln}(2\pi N\stackrel{~}{n}_b(1+\stackrel{~}{n}_b))`$ (46)
$`=`$ $`NT\mathrm{ln}(1+\stackrel{~}{n}_b)NqT\mathrm{ln}({\displaystyle \frac{1+\stackrel{~}{n}_b}{\stackrel{~}{n}_b}})+{\displaystyle \frac{T}{2}}\mathrm{ln}(2\pi N\stackrel{~}{n}_b(1+\stackrel{~}{n}_b))`$ (47)
For the usual mean field saddle-point, we obtain $`\stackrel{~}{n}_b=q`$ so the entropy is found to be:

$$S_{mfc}=N(1+q)\mathrm{ln}(1+q)Nq\mathrm{ln}(q)\frac{1}{2}\mathrm{ln}(2\pi Nq(1+q))$$
(48)
which reduces to

$$\frac{2S_{mfc}}{N}=3\mathrm{ln}(\frac{3}{2})+\mathrm{ln}2\frac{1}{N}\mathrm{ln}(\frac{3\pi N}{2})$$
(49)
for the special case $`q=\frac{1}{2}`$, agreeing with the leading approximation to Stirling’s approximation to the exact entropy. For a fully ’self-consistent’ saddle-point the consistency conditions would yield $`\stackrel{~}{n}=\frac{q+\frac{1}{2N}}{1\frac{1}{N}}`$ which is less close to the known exact result.

B. Comparison of Accuracy

1. Susy

At N=2, the susy approximation reverts to either Abrikosov fermions or Schwinger bosons (see Fig 3(b)) recovering the entropy of these approximations if one is willing to remove a divergent term by hand: $`\frac{1}{N}ln(2\pi N(0))`$. It is interesting to speculate that a correction of the form $`\frac{2}{N^2}log(\frac{\stackrel{~}{n}_b(1+\stackrel{~}{n}_b)+\stackrel{~}{n}_f(1\stackrel{~}{n}_f)}{2\pi \stackrel{~}{n}_b(1+\stackrel{~}{n}_b)\stackrel{~}{n}_f(1\stackrel{~}{n}_f)})`$ may be necessary to remove this divergence in these fluctuations and might arise in higher order diagrams, but such a correction is beyond the scope of this paper. Nonetheless, the agreement for N$`>`$2 is quite reasonable (see Fig. 3(a)).

2. Abrikosov fermions

To compare these approximate solutions with the exact solution, by $`N=10`$ a graph is indistinguishable by the eye (see Fig. 4(a)). At $`N=2`$, $`\frac{2S_{exact}}{2}=\mathrm{ln}2`$, so the percentage error is given by

$$\frac{\frac{2S_{Stirling^{}s,mfc}}{2}\frac{2S_{exact}}{2}}{\frac{2S_{exact}}{2}}100=17\%$$
(50)
which agrees with the fully self-consistent version for $`q=\frac{1}{2}`$. So there seems to be no advantage here to fixing $`\frac{F_{fsc}}{\stackrel{~}{n}_f}=0`$ over $`\frac{F_{mft}}{\stackrel{~}{n}_f}=0`$.

3. Schwinger bosons

The self-consistent approximation ($`\frac{F_{fsc}}{n_b}=0`$) does dramatically worse ($`40\%`$ error at N =2) than mean field solution ($`\frac{F_{mft}}{n_b}=0`$). At $`N=2`$,

$$\frac{\frac{2S_{mfc}}{2}\frac{2S_{exact}}{2}}{\frac{2S_{exact}}{2}}100=14\%$$
(51)
whereas

$$\frac{\frac{2S_{Stirling^{}s}}{2}\frac{2S_{exact}}{2}}{\frac{2S_{exact}}{2}}100=17\%$$
(52)
while for large N, both approximations approach the correct result, above $`N=12`$ being indistinguishable (see Fig.4(b)). It is interesting to note that the leading $`\frac{1}{N}`$ approximation to Stirling’s formula does a better job than the formula itself, and the mfc approach reproduces this nice feature at N=2.

3. ATOMIC MODEL

Next, we consider a simple atomic model. In such a model there is an excitation energy for a spin to leave its localized site, creating an empty or vacuum site. The Hamiltonian for such a system is particularly simple, being given by simply the energy difference between the unoccupied site and the occupied spin-ful site. In terms of Hubbard operators, this can be written:

$$H=E_dX_{00}.$$
(53)
A. Entropies, Heat Capacities and Wavefunctions

1. Exact Results

(a) N=2. For a single SU(2) spin with no double occupancy, we know that the exact partition function is given by

$$Z=2+e^{\beta E_d}$$
(54)
which means that the free energy is given by

$$F=T\mathrm{ln}(2+e^{\beta E_d})$$
(55)
the entropy is

$$S=\mathrm{ln}(2+e^{\beta E_d})+\frac{E_de^{\beta E_d}}{T(2+e^{\beta E_d})}$$
(56)
and the heat capacity is given by

$$C_v=\frac{2(\beta E_d)^2e^{\beta E_d}}{(2e^{\beta E}+1)^2}$$
(57)
By writing out the possible states of the system and their energies, we can count the total number of states available per site, as done for the Q = 1, Y = 0, N = 2 case above. This procedure allows us to compute exact results for thermodynamic quantities for larger representations. In Table I, we list the results for Q=2,3 calculated in Appendices B, C.

(b) Exact Wavefunctions. For the simple case (Q=1, N=2) we can choose to represent the spin as either a fermion or a boson, as in reality a single spin is neither. This also applies to the empty state. It can be described by a grassman variable or a spinless boson. This dictates the possible wavefunctions for Q = 1, Y = 0, N = 2 to be:

$$f_{}^{}|0>,f_{}^{}|0>,b_{}^{}|0>,b_{}^{}|0>,\varphi ^{}|0>,\chi ^{}|0>$$
(58)
so we have a double-counting of the available states of the system. We have four states of energy 0, and two of energy -Ed. To agree with physical results, we remove this double counting, so for each partition function written we divide the number of available states to the system by two. That is, for the single box case rather than $`Z=4+2e^{\beta E_d}`$ as one might expect from the degeneracies of the wavefunctions written, we will have $`Z=2+e^{\beta E_d}`$ which agrees with our prior knowledge of this system.

The rules for writing wavefunctions are fairly simple. Since

$$[\theta ^{()},Q]=[\theta ^{()},Y]=0|\nu >\{Q=Q_0,Y=Y_0\}\theta ^{()}|\nu >\{Q=Q_0,Y=Y_0\}$$
(59)
$`[\theta ^{()},X_{\sigma \sigma ^{}}]=0`$ $``$ $`[{\displaystyle \frac{1}{Q}}\{\theta \theta ^{}\theta ^{}\theta \},X_{\sigma \sigma ^{}}]=0`$ (60)
$``$ $`|\nu >\{Q=Q_0,Y=Y_0\}X_{\sigma \sigma ^{}}|\nu >\{Q=Q_0,Y=Y_0\}`$ (61)
$`\{\theta ^{()},X_{0\sigma }\}=0`$ $``$ $`[{\displaystyle \frac{1}{Q}}\{\theta \theta ^{}\theta ^{}\theta \},X_{0\sigma }]=0`$ (62)
$``$ $`|\nu >\{Q=Q_0,Y=Y_0\}X_{0\sigma }|\nu >\{Q=Q_0,Y=Y_0\}`$ (63)
Where we have used the fact that the Hubbard operators conserve particle number and slave representation in the later two statements. This means that by finding a ‘highest weight’ state and applying these operators we can generate a basis of states consistent with $`Q=Q_0,Y=Y_0,N=N_0`$. An additional subtlety to be aware of is that the Hubbard commutation relations mean that for blind application of this method we will end up with a non-orthogonal basis since:

$$[X_{},X_0]=X_0$$
(64)
so we will intentionally discard one of the three to form a basis.

(c) Exact results at general N. To calculate the general partition function for the atomic model, we need to simply count the number of possible states of the system at each energy. This can be simply performed by writing the corresponding Young Tableau for each energy level and using the counting shown above in Eq. (8) (see Fig. 5). As one can only have one slave fermion (as it is a grassman variable), but up to h slave bosons, an L-shaped tableau of energy 0 and dimensions (h,w) allows two L-shaped energy $`E_d`$ states to arise after application of the Hubbard operators with corresponding dimensions (h, w-1) and (h-1, w). Application of Eq. (8) yields the coefficient of the $`e^{\beta E_d}`$ term. At energy $`E_d`$ we have Young Tableaux with dimension (h-1,w-1) and (h-2,w), and one can continue this procedure until energy $`hE_d`$ when only one diagram (1,w-1) remains. This procedure allows us to generate exact results for the free energy in cases that would otherwise prove unaccessible to us.

Using the above algorithm, it is simple to generate exact results for the free energy of an arbitrarily shaped Young Tableau (see Table II).

2. Mean Field + Gaussian corrections

(a) Susy. For the full atomic model, the (mean field) Hamiltonian reads

$`H=E_d(\widehat{n}_\varphi +\widehat{n}_\chi )+\lambda (\widehat{n}_b+\widehat{n}_f+\widehat{n}_\chi +\widehat{n}_\varphi Q_0)+\zeta (\widehat{n}_f+\widehat{n}_\varphi (\widehat{n}_b+\widehat{n}_\chi )+{\displaystyle \frac{1}{Q_0}}<[\theta ,\theta ^{}]>Y_0)`$ (65)
where $`n_\chi `$ and $`n_\varphi `$ are slave fermion and slave boson fields respectively. As we saw in the free spin limit, it is useful to express Q0 and Y0 in terms of the height and width of the Young tableau (h,w), which slightly modifies the constraint on Q0, as we absorb a term $`\frac{\lambda }{N}`$ into the gaussian fluctuations. Then the free energy is

$$F=F_{mft}+\frac{1}{N}F_{g.c}$$
(66)
where

$`F_{mft}`$ $`=`$ $`T\mathrm{ln}(1+e^{\beta (\lambda _bE_d)})NT\mathrm{ln}(1+e^{\beta \lambda _f})+T\mathrm{ln}(1e^{\beta (\lambda _fE_d)})+NT\mathrm{ln}(1e^{\beta \lambda _b})+T\mathrm{ln}(e^{\beta \zeta }+e^{\beta \zeta })`$ (68)
$`\lambda (h+w1)\zeta (hw)`$
$`=`$ $`NT(\stackrel{~}{h}\mathrm{ln}(\stackrel{~}{n}_f)+(1\stackrel{~}{h})\mathrm{ln}(1\stackrel{~}{n}_f)+\stackrel{~}{w}\mathrm{ln}(\stackrel{~}{n}_b)(1+\stackrel{~}{w})\mathrm{ln}(1+\stackrel{~}{n}_b))`$ (70)
$`+T\mathrm{ln}(1n_\chi )T\mathrm{ln}(1+n_\varphi )T\mathrm{ln}({\displaystyle \frac{\stackrel{~}{n}_b\stackrel{~}{n}_f}{\stackrel{~}{n}_b+\stackrel{~}{n}_f}})`$
and

$`{\displaystyle \frac{1}{N}}F_{g.c}`$ $`=`$ $`F_{\delta \lambda _b}+F_{\delta \lambda _f}`$ (71)
$`=`$ $`+{\displaystyle \frac{T}{2}}\mathrm{ln}((2\pi )^2((N\stackrel{~}{n}_b(1+\stackrel{~}{n}_b)+n_\chi (1n_\chi )n_\alpha (1n_\alpha ))(N\stackrel{~}{n}_f(1\stackrel{~}{n}_f)+n_\varphi (1+n_\varphi )n_\alpha (1n_\alpha ))`$ (73)
$`n_\alpha ^2(1n_\alpha )^2))`$
where $`\stackrel{~}{n}_b=\frac{1}{e^{\beta \lambda _b}1}`$ and $`\stackrel{~}{n}_f=\frac{1}{e^{\beta \lambda _f}+1}`$. We have again dropped the $`T\mathrm{ln}2`$ from the double counting of the corner box and the familiar $`T\mathrm{ln}(\frac{\stackrel{~}{n}_b\stackrel{~}{n}_f}{\stackrel{~}{n}_b+\stackrel{~}{n}_f})`$ has again been included in the mean field to get us as close as possible to the correct saddle-point. The saddle-point solution $`\frac{F_{mft}}{\lambda _f}=\frac{F_{mft}}{\lambda _b}=0`$ now yields:

$`\stackrel{~}{n}_f+{\displaystyle \frac{1}{N}}n_\varphi +{\displaystyle \frac{1}{N}}(1n_\alpha )=\stackrel{~}{h}`$ (74)
$`\stackrel{~}{n}_b+{\displaystyle \frac{1}{N}}n_\chi +{\displaystyle \frac{1}{N}}n_\alpha =\stackrel{~}{w}`$ (75)
where $`n_\alpha =\frac{1}{e^{2\beta \zeta }+1}=\frac{1}{\frac{\frac{1}{\stackrel{~}{n}_f}1}{\frac{1}{\stackrel{~}{n}_b}+1}+1}`$ gauges the degree by which a representation is fermionic or bosonic (refer to Appendix D for a discussion), $`n_\varphi =\frac{1}{e^{\beta E_d}(\frac{1}{\stackrel{~}{n}_f}1)1}`$ and $`n_\chi =\frac{1}{e^{\beta E_d}(\frac{1}{\stackrel{~}{n}_b}+1)+1}`$. Hence as in the free spin limit, we include as much as possible in the mean field saddle-point before treating the gaussian fluctuations.

As in the free spin case, the saddle-point has an analytic solution for columns or rows. The special case Q =1, N=2 is of interest despite the non-physical divergence of the entropy noted in the free spin limit (Fig.3). The heat capacity does not suffer from this malady, and is accessible to us as we can analytically (see Table III) remove the divergent terms in the entropy (which are independent of temperature–at least in the slave boson/slave fermion limits) by hand. Furthermore, the phase diagram generated by the saddle point solutions is interesting in its own right, with three possible solutions arising: slave bosons, slave fermions or mixed solutions. In Fig.6 we sketch the simple picture arising. By minimizing the free energy (see Fig. 6(a)), one is driven to the cross-over from a slave fermion description to a slave boson description as one increases the value of $`E_d`$. This appears to be the physically realized situation. A more compelling phase diagram corresponds to maximizing the free energy of the saddle-point solutions, yielding a localized spin slave boson description, a small window of mixed solution followed by a slave fermion description at high temperatures which looks superficially similar to that seen close to quantum critical points (see Fig. 6(b)). Unfortunately it is not as simple to remove the non-physical divergence from this term, and it appears that simple procedures would lead to a non-physical heat capacity, such that the mixed solution likely is not chosen in the atomic model. It is interesting to note that at the mean field level the entropies match at these phase boundaries, and further that the difference between the entropy difference calculated from these approximations at the level of gaussian fluctuations would lead to a weakly positive heat capacity, leading one to believe that if the gaussian fluctuations were properly controlled in the supersymmetric approximation, one might in fact see a second order transition. Unfortunately as yet we have not discovered a way to realize this goal. It would be interesting however to imagine turning on a hybridization to a conduction sea and investigating if beyond some hybridization this state becomes realizable. For general Young Tableau of the row type we see that an analytic slave fermion solution exists, the column type has an analytic slave boson solution arising from the mean field consistency conditions. For more general L-shapes a numerical root-finding is necessary.

(b) Slave fermions. For slave fermions, the Hamiltonian is

$$H_{sf}=E_d\widehat{n}_\chi +\lambda ((\widehat{n}_\chi +\widehat{n}_b)Q_0)$$
(76)
such that the mean field free energy is:

$$F_{mft}=T\mathrm{ln}(1+e^{\beta (\lambda E_d)})+NTq\mathrm{ln}(\stackrel{~}{n}_b)NT(1+q)\mathrm{ln}(1+\stackrel{~}{n}_b)$$
(77)
and the gaussian corrections to this are

$$\frac{1}{N}F_{g.c}=\frac{T}{2}\mathrm{ln}(2\pi (N\stackrel{~}{n}_b(1+\stackrel{~}{n}_b)+n_\chi (1n_\chi )))$$
(78)
where $`n_\chi =\frac{1}{e^{\beta E_d}(1+\frac{1}{\stackrel{~}{n}_b})+1}`$ is the number of slave fermions, and the particles have been constrained such that

$$\stackrel{~}{n}_b+\frac{1}{N}n_\chi =q$$
(79)
For this case, we see that as $`T0`$ for $`E_d<0`$, $`n_\chi 0`$ and $`\stackrel{~}{n}_bq`$, while for $`E_d>0`$, we have $`n_\chi 1`$ and $`\stackrel{~}{n}_bq\frac{1}{N}n_\chi `$ (see Appendix D).

(c) Slave bosons. For slave bosons, the Hamiltonian is given by

$$H_{sb}=E_d\widehat{n}_\varphi +\lambda (\widehat{n}_f+\widehat{n}_\varphi Q_0)$$
(80)
which means that the mean field free energy is

$$F_{mft}=T\mathrm{ln}(1e^{\beta (\lambda E_d)})+NT(1q)\mathrm{ln}(1\stackrel{~}{n}_f)+NTq\mathrm{ln}(\stackrel{~}{n}_f)$$
(81)
and the gaussian corrections are

$$\frac{1}{N}F_{g.c}=\frac{T}{2}\mathrm{ln}(2\pi (N\stackrel{~}{n}_f(1\stackrel{~}{n}_f)+n_\varphi (1+n_\varphi )))$$
(82)
where $`n_\varphi =\frac{1}{e^{\beta E_d}(\frac{1}{\stackrel{~}{n}_f}1)1}`$ and the constraint maintained is:

$$\stackrel{~}{n}_f+\frac{1}{N}n_\varphi =q$$
(83)
For $`E_d>0`$, as $`T0`$, $`\stackrel{~}{n}_f0`$ and $`\frac{1}{N}n_\varphi q`$, while for $`E_d<0`$, $`n_\varphi 0`$ and $`\stackrel{~}{n}_fq`$. Despite the complete transition to the bosonic vacuum state for $`E<0`$, at low T some entropy does remain in the system. Both entropies remain positive in this approximation, while the entropy of the system with $`E>0`$ is much larger.

B. Comparison of accuracy

1. N=2, q = $`\frac{1}{2}`$

As these are all large N techniques, we note that this is an especially cruel limit. Nonetheless, a comparison shows that the heat capacity of slave bosons does the best job for $`E_d>0`$ (empty state favored) while the heat capacity for slave fermions does the best job for $`E_d<0`$ (spin state favored). Choosing the lowest free energy saddlepoint corresponding to Fig 6 (a), the susy approximation agrees with the slave boson result for $`E_d>0`$ and slave fermion for $`E_d<0`$ at low temperatures (see Fig. 7).

2. The Y=0 series – results and errors

In Fig. 8 we compare the heat capacity of our approximation with the exact result. We see that as we increase N while holding $`\stackrel{~}{h}=\frac{1}{2}`$ and Y=0 constant the fit improves, although it is surprisingly good already at N=4. This is partly because the non-physical divergence of the susy approximation as $`N2`$ does not affect the heat capacity as noted in section A2 (a).

In Fig. 9 we see that the entropy being an extrinsic quantity improves dramatically in its description from N=3 to N=20 for the same Y=0 series.

But the eye is somewhat deceiving. In Fig. 10 and 11 we show the absolute errors the entropy and heat capacity respectively as a function of temperature. Both results seem to indicate that there is an error inherent to the approximation of the order of $`\frac{1}{N}`$, which may also be present for the slave boson and slave fermion approximations and should be investigated.

4. CONCLUSIONS

We have presented a simple strongly correlated electron model to test three representations of Hubbard operators. For the limiting case of spin fluctuations, we have demonstrated that mean field theory plus gaussian corrections is capable of capturing the leading order O($`\frac{1}{N}`$) contributions for our model by comparison with Stirling’s approximation expansions of exact results. Also in this case, by considering a mean field treatment plus gaussian corrections we have shown that a Schwinger boson description produces the most accurate results of the three large N methods. By introducing charge fluctuations, we have shown that this method of comparison can be extended to calculate heat capacities, and have found that spin states within this model generally are better described by bosonic descriptions (slave fermion) while empty favored by fermionic descriptions (slave boson), in the dangerous limit of N=2. The lowest energy saddle-point solution arising from the susy representation automatically makes this choice at the mean field level, but encounters a non-physical divergence in its gaussian fluctuations at N=2. The highest energy saddle-point of the susy representation on the other hand exhibits a qualitative similarity to that seen experimentally by systems near an antiferromagnetic quantum critical point, and it will be interesting to see if adding interactions to this simple model might drive us towards the realization of this type of phase diagram which might be expected to have two distinct bosonic order parameters: charge $`\varphi _0^2`$ which might couple to the hybridization and spin $`<b>`$ which might allow the realization of a magnetically ordered state. Furthermore, we have shown that these results can be extended to higher N by the application of simple rules. We have seen that common to these approaches, the results become exact in the limit of large N, implying that if magnetism and paramagnetism co-exist in the ground-state of a material, at large N one might be justified in using the supersymmetric operators.

ACKNOWLEDGEMENTS

Much of the early work on this model was carried out in collaboration with C. Pépin. In addition we gratefully acknowledge discussions on this work with C. Pépin, A. Katanin and O. Parcollet who pointed us towards ref. (22). This work was supported by National Science Foundation grant DMR 9983156.

APPENDIX A: GAUSSIAN APPROXIMATION TO THE ENTROPY OF FREE SPINS

The first term of Eqn (7) corresponds to the mean field theory in our approximate large N treatment. Since we cannot enforce the constraints $`Q=Q_0`$ and $`Y=Y_0`$ rigorously at the mean field level, (Y has a term of order $`\frac{1}{N}`$), we would like to see how inclusion of these terms consistently would affect the saddle point approximation we have made. The total free energy then is made up of five terms to order $`\frac{1}{N}`$,

$$F_{tot}=NF_f+NF_b+F_\eta \lambda Q_0\zeta Y_0+F_{\delta \lambda _f}+F_{\delta \lambda _b}=F_{mft}+F_{\frac{1}{N}g.c.}$$
(A1)
where the two terms of order N define the location of the large N saddle point, and the two latter terms are evaluated at the location of the mean field saddle point. The third term is the expectation value of the fermionic gauge field $`\theta `$ which we calculate at the level of gaussian fluctuations and then include in the saddle point, while fourth and fifth terms are corrections to the fermionic and bosonic free energies respectively. The first and second terms are simply the free fermionic and bosonic free energies:

$$NF_f=NT\mathrm{ln}(1+e^{\beta \lambda _f})NF_b=NT\mathrm{ln}(1e^{\beta \lambda _b})$$
(A2)
The third term is the constraint term that we can not strictly handle at the mean field level, $`\frac{\zeta }{Q}<[\theta ,\theta ^{}]>=\frac{2\zeta }{Q}\theta ^{}\theta +\zeta `$. To accomplish this, we perform a Hubbard-Stratonovich transformation in the following way:

$$\frac{2\zeta }{Q_0}(\theta ^{}+\frac{Q_0}{2\zeta }\overline{\eta })(\theta +\frac{Q_0}{2\zeta }\eta )\frac{2\zeta }{Q}\theta ^{}\theta =\frac{Q_0}{2\zeta }\overline{\eta }\eta +\theta ^{}\eta +\overline{\eta }\theta $$
(A3)
Since with the introduction of this grassman field we no longer have a quartic term in fields, we can treat these fluctuations at the gaussian level. This corresponds to calculating

$$\mathrm{\Pi }_\eta =\frac{Q_0}{2\zeta }+\text{}=\frac{Q_0}{2\zeta }+\varphi (\omega )$$
(A4)
where

$`\varphi (i\omega _n)=T{\displaystyle \underset{\nu }{}}G_f(\omega +\nu )G_b(\nu )={\displaystyle \frac{dzn(z)}{2\pi i}\frac{1}{z+i\omega _n+\lambda _f}\frac{1}{z\lambda _b}}`$ (A5)
$`={\displaystyle \frac{n(i\omega _n+\lambda _f)}{i\omega _n+2\zeta }}+{\displaystyle \frac{n(\lambda _b)}{i\omega _n2\zeta }}={\displaystyle \frac{f(\lambda _f)+n(\lambda _b)}{\omega 2\zeta }}={\displaystyle \frac{Q_0}{\omega 2\zeta }}`$ (A6)
Which means that

$$\mathrm{\Pi }_\eta =\frac{Q_0}{2\zeta }+\frac{Q_0}{\omega 2\zeta }=\frac{Q_0\omega }{2\zeta (\omega 2\zeta )}$$
(A7)
So the contribution to the free energy from these fluctuations is given by

$`F_\eta ={\displaystyle \frac{1}{\beta }}{\displaystyle \underset{i\omega _n}{}}\mathrm{ln}(\mathrm{\Pi }(i\omega _n))={\displaystyle \frac{1}{\beta }}{\displaystyle \underset{i\omega _n}{}}\mathrm{ln}({\displaystyle \frac{Q_0}{2\zeta }})+{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{i\omega _n}{}}\mathrm{ln}(i\omega _n2\zeta ){\displaystyle \frac{1}{\beta }}{\displaystyle \underset{i\omega _n}{}}\mathrm{ln}(i\omega _n)`$ (A8)
$`=0{\displaystyle \frac{dzf(z)}{2\pi i}\mathrm{ln}(z2\zeta )}+{\displaystyle \frac{dzf(z)}{2\pi i}\mathrm{ln}(z+0)}`$ (A9)
$`={\displaystyle \frac{1}{\beta }}\mathrm{ln}(1+e^{2\beta \zeta }{\displaystyle \frac{1}{\beta }}\mathrm{ln}(2))={\displaystyle \frac{1}{\beta }}\mathrm{ln}({\displaystyle \frac{1+e^{2\beta \zeta }}{2}})`$ (A10)
Now, adding back in the additional $`\zeta `$ from the constraint (as $`e^{\mathrm{ln}(\beta \zeta )}`$) which was not done for the susy Kondo model, this becomes $`F_\eta =T\mathrm{ln}(\frac{e^{\beta \zeta }+e^{\beta \zeta }}{2})`$. At the mean field level the constraint term reads $`\lambda Q=\lambda (n_b+n_f1)`$ which implies an extra contribution of $`\lambda `$ to the free energy hence this fluctuation term becomes $`F_\eta =T\mathrm{ln}(\frac{e^{\beta \lambda _f}+e^{\beta \lambda _b}}{2})=T\mathrm{ln}(\frac{\stackrel{~}{n}_b+\stackrel{~}{n}_f}{2\stackrel{~}{n}_b\stackrel{~}{n}_f})`$.

The corrections coming from varying the Lagrange multiplier about its mean field value and from the new gauge field arise from an expansion of the Lagrangian as $`Z=𝒟(f_\sigma ,b_\sigma ,\delta \lambda _f,\delta \lambda _b,\eta )e^{{\scriptscriptstyle 𝑑\tau _0}}(1+\delta +\frac{1}{2}(\delta )^2)`$ where $`L_0`$ represents the mean field Lagrangian, and we perform a saddle point evaluation. For a grassman gauge field, the expansion terminates at second order and the first order terms vanish so we have therefore $`e^{(\delta )^2}`$, a grassman gaussian integral. For the Lagrange multiplier corrections, the first non-zero term is the second order diagram and in the gaussian approximation, we treat this as an exponential correction analogous to the fermionic case. This being the case, the corrections are as follows. The diagram of interest for $`F_{\delta \lambda _f}`$ is :

(A11)
This is given by:

$`T{\displaystyle \underset{i\omega _n\alpha =1..N}{}}G_f(i\omega _n)G_f(i\omega _n)={\displaystyle \underset{\alpha }{}}{\displaystyle \frac{dz}{2\pi i}\frac{}{z}f(z)\frac{1}{z\lambda _f}}`$ (A12)
$`=N(\beta \stackrel{~}{n}_f(\lambda _f)(1\stackrel{~}{n}_f(\lambda _f)))`$ (A13)
For a real gaussian variable, we have $`Z=\frac{d(\delta \lambda _f\beta )}{\pi }e^{N_o^\beta 𝑑\tau \beta \stackrel{~}{n}_f(\lambda _f)(1\stackrel{~}{n}_f(\lambda _f))(\delta \lambda _f)^2}=\frac{1}{\pi }\sqrt{\frac{\pi }{2N\stackrel{~}{n}_f(\lambda _f)(1\stackrel{~}{n}_f(\lambda _f))}}`$ where the integration is a sum over imaginary time (so it doesn’t affect $`\stackrel{~}{n}_f`$), which totals $`\beta `$, the inverse temperature. The integration variable is normalized by $`\beta `$ to make it dimensionless. The corresponding free energy contribution from these fluctuations is then:

$$F_{\delta \lambda _f}=T\mathrm{ln}Z=\frac{T}{2}(\mathrm{ln}(2\pi N\stackrel{~}{n}_f(1\stackrel{~}{n}_f)))$$
(A14)
which is familiar as the gaussian fluctuations about the mean field for Abrikosov fermions. The diagram of interest for $`F_{\delta \lambda _b}`$ is:

(A15)
Which is given by:

$`T{\displaystyle \underset{i\nu _n\alpha =1..N}{}}G_b(i\nu _n)G_b(i\nu _n)={\displaystyle \underset{\alpha =1..N}{}}{\displaystyle \frac{dz}{2\pi i}n(z)\frac{1}{(z\lambda _b)^2}}`$ (A16)
$`={\displaystyle \underset{\alpha =1..N}{}}{\displaystyle \frac{dz}{2\pi i}\frac{n(z)}{z}\frac{1}{z\lambda _b}}=N\beta (\stackrel{~}{n}_b(\lambda _b)(1+\stackrel{~}{n}_b(\lambda _b)))`$ (A17)
Following the same procedure to find the partition function factor, we get $`Z=\frac{1}{\pi }\sqrt{\frac{\pi }{2N\stackrel{~}{n}_b(\lambda _b)(1+\stackrel{~}{n}_b(\lambda _b))}}`$ so the free energy contribution from these fluctuations is:

$$F_{\delta \lambda _b}=\frac{T}{2}(\mathrm{ln}(2\pi N\stackrel{~}{n}_b(1+\stackrel{~}{n}_b)))$$
(A18)
which is familiar as the gaussian fluctuations about the mean field for Schwinger bosons. Enforcing the constraint on Y at the mean field level, couples these two contributions as additional terms arise as $`\frac{^2F_\eta }{\lambda _f\lambda _b}`$=$`\frac{^2F_\eta }{\lambda _b\lambda _f}`$=-$`\frac{^2F_\eta }{\lambda _b^2}`$=-$`\frac{^2F_\eta }{\lambda _f^2}`$=-$`\beta n_\alpha (1n_\alpha )`$ leading to the combined contribution to the gaussian fluctuations as

$$F_{\delta \lambda _f,\delta \lambda _b}=\frac{T}{2}\mathrm{ln}((2\pi )^2Det\left[\begin{array}{cc}N\stackrel{~}{n}_f(1\stackrel{~}{n}_f)n_\alpha (1n_\alpha )& n_\alpha (1n_\alpha )\\ n_\alpha (1n_\alpha )& N\stackrel{~}{n}_b(1+\stackrel{~}{n}_b)n_\alpha (1n_\alpha )\end{array}\right])$$
(A19)
The generalization of this term to the atomic model is straightforward as one incorporates fluctuations in the number of slave partners as well.

APPENDIX B: Q=2, N=2 ATOMIC WAVEFUNCTIONS

For Q = 2, N = 2, Y = 1, we have a totally antisymmetric wavefunction, so to write wavefunctions for these states we start with a simple two-particle anti-symmetric wavefunction and operate with lowering operators ($`X_{0\sigma }`$, $`\mathrm{\Theta }`$). The spin-spin Hubbard operator does not create new states because the initial wavefunctions are totally antisymmetric. We have with no slaves :

$$|\widehat{A}>=f_{}^{}f_{}^{},\mathrm{\Theta }|\widehat{A}>$$
(B1)
whose wavefunctions associated are:

$$f_{}^{}f_{}^{}|0>,b_{}^{}f_{}^{}|0>b_{}^{}f_{}^{}|0>$$
(B2)
With one vacuum state we have four possibilities:

$$X_0|\widehat{A}>,X_0|\widehat{A}>,\mathrm{\Theta }X_0|\widehat{A}>,\mathrm{\Theta }X_0|\widehat{A}>$$
(B3)
$$f_{}^{}\varphi ^{}|0>,f_{}^{}\varphi ^{}|0>,(b_{}^{}\varphi ^{}\chi ^{}f_{}^{})|0>,(b_{}^{}\varphi ^{}\chi ^{}f_{}^{})|0>$$
(B4)
While with two empty states,

$$X_0X_0|\widehat{A}>,\mathrm{\Theta }X_0X_0|\widehat{A}>$$
(B5)
are the remaining possibilities:

$$(\varphi ^{})^2|0>,\varphi ^{}\chi ^{}|0>$$
(B6)
For Q = 2, Y = -1, N = 2, we have a completely symmetric wavefunction, so we naturally start with a bosonic description as the highest weight state and apply boson lowering operators ($`X_{0\sigma },\mathrm{\Theta }^{}`$) and spin exchange operators ($`X_{\sigma \sigma }`$) to generate the complete set. The fully particle description is comprised of

$$|B>=b_{}^{}b_{}^{}|0>,X_{}|B>,X_{}X_{}|B>,\mathrm{\Theta }^{}|B>,\mathrm{\Theta }^{}X_{}|B>,\mathrm{\Theta }^{}X_{}X_{}|B>$$
(B7)
which in terms of fermionic and bosonic operators appears as:

$$b_{}^{}b_{}^{}|0>,b_{}^{}b_{}^{}|0>,b_{}^{}b_{}^{}|0>,f_{}^{}b_{}^{}|0>,f_{}^{}b_{}^{}|0>+f_{}^{}b_{}^{}|0>,f_{}^{}b_{}^{}|0>$$
(B8)
Only fluctuations to one empty site are possible ($`X_{0\sigma }^2=0`$)within the constraints for these quantum numbers. Thus the other possible states are categorized by:

$$X_0|B>,X_0|B>,X_0\mathrm{\Theta }^{}|B>,X_0\mathrm{\Theta }^{}X_{}X_{}|B>$$
(B9)
$$b_{}^{}\chi ^{}|0>,b_{}^{}\chi ^{}|0>,\varphi ^{}b_{}^{}|0>+\chi ^{}f_{}^{}|0>,\varphi ^{}b_{}^{}|0>+\chi ^{}f_{}^{}|0>$$
(B10)
APPENDIX C: Q=3, N=2 ATOMIC WAVEFUNCTIONS

For the (mixed) atomic model we can treat Q = 3 exactly (without comparison) for all three possible cases (Y = -2,0,2). Here for comparison’s sake we list the wavefunctions of the N = 2 result, so that we can compare with a fairly simple calculation where the corner particle gets to choose its type in a non-trivial way. Since the degeneracy of the level gives us a value for the partition function, we can easily calculate a free energy and a heat capacity for this atomic model, and compare these results to the algorithm of section II i) d). For $`Q=3,Y=0,N=2`$, our basis consists of four states with no slave particles, $`8=2^3`$ distinct states with one slave particle and four states with two slave particles. For completeness I will list them here, and how they were derived. With no slave particles we have:

$$|A>=b_{}^{}f_{}^{}f_{}^{}|0>,\mathrm{\Theta }|A>,X_{}|A>,\mathrm{\Theta }X_{}|A>$$
(C1)
which are respectively:

$$b_{}^{}f_{}^{}f_{}^{}|0>,b_{}^{}(b_{}^{}f_{}^{}|0>b_{}^{}f_{}^{}|0>),b_{}^{}f_{}^{}f_{}^{}|0>,b_{}^{}(b_{}^{}f_{}^{}|0>b_{}^{}f_{}^{}|0>).$$
(C2)
With one slave particle:

$`X_0`$ $`|A>,X_0|A>,\mathrm{\Theta }X_0|A>,\mathrm{\Theta }X_0|A>,`$ (C3)
$`X_0`$ $`X_{}|A>,\mathrm{\Theta }X_0X_{}|A>,X_0X_{}|A>,\mathrm{\Theta }X_0X_{}|A>`$ (C4)
which are:

$`\varphi ^{}b_{}^{}f_{}^{},\chi ^{}f_{}^{}f_{}^{}|0>\varphi ^{}b_{}^{}f_{}^{}|0>,b_{}^{}(\varphi ^{}b_{}^{}|0>\chi ^{}f_{}^{}|0>),`$ (C5)
$`2b_{}^{}\chi ^{}f_{}^{}|0>b_{}^{}(\chi ^{}f_{}^{}|0>+\varphi ^{}b_{}^{}|0>),\chi ^{}f_{}^{}f_{}^{}|0>+\varphi ^{}b_{}^{}f_{}^{}|0>,`$ (C6)
$`b_{}^{}(\chi ^{}f_{}^{}|0>+\varphi ^{}b_{}^{}|0>)2b_{}^{}\chi ^{}f_{}^{}|0>,\varphi ^{}b_{}^{}f_{}^{}|0>,b_{}^{}(\chi ^{}f_{}^{}|0>\varphi ^{}b_{}^{}|0>).`$ (C7)
Here we have not written two equivalent eigenvectors. While with two slave particles:

$$X_0X_0|A>,\mathrm{\Theta }X_0X_0|A>,X_{}X_0X_0|A>,\mathrm{\Theta }X_{}X_0X_0|A>$$
(C8)
whose wavefunctions are:

$$\varphi ^{}(\chi ^{}f_{}^{}|0>+\varphi ^{}b_{}^{}|0>),b_{}^{}\varphi ^{}\chi ^{}|0>,\varphi ^{}(\chi ^{}f_{}^{}|0>+\varphi ^{}b_{}^{}),b_{}^{}\varphi ^{}\chi ^{}|0>$$
(C9)
We can also treat the fully symmetric case Q = 3, Y = -2, N = 2. As in the case with one less particle, we have spin fluctuations, wavefunction character fluctuations and limited charge fluctuations. The filled particle levels are described by:

$`|C>=(b_{}^{})^3|0>,X_{}|C>,X_{}X_{}|C>,X_{}^3|0>,`$ (C10)
$`\mathrm{\Theta }^{}|C>,\mathrm{\Theta }^{}X_{}|C>,\mathrm{\Theta }^{}X_{}X_{}|C>,\mathrm{\Theta }^{}X_{}^3|0>`$ (C11)
which are:

$`(b_{}^{})^3|0>,(b_{}^{})^2b_{}^{}|0>,(b_{}^{})^2b_{}^{}|0>,(b_{}^{})^3|0>`$ (C12)
$`f_{}^{}(b_{}^{})^2|0>,b_{}^{}(f_{}^{}b_{}^{}|0>+2f_{}^{}b_{}^{}|0>),b_{}^{}(f_{}^{}b_{}^{}|0>+2f_{}^{}b_{}^{}|0>),f_{}^{}(b_{}^{})^2|0>`$ (C13)
and with one less:

$`X_0|C>,X_0X_{}|C>,X_0X_{}^2|C>,`$ (C14)
$`X_0\mathrm{\Theta }^{}|C>,X_0\mathrm{\Theta }^{}X_{}|C>,X_0\mathrm{\Theta }^{}X_{}^2|C>`$ (C15)
$`\chi ^{}(b_{}^{})^2|0>,\chi ^{}b_{}^{}b_{}^{}|0>,\chi ^{}(b_{}^{})^2|0>,`$ (C16)
$`\varphi ^{}(b_{}^{})^2|0>+2\chi ^{}f_{}^{}b_{}^{}|0>,\varphi ^{}b_{}^{}b_{}^{}|0>+\chi ^{}(f_{}^{}b_{}^{}|0>+f_{}^{}b_{}^{}|0>),\varphi ^{}(b_{}^{})^2|0>+2\chi ^{}f_{}^{}b_{}^{}|0>`$ (C17)
And of course we can treat the fully antisymmetric representation (Q=3,Y=2,N=2). Note that since N = 2 we can no longer have 3 spins! Here’s what we have (totally unoccupied):

$$|D>=(\varphi ^{})^3|0>,\mathrm{\Theta }|D>$$
(C18)
which are

$$(\varphi ^{})^3|0>,\chi ^{}(\varphi ^{})^2|0>$$
(C19)
one spin

$$X_0|D>,X_0|D>,X_0\mathrm{\Theta }|D>,X_0\mathrm{\Theta }|D>$$
(C20)
which look like

$$f_{}^{}(\varphi ^{})^2|0>,f_{}^{}(\varphi ^{})^2|0>,b_{}^{}(\varphi ^{})^2|0>+2f_{}^{}\chi ^{}\varphi ^{}|0>,b_{}^{}(\varphi ^{})^2|0>+2f_{}^{}\chi ^{}\varphi ^{}|0>$$
(C21)
two spins:

$$X_0X_0|D>,X_0X_0\mathrm{\Theta }|D>$$
(C22)
and these are:

$$f_{}^{}f_{}^{}\varphi ^{}|0>,\varphi ^{}(f_{}^{}b_{}^{}|0>f_{}^{}b_{}^{}|0>)+\chi ^{}f_{}^{}f_{}^{}|0>$$
(C23)
APPENDIX D: MEAN FIELD OCCUPATION

To gain an appreciation of how the character of the spin system changes as a function of temperature, in Fig. 12 we show the mean field results of the slave fermion approximation for Q = 1, N=2. We see (Fig. 12 (a)) that in this approximation when $`E_d>0`$ the Hubbard operators describe a charged hole which obeys grassman statistics. As the temperature is raised the bosonic description of the spin begins to carry more weight. When $`E_d<0`$, the ground state at T=0 is a pure spin state which acquires a bosonic description (Schwinger bosons) (Fig. 12 (b)).

In Fig. 13 we show the mean field results of the slave boson approximation for Q=1, N=2. In Fig. 13 (a) we see that for $`E_d>0`$ the zero temperature ground state of the system is a positively charged bosonic hole. As the temperature is increased some fermionic spinorial character of the electron returns. When $`E_d<0`$ the system enters a fermionic ground state at low temperatures which is just the free spin case treated in section 2 (Abrikosov fermions). (see Fig. 13 (b)).

For the case of physical interest (Q=1,N=2,Y=0), we saw that the supersymmetric description has three possible physical analytic solutions (see Table III). When $`n_\alpha =0`$ one recovers the slave fermion mean field results (Fig. 12), when $`n_\alpha =1`$ one recovers the slave boson results (Fig. 13), and a mixed solution exists in the small region $`E_d=0`$ to $`E_d=T\mathrm{ln}(2)`$ as shown in the inset to Fig. 14. If one follows the Y=0 series by increasing N while holding $`\stackrel{~}{h}=\frac{h}{N}=\frac{1}{2}`$ constant (see Fig.14), one sees that one generically has a mixed solution – neither slave bosons nor slave fermions can describe the physics of systems exhibiting this symmetry. As approximated methods become asymptotically exact as N $`\mathrm{}`$, if the phase diagram of Fig. 6 (b) were to become the physically realized state, one might be able to justify a large N expansion to treat the physics of the mixed region using susy Hubbard operators at Y=0.


text:
# Including Systematic Uncertainties in Confidence Interval Construction for Poisson Statistics

## I  Introduction

A limit on , or a measurement of , a physical quantity at a given confidence level is usually set by comparing a number of detected events, $`n_o`$, with the number of expected events from the known background sources contributing to the physical process in question, $`n_b`$. How ’compatible’ these numbers are determines how much room there is for new processes, ie., for a signal. How well do the number of observed events and expected background compare, strongly depends on the systematic uncertainties present in the measurement. Systematic uncertainties must, therefore, be taken into account in the limit or confidence belt calculation that is finally published.

Traditionally, confidence limits are set using a Neyman construction Neyman:1937a . This is a purely frequentist method. G. Feldman and R. Cousins Feldman:1998a  have proposed an improved method to construct confidence intervals based on likelihood ratios, a method allready known in statistics and originally described in Kendall:91 . Still, this method is based on the original Neyman construction, and needs to be extended to incorporate systematic uncertainties in the measurement. Along this line, a modification of the Neyman method that incorporates systematic uncertainties in the experimental signal efficiency has been proposed by V. Highland & R. Cousins Cousins:1992a . These authors use a “semi”- Bayesian approach where an average over the probability distribution of the experimental sensitivity (and its uncertainty) is performed. By construction, the method is of limited accuracy in the limit of high relative systematic uncertainties.

Recently, an entirely frequentist approach has been proposed for uncertainty in the background rate prediction Rolke:2001a . That approach is based on a two-dimensional confidence belt construction and likelihood ratio hypothesis testing and treats the uncertainty in the background as a statistical uncertainty rather than as a systematic one.

The interest aroused recently in the High Energy Physics community about the many open issues on setting limits and quoting confidence levels is stresseded by the organization of devoted workshops on the subject. We refer the reader to the proceedings of the recent workshops at CERN CERN:00a , FERMILAB FERM:00a  and Durham DURHAM:02  for a review of the status of the field.

In this paper we extend the method of confidence belt construction proposed in Cousins:1992a  to include systematic uncertainties both in the signal and background efficiencies as well as theoretical uncertainties in the background prediction. The proposed method allows as well to use newer ordering schemes. A recent attempt to include systematic uncertainty in the background prediction in a similar manner has been presented in Giunti:1998xv . The paper is organized as follows. In section II we give a short review of the confidence belt construction schemes which we will use. In section III we describe how to include the systematic uncertainties, in section IV we discuss how the confidence belt construction is performed and present some selected results. We compare the results of this method with other methods to include systematics in section V.

We introduce the tests of coverage performed in section VI and present an example based on data from the AMANDA neutrino experiment in section VII.

## II  The Construction of Confidence Intervals

The frequentist construction of confidence intervals is described in detail elsewhere Eadie:1982a . Here we will give just a short review.

Let us consider a Poissonian probability density function (PDF), $`p(n)_{s+b}`$, for a fixed but unknown signal, $`s`$, in the presence of a known background with mean $`b`$. For every value of $`s`$ we can find two values $`n_1`$ and $`n_2`$ such that

$$\underset{n^{}=n_1}{\overset{n_2}{}}p(n^{})_{s+b}=1\alpha $$
(1)
where $`1\alpha `$ denotes the confidence level (usually quoted as a 100(1-$`\alpha `$)% confidence interval). Since we assume a Poisson distribution, the equality will generally not be fulfilled exactly. A set of intervals $`[n_1(s+b,\alpha ),n_2(s+b,\alpha )]`$ is called a confidence belt. Graphically, upon a measurement, $`n_o`$, the confidence interval $`[s_1,s_2]`$ is determined by the intersection of the vertical line drawn from the measured value $`n_o`$ and the boundary of the confidence belt. This is illustrated in figure 1. The probability that the confidence interval will contain the true value $`s`$ is $`1\alpha `$, since this is true for all $`s`$ per construction.

The choice of the $`n_1`$ and $`n_2`$ is, however, not unique to define the confidence belt. An additional criterion has to be applied. The choices originally proposed by Neyman Neyman:1937a  are

$$\underset{n^{}=0}{\overset{n_1}{}}p(n^{})_{s+b}=\underset{n^{}=n_2}{\overset{\mathrm{}}{}}p(n^{})_{s+b}=\frac{1\alpha }{2}$$
(2)
for central confidence intervals, and

$$\underset{n^{}=0}{\overset{n_1}{}}p(n^{})_{s+b}=1\alpha $$
(3)
for upper confidence limits. This method presents certain drawbacks in the case of small samples and, in particular, can yield null results (in the sense that the algorithm gives no answer) in the case when no events have been observed. Also, the decision of quoting a measurement (that is, a central confidence interval) or an upper limit might not be straightforward before performing an experiment.

### II.1 Likelihood ratio ordering

To solve this problem, a modification of the Neyman method has been proposed Kendall:91  Feldman:1998a  that is based on a more rationalized ordering scheme of the elements in the sum in equation (1), based on likelihood ratios. This approach automatically provides central confidence intervals when motivated and upper limits when necessary, therefore it is often denoted as the “unified approach”. Instead of using the choices given in the previous section, the following ordering scheme is applied in solving equation (1):
For each $`n`$ the $`s_{best}`$ is found which maximizes the likelihood $`(n)_{s+b}`$. In case of a simple Poissonian distribution with known background, $`s_{best}`$ is given by $`max(0,nb)`$. Then for a fixed $`s`$ the ratio

$$R(s,n)_{}=\frac{_{s+b}(n)}{_{s_{best}+b}(n)}$$
(4)
is computed for each $`n`$, and all $`n`$’s are consequently ranked according to the value of this ratio. Values of $`n`$ are included in the confidence belt starting with the $`n`$ with the highest rank (largest $`R_{}`$) and then decreasing rank until $`_{n=n_1}^{n_2}p(n)_{s+b}=1\alpha `$. After the confidence belt has been constructed in this way, the confidence interval $`[s_1,s_2]`$ is found as described in the previous section. Note that this ordering principle is a standard method within the theory of likelihood ratio tests Kendall:91 .

This approach has some undesired features as well. There is a background dependence of the upper limit in case of less events observed than expected from background. This can lead to situations where measurements with higher background give a better limit, a clearly undesireable effect. B. Roe & M. Woodroofe Roe:1999a  proposed a solution to this problem which we briefly describe next.

### II.2 Conditioning

A variation of the classical method of constructing confidence belts is to use the fact that, given an observation $`n_o`$, it is known that the background can not have been larger than $`n_o`$ itself. To incorporate this knowledge into the PDF, the authors in Roe:1999a  have proposed the following modification:

$$q_{s+b}^{n_o}(n)=\{\begin{array}{cc}\frac{p\left(n\right)_{s+b}}{{\displaystyle \underset{n^{}=0}{\overset{n_o}{}}}p\left(n^{}\right)_b}\hfill & \text{if }nn_o\hfill \\ & \\ \frac{{\displaystyle \underset{n^{}=0}{\overset{n_o}{}}}p\left(n^{}\right)_bp\left(nn^{}\right)_s}{{\displaystyle \underset{n^{}=0}{\overset{n_o}{}}}p\left(n^{}\right)_b}\hfill & \text{if }n>n_o\hfill \end{array}$$
(5)
The likelihood ratio ordering can then be applied with this new PDF. Note that in this case the PDF is dependent on the number of observed events. This approach solves the background dependence of the upper limit: a limit set when no events are observed stays constant at a value of 2.44 independent of the expected background (which agrees with the result of the original likelihood ordering for no events observed and no expected background). However, this method does not satisfy all the requirements of proper coverage Zech:2001a  and has problems when applied to the case of a Gaussian distribution with boundaries Cousins:2000a . An extension based on a Bayesian approach with tests of coverage can be found in Roe:2001a .

## III  The inclusion of systematic uncertainties

The way of incorporating systematic uncertainties into the confidence belt construction presented in this paper does not affect the particular ordering scheme. Instead, it takes into account the systematic uncertainties by assuming (or if possible determining) a PDF which parameterizes our knowledge about the uncertainties and integrating over this PDF. It has been noted that averaging over systematic uncertainties in itself is a Bayesian approach Cousins:1992a  . Therefore the method presented is refered to as “semi-Bayesian”, combining classical and Bayesian elements. We return to this point in section VI. Usually, uncertainties are assumed to be described by a Gaussian distribution, which we will adopt for the remainder of this paper. The implementation, however, makes it easy to use other parametrizations for the uncertainties.

We will refer in the following to the parameters with systematic uncertainties also as nuisance parameters.

Two examples of how the PDF modifies if systematic uncertainties are present are the following. In the case that the only uncertainty present is a theoretical uncertainty of the background process the PDF is modified to:

$$q(n)_{s+b}=\frac{1}{\sqrt{2\pi }\sigma _b}\underset{0}{\overset{\mathrm{}}{}}p(n)_{s+b^{}}e^{\frac{(bb^{})^2}{2\sigma _b^2}}𝑑b^{}$$
(6)
Here $`b`$ is the estimated background level, and $`\sigma _b`$ is the uncertainty in the background estimation. If, in addition to the theoretical uncertainty for background, there is the need to include the uncertainty in the signal detection efficiency the expression for $`q(n)_{s+b}`$ might be extended to:

$$\begin{array}{cc}q(n)_{s+b}=\frac{1}{2\pi \sigma _b\sigma _ϵ}\times \hfill & \\ _0^{\mathrm{}}_0^{\mathrm{}}p(n)_{b^{}+ϵ^{}s}e^{\frac{(bb^{})^2}{2\sigma _b^2}}e^{\frac{(1ϵ^{})^2}{2\sigma _ϵ^2}}𝑑b^{}𝑑ϵ^{}\hfill & \end{array}$$
(7)
where $`\sigma _ϵ`$ is the uncertainty in the detection efficiency expressed in relative terms with respect to the nominal efficiency. It is important to realize that the integration variables, here $`ϵ^{}`$ and $`b^{}`$, are the possible “true” (but unknown) values of nuisance parameter. This indicates that this method is based on Bayesian statistics.

## IV  POLE: A General Algorithm for Confidence Belt Construction

The integrals (6) and (7) can be solved using different methods. We note, however, that they are examples of simplified cases. The most general experimental situation involves both an uncertainty in signal efficiency as well as in the background detection efficiency, which are usually correlated, and possibly an additional theoretical uncertainty in the background process prediction. We have developed an algorithm that takes these effects into account with the proper correlations between them. The algorithm performs a Monte Carlo integration over the systematic uncertainties. It has been implemented as a FORTRAN program, POLE (POissonian Limit Estimator) Conrad:2001a . In the examples used in this section a Gaussian distribution of uncertainties is assumed, but the algorithm makes it easy to implement PDFs other than Gaussian (see the next section for an example of using a different distribution). For the moment the code supports a Gaussian, flat and log-normal parametrization of the uncertainties. After determining the PDF through evaluation of the integrals, different ordering schemes can be applied for the final calculation of the confidence belt. The results presented here are mainly for the likelihood ratio ordering scheme with and without conditioning. We restrict ourselves to present systematic uncertainties of signal and background efficiencies separately to give a clear idea of the effect of varying a single variable at a time. Real applications usually combine those uncertainties.

The confidence belt constructions have been performed using steps of 0.05 in signal expectation and performing the construction up to a maximal signal expectation of 50 and a maximal number of detected events of 100. Including systematic uncertainties generally leads to a widening of the confidence belt. Figure 2 shows an example of a likelihood ratio confidence belt construction with and without uncertainty in the signal efficiency, where a background expectation $`b=2`$ has been assumed.

Examples for some resulting intervals are given in tables 1 and 2. Different combinations of number of observed events, $`n_0`$, and expected background, $`b`$, are given for different uncertainties in the signal and background efficiency.
The width of the interval for two particular examples of observed events and expected background as function of signal efficiency uncertainty and background uncertainty is shown in figure 3. Note that for low background expectation, the uncertainties in the background can be neglected (see also table 2). Figures 4 and 5 give more extended information on resulting intervals.

An interesting case arises when there are significantly less events observed than expected from background and there is an uncertainty in the signal efficiency. In this case, the width of the confidence interval does not increase (see table 3). Note that if we use conditioning the effect disappears. The same can not be observed in the case where we only consider an increasing background uncertainty (see table 4).

### IV.1 Negative values of the nuisance parameters: Using a log - normal distribution.

In experimental situations where the systematic uncertainties are high, a problem might arise due to the fact that sampling from a Gaussian PDF allows negative values. POLE is dealing with these cases truncating the Gaussian distribution and renormalizing the part above zero.
We examine the effect of truncating the Gauss distribution by calculating the confidence interval for different values of truncation point (see figure 6).

Considering a Gaussian distribution centered on one widh $`\sigma `$ = 40 %, a truncation at zero removes only 0.7 %. Figure 6 therefore indicates that effects on the confidence interval due to the truncation are negligible for all cases considered in this paper.
A PDF for the nuisance parameters extending to negative values or which falls off to zero discontinuously is certainly undesired from a conceptual point of view. We therefore test the behavior of the confidence interval if we replace the Gaussian distribution with a log-normal distribution, which in the general form is given by:

$$q(x)_{\mu ,\sigma }=\frac{1}{\sqrt{2\pi }x\sigma }e^{\frac{(\mathrm{ln}x\mu )^2}{2\sigma ^2}}$$
(8)
We require the mean of the log-normal distribution to be the nominal value of the nuisance parameter and use the Gaussian standard deviation as before (the variance of the log normal distribution will then be approximately the same). The confidence interval for Neyman and likelihood ratio ordering under these assumptions are shown for one particular example of number of observed events and expected background as a function of signal efficiency in figure 7. The differences between using a Gaussian distribution and using a log-normal distribution are generally small, in our example less than $``$ 2 %. The use of a log-normal distribution is implemented as an option in POLE.

## V  Comparison with the $`\chi ^2`$ method.

Since the ratio of the likelihoods is asymptotically $`\chi ^2`$ distributed the approximation:

$$\mathrm{\Delta }\chi ^2=2\times \underset{ϵ}{\mathrm{min}}\frac{(n_o,s,ϵ)}{(n_o,s_{best},ϵ_{best})}$$
(9)
is often used, see e.g. Feldman:2000 .
Here for a given observation,$`n_0`$, the $`\mathrm{\Delta }\chi ^2`$ is calculated as a function of $`s`$ and a cut on $`\mathrm{\Delta }\chi ^2`$ is performed to obtain the confidence interval (e.g. $`\mathrm{\Delta }\chi ^2`$ = 2.71 corresponds to 90 % confidence level for one degree of freedom).

Figure 8 illustrates the effect of including uncertainties on the resulting confidence intervals for the $`\chi ^2`$ approximation as compared to the method proposed here.
Generally, the $`\chi ^2`$ approximation gives more conservative results than the pole method. Since - as we will see in the following sections - using the pole method leads to some over-coverage, this is clearly undesirable.

## VI  Tests of Coverage

From a frequentist point of view, an algorithm is said to have the correct coverage if, given a confidence level $`1\alpha `$ and a large number of repeated identical experiments, it provides correct answers in a fraction $`1\alpha `$ of the cases, independent of the value of $`s`$. To test the coverage of the algorithm proposed in this paper, we perform the construction described in the previous sections for a large number of simulated experiments, where we predefine the true signal and background and then determine $`n_o`$ by random sampling from a Poisson distribution. We then calculate how often the obtained confidence interval does not contain the predefined $`s`$. We define the coverage ratio:

$$R=\frac{n_{false}}{n_{tot}}.$$
(10)
Here, $`n_{false}`$ denotes the number of simulated experiments in which the result of the algorithm does not contain the predefined $`s`$, and $`n_{tot}`$ denotes the number of simulated experiments performed. If we choose $`1\alpha `$ to be 0.9, perfect coverage would mean $`R=0.1`$, independent of signal expectation assumption.
A value or $`R`$ smaller than $`0.1`$ means that the method over-covers. Expected coverage was studied mostly in the context of Bayesian intervals, small number of events or including conditioning Narsky:2000  Mandelkern:2001  Helene:99  and without taking into account systematic uncertainties. Very recently, a study was presented considering coverage including systematic uncertainties Conrad:2002ur .

In the next sections, besides presenting coverage tests done for higher signal expectations without uncertainties, we will focus on the coverage of the methods if systematic uncertainties are included.

### VI.1 Coverage without systematic uncertainties

We show an example of a plot of the coverage ratio (here using steps of 0.1 in signal space) in figure 9 for Neyman and likelihood ratio ordering. Both methods seem to over-cover for almost all cases (which is expected because of the discreteness of the Poisson distribution). There is no signal expectation dependence of the coverage ratio except for the “see-saw” behavior which again reflects the discreteness of the Poisson distribution.

### VI.2 Coverage with systematic uncertainties

Introducing systematic uncertainties in the calculation of confidence intervals and tests of coverage leads to the question of what is meant by a repeated experiment. If we adhere to the traditional definition in which an experiment is repeated with fixed parameters such as efficiency or background rate, the algorithm presented here will inevitably yield over-coverage. Figure 10 shows the mean coverage ratio (mean here taken over six different signal expectation assumptions) as a function of different systematic uncertainties. The over-coverage will not only be increasing with increasing uncertainties but also be dependent on the signal expectation (figure 12).

#### VI.2.1 Bayesian Coverage.

The over-coverage described in the previous section is a consequence of the fact that efficiencies and background are not random variables (there is a true but unknown fixed efficiency and background rate) but they are treated as random variables in the construction of the confidence belt (equations (6) and (7)).
Thus, while in the construction we are using a PDF which is a convolution of a Poisson distribution with a Gaussian distribution, repeated measurements (with parameters fixed) will produce a Poisson distribution.

However, one has to keep in mind, that the distribution obtained in this way is not the underlying one. To infer from the measured Poisson distribution the underlying one, the signal efficiency and the background have to be taken into account. In particular, if these parameters are uncertain, there will not be a single underlying Poisson distribution, but a set of distributions that are weighted with the probabilities of the possible different efficiencies and backgrounds. In a way, we thus give different hypotheses different weights. To take this into account we modify the coverage test described in the previous section. Instead of drawing a measurement from Poisson distributions with predefined signal expectation and background, we draw the signal expectation and background prediction used in each simulated experiment from Gaussian distributions centered on the predefined true signal and background, where the width of the Gaussian is the associated systematic uncertainty. The measurement is then produced by taking these new values as means for the final Poisson distributions.
Since, by using this approach, we give different weights to different hypotheses, we call this modified coverage test Bayesian Coverage. In this way, the PDF used in the construction and in the coverage test are consistent with each other, and the algorithm should, per construction, give the correct coverage (except for discreteness effects). In particular, the coverage defined in this way should be independent of the magnitude of the uncertainties present in the experiment. Figure 12 shows the mean Bayesian coverage for different uncertainties in the signal efficiency together with the frequentist coverage. The mean is here taken over the 29 points in signal expectation space which were tested. As expected the Bayesian coverage ratio is nearly constant.
Thus, if we loosen the criteria on the definition of “repeated experiment”, allowing the “true” (unknown) efficiencies vary for each experiment repitition , the method has the desired statistical property of correct coverage.

#### VI.2.2 Remark on the choice of ensemble.

In the previous subsection we consider an ensemble in which the true value of the nuisance parameter, corresponding to $`ϵ^{}`$ in equation 7, is varied in each of the members of the ensemble. We find, as expected, that the POLE\- method fulfills the requirement of correct coverage with respect to this ensemble.
However, it can be argued, that this ensemble does not describe the situation usually encountered in experimental physics. The systematic uncertainty is a measurement uncertainty, i.e. not the true value of the nuisance parameter changes in each experiment but the measured one.
Studies for a few cases indicate, that with respect to such an ensemble the POLE method leads to moderate over-coverage Conrad:2002ur , Cousins:2002 .

## VII Limits on High Energy Cosmic Neutrinos.

In experimental situations where the systematic uncertainties are negligible small, limits or central confidence intervals can be calculated without evaluating the effects of the former. In more general situations, including systematic uncertainties in the calculation of the limits is essential, since it is a way of incorporating the real sensitivity of a given experiment to the quantity being measured.
In this section we will consider two real examples taken from published results of the AMANDA neutrino telescope, where the program POLE was used to include the systematics in the final results.

The AMANDA collaboration has published recently 90% confidence limits on the difuse flux of cosmic electron neutrinos and of neutrinos all flavours in the energy range between 5 TeV and 300 TeV Ahrens:2002wz . The analysis revealed zero events with an expected background from atmospheric neutrinos of 0.01 events. The systematic uncertainty in the signal efficiency for this analysis was $``$25%, determined from Monte Carlo studies. On top of that, the current theoretical systematic uncertainty in the atmospheric neutrino flux prediction in the energy range relevant to the analysis is about 30% Gaisser:2001a , which has to be taken into account as well. The effect of including the systematic uncertainties is shown in figure 12. The effect is to worsen the limits by about 10%.

Another example worth noticing are the results published by the same collaboration on searches for supersymmetric dark matter in the form of weakly interacting massive particles, WIMPs Ahrens:2002eb . In this analysis the uncertainties in signal efficiency range from 10% to 25% depending on the assumed signal spectrum and, additionally, an uncertainty in the background detection efficiency, estimated to be 20%, has to be taken into account in this case. A further complication to include these uncertainties in the calculation of the final limits arises since the efficiencies in signal and background detection are correlated. Moreover, the mentioned theoretical uncertainty on the overall normalization of the atmospheric neutrino flux has to be added. In figure 13 we show the limits to the muon flux from the center of the earth as a function of WIMP mass. The full lines show the limits for two diferent assumptions on the signal spectra, and the dashed lines the corresponding limits without including systematic uncertainties in the calculations.

For the purpose of illustration, we show in figure 14 the effect of including systematic uncertainties in the limit calculation for five different values (dashed lines). From bottom to top: a 20% uncertainty in background expectation only, and a 10%, 20% 30% and 40% uncertainty in signal efficiency (on top of the mentioned 20% background uncertainty). The absolute scale of the plot is arbitrary since we are just interested in showing the relative effect of the inclusion of systematic uncertainties in the limit calculation, with respect to the no-systematics case (full line). The figure shows the importance of correctly evaluating systematics and including them in the final result, since the effects can be important.

## VIII Conclusions

In this note we present a Monte Carlo algorithm for introducing systematic uncertainties in the evaluation of classical confidence intervals which allows to include uncertainties in the background prediction, in the background detection efficiency and in the signal detection efficiency, and correlations between them, by integrating over the (assumed) PDFs of these parameters. We apply the method for a Poisson process with background under the assumption of a Gaussian PDF describing the uncertainties. We present results where the construction has been performed using likelihood ratio ordering with and without conditioning.
Generally, the introduction of systematic uncertainties leads to an increase in confidence interval width. However, an interesting result is that likelihood ratio (as well as Neyman) confidence intervals which take into account the systematic uncertainty in the signal efficiency do not become larger with larger uncertainty, in the case that significantly less events have been observed than expected background. With respect to an ensemble with strictly identical experiments, introducing systematic uncertainties in the presented manner inevitably leads to over-coverage, increasing with the magnitude of systematic uncertainties. However, we show that with respect to an ensemble where the systematic uncertainties are taken into account in the coverage test by varying their true assumed values in each pseudoexperiment, the method presented here provides over-coverage only on the level already present due to the discreteness of the Poisson distribution.Both ensembles are ideal ensembles and have to be seen as approximations to the ensemble encountered in experimental physics. In summary, the algorithm presented here provides a practical and flexible way to quantitatively take into account systematic uncertainties present in experimental situations in the calculation of confidence intervals.

###### Acknowledgements.

We are thankful to Robert Cousins for valuable comments on the manuscript. We thank John Conway for making information available in the early stages of this work. Alexander Biron is acknowledged for careful reading of the manuscript in an earlier version.