text:
# Surface Superconductivity of Dirty Two-Band Superconductors: Applications to MgB₂.

## Abstract

The minimal magnetic field $`H_{c2}`$ destroying superconductivity in the bulk of a superconductor is smaller than the magnetic field $`H_{c3}`$ needed to destroy surface superconductivity if the surface of the superconductor coincides with one of the crystallographic planes and is parallel to the external magnetic field. While for a dirty single-band superconductor the ratio of $`H_{c3}`$ to $`H_{c2}`$ is a universal temperature-independent constant 1.6946, for dirty two-band superconductors this is not the case. I show that in the latter case the interaction of the two bands leads to a novel scenario with the ratio $`H_{c3}/H_{c2}`$ varying with temperature and taking values larger and smaller than 1.6946. The results are applied to $`\mathrm{MgB}_2`$ and are in agreement with recent experiments \[A. Rydh et al., cond-mat/0307445\].

Introduction. It is well-established that strong magnetic field destroys superconductivity. If an external field $`H`$ applied to a type-II superconductor exceeds the second critical field $`H_{c2}`$, the bulk order parameter in the superconductor vanishes. However, even for $`H>H_{c2}`$ superconductivity might still exist in a thin layer close to the surface if $`H`$ is smaller than the third critical field $`H_{c3}>H_{c2}`$Saintjames . In this paper I investigate the onset of superconductivity via surface nucleation for the field $`H`$ slightly below the threshold $`H_{c3}`$.

In their pioneering workSaintjames  Saint-James and de Gennes have shown that if the external magnetic field is applied parallel to the surface of an isotropic single-band superconductorfootnote\_1  with a temperature close to the transition temperature $`T_c`$, the ratio $`H_{c3}/H_{c2}`$ takes the universal value $`\eta =1.6946`$ independently of the superconducting material. It turns outSaintjames  that for $`H`$ slightly below $`H_{c3}`$ the superconducting order parameter exists within the distance $`\zeta (T)`$ (the coherence length of the superconductor) from the surface. For distances exceeding $`\zeta (T)`$ the order parameter approaches zero rapidly. The dependence of the ratio $`H_{c3}/H_{c2}`$ on the material propertiesSaintjames\_book ; Degennes\_book ; Strongin , sample geometry and topologyKogan ; Saintjames\_2 ; Moshchalkov , and temperatureAbrikosov ; Kulik ; Keller  has become a subject of intensive investigations.

A novel window for investigating surface superconductivity was opened after the discovery of the two-band superconductor $`\mathrm{MgB}_2`$Nagamatsu . Not only has it a relatively high ($`40`$ K) $`T_c`$ but also there exist two different superconducting gaps. As the consequence of this fact, various properties of $`\mathrm{MgB}_2`$ are quite different from those of single-band superconductors. For example, the anisotropy $`\gamma (T)=H_{c2}^{(ab)}(T)/H_{c2}^{(c)}(T)`$ (here, $`H_{c2}^{(ab)}(T)`$ and $`H_{c2}^{(c)}(T)`$ stand for the second critical fields in the $`ab`$ and $`c`$-directions respectively; note that the crystal of $`\mathrm{MgB}_2`$ is uniaxial) of the second critical field $`\mathrm{H}_{c2}`$ exhibits strong dependence on temperature, see e.g. Ref. Angst . For single-band superconductors this ratio is constant.

Another puzzle is that $`\gamma (T)`$ varies widely in different experimentsWelp . This can be attributed to the existence of surface superconductivity which might affect the observable values of $`H_{c2}`$ and, hence, the anisotropy. Consequently, the determination of the third critical field $`H_{c3}`$ is a very important problem. In a recent experimentRydh  it has been shown that $`H_{c3}/H_{c2}`$ for $`\mathrm{MgB}_2`$ might be reduced.

In the present paper I investigate the ratio $`H_{c3}/H_{c2}`$ for a dirty $`\mathrm{MgB}_2`$ crystal. The existence of two different gaps manifests itself through the remarkable dependence of $`H_{c3}/H_{c2}`$ on temperature. This is in sharp contrast with the case of a dirty single-gap superconductor where $`H_{c3}/H_{c2}=\eta `$ in the whole temperature range. For a magnetic field lying in the $`ab`$-plane of the $`\mathrm{MgB}_2`$ crystal I find that if one starts decreasing temperature, $`H_{c3}/H_{c2}`$ first exhibits a maximum at $`T0.99T_c`$ and then a minimum at $`T0.9T_c`$. As temperature decreases further, $`H_{c3}/H_{c2}`$ increases and tends to a value slightly below $`\eta `$, see Fig. 1. Naively, one could try to use the Ginzburg-Landau theory (GLT) in order to find $`H_{c3}/H_{c2}`$. However, as it will be explained below, this would lead to the ratio $`H_{c3}/H_{c2}=\eta `$ in the whole temperature range, i.e. one needs a more rigorous approach in order to explain the deviation of $`H_{c3}/H_{c2}`$ from $`\eta `$.

General formalism. An appropriate tool to investigate magnetic properties of dirty superconductors is the Usadel equationsUsadel . For two-band superconductors they have been derived by Koshelev and GolubovKoshelev  and by GurevichGurevich . Since I investigate the onset of superconductivity near $`H_{c2}`$ or $`H_{c3}`$, it is possible to write the Usadel equations in the linearized form

$$\omega f_\alpha +\left(\underset{j}{}\frac{D_{\alpha ,j}}{2}\left(_j\frac{2\pi i}{\mathrm{\Phi }_0}A_j\right)^2\right)f_\alpha =\mathrm{\Delta }_\alpha $$
(1)
$$\mathrm{\Delta }_\alpha =2\pi T\underset{\omega >0}{\overset{\omega _D}{}}\mathrm{\Lambda }_{\alpha \beta }f_\beta $$
(2)
Here, $`\omega =2\pi T(n+1/2),n=0,1,\mathrm{}`$ and $`\omega _D`$ are the Matsubara and cutoff phonon frequencies. $`D_{\alpha ,j}`$ is the diffusion coefficient of the band $`\alpha =1,2`$ along the direction $`j=x,y,z`$. The indices 1 and 2 correspond to the $`\sigma `$\- and $`\pi `$-bands respectively. $`𝐀`$ is the vector potential. $`\mathrm{\Delta }_\alpha `$ and $`f_\alpha `$ are the superconducting gap and anomalous green function for the band $`\alpha `$. The matrix $`\widehat{\mathrm{\Lambda }}`$ represents the strength of the coupling parameters and has the values $`\lambda _{11}0.81`$, $`\lambda _{22}=0.285`$, $`\lambda _{12}0.119`$, and $`\lambda _{21}0.09`$, see Ref. Mazin . In the present paper I will concentrate on two geometries: i) the magnetic field is parallel to the $`c`$-axis and the surface of the crystal; ii) the surface of the superconductor coincides with the $`ab`$-plane and the field H lies in the $`ab`$-plane. As I will show, in case i $`H_{c3}/H_{c2}=\eta `$ at any temperature. In case ii, $`H_{c3}/H_{c2}`$ is shown in Fig. 1. I assume that $`𝐇`$ is parallel to the $`z`$-axis.

Choosing the gauge as $`A_y=Hx`$ I look for the solution of the form $`f_\alpha f_\alpha (\omega ,x)\mathrm{exp}\left(ik_yy+ik_zz\right)`$ and $`\mathrm{\Delta }_\alpha \mathrm{\Delta }_\alpha (x)\mathrm{exp}\left(ik_yy+ik_zz\right)`$. In general, Eqs. (1) and (2) define a sequence of solutions corresponding to different eigenvalues $`H=H_{c2}`$ or $`H=H_{c3}`$. One should look for the maximal possible values of $`H_{c2}`$ or $`H_{c3}`$. This corresponds to the case $`k_z=0`$. Substituting the Ansatz for $`f_\alpha `$ and $`\mathrm{\Delta }_\alpha `$ into (1) and (2) I obtain the system of equations

$$\widehat{\mathrm{\Lambda }}^1\left(\begin{array}{c}\mathrm{\Delta }_1(x)\\ \mathrm{\Delta }_2(x)\end{array}\right)=\left(\begin{array}{cc}2\pi T\underset{\omega >0}{\overset{\omega _D}{}}\frac{1}{\omega +\widehat{H}_1(x_0)}& 0\\ 0& 2\pi T\underset{\omega >0}{\overset{\omega _D}{}}\frac{1}{\omega +\widehat{H}_2(x_0)}\end{array}\right)\left(\begin{array}{c}\mathrm{\Delta }_1(x)\\ \mathrm{\Delta }_2(x)\end{array}\right),$$
(3)
where

$$\widehat{H}_\alpha (x_0)=\frac{\stackrel{~}{D}_{\alpha ,x}}{2}\frac{^2}{x^2}+\frac{\stackrel{~}{D}_{\alpha ,z}}{2}\left(\frac{2\pi H}{\mathrm{\Phi }_0}\right)^2\left(xx_0\right)^2,$$
(4)
with $`\stackrel{~}{D}_{\alpha ,x}=\stackrel{~}{D}_{\alpha ,z}=D_{\alpha ,a}`$ for case i and $`\stackrel{~}{D}_{\alpha ,x}=D_{\alpha ,a}`$ and $`\stackrel{~}{D}_{\alpha ,z}=D_{\alpha ,c}`$ for case ii ($`D_{\alpha ,a}=D_{\alpha ,b}D_{\alpha ,c}`$ are the diffusion coefficients along the crystallographic axes). $`x_0=k_y\mathrm{\Phi }_0/2\pi H`$ is the parameter characterizing how far away the superconducting nucleus is situated from the surface. Note that $`x_0`$ is the same for the both bands. The operator (4) can be rewritten in the form $`\left(\pi H/\mathrm{\Phi }_0\right)\sqrt{\stackrel{~}{D}_{\alpha ,x}\stackrel{~}{D}_{\alpha ,z}}\widehat{h}_\alpha (x_0^{})`$, with

$$\widehat{h}_\alpha (x_{\alpha ,0}^{})=\frac{^2}{x_{\alpha }^{}{}_{}{}^{2}}+\left(x_\alpha ^{}x_{\alpha ,0}^{}\right)^2,$$
(5)
where I have made the variable substitution $`x=\beta _\alpha x_\alpha ^{}`$ and $`x_0=\beta _\alpha x_{\alpha ,0}^{}`$, with $`\beta _\alpha =\left(\stackrel{~}{D}_{\alpha ,x}/\stackrel{~}{D}_{\alpha ,z}\right)^{1/4}\left(\mathrm{\Phi }_0/2\pi H\right)^{1/2}`$.

The system (3) should be solved with the boundary conditions (BC) $`\mathrm{\Delta }_\alpha /x|_{x=0}=0`$, and $`\mathrm{\Delta }_\alpha (x+\mathrm{})0`$, $`\alpha =1,2`$ valid for geometries i and ii, see above. For $`H_{c3}`$ the BC are well-established for dirty superconductorsDegennes\_book . For $`H_{c2}`$ the application of these BC gives the same result as the BC requiring the appearance of a superconducting nucleus in the bulk. The procedure of finding $`H_{c3}/H_{c2}`$ is following: first, set $`x_0=0`$ in (3) and find the maximal possible field $`H`$ for which the solution satisfying the BC exists. This gives $`H_{c2}`$. Next, for $`x_00`$ find the maximal field $`H=H(x_0)`$ for which the solution of (3) exists. Then, $`H_{c3}=\mathrm{max}_{x_0}\left\{H(x_0)\right\}`$. I would like to mention that there are complementary approaches for calculating $`H_{c2}`$ based on macroscopic theoryMiranovic ; Dahm  and GLTAskerzade ; Zhitomirsky .

Here, it is instructive to study briefly the case of a single-gap superconductor. This corresponds to $`\lambda _{12}=\lambda _{21}=0`$. $`H_{c2}`$ and $`H_{c3}`$ are then determined by those for band 1 (as $`\lambda _{11}>\lambda _{22}`$). The solution for $`\mathrm{\Delta }_1(x)`$ is proportional to the ground state wave function of the operator $`\widehat{H}_1(x_0)`$ and $`\mathrm{\Delta }_2(x)=0`$. Substituting this Ansatz into (3) I obtain the transcendental equation of the form

$$1\lambda _{11}\underset{\omega >0}{\overset{\omega _D}{}}\frac{1}{\omega +\left(\pi H/\mathrm{\Phi }_0\right)\sqrt{\stackrel{~}{D}_{1,x}\stackrel{~}{D}_{1,z}}ϵ_0(x_{1,0}^{})}=0,$$
(6)
with $`ϵ_0(x_{\alpha ,0}^{})`$ the lowest eigenvalue of the operator $`\widehat{h}_\alpha (x_{\alpha ,0}^{})`$. The field $`H_{c2}`$ can be found as the solution of the above equation for $`x_{1,0}^{}=0`$; note that $`ϵ_0(0)=1`$. Assume, a certain value of the magnetic field $`H_{c2}`$ is found; let us change the parameter $`x_{1,0}^{}`$. This leads to the decrease of the eigenvalue $`ϵ_0(x_{1,0}^{})`$Saintjames . In order to satisfy Eq. (6) one has to increase the field $`H`$; that is why $`H_{c3}>H_{c2}`$. The minimal $`ϵ_0(x_{1,0}^{})`$ can be realized for $`x_{1,0}^{}=0.7618`$Saintjames  and is equal to $`0.5901`$Saintjames . This means that $`H_{c3}/H_{c2}=1/0.5901=\eta `$ for any temperature $`T`$. Remarkably, it is not necessary to solve Eq. (6) in order to find the ratio $`H_{c3}/H_{c2}`$, although the determination of $`H_{c2}`$ or $`H_{c3}`$ alone would require the complete analysis.

Case i. In this case, the ratio $`H_{c3}/H_{c2}=\eta `$ at all temperatures. This is a consequence of the fact that the operators $`\widehat{H}_1(x_0)`$ and $`\widehat{H}_2(x_0)`$ have identical eigenfunctions (since $`\stackrel{~}{D}_{1,x}/\stackrel{~}{D}_{1,z}=\stackrel{~}{D}_{2,x}/\stackrel{~}{D}_{2,z}=1`$). The functions $`\mathrm{\Delta }_1(x)`$ and $`\mathrm{\Delta }_2(x)`$ are proportional to the ground state wave function of the operator $`\widehat{H}_1(x_0)`$ (or $`\widehat{H}_2(x_0)`$). The equation determining the critical fields $`H_{c2}`$ and $`H_{c3}`$ has the form $`F(\pi D_{1,a}Hϵ_0(x_{1,0}^{})/\mathrm{\Phi }_0,\pi D_{2,a}Hϵ_0(x_{2,0}^{})/\mathrm{\Phi }_0)=0`$, with $`F(y_1,y_2)`$ a certain function of two arguments. Note that in the present case $`\beta _1=\beta _2`$ and, consequently, $`x_{1,0}^{}=x_{2,0}^{}`$ and $`ϵ_0(x_{1,0}^{})=ϵ_0(x_{2,0}^{})`$. The maximal value of $`H_{c3}`$ can be realized for $`x_{1,0}^{}=0.7618`$ and is equal to $`\eta H_{c2}`$ for all temperatures.

Case ii. If the magnetic field lies in the $`ab`$-plane, the operators $`\widehat{H}_1(x_0)`$ and $`\widehat{H}_2(x_0)`$ have different eigenfunctions. This leads to a complicated transcendental equation depending on all eigenvalues of the operators $`\widehat{H}_1(x_0)`$ and $`\widehat{H}_2(x_0)`$ and not only on the ground state ones. Eqs. (3) can be solved via expanding functions $`\mathrm{\Delta }_1(x)`$ and $`\mathrm{\Delta }_2(x)`$ over the eigenfunctions of the operators $`\widehat{H}_1(x_0)`$ and $`\widehat{H}_2(x_0)`$. I have truncated the basis of the operators $`\widehat{H}_1(x_0)`$ and $`\widehat{H}_2(x_0)`$ to subspaces consisting of 70 eigenfunctions and solved the system (3) numerically. In Fig. 1 I show the results of the numerics. Here, I take $`D_{1,a}/D_{1,c}=40.0`$, $`D_{2,a}/D_{2,c}=0.665`$. These ratios can be obtained using the results for the average velocity on the $`\mathrm{MgB}_2`$ Fermi surfaces and assuming isotropic scattering, see Refs. Golubov ; Brinkman . The ratio $`D_{2,c}/D_{1,c}`$ takes three values: 100, 300, and 600. This choice is motivated by the facts that the ratio $`D_{2,c}/D_{1,c}100`$ can be obtained assuming that the scattering rate of electrons is the same in both bands. On the other hand $`R=D_{2,c}/D_{1,c}=600`$ gives a better fit with experiments on the anisotropy measurementsGolubov .

The results are shown in Fig. 1. For temperatures $`T0.6T_c`$ I have found that the ratio $`H_{c3}/H_{c2}`$ is nearly constant and has a value slightly below $`\eta `$. This can be explained by the fact that at low temperatures the fields $`H_{c2}`$ and $`H_{c3}`$ are determined mostly by the $`\sigma `$-band whose coherence length is much smaller than that of the $`\pi `$-band. At low $`T`$ the magnetic field $`H_{c2}`$ depends on the ground state eigenvalues of the operators $`\widehat{H}_1(x_0)`$ and $`\widehat{H}_2(x_0)`$ and the contribution of excited states is negligibleGurevich ; Koshelev . The ratios $`\left(\stackrel{~}{D}_{1,x}/\stackrel{~}{D}_{1,z}\right)^{1/4}`$ and $`\left(\stackrel{~}{D}_{2,x}/\stackrel{~}{D}_{2,z}\right)^{1/4}`$ determining the length $`x_0`$ are equal to $`2.5`$ and $`0.9`$ respectively. This means that one can maximize the field $`H_{c3}`$ by choosing $`x_{1,0}^{}=0.7618`$. The length $`x_{2,0}^{}`$ then is large and the ground-state eigenvalue of the operator $`\widehat{h}_2(x_{2,0}^{})`$ is close to $`1`$. The ratio $`H_{c3}/H_{c2}`$ then can be calculated as follows: take the zero temperature expression for $`H_{c2}`$Gurevich ; Golubov  and make there a substitution $`D_{1,j}D_{1,j}/\eta `$. At low temperaturesGurevich ; Golubov , $`H_{c2}=\mathrm{\Phi }_0T_c\mathrm{exp}\left(g/2\right)/2\gamma \left(D_{1,a}D_{1,c}D_{2,a}D_{2,c}\right)^{1/4},`$ with $`g=\left(\lambda _0^2/w^2+\mathrm{ln}^2\kappa /4+2\lambda _{}\mathrm{ln}\kappa /w\right)^{1/2}\lambda _0/w,`$ $`\kappa =D_{2,a}D_{2,c}/D_{1,a}D_{1,c}`$, $`\lambda _{}=\lambda _{11}\lambda _{22}`$, $`\lambda _0=\left(\lambda _{}^2+4\lambda _{12}\lambda _{21}\right)^{1/2}`$, $`w=\lambda _{11}\lambda _{22}\lambda _{12}\lambda _{21}`$, and $`\mathrm{ln}\gamma =0.5772`$. This procedure yields $`H_{c3}/H_{c2}=1.688,1.690,`$ and $`1.691`$ for $`D_{2,c}/D_{1,c}=100,300,`$ and $`600`$ respectively. The values obtained in the numerics are slightly larger (but still smaller than $`\eta `$) due to a small contribution to the adjustment of $`x_0`$ from the $`\pi `$-band.

If one increases temperature, the ratio $`H_{c3}/H_{c2}`$ decreases and exhibits a minimum at $`T0.9T_c`$. Then, the value of $`H_{c3}/H_{c2}`$ goes up and takes a maximum at $`T0.99T_c`$. At $`T=T_c`$, $`H_{c3}/H_{c2}=\eta `$. The nontrivial behavior of the ratio $`H_{c3}/H_{c2}`$ in $`\mathrm{MgB}_2`$ is due to the changing relative importance of the $`\pi `$-band. While at low $`T`$ it is unimportant, at high $`T`$ it gives a comparable with the $`\sigma `$-band contribution to the fields $`H_{c2}`$ and $`H_{c3}`$.

Let us analyze the situation close to $`T_c`$ in more detail. In particular, let us explain why at $`T=T_c`$, $`H_{c3}/H_{c2}=\eta .`$ For $`T_cTT_c`$ the field $`H_{c3}`$ is small and so are the eigenvalues of the operators $`\widehat{H}_1(x_0)`$ and $`\widehat{H}_2(x_0)`$, i.e. one can use the expansion $`_{\omega >0}^{\omega _D}1/\left(\omega +\widehat{H}_\alpha (x_0)\right)_{\omega >0}^{\omega _D}1/\omega _{\omega >0}^{\omega _D}\widehat{H}_\alpha (x_0)/\omega ^2+\mathrm{}`$ and Eqs. (3) can be rewritten in the form

$$\widehat{W}\left(\begin{array}{c}\mathrm{\Delta }_1(x)\\ \mathrm{\Delta }_2(x)\end{array}\right)=\left(\begin{array}{cc}\widehat{R}_1(x_0)& 0\\ 0& \widehat{R}_2(x_0)\end{array}\right)\left(\begin{array}{c}\mathrm{\Delta }_1(x)\\ \mathrm{\Delta }_2(x)\end{array}\right),$$
(7)
with $`\stackrel{~}{W}=\widehat{\mathrm{\Lambda }}^1\mathrm{ln}\left(2\gamma \omega _D/\pi T\right)𝟙_\mathrm{𝟚}`$ and $`\widehat{R}_\alpha (x_0)=2\pi T_{\omega >0}^{\omega _D}1/\left(\omega +\widehat{H}_\alpha (x_0)\right)2\pi T_{\omega >0}^{\omega _D}1/\omega `$. $`T_c`$ is determined by $`\mathrm{det}\widehat{W}=0`$. Solving the system (7) for $`\mathrm{\Delta }_1(x)`$ I obtain

$$\left(W_{11}\widehat{R}_2(x_0)+W_{22}\widehat{R}_1(x_0)\widehat{R}_2(x_0)\widehat{R}_1(x_0)\right)\mathrm{\Delta }_1(x)=0.$$
(8)
Eq. (8) allows to determine the fields $`H_{c2}`$ and $`H_{c3}`$ in a regular way. To lowest order, one can neglect the term $`\widehat{R}_2(x_0)\widehat{R}_1(x_0)`$.

The equation $`\left(W_{11}\widehat{R}_2(x_0)+W_{22}\widehat{R}_1(x_0)\right)\mathrm{\Delta }_1(x)=0`$ has the ground-state solution of the same form as Eq. (3) for a single-gap superconductor (the case $`\lambda _{12}=\lambda _{21}=0`$) with $`\stackrel{~}{D}_{1,x}D_X=W_{22}\stackrel{~}{D}_{1,x}+W_{11}\stackrel{~}{D}_{2,x}`$ and $`\stackrel{~}{D}_{1,z}D_Z=W_{22}\stackrel{~}{D}_{1,z}+W_{11}\stackrel{~}{D}_{2,z}`$ and the problem of finding $`H_{c3}`$ becomes equivalent to the original one considered by Saint-James and de GennesSaintjames . Consequently, to lowest order in $`T_c`$ the ratio $`H_{c3}/H_{c2}`$ has the same value $`\eta `$ as in the case of a single-gap superconductor. The approximation described above is equivalent to the GLT. Consequently, the GLT is unable to explain deviations of $`H_{c3}/H_{c2}`$ from $`\eta `$.

The maximum of $`H_{c3}/H_{c2}`$ takes place very close to $`T_c`$ and is at the boundary of the accuracy of the present numerical calculations. Hence, an analytical approach would be useful. Temperature corrections to $`H_{c3}/H_{c2}`$ can be found by expanding operators $`\widehat{R}_1(x_0)`$ and $`\widehat{R}_2(x_0)`$ to second order in $`(T_cT)`$. One can decompose the operator in the left-hand side of (8) as a sum $`\widehat{L}_1(H,x_0)+\widehat{L}_2(H,x_0)`$, with $`\widehat{L}_1(H,x_0)T_cT`$ and $`\widehat{L}_2(H,x_0)\left(T_cT\right)^2`$. Let $`|\varphi _0`$ be the solution of the equation $`\widehat{L}_1(H,x_0)|\varphi _0=0`$ and $`H=H^{(0)}(x_0)`$ the critical field to this order. The correction to the eigenvalue can be found perturbatively and are determined by the implicit relation

$`\varphi _0|\widehat{L}_1(H,x_0)|\varphi _0+\varphi _0|\widehat{L}_2(H^{(0)},x_0)|\varphi _0=0.`$ (9)
To this order $`x_0=0.7618\left(D_X/D_Z\right)^{1/4}\left(1/2\pi H^{(0)}\right)^{1/2}`$. Temperature corrections to $`H_{c3}`$ due to change in $`x_0`$ are proportional to $`\left(T_cT\right)^3`$ and can be neglected. Since $`|\varphi _0`$, $`\widehat{L}_1(H,x_0)`$, and $`\widehat{L}_2(H,x_0)`$ are known, one can find $`H_{c3}/H_{c2}`$ analytically. Straightforward but quite cumbersome calculationsGorokhov  show that near $`T_c`$, $`H_{c3}/H_{c2}\eta +b(T_cT)/T_c`$, with $`b1`$ for $`D_{2,c}/D_{1,c}`$ in the range from 100 to 600, in accordance with the numerics, see inset in Fig. 1.

Experiment. Recent experimentsRydh  show that that the ratio $`H_{c3}/H_{c2}`$ is reduced in the case ii. The values $`H_{c3}/H_{c2}1.5`$ in the temperature range 20–30 K have been reported. For the case i $`H_{c3}/H_{c2}1.7`$Welp ; Rydh . The present theory gives that $`H_{c3}/H_{c2}=\eta `$ for case i and $`H_{c3}/H_{c2}<\eta `$ for case ii, in agreement with Welp ; Rydh . Theoretical calculations showGolubov  that the anisotropy of $`\mathrm{MgB}_2`$ is distributed over a wider temperature range in an experiment that theory suggests. A somewhat similar situation takes place in the present work for the ratio $`H_{c3}/H_{c2}`$. There are two main sources of deviation between theory and experiment. First, surface qualityHart  might affect $`H_{c3}/H_{c2}`$. Second, $`\mathrm{MgB}_2`$ is situated somewhere at the boundary of the applicability of the weak-coupling BCS-theory. It would be very interesting to repeat the calculation done in the present paper starting from the Eliashberg equations.

The method described above can be generalized for an arbitrary direction of crystallographic axes with respect to the surface of a superconductor. For strongly anisotropic superconductors, surface superconductivity might disappear if the surface does not coincide with crystallographic planesKogan . EstimatesGorokhov  show that $`\mathrm{MgB}_2`$ is sufficiently anisotropic in order to observe this kind of effects. The detailed analysis of the onset of surface superconductivity in this case is challenging for both theorists and experimentalists.

In conclusion, I have presented the calculation of the ratio $`H_{c3}/H_{c2}`$ for the two-band superconductor $`\mathrm{MgB}_2`$ in the dirty limit. Remarkably, in contrast to the case of a single-gap superconductor, the above ratio is temperature-dependent. The Ginzburg-Landau theory is unable to explain deviations of $`H_{c3}/H_{c2}`$ from $`\eta `$.

The present work is supported by the Packard Foundation. D. A. G. thanks M. Angst and A. E. Koshelev for helpful discussions.


format: text

text:
# Evidence for entangled states of two coupled flux qubits

## Abstract

We have studied the low-frequency magnetic susceptibility of two inductively coupled flux qubits using the impedance measurement technique (IMT), through their influence on the resonant properties of a weakly coupled high-quality tank circuit. In a single qubit, an IMT dip in the tank’s current–voltage phase angle at the level anticrossing yields the amplitude of coherent flux tunneling. For two qubits, the difference (IMT deficit) between the sum of single-qubit dips and the dip amplitude when both qubits are at degeneracy shows that the system is in a mixture of entangled states (a necessary condition for entanglement). The dependence on temperature and relative bias between the qubits allows one to determine all the parameters of the effective Hamiltonian and equilibrium density matrix, and confirms the formation of entangled eigenstates.

While its exact role in the power of quantum computers can still be debated, there is no doubt that entanglement is a key feature of quantum registers, absent in their classical counterparts N&C . From the point of view of scalability, implementing quantum registers as integrated circuits is attractive. However, demonstrating entanglement in macroscopic solid-state systems remains a daunting task, even after the observation of quantum coherence in many types of superconducting qubits Nakamura99 ; Lukens ; VdW ; Martinis ; Vion  and also entanglement in charge Pashkin  and current-biased Josephson-junction Berkley  qubits, using pulse and spectroscopic techniques respectively. In this paper we investigate two coupled three-junction flux qubits Mooij ; Majer  using a different method, providing a simple criterion for eigenstate entanglement, by studying the latter’s influence on the two-qubit magnetic susceptibility.

The impedance measurement technique (IMT) Greenberg:2002  relies on monitoring the current and voltage in a high-quality, low-frequency resonant tank circuit inductively coupled to the qubit. In various modifications, this method was applied to map the Josephson potential profile APL\_IMT , and to observe Landau–Zener transitions Landau-Zener  and Rabi oscillations Rabi  in a three-junction flux qubit. Due to the tank’s weak coupling to the qubit and high quality factor, decoherence times can be as high as 2.5 $`\mu `$Rabi . Note that the IMT approach requires a small decoherence rate compared to only the tunneling amplitudes (0.5$``$1 GHz), not the tank frequency (10$``$20 MHz). This paper generalizes to two qubits an IMT study of coherent tunneling IMT:2003 . The tank is fed a small-amplitude ac signal; its effective impedance, and therefore the current–voltage phase angle, are sensitive to the qubit susceptibility $`\chi `$. For, say, one qubit at temperature $`T=0`$, $`\chi _\mathrm{\Phi }^2E_{}`$, the curvature of the ground-state energy $`E_{}`$ vs external flux $`\mathrm{\Phi }`$ \[a different expression for $`\chi `$ is used in Eq. (4) below\]. Thus, the level anticrossing due to flux tunneling is revealed by a large peak in $`|\chi |`$.

Our system of two flux qubits inductively coupled to each other and to the tank is shown in Fig. 1. The qubits were fabricated out of aluminum by conventional shadow evaporation, nominally 1 $`\mu `$m apart, at the centre of a niobium pickup coil. The area of each qubit was $`80`$ $`\mu `$m2, with self-inductance $`L_{a/b}39`$ pH; the two biggest junctions in each qubit had critical current $`400`$ nA and Coulomb energy $`e^2/2Ch3.2`$ GHz (for the third junction, these values are 10$``$20% smaller and larger, respectively). The mutual inductance between the qubits $`M_{ab}=2.7`$ pH was estimated numerically from the electron micrograph. The Nb coil, which together with an external capacitor of $`C_\mathrm{T}470`$ pF forms a parallel tank circuit with inductance $`L_\mathrm{T}130`$ nH, resonant frequency $`\omega _\mathrm{T}/2\pi =20.139`$ MHz and quality factor $`Q_\mathrm{T}\omega _\mathrm{T}R_\mathrm{T}C_\mathrm{T}=1680`$ (at 10 mK; $`R_\mathrm{T}`$ is an effective tank resistance) rectify , was fabricated using e-beam lithography. The external magnetic flux through the qubits was created by the dc component of the current in the coil $`I_{\mathrm{dc1}}`$, and by the bias current $`I_{\mathrm{dc2}}`$ through a wire close to one of the qubits. This allowed independent control of the bias in each qubit.

The system of Fig. 1 is described by the Hamiltonian $`H=H_0+H_\mathrm{T}+H_{\mathrm{int}}+H_{\mathrm{diss}}`$, where the two-qubit Hamiltonian in the two-state approximation is expressed through Pauli matrices as Alec03

$$H_0=\mathrm{\Delta }_a\sigma _x^{(a)}\mathrm{\Delta }_b\sigma _x^{(b)}+ϵ_a\sigma _z^{(a)}+ϵ_b\sigma _z^{(b)}+J\sigma _z^{(a)}\sigma _z^{(b)},$$
(1)
$`H_\mathrm{T}`$ is the tank Hamiltonian (a harmonic oscillator), the qubit–tank interaction is

$$H_{\mathrm{int}}=(\lambda _a\sigma _z^{(a)}+\lambda _b\sigma _z^{(b)})I_\mathrm{T},$$
(2)
and $`H_{\mathrm{diss}}`$ describes the standard weak coupling of the qubits to a dissipative bath Weiss . Here $`I_\mathrm{T}`$ is the current through $`L_\mathrm{T}`$. The coefficients are $`\lambda _{a/b}=M_{a/b,\mathrm{T}}I_{a/b}`$, where $`M_{a/b,\mathrm{T}}`$ is the qubit–tank mutual inductance and $`I_{a/b}`$ is the magnitude of the persistent current in the corresponding qubit. The qubit biases are given by $`ϵ_a=I_a\mathrm{\Phi }_0(f_\mathrm{x}\frac{1}{2}+f_{\mathrm{shift}})`$, $`ϵ_b=I_b\mathrm{\Phi }_0(f_\mathrm{x}\frac{1}{2}+\eta f_{\mathrm{shift}})`$, where $`f_\mathrm{x}=\mathrm{\Phi }/\mathrm{\Phi }_0`$ accounts for the external flux $`\mathrm{\Phi }I_{\mathrm{dc1}}`$ created by the Nb coil in both qubits, while the parameters $`f_{\mathrm{shift}}I_{\mathrm{dc2}}`$ and $`\eta =M_{b\mathrm{w}}/M_{a\mathrm{w}}<1`$ describe the bias difference between the qubits created by the additional wire. Here $`M_{a\mathrm{w}}`$ ($`M_{b\mathrm{w}}`$) are the mutual inductances between the $`a`$ ($`b`$) qubit and the additional dc wire (for our sample, $`M_{a\mathrm{w}}`$ and $`M_{b\mathrm{w}}`$ were calculated numerically, yielding $`\eta =0.32`$). The qubit–qubit coupling $`J=M_{ab}I_aI_b`$ is positive because the two qubits are in the same plane side by side, leading to antiferromagnetic coupling (according to the north-to-south attraction law).

The effective inductance of the tank, and therefore its resonant frequency, depend on the state of the qubits. In the IMT method Greenberg:2002 , this dependence is observed by continuously monitoring the phase angle $`\mathrm{\Theta }`$ between the average tank voltage and bias current $`I_{\mathrm{bias}}(t)`$ fed into the tank (Fig. 1). For small enough $`I_{\mathrm{bias}}`$, $`\mathrm{tan}\mathrm{\Theta }`$ is determined by the real part of the two-qubit contribution to the tank susceptibility $`\chi ^{}`$ at $`\omega _\mathrm{T}`$ SmirnovR2003 , viz.,

$$\mathrm{tan}\mathrm{\Theta }=\frac{Q_\mathrm{T}}{L_\mathrm{T}}\chi ^{}(\omega _\mathrm{T}),$$
(3)
where $`\chi (\omega )`$ has a Kubo-type linear-response expression through retarded Green’s functions of the qubit operators $`\sigma _z^{(a/b)}`$. For weak damping $`H_{\mathrm{diss}}`$, the latter can be calculated with the equilibrium density matrix corresponding to $`H_0`$ in Eq. (1). It can be generally assumed that the eigenvalues $`E_\mu `$ of $`H_0`$, $`\mu =1,2,3,4`$, are non-degenerate and the eigenstates normalized, $`\nu |\mu =\delta _{\mu \nu }`$. Taking into account the qubits’ interaction with a dissipative environment SmirnovR2003 ; Blum , a standard calculation gives SmirnovUnpub

$$\chi (\omega )=\underset{\mu \nu }{}\frac{\rho _\mu \rho _\nu }{\mathrm{}\omega +E_\mu E_\nu +i\mathrm{}\mathrm{\Gamma }_{\mu \nu }}R_{\mu \nu },$$
(4)
where $`\rho _\mu =e^{E_\mu /T}/[_\nu e^{E_\nu /T}]`$ is the population of energy level $`\mu `$, $`\mathrm{\Gamma }_{\mu \nu }`$ in the energy denominators are decoherence rates, and the real matrix elements are

$`R_{\mu \nu }`$ $`=`$ $`\lambda _a^2\mu |\sigma _z^{(a)}|\nu \nu |\sigma _z^{(a)}|\mu +\lambda _b^2\mu |\sigma _z^{(b)}|\nu \nu |\sigma _z^{(b)}|\mu `$ (5)
$`+\lambda _a\lambda _b\mu |\sigma _z^{(a)}|\nu \nu |\sigma _z^{(b)}|\mu `$
$`+\lambda _a\lambda _b\mu |\sigma _z^{(b)}|\nu \nu |\sigma _z^{(a)}|\mu .`$
Substitution into Eq. (3) yields

$$\mathrm{tan}\mathrm{\Theta }=2\frac{Q_\mathrm{T}}{L_\mathrm{T}}\underset{\mu <\nu }{}\frac{\rho _\mu \rho _\nu }{E_\nu E_\mu }R_{\mu \nu }.$$
(6)
At low frequencies $`\omega _\mathrm{T}|E_\mu E_\nu |/\mathrm{}`$ and for weak damping $`\mathrm{\Gamma }_{\mu \nu }|E_\mu E_\nu |/\mathrm{}`$, the $`\mathrm{\Gamma }_{\mu \nu }`$ do not affect $`\mathrm{tan}\mathrm{\Theta }`$, but are responsible for establishing the equilibrium distribution.

The first two terms in Eq. (5) are clearly positive, and non-zero even if the two-qubit states are factorized. The first (second) term corresponds to the contribution of qubit $`a`$ ($`b`$), and peaks near that qubit’s degeneracy point. These contributions are practically independent of whether the qubits’ degeneracy points coincide or not.

The last two terms in Eq. (5) describe coherent flipping of both qubits, which is *only* possible for non-factorizable (entangled) eigenstates $`|\mu `$, $`|\nu `$. For $`J>0`$, these terms are found to be negative. Therefore the difference (*IMT deficit*) between the coinciding two-qubit IMT dip and the sum of two single-qubit dips provides a measure of how coherent the two-qubit dynamics is, that is, whether entangled eigenstates of $`H_0`$ are formed SmirnovUnpub . Eigenstate entanglement is a necessary condition for the equilibrium state to be entangled; of course, the eigenstate populations play a role as well. However, it is already a sufficient condition for performing quantum gate operations (if decoherence times are long enough), since one can also initialize the system in a nonequilibrium state, say, by letting it relax at a point in parameter space for which there is a large gap above the ground state, and subsequently returning to the working point Nakamura99 ; Lukens .

The measurements on the system of Fig. 1 allowed to determine all parameters of $`H_0`$. The data are presented in Fig. 2 ($`\mathrm{tan}\mathrm{\Theta }`$ vs flux bias) and Fig. 4 ($`T`$-dependence of dip amplitudes).

First, consider $`I_{\mathrm{dc2}}=27.3`$ ($`32.7`$$`\mu `$A, i.e., unequally biased qubits. In this case, whenever one qubit is near-degenerate, the other is strongly biased, putting it into a unique classical ground state (at the low end of our $`T`$-range) in which it behaves trivially. Hence, one essentially observes single-qubit properties. We match the corresponding IMT dips of Fig. 2 with the predictions of Eqs. (6), (5) for both the width and amplitude; see Fig. 3. This yields $`\mathrm{\Delta }_a/h=550`$ MHz, $`\mathrm{\Delta }_b/h=450`$ MHz, and the persistent currents $`I_aI_bI_\mathrm{p}=320`$ nA, for a flux-to-energy bias conversion factor $`I_\mathrm{p}\mathrm{\Phi }_0/h=990`$ GHz Delta . The width of the dips is practically $`T`$-independent, as expected IMT:2003 . The $`T`$-dependence of the amplitudes (Fig. 4, squares and triangles) agrees with these values for $`\mathrm{\Delta }_{a/b}`$. (The saturation below 35–40 mK IMT:2003  is likely due to a discrepancy between the mixing-chamber and sample temperatures.)

We can determine the coupling constant in Eq. (1) from the *single*-qubit measurements, as $`J/h=M_{ab}I_\mathrm{p}^2/h=410`$ MHz. This agrees well with the value $`J/h=420`$ MHz obtained by fitting the $`T`$-dependence of the coincident IMT dip amplitude above the saturation temperature (Fig. 4, circles) with the theoretical curve.

Comparison of the single-qubit and coincident IMT dips shows that the contribution to $`\mathrm{tan}\mathrm{\Theta }`$ due to eigenstate entanglement is significant. Indeed, the amplitude of the central dip in Fig. 2 at $`T=50`$ mK is 1.12, compared to the sum 1.69 of the single-qubit dips. This means that the entangled terms \[last two in Eq. (5)\] are responsible for a contribution $`0.57`$. As the dashed line in Fig. 3 shows, dropping these terms would yield the coincident IMT dip in marked deviation from both the full theory and experiment.

The IMT deficit also confirms the sign of the qubit–qubit interaction: in the case of ferromagnetic coupling, the coincident dip should be *larger* than the sum of the single-qubit ones.

Plotting the ratio of the coincident dip to the sum of single-qubit ones (Fig. 4, inset), we see that it grows with $`T`$, approaching $`1`$. This is to be expected: thermal excitations tend to destroy coherent correlations between the qubits (which will then behave as independent quantum systems with zero IMT deficit)—eventually, of course, destroying the single-qubit IMT dips as well.

At 50 mK, $`T`$ is comparable to the characteristic qubit energies (at the two-qubit degeneracy point, the top excited state is $`100`$ mK above the ground state). The measurement time equals at least the tank saturation time $`Q_\mathrm{T}/\omega _\mathrm{T}`$ and is longer in practice, 1$``$10 ms. Since this exceeds any conceivable qubit relaxation time $`\mathrm{\Gamma }_{\mu \nu }^1`$, the system has time to equilibrate. Indeed, the good quantitative agreement between experiment (Figs. 2, 4) and theory \[Eq. (6)\] in a wide range of $`T`$ confirms that our system is described by the equilibrium density matrix with the Hamiltonian (1). In other words, it is in a mixture of entangled two-qubit eigenstates.

A measure of the entanglement of an arbitrary state is the concurrence $`𝒞`$ Wootters , changing monotonically from 0 for non-entangled to 1 for maximally entangled states. For a normalized pure state $`\alpha |00+\beta |01+\gamma |10+\delta |11`$, one has $`𝒞=2|\alpha \delta \beta \gamma |`$; in general, $`𝒞`$ depends both on the existence of entangled eigenstates and on their population, and can be erased by temperature faster than quantum coherence in individual qubits Yu\_Eberly . Using the experimentally established parameters of the Hamiltonian (1), the concurrences of its eigenstates at the two-qubit degeneracy point are immediately calculated: $`𝒞_1=𝒞_4=0.39`$, $`𝒞_2=𝒞_3=0.97`$. For the equilibrium density matrix of our system at 10 mK, these give $`𝒞_{\mathrm{eq}}=0.33`$. This value will decrease with rising $`T`$, vanishing around 30 mK. However, from Eq. (6) it can be seen that $`\mathrm{\Theta }`$ depends on $`T`$ much more weakly (through the Boltzmann factors $`\rho _\mu `$), so that the IMT deficit can be reliably extracted even at temperatures where $`𝒞_{\mathrm{eq}}`$ is small or even zero.

We stress that for quantum computing, it is the existence of entangled eigenstates that matters, not the measure of equilibrium entanglement \[cf. the discussion a bit below Eq. (6)\], since the system anyway must operate on a faster scale than the equilibration time. The crucial requirement is that the decoherence time $`\tau _{2\mathrm{q}}`$ of the entangled eigenstates exceeds the operation time. The presence of an IMT deficit shows that in our system $`\tau _{2\mathrm{q}}`$ is larger than $`\mathrm{min}(\mathrm{}/\mathrm{\Delta }_{a/b},\mathrm{}/J)1`$ ns; its actual value can be determined using, e.g., the approach of Rabi .

In conclusion, we have investigated two inductively coupled aluminum flux qubits with independently controlled bias fluxes. The impedance measurement technique directly measures the low-frequency part of the system’s magnetic susceptibility. The quantitative agreement between theory and experiment confirms that the system is in an equilibrium mixture of entangled states. Besides trying to scale the approach to larger qubit systems, it will next be important to study the system’s coherence properties in detail, and to establish time-domain control enabling actual quantum gate operations.

###### Acknowledgements.

We thank N. Oukhanski, H. Mühlig, and D. Born for technical assistance, U. Hübner and T. May for sample fabrication, and I. Affleck, D. Bonn, M. Franz, W. Hardy, H.E. Hoenig, P.C.E. Stamp, and M.F.H. Steininger for fruitful discussions. MG acknowledges partial support by the Slovak Grant Agency VEGA (Grant No. 1/9177/02).


text:
# Cosmic Ray Production of 6Li by Virialization Shocks in the Early Milky Way

## 1 Introduction

Apart from 7Li, the bulk of the light elements Li, Be and B are believed to arise from nonthermal nuclear reactions induced by cosmic rays (CRs) (Meneguzzi , Audouze & Reeves 1971). In the last decade, extensive observations of LiBeB in population II, metal-poor halo stars (MPHS) have turned up new and unexpected results, spurring controversy as to what type of CR sources and production mechanisms were operating in the halo of the early, forming Galaxy (e.g. Vangioni-Flam, Cassé & Audouze 2000). To date, most models of light element evolution in the early Galaxy have focused on strong shocks driven by supernovae (SNe) as the principal sources of CRs. A general consensus is that Be and B in MPHS mainly originate from the “inverse” spallation process, whereby CR CNO particles are transformed in flight into LiBeB by impinging on ISM H or He atoms (Duncan, Lambert & Lemke 1992). This can be realized if a sizable fraction of the CRs responsible for spallation comprise fresh, CNO-rich SN ejecta (e.g. Vangioni-Flam et al. 2000, Ramaty et al. 2000, Suzuki, Yoshii & Kajino 1999, Suzuki & Yoshii 2001,hereafter SY), as opposed to CRs injected from the average ISM (e.g. Fields & Olive 1999).

The origin of 6Li in MPHS (Hobbs 2000 and references therein) is more mysterious, as current models involving SN CRs face some difficulties. A peculiar aspect of Li is that in addition to spallation, the fusion process of CR $`\alpha `$ particles with ambient He atoms can be effective, and should actually dominate Li production at low metallicities. (Note that while both 7Li and 6Li are synthesized in comparable amounts, the CR-produced 7Li component is generally overwhelmed by the “Spite plateau” from primordial nucleosynthesis in the metallicity range under consideration; e.g. Ryan et al. 2001.) If the CR energy spectrum is taken to be a standard power-law distribution in momentum (§3), one requires a CR injection efficiency much higher than normally inferred to reproduce the 6Li observations, whether the CR composition is metal-enriched or not (Ramaty et al. 2000; SY, see their fig.2). This raises the question of whether there may have been other sources for 6Li.

Suzuki & Inoue (2002; hereafter SI) proposed a new and more natural 6Li production scenario based on a previously unconsidered CR source: CRs accelerated at virialization shocks (VSs), i.e. shocks driven by the gravitational infall and merging of sub-Galactic gas clumps during the hierarchical build-up of structure in the early Galaxy. Such shocks are inevitable consequences in the currently standard theory of hierarchical structure formation (see however §5). We further discuss a number of testable predictions and important implications expected in our scenario for understanding the formation of our Galaxy. The potential relevance to some fundamental issues in cosmology and structure formation on sub-galactic scales is also mentioned.

## 2 Cosmic Ray Sources in the Early Galaxy

SNe are known to release $`E_{\mathrm{SN}}10^{51}`$ erg of kinetic energy in each explosion, driving strong shocks into the ambient medium. These SN shocks are favorable sites for efficient CR acceleration through the first order Fermi mechanism. A plausible value for the injection efficiency $`\xi _{\mathrm{SN}}`$, i.e. the fraction of the SN kinetic energy imparted to CRs, is $`1020\%`$, deduced from comparison of the SN rate and the energy content of CRs currently observed in the Galactic disk. The global energetics of SNe for the early Galactic halo can be estimated from the cosmic star formation history or the total amount of heavy elements ejected by halo SNe (see SI for details of the estimation). Both methods give $`ϵ_{\mathrm{SN}}_{\mathrm{SN}}\mu m_p/M_g0.15`$ keV per particle as the average specific energy input from SNe.

SNe are not the only sources of mechanical energy (and hence CRs through shock acceleration) that may have been active in the early Galaxy. In the framework of the currently successful picture of hierarchical structure formation in the universe, large scale objects are formed through the merging and virialization of smaller subsystems, driven by gravitational forces acting on the dark matter. For each merging hierarchy, shocks should inevitably arise in the associated baryonic gas component, whereby the kinetic energy of infall is dissipated and the gas heated to the virial temperature of the merged halo (e.g. White & Rees 1978). It is quite plausible that such virialization shocks (VSs) also accelerate CRs (e.g. Miniati et al. 2001).

A guide to how structure formation may have proceeded in our Galaxy, particularly for the Galactic halo, may be offered by the recent numerical simulations of Galaxy formation by Bekki & Chiba (2001, hereafter BC; see also Samland & Gerhard 2003, Abadi et al. 2003). It is expected that sub-Galactic structure eventually merge into a single entity in the central region at redshift $`z2`$, whereby the majority of the infall kinetic energy is virialized. The total energy dissipated at this main VS can be evaluated from the virial temperature of the merged system. The virial temperature $`T_v`$ for a halo of total mass $`M_t`$ virializing at redshift $`z2`$ is $`k_BT_v=\mu m_pGM_t/2r_v0.26\mathrm{keV}\left(M_t/3\times 10^{12}\mathrm{M}_{}\right)^{2/3}`$, where $`r_v`$ is the virial radius, and we have assumed a cosmology with $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`h=0.7`$ (see SI). The specific energy dissipated at the VS is thus $`ϵ_{\mathrm{VS}}3k_BT_v/20.4`$ keV per particle, higher than the above estimate for SNe by a factor of $`2.6`$, if we adopt $`M_t=3\times 10^{12}\mathrm{M}_{}`$ (Sakamoto, Chiba & Beers 2003). The CR contribution from VSs compared to SNe should similarly be higher, as the CR injection efficiency $`\xi _{\mathrm{VS}}`$ should not be too different from that for SNe (Miniati et al. 2001). CRs accelerated by VSs should therefore be at least as important as SN CRs, and may well dominate at early epochs.

Another, more speculative but potentially interesting possibility is outflows powered by massive black hole(s), which may have been active in the early Galaxy. However, there are numerous ambiguities with such a picture, and here we will concentrate on CRs from VSs and SNe (see SI for more discussion).

## 3 Model

We employ assumptions and parameters deemed most plausible for the SN CRs. The CR injection efficiency is taken to be $`\xi _{\mathrm{SN}}=0.15`$ for each SN of kinetic energy $`E_{\mathrm{SN}}=10^{51}\mathrm{erg}`$. A standard, single power-law distribution in particle momentum is adopted for the CR spectrum (see SI). The injection spectral index is chosen to be $`\gamma _{\mathrm{SN}}=2.1`$, appropriate for strong SN shocks, and consistent with the source spectrum inferred for present-day CRs. As assumed in SY, the composition of SN CRs is a mixture of SN ejecta containing freshly synthesized CNO and Fe and the ambient ISM swept up by the SN blastwave.

CRs from VSs are markedly different from SN CRs in a number of important ways. First and foremost, VSs do not synthesize fresh CNO nor Fe, so that the composition of these CRs is completely ascertained by the pre-existing ISM. When the ISM is metal-poor, these shocks induce very little Be or B production through inverse spallation, and are only efficient at spawning Li via $`\alpha \alpha `$ fusion. Second, VSs are not necessarily strong ones, particularly for major mergers of systems with comparable masses (Miniati et al. 2001), provided that the pre-shock gas has not cooled significantly. Shock acceleration should then lead to injection indices $`\gamma _{\mathrm{VS}}`$ steeper than the strong shock limit value of 2, which works in favor of Li production (§4). Major merger shocks may possess Mach numbers as low as $`23`$, corresponding to $`\gamma _{\mathrm{VS}}2.53.3`$; $`\gamma _{\mathrm{VS}}=3`$ is generally chosen below. As with SN CRs, we take the spectral shape to be a momentum power-law distribution and the injection efficiency to be $`\xi _{\mathrm{VS}}=0.15`$.

A further distinction from SNe is that the VS CR flux should not entail any direct dependence on the metallicity, which is in fact an obstacle to predictive modeling. While the hierarchical growth of structure with respect to cosmic time may be evaluated in concrete ways using e.g. extended Press-Schechter formalisms, relating this to \[Fe/H\] requires additional knowledge of how the combination of star formation, SN nucleosynthesis and chemical evolution in our Galaxy proceeded with cosmic time. This involves large uncertainties, and is not specified in our chemical evolution model. As a first step, we choose to describe the time evolution of VS CRs in a simple, parameterized way, assuming a ‘step function’ behavior: VS CRs begin to be injected from a certain time $`t_{\mathrm{VS}}`$, maintains a constant flux for a duration $`\tau _{\mathrm{VS}}`$, and then returns to zero. We take the injection duration $`\tau _{\mathrm{VS}}`$ to be roughly the dynamical time of the major merger, $`3\times 10^8`$ yr. The injection is also assumed to be uniform, since the effect of the main VS should be global throughout the gas under consideration, unlike SNe. The injected VS CR flux integrated over $`\tau _{\mathrm{VS}}`$ is normalized to the above values for $`\xi _{\mathrm{VS}}`$ and $`ϵ_{\mathrm{VS}}`$. The true evolutionary behavior of VS activity relative to metallicity should actually be probed through future observations of 6Li at low \[Fe/H\] (§4, §5).

After injection by either SNe or VSs, the spectral flux $`F_i(E,t)`$ for each CR element $`i`$ evolves with time during subsequent interstellar propagation. This is obtained from time-dependent solutions of the CR transport equation for a leaky box propagation model (SY). Using the transported spectra, we calculate the CR production of LiBeB in the ISM including all three types of reactions: forward spallation of ISM CNO by CR protons and $`\alpha `$’s, inverse spallation of CR CNOs by ISM H and He, and the fusion of CR $`\alpha `$’s with ISM He. For more details, consult SY. We do not consider here the potentially complicating effects of stellar depletion, which are highly uncertain at the moment.

## 4 Results

For our calculations, we have selected the following sets of parameters for $`t_{\mathrm{VS}}`$, $`\tau _{\mathrm{VS}}`$ and $`\gamma _{\mathrm{VS}}`$, respectively, labeled models I - IV: I (0.22, 0.1, 3), II (0.22, 0.1, 2), III (0.32, 0.1, 3), and IV (0.1, 0.5, 3), where $`t_{\mathrm{VS}}`$ and $`\tau _{\mathrm{VS}}`$ are in units of Gyr. These were chosen to provide results exemplary of light element production by VS CRs, in contradistinction to that by SN CRs. The evolution of 6Li and Be vs. metallicity calculated for each model until the end of halo chemical evolution (\[Fe/H\]$`1.5`$) is shown in Fig.1, along with the current observational data in MPHS.

We discuss some salient points regarding these results. First, it is confirmed that with our fiducial parameters, production by SN CRs alone (dashed) works very well for the observed Be (and B, not shown), yet falls short of the observed 6Li. In contrast, with reasonable values for $`ϵ_{\mathrm{VS}}`$, $`\xi _{\mathrm{VS}}`$ and $`\gamma _{\mathrm{VS}}`$, production by VS CRs is capable of explaining the current 6Li data quite adequately. This mainly owes to two facts: 1) VSs are more energetic than (or at least as energetic as) SN shocks, as estimated in §2, and 2) VS CRs can generate 6Li at early epochs independently of the metallicity. Regardless of the early evolutionary behavior, identical 6Li abundances are attained at the end of the halo phase for a given $`\gamma _{\mathrm{VS}}`$ (I, III, & IV), since this is determined by the time-integrated CR flux, for which we had assumed a fixed value. Compared to a flat spectral index of $`\gamma _{\mathrm{VS}}=2`$ (II), a steeper one of $`\gamma _{\mathrm{VS}}=3`$, more appropriate for low Mach number VS shocks (§3), results in a larger 6Li yield, by about factor of 3. This is because with a constant total CR energy, a steeper index implies a larger CR flux in the subrelativistic energy range $`E10`$$`100`$ MeV, where the 6Li production cross section peaks (see fig.1 in SY). Note that a conservatively lower specific energy for VSs, $`ϵ_{\mathrm{VS}}0.15`$ keV/particle (i.e. comparable to SNe), can still be consistent with the available 6Li data provided that $`\gamma _{\mathrm{VS}}3`$.

As already stressed, 2) is a consequence of 6Li synthesis being dominated by $`\alpha \alpha `$ fusion, and VSs not creating any new Fe. Depending on the onset time and duration of the VS, the 6Li abundance may potentially reach large values quickly at very low metallicity, which can be followed by a plateau or a very slow rise. On the other hand, SNe unavoidably give forth to freshly synthesized Fe, so a correlation with Fe/H must arise; in fact 6Li/H vs. Fe/H for SN CRs can never be much flatter than linear (e.g. Ramaty et al. 2000). Moreover, since VSs do not eject fresh CNO, they produce very little Be or B through spallation; only a minuscule contribution can appear as the ISM becomes metal-enriched. This may allow an extremely large 6Li/Be ratio at low \[Fe/H\]. Conversely, we see that SN CRs must play an indispensable role in generating the Be and B observed in MPHS.

Independent of the particular evolutionary parameters, the following abundance trends are characteristic of VS CR production and should serve as distinguishing properties of the scenario for future observations. Going from high to low metallicity: a plateau or a very slow decrease in $`\mathrm{log}^6`$Li/H vs. \[Fe/H\], followed by a steeper decline in some range of \[Fe/H\] corresponding to the main epoch of VS; a steady increase in 6Li/Be, possibly up to values exceeding $`100`$, also followed by a downturn. These traits are very distinctive and not expected in SN CR models, for which the slope of $`\mathrm{log}^6`$Li/H - \[Fe/H\] must be $`1`$ or greater, and the 6Li/Be ratio constant at sufficiently low \[Fe/H\]. Distinction from any production processes in the early universe (e.g. Jedamzik 2000) should also be straightforward, as they predict a true plateau down to the lowest \[Fe/H\], in contrast to an eventual decrease for VS CR models. Further, unique diagnostic features are discussed in §5.

## 5 6Li as Fossil Record of Dissipative Processes during Galaxy Formation

A truly unique and intriguing aspect of the VS CR picture is that 6Li in MPHS can be interpreted and utilized as a fossil record of dissipative gas dynamical processes in the early Galaxy. Of particular interest are various correlations expected between the 6Li abundance and the kinematic properties of the stars. On the one hand, 6Li arises as a consequence of gaseous dissipation through gravitationally-driven shocks, and survives to this day as signatures of the dynamical history of hierarchical structure formation in the early Galaxy. On the other, the kinematic characteristics of stars presently observed should reflect the past dynamical state of their parent gas systems, because once stars form, they become collisionless and have long timescales for phase space mixing (e.g. Chiba & Beers 2000). Interesting relationships may then exist among the two observables.

For example, recent intensive studies of the structure and kinematics of MPHS in our Galaxy based on Hipparcos data (Chiba & Beers 2000) have elucidated the detailed characteristics of our Galaxy’s halo, such as its two-component nature: an inner halo which is flattened and rotating, and an outer halo which is spherical and non-rotating. This dichotomy has been suggested to result from differences in the physical processes responsible for their formation, dissipative gas dynamics being crucial for the former, and dissipationless stellar dynamics determining the latter. The numerical simulations of Galaxy formation by BC support this conjecture: the outer halo forms through dissipationless merging of small sub-Galactic clumps that have already turned into stars (c.f. Searle & Zinn, 1978), whereas the inner halo mainly forms through dissipative merging and accretion of larger clumps that are still gas rich (c.f. Eggen, Lynden-Bell & Sandage, 1962). In our scenario, 6Li production is a direct outcome of the principal gas dissipation mechanism of gravitational shock heating. If the above inferences on the formation of halo structure are correct, 6Li should be systematically more abundant in stars belonging to the inner halo compared to those of the outer halo, a clearly testable prediction. An important prospect is that 6Li may provide a quantitative measure of the effectiveness of gas dynamical processes during formation of halo structure, rather than just the qualitative deductions allowed by kinematic studies.

Another possibility regards the main epoch of VS with respect to metallicity (§3). As already mentioned, the relation between 6Li/H and Fe/H should mirror the time evolution of dissipative energy release through VSs, but is complicated by being convolved with the uncertain ingredients of star formation and chemical evolution. The Fe abundance can be a bad tracer of time, especially at low \[Fe/H\] where the effects of dispersion in SN yields can be extremely large (SY). Stellar kinematics information may offer a handle on this problem. The observed relation between \[Fe/H\] and $`<V_\varphi >`$, the mean azimuthal rotation velocity of MPHS, seems to manifest a distinctive kink around \[Fe/H\] $``$ -2 (Chiba & Beers 2000). Through chemo-dynamical modeling of the early Galaxy, BC have proposed that this kink may correspond to the epoch of the major merger (§3). If this was true, a simple expectation in the context of the VS CR model is that 6Li/H should be just rising near this value of \[Fe/H\], which is the range occupied by the currently 6Li-detected stars; also expected are a steep decline at lower \[Fe/H\], as well as a plateau or slow rise at higher \[Fe/H\] (i.e. close to model III in fig.1). However, any inferences related to \[Fe/H\] are always subject to the chemical evolution ambiguities. A more reliable and quantitative answer may be achieved by looking for correlations between 6Li/H and $`<V_\varphi >`$ without recourse to \[Fe/H\], as 6Li is a direct and pure indicator of dynamical evolution in the early Galaxy.

Thus the VS CR model for 6Li bears important implications for understanding how our Galaxy formed. If the above mentioned trends are indeed observed, it would not only confirm the virialization shock origin of 6Li, but may potentially point to new studies of “6Li Galactic archaeology”, whereby extensive observations of 6Li in MPHS coupled with detailed chemo-dynamical models can be exploited as a robust and clear-cut probe of dissipative dynamical processes that were essential for the formation of the Galaxy.

Furthermore, 6Li in Galactic MPHSs may potentially offer interesting clues to a number of outstanding current problems in cosmology and structure formation theory, all involving physics on sub-galactic scales. 1) The global importance of dynamical feedback by SNe has been a long-standing uncertainty in galaxy formation theory (White & Rees 1978, Abadi et al. 2003). Further observations of 6Li vs. Fe or possibly Be and B at low metallicity and comparison with detailed theoretical models may constrain this crucial unknown. 2) The efficiency and location of VSs as commonly assumed on sub-galactic scales has recently been brought into question (Katz et al. 2002, Birnboim & Dekel 2003), with radical implications for how stars and galaxies form and how the galaxy luminosity function is shaped (Binney 2003). The early evolution of 6Li may provide a direct probe of this presently speculative but important suggestion. 3) Comparison of theoretical simulations with observations of dwarf galaxy cores and of satellite galaxies of our Galaxy and within the Local Group indicate that standard cold dark matter (CDM) produces much more substructure than is actually seen. This “CDM crisis” may point to non-standard dark matter properties, such as warm dark matter, self-interacting dark matter or even more exotic proposals (e.g. Ostriker & Steinhardt 2003, Madau & Kuhlen 2003). On the other hand, the apparent discrepancy may be the result of strong feedback by SNe or a UV background. These different possibilities should result in differences in the early evolution of 6Li, potentially constituting a unique probe. 4) From a combined analysis of measurements of cosmic microwave background anisotropies by WMAP and of the power spectrum on galactic and sub-galactic scales, it has been suggested that the spectrum of primordial density fluctuations deviates from a standard, single power-law (Spergel et al. 2003). Although less drastic than non-standard dark matter, this would also modify the growth of structure on subgalactic scales (Madau & Kuhlen 2003), and the consequent 6Li production at early epochs. More detailed investigations of these intriguing prospects are certainly necessary, but in principle, 6Li in our Galactic halo may shed light on these issues of paramount importance for cosmology and the physics of the early Universe.

## References

Abadi, M. G. et al. 2003, ApJ, 591, 499

Aoki, W. et al. 2003, A&A, submitted

Asplund, M. et al. 2001, in “Cosmic Evolution” Eds. Vangioni-Flam et al., World Scientific

Bekki, K. & Chiba, M. 2001, ApJ, 558, 666 (BC)

Binney, J. 2003, MNRAS, submitted (astro-ph/0308172)

Birnboim, Y., & Dekel, A. 2003, MNRAS, in press (astro-ph/0302161)

Boesgaard et al. 1999, AJ, 117, 1549

Chiba, M. & Beers, T. C. 2000, AJ, 119, 2843

Duncan, D. K., Lambert, D. L., Lemke, M. 1992, ApJ, 401, 584

Eggen, O. J., Lynden-Bell, D. & Sandage, A. R. 1962, ApJ, 136, 748

Fields, B. D. & Olive, K. A. 1999, New Ast., 4, 255

Hobbs, L. M. 2000, Phys. Rep., 333, 449

Jedamzik, K. 2000, Phys.Rev.Lett., 84, 3248

Katz, N. et al. 2002, astro-ph/0209279

Madau, P. & Kuhlen, M. 2003, astro-ph/0303584

Meneguzzi, M., Audouze, J. & Reeves, H. 1971, A&A, 15, 337

Miniati, F. et al. 2001, ApJ, 559, 59

Nissen, P. E. et al. 2000, A&A, 357, L49

Ostriker, J. P. & Steinhardt, P. 2003, Science, 300, 190

Ramaty, R. et al. 2000, ApJ, 534, 747

Ryan, S. G. et al. 2001, ApJ, 549, 55

Sakamoto, T., Chiba, M. & Beers, T. C. 2003, A&A, 397, 899

Samland, M. & Gerhard, O. E. 2003, A&A, 399, 961

Searle, L. & Zinn, R. 1978, ApJ, 225, 357

Smith, V. V., Lambert, D. L., & Nissen, P. E. 1998, ApJ, 506, 405

Spergel, D. N. et al. 2003, ApJS, 148, 175

Suzuki, T. K., Yoshii, Y. & Kajino, T. 1999, ApJ, 522, L125

Suzuki, T. K. & Yoshii, Y. 2001, ApJ, 549, 303 (SY)

Suzuki, T. K. & Inoue, S. 2002, ApJ, 573, 168

Vangioni-Flam, E., Cassé, M. & Audouze, J. 2000, Phys. Rep., 333, 365

White, S. D. M. & Rees, M. J. 1978, MNRAS, 183, 341


text:
# Baryon operators and spectroscopy in lattice QCD

## 1 INTRODUCTION

The determination of the excited baryon spectrum is an important step on the road to understanding the dynamics of QCD, and it is, for example, a vital element of the Jefferson Laboratory experimental program. Lattice gauge theory has a pivotal role not only in providing ab initio calculations of the masses of the lowest lying states, but also in establishing the quark and gluon content of the states. There has thus been increasing activity amongst the lattice community that has demonstrated that the masses of the lowest-lying baryon states of both parities for spin-$`1/2`$ and spin-$`3/2`$ are accessible to lattice calculation, and there have been calculations of the mass of the first positive-parity, spin-$`1/2`$ excitation of the nucleon, the so-called Roper resonance.

Most of these calculations have been made in the quenched approximation to QCD, and an important component of future lattice studies will be a careful analysis of the systematic uncertainties on these results. However, the confrontation of lattice calculations with experiment will require a more complete understanding of the spectrum, obtaining both the masses of the higher spins and of the radial excitations. Furthermore, the discovery of the quark and gluon structure of the excited states will require a more extensive basis of interpolating operators than used so far, allowing both for excited glue and for multi-quark operators. The aim of this talk is to describe the design of such operators, and in particular to construct operators that transform irreducibly under the symmetries of the lattice.

The remainder of this talk is laid out as follows. In the next section, the methodology of computing the spectrum is described, and the symmetries used to classify baryon states are detailed. Gauge-invariant “elemental” operators that have the correct flavor structure and that explore the spatial structure of the baryons are then introduced. The reduction of these operators to ones that transform irreducibly under the cubic symmetry of the lattice is then performed. Finally, the reduction is illustrated for the case of point-like elemental operators obtained using domain-wall fermion (DWF) valence propagators computed on a dynamical asqtad-fermion background.

## 2 THE SPECTRUM FROM LATTICE CALCULATIONS

The computation of the spectrum of states in lattice QCD is in principle straightforward.

1. Choose an interpolating operator $`𝒪`$ that has a good overlap with $`P`$, the state of interest,

$$0𝒪P0,$$

and ideally a small overlap with other states having the same quantum numbers.
2. Form the time-sliced correlation function

$$C(t)=\underset{\stackrel{}{x}}{}𝒪(\stackrel{}{x},t)𝒪^{}(\stackrel{}{0},0).$$
3. Insert a complete set of states between $`𝒪`$ and $`𝒪^{}`$. The time-sliced sum puts the intermediate states at definite momentum, and we find

$`C(t)={\displaystyle \underset{\stackrel{}{x}}{}}{\displaystyle \underset{P}{}}{\displaystyle \frac{d^3k}{(2\pi )^32E(\stackrel{}{k})}}`$
$`\times 0|𝒪(\stackrel{}{x},t)|P(\stackrel{}{k})P(\stackrel{}{k})|𝒪^{}(\stackrel{}{0},0)|0`$
$`={\displaystyle \underset{P}{}}{\displaystyle \frac{0𝒪P^2}{2E_P(\stackrel{}{0})}}e^{iE_P(\stackrel{}{0})t}.`$
where the sum over $`P`$ includes the contributions from two-particle and higher states.
4. Continue to Euclidean space $`tit`$, yielding

$$C(t)=\underset{P}{}\frac{0𝒪P^2}{2m_P}e^{m_Pt},$$
(1)
so that the correlator falls off exponentially with the mass of the lightest state at large times.

In order to obtain the mass of the higher resonances in Eqn. 1, the most appealing approach is to apply variational methods, beginning with the computation of a matrix of correlators

$$C_{ij}(t)=\underset{\stackrel{}{x}}{}𝒪_i(\stackrel{}{x},t)𝒪_j^{}(0),$$
(2)
where the $`\{𝒪_i\}`$ form a basis of operators of definite quantum numbers. In the usual application of the variational technique, we find a basis of eigenvectors $`V(t_0)`$ that diagonalizes the transfer matrix $`C(t_0)^1C(t)`$ for $`t=t_0+1`$ for some $`t_0`$ close to the source, and determine the spectrum from the eigenvalues of

$$V(t_0+1)^1C(t_0)^1C(t)V(t_0+1)$$
(3)
at subsequent times. The method relies on constructing a basis of operators $`\{𝒪_i\}`$ that provides a good description of the states of interest, and the remainder of this talk will focus on the construction of that basis.

We classify the states of the continuum by their flavor structure, $`F`$, and by their parity $`P`$, total spin $`J`$ and “helicity” $`J_z`$. Whilst the flavor structure can be faithfully preserved in a lattice calculation, the parity, spin and helicity labels arise from looking at the symmetry properties of states under rotations and reflections. In a lattice calculation, we have replaced continuum space time by a hypercubic lattice, and thus must classify states according to the irreducible representations $`\mathrm{\Lambda }`$ of the cubic group (for particles at rest); the row $`\lambda `$ within the representation is then the analogue of helicity.

In the following section, we will construct a set of elemental baryon operators $`\{B_i^F(x)\}`$ having the correct flavor properties $`F`$. We will then use group theory to generate a set of baryon operators $`\{B_i^{F\mathrm{\Lambda }\lambda }(x)\}`$ that transform irreducibly under the rotation and reflection symmetries of the lattice, and from which we will construct our correlation matrix

$$C_{ij}^{\mathrm{\Lambda }\lambda F}(t)=\underset{\stackrel{}{x}}{}B_i^{\mathrm{\Lambda }\lambda F}(x)\overline{B}_j^{\mathrm{\Lambda }\lambda F}(0),$$
(4)
and hence extract the spectrum of states corresponding to $`\mathrm{\Lambda },\lambda `$.

## 3 ELEMENTAL BARYON OPERATORS

The color and flavor structure of the operators is dictated by the requirements of gauge invariance and of isospin respectively; it is natural in a lattice calculation to assume exact isospin symmetry, and to classify states and operators by their isospin and strangeness, rather than according to $`SU(3)`$ flavor.

The use of “smearing” has proven essential to reduce the coupling of operators to the higher excited states, and thus to ensure that the correlators Eqn. 2 are dominated by only a few states even close to the source. In the example given here, we will adopt a gauge-covariant spatial Laplacian to smear the fields:

$$\stackrel{~}{\psi }(x)=(1+\sigma ^2\stackrel{~}{\mathrm{\Delta }}/4N)^N\psi (x)$$
(5)
with the three-dimensional Laplacian defined by

$$\stackrel{~}{\mathrm{\Delta }}\psi (x)=\underset{k=\pm 1,\pm 2,\pm 3}{}\left(\stackrel{~}{U}_k(x)\psi (x+\widehat{k})\psi (x)\right),$$
(6)
where $`\stackrel{~}{U}`$ denotes a link variable that can be smeared according to, say, the APE prescription. Both $`N`$ and $`\sigma `$ are tunable parameters, and in the limit $`N\mathrm{}`$, Eqn. 5 reduces to Gaussian smearing of width $`\sigma `$. Note that the square of the smeared field vanishes, in the manner of a simple Grassmann field.

An important consideration is to construct operators that enable us to build up the radial and orbital structure of the states. This is accomplished by additional powers of the Laplacian, and by spatial derivatives, respectively, so that the quark building blocks become:

$`\chi _{Aa\alpha }^n(x)`$ $``$ $`\left(\stackrel{~}{\mathrm{\Delta }}^n\stackrel{~}{\psi }(x)\right)_{Aa\alpha },`$ (7)
$`\xi _{Aa\alpha }^{npj}(x)`$ $``$ $`\left(\stackrel{~}{D}_j^{(p)}\stackrel{~}{\mathrm{\Delta }}^n\stackrel{~}{\psi }(x)\right)_{Aa\alpha },`$ (8)
where $`A`$ is a flavor index, $`a`$ is a color index, $`\alpha `$ a Dirac spin index, and $`n`$ is an integer. The $`p`$-link gauge-covariant forward displacement operator is defined by

$`\stackrel{~}{D}_j^{(p)}O(x)=`$ (9)
$`\stackrel{~}{U}_j(x)\mathrm{}\stackrel{~}{U}_j(x+(p1)\widehat{j})O(x+p\widehat{j}),`$
where $`j=1,2,3`$, and $`p`$ is an integer specifying the length of the displacement.

From these quark building blocks, we construct our elemental three-quark operators:

$`\varphi _{ABC}^F\epsilon _{abc}\chi _{Aa\alpha }^{n_1}\chi _{Bb\beta }^{n_2}\chi _{Cc\gamma }^{n_3},`$ (10)
$`\varphi _{ABC}^F\epsilon _{abc}\chi _{Aa\alpha }^{n_1}\chi _{Bb\beta }^{n_2}\xi _{Cc\gamma }^{n_3pj},`$ (11)
$`\varphi _{ABC}^F\epsilon _{abc}\chi _{Aa\alpha }^{n_1}\xi _{Bb\beta }^{n_2p_1j}\xi _{Cc\gamma }^{n_3p_2k}.`$ (12)
The flavor structure is specified by $`\varphi _{ABC}^F`$, whilst the Levi-Civita symbol ensures that these are color singlets. The radial and orbital structures are explored by varying $`n_1,n_2,n_3`$, and through the choices of $`p_1,p_2,j,k`$ respectively. A large number of baryons can then be studied using a somewhat small number of quark propagator sources.

## 4 SYMMETRIES OF THE LATTICE

Proper rotations restricted to an isotropic cubic lattice form the cubic group $`O`$. It has 24 elements, and five conjugacy classes and thus five single-valued irreducible representations: $`A_1,A_2,E,T_1,T_2`$, of dimensions $`1,1,2,3\text{and}3`$ respectively. If we allow for spatial inversions, corresponding to $`P=\pm 1`$, we obtain the group $`O_h`$, and the irreducible representations acquire a further label, $`g`$ or $`u`$, corresponding to positive and negative parity states, respectively. The construction of operators transforming irreducibly according to the single-valued representations is well known, and has been crucial to identifying the states in lattice calculations of the glueball spectrum.

There are four two-dimensional spinorial representations of $`O_h`$: $`G_{1g},G_{1u},G_{2g}\text{and}G_{2u}`$, and two four-dimensional representatons, $`H_g\text{and}H_u`$. The irreducible representations $`J`$ of the continuum group $`SU(2)`$ are reducible under the cubic group $`O`$; the number of times $`n_\mathrm{\Gamma }^J`$ that each of these reducible representations occurs in the irreducible representation $`\mathrm{\Gamma }`$ of $`O`$ is shown in Table 1; the extension to $`O_h`$ is straightforward.

Whilst full rotational symmetry is restored in the continuum limit, the irreducible representations of $`O_h`$ contain states of many representations of $`SU(2)`$, and in general the degrees of freedom, corresponding to different helicities, lie in different irreducible representations. Thus, for example, a state of spin $`5/2`$ has four degrees of freedom in the four-dimensional representation $`H`$, and two degrees of freedom in the two-dimensional representation $`G_2`$; the identification of the spin of a state is accomplished by looking for the approach to a common mass across the different irreducible representations in the continuum limit.

Central to the task of finding the $`B_i^{\mathrm{\Lambda }\lambda F}(x)`$ is the projection formula

$$B_i^{\mathrm{\Lambda }\lambda F}(\stackrel{}{x})=\frac{d_\mathrm{\Lambda }}{g_{O_h}}\underset{RO_h}{}D_{\lambda \lambda }^{(\mathrm{\Lambda })}(R)U_RB_i^F(\stackrel{}{x})U_R^{},$$
(13)
where $`\mathrm{\Lambda }`$ refers to an $`O_h`$ irrep, $`\lambda `$ is the irrep row, $`g_{O_h}`$ is the number of elements in $`O_h`$, $`d_\mathrm{\Lambda }`$ is the dimension of the $`\mathrm{\Lambda }`$ irrep, $`D_{mn}^{(\mathrm{\Lambda })}(R)`$ is a $`\mathrm{\Lambda }`$ representation matrix corresponding to group element $`R`$, and $`U_R`$ is the quantum operator which implements the symmetry operations; the temporal argument is suppressed.

Application of this formula requires explicit representation matrices for every group element. Representation matrices for all allowed proper rotations can be generated from the matrices for $`C_{4y}`$ and $`C_{4z}`$, the rotations by $`\pi /2`$ about the $`y`$\- and $`z`$-axes, respectively. The explicit matrices we use for these generators are given by

$$D^{(G_1)}(C_{4y})=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}\hfill 1& \hfill 1\\ \hfill 1& \hfill 1\end{array}\right]=D^{(G_2)}(C_{4y}),$$

$$D^{(G_1)}(C_{4z})=\frac{1}{\sqrt{2}}[\begin{array}{cc}1i& 0\\ 0& 1+i\end{array}]=D^{(G_2)}(C_{4z}),$$

$$D^{(H)}(C_{4y})=\frac{1}{2\sqrt{2}}\left[\begin{array}{cccc}\hfill 1& \hfill \sqrt{3}& \hfill \sqrt{3}& \hfill 1\\ \hfill \sqrt{3}& \hfill 1& \hfill 1& \hfill \sqrt{3}\\ \hfill \sqrt{3}& \hfill 1& \hfill 1& \hfill \sqrt{3}\\ \hfill 1& \hfill \sqrt{3}& \hfill \sqrt{3}& \hfill 1\end{array}\right],$$

$$D^{(H)}(C_{4z})=\frac{1}{\sqrt{2}}\left[\begin{array}{cccc}1i& 0& 0& 0\\ 0& 1i& 0& 0\\ 0& 0& 1+i& 0\\ 0& 0& 0& 1+i\end{array}\right].$$

## 5 EXAMPLE: LOCAL NUCLEON INTERPOLATING OPERATOR

We will illustrate the procedure above by examining the reduction of the local proton interpolating operator, with $`I=1/2,I_3=1/2`$. The elemental operators that are gauge invariant and have the correct isospin are

$$\mathrm{\Phi }_{\alpha \beta \gamma }=ϵ^{abc}(u_\alpha ^ad_\gamma ^bu_\beta ^cd_\alpha ^au_\gamma ^bu_\beta ^c).$$
(14)
Examination of Eqn. 14 reveals the constraints $`\mathrm{\Phi }_{\alpha \beta \gamma }+\mathrm{\Phi }_{\gamma \beta \alpha }=0`$ and $`\mathrm{\Phi }_{\alpha \beta \gamma }+\mathrm{\Phi }_{\beta \gamma \alpha }+\mathrm{\Phi }_{\gamma \alpha \beta }=0`$ so that there are only 20 independent operators. The projection operation Eqn. 13 is applied to each of the twenty linearly independent operators in Eqn. 14, and linearly independent operators identified that transform irreducibly under $`O_h`$. The procedure depends on the Dirac basis adopted, and therefore the reduction must be performed using the basis employed to compute the quark propagators. The reduction in the DeGrand-Rossi basis is shown in Table 2; the corresponding reduction in the Dirac basis is provided in ref. .

There are three embeddings of $`G_1`$ and a single embedding of $`H`$, of both parities. This is in accord with expectations from an analysis of the continuum operators, where it is possible to construct three linearly independent nucleon interpolating currents, and a single spin-$`3/2`$ interpolating current.

The computation of the correlator matrix Eqn. 4 is straightforward. We compute the gauge-invariant “generalised” baryon correlator

$`B_{\alpha \beta \gamma ,\overline{\alpha }\overline{\beta }\overline{\gamma }}(t;0)=`$ (15)
$`{\displaystyle \underset{\stackrel{}{x}}{}}ϵ^{abc}ϵ^{\overline{a}\overline{b}\overline{c}}Q_{\alpha \overline{\alpha }}^{a\overline{a}}(x;0)Q_{\beta \overline{\beta }}^{b\overline{b}}(x;0)Q_{\gamma \overline{\gamma }}^{c\overline{c}}(x;0)`$
where $`Q`$ is the quark propagator, whence

$$\underset{\stackrel{}{x}}{}\mathrm{\Phi }_{\alpha \beta \gamma }(\stackrel{}{x},t)\overline{\mathrm{\Phi }}_{\overline{\alpha }\overline{\beta }\overline{\gamma }}=\stackrel{~}{B}_{\alpha \beta \gamma ,\overline{\alpha }\overline{\beta }\overline{\gamma }},$$
(16)
where $`\stackrel{~}{B}`$ is obtained from $`B`$ by symmetrizing over $`\alpha `$ and $`\beta `$ to account for the two possible Wick contractions, and antisymmetrizing over $`\alpha ,\gamma `$ and over $`\overline{\alpha },\overline{\gamma }`$ to take account of the two terms in Eqn. 14. Finally, we form the correlators

$$C_{ij}(t)=\underset{\stackrel{}{x}}{}\mathrm{\Gamma }_i^{\alpha \beta \gamma }\overline{\mathrm{\Gamma }}_j^{\overline{\alpha }\overline{\beta }\overline{\gamma }}\mathrm{\Phi }_{\alpha \beta \gamma }(\stackrel{}{x},t)\overline{\mathrm{\Phi }}_{\overline{\alpha }\overline{\beta }\overline{\gamma }}(0),$$
(17)
where the $`\mathrm{\Gamma }_i,\overline{\mathrm{\Gamma }}_j`$ are given in Table 2. While the computation of the generalized baryon correlators is computationally quite demanding, the subsequent calculation of the $`20\times 20`$ matrix of correlators of Eqn. 17 can be performed in a matter of minutes on a workstation, even for several hundred configurations.

It is instructive to list some of the properties of the correlation matrix Eqn. 17:

* Cross-correlations between operators within different irreducible representations vanish.
* Cross-correlations between the operators of different rows within the same representation vanish.
* The different rows within a representation can be chosen to have the same normalisation, so that correlators thus constructed have the same expectation values.
* Cross-correlations between the same row of different embeddings of the same irreducible representation do not vanish; they can be constructed to vanish only on a particular time slice.

Thus the strategy for extracting the mass spectrum from the data is now clear. Firstly, choose a particular irreducible representation, possibly averaging over the rows within that representation. Then form the matrix of correlators between all the averaged operators within that representation, and extract the spectrum using, say, the method of Eqn. 3.

The talk concludes with an examination of the correlators on an ensemble of lattice data, and in particular a verification of the properties above. To perform this study, we employ an ensemble of around 300 configurations computed on a $`20^3\times 64`$ lattice with $`N_f=2+1`$ flavors of asqtad staggered quarks having $`am_u/am_s=0.01/0.05`$. The lattice spacing $`a=0.124\mathrm{fm}`$ is obtained from the 1P-1S splitting in the Upsilon spectrum.

We employ domain-wall fermions at a single mass $`m=0.01`$ for the valence quarks, which are computed on “chopped” $`20^3\times 32`$ lattices, with $`N_5=16`$; to reduce the residual mass, the configurations are first blocked using a non-perturbative HYP blocking scheme. Smeared sources are used, using $`\sigma =4.35`$ and $`N=30`$ in Eqn. 5, but local sinks, and thus the correlators are not positive definite; full details of the calculation will appear in a subsequent paper .

The effective masses obtained from the diagonal correlators using the three embeddings of $`G_{1g}`$ and $`G_{1u}`$ are shown in figure 1. The plot reveals a clear signal for both positive- and negative parity states from each of the three embeddings. Note that the three usual baryon interpolating operators, for example in ref. , correspond to linear combinations of the operators used here; the identification of a “good” baryon operator corresponding to the nucleon, and a baryon operator having a good overlap onto the Roper resonance, should appear through the variational analysis. The effective masses for $`H_g`$ and $`H_u`$ are shown as the lower plot in figure 1; there is a clear separation from the mass of the ground-state nucleon, with the negative-parity state being lower in mass, as observed experimentally.

We illustrate the orthogonality properties in two ways in figure 2: firstly, by confirming that the cross-correlators between different irreducible representations are indeed zero, and secondly by verifying that those between different rows of the same irreducible representation vanish. The orthogonality properties are clearly well satisfied in these representative examples, indicating that the numerical implementation is under control.

## 6 CONCLUSIONS

We have detailed the design of operators needed to determine the excited baryon spectrum. The method used to construct operators transforming irreducibly under the symmetries of the lattice has been outlined, and the efficacy of the method illustrated. The calculation of “generalised” baryon correlators followed by group-theory projections is computationally efficient, and allows straightforward identification of the quantum numbers of the correlators. Furthermore, the methodology is readily extendible to more complicated operators, such as those for pentaquarks and those having excited glue.

This work was supported by the U.S. National Science Foundation under Awards PHY-0099450 and PHY-0300065, and by the U.S. Department of Energy under contracts DE-AC05-84ER40150 and DE-FG02-93ER-40762. We are indebted to members of the MILC collaboration both for making available their configurations and for many valuable insights, and to Anna Hasenfratz and John Negele for many helpful discussions. Computations were performed on the 128-node Pentium IV cluster at JLab, and at ORNL, under the auspices of the US DOE’s SciDAC Initiative.