text:
# Replica Selection in the Globus Data Grid

## 1 Introduction

An increasing number of scientific applications ranging from high-energy physics to computational genomics require access to large amounts of data—currently terabytes and soon petabytes—with varied quality of service requirements. This diverse demand has contributed to the proliferation of storage system capabilities, thus making storage devices an integral part of the Grid environment and thereby constituting the Data Grid. The Globus  Data Grid architecture is an effort to standardize access to the multitude of storage systems spread across the Grid environment. It attempts to abstract these diverse elements of the data grid by providing a set of core services that can be used to construct a variety of higher-level services .

In this paper, we first describe the basic Globus Data Grid architecture, explaining briefly the various components. We then describe the process of identifying characteristics of interest about a storage replica resource, and we illustrate the mechanism with which these characteristics can be published. Next, we explain the architecture of replica selection, a higher-level service, built by using core services provided by the Globus data grid. We conclude by presenting the use of Condor’s classified advertisements (ClassAds) and matchmaking mechanisms  as an elegant matching and ranking tool in a storage context.

## 2 Architecture of Data Grid

The Globus data grid is organized into two layers, namely, core services and higher-level services that are built using the core services. The working hypothesis is that this hierarchical organization will make it possible to reuse services and code across a variety of applications and tools, so that, for example, both application-specific replica management solutions and sophisticated storage management systems such as the Storage Request Broker (SRB)  can share common low-level mechanisms.

### 2.1 Core Services

Core Data Grid services (Figure 1) seek to abstract the multitude of storage systems that exist in Grid environments, so that higher-level services—and applications—can access those storage systems uniformly. These core services include the following:

* Storage Systems and Data Access \- Storage system and data access services provide basic mechanisms for accessing and managing the data located in storage systems. These mechanisms provide abstractions for uniformly creating, deleting, accessing and modifying file instances across storage systems, regardless of their physical location. They target secondary and archival storage systems such as Unix file systems, High Performance Storage System (HPSS) , and Unitree and can also be used to access more sophisticated systems such as SRB.
* Metadata Access \- The metadata access service provides mechanisms to access and manage information about the data stored in storage systems. Various kinds of metadata are provided by this service:

1. Application Metadata describes the contents of the file, the circumstances under which the data was collected, and various details used by the application.
    2. Replica Metadata describes the mapping between file instances and particular replica locations.
    3. System Configuration Metadata describes the capabilities of storage systems, thus providing information on the fabric of the Data Grid.

The metadata service further provides a uniform method with which the various metadata can be published and accessed .

These services themselves build on basic security and information services provided by the Globus Toolkit .

### 2.2 Higher-Level Services

The core Data Grid services can be used to construct a variety of higher-level services (Figure 1). For example:

* Replica Management \- Replica management is the process of creating or deleting replicas at a storage site. Most often, these replicas are exact copies of the original files, created only to harness certain performance benefits. A replica manager typically maintains a replica catalog containing replica site addresses and the file instances.
* Replica Selection \- Replica selection is the process of choosing a replica from among those spread across the Grid, based on some characteristics specified by the application. One common selection criteria would be access speed .

## 3 Storage System Functionality and Replication

Intelligent replica selection requires information about the capabilities and performance characteristics of a storage system . We address this information discovery problem by leveraging machinery provided by the Globus Metacomputing Directory Service (MDS) , an information collection, publication, and access service for Grid resources.

The Globus MDS uses the Lightweight Directory Access Protocol (LDAP)  as its access protocol and LDAP object classes as its data representation, but it adopts innovative approaches to the problems of resource registration and discovery. Information about an individual resource or set of resources is collected and maintained by a Grid Resource Information Service (GRIS) daemon, which responds to LDAP requests with dynamically generated information and can be configured to register with one or more Grid Index Information Services (GIISs). Users will typically direct broad queries to GIIS to discover resources and then drill down with direct queries to GRIS to get up-to-date, detailed information about individual resources.

Information in LDAP is organized in a tree structure referred to as a Directory Information Tree (DIT), with nodes in a DIT tree corresponding to the LDAP structured data types called object classes. Figure 3 depicts the DIT hierarchy of object classes that we have defined to describe storage systems in a Data Grid environment. We describe these object classes in the following section.

### 3.1 Storage GRIS

Each storage resource in the Globus Data Grid incorporates a Grid Resource Information Server, configured to collect and publish system configuration metadata describing that storage system. This information typically includes attributes such as storage capacity and seek times, and descriptions of site-specific policies governing storage system usage. Figure 2 shows the object class definition that we have developed for this data. Attributes are labeled as either “MUST CONTAIN” or “MAY CONTAIN,” indicating whether or not they are mandatory.

This object class contains both dynamic and static attributes: attributes such as totalSpace, availableSpace, and mountPoint are dynamic, varying with various frequencies; attributes such as diskTransferRate, drdTime, dwrTime, and requirements are more static.

The requirements attribute is particularly interesting: it allows an administrator to specify the conditions under which the specified storage replica can be used, based, for example, on device utilization. For example, the requirements attribute could be a Boolean expression of maximum allowable storage and transfer bandwidth and can be specified using the ClassAd  mechanism. We discuss ClassAds in detail in a later section.

The data shown in Figure 2 is gathered in a variety of ways. We base our general approach on a generic LDAP-based GRIS server developed within the Globus project. The OpenLDAP server has a feature by which shell scripts (“shell-backends”) can be executed at the back-end in response to search queries; we use such scripts to gather dynamic attributes such as availableSpace, mountPoint, and totalSpace. Static attributes such as site usage policies and seek times can be specified by the administrator in a configuration file. Data collected in this manner is published in a suitable format (for example, LDIF ).

### 3.2 Data Access Service

We are also interested in the speed of a storage system, or, rather, in the time that the storage system can be expected to take to deliver a replica. One approach to determining this information is to construct a performance model of the relevant components (e.g., see ). We favor an alternative approach in which historical information concerning data transfer rates is used as a predictor of future transfer times. In brief, storage systems are configured to provide information on their own behavior and performance. Attributes such as maximum achievable read and write transfer bandwidths across networks can help an application choose one replica over another. Such data can be obtained by the storage replica by monitoring their own performance. This feature can be extended further, to obtain statistical information based on the performance data, such as average transfer bandwidths and their standard deviations, that can help predict the behavior of a particular replica. Further, we expect that there will be significant reuse of storage servers by clients, thereby justifying performance information on a per source basis, which provides a client useful information on end-to-end transfer performance between the server and the client. For example, a simple heuristic of combining past observed performance with current load of server might give a client a reasonably good choice of server.

Figures 4 and 5, along with Figure 3, depict the object classes used to record performance data associated with a replica location. The object class in Figure 4 specifies a summary of transfer bandwidth performance for all replica transfers, and the object class in Figure 5 specifies the performance details from the replica location to particular source sites. We gather this performance data by using instrumentation incorporated in the GridFTP server .

## 4 Classified Advertisements

Classified Advertisements (ClassAds)  are an elegant matching mechanism used in the Condor high-throughput computing system to map resource capabilities against job requirements. ClassAds allow resources and jobs to advertise their capabilities and requirements as attribute-value pairs that can be matched. Two ClassAds match if the logical expressions contained in the “requirements” attribute in both of them is satisfied. This requirements attribute itself can be represented in terms of other attributes. Further, ClassAds provide a mechanism for ranking the matches based on some attribute value, or value computed from multiple attributes. Up until now, ClassAds have been used extensively in Condor for job placements. In this section, we present the use of ClassAds in a storage context.

The various storage attributes described in the preceding section can be represented as attribute-value pairs in a ClassAd. For example, the following is a simple ClassAd describing the capabilities of a storage resource:

* hostname = ‘‘hugo.mcs.anl.gov’’;

volume = ‘‘/dev/sandbox’’;

availableSpace = 50G;

MaxRDBandwidth = 75K/Sec;

requirement = other.reqdSpace \< 10G

&& other.reqdRDBandwidth \< 75K/Sec;

This ClassAd describes a volume of a storage resource by specifying its attributes. The ClassAd also specifies a usage policy enforced by the resource, whereby only applications requiring storage less than 10 GB and transfer bandwidths less than 75 KB/sec are granted access. When two ClassAds are being matched, a MatchClassAd is created that contains both ClassAds. Each ClassAd can refer to the other ClassAd by using the “other” keyword.

A ClassAds representation of storage capabilities provides an efficient environment for matching, querying, and ranking requests. We will discuss the specifics of a request ClassAd and the matching process once we describe the replica selection mechanism.

## 5 Replica Selection

As noted above, a Data Grid provides a convenient environment for a community of researchers, interested in particular data sets, to maintain replicas of the data sets at their respective sites. Such an environment would provide both faster access and better performance characteristics. Replica selection is a high-level service provided by the Data Grid based on the core services. The replica selection process allows an application to choose a replica, from among those in a replica catalog, based on its performance and data access features.

An application that requires access to replicated data begins by querying an application specific metadata repository, specifying the characteristics of the desired data. The metadata repository maintains associations between representative characteristics and logical files, thus enabling the application to identify logical files based on application requirements. Once the logical file has been identified, the application uses the replica catalog to locate all replica locations containing physical file instances of this logical file, from which it can choose a suitable instance for retrieval.

The entity that identifies the suitable instance of a replicated file based on application requirements is referred to as a broker. In effect, the responsibility of the broker is to map application requirements against storage resource capabilities. In the following section, we discuss the storage broker in detail.

### 5.1 The Storage Broker

In this section, we first delve into the details of the storage broker architecture. We then discuss the overall architecture, a decentralized selection process, and a few issues involved in its realization.

#### 5.1.1 Decentralized Selection Process

Traditionally, resource brokers have adopted a centralized approach to resource management, wherein a single node is responsible for decision making. An example of such an environment is the Condor  high-throughput computing platform, wherein a central manager is responsible for matching resources against jobs. Obvious disadvantages to this approach are scalability and a single point of failure. Of course, Condor has an efficient recovery mechanism to address failure and has been proven to scale to thousands of resources and users.

But there is a more fundamental problem with this centralized approach when applied to Grids. In these highly distributed environments, there are numerous user communities and shared resources, each with distinct security requirements. No single resource broker is likely to be trusted by all of these communities and resources with the necessary information to make decisions. At the extreme, each user may need his or her own resource broker, because only that user has the authorization to gather all of the information necessary to make brokering decisions.

For this reason, we have designed a decentralized storage brokering strategy wherein every client that requires access to a replica performs the selection process rather than a central manager performing matches against clients and replicas.

#### 5.1.2 Architecture

Figure 6 presents a snapshot of the Grid environment where storage resources are scattered and each client requiring access to a replica initiates a decentralized replica selection mechanism. As can be seen, there is no central point of control, and decision making is delegated to each and every client.

An application requiring access to a file presents its requirements as a ClassAd to the broker. These requirements might be as simple as a Boolean expression indicating the storage and transfer bandwidth required or may indicate more complex constraints on storage system type, state, and policy. The broker then performs the following sequence of actions:

* Search Phase

1. The broker attempts to find a suitable replica matching the application’s requirements. To achieve this, it queries the replica catalog, which contains addresses of all replicas for each logical file.
    2. The next step is to query each replica location. As we have seen before, this involves using LDAP searches to query GRIS servers associated with storage systems.
    3. In response to the LDAP search, each storage system returns its capabilities and usage policies in the LDAP Information Format.
    4. The broker collects the capabilities of all replica resources and proceeds to the matching phase.
* Match Phase

1. The replica capability data is converted into ClassAd format in preparation for invoking the Condor matchmaking mechanisms. We thus obtain a list of classified advertisements, representing the various replica sites.
    2. The broker then performs a match of the application’s requirement ClassAd against the list of replica capability ClassAds, obtaining a set of replica locations that satisfy the criterion.
    3. The ClassAd ranking feature can be used to prioritize successful matches based on some attribute, specified by the application.
* Access Phase

1. Once a suitable replica has been identified, the file is accessed using a high-speed file transfer protocol, for example the GridFTP tools provided within the Globus Toolkit.

### 5.2 An Example Application Request

In Section 4, we described a storage ClassAd that indicates its willingness to accept application requests that require space less than 10 GB and require a transfer bandwidth less than 75 KB/sec. In this section, we look at a user request ClassAd and how it is matched against the storage ClassAd.

An application might advertise its request to the broker as follows:

* hostname = ‘‘comet.xyz.com’’;

reqdSpace = 5G;

reqdRDBandwidth = 50K/Sec;

rank = other.availableSpace;

requirement = other.availableSpace \> 5G && other.MaxRDBandwidth \> 50K/Sec;

The application indicates its preference for a storage resource that has more than 5 GB available space and a maximum transfer bandwidth greater that 50 KB/sec. The broker, upon receiving such a request, attempts to contact all replica locations identified by the replica catalog, using LDAP search queries to request the attributes of interest: in this case, availableSpace and MaxRDBandwidth. The broker thus uses the application ClassAd to build specialized LDAP search queries.

The results obtained via these LDAP queries are converted from LDAP Interchange Format (LDIF) into ClassAds and standard Condor mechanisms are invoked to match the application ClassAd with each ClassAd in the list of storage system ClassAds. (Note that the storage system ClassAds comprise their site specific usage policy, as mentioned in the sample storage ClassAd in Section 4.) Any matched ClassAd can be further ranked by querying their rank attribute. In our case, we rank the replica servers based on their available space, thus obtaining the “best” match .

## 6 Results

We summarize our initial results based on our experience in building this prototype.

We have succeeded in building a replica selection broker and have demonstrated that this high-level Data Grid service can be built using the services provided by the Globus Data Grid Toolkit, such as: Storage GRIS, GridFTP and Replica Catalog.

We have further demonstrated the process of: Identifying characteristics of interest on a storage resource, based on Storage GRIS and GridFTP protocols; Classifying them as belonging to appropriate object classes according to LDAP; Gathering these attributes, via such mechanisms as shell backend scripts, tuning FTP servers and configuration files; Publishing the features in a suitable, efficiently queriable format using LDAP protocol.

We have explored the use of Condor ClassAds as a mechanism for expressing storage resource capabilities in the Grid environment. Although LDAP provides an equivalent method of publishing characteristics in attribute-value pairs, ClassAds provides a richer matching and ranking environment. Further, in our prototype implementation, we have demonstrated that the process of converting data, represented in LDAP format, into ClassAds is not cumbersome and is worth the effort. We have, in fact, developed primitive libraries to achieve the conversion of this attribute set.

## 7 Conclusions and Future Directions

Our current prototype implementation is primarily a proof of concept of the Globus Data Grid services. As the next logical step in its development, we need to demonstrate its applicability in a real application, thereby corroborating its usefulness. Further improvements can be made with regards to the information published in the storage GRIS. Finally, the statistical information published by the storage resource can be fed to an information service, such as the Network Weather Service , to perform predictive analysis of the behavior of storage resources.

## Acknowledgments

This work was supported in part by the Mathematical, Information, and Computational Sciences Division subprogram of the Office of Advanced Scientific Computing Research, U.S. Department of Energy, under Contract W-31-109-Eng-38; by the National Science Foundation; and by the NASA Information Power Grid program.


format: text

text:
# The QCD spectrum with three quark flavors

## I Introduction

Computation of the properties of hadrons beginning from the QCD Lagrangian is a major goal of lattice gauge theory, and steady progress has been made. The computational burden of including dynamical quarks is a major obstacle in the use of lattice QCD to compute hadronic properties. As a result, many quantities are much better determined in the quenched approximation than in full QCD. One way of studying the effects of dynamical quarks is to calculate in quenched and full QCD, using the same valence quark action in both cases, and matching the lattice spacings and physical sizes of the lattices, so that any differences that are found can convincingly be ascribed to the dynamical quarks. Here we present a calculation of the hadron spectrum at a lattice spacing of about $`a=0.13`$ fm, using quenched and full QCD lattices at the same lattice spacing. The lattice spacing was tuned by making short runs on smaller lattices, adjusting the parameters to match the lattice spacing of an initial quenched run at $`10/g^2=8.0`$.

We use an improved Kogut-Susskind quark action which removes tree level order $`a^2`$ lattice artifacts. The gauge action is a one-loop Symanzik improved action. This action has been shown to reduce flavor symmetry breaking and to improve rotational symmetry of the hadron spectrum, and to give improved scaling of hadron masses as a function of lattice spacing.

Another important improvement of these calculations over previous generations is that we use three flavors of dynamical quarks. For quark masses larger than the strange quark mass we use three degenerate flavors, and for light quark masses less than $`m_s`$ we use two light flavors, keeping the third quark mass at about the strange quark mass. We have also done one two-flavor simulation on a matched lattice to check for effects of the dynamical strange quark. For the runs with $`2+1`$ dynamical flavors, we computed hadron spectra using valence quark masses equal to the sea quark masses. In the quenched run we computed hadron masses with valence quarks with the same masses, and nondegenerate propagators using a strange quark mass of $`am_s=0.05`$. Finally, in the two dynamical flavor run we computed hadron propagators using light quark masses equal to the dynamical mass, $`am_{u,d}=0.02`$ and a non-dynamical strange quark with $`am_s=0.05`$.

Two issues that we do not completely address are extrapolation to zero lattice spacing and extrapolation to the physical light quark mass. Using matched quenched and full QCD lattices allows us to draw conclusions about the effects of dynamical quarks without explicit continuum extrapolation. Although it is in principle possible that the discretization errors in quenched and full QCD are very different, we expect that the differences in lattice spacing dependence are in fact comparable to the differences in physical quantities themselves. Since corrections to scaling are in any case quite small with our current improved action (see also below), we confidently expect conclusions drawn at fixed lattice spacing to survive in the continuum limit. We are beginning a series of simulations at a smaller lattice spacing which will eventually allow us to make the continuum extrapolation. A few preliminary quenched points from these finer lattice runs are included here to give an idea of the size of these effects. A complete chiral extrapolation will be more difficult. In this work, we attempt an explicit chiral extrapolation only for the shape of the static quark potential, and show other quantities as functions of the quark mass. Some quantities such as “$`J`$ are only minimally sensitive to chiral extrapolation, and such quantites provide immediately useful tests of dynamical quark effects.

In addition to presenting hadron spectra, we update our computation of the static quark potential. For this quantity we have sufficiently accurate data that we can hazard an extrapolation to zero quark mass to produce numbers that can be compared with phenomenological potential models. Because the static potential is determined very accurately, it clearly shows differences between quenched and full QCD. In fact, one can even see the differences between two and three dynamical flavors, and a “kink” in the plots where we change from three degenerate flavors to two light and one heavy indicates a noticeable difference between two light plus one heavy and three light dynamical flavors.

In the meson sector we find differences between full and quenched QCD. A nice way of exposing these differences is the ratio $`J`$ proposed by Lacock and Michael. We find that this quantity is in better agreement with experiment in full QCD than in the quenched approximation, as predicted in . This is consistent with results of the CP-PACS and JLQCD collaborations, who also conclude that inclusion of two flavors of dynamical quarks improves agreement of the lattice spectrum with the real world.

In the isovector $`J^{PC}=0^{++}`$ ($`a_0`$) channel we find a large difference between quenched and three-flavor results. We interpret the three flavor results as an avoided level crossing between the $`0^{++}`$ meson and a two pseudoscalar state. In other words, we see the $`a_0`$ decay to two mesons.

We include tables of the mass fits we have chosen, so the reader can compute his or her own favorite mass ratios.

## II Simulation parameters

For our two and three flavor simulations we used the standard hybrid-molecular dynamics “R algorithm”,  with one pseudofermion field for runs with degenerate quarks, and two pseudofermion fields for runs with different strange and up and down quark masses. In all cases, we used trajectories with unit length in the simulation time. Basic parameters of these runs are summarized in Table I.

Two sources of systematic error in this method are the accuracy of the approximate sparse matrix solution required at each time step in the integration of the molecular dynamics equations, and the effect of the nonzero time step used in the integration. We have investigated these effects in the 2+1 flavor simulations at $`10/g^2=6.8`$ with two flavors of quarks at mass $`am_{u,d}=0.02`$ and one flavor with $`am_s=0.05`$ on a $`16^3\times 48`$ lattice. (These masses are approximately 0.4 times the physical strange quark mass and the physical strange quark mass, respectively). Figures 1 and 2 show the plaquette and $`\overline{\psi }\psi `$ as a function of the conjugate gradient residual used in the updating. Finally, to see how this effect propagates into the hadron masses, we show the Goldstone pion mass in these same runs in Fig. 3. Since the effect of this error is poorly understood, we adopted a conservative choice of $`1\times 10^4`$ or $`5\times 10^5`$ in most of our runs, and used $`2\times 10^5`$ for the heaviest quark mass ($`am_q=0.40`$) where convergence was very fast.

The effect of integration step size is better understood. We verified the expected quadratic dependence of the error on the step size in our three flavor code. The first panel in Fig. 4 shows the result of one such test, where we ran on a quite coarse lattice using the conventional action, comparing the old code with three degenerate flavors to the two-plus-one flavor code with $`m_{u,d}=m_s=0.02/a`$. (The “old code” uses a single pseudofermion vector, with a weight of $`3/4`$ in the fermion force, while the “two-plus-one flavor” code uses two pseudofermion vectors, one with a weight of $`2/4`$ and the other with weight $`1/4`$. In this test, both pseudofermion vectors had the same mass. In each case, the multiplication of the Gaussian random vector by $`M^{}`$ to produce the pseudofermion vector was done at the appropriate point in the time step to make the error in a single time step order $`ϵ^3`$, making the accumulated error over a trajectory order $`ϵ^2`$..) It can be seen that both sets of points extrapolate to the same limit, although the size of the effect is quite different. The right hand panel of Fig. 4 shows a similar plot from one of our pre-production tuning runs with the improved gauge and quark action, at $`10/g^2=6.80`$ and $`am_q=0.05`$, on a $`16^3\times 48`$ lattice. This is close to the value of $`10/g^2`$ used in the production run at this quark mass. Note that the finite step size corrections are quite sensitive to the action being used – even the sign of the effect differs between the two actions. Based on these tests and previous experience, we used a step size of no more than two-thirds of the light quark mass or 0.03, whichever was smaller, in our production runs. In the production run at $`am_q=0.05`$ we used a step size of $`ϵ=0.02`$. From the slope in Fig. 4 we can infer that this caused a fractional error in $`\overline{\psi }\psi `$ of about 0.004, and a similar analysis for the plaquette suggests a fractional error of about 0.0006. We also looked at pion masses and some of the Wilson loops involved in the computation of the static quark potential, and we were unable to resolve statistically significant effects on these quantities in our tuning runs.

In our production runs with dynamical quarks we measured the potential and the spectrum at intervals of six simulation time units, and archived these lattices. The autocorrelation of the plaquette at six time units, or successive measurements, was generally about 0.1. We investigated the effect of autocorrelations on the potential and spectrum by blocking together different numbers of measurements before doing the fitting. For the potential measurements we chose to block five measurments (30 time units). For the hadron spectrum most of the particles showed no systematic effects of blocking measurements together. The exception was the pseudoscalar mesons, where blocking gave significantly better confidence levels and larger error estimates. For the pseudoscalars we chose to block together four measurements, or 24 time units.

The physical size of our lattices is $`La=20a2.6`$ fm, which is similar to or larger than other recent full QCD simulations. Basic parameters of several of these calculations are summarized in Ref. .

## III Length scales from the static potential

The static quark potential is often used to define the length scale in lattice simulations. Advantages of using the potential include ease and accuracy of its computation, and its lack of dependence on the valence quark mass. In comparing quenched and full simulations, subtleties arise because the potential does depend on the masses of the sea quarks. In Ref.  we demonstrated the effects of sea quarks on the potential using this improved action. Because these effects are important in our analysis of the hadron spectrum, we update and extend these results here. Our methods for computing the potential and our reasons for using $`r_1`$, a variant of the conventional $`r_0`$ , are described in . $`r_0`$ is conventionally defined by $`r_0^2F(r_0)=1.65`$, and $`r_1`$ by $`r_1^2F(r_1)=1.00`$.

The fitting form used here is slightly more complicated than the form used in Ref. , with an extra term to take into account lattice artifacts at the shortest distances. Following a procedure used in Ref. ,

$$V(\stackrel{}{r})=C+\sigma r\alpha /r+\lambda \left(V_{free}(\stackrel{}{r})1/r\right).$$
(1)
The last term, used for $`r<2.5`$, approximately compensates for remaining lattice artifacts. Here $`V_{free}(\stackrel{}{r})`$ is the potential calculated in free field theory, using the improved gauge action. Adding this term to the fits significantly improves the goodness of fit and makes the fit parameters less sensitive to the choice of distance range. For the $`a0.13`$ fm runs we typically find $`\lambda 0.30.4`$.

In Figs. 5 and 6 we show the dimensionless quantities $`r_0\sqrt{\sigma }`$ and $`r_1\sqrt{\sigma }`$ respectively as functions of the quark mass, represented by $`(m_\pi /m_\rho )^2`$. This places the quenched approximation at $`(m_\pi /m_\rho )^2=1`$, and the chiral limit at the left side of the graph. In these plots the octagons are runs with three degenerate sea quarks, except for the rightmost point which is the quenched limit. Squares are runs with $`am_s=0.05`$, its approximate physical value, and $`am_{u,d}<0.05`$. The isolated diamond is our two flavor run. Finally, the cross at $`(m_\pi /m_\rho )^2=1`$ is the finer lattice quenched run. From the two quenched points we see that remaining lattice artifacts are small compared with the effects of the sea quarks. In particular, the central values for $`r_0\sqrt{\sigma }`$ and $`r_1\sqrt{\sigma }`$ changed by less than 1% when the lattice spacing was reduced from 0.13 fm to 0.09 fm, a change of about 35%. The kink in the plots at $`(m_\pi /m_\rho )^20.46`$ ($`am_{u,d}=0.05`$) shows the transition between three degenerate flavors and “$`2+1`$” flavors. We can clearly see the distinction between two and three flavors, as well as the effect of using two light and one heavy flavor rather than three degenerate flavors (the “kink’ at $`(m_\pi /m_\rho )^20.46`$).

If we extrapolate $`r_0\sqrt{\sigma }`$ to the physical quark mass, as shown in Fig. 5, we find $`r_0\sqrt{\sigma }=1.114(4)`$ (statistical errors only), a number which can be compared with phenomenological potential models. The two quenched points give some idea of the possible systematic error. Since the squared lattice spacing in the finer lattice is about one half that of the coarser, we might expect a shift of about twice the separation of these points in the continuum limit. Since the error bars on these points overlap, we don’t know this systematic error well enough to justify such an extrapolation at this point. We expect that the effect of this systematic error is mostly an overall shift of the graphs, but the next generation of simulations should clarify this.

While using $`r_1`$ to define the length scale has the advantage that it can be done more accurately, $`r_0`$ has the advantage that it has been related to phenomenological potential models, which consistently place it around 0.5 fm. Therefore, to estimate $`r_1`$ in physical units, we plot $`r_0/r_1`$ in Fig. 7. Extrapolating linearly in $`(m_\pi /m_\rho )^2`$ to the physical value gives $`r_0/r_1=1.449(5)`$ (statistical error only), or with $`r_0=0.5`$ fm, $`r_1=0.35`$ fm.

## IV Procedures for determining masses

All of our hadron propagators used wall sources and local sink operators. Several different wall sources were used. For the “pointlike” hadrons, for which all the quarks can be on a single corner of the hypercube, a “corner wall” source gave the best results. This source is simply a 1 on each $`(0,0,0)`$ corner of the $`2^3`$ cubes on the chosen time slice. However, to isolate the decuplet baryons a non-pointlike source is essential, and we used “even” and “odd” wall sources, where $`1`$ or $`(1)^{x+y+z}`$ is placed on each site, respectively. This set of sources was developed in Ref. . For the nonlocal pseudoscalar and vector mesons we used two wall sources made from empirically determined linear combinations of the nonlocal pion operators. Finally, for the nonzero momentum mesons we used a quark source with $`1`$ on each site, and an antiquark source with $`e^{i\stackrel{}{k}\stackrel{}{x}}`$ on each site. All of the configurations were gauge fixed to the Coulomb gauge before computing the propagators.

In most cases we computed propagators from four source times evenly spread through the lattice (only one source slice was “turned on” at a time). For the corner source we used eight source time slices for the light quark particles with $`am_q0.04`$, and we used eight source time slices for the “even” and “odd” source baryon propagators. (The $`\mathrm{\Delta }`$ propagator is very noisy, and propagators computed from source times separated by $`8\times 0.13`$ fm were basically independent.)

For Kogut-Susskind quarks the meson propagators have the generic form

$$(t)=\underset{i}{}A_i\left(e^{m_it}+e^{m_i(N_tt)}\right)+\underset{i}{}A_i^{}(1)^t\left(e^{m_i^{}t}+e^{m_i^{}(N_tt)}\right).$$
(2)
Here the oscillating terms correspond to particles with opposite parity from the ordinary exponential terms. Baryon propagators are similar, but have antiperiodic boundary conditions and the “backwards” terms include an extra factor of $`(1)^t`$. In most cases only one mass with each parity is included in the fits, but for half of the pseudoscalar meson operators the opposite parity terms are not present, and for the P-wave mesons we found it necessary to keep two simple exponentials. The quantum numbers for the various operators are tabulated in Ref. .

Hadron masses were determined from fits to propagators, using the full covariance matrix to estimate errors. The maximum time distance used in the fits was chosen to include points with fractional error less than 0.3. Because of the oscillating components in the staggered quark propagators, it sometimes happens that the fractional error exceeds the threshold at one distance but is smaller at larger distances, so the complete criterion for the maximum distance included is the largest distance such that the fractional error on each of the next two points exceeds 0.3. Since the points at largest distance contribute little information, the exact large distance cutoff is not critical. To choose the minimum distance included, we first went through the fits and chose a minimum distance for each hadron in each dataset, choosing a distance where the confidence level was reasonable and where the mass appeared to reach a plateau. As expected, some propagators had larger(smaller) fluctuations than other similar propagators, resulting in the choice of a larger(smaller) minimum distance. For the particles of greatest interest, to reduce this effect we then “smoothed” these minimum distances, requiring that the minimum distances be smooth functions of quark mass and be the same for the quenched and dynamical runs. The resulting minimum distances are strongly dependent on the quark masses, with smaller minimum distances for smaller masses. Most of this is due to the larger statistical errors at smaller quark mass, which result in the excited state contributions disappearing into the noise at shorter distance. However, it is in part physical, since splittings between the ground and excited states are larger for smaller quark masses. Table II shows the minimum distances that we chose, and the number of particles with each parity.

## V Results

### A Pseudoscalar mesons

We calculated propagators for all eight flavor combinations of the staggered quark pseudoscalar mesons. These masses obey the “partial flavor symmetry restoration” predicted by Lee and Sharpe to very good accuracy. Specifically, Ref.  predicts that the leading order flavor symmetry breaking effects, which are order $`a^2`$, leave degeneracies between pairs of pseudoscalar mesons for which $`\gamma _0`$ is replaced by $`\gamma _i`$ in the flavor structure. For example, the local non-Goldstone pion, $`\gamma _0\gamma _5\gamma _0\gamma _5`$ in the “spin $``$ flavor” notation, is degenerate with the distance one pion, $`\gamma _5\gamma _i\gamma _5`$, to this order. Moreover, all of the squared pion masses should depend linearly on the quark mass with the same slope to lowest order. Figure 8 shows this behavior for the quenched pion masses for $`am_{u,d}0.05`$. The results for the full QCD runs are similar. Since the local $`\gamma _5\gamma _5`$ pseudoscalar has the correct chiral behavior (and the best scaling behavior), we will use this pseudoscalar in the rest of the analysis unless we specifically indicate otherwise. Results for the full QCD runs are similar, but the flavor symmetry breaking is somewhat larger.

In Fig. 8 the relation between the squared pseudoscalar mass and the quark mass is clearly nearly linear. The deviations from linearity and the effect of the dynamical quarks can be exposed by plotting the squared pseudoscalar masses divided by the quark mass, in Fig. 9. This is essentially $`\psi \overline{\psi }r_1/f_\pi ^2`$ with an (unknown) renormalization factor. This plot contains pseudoscalar mesons with both light and strange valence quarks (pions, kaons and “unmixed $`s\overline{s}`$’s”). There is clearly a systematic difference between quenched and full QCD. This difference increases with decreasing quark mass, and the two-flavor point falls in between the quenched and three-flavor points. The bursts among the quenched points are from the $`10/g^2=8.4`$, $`a0.09`$ fm run, showing gratifying agreement with the $`a0.13`$ fm points. Unfortunately, a coarser three flavor lattice, $`a0.2`$ fm, shows a large effect, so we would not want to use much coarser lattices in studying this effect. We do note that we expect scaling violations to be similar for the quenched and dynamical theories, so it is an advantage to have runs with matched lattice spacings. The deviations from linearity of $`m_\pi ^2`$ are similar in quenched and full QCD. The upturn for larger quark masses signals the beginning of the transition to the heavy quark regime, where $`m_\pi ^2m_q^2`$. We do not fully understand the shape of this plot for small quark mass. There are several ways to interpret the difference between quenched and dynamical results. One could say that $`\overline{\psi }\psi /f_\pi ^2`$ is too small in the quenched approximation, or one could say that the quark mass at which a desired $`m_\pi /m_\rho `$ is reached is larger in the quenched approximation than in full QCD. This second interpretation is consistent with CP-PACS results on the quark masses, in which they find that the quark mass needed to reach a given value of $`m_\pi /m_\rho `$ is smaller in two flavor QCD than in quenched QCD. (Indeed, one could even use this quantity as a length scale, and conclude that $`r_1`$ is different in quenched and full QCD.)

The largest part of the error bars in Fig. 9 come from the uncertainty in $`r_1`$. However, this uncertainty is common to all of the points coming from the same set of lattices. In particular, all the $`a0.13`$ fm quenched points are correlated in this respect, as are the three two-flavor points. If we are interested in the dependence of the pseudoscalar mass on the quark mass on a fixed lattice, we may want to consider only the error from the determination of the meson mass in units of $`a`$. The left hand panel in Fig. 10 shows $`m_\pi ^2r_1^2/(r_1(m_1+m_2))`$ for the quenched calculation, including only the error from $`am_{PS}`$ and showing only the reasonably light mass points. In this panel the octagons are “pions”, with $`m_1=m_2=am_{u,d}`$, and the bursts are “kaons”, with $`am_2`$ fixed at 0.05, which is approximately the physical value of the strange quark mass. We see that this quantity is dependent only on the sum of the quark masses to very good accuracy. The center panel contains the same plot for the three-flavor runs, where now $`r_1`$ is determined independently in each run. In this panel the octagons are runs with $`am_{u,d}=am_s0.05`$, the squares the run with $`am_{u,d}=0.04`$, the diamonds from $`am_{u,d}=0.03`$, the crosses from $`am_{u,d}=0.02`$ and the bursts from $`am_{u,d}=0.01`$. The three symbols for each of the runs with $`am_{u,d}<0.05`$ correspond to the “pion”, with both valence quarks light, the “kaon”, with one light and one strange valence quark, and an “unmixed $`s\overline{s}`$”, with two valence quarks of mass approximately equal to that of the strange quark, but no $`q\overline{q}`$ annihilation. This graph is far from smooth, but most of the scatter comes from the fact that each set of dynamical quark masses has an independent uncertainty in $`r_1`$. Note that each dynamical run (for example, the three bursts) shows qualitatively the same behavior as the quenched case, with the light-light pseudoscalar tending to a larger value. An interesting question is how the sea quark mass affects the pseudoscalar mass. This can be investigated by looking at the points at $`a(m_1+m_2)=0.1`$, which are the “unmixed $`s\overline{s}`$” points with both valence quark masses equal to $`0.05/a`$. In the right hand panel of Fig. 10 we plot these points as a function of the light quark mass. There is a noticeable effect, with smaller light quark mass producing larger $`s\overline{s}`$ mass. The direction of this effect is consistent with the smaller pseudoscalar masses in the two flavor and quenched calculations seen in Fig. 9.

The selected pseudoscalar meson mass fits in units of the lattice spacing are tabulated in table III. In addition to the “pions”, the table also contains fits with one quark at about the strange quark mass and one lighter quark, or “kaons”. For the two and three flavor runs we also tabulate “unmixed $`s\overline{s}`$ mesons”, with two valence quarks with $`am_v=0.05`$.

We have attempted to fit the results in Fig. 9 to the forms predicted by chiral perturbation theory. In the quenched case, the behavior of $`m_\pi ^2`$ as a function of quark mass is derived in Refs. . We use eqn. (9) in , with the parameter $`\alpha `$, which is believed to be small, set equal to zero, and the analytic correction term added to the chiral log. For pions, then one has

$$\frac{m_\pi ^2}{m}=C(1\delta \mathrm{log}\left(\frac{Cm}{\mathrm{\Lambda }^2}\right)+Km),$$
(3)
where $`m_1=m_2m`$, $`C`$ and $`K`$ are constants, and the chiral scale $`\mathrm{\Lambda }`$ may be taken as the $`\eta `$ mass. As shown in Fig. 9, the fit to the $`a=.13`$ fm quenched data is good. It gives $`\delta =0.061(3)`$ (statistical error only), which is on the low side but compatible with the range reported by CP-PACS  and is in excellent agreement with the result of Bardeen et al. , $`\delta =0.065(13)`$.

Unfortunately, our attempts to fit the 3-flavor pion data in Fig. 9 to the corresponding full QCD chiral form  have been unsuccessful to date. We can fit the 5 lowest-mass pions with reasonable confidence level, but the coefficients of the analytic terms are unreasonably large, and the fit misses the next lightest pion by a wide margin. If we try to fit the 6 lowest-mass pions, the fit has terrible confidence level. Finally, good fits can be obtained by introducing, as an additional free parameter, an overall coefficient in front of the chiral logs. However, the value of that coefficient in the fit is much smaller than its predicted  value. We are continuing to study this puzzling situation. Our current running at smaller lattice spacing may provide additional insight here.

### B Vector mesons

We calculated propagators for the two local vector mesons, $`\gamma _i\gamma _i`$ (VT) and $`\gamma _0\gamma _i\gamma _0\gamma _i`$ (PV), and two distance one vector mesons, $`\gamma _i\mathrm{𝟏}`$ and $`\gamma _0\gamma _i\gamma _0`$. Any flavor symmetry breaking among these mesons is smaller than the statistical errors, so we simply quote results for the local $`\gamma _i\gamma _i`$, or “VT” mesons. Table IV contains these masses in units of the lattice spacing.

In Figure 11 we plot the vector meson mass in units of $`r_1`$ versus the squared pion/rho mass ratio. In this plot there is a clear difference between the quenched and dynamical masses, with the full QCD vector mesons lying lower. However, the size of this effect depends on the length standard chosen, as illustrated in Fig. 12, where the same quantity is plotted in units of the string tension. Of course, the difference between these two plots simply arises from the difference in $`r_1\sqrt{\sigma }`$ plotted in Fig. 6.

In our calculations we used a mass of $`am_s=0.05`$ for the strange quark mass. With the meson spectrum in hand, we can now go back and ask whether this choice was exactly correct. There are several quantities that we could choose to define the strange quark mass, and, especially in the quenched calculation, they will in general give different results. Perhaps the most accessible quantity is the ratio of the $`s\overline{s}`$ pseudoscalar to vector meson mass. (Our “$`s\overline{s}`$ pseudoscalar” does not include $`q\overline{q}`$ annihilation diagrams, so it is not the $`\eta `$ or $`\eta ^{}`$.) We therefore tune this quantity using $`m_{s\overline{s}}=\sqrt{2m_K^2m_\pi ^2}=686`$ MeV and $`m_\varphi =1020`$ MeV, or $`m_{PS}/m_V=0.673`$.

In the quenched spectrum with $`am_q=0.05`$ we have $`m_{PS}/m_V=0.643`$, indicating that $`am_q=0.05`$ is smaller than the desired strange quark mass. To compute the mass at which $`m_{PS}/m_V=0.673`$, we can do a linear fit to the vector meson mass as a function of quark mass, and combine that with a squared pseudoscalar mass proportional to the quark mass to conclude that the quenched strange quark mass defined by $`m_{PS}/m_V`$ at this lattice spacing is $`am_s=0.058`$. In contrast, for the dynamical runs $`am_s=0.05`$ is a fairly good estimate of the strange quark mass, with $`m_{PS}/m_V`$ taking values between 0.670(2) and 0.687(1).

Differences between the quenched meson spectrum and the real world have been observed by the UKQCD collaboration, and improvements of the spectrum when dynamical quarks are included have been reported by the CPPACS and JLQCD collaborations. In particular, the UKQCD collaboration studied the quantity

$$J=m_K^{}\frac{m_V}{m_{PS}^2},$$
(4)
where $`m_V`$ and $`m_{PS}`$ are the vector and pseudoscalar meson masses. This quantity has the advantage of being relatively insensitive to the quark masses, so that accurate tuning of the strange quark mass or extrapolation of the masses to the chiral limit is not essential. Of course, to compare to experiment the derivative in this expression must be replaced by a ratio of mass differences, and we choose

$$J=m_K^{}\frac{\left(m_\varphi m_\rho \right)}{2\left(m_K^2m_\pi ^2\right)}.$$
(5)
Here $`m_\rho `$ is the mass of the vector meson including two light quarks, etc. We choose the $`m_\varphi m_\rho `$ mass difference because the statistical error in $`J`$ is dominated by the error in the vector meson masses, and the larger difference in $`m_\varphi m_\rho `$ relative to, say, $`m_K^{}m_\rho `$ leads to smaller statistical errors on $`J`$. Because all of the masses in Eq. 5 are correlated, we use a jackknife analysis to compute the error on $`J`$. Figure 13 shows the results for $`J`$ in quenched and three flavor QCD. Following UKQCD, we plot this versus $`m_K^{}/m_K`$, for which the real world value is 1.8. The burst is the real world value of this definition of $`J`$ (0.49), and the cross is the value of $`J`$ found in the UKQCD quenched simulations. We see a clear effect of the sea quarks on this quantity. Indeed, any reasonable extrapolation of our data in $`m_K^{}/m_K`$ would pass near the real world point. Figure 13 also contains one point with two dynamical flavors. This point falls near the three flavor points, indicating that the dynamical strange quark is less important than the two light quarks. Although our quenched results are somewhat higher than the UKQCD value, they are significantly below the experimental value. The fact that the quenched points in this plot are to the right of the full QCD points is largely due to the fact that the mass of $`am_q=0.05`$ used for the quenched strange quark was too small, as discussed above. If we use the observed slopes of the quenched vector meson masses and squared pseudoscalar masses as functions of quark mass to adjust these points to a quenched strange quark mass of $`am_s=0.058`$, the main effect is to shift the quenched points to the left. In particular, the rightmost quenched point moves to $`m_K^{}/m_K=1.78`$, but moves up by only 0.004.

### C P-wave mesons

Mesons with quantum numbers $`J^{PC}=0^{++}`$, $`1^{++}`$ and $`1^+`$, which are P-wave mesons in the nonrelativistic limit, are found as the oscillating “parity partners” of the $`\pi `$ and $`\rho `$ propagators with Kogut-Susskind quarks. Fitting these particles is difficult because of their larger mass. In order to get good fits, we need to allow two particles in the non-oscillating components, since we typically find that the excited “$`\pi `$” or “$`\rho `$” state is comparable in mass to the lowest oscillating (P-wave) state. Thus we find large errors, and must use small minimum distances to get good fits. In addition, plateaus in the effective mass are short. Figure 14 illustrates an example of a difficult but crucial case. This figure shows fits to the scalar ($`a_0`$) mass as a function of the minimum distance included in the fit. With these caveats, selected fits for the P-wave mesons are shown in Fig. 15 and Tables V, VI and VII. This figure and these tables contain several interesting features. For the $`a_1`$ and $`b_1`$ the full QCD runs give consistently better fits than the quenched run, and smaller masses for the lighter quarks. The $`b_1`$ is consistently slightly heavier than the $`a_1`$, although with the difficulties in extracting these particles we would not want to make too much of this. The diamond at the left in Fig. 15 is the experimental value for the $`a_1`$ and $`b_1`$ masses.

The scalar channel, “$`a_0`$”, is clearly very different in the quenched and full QCD runs. For large quark masses there is no visible difference, but as the quark mass is decreased the full QCD $`0^{++}`$ mass drops below all the other masses. For all but the lowest quark mass, the quenched $`0^{++}`$ is close to the other P-wave meson masses. We ascribe the behavior of the full QCD mass to the decay of the $`a_0`$ into $`\pi +\eta `$. (Bose symmetry plus isospin forbids decay into two pions.) Figure 16 illustrates this interpretation. In the figure we plot the quenched and full $`0^{++}`$ masses versus quark mass. The straight line in the graph is a fit to the quenched mass for the heavier quarks, and represents the mass of a $`q\overline{q}`$ state. The curved line with the kink at $`am_q=0.05`$ represents the mass of $`\pi +\eta `$. For $`am_q0.05`$ we used three degenerate quark flavors, so the $`\eta `$ and $`\pi `$ are degenerate and this line is simply twice the pion mass. For $`am_q<0.05`$ we don’t have direct information on the $`\eta `$ mass, so we use the Gell-man–Okubo formula written in terms of an “unmixed $`s\overline{s}`$” mass (just our pseudoscalar mass at $`am_q=0.05`$).

$$m_\eta ^2=\left(m_\pi ^2+2m_{s\overline{s}}^2\right)/3$$
(6)
In the quenched case the $`a_0`$ mesons can couple to two-meson states through a ”hairpin diagram” on one of the meson lines. Such diagrams, like Fig. 1(b) in Ref. , can behave like powers of $`t`$ times $`e^{2m_\pi t}`$ and therefore masquerade as a light $`a_0`$ when $`2m_\pi <m_{a_0}`$. This may explain the lightest quark mass quenched point. In this analysis we used the $`a_0`$ from the local source, or the “$`\mathrm{𝟏}\mathrm{𝟏}`$” operator, which gave the best signal. Within the very large statistical errors, we saw relatively large breaking of flavor symmetry among the different $`a_0`$ channels, with some evidence that this reflects the masses of the different lattice pseudoscalars to which the various $`a_0`$’s should couple.

### D Baryons

We have evaluated propagators for baryons using the “corner wall” source for both degenerate and nondegenerate quarks, using a pointlike sink operator with all three quarks on the same lattice site. With nondegenerate quarks, the lightest states in this channel for zero, one or two strange quarks are the $`N`$, the $`\mathrm{\Lambda }`$ and the $`\mathrm{\Xi }`$ respectively. However, since we took no special measures to make the operator with one strange quark and two light quarks orthogonal to the $`\mathrm{\Sigma }`$, these propagators undoubtedly contain contamination from a nearby $`\mathrm{\Sigma }`$. In order to get the decuplet baryons, we followed Ref. , using a wall source on every spatial site, and using the operator in Eq. 6.3 of Ref.  for the $`\mathrm{\Delta }`$. This is necessary because the corner wall source does not overlap this $`\mathrm{\Delta }`$ operator. As a byproduct of the calculation of the decuplet mass, we obtain nucleon propagators from even site wall sources and a wall source containing all sites. These propagators are generally noisier than the corner wall source propagators, and the plateau in the effective mass occurs at larger distances. Therefore we generally use the corner source propagators. However, if one uses our fitting procedures on the even-wall or full-wall propagators, one invariably selects a smaller mass than from the corner-wall source propagators. This situation is illustrated by Fig. 17, which show mass fits to these two propagators as a function of the minimum distance included in the fit. In this figure the corner source propagators reach a plateau earlier and with smaller error bars. However, there are perfectly acceptable fits to the full wall source propagators giving masses significantly smaller than the corner source values. Perhaps the only good thing we can say about this situation is that the effect is similar in the quenched and dynamical runs, so as long as we are careful to make the same choices in both cases, we can investigate the effects of sea quarks on the spectrum while taking the statistical errors at face value. The results of the fits that we selected are listed in Tables VIII and IX.

The nucleon to rho mass ratio, or “Edinburgh plot”, has long been used as a way of displaying lattice spectrum results. This ratio is known to be sensitive to lattice spacing, lattice volume and quark masses. As mentioned above, we do not address the issues of continuum extrapolation and chiral extrapolation in this paper. However, since our quenched and dynamical lattices are matched in lattice spacing and physical size, we are well positioned to ask if effects of dynamical quarks show up in this ratio. In Fig. 18 we show a variant of the Edinburgh plot, the “APE plot”, where $`(m_\pi /m_\rho )^2`$ is used as the abscissa, so that for small quark masses the abscissa is proportional to the quark mass. Most of the points on this plot are from the $`a0.13`$ fm matched lattice runs. It can be seen that there are no significant differences between the quenched and three flavor runs. The single two flavor point lies slightly above the trend, although this is probably not significant. This agreement between $`m_N/m_\rho `$ for quenched and full QCD is in apparent conflict with our extrapolations of the conventional action. The discrepancy, which may be due to residual discretization effects in one or both calculations, is under study. We are hopeful that currently running three flavor simulations at finer lattice spacing will shed new light here. Figure 18 also contains a point from a coarser lattice three flavor run and a preliminary point from a finer lattice quenched run ($`10/g^2=8.4`$, $`28^3\times 96`$ lattice, $`a0.09`$ fm) These two points suggest that when we are in a position to do a continuum extrapolation the continuum results will be lower.

It is interesting to compare these results to the conventional Kogut-Susskind quark action. Figure 19 shows the improved action results together with conventional action results at $`6/g^2=5.7`$, $`5.85`$. $`6.15`$ and $`6.5`$, which correspond to lattice spacings of about $`0.16`$, $`0.12`$, $`0.07`$ and 0.043 fm respectively. While a continuum extrapolation will be deferred until the 0.1 fm runs are completed, we can see in this plot that the improved action at 0.13 fm gives results similar to the conventional action at $`a=0.07`$ fm.

Just as for the mesons, one of the biggest problems in comparing baryon masses to the real world is the need for an extrapolation in quark mass. For the mesons the quantity $`J`$ has the nice feature that it is only minimally sensitive to this extrapolation. It is tempting to try to construct similar quantities for the baryons. This suggests looking at the $`\mathrm{\Omega }^{}`$, the $`sss`$ decuplet baryon. More generally, we could plot the mass of the decuplet baryon as a function of $`(m_{PS}/m_V)^2`$ to produce a variant of the Edinburgh plot which has the nice feature that there are two experimental points, one for $`m_\mathrm{\Omega }^{}/m_\varphi `$ and another for $`m_\mathrm{\Delta }/m_\rho `$. Unfortunately, because the decuplet masses are hard to determine on the lattice, the errors are rather large. The result of this exercise is in Fig. 20. Although the error bars and the scatter among the points are large, the overall trend of this plot is encouraging. It is premature to say whether this plot shows real differences between full and quenched QCD.

### E Nonzero momentum

As a check on the quality of the dynamics in our simulations, we calculated the energies of a few nonzero momentum mesons, namely the Goldstone pion with momenta $`\frac{2\pi }{L}(0,0,1)`$ and $`\frac{2\pi }{L}(0,1,1)`$ and the $`\gamma _z\gamma _z`$ rho with momentum $`\frac{2\pi }{L}(0,0,1)`$. We compare these energies to the ideal dispersion relation by tabulating the “speed of light” in Table X.

$$c^2=\frac{E(\stackrel{}{k})E(\stackrel{}{0})}{k^2}$$
(7)
We see that the dispersion relation is generally very good, but with noticeable deviations from one for the heavier mesons. The results are similar in the quenched and full QCD calculations, as can be seen by comparing lines with the same valence quark mass in Table X.

## VI Conclusion

We have used simulations with three flavors of dynamical quarks and a quenched simulation on lattices with matched lattice spacing and physical size to isolate effects of the sea quarks on the hadron spectrum and on the static quark potential. This was done with a Symanzik improved gauge action and an improved Kogut-Susskind quark action to make the effects of the nonzero lattice spacing as small as practical. Effects of the sea quarks are clearly visible in the static potential. The mesonic mass ratio $`J`$ is much closer to the experimental value when dynamical quarks are included. The $`a_0`$ meson couples strongly to two meson states, as expected when sea quarks are included.

Several aspects need further study. While the quenched pion mass has the expected form, we have not understood the dependence of the three flavor pion mass on the quark masses. It would be very nice to be able to extract excited state masses, especially in the avoided level crossing of the $`a_0`$, and to extract decay rates to be compared to experiment. Although we have minimized the lattice artifacts by using an improved action, an empirical investigation of these effects is necessary. We are beginning a series of simulations with $`a0.09`$ fm to investigate these effects.

## acknowledgements

Computations were done on the Origin 2000 clusters at LANL, on the T3E and SP3 at NERSC, on the T3E, SP2 and SP3 at SDSC, on the NT and Linux clusters and Origin 2000 at NCSA, on the Origin 2000 at BU, on the Roadrunner and Los Lobos clusters at UNM and on the PC clusters at NERSC, Indiana University and the University of Utah. This work was supported by the U.S. Department of Energy under contracts DOE – DE-FG02-91ER-40628, DOE – DE-FG03-95ER-40894, DOE – DE-FG02-91ER-40661, DOE – DE-FG02-97ER-41022 and DOE – DE-FG03-95ER-40906 and National Science Foundation grants NSF – PHY99-70701 and NSF – PHY97–22022.


text:
# Asymmetries in ϕ photoproduction and the Okubo-Zweig-Iizuka violation

## Abstract

We study the vector meson density matrix and the associated polarization asymmetries in $`\varphi `$ photoproduction near threshold. In order to study the OZI violating processes, we consider the direct $`\varphi NN`$ coupling as well as the knockout processes arising from the non-vanishing $`s\overline{s}`$ sea quarks in the nucleon. The polarization asymmetries in $`\varphi `$ meson decays are found to be very useful to constrain the $`\varphi NN`$ coupling constant. The hidden strangeness of the nucleon can also be studied by these asymmetries, but mostly at large $`|t|`$ region.

preprint: nucl-th/0104068

The electromagnetic production of $`\varphi `$ meson from the nucleon has been suggested as a sensitive probe to study the hidden strangeness of the nucleon . This is because the $`\varphi `$ is nearly pure $`s\overline{s}`$ state so that its direct coupling to the nucleon is suppressed by the OZI rule. However, if there exists a non-vanishing $`s\overline{s}`$ sea quark component in the nucleon, the strange sea quarks can contribute to the $`\varphi `$ production via the OZI evasion processes. Investigation of such processes can then be expected to shed light on the nucleon strangeness content, if any .

In a series of publication, it has been argued that the double polarization asymmetries of $`\varphi `$ photo- and electroproduction from the nucleon targets at threshold energy region can be used to probe the hidden strangeness content of the nucleon . In addition to the direct knockout production mechanisms, the direct $`\varphi NN`$ coupling can also induce an OZI evading process . Even if we assume the ideal mixing of the $`\varphi `$ with the $`\omega `$, the effective $`\varphi NN`$ coupling is allowed by the kaon loops and hyperon excitations. But the smallness of the effective $`\varphi NN`$ coupling constants estimated by this approach reveals the validity of the OZI rule . However the analyses of electromagnetic nucleon form factors  and the baryon-baryon scattering  suggest strongly a fairly large $`\varphi NN`$ vector coupling constant which violates the OZI rule. The recent CLAS experiments on $`\varphi `$ photoproduction at large momentum transfers , therefore, open another way to constrain the $`\varphi NN`$ coupling. The data show increasing differential cross sections at large $`|t|`$, which has been interpreted as the contributions from the $`u`$-channel nucleon exchanges . It has been recently re-analysed to probe the nucleon resonances which couple to the $`\varphi N`$ channel in an SU(3) quark model . The recent CLAS data on $`\varphi `$ photo- and electroproduction  also lead to a discussion on the role of the axial ghost pole in vector meson electromagnetic production at low energies and large momentum transfers .

In this paper, we re-analyse the CLAS data by extending the work of Refs. . As the production models, we consider the Pomeron exchange, $`\pi `$ and $`\eta `$ exchange, nucleon pole terms with adjustable $`\varphi NN`$ couplings, and the direct knockouts. We study the vector meson density matrix and the associated polarization asymmetries in $`\varphi `$ photoproduction and find that they are useful to investigate the $`\varphi NN`$ coupling as well as the knockout mechanisms.

We start with defining the kinematic variables. The four-momenta of the incoming photon, outgoing $`\varphi `$ meson, and initial/final proton are denoted by $`k`$, $`q`$, $`p`$, and $`p^{}`$, respectively. The Mandelstam variables are defined as $`W^2=(p+k)^2`$ and $`t=(pp^{})^2`$. The photon laboratory energy is denoted by $`E_\gamma `$. We assume that the production amplitudes include the Pomeron exchange \[Fig. 1(a)\] and the pseudoscalar-meson $`(\pi ,\eta )`$ exchange \[Fig. 1(b)\]. With the above $`t`$-channel background processes we study the direct and crossed nucleon pole terms \[Fig. 1(c,d)\] as well as the direct knockout mechanisms \[Fig. 1(e,f)\].

For the Pomeron exchange, which governs the total cross sections and differential cross sections at low $`|t|`$, we follow the Donnachie-Landshoff model , which gives

$$T_{fi}^P=iM_0(s,t)\overline{u}_{m_f}(p^{})\epsilon _\mu ^{}(\varphi )\left\{k/g^{\mu \nu }k^\mu \gamma ^\nu \right\}\epsilon _\nu (\gamma )u_{m_i}(p),$$
(1)
where $`\epsilon _\mu (\varphi )`$ and $`\epsilon _\nu (\gamma )`$ are the polarization vectors of the $`\varphi `$ meson and photon, respectively, and $`u_m(p)`$ is the Dirac spinor of the nucleon with momentum $`p`$ and spin projection $`m`$. The Pomeron exchange is described by the following Regge parameterization:

$$M_0(s,t)=C_VF_1(t)F_V(t)\left(\frac{s}{s_0}\right)^{\alpha _P(t)1}\mathrm{exp}\left\{\frac{i\pi }{2}[\alpha _P(t)1]\right\},$$
(2)
where $`F_1(t)`$ is the isoscalar electromagnetic form factor of the nucleon and $`F_V(t)`$ is the form factor for the vector-meson–photon–Pomeron coupling:

$$F_1(t)=\frac{4M_N^22.8t}{(4M_N^2t)(1t/t_0)^2},F_V(t)=\frac{1}{1t/M_V^2}\frac{2\mu _0^2}{2\mu _0^2+M_V^2t},$$
(3)
where $`t_0=0.7`$ GeV2 and $`M_N`$ ($`M_V`$) is the nucleon ($`\varphi `$ meson) mass. The Pomeron trajectory is known to be $`\alpha _P(t)=1.08+0.25t`$. The strength factor $`C_V`$ reads $`C_V=12\sqrt{4\pi \alpha _{\mathrm{em}}}\beta _0\beta _s/f_V`$ with the vector meson decay constant $`f_V`$ ($`=13.13`$ for the $`\varphi `$ meson) and $`\alpha _{\mathrm{em}}=e^2/4\pi `$. By fitting all of the total cross section data for $`\omega `$, $`\rho `$, and $`\varphi `$ photoproduction at high energies, the remaining parameters of the model are determined: $`\mu _0^2=1.1`$ GeV2, $`\beta _0=2.05`$ GeV-1, $`\beta _s=1.60`$ GeV-1, and $`s_0=4`$ GeV2.

The pseudoscalar-meson exchange is known to be an important correction to the Pomeron exchange at low energies. The amplitude can be calculated from the following effective Lagrangians,

$`_{\varphi \gamma \mathrm{\Pi }}`$ $`=`$ $`{\displaystyle \frac{eg_{\varphi \gamma \mathrm{\Pi }}}{M_V}}ϵ^{\mu \nu \alpha \beta }_\mu \varphi _\nu _\alpha A_\beta \mathrm{\Pi },`$ (4)
$`_{\mathrm{\Pi }NN}`$ $`=`$ $`{\displaystyle \frac{g_{\pi NN}}{2M_N}}\overline{N}\gamma ^\mu \gamma _5\tau _3N_\mu \pi ^0+{\displaystyle \frac{g_{\eta NN}}{2M_N}}\overline{N}\gamma ^\mu \gamma _5N_\mu \eta ,`$ (5)
where $`\mathrm{\Pi }=(\pi ^0,\eta )`$ and $`A_\beta `$ is the photon field. The resulting invariant amplitude reads

$`T_{fi}^{ps}={\displaystyle \underset{\mathrm{\Pi }=\pi ,\eta }{}}{\displaystyle \frac{iF_{\mathrm{\Pi }NN}(t)F_{\varphi \gamma \mathrm{\Pi }}(t)}{tM_\mathrm{\Pi }^2}}{\displaystyle \frac{eg_{\varphi \gamma \mathrm{\Pi }}g_{\mathrm{\Pi }NN}}{M_V}}\overline{u}_{m_f}(p^{})\gamma _5u_{m_i}(p)\epsilon ^{\mu \nu \alpha \beta }q_\mu k_\alpha \epsilon _\nu ^{}(\varphi )\epsilon _\beta (\gamma ).`$ (6)
In the above, we have followed Ref.  to include the form factors to dress the $`\mathrm{\Pi }NN`$ and $`\varphi \gamma \mathrm{\Pi }`$ vertices as

$$F_{\mathrm{\Pi }NN}(t)=\frac{\mathrm{\Lambda }_\mathrm{\Pi }^2M_\mathrm{\Pi }^2}{\mathrm{\Lambda }_\mathrm{\Pi }^2t},F_{\varphi \gamma \mathrm{\Pi }}(t)=\frac{\mathrm{\Lambda }_{\varphi \gamma \mathrm{\Pi }}^2M_\mathrm{\Pi }^2}{\mathrm{\Lambda }_{\varphi \gamma \mathrm{\Pi }}^2t}.$$
(7)
We use $`g_{\pi NN}^2/4\pi =14`$ for the $`\pi NN`$ coupling constant. The $`\eta NN`$ coupling constant is not well determined . In the case of $`\omega `$ photoproduction, the higher mass of the $`\eta `$ and the smallness of the $`\omega \gamma \eta `$ coupling constant suppress the contribution of the $`\eta `$ exchange compared with that of the $`\pi `$ exchange . However in the case of $`\varphi `$ photoproduction, the $`\eta `$ exchange contribution may be larger than the $`\pi `$ exchange at large $`|t|`$ depending on the $`\eta NN`$ coupling constant because of the large value of the $`\varphi \gamma \eta `$ coupling. The dependence of the cross sections on $`g_{\eta NN}`$ was discussed in detail by Ref.  and in this study we use $`g_{\eta NN}^2/4\pi =0.99`$ which is obtained from making use of the SU(3) symmetry relation together with a recent value of $`F/D=0.575`$. The coupling constants $`g_{\varphi \gamma \mathrm{\Pi }}`$ can be estimated through the decay widths of $`\varphi \gamma \pi `$ and $`\varphi \gamma \eta `$ , which lead to $`g_{\varphi \gamma \pi }=0.14`$ and $`g_{\varphi \gamma \eta }=0.704`$. The amplitudes are also dependent of the cutoff parameters $`\mathrm{\Lambda }_\mathrm{\Pi }`$ and $`\mathrm{\Lambda }_{\varphi \gamma \mathrm{\Pi }}`$ in Eq. (7), which are fitted by the experimental data  at small and medium $`|t|`$ region and from Ref. : $`\mathrm{\Lambda }_\pi =0.9`$ GeV, $`\mathrm{\Lambda }_{\varphi \gamma \pi }=0.9`$ GeV, $`\mathrm{\Lambda }_\eta =1.0`$ GeV and $`\mathrm{\Lambda }_{\varphi \gamma \eta }=0.9`$ GeV.

We evaluate the nucleon pole terms shown in Fig. 1(c,d) from the following interaction Lagrangians,

$`_{\gamma NN}`$ $`=`$ $`e\overline{N}\left(\gamma _\mu {\displaystyle \frac{1+\tau _3}{2}}A^\mu {\displaystyle \frac{\kappa _N}{2M_N}}\sigma ^{\mu \nu }_\nu A_\mu \right)N,`$ (8)
$`_{\varphi NN}`$ $`=`$ $`g_{\varphi NN}\overline{N}\left(\gamma _\mu \varphi ^\mu {\displaystyle \frac{\kappa _\varphi }{2M_N}}\sigma ^{\mu \nu }_\nu \varphi _\mu \right)N,`$ (9)
with the anomalous magnetic moment of the nucleon $`\kappa _{p(n)}=1.79`$ $`(1.91)`$. There are big uncertainties with the $`\varphi NN`$ coupling constants. First, the effective $`\varphi NN`$ coupling through kaon loops and hyperon excitations was estimated to be very small, $`g_{\varphi NN}0.24`$ , which is consistent with the OZI rule prediction. However the analyses on the electromagnetic nucleon form factors  and the baryon-baryon scattering  favor a large $`\varphi NN`$ vector coupling constant which strongly violates the OZI rule: $`g_{\varphi NN}/g_{\omega NN}=0.30.43`$, therefore $`g_{\varphi NN}=2.34.7`$ with $`g_{\omega NN}=7.011.0`$ . In our study the nucleon pole terms are responsible to the recent CLAS data on $`\varphi `$ photoproduction at large $`|t|`$ and we fit the data with the $`g_{\varphi NN}`$ and find that $`|g_{\varphi NN}|=3.0`$ can reproduce the data, which is consistent with the Regge trajectory study of Ref. . The tensor coupling constant $`\kappa _\varphi `$ has been estimated to be $`\kappa _\varphi /\kappa _\omega =0.52.1`$  or $`\kappa _\varphi \pm 0.2`$ . However, it is known that the value of the tensor coupling $`\kappa _\omega `$ is very small as discussed in a study of $`\pi N`$ scattering and pion photoproduction . Thus we use $`\kappa _\varphi =0`$ in our calculation. We also found that the results are not dependent on $`\kappa _\varphi `$ if $`\kappa _\varphi <1`$. The resulting invariant amplitude reads

$$T_{fi}^N=\overline{u}_{m_f}(p^{})\epsilon ^\mu (\varphi )M_{\mu \nu }\epsilon ^\nu (\gamma )u_{m_i}(p),$$
(10)
where

$`M_{\mu \nu }=eg_{\varphi NN}\left[\mathrm{\Gamma }_\mu ^\varphi (q){\displaystyle \frac{\overline{)}p+\overline{)}k+M_N}{sM_N^2}}\mathrm{\Gamma }_\nu ^\gamma (k)F_N(s)+\mathrm{\Gamma }_\nu ^\gamma (k){\displaystyle \frac{\overline{)}p\overline{)}q+M_N}{uM_N^2}}\mathrm{\Gamma }_\mu ^\varphi (q)F_N(u)\right]`$ (11)
with

$$\mathrm{\Gamma }_\mu ^\varphi (q)=\gamma _\mu i\frac{\kappa _\varphi }{2M_N}\sigma _{\mu \alpha }q^\alpha ,\mathrm{\Gamma }_\nu ^\gamma (k)=\gamma _\nu +i\frac{\kappa _p}{2M_N}\sigma _{\nu \beta }k^\beta ,$$
(12)
and $`s=(p+k)^2`$, $`u=(pq)^2`$. Here we follow Ref.  to include a form factor

$`F_N(r)={\displaystyle \frac{\mathrm{\Lambda }_N^4}{\mathrm{\Lambda }_N^4(rM_N^2)^2}}`$ (13)
with $`r=s`$ or $`t`$. The gauge invariance is restored by projecting out the gauge non-invariant parts  and we use the cutoff parameter $`\mathrm{\Lambda }_N=0.5`$ GeV as in the study of $`\omega `$ photoproduction .

If we assume a non-vanishing strange sea quarks in the nucleon wavefunction, it can contribute to $`\varphi `$ photoproduction through direct knockout mechanisms  \[Fig. 1(e,f)\].The direct knockout mechanisms may be related to the loop diagrams of the effective $`\varphi NN`$ vertex considered in Ref. , i.e., the meson clouds of the nucleon or the intermediate baryons may interact with the photon to produce the $`\varphi `$ meson. Since we are considering tree diagrams only, however, clarifying that issue is beyond the scope of this work. In order to investigate the effects from the hidden strangeness content of the proton in $`\varphi `$ photoproduction, we parameterize the proton wavefunction in Fock space as

$$|p=A_0|uud+\underset{X}{}A_X|uudX+\underset{X}{}B_X|uuds\overline{s}X,$$
(14)
where $`X`$ denotes any combination of gluons and light quark pairs of $`u`$ and $`d`$ quarks. Ellis et al.  estimated it to be 1–19% from an analysis of $`p\overline{p}`$ annihilation. From the $`\varphi `$ electroproduction process, Henley et al.  claimed that its theoretical upper-bound would be 10–20%. By employing a relativistic quark model  it was shown that the upper-bound could be lowered to 3–5%. For simplicity and for our qualitative study, we approximate the proton wavefunction (14) as

$`|p`$ $`=`$ $`A|uud+B|uuds\overline{s}`$ (15)
$`=`$ $`A|[uud]^{1/2}+{\displaystyle \underset{j_{s\overline{s}}=0,1;j_c}{}}b_{j_{s\overline{s}}}|\left[\mathbf{[}[uud]^{1/2}[𝐋]\mathbf{]}^{j_c}[s\overline{s}]^{j_{s\overline{s}}}\right]^{1/2},`$ (16)
where $`𝐋`$ is the relative angular momentum of the two clusters. This parameterization of the nucleon wavefunction can be justified in $`\varphi `$ production as argued in Refs. . The strangeness content $`B^2`$ should be determined by further theoretical and experimental studies and we present our results with varying this value. The production amplitude of this mechanism has been discussed extensively in Refs.  and will not be repeated here.

Within the models discussed above, we present the differential cross sections of $`\varphi `$ photoproduction at $`E_\gamma =2.0`$ GeV and $`3.6`$ GeV in Fig. 2. In the considered energy region, it can be seen that the Pomeron exchange (dashed lines) gives the dominant contribution at small $`|t|`$ region. The pseudoscalar-meson exchange (dotted line) is important at medium $`|t|`$ region while the nucleon pole terms are responsible for the increase of the differential cross sections at large $`|t|`$. With $`|g_{\varphi NN}|=3.0`$, we could reproduce the recent CLAS data on $`\varphi `$ photoproduction . This is consistent with the results of Ref.  based on Regge trajectory. The direct knockout mechanism is suppressed and thus cannot be distinguished from the other mechanisms unless we study other quantities such as the polarization asymmetries .

Although the analyses on electromagnetic nucleon form factor and baryon-baryon scattering supports the large value of the $`\varphi NN`$ coupling, one can further constrain the coupling in $`\varphi `$ photoproduction. The electromagnetic nucleon form factor study based on vector meson dominance only tells that the relative phase between the $`\varphi \gamma `$ coupling and the $`\varphi NN`$ coupling is negative. We found that the asymmetries in $`\varphi `$ photoproduction can give us a clue on the phase of the $`\varphi NN`$ coupling and therefore the vector meson dominance in electromagnetic nucleon form factor can be tested directly. (In our convention, the $`\varphi \gamma `$ coupling has positive sign.) For this purpose, we compute the $`\varphi `$ meson density matrix elements within our model. The results at $`E_\gamma =2.8`$ GeV in helicity frame  are shown in Fig. 3. The solid lines are obtained with $`g_{\varphi NN}=3.0`$ while the dashed lines are with $`g_{\varphi NN}=+3.0`$. The contributions from the direct knockout mechanisms are small and cannot be seen here. Figure 3 shows that the vector meson density matrix is strongly dependent on the phase of the $`g_{\varphi NN}`$ at large $`|t|`$.

This difference can be readily seen in the photon polarization asymmetry $`\mathrm{\Sigma }_\varphi `$ and the parity asymmetry $`P_\sigma `$ which are defined as

$`\mathrm{\Sigma }_\varphi `$ $``$ $`{\displaystyle \frac{\sigma _{}\sigma _{}}{\sigma _{}+\sigma _{}}}={\displaystyle \frac{1}{p_\gamma }}{\displaystyle \frac{W^L(0,\frac{\pi }{2},\frac{\pi }{2})W^L(0,\frac{\pi }{2},0)}{W^L(0,\frac{\pi }{2},\frac{\pi }{2})+W^L(0,\frac{\pi }{2},0)}},`$ (17)
$`P_\sigma `$ $``$ $`{\displaystyle \frac{\sigma ^N\sigma ^U}{\sigma ^N+\sigma ^U}}=2\rho _{11}^1\rho _{00}^1,`$ (18)
where $`W^L`$ is the decay angular distribution of the $`\varphi `$ meson, $`\sigma _{}`$ ($`\sigma _{}`$) is the cross sections for symmetric decay particle pairs produced parallel (normal) to the photon polarization plane and $`\sigma ^N`$ ($`\sigma ^U`$) is the contribution of natural (unnatural) parity exchanges to the cross section. The details on the definitions and physical meanings of those asymmetries can be found, e.g., in Refs. . Shown in Fig. 4 are our results for the asymmetries at $`E_\gamma =3.6`$ GeV. As can be seen in Fig. 4, measuring those asymmetries can be used to constrain the $`\varphi NN`$ coupling and can test the vector meson dominance hypothesis.

As discussed in Refs.  and can be verified from Fig. 2, however, the direct knockout mechanisms could be investigated only at lower energies, i.e., near the threshold. Therefore, we compute the polarization asymmetries at $`E_\gamma =2.0`$ GeV in Fig. 5. The dashed lines are obtained without the knockout mechanisms, i.e., with $`B^2=0`$, while the solid and dot-dashed lines are with $`B^2=1`$ % and $`2`$ %, respectively. The results are dependent on the phases of $`b_{0,1}`$ and the phases of $`b_0`$ and $`b_1`$ are assumed to be positive in Fig. 5 for simplicity. We have used $`g_{\varphi NN}=3.0`$ in the calculation. The difference of the results depending on $`B^2`$ can be seen especially at large $`|t|`$ region. Although the difference is not as large as in Fig. 4, it can be used for to study the $`uud`$ knockout process \[Fig. 1(f)\] and might be tested experimentally.

In summary, we have investigated the vector meson density matrix and the associated polarization asymmetries in $`\varphi `$ photoproduction at low energies. As the models for $`\varphi `$ photoproduction, we have considered the Pomeron and pseudoscalar $`(\pi ,\eta )`$ exchanges, nucleon pole terms, and the direct knockouts. We found that the recent CLAS data on $`\varphi `$ photoproduction can be explained by a rather large $`\varphi NN`$ coupling, i.e. $`g_{\varphi NN}=3.0`$, which supports the result of Ref. . This value is larger than the expectations based on the OZI rule, but is consistent with the values estimated in electromagnetic nucleon form factors and baryon-baryon scattering. As another test for constraining $`g_{\varphi NN}`$ we computed the vector meson density matrix and the polarization asymmetries. It was shown that they are strongly dependent on the $`\varphi NN`$ coupling constant. The results are different from the predictions of Ref.  which includes the nucleon resonances decaying into $`\varphi N`$ channel within an SU(3) quark model. Thus the polarization asymmetries can be used to test the models for the nucleon pole terms and the nucleon resonances. Finally, we have shown that they can be used to probe the hidden strangeness content of the nucleon. Those observables could be measured at current experimental facilities such as Thomas Jefferson National Accelerator Facility and LEPS of SPring-8. On the theoretical side, further investigations are needed for clarifying the model dependence of our results and understanding the contributions from the loop diagrams which give an effective $`\varphi NN`$ coupling through the rescattering processes.

###### Acknowledgements.

We are grateful to V. Burkert, M. Fujiwara, T.-S. H. Lee, and E. Smith for fruitful discussions. This work was supported in part by the Brain Korea 21 project of Korean Ministry of Education and the International Collaboration Program of KOSEF under Grant No. 20006-111-01-2.


text:
# REFERENCES

Reply: Van den Brink, Khaliullin and Khomskii

In the preceding Comment , Shen points out that the on-site Coulomb interaction, that can cause charge order in half-doped manganites , also destabilizes the magnetic CE-phase observed in these systems.

This is a valid observation, but it is not a priori clear whether in the relevant parameter regime the C-phase is indeed lower in energy then the CE-phase within our model . From the exact diagonalization of a 16-site ring, Shen finds a groundstate energy for the CE-phase as a function of $`U`$, which is compared to the energy of an infinite ring in the C-phase, in order to find a critical coupling $`U_c=5t`$ above which the C-phase is stable. In order to determine the critical coupling, however, one should handle finite size effects in both phases on equal footing and compare the CE-phase energy with the energy of a 16-site C-phase ring, which is $`E_{16}^C=0.628t`$. Using this value for the energy, we can estimate from Fig.1 in Ref.  that $`U_c10t`$, which is in the relevant parameter regime for the half-doped manganites.

As Shen states, other interactions which are neglected in the model could very well compensate for the small energy differences between the C- and CE-phase. We wish to point out in particular that coherent and incoherent hopping between two sites with opposite spin orientation lowers the kinetic energy of the CE-phase appreciably.

The main aim of our model is to explain the observed charge and orbital order, given the CE magnetic structure. As to the stability of the phases, indeed not all physical relevant effects are taken into account. Especially the contribution of lattice distortions, which we neglected, are important . Also, for instance, anharmonic effects that favor occupation of elongated orbitals  that may stabilize the CE-phase, are not incorporated in the model.

Moreover, one should also take into account that the hopping between two sites with antiferromagnetically coupled core spins is not completely blocked. For a finite Hund’s rule coupling ($`J_H`$) an electron can hop coherently from a high spin state to lower spin states on a neighboring site, but even for a infinitely large $`J_H`$, there exist quantum effects which allow hopping that preserve the total high spin state $`S+1/2`$, but with different projections $`S^z=S+1/2`$ and $`S1/2`$ at different sublattices . Such electron hopping processes between C-type chains or between CE zigzag chains modify the total energy. One can easily estimate the coherent contribution to the kinetic energy in second order perturbation theory. The contribution from hopping between two chains in the C-phase is then $`\delta E_C=\frac{3}{8}\frac{t^2}{\mathrm{\Delta }}`$ and in the CE-phase $`\delta E_{CE}=\frac{5}{8}\frac{t^2}{\mathrm{\Delta }}`$. This difference is due to the fact that in the CE-phase also the $`x^2y^2`$ orbital on the corner site is partially occupied. As typically $`\mathrm{\Delta }J_H`$ and $`J_H`$ is in the order of $`5t`$, the energy lowering of the CE-phase due to this effect would be larger than the energy loss due to the on-site Coulomb interactions in the CE-phase. The incoherent contributions that change the spin projections on different sublattices are more involved to calculate, but they would certainly lower the energy of the CE-phase with respect to the C-phase as the inter-chain hopping integrals are larger in the former.

We therefore conclude that the model of Ref. , which correctly captures the interplay of spin, charge and orbital degrees of freedom in the half-doped manganites and gives a reasonable description of their electronic structure, is by itself not sufficient for the precise determination of the regions of stabilities of different phases. For this several other factors, in particular those mentioned above, should be taken into account.

Jeroen van den Brink
Faculty of Applied Physics and MESA<sup>+</sup> Institute
University of Twente
Box 217, 7500 AE Enschede, The Netherlands

Giniyat Khaliullin
Max-Plack-Institute für Festkörperforschung
Heisenbergstrasse 1, 70569, Stuttgart, Germany

Daniel Khomskii
Laboratory of Solid State Physics
University of Groningen
Nijenborgh 4, 9747 AG Groningen, The Netherlands