Group-theoretical construction of extended baryon operators in lattice QCD

The design and implementation of large sets of spatially extended, gauge-invariant operators for use in determining the spectrum of baryons in lattice QCD computations are described. Group-theoretical projections onto the irreducible representations of the symmetry group of a cubic spatial lattice are used in all isospin channels. The operators are constructed to maximize overlaps with the low-lying states of interest, while minimizing the number of sources needed in computing the required quark propagators. Issues related to the identiﬁcation of the spin quantum numbers of the states in the continuum limit are addressed.


I. INTRODUCTION
Spectroscopy is a powerful tool for uncovering the important degrees of freedom of a physical system and the interaction forces between them. The spectrum of quantum chromodynamics (QCD) is indeed very rich: conventional baryons (nucleons, , , , , . . . ) and mesons (, K, , , . . . ) have been known for nearly half a century, but other higher-lying exotic states, such as glueballs, four-quark states, and so-called hybrid mesons and hybrid baryons bound by an excited gluon field, have proved more elusive, mainly because our theoretical understanding of such states is insufficient, making their identification problematical.
Recent discoveries of several new hadronic resonances have generated much excitement in the field of hadron spectroscopy. The E852 collaboration [1] reported a signal for a 1 ÿ hybrid meson decaying into with a mass near 1.6 GeV, though the significance of this result has been questioned [2]. Another exotic 1 ÿ candidate at 1.4 GeV has been tentatively identified in the channel by E852 [3], VES [4], KEK [5], and Crystal Barrel [6]. The first observation of a doubly charmed baryon has been reported [7], and evidence for the possible existence of a strangeness S 1 pentaquark state has emerged at SPring-8, JLab, and elsewhere [8][9][10][11]. Interest in excited baryon resonances has also been raised by experiments dedicated to mapping out the N spectrum in Hall B at the Thomas Jefferson National Accelerator Facility (JLab), and the search for hybrid mesons and glueballs is intensifying due to the Hall D initiative at JLab and experiments at CLEO-c.
Unfortunately, our understanding of many conventional and excited hadron resonances is vestigial and comes only from QCD-inspired phenomenological models, such as the bag model, the nonrelativistic quark model, and quarkdiquark models, or from approaches, such as QCD sum rules and methods based on Schwinger-Dyson equations, which use approximations whose justifications are unclear. There are a growing number of resonances which cannot be easily accommodated within quark models. States bound by an excited gluon field, such as hybrid mesons and baryons, are still poorly understood. The natures of the Roper resonance and the 1405 remain controversial. Experiment shows that the first excited positive-parity spin-1=2 baryon lies below the lowest-lying negativeparity spin-1=2 resonance, a fact which is difficult to reconcile in quark models. Given the current surge in experimental activity, the need for an understanding of such states from QCD itself has never been greater.
This need has motivated us to undertake a comprehensive ab initio study of the hadron spectrum in QCD. Presently, Monte Carlo estimates of QCD path integrals defined on a space-time lattice offer the best way to make progress in this regard, so this is the calculational approach we have adopted. The Monte Carlo approach has been used to investigate hadrons throughout the past two decades, but the number of states studied to date has been somewhat limited. Also, prior works have not strived to identify the continuum spin J of the states studied, simply assuming that the lowest allowed J would be the lowest-lying state. The novel features of our approach are its comprehensiveness and its use of techniques to identify spin. Given the vast amount of experimental data being generated at accelerator facilities, such as Jefferson Lab, there is an urgency to investigate the spectrum (masses, widths, transition rates, form factors, and so on) as completely as possible. The number of hadron eigenstates which can be reliably extracted in Monte Carlo computations is not currently known, so our undertaking will be partly an exploration of the limits of possibility. Another aim is to discover whether a very large number of interpolating operators are needed to extract the low-lying spectrum, or whether a handful of carefully chosen ones is sufficient, and this work outlines a systematic means of finding such operators.
Our first goal is to calculate the masses of as many lowlying hadron resonances as possible in QCD using Monte Carlo techniques. The masses and widths of resonances (unstable hadrons) cannot be calculated directly in finitevolume Monte Carlo computations, but must be deduced from the discrete spectrum of finite-volume stationary states for a range of box sizes [12 -15]. The rigorous application of such techniques to obtain the resonance parameters to high accuracy would require vast computational resources, but our goal here is merely to obtain a first exploratory scan of the spectrum of QCD, not to pin down each mass to very high precision. Hence, simply obtaining the finite-volume spectrum for a few judiciously chosen volumes should suffice for ferreting out the hadron resonances from the less interesting scattering states and may even give qualitative information about preferred decay modes.
To compute the finite-volume stationary-state energies, the temporal correlations C ij t h0jTO i tO j 0j0i, where T denotes time-ordering, of a set of operators O j 0 which create the states of interest at an initial time t 0 with a corresponding set of operators O i t which annihilate the states of interest at a later time t must be determined. The correlation functions C ij t can be expressed in terms of path integrals over the quark and gluon fields, and when formulated on a Euclidean space-time lattice, such path integrals can be estimated using the Monte Carlo method with Markov-chain importance sampling. Incorporating the quark-field effects into the Monte Carlo updating for realistically light quark masses remains a challenge, but there is steady progress with improving algorithms and increasing computational power.
The procedure for extracting the lowest stationary-state energies E 0 , E 1 , E 2 , . . . from the Hermitian matrix of correlation functions C ij t is well known [16,17]. Let n t; t 0 denote the eigenvalues of the Hermitian matrix Ct 0 ÿ1=2 CtCt 0 ÿ1=2 , where t 0 is some fixed reference time (typically small) and the eigenvalues, also known as the principal correlation functions, are ordered such that 0 1 as t becomes large. Then one can show that Determinations of the principal correlators n t; t 0 for sufficiently large temporal separations t yield the desired energies E n . The above equations illustrate the difficulties which must be faced in order to extract the stationary-state energies E 0 , E 1 , . . . from the temporal correlations of the hadronic operators. In each principal correlator, there are contaminating contributions from all other states which can be created and annihilated by the operators used. In order to reliably extract the single decaying exponential of interest, the obscuring contributions from all of these other states must somehow be suppressed.
There are two crucial ingredients in reducing the unwanted contributions in n t; t 0 . The first is to use sufficiently large values of t. However, there are often practical considerations which limit how large t can be, and the statistical uncertainties in the Monte Carlo estimates of the correlation functions generally increase with t. The second, and more important, consideration in suppressing the contamination in n t; t 0 is to use cleverly devised operators which couple minimally to the unwanted states.
Successfully extracting the spectrum of QCD in our Monte Carlo computations will hinge crucially on using carefully designed hadronic operators. Excited meson and baryon resonances are expected to be large objects, so the use of spatially extended operators is important. Since our calculations will be carried out on hypercubic space-time lattices, the energies of states in all irreducible representations (irreps) of the O h cubic point group must be determined in order to identify the continuum-limit spin J of each physical state. Because determining the mass of a particular resonance requires determining the energies of all lower-lying stationary states, including scattering states, the set of operators we use must include not only meson and baryon operators, but also multihadron operators. Another very important fact to keep in mind is the computational cost of evaluating quark propagators, especially for light quark masses. Hence, our operators must be devised with an eye towards minimizing the number of sources needed to calculate the required quark propagators.
Designing the hadronic operators is an important first step in our comprehensive study of the mass spectrum in lattice QCD. From the above considerations, the guiding principles in devising our operators are maximizing overlaps with the states of interest while minimizing the number of quark-propagator sources. Although operators for baryons, mesons, and their scattering states will be needed, we restrict our attention to the construction of the threequark baryon operators in this first paper. Because of the complexity of these calculations and the importance of providing checks on our final results, we have been pursuing two different approaches to constructing the baryon interpolating field operators. Only one approach is described here; an alternative approach based on Clebsch-Gordan techniques will appear elsewhere [18]. Meson and multihadron operators will be detailed in subsequent works. Furthermore, this paper deals only with issues related to the construction and utilization of these operators. Results from Monte Carlo calculations using these operators will be presented in later publications. This paper is organized as follows. An overview of our approach to constructing the hadronic operators is first outlined in Sec. II. Our operators are assemblages of basic building blocks, which are described in Sec. III, along with the conventions we use for the Dirac -matrices and in Wick rotating into Euclidean space-time. We use quark fields which are Dirac spinors, so our hadron operators apply to Monte Carlo calculations involving Wilson, domain-wall, and overlap fermion actions, but not to computations involving staggered fermions. The basic building blocks are then combined into gauge-invariant three-quark operators (referred to as elemental operators) having appropriate flavor structure in Sec. IV. Projections onto the rotation-reflection symmetry sectors produce the final operators in Sec. V. Issues related to the use of these projections in constructing the baryon propagators are discussed in Sec. VI. Concluding remarks and plans for future work are outlined in Sec. VII.

II. OVERVIEW OF OPERATOR CONSTRUCTION
Devising relativistic hadronic operators in continuous space-time usually involves combining Dirac spinors to form Lorentz scalars, pseudoscalars, vectors, axial-vectors, and so on, using the Dirac -matrices and the charge conjugation matrix C satisfying C C ÿ1 ÿ T . It has been common practice in lattice QCD simulations to use hadron operators built in a similar fashion through simply discretizing their continuum analogs. However, such an approach becomes very cumbersome when constructing higher spin states or complicated extended operators. Also, the above operators generally couple to states belonging to different J P (spin-parity) symmetry sectors. Since the hypercubic lattice breaks Lorentz covariance, there is really no reason to construct operators according to Lorentz symmetries. Finally, we wish to extract a large portion of the low-lying spectrum, which means that large sets of operators will be needed to compute complete correlation matrices. Hence, the usual approach which mimics that used in continuous space-time is not feasible for our purposes.
Instead, we advocate an approach which more directly combines the physical characteristics of baryons with the symmetries of the lattice regularization of QCD. Recall that baryon states are characterized by their total momentum p, their total (half-integral) spin J, a projection of this spin onto some axis (the z-axis or the momentum, say), their parity P 1, and their quark flavor content. The masses of the light u and d quarks are very nearly equal, and therefore, we work in the approximation that m u m d . In this approximation, the theory has an exact isotopic spin symmetry, and states carry two more labels, total isospin I, and its projection I 3 onto a given axis. The other flavor quantum numbers which label the states are strangeness S, charm C, and bottomness B (we do not consider the top quark). In our calculations, isospin remains an exact symmetry since we neglect electromagnetic interactions.
Since our simulations will be performed using a hypercubic space-time lattice, our operators should be classified according to the symmetries of the lattice, rather than the full rotational symmetries of continuous space-time. Of course, we expect to recover the symmetries of the space-time continuum as the lattice spacing is made small. Since we are interested only in determining the masses of the baryon states, we restrict our attention to representations corresponding to zero total three-momentum p 0. Hence, our operators must be invariant under all allowed spatial translations and we require that they transform under spatial rotation-reflection symmetry operations according to the irreducible representations of the octahedral point group O h . These irreducible representations are the lattice analogs of the continuum J P labels, and the row of the representation is the analog of the spin projection . Thus, our (annihilation) operators can be written where indicates the irreducible representation of O h , is the row of the representation, F denotes all of the quantum numbers associated with the flavor content of the operator, and i labels the different operators in the F symmetry sector. Under a symmetry operation R, these operators transform according to where U R denotes the quantum operator which effects the symmetry operation corresponding to group element R (not to be confused with the gauge link variables), and ÿ R are the elements of the representation matrix corresponding to group element R. For baryons, we are only interested in states corresponding to half-integral spin J in the continuum limit, so we can restrict our attention to the six spinorial representations of O h . There are four twodimensional irreducible representations G 1g , G 1u , G 2g , and G 2u (adopting a Mulliken-like naming convention), and two four-dimensional representations H g and H u . These representations will be discussed in greater detail later.
Our general approach to constructing the B F i t operators is to (1) first identify appropriate basic ''building blocks'' to use in constructing all baryon operators, (2) devise simple elemental operators containing the appropriate flavor and color structure, then (3) apply appropriate group-theoretical projection operators to find linear combinations of the elemental operators with the desired transformation properties under the symmetry group of a spatial cubic lattice. Let B F i t P x B F i x; t denote a gaugeinvariant elemental operator with the appropriate quark flavor content and which is invariant under allowed spatial translations, then an operator which transforms according to the row of the irreducible representation is obtained Note that the expansion coefficients in the B F i operators are the complex conjugates of those in

III. THE BASIC BUILDING BLOCKS
The oscillating path integral weight e iS M in quantum field theory, where S M is the action defined in Minkowski space-time and using natural units @ c 1, is unsuitable for applying the Monte Carlo method to evaluate the correlation functions of the theory via Feynman path integrals. A rotation to imaginary time t ! ÿi, where is real, leads to path integrals with weight e ÿS , where S is the action defined in Euclidean space-time. If S is real, the path integral weight is real and positive and, hence, can be interpreted as a probability, allowing the application of Monte Carlo methods with suitable importance sampling. The Euclidean action S is defined such that S is invariant under all symmetries of Euclidean space-time and all Green's functions of the theory are identical to the Green's functions of the Minkowski theory, analytically continued to imaginary time t ! ÿi. Although our simulations employ a path integral quantization of the field theory, a canonical quantization viewpoint can be adopted when discussing the quantum operators.
Our conventions for the continuation from Minkowski space-time with metric g diag1; ÿ1; ÿ1; ÿ1 into Euclidean space-time (imaginary time) are as follows. We define the following Euclidean space-time coordinates and derivatives (a subscript or superscript M indicates a Minkowski space-time quantity): @ 4 @ 4 ÿi@ M 0 ; @ j @ j ÿ@ j M @ M j ; (10) for spatial directions j 1, 2, 3. The metric in Euclidean space-time is so there is no distinction between covariant and contravariant indices. Our Monte Carlo calculations will be carried out using a hypercubic space-time lattice, and we require that the lattice spacings in the three spatial directions are the same, denoted by a s ; the temporal lattice spacing a t may differ from a s , allowing us to exploit the known benefits of anisotropic lattices [19]. Throughout this paper, we set a s 1 to simplify the notation. Four vectors of unit length pointing along the spatial axes of the lattice with a vanishing temporal component will be denoted byĵ,k, and so on, for j; k 1; 2; 3.
As usual in lattice gauge theory, the gluon field is introduced using the parallel transporter U x given by the path-ordered exponential of the gauge field along each link connecting neighboring sites of the hypercubic lattice. We also introduce the Dirac spinor field A a x which annihilates a quark and creates an antiquark, where A refers to the quark flavor, a refers to color, and is the Dirac spin index, and the field A a x which annihilates an antiquark and creates a quark. Unlike in Minkowski space-time, and must be treated as independent fields, so we emphasize that Þ y 4 . This is required in order to simultaneously satisfy Euclidean covariance of the fields, the canonical anticommutation relations, and the equality of the Euclidean two-point function with the relativistic Feynman propagator continued to imaginary times [20,21]. Our Euclidean space-time Dirac-matrices are related to their Minkowski counterparts by Throughout this paper, we use the standard Dirac-Pauli representation for the -matrices: where the Pauli spin matrices are given by It has long been known that operators constructed out of smeared fields have dramatically reduced mixings with the high frequency modes of the theory. Thus, our operators are constructed using spatially smoothed link variables Ũ j x and spatially-smeared quark fields~ x. The spatial links are smeared using either the stout-link procedure described in Ref. [22] or the method introduced in Ref. [23]. Note that only spatial staples are used in the link smoothening; no temporal staples are used, and the temporal link variables are not smeared. The smeared quark fields are defined by [24] x where s and n are tunable parameters (n is a positive integer) and the three-dimensional covariant Laplacian operators are defined in terms of the smeared link variables U j x as follows: where Ox and Ox are operators defined at lattice site x with appropriate color structure, and noting thatŨ ÿk x U y k x ÿk. The smeared fields~ and~ are Grassmannvalued; in particular, these fields anticommute in the same way that the original fields do, and the square of each smeared field vanishes.
Our operators are designed with one eye on capturing the states of interest and the other eye on facilitating the efficient computation of the hadron correlation functions.
Since the baryon resonances are expected to be large objects, the use of extended operators is crucial. Our choice of basic building blocks is motivated by the need to incorporate spatial extensions so as to facilitate the efficient and gauge-invariant assembly of many large hadron operators. To capture orbital structure, the quarks must be displaced in different directions, and to capture radial structure, the quarks must be displaced by different distances. To maintain gauge invariance, covariant displacements in terms of the gauge-field parallel transporters must be used. To simplify matters, we consider only straight-path displacements in the six directions along the axes of the spatial cubic lattice: j 1; 2; 3. Hence, we shall assemble our baryon (and later, meson) operators using the following basic building blocks: where A is a flavor index, a is a color index, is a Dirac spin index, and the p-link gauge-covariant displacement operator in the jth direction is defined byD for j 1; 2; 3, and p 1. In what follows, we can achieve a significant economy in the notation by defining a zero-displacement operator to indicate no displacement. In this way, our basic building blocks can be listed more succinctly: Sometimes it will be convenient to make the flavor quantum number more apparent by writing and similarly for the s, c, b quarks. It is important to remember that theũ,d,s, etc. operators so defined refer to smeared quark fields.

IV. THREE-QUARK ELEMENTAL OPERATORS
Having chosen the basic building blocks for our hadron operators listed in Eq. (22), the next step is to devise elemental baryon operators B F i t having the appropriate color and flavor structure. These operators are chosen such that they are gauge-invariant and transform irreducibly under the isotopic spin symmetry. Explicitly dealing with SU3 flavor symmetry is not necessary, as will be discussed below.
Gauge invariance is easily handled. Given three quarks with color indices a, b, c associated with the same lattice site x, there is only one way of combining the color indices to arrive at a locally gauge-invariant object: the use of the antisymmetric Levi-Civita symbol " abc . Similarly, covariantly displaced quark fields must also be connected by an " abc coupling at a single lattice site. To simplify the operator construction, we consider only combinations of displaced quarks having the same displacement length p. Thus, all of our three-quark baryon operators are linear superpositions of gauge-invariant terms of the form The ''barred'' operators have the form Note that the barred composite operator is defined differently than the barred fermion field operator due to the presence of the 4 spin matrices in Eq. (25). Their presence is needed to ensure that the resulting correlation matrices satisfy the desirable property of hermiticity. These 4 matrices do not affect the transformation properties of these operators under proper spatial rotations and parity, although they do affect the Lorentz boost properties.
The simplest way to combine the basic building blocks previously described is to combine three-quark fields at a single lattice site, corresponding to i j k 0 in Eq. (24). We refer to these operators as single-site operators. Next, one of the three quarks can be displaced; such singly displaced operators correspond to i j 0, k Þ 0 in Eq. (24). In these operators, the two quarks which are not displaced may be viewed as forming a localized diquark, so such operators may be important if baryon formation is dominated by a quark-diquark mechanism.
Two of the quarks can be displaced; they can be displaced either in opposite directions (doubly displaced-I), or in orthogonal directions (doubly displaced-L). In such operators, one can in general choose different lengths for the two displacements, but to simplify matters, we restrict our attention to the case in which both displacements have the same length p. Since the displacement of two quark fields essentially results in an object in which all three quarks are at different lattice sites, one may be tempted to exclude these operators in favor of triply displaced operators. However, the relative importance of so-called Y-flux and -flux formation of the gluon field in threequark systems (the is actually a quantum superposition of V-flux forms) is a much-discussed issue (see, for example, Refs. [25][26][27][28][29][30][31]). Hence, it is important to include some operators with significant overlap with a -flux configuration, as well as operators with strong mixing with a Y-flux configuration. The doubly displaced operators above are suitable for -flux formation.
Lastly, all three quarks can be displaced (again, restricting attention to the case of equal distances). If two of the quarks are displaced in opposite directions, this produces a coplanar T-shape (triply displaced-T); alternatively, the three quarks can be displaced in mutually orthogonal directions (triply displaced-O). These operators are suitable for producing Y-flux configurations.
These operators, summarized and illustrated in Table I, allow a large number of baryon operators to be constructed using a relatively small number of quark-propagator sources. For a given reference source site x, quark propagators must be evaluated using only a handful of different sources: for each quark mass value, we need a local source and displaced sources in at most each of the six directions. However, rotational invariance of the baryon correlation functions can be exploited to reduce the number of source displacement directions. For singly displaced operators at the source, a simultaneous rotation of the source and sink can always be used to align the source displacement along the z direction, say. For doubly displaced-I sources, a rotation can always be done to align one displacement along the z direction and the other along the ÿz direction. Doubly displaced-L sources can be rotated so the source displacements are along the y and z directions, triply displaced-T sources can be rotated to align the displacements along the y, z, and ÿz directions, and triply displaced-O sources can always be rotated so the displacements are along the x, y, z directions. In total, only source displacements along the x, y, z, ÿz directions (four directions) are required. Hence, the number of conjugate gradient inversions needed is where N c 3 is the number of colors, N sp 4 is the number of Dirac spin components, N p is the number of displacement lengths p, and N is the number of quark masses to be used. Incorporating the isospin symmetry is also straightforwardly done. Let 1 , 2 , 3 denote the three Hermitian Because of the isospin symmetry of QCD in the m u m d approximation, the particle masses do not depend on I 3 , so we limit our attention to only one value of I 3 in each isospin I, strangeness S sector, choosing the maximal I 3 I value. Using the above isospin relations, we first write down all possible flavor combinations appropriate for each  The linearly independent operators we chose are described in the second column, and the numbers of independent operators of each type are listed in the third column. Replacing each u quark by an s quark in the operators below yields ÿ elemental operators.

Operator type
Restrictions Multiplicity  (24). The linearly independent operators we chose are described in the second column, and the numbers of independent operators of each type are listed in the third column. Interchanging the u and s quarks yields the 0 elemental operators.

Operator type
Restrictions Multiplicity isospin channel for the three-quark elemental operators from Table I. These are listed in Table II. For each flavor sector and quark-displacement type, we then determine a maximal set of linearly independent operators, giving us the final set of elemental operators we use. A symbolic computer program capable of manipulating Grassmann fields was written using MAPLE 9.5 and utilized to identify linearly independent operators. The independent elemental operators we chose in the different flavor sections are described in Tables III, IV Lastly, note that the baryons c , c , cc , and ccc can be studied using the operators presented above if the s quark is replaced with a c quark. Similarly, replacing the s quark by a b quark in the above operators allows us to investigate the b , b , bb , and bbb baryons. Some other baryons, such as cc , can be studied using other suitable flavor replacements in the above operators, whereas the investigations of baryons such as c containing usc quarks cannot directly exploit the above tables, requiring slight modifications.

V. PROJECTIONS ONTO SYMMETRY SECTORS
The final step in our operator construction is to apply group-theoretical projections to obtain operators which transform irreducibly under all lattice rotation and reflection symmetries. The basic building blocks used to assemble our baryon operators transform under the allowed spatial rotations and reflections of the point group O h according to (34) where the transformation matrices for spatial inversion I s and proper rotations C nj through angle 2=n about axis Oj are given by  (24)). The linearly independent operators we chose are described in the second column, and the numbers of independent operators of each type are listed in the third column.

Operator type
Restrictions Multiplicity  (24)). The linearly independent operators we chose are described in the second column, and the numbers of independent operators of each type are listed in the third column.

Operator type
Restrictions Multiplicity with ! kl ÿ2" jkl =n and ! 4k ! k4 0 (! is an antisymmetric tensor which parametrizes rotations and boosts). A rotation by =2 about the y-axis is conventionally denoted by C 4y , and C 4z denotes a rotation by =2 about the z-axis. These particular group elements are given by The allowed rotations on a three-dimensional spatially isotropic cubic lattice form the octahedral group O which has 24 elements. Inclusion of spatial inversion yields the point group O h which has 48 elements occurring in ten conjugacy classes. All elements of O h can be generated from appropriate products of only C 4y , C 4z , and I s . Charge conjugation is another symmetry of our theory. Under charge conjugation C, the link variables U ! U and our basic building blocks transform according to where the charge conjugation matrix C must be unitary, antisymmetric, and satisfy C C y ÿ T . Our choice for C in the Dirac-Pauli representation is C 4 2 .
Operators which transform according to the irreducible representations of O h can then be constructed using the well-known group-theoretical projections given in Eq. (6). Since baryons are fermions, we need only be concerned with the six double-valued irreps of O h . There are four two-dimensional irreps G 1g , G 1u , G 2g , and G 2u , and two four-dimensional irreps H g and H u . The subscript g refers to even-parity states, whereas the subscript u refers to oddparity states. The irreps G 1g and G 1u contain the spin-1=2 states, spin-3=2 states reside in the H g and H u , and two of the spin projections of the spin-5=2 states occur in the G 2g and G 2u irreps, while the remaining four projections reside in the H g and H u irreps. The spin content of each O h irrep in the continuum limit is summarized in Table VII To carry out the projections in Eq. (6), explicit representation matrices are needed. Our choice of representation matrices is summarized in Table VIII. Matrices for only the group elements C 4y , C 4z , and I s are given in Table VIII   TABLE VII. Continuum-limit spin identification: the number n J of times that the irrep of the octahedral point group O h occurs in the (reducible) subduction of the J irrep of SU2. The numbers for G 1u , G 2u , H u are the same as for G 1g , G 2g , H g , respectively. The G 1u , G 2u , H u matrices for the rotations C 4y , C 4z are the same as the G 1g , G 2g , H g matrices, respectively, given below. Each of the G 1g , G 2g , H g matrices for spatial inversion I s is the identity matrix, whereas each of the G 1u , G 2u , H u matrices for I s is ÿ1 times the identity matrix.
The matrices for all other group elements can be obtained from appropriate multiplications of the C 4y , C 4z , and I s matrices. ÿ C 4y ÿ C 4z since the representation matrices for all other group elements can be obtained by suitable multiplications of the matrices for the three generating elements. For baryons, the representation matrix for E in each of the O D h extra irreps is ÿ1 times the identity matrix.
Note that Table VII is the key to identifying the continuum-limit spin J corresponding to the masses extracted in our Monte Carlo calculations. For example, to identify an even-parity baryon as having J 1=2, a level must be observed in the G 1g channel, and there must be no degenerate partners in either of the G 2g or H g . A level observed in the H g channel with no degenerate partners in the G 1g and G 2g channels (in the continuum limit) is a J 3=2 state. Degenerate partners observed in the G 2g and H g channels with no partner in the G 1g channel indicates a J 5=2 baryon. In other words, Table VII details the patterns of continuum-limit degeneracies corresponding to each half-integral J value.
Our operators are constructed using fermion fields x which annihilate a quark and create an antiquark. Hence, each of our baryon operators annihilates a three-quark system of a given parity P and creates a three-antiquark system of the same parity P. This means that in the baryon propagator, a baryon of parity P propagates forward in time while an antibaryon of parity P propagates backwards in time. Unlike boson fields, a fermion and its antifermion have opposite intrinsic parities, so that the antibaryon propagating backwards in time is the antiparticle of the parity partner of the baryon propagating forwards in time. Since chiral symmetry is spontaneously broken, the masses of parity partners may differ. The forward propagating baryon will have a mass different from that of the antibaryon propagating backwards in time. If the even and odd-parity baryon operators are carefully designed with respect to one another, it is possible to arrange a definite relationship between the correlation matrix elements of one parity for t > 0 and the opposite-parity matrix elements for t < 0, allowing an increase in statistics. Our operators are designed to take advantage of this symmetry (see below).
Our construction of the irreducible B F i t baryon operators is done in the following sequence of steps.
(a) A set of M B linearly independent elemental operators B F j t that transform among one another under O D h is identified with the help of the computer software mentioned earlier. This is done by starting with all possible operators of a given type (single-site, singly displaced, and so on), then using the MAPLE program to detect dependencies between the operators. To find such dependencies, the computer program expresses each B F j t operator as a sum of products of Grassmann fields (or gauge-covariantly displaced Grassmann fields) with explicit color, flavor, and spin indices: . Since each O k operator is a single product of three displaced Grassmann fields, these operators are linearly independent, so the linear independence of the B F j operators boils down to the linear independence of the rows of the g F jk matrix. In this way, the set of all possible operators of a given type is easily reduced to a set containing only linearly independent operators.
(b) We obtain the M B M B representation matrices W ij R which satisfy Our MAPLE program determines the ith column of the W ji R matrix for a given symmetry transformation R as follows. First, the B F i t operator is expressed as a sum of products of displaced Grassmann fields with explicit color, flavor, and spin indices as in the previous step: Next, the R symmetry transformation is applied to the displaced Grassmann fields in each O k term in this sum of products using Eq. (33). The resulting sum of terms U R B F i tU y R P k h F ik O k is then expressed as a linear superposition of the original M B operators using the Moore-Penrose pseudoinverse [32] of the g F ik matrix. If the transformed operator contains Grassmann products which are not in the original set of O k operators, then this signals that the starting basis of B F j operators is incomplete, but with our method in the previous step of choosing linearly independent operators, this did not occur. In this way, the W ji R matrices for the generating group elements C 4y , C 4z , and I s are obtained. The matrices for all other group elements are then determined by appropriate products of these three matrices. At this point, we have a set of M B operators B F i t which form the basis of a reducible representation given by the W ij R matrices. Our remaining task is to find a change of basis such that the resulting representation matrices are block-diagonal with the blocks given by the irreducible representations of O h .
(c) Since the resulting W ji R matrices may not be unitary, we compute the Hermitian metric matrix M This matrix will be needed in a later step where it will facilitate the full block-diagonalization of the W ij R matrices.
(d) For each even-parity irrep , we compute the large M B M B projection matrix for row 1: This is one of the most important steps in our operator construction. Applying the group-theoretical projection of Eq. (6) to the operator B F i , then utilizing Eq. (40), produces a new operator B F i t P j P F ij B F j t which resides in the subspace of operators which transform according to the row 1 of the given irrep . In other words, the ith row of the projection matrix P contains the superposition coefficients of the projected B F i operator.
(e) Although group theory guarantees that the resulting projected B F i operators reside in the subspace associated with row of the irrep, it does not guarantee that all of the resulting operators are linearly independent. The maximum number of independent operators in the projected subspace is given by the rank r of the projection matrix. Hence, the next step is to form r superpositions of the B F i operators such that the resulting r operators are linearly independent. The choice of these operators is not unique. In practice, these linear combinations are obtained using the well-known Gram-Schmidt procedure, but with a modified inner product to incorporate the metric matrix M. The use of the metric matrix M ensures full blockdiagonalization of the original W ij R matrices. Using such a procedure, the final operators, expressed in terms of the original set of B F i operators by The use of operators belonging to other rows will be important for increasing the statistics of our Monte Carlo calculations, as will be discussed below.
(g) Although the odd-parity operators can be obtained using the same procedure described above, we instead utilize charge conjugation to construct the odd-parity operators. Consider the correlation matrix element of evenparity operators for t 0. Suppressing flavor and displacement indices, one sees that invariance under charge conjugation implies that using invariance under time translations of the above expectation value and invariance of the vacuum under charge conjugation. The last line above represents the correlation of odd-parity operators propagating temporally backwards. Hence, for a given even-parity operator B g i t, we can define an odd-parity operator B u i t by rotating the three Dirac indices using the 2 matrix and replacing the expansion coefficients by their complex conjugates such that the correlation matrices of the even and odd-parity operators are related by relate operators of different types, such as singly displaced operators with single-site operators, identifying these relationships must be done at a later stage in the calculations.

VI. BARYON PROPAGATORS
To extract the baryon masses, we need to compute the correlations C F ij t h0jTB F i tB F j 0j0i, using the operators constructed as described above. In this section, several issues related to computing these correlation matrix elements are discussed. In particular, we discuss the use of symmetry to minimize the number of quark-propagator sources, the use of gauge-invariant three-quark propagators as an intermediate step in the baryon propagator determinations, and detail the application of Wick's theorem.
The baryon propagators may be expressed in terms of the correlations of the elemental operators by Note that correlations between operators in different rows of the same irrep vanish. Since the number of elemental operators is large and the quark propagators are rather expensive to compute, it is very important to use symmetry to reduce the number of quark-propagator sources. Given the invariance of the vacuum and the unitarity of the symmetry transformation operators, we know that for any group element R of O h . Hence, for each source B F l 0, we can choose a group element R l such that we minimize the total number of source elemental operators which must be considered. For example, consider the singly-displaced operators. We can choose an R l such that the displaced quark in the source is always displaced in the z direction. Similarly, a group element R l can always be chosen to rotate each of the other types of operators into a specific orientation.
The coefficients c F ij in the baryon operators involve only the Dirac spin components and the quarkdisplacement directions and are independent of the color indices and spatial sites. Thus, in calculating the baryon correlators, it is convenient to first calculate gaugeinvariant three-quark propagators in which all summations over color indices and spatial sites have been done. A TABLE XIII. The numbers of operators of each type which project into each row of the G 1g , G 2g , and H g irreps for the , , N , and 0 baryons. The numbers for the G 1u , G 2u , and H u irreps are the same as for the G 1g , G 2g , and H g , respectively.
whereQ A aipja i p x; tjx 0 ; 0 denotes the propagator for a single smeared quark field of flavor A from source site x 0 at time t 0 to sink site x at time t. At the sink, a denotes color, is the Dirac spin index, i is the displacement direction, and p is the displacement length, and similarly at the source for a, , i, and p, respectively. Notice that the three-quark propagator is symmetric under interchange of all indices associated with the same flavor. As usual, translation invariance is invoked at the source so that summation over spatial sites is done only at the sink. These threequark propagators are computed for all possible values of the six Dirac spin indices.
Each baryon correlator is simply a linear superposition of elements of the three-quark propagators. These superposition coefficients are calculated as follows: first, the baryon operators at the source and sink are expressed in terms of the elemental operators; next, Wick's theorem is applied to express the correlator as a large sum of threequark propagator components; finally, symmetry operations are applied to minimize the number of source orientations, and the results are averaged over the rows of the representations. C++ code was written to perform these computations, and the resulting superposition coefficients are stored in computer files which are subsequently used as input to the Monte Carlo runs.
Wick's theorem is an important part of expressing the baryon correlators in terms of the three-quark propagators. To simplify the notation in the following, let the indices , , each represent a Dirac spin index and a displacement direction, and suppress the displacement lengths. Define c i c i 0 0 0 4 0 4 0 4 0 ; then the elements of the baryon correlation matrix in the channel are given in terms of three-quark propagator components (before source-minimizing rotations) by C ij t c i c j fG uuu jjj t G uuu jjj t G uuu jjj t G uuu jjj t G uuu jjj t G uuu jjj tg: The N correlators are expressed in terms of components of three-quark propagators by C N ij t c i c j fG uud jjj G uud jjj ÿG uud jjj ÿG uud jjj ÿG uud jjj ÿG uud jjj G uud jjj G uud jjj g; and for the and 0 channels, one finds C ij t c i c j fG uus jjj t G uus jjj tg; (53) C ij t c i c j fG uds jjj ÿG uds jjj ÿG uds jjj G uds jjj g:

VII. CONCLUSION
We plan to undertake a comprehensive study of the spectrum of QCD using Monte Carlo computations. Our first goal in this study is to calculate the masses of as many low-lying hadron resonances as possible. Successfully extracting these masses will depend crucially on using carefully designed spatially extended hadronic operators. In this first paper, the construction of three-quark , N, , , , baryon operators using group-theoretical projections was detailed. The operators were assembled out of gaugecovariantly displaced quark fields and transform according to the irreducible representations of the symmetry group of a spatial simple cubic lattice. Single-site, singly displaced, doubly displaced, and triply displaced three-quark operators were considered. The guiding principles in devising our operators were maximizing overlaps with the states of interest while minimizing the number of quark-propagator sources. Identifying the continuum-limit spins J of the states was addressed, and various issues related to computing the correlation matrix elements of the baryon operators were discussed.
Because of the complexity of these calculations and the importance of providing checks on our final results, we have been pursuing two different approaches to constructing the baryon operators. An alternative method of building the baryon operators based on Clebsch-Gordan techniques will be presented elsewhere [18]. The construction of meson and multihadron operators will be described in subsequent papers.
The Monte Carlo software to evaluate the correlation matrix elements of these baryons operators has been written and thoroughly tested using a large variety of checks, including comparison with known results in the case of a uniform constant background gauge field. This software uses the Chroma Software System for Lattice QCD [33] with QDP++ and QMP, developed under the Scientific Discovery through Advanced Computing initiative of the U.S. Department of Energy. Preliminary results have al-ready been presented in Refs. [34,35]. Although a very large number of baryon operators have been devised, it is not our intent to evaluate correlation matrices using all of these operators. Such calculations would not be feasible. The next important step in our study is to remove dynamically redundant and ineffective operators using lowstatistics Monte Carlo calculations, with the goal of finding some reasonably small subsets of operators adequate for extracting the low-lying masses of interest. Such calculations are currently in progress.