text:
# Freezing of a Stripe Liquid

## Abstract

The existence of a stripe-liquid phase in a layered nickelate, La1.725Sr0.275NiO4, is demonstrated through neutron scattering measurements. We show that incommensurate magnetic fluctuations evolve continuously through the charge-ordering temperature, although an abrupt decrease in the effective damping energy is observed on cooling through the transition. The energy and momentum dependence of the magnetic scattering are parametrized with a damped-harmonic-oscillator model describing overdamped spin-waves in the antiferromagnetic domains defined instantaneously by charge stripes.

One of the key issues in current debates over copper-oxide superconductors concerns the nature and relevance of charge stripes Emery et al. (1999). It has been demonstrated in one cuprate family that the holes doped into the CuO2 planes can order in an array of periodically-spaced stripes, separating antiferromagnetic domains Ichikawa et al. (2000). For charge stripes to be relevant to superconductivity, they must be ubiquitous among the cuprates, and for the latter to be true, they must be able to exist in a liquid state. Indeed, theoretical models for such electronic liquid crystal phases have been proposed Zaanen et al. (1996); Kivelson et al. (1998). The inelastic incommensurate magnetic scattering observed in La2-xSrxCuO4 Cheong et al. (1991); Yamada et al. (1998) and YBa2Cu3O6+x Mook et al. (1998) has sometimes been interpreted as evidence for dynamic stripes; however, this interpretation has been controversial Bourges et al. (2000).

Here we present evidence for a stripe liquid phase in a related system, La2-xSrxNiO4. Ordered stripes have been observed in this system over a large range of hole concentrations Yoshizawa et al. (2000). Although the maximum charge and spin ordering temperatures occur for $`x=0.33`$ Lee and Cheong (1997), we have chosen to study crystals of $`x=0.275`$, which have the advantage that the spin and charge-ordering wave vectors do not coincide. Using neutron scattering, we follow the magnetic inelastic scattering to temperatures as high as 1.7 times the charge-ordering transition, $`T_{\mathrm{co}}`$, and demonstrate the existence of a stripe liquid. (La2-xSrxNiO4 is nonmetallic even at $`TT_{\mathrm{co}}`$, as indicated by a peak in its optical conductivity at $`0.5`$ eV Katsufuji et al. (1996); Pashkevich et al. (2000), so there can be little controversy over possible alternative interpretations involving Fermi-surface nesting.) We show that the $`\omega `$ and $`𝐐`$ dependence of the data can be effectively parametrized with a damped-harmonic-oscillator (DHO) model describing overdamped spin waves associated with the antiferromagnetic domains defined instantaneously by the charge stripes. Using a spin-wave energy dispersion with an effective gap energy, we find that the gap energy is proportional to the inverse correlation length and varies linearly with temperature through $`T_{\mathrm{co}}`$. On the other hand, the behavior of the damping parameter changes abruptly near $`T_{\mathrm{co}}`$.

The neutron scattering measurements were performed on the SPINS triple-axis spectrometer in the cold-neutron guide hall at the NIST Center for Neutron Research (NCNR). Initially, a single crystal of 4.4 g was used; later, a second crystal of 6.3 g was mounted along side the first. Both crystals were grown at Kyoto University, and annealed at the University of Delaware in order to achieve a stoichiometric oxygen concentration. (The first crystal was used previously for the work reported in Lee et al. (2001), but without annealing.) The neutron spectrometer was equipped with a vertically-focusing monochromator and horizontally-focusing analyzer, both utilizing PG (002) reflections. Elastic scans of superlattice peaks were performed with 5-meV neutrons. Inelastic measurements were done in energy-gain mode, with incident energies of either 5 or 13.7 meV. The sample temperature was controlled with a displex refrigerator.

The stripe order in this sample is consistent with that observed in previous studies of La2-xSrxNiO4 Yoshizawa et al. (2000). Figure 1(a) shows the temperature dependence of the spin and charge-order superlattice peaks. The spin order approaches zero at $`T_{\mathrm{so}}120`$ K, while the charge order disappears by $`T_{\mathrm{co}}190`$ K.

Our focus here is on the dynamic correlations, especially those at $`T>T_{\mathrm{co}}`$. Figure 2 shows examples of constant-energy-transfer scans ($`\mathrm{}\omega =4`$ meV) through the incommensurate magnetic peak positions $`(1,\pm ϵ,0)`$ at temperatures (a) $`T>T_{\mathrm{co}}`$, (b) $`T_{\mathrm{co}}>T>T_{\mathrm{so}}`$, and (c) $`T_{\mathrm{so}}>T`$. We clearly observe well-resolved peaks at all temperatures. Given that the ordered state at $`T<T_{\mathrm{co}}`$ is already well characterized in terms of periodically spaced charge stripes separating antiferromagnetic domains, the continuous evolution of the inelastic magnetic scattering \[see Fig. 1(b)\] provides direct evidence for the existence of instantaneously correlated magnetic domains in the disordered state. The fact that the scattering at the commensurate antiferromagnetic wave vector (1,0,0) is always a local minimum indicates that neighboring antiferromagnetic domains maintain their antiphase relationship, which is caused by the segregation of the doped holes to the domain walls. Thus, we feel that Fig. 2(a) is firm evidence for a stripe-liquid phase.

The variation of the $`Q`$-widths of the peaks in the constant-$`E`$ scans for $`\mathrm{}\omega =0`$, $`4`$, and $`8`$ meV is plotted in Fig. 3. For $`TT_{\mathrm{so}}`$, the widths are roughly temperature independent but vary substantially with energy. The variation with energy is similar to what one might expect from spin-wave dispersion, as observed previously Tranquada et al. (1997). At higher temperatures, the widths appear to vary linearly with temperature, and the dependence on frequency is reduced.

The distribution of the magnetic-scattering strength with frequency also evolves with temperature. Figure 2(d)-(f) show the variation of the imaginary part of the dynamic susceptibility, $`\chi ^{\prime \prime }(𝐐,\omega )`$, at several temperatures; here, $`𝐐=𝐐_0`$, where $`𝐐_0`$ stands for an incommensurate magnetic wave vector. $`\chi ^{\prime \prime }`$ is related to the experimentally measured $`S(𝐐,\omega )`$ via

$$S(𝐐,\omega )=\left(1e^{\mathrm{}\omega /k_\mathrm{B}T}\right)^1\chi ^{\prime \prime }(𝐐,\omega );$$
(1)
the correction for the detailed-balance factor was applied after subtracting the background, measured at $`𝐐=(1,0.9,0)`$. The position of the maximum of $`\chi ^{\prime \prime }`$ as a function of energy corresponds to a characteristic damping energy $`\mathrm{\Gamma }`$, and one can see that $`\mathrm{\Gamma }`$ increases with temperature.

To describe the data more quantitatively, we have chosen to use the damped-harmonic-oscillator (DHO) model:

$$\chi ^{\prime \prime }(𝐐,\omega )=\underset{𝐐_0}{}\frac{2\omega \gamma \chi _0}{(\omega ^2\omega _𝐐^2)^2+(2\omega \gamma )^2},$$
(2)
where the excitations are assumed to have the dispersive form

$$\omega _𝐐^2=\omega _0^2+c^2(𝐐𝐐_0)^2.$$
(3)
(Direct evidence for dispersive modes at high energies will be presented in a complementary study Bourges et al. .) To fit the data, the model $`S(𝐐,\omega )`$ was convolved with the spectrometer resolution function. The fitting parameters are $`\omega _0`$, $`\gamma `$, and $`\chi _0`$, with the effective spin-wave velocity held fixed at $`\mathrm{}c=300`$ meV Å, the approximate value determined for pure La2NiO4 Yamada et al. (1991). To describe the $`𝐐`$ dependence of the data, the relevant combinations are $`\omega _0/c`$ and $`\gamma /c`$; for the energy dependence, as we will discuss, it is only the ratio $`\omega _0^2/2\gamma `$ that matters.

The solid curves in Fig. 2 represent the fits to the data. The parameter values obtained from such fits are displayed in Fig. 4. For $`T>T_{\mathrm{so}}`$, $`\omega _0`$ varies linearly with temperature. In contrast, $`\gamma `$ is roughly constant up to $`TT_{\mathrm{co}}`$ where it begins to increase rapidly. $`\chi _0`$ has a fairly weak temperature dependence which is described approximately by $`A/(1+\gamma /B)`$, where $`B120`$ meV. The temperature dependence of $`\gamma `$ suggests that a major contribution to this damping factor comes from the fluctuations of the charge stripes Zaanen and van Saarloos (1997); it is consistent with the growth in the inverse correlation length for charge stripes observed by x-ray scattering in La2-xSrxNiO4 with $`x=\frac{1}{3}`$ Du et al. (2000). The saturation of $`\gamma `$ below $`T_{\mathrm{co}}`$ might be due in part to the quenched disorder associated with the randomly-distributed Sr dopants.

Our scans versus frequency (Fig. 2) are measured at $`𝐐=𝐐_0`$, and most of the data correspond to $`\omega ^2\omega _0^2`$, in which case Eq. (2) simplifies to

$$\chi ^{\prime \prime }(𝐐_0,\omega )\frac{\chi _0}{2\gamma }\frac{\omega }{\omega ^2+\mathrm{\Gamma }^2},$$
(4)
where $`\mathrm{\Gamma }=\omega _0^2/2\gamma `$. The values of $`\mathrm{\Gamma }`$ obtained from the values of the fitted parameters are plotted in Fig. 3(b). At $`T>T_{\mathrm{co}}`$, $`\mathrm{\Gamma }`$ changes slowly and has a value of roughly 7 meV. On cooling through $`T_{\mathrm{co}}`$, $`\mathrm{\Gamma }`$ abruptly starts to decrease, and becomes rather small by $`T_{\mathrm{so}}`$. Since $`\omega _0`$ varies smoothly through this regime, the change in behavior at $`T_{\mathrm{co}}`$ is controlled by $`\gamma `$. If $`\gamma `$ is a measure of charge-stripe fluctuations, as discussed above, then it appears that the abrupt decrease in $`\mathrm{\Gamma }`$ below $`T_{\mathrm{co}}`$ may be a direct result of charge order. Low-frequency probes such as nuclear magnetic resonance Abu-Shiekah et al. (1999); Yoshinari et al. (1999) and muon-spin relaxation Chow et al. (1996); Jestädt et al. (1999) should be sensitive to the variations in $`\mathrm{\Gamma }`$, and it appears that such variations above $`T_{\mathrm{so}}`$ have been detected in studies of related nickelates Abu-Shiekah et al. (1999); Chow et al. (1996); Jestädt et al. (1999). The sensitivity of magnetic damping to charge order may also be relevant to the mechanism of the “wipe-out” effect observed in nuclear-quadrupole-resonance studies in cuprates Hunt et al. (2001); Curro et al. (2000); Julien et al. (2001).

In order to get further insight into the significance of the model parameters, we numerically integrated $`S(𝐐,\omega )`$ over frequency. The result corresponds to $`𝐐`$-dependent peaks with a line shape that is approximately Lorentzian. The half-width-at-half-maximum of this function should correspond to the instantaneous inverse correlation length, $`\xi ^1\kappa `$, and we find that, to within 1%, $`\kappa =0.86\omega _0/c`$. This result, that the effective spin-excitation gap $`\omega _0`$ is approximately equal to $`\kappa c`$, is equivalent to the form proposed for the paramagnetic state of the undoped Heisenberg antiferromagnet Auerbach and Arovas (1988); Mata and Arnold (1989).

It is of interest to consider the transport properties of La2-xSrxNiO4 in the regime $`T>T_{\mathrm{co}}`$ which we now associate with the stripe-liquid phase; numerous single-crystal studies of transport and optical properties have been reported previously (see, e.g., Katsufuji et al. (1996); Pashkevich et al. (2000); Katsufuji et al. (1999)). Figure 5 compares the temperature-dependence of the in-plane resistivity measured on our crystal with similar measurements on two La2-xSrxCuO4 samples Ichikawa et al. (2000); Wakimoto et al. (2001). We observe that the resistivity of our nickelate crystal is only an order of magnitude greater than that of underdoped cuprates at 300 K, and we expect that it will develop a “metallic” temperature derivative at higher temperatures, as observed above 400 K in La2-xSrxNiO4 with $`x=0.33`$ Katsufuji et al. (1996). Furthermore, at high temperatures the resistivities of all three samples exceed a critical limit; to the extent that $`d\rho /dT`$ is positive, they qualify as “bad” metals Emery and Kivelson (1995). We believe it is plausible that the “bad” metallic behavior of the cuprates might be associated with a stripe-liquid phase. (Resistivity in stripe-ordered cuprates can be low Ichikawa et al. (2000).)

The comparable magnitudes of room-temperature resistivity in the nickelates and cuprates make it plausible that a stripe-liquid phase could be relevant to transport properties in both materials. There are certainly differences in the two systems: the stripes inevitably order in the nickelates, whereas stripe order is generally avoided in the cuprates. However, these differences might be associated more with the magnitude of fluctuations rather than the nature of instantaneous correlations. In La2-xSrxCuO4 with $`x=0.14`$, the measurements of Aeppli and coworkers Aeppli et al. (1997) indicate that, using the present parametrization, $`\mathrm{}\mathrm{\Gamma }10`$ meV at $`T=35`$ K, with $`\mathrm{\Gamma }`$ increasing substantially at higher temperatures. In that case, the relevance of quantum critical phenonema has been proposed Castellani et al. (1997); Sachdev (2000).

To conclude, we have presented experimental evidence for the existence of a liquid phase of charge stripes in a hole-doped nickelate. The magnetic correlations evolve smoothly through the freezing transition, although the damping of the fluctuations does not. The in-plane resistivity, though non-metallic, is relatively low in the stripe-liquid phase.

We gratefully acknowledge helpful comments from S. A. Kivelson and J. Zaanen. KY acknowledges support from the Japan Science and Technology Corporation, the Core Research for Evolutional Science and Technology Project (CREST), and Grants-in-Aid for Scientific Research on Priority Areas, 12046239, 2001 and Research (A), 10304026, 2001 and for Creative Scientific Research (13NP0201) from the Japanese Ministry of Education, Culture, Sports, Science, and Technology. Work at SPINS is based upon activities supported by the National Science Foundation under Agreement No. DMR-9986442. Research at Brookhaven is supported by the Department of Energy’s (DOE) Office of Science under Contract No. DE-AC02-98CH10886. DJB acknowledges support from DOE under Contract No. DE-FG02-00ER45800.


format: text

text:
# Complete Factorization of Equations of Motion in Wess-Zumino Theory

## Abstract

We prove that the equations of motion describing domain walls in a Wess-Zumino theory involving only one chiral matter multiplet can be factorized into first order Bogomol’nyi equations, so that all the topological defects are of the Bogomol’nyi-Prasad-Sommerfield type.

PACS numbers: 11.27.+d, 11.30.Pb, 11.30.Er

Domain walls are defect structures that appear in diverse branches of physics. They usually live in three spatial dimensions as bidimensional objects that arise in systems described by potentials that contain at least two isolated degenerate minima. They involve energy scales as different as the ones that appear in Condensed Matter  and in Cosmology .

A lot of attention has been drawn recently to domain walls in field theories, in models that have been investigated under several distinct motivations . A very specific motivation concerns the presence of domain walls arising in between non zero vacuum expectation values of scalar fields in supergravity . Another line deals with the formation of defects inside domain walls . A great deal of attention has also been drawn to $`SU(N)`$ supersymmetric gluodynamics, where nonperturbative effects give rise to gluino condensates that may form according to a set of $`N`$ isolated degenerate chirally asymmetric vacua, from where domain walls spring interpolating between pairs of vacua . The interest in domain walls in general widens because of the interplay between Field Theory and the low energy world volume dynamics of branes in String Theory .

In the present work we examine the bosonic portion of a Wess-Zumino theory that involves only one chiral matter multiplet. We write the Lagrangian density in the standard form

$`={\displaystyle \frac{1}{2}}_\alpha \overline{\phi }^\alpha \phi {\displaystyle \frac{1}{2}}|W^{}(\phi )|^2.`$ (1)
We assume that the superpotential $`W(\phi )`$ is a holomorphic function in the complex variable $`\phi `$, and $`W^{}(\phi )`$ stands for the complex derivative of $`W`$ with respect to $`\phi `$. The specific form of the potential shows that its minima are critical values of the superpotential.

We search for defect structures, for static solutions of the equation of motion. We suppose the static solutions depend only on $`x`$, so the equation of motion becomes

$$\frac{d^2\phi }{dx^2}=W^{}(\phi )\overline{W^{\prime \prime }(\phi )}.$$
(2)
We examine the energy of the static configurations to show that the second order ordinary differential equation is equivalent to a family of first order ordinary differential equations of Bogomol’nyi type, namely

$$\frac{d\phi }{dx}=\overline{W^{}(\phi )}e^{i\xi },$$
(3)
where $`\xi `$ is the family parameter ranging in the set of real numbers, constant with respect to $`x`$. Thus, we say that the equation of motion (2) completely factorizes into a family of first order equations. This is our main result. It is valid under the boundary conditions $`lim_x\mathrm{}\phi (x)=v^k`$, $`lim_x\mathrm{}(d\phi /dx)=0`$, as required by the topological solutions, where $`v^k`$ is some critical point of the superpotential $`W`$. A direct consequence of this result is that in such models all the topological solutions are of the BPS type, that is, these models do not support non-BPS states.

Next, we proceed to prove our main result, that is, the equivalence of the equation of motion (2) and the family of first order equations (3). It is obvious that (3) implies (2), for any $`\xi `$, as one can see just by differentiating equation (3) with respect to $`x`$. To prove the reverse, that is, to prove that any solution of (2) under the aforementioned boundary conditions is a solution of (3) for some $`\xi `$, we introduce the ratio

$`R(\phi )={\displaystyle \frac{1}{\overline{W^{}(\phi )}}}{\displaystyle \frac{d\phi }{dx}}.`$ (4)
Differentiating $`R(\phi )`$ with respect to $`x`$ and using the equation of motion (2), we obtain

$`{\displaystyle \frac{dR(\phi )}{dx}}`$ $`=`$ $`\left({\displaystyle \frac{1}{\overline{W^{}(\phi )}}}\right)^2\left[\overline{W^{}(\phi )}{\displaystyle \frac{d^2\phi }{dx^2}}{\displaystyle \frac{d\phi }{dx}}\overline{W^{\prime \prime }(\phi )}{\displaystyle \frac{d\overline{\phi }}{dx}}\right]`$ (5)
$`=`$ $`\left({\displaystyle \frac{1}{\overline{W^{}(\phi )}}}\right)^2\left[\left|W^{}(\phi )\right|^2\left|{\displaystyle \frac{d\phi }{dx}}\right|^2\right]\overline{W^{\prime \prime }(\phi )}.`$ (6)
We notice that

$`{\displaystyle \frac{d}{dx}}|W^{}(\phi )|^2=\overline{W^{}(\phi )}W^{\prime \prime }(\phi ){\displaystyle \frac{d\phi }{dx}}+W^{}(\phi )\overline{W^{\prime \prime }(\phi )}{\displaystyle \frac{d\overline{\phi }}{dx}}={\displaystyle \frac{d}{dx}}\left|{\displaystyle \frac{d\phi }{dx}}\right|^2.`$ (7)
Besides, due to the boundary conditions we also have $`lim_x\mathrm{}W^{}(\phi )=0`$. Thus

$$|W^{}(\phi )|^2\left|\frac{d\phi }{dx}\right|^2=0,$$
(8)
and so, $`R(\phi )`$ does not depend on $`x`$. It remains only to show that $`R(\phi )`$ has modulus one. This is a consequence of (8) again. Indeed,

$`|R(\phi )|^2=R(\phi )\overline{R(\phi )}={\displaystyle \frac{1}{|W^{}(\phi )|^2}}\left|{\displaystyle \frac{d\phi }{dx}}\right|^2=1.`$ (9)
We now present our last remarks. First, our result is valid in the case of Wess-Zumino theory involving only one chiral matter multiplet, so it does not work for models with more than one multiplet. Second, in systems of scalar fields one may occasionally run into solutions that engender no topological feature. They are named nontopological solutions, and have been found for instance in Ref. . We first notice that in systems of real scalar fields described by some superpotential, the nontopological solutions cannot appear as solutions of the first order Bogomol’nyi equations, because they should have zero energy, and this is the energy of the vacuum states. We add this to our former result to obtain another result: in the bosonic sector of a Wess-Zumino theory involving only one chiral matter multiplet there is no room for non topological solutions. Third, if we deal with a single real scalar field $`\varphi `$, we have to consider potentials in the form $`V(\varphi )=(1/2)W_\varphi ^2`$. In this case the equation of motion for $`\varphi =\varphi (x)`$ becomes $`d^2\varphi /dx^2=W_\varphi W_{\varphi \varphi }`$. It also factorizes into the two first order Bogomol’nyi equations $`d\varphi /dx=\pm W_\varphi `$, for solutions that obey $`lim_x\mathrm{}\varphi (x)v^k`$ and $`lim_x\mathrm{}(d\varphi /dx)=0`$, where $`v^k`$ is a minimum of $`V(\varphi )`$, a critical point of $`W(\varphi )`$.

In summary, the new idea of factorizing equations of motion into first order Bogomol’nyi equations is of direct interest to supersymmetry. It can be used in diverse applications, in particular to investigate the presence of BPS walls interpolating between distinct pairs of vacua, and exact integrability. The complete factorization of equations of motion in systems described by $`W(\phi )`$ holomorphic ensures BPS feature to domain walls that spring in such systems. Because of the associated BPS character, these domain walls partially preserve the supersymmetry.

We thank F.A. Brito, J.R. Morris, J.R.S. Nascimento, and R.F. Ribeiro for discussions, and CAPES, CNPq and PRONEX for partial support.


text:
# A new fermion Hamiltonian for lattice gauge theory

## Abstract

We formulate Hamiltonian vector-like lattice gauge theory using the overlap formula for the spatial fermionic part, $`H_f`$. We define a chiral charge, $`Q_5`$ which commutes with $`H_f`$, but not with the electric field term. There is an interesting relation between the chiral charge and the fermion energy with consequences for chiral anomalies.

Lattice gauge theory and chiral symmetry represent two venerable non-perturbative approaches in particle theory. Historically, however, they have not easily meshed. In particular, the old species doubling problem forces one to insert chiral symmetry breaking somewhere into the regulator. The issues involved are deeply entwined with anomalies . Stimulated by the domain wall fermion idea of Kaplan  and the overlap formulation of Neuberger and Narayanan , this topic has recently seen renewed attention and dramatic progress.

Here we explore an adaptation of the overlap idea to the Hamiltonian formalism. We find an interesting operator structure and a simple intuitive picture of how anomalies work in terms of a fermion eigenvalue flow.

To start, consider the conventional “continuum” Hamiltonian for a gauge theory

$$\begin{array}{c}H=H_f+H_g\\ H_f=\overline{\psi }D_c\psi =\psi ^{}\gamma _0D_c\psi \\ H_g=E^2+B^2\end{array}$$
(1)
where the fermion fields $`\psi `$ anticommute and

$$D_c=\stackrel{}{\gamma }(\stackrel{}{}+ig\stackrel{}{A})+m.$$
(2)
The operator $`D_c`$ consists of an antihermitean kinetic piece and the hermitean mass term (we use hermitean gamma matrices). The antihermitean part anticommutes with both $`\gamma _0`$ and $`\gamma _5`$, giving the properties

$$\gamma _5D_c=D_c^{}\gamma _5,\gamma _0D_c=D_c^{}\gamma _0.$$
(3)
The eigenvalues of $`D_c`$ lie along the line $`\mathrm{Re}\lambda =m`$ and occur in complex conjugate pairs; i.e. if we have $`D_c\chi =\lambda \chi `$, then $`D_c\gamma _5\chi =\lambda ^{}\gamma _5\chi `$.

The free theory is particularly simple in momentum space, with $`D_c=i\stackrel{}{p}\stackrel{}{\gamma }+m`$. The eigenvalues are $`\lambda =\pm i|\stackrel{}{p}|+m`$. This spectrum is discrete in finite volume, with momentum quantized in units of $`\frac{2\pi }{L}`$, where $`L`$ is the size of our box.

Now for the lattice, “naive” fermions replace the momentum with a trigonometric function $`p_i\mathrm{sin}(p_ia)/a`$. This causes the well known doubling, with extra low energy states when $`p_i\frac{\pi }{a}`$. The Wilson solution to this problem gives the mass a momentum dependence $`mm+\frac{1}{a}_i(1\mathrm{cos}(p_ia))`$ so that the doubler mass becomes of $`O(1/a)`$. The free Wilson-Dirac operator in the Hamiltonian formulation is thus

$$D_w=m+\frac{1}{a}\underset{i}{}\left(i\mathrm{sin}(p_ia)\gamma _i+1\mathrm{cos}(p_ia)\right)$$
(4)
The eigenvalues of $`D_w`$ lie on a superposition of ellipses, as sketched here

To simplify notation, from now on we work in lattice units with $`a=1`$. For the continuum limit we are interested in small eigenvalues $`\lambda 0`$.

When gauge fields are added, several things happen. First the fermi eigenvalues move around, starting to fill the holes in the above figure. In the process the operator $`D_w`$ ceases to be normal, i.e $`[D_w,D_w^{}]0`$. The eigenvectors are no longer orthogonal. Also the pairing due to the $`\gamma _5`$ symmetry is lost. This is directly related to the chiral symmetry violation in the Wilson term.

To get a nicer behavior, we mimic ref.  and project $`D_w`$ onto a unitary matrix

$$V=D_w(D_w^{}D_w)^{1/2}$$
(5)
Being unitary, $`V^{}V=1`$, this is a normal matrix. From this we construct the operator

$$D=1+V$$
(6)
and our new fermion Hamiltonian

$$\begin{array}{c}H_f=\psi ^{}h\psi \\ h=\gamma _0D=\gamma _0(1+V)\end{array}$$
(7)
We recover the continuum properties

$$[D,D^{}]=0,\gamma _5D=D^{}\gamma _5$$
(8)
and the eigenvalues again lie in complex conjugate pairs. By construction, the eigenvalues lie on a circle, as sketched here

As the above figure is drawn, there are no low energy states. For such, we must make the starting Wilson mass $`m`$ negative, keeping the doubler masses positive. Then the figure is more like this

As with the euclidian overlap operator, combining the unitarity condition $`V^{}V=1=D^{}DDD^{}+1`$ with the hermiticity condition $`D^{}=\gamma _5D\gamma _5`$ gives rise to the Ginsparg-Wilson  relation

$$\gamma _5D+D\gamma _5D\gamma _5D=0$$
(9)
However we now have another variation on this following from $`D^{}=\gamma _0D\gamma _0`$

$$\gamma _0D+D\gamma _0D\gamma _0D=0$$
(10)
By either multiplying the first of these two Ginsparg-Wilson relations by $`\gamma _0`$ or the second by $`\gamma _5`$ we obtain the exact commutation relation

$$[\gamma _5(1\frac{D}{2}),h]=0$$
(11)
This suggests defining an axial charge

$$\begin{array}{c}Q_5\psi ^{}q_5\psi \\ q_5\gamma _5(1\frac{D}{2})=\gamma _5\frac{1V}{2}\end{array}$$
(12)
This charge commutes with the fermion part of our Hamiltonian (see also Ref. )

$$[Q_5,H_f]=0.$$
(13)
Since $`D`$ can be reconstructed from either $`h`$ or $`q_5`$, the matrices $`q_5`$ and $`h`$ are closely correlated

$$\begin{array}{c}D=\gamma _0h=22\gamma _5q_5\\ \gamma _5q_5+\gamma _0\frac{h}{2}=1\end{array}$$
(14)
Squaring this gives

$$q_5^2+\frac{h^2}{4}=1.$$
(15)
Since $`q_5`$ and $`h`$ commute they can be simultaneously diagonalized. The above equation states that these eigenvalues lie on a circle.

Low-energy states have a well defined chirality; i.e. $`h0q_5\pm 1`$. In contrast, high-energy states all have $`|q_5|<1`$. The combination $`\gamma _0\gamma _5`$ flips the sign of both eigenvalues, which are thus paired on opposite sides of the the circle.

This formulation gives a simple understanding of anomalies, analogous to the domain-wall discussion in . An adiabatic change of gauge fields shifts eigenvalues continuously while the pairing is preserved. A topologically non-trivial shift moves an $`h>0`$ eigenvalue to $`h<0`$ while its paired state goes the other way. In this way a “left” hole and “right” particle (or vice versa) are generated out of the Dirac sea, as sketched here

We have been considering the fermion states in a given background gauge field. For the full coupled quantum theory theory we add the gauge field Hamiltonian

$$H=H_f+E^2+B^2$$
(16)
The operator $`D`$ involves link variables, which do not commute with the electric field term, $`[Q_5,E^2]0`$. This is essential for the $`U(1)`$ anomaly to generate the $`\eta ^{}`$ mass via mixing with gluons.

There is a close connection between these zero crossings in energy and zero modes in Euclidean space. Adiabatically change the gauge field with time, $`\frac{d}{dt}h(t)=O(1/T)`$ with $`\frac{T}{2}<t<\frac{T}{2}`$. Then consider an eigenvalue $`h(t)\chi (t)=E_i(t)\chi (t)`$ which changes in sign, $`E(\frac{T}{2})<0`$, $`E(\frac{T}{2})>0`$. From this construct

$$\varphi (t)=e^{E(t)t}\chi (t)$$
(17)
This satisfies

$$D_4\varphi \gamma _0\left(_0+h(t)\right)\varphi (t)=0+O(1/T)$$
(18)
This wave function is normalizable if the eigenvalue rises through 0.

Because of $`\gamma _5`$ hermiticity, $`\gamma _5D_4\gamma _5=D_4^{}`$, complex eigenvalues of $`D_4`$ are paired; thus, unpaired zero modes are robust. This is the lattice version of the index theorem.

A variety of questions remain. One involves flavored axial charges such as $`Q_5^\alpha =\psi ^{}\tau ^\alpha \gamma _5(1\frac{D}{2})\psi `$. It appears that these also do not commute with the electric field term of the Hamiltonian, $`[Q_5^\alpha ,E^2]0`$. Why is this so, despite the euclidean overlap formulation having an exact flavored chiral symmetry? Another question arises at the level of currents. Since $`Q_5`$ is not ultra locally defined , what is the natural associated $`\stackrel{}{J}_5`$? Finally, our construction was for vectorlike theories. Can something similar be done for chiral gauge theories? In these cases anomalies change species as in the t’Hooft  baryon decay process. For the standard model, how do the quark and lepton “circles” interact to give this phenomenon?


text:
# Lattice study of nucleon properties with domain wall fermions

## Abstract

Domain wall fermions (DWF) are a new fermion discretization scheme with greatly improved chiral symmetry. Our final goal is to study the nucleon spin structure through lattice simulation using DWF. In this paper, we present our current progress on two topics toward this goal: 1) the mass spectrum of the nucleon excited states and 2) the iso-vector vector and axial charges, $`g_V`$ and $`g_A`$, of the nucleon.

1. Introduction

The RIKEN-BNL-Columbia-KEK QCD Collaboration has been pursuing the domain-wall fermion (DWF) method in lattice quantum chromodynamics (QCD). In DWF an extra fifth dimension is added to the lattice. By manipulating the domain-wall structure of the fermion mass in this fifth dimension, we control the number of light fermion species in the other four space-time dimensions. These light fermions possess exact chiral symmetry in the limit of an infinite fifth dimension. In particular: 1) fermion near-zero mode effects are well understood, 2) explicit chiral symmetry breaking induced by a finite extra dimension is described by a single residual mass parameter, which is very small in the present calculation, in the low-energy effective lagrangian, and 3) non-perturbative renormalization works well.

DWF is a promising new approach for treating fermions on the lattice. However, we need several tests of DWF in the baryon sector to reach our final goal of establishing the spin structure of the nucleon from first principles. Here we report our recent studies of the mass spectrum of the nucleon and its excited states and the nucleon matrix elements of the iso-vector vector and axial charges in quenched lattice QCD with DWF. Although most of the latter results are preliminary, the conclusive results in the former subject have been reported in Ref..

2. Nucleon excited states

First, we discuss the mass spectrum of the nucleon $`N`$ and its excited states (the negative-parity nucleon $`N^{}`$ and the positive-parity first excited nucleon $`N^{}`$) by means of a systematic investigation utilizing two distinct interpolating operators $`B_1^\pm `$ and $`B_2^\pm `$. For an explanation of those operators, see Ref.. Our quenched DWF calculation was employed on lattice with size $`16^3\times 32\times 16`$, gauge coupling $`\beta =6/g^2=6.0`$, and domain wall height $`M_5=1.8`$. Additional details of our simulation can be found in Ref..

In Fig.1 we show the low-lying nucleon spectrum as a function of the quark mass, $`m_f`$ in lattice units ($`a^11.9`$ GeV set from $`aM_\rho `$=0.400(8) in the chiral limit). $`B_1^+`$ gives the ground-state nucleon mass $`N`$ (cross). The $`N^{}`$ mass estimates (square and diamond) are extracted from both $`B_1^{}`$ and $`B_2^{}`$. The corresponding experimental values for $`N`$ and $`N^{}`$ are marked with lower and upper stars. Both $`N^{}`$ mass estimates extracted from two distinct operators agree with each other. The large $`N`$-$`N^{}`$ mass splitting is clearly evident.

In contrast to the negative parity operators, we find that the mass estimates from a second operator, $`B_2^+`$, are considerably larger than the ground state obtained from $`B_1^+`$. This suggests that $`B_2^+`$ has negligible overlap with the nucleon ground state and provides a signal for the positive-parity excited nucleon $`N^{}`$. To justify this possibility, we employ a sophisticated approach which utilizes the transfer matrix of a $`2\times 2`$ correlation function constructed from both $`B_1^+`$ and $`B_2^+`$. The diagonalization of the transfer matrix yields the excited state. Fig.2 shows a comparison of the fitted mass from $`B_2^+(t)\overline{B}_2^+(0)`$ (circle) and the estimated mass from the average effective mass given by the smaller eigenvalue of the transfer matrix (bullet). The cross symbol corresponds to the nucleon ground state mass evaluated from the larger eigenvalue of the transfer matrix, which is quite consistent with the fitted mass from $`B_1^+(t)\overline{B}_1^+(0)`$

Another important conclusion can be drawn from Fig.1 and Fig.2. In the heavy quark mass region, the ordering of the negative-parity nucleon ($`N^{}`$) and the positive-parity excited nucleon ($`N^{}`$) is inverted relative to experiment. This remarkable result was originally reported in our early paper and subsequently confirmed in Ref.. Further systematic calculation is required to determine whether this ordering switches to the observed ordering as one approaches the chiral limit.

3. Nucleon matrix elements

The nucleon (iso-vector) axial charge $`g__A`$ is a particularly interesting quantity. We know precisely the experimental value $`g__A=1.2670(35)`$ from neutron beta decay. Why does $`g__A`$ deviate from unity in contrast to vector charge, $`g__V=1`$? The simple explanation is given by the fact that the axial current is only partially conserved in the strong interaction while the vector current is exactly conserved. Thus, the calculation of $`g__A`$ is an especially relevant test of the chiral properties of DWF in the baryon sector. In addition, calculation of $`g__A`$ is an important first step in studying polarized nucleon structure functions since $`g__A=\mathrm{\Delta }u\mathrm{\Delta }d`$ where $`p,s|\overline{q}_f\gamma _5\gamma _\mu q_f|p,s=2s_\mu \mathrm{\Delta }q_f`$ with $`s^2=1`$ and $`sp=0`$.

We follow the standard practice for the calculation of $`g__V`$ and $`g__A`$. The three-point functions for the local vector current $`J__V^f=\overline{q}_f\gamma _4q_f`$ and axial current $`J__A^f=\overline{q}_f\gamma _5\gamma _iq_f`$ ($`i=1,2,3`$) are defined by

$$G_\mathrm{\Gamma }^{u,d}(t,t^{})=\mathrm{Tr}[P__\mathrm{\Gamma }\underset{\stackrel{}{x},\stackrel{}{x}^{}}{}TB_1(\stackrel{}{x},t)J__\mathrm{\Gamma }^{u,d}(\stackrel{}{x}^{},t^{})\overline{B}_1(0,0)]$$
(1)
with $`\mathrm{\Gamma }`$ = $`V`$ (vector) or $`A`$ (axial) where $`P__V=P_+=(1+\gamma _4)/2`$ and $`P__A=P_+\gamma _i\gamma _5`$. For the axial current, the three-point function is averaged over $`i=1,2,3`$. The lattice estimates of vector and axial charges can be derived from the ratio between two- and three-point functions

$$g__\mathrm{\Gamma }^{\mathrm{lattice}}=\frac{G_\mathrm{\Gamma }^u(t,t^{})G_\mathrm{\Gamma }^d(t,t^{})}{G__N(t)},$$
(2)
where $`G__N(t)=\mathrm{Tr}[P_+_\stackrel{}{x}TB_1(\stackrel{}{x},t)\overline{B}_1(0,0)]`$.

Recall that in general lattice operators $`𝒪_{\mathrm{lat}}`$ and continuum operator $`𝒪_{\mathrm{con}}`$ are regularized in different schemes. The operators are related by a renormalization factor $`Z_𝒪`$: $`𝒪_{\mathrm{con}}(\mu )=Z_𝒪(a\mu )𝒪_{\mathrm{lat}}(a)`$. This implies that the continuum value of vector and axial charges are given by $`g__\mathrm{\Gamma }=Z__\mathrm{\Gamma }g__\mathrm{\Gamma }^{\mathrm{lattice}}`$. In the case of conventional Wilson fermions, the renormalization factor $`Z__A`$ is usually estimated in perturbation theory ($`Z__A`$ differs from unity because of explicit symmetry breaking). For DWF, the conserved axial current receives no renormalization. This is not true for the lattice local current. An important advantage with DWF, however, is that the lattice renormalizations, $`Z__V`$ and $`Z__A`$, of the local currents are the same so that the ratio $`(g__A/g__V)^{\mathrm{lattice}}`$ directly yields the continuum value $`g__A`$.

Our preliminary results are analyzed on 200 quenched gauge configurations at $`\beta =6.0`$ on a $`16^3\times 32\times 16`$ lattice with $`M_5=1.8`$. We choose a fixed separation in time of the nucleon interpolating operators, $`t=t_{\mathrm{source}}t_{\mathrm{sink}}`$ and $`t^{}<t`$, with currents inserted in between. We take $`t_{\mathrm{source}}=5`$ and $`t_{\mathrm{sink}}=21`$.

In Fig.3 we show the dependence of the vector renormalization, $`Z__V=1/g__V^{\mathrm{lattice}}`$ on the location of current insertions. A good plateau is observed so that $`Z__V`$ is certainly well behaved. The value 0.763(5) at $`m_f=0.02`$ (obtained by averaging over time slices denoted by the solid line in Fig.3) agrees well with $`Z__A=0.7555(3)`$, which was obtained from a completely different calculation of meson two-point correlation functions based on the relation $`A_\mu ^{\mathrm{conserved}}(t)\overline{q}\gamma _5q(0)=Z__AA_\mu ^{\mathrm{local}}(t)\overline{q}\gamma _5q(0)`$.

For the axial charge, $`g__A^{\mathrm{lattice}}`$, plateaus are seen for $`10t16`$ in Fig.4, so the charge ratios $`(g__A/g__V)^{\mathrm{lattice}}`$ at each $`m_f`$ are averaged in this time slice range. We find that there is a fairly strong dependence of $`(g__A/g__V)^{\mathrm{lattice}}`$ on $`m_f`$ as be shown in Fig.5. A simple linear extrapolation to $`m_f=0`$ yields $`g__A=0.62(13)`$, which is roughly 1/2 the experimental value. However, a simple linear ansatz may not describe the data in which case the result in the chiral limit may be even smaller.

To compare to results using Wilson fermions, we plot $`g__A`$ versus the square of the $`\pi `$-$`\rho `$ mass ratio in Fig.6. Our result is given by the bullets. We also include two heavier mass points (with large errors) from an earlier simulation using a larger separation between the source and the sink. The triangles and diamonds are from Wilson simulations at relatively strong and weak coupling. At first glance, the DWF and Wilson fermions seem to be in rough agreement except for the lightest point. However, our DWF results have a strong mass dependence. This may be a finite volume effect; the Wilson results at strong coupling were simulated on a lattice with a volume which is twice ours.

A couple of comments on the mass dependence of our data are in order. First, it is interesting to note that our results appear consistent with the value 5/3 (marked as star) in the heavy quark limit, while the others seem inconsistent with this limit. Second, our results may also be consistent with vanishing in the chiral limit. This can be explained through the axial Ward-Takahashi identity which governs $`g__A`$. If the PCAC relation $`m_\pi ^2m_f`$ is modified, for example by chiral logarithms, the nucleon matrix element of the pseudoscalar density does not develop a pole as $`m_f0`$, and the r.h.s of the identity vanishes in the chiral limit. Thus, $`g__A`$ must also vanish as $`m_f0`$. Indeed, we already know that the PCAC relation for the pion mass is modified in the quenched approximation by two effects: zero-modes of the Dirac operator and the quenched chiral logarithm. Further investigation of the Ward-Takahashi identity is under way.

4. Conclusions

We have explored several nucleon properties in quenched lattice QCD using domain wall fermions toward our final goal of studying the spin structure of the nucleon from first principles.

Our quenched DWF calculation reproduces very well the large mass splitting between the nucleon $`N(939)`$ and its parity partner $`N^{}(1535)`$. We have also calculated the mass of the first positive-parity excited state $`N^{}(1440)`$ by the diagonalization of a $`2\times 2`$ matrix correlator and confirmed that it is heavier than the negative-parity excited state $`N^{}(1535)`$. A remaining puzzle is whether or not a switching of $`N^{}`$ and $`N^{}`$ occurs close to the chiral limit.

A preliminary calculation of iso-vector vector and axial charges shows that all the relevant three-point functions are well behaved. $`Z__V`$ determined from the nucleon matrix element of the vector current agrees closely with that from an NPR study of quark bilinears and a direct calculation using meson correlation functions. This indicates that $`g__V=1`$ and $`Z__V=Z__A`$ are mutually satisfied in our quenched DWF calculation. However, a linear extrapolation of $`g__A`$ to the chiral limit yields a value which is a factor of two smaller than the experimental value. We are currently investigating the Ward-Takahashi identity that governs $`g__A`$ to shed light on this behavior. We also plan to check related systematic effects arising from finite lattice volume and quenching (for example quenched chiral logarithms, zero modes, and the absence of the full pion cloud), especially in the lighter quark mass region.

Acknowledgments

The author would like to thank the organizers of NSTAR2001 for an invitation. This work was done in collaboration with Tom Blum and Shigemi Ohta as a part of the RIKEN-BNL-Columbia-KEK QCD Collaboration. We thank RIKEN, Brookhaven National Laboratory, and the U.S. Department of Energy for providing the facilities essential for the completion of this work.


text:
# The Mauna Kea Observatories Near-Infrared Filter Set. II. Specifications for a New 𝐽𝐻𝐾𝐿'𝑀' Filter Set for Infrared Astronomy

## 1 Introduction

The rationale for a new set of infrared filters was presented by Simons & Tokunaga (2001; hereafter Paper I). The goals of the design of this new filter set are similar to those of Young, Milone, & Stagg (1994), namely to construct a filter set that minimized color terms in the transformations between photometric systems and that reduced the uncertainty in determining the absolute calibration of photometric systems. However, in contrast to the filters proposed by Young et al., the filters in this new set were also designed to maximize throughput, in addition to minimizing the effects of atmospheric absorption. The reduced dependence on atmospheric absorption also permits photometry to be less sensitive to water vapor variations and to the altitude of the observatory. Most importantly, the new filters permit extrapolation to zero air mass with small errors.

We required large-size, high-quality filters for several facility instruments for the Gemini and Subaru Telescopes. Since the cost of each filter is dominated by the technical difficulties of the coating process it is desirable to place as many substrates as possible into the coating chamber in order to reduce the cost per filter. Hence there is a strong economic driver to produce custom filters in a consortium production run. This has the added benefit that with more observatories involved there would be greater standardization among the observatories.

Although this filter set was designed for the Gemini and Subaru Telescopes, all of the optical/infrared observatories at Mauna Kea are presently using these filters (NASA Infrared Telescope Facility, United Kingdom Infrared Telescope, Canada-France-Hawaii Telescope, Keck, Gemini, Subaru). In addition, the filter set was discussed informally by the IAU Working Groups on IR Photometry and on Standard Stars at the 2000 General Assembly, and endorsed as the preferred “standard” near-infrared photometric system, to be known as the Mauna Kea Observatories Near-Infrared (MKO-NIR) photometric system.

In this paper we describe a set of filters that were fabricated according to the definitions presented in Paper I. We include a list of specifications so that similar filters can be produced in the future. It is our hope that if enough observatories adopt these filters, we will have greater uniformity among photometric systems, as well as reduced systematic errors when comparing observations from different observatories.

## 2 Filter Specifications

We present here the list of technical specifications that were required to be satisfied by the filter manufacturer. The center, cut-on, and cut-off wavelengths are given in Table 1, and they follow the filter definitions given in Paper I. The specifications for substrate flatness, parallelism of the filters, and use of a single substrate are required for use with adaptive optics systems.

1. Out-of-band transmission: $`<`$10-4 out to 5.6 $`\mu `$m.

This specification is required for use with InSb detectors. For HgCdTe 2.5 $`\mu `$m cut-off detectors, the out-of-band blocking can be specified for wavelengths less than 3.0 $`\mu `$m instead of 5.6 $`\mu `$m. The desired blocking is better than 10-4, but practical considerations such as cost and manufacturing difficulty make it impossible to go lower. In cases where blocking to this level or better is not possible, a separate blocking filter should be used. PK-50 is a suitable blocker for wavelengths less than 2.0 $`\mu `$m. For longer wavelengths, an interference blocking filter may be required for InSb detectors.
2. Operating temperature: 65 K. Cold filter scans of witness samples to be provided, together with prediction of wavelength shift with temperature.

The filters are expected to be used at 50–77 K, but the specified filter temperature was set at 65 K primarily because the Gemini and Subaru instruments were designed for use with InSb detectors and with cryogenic motors inside the cryostat. This requires cooling to 65 K to avoid excess thermal emission from the motors within the cryostat.
3. Average transmission: $`>`$80% (goal $`>`$90%).
4. Transmission between the 80% points to have a ripple of less than $`\pm `$5%.
5. Cut-on and cut-off: $`\pm `$0.5% wavelength error.
6. Roll-off: %slope $``$2.5%.

The %slope $`=`$ $`[\lambda (80\%)\lambda (5\%)]/\lambda (5\%)\times 100`$, where $`\lambda (80\%)`$ is the wavelength at 80% transmission and $`\lambda (5\%)`$ the wavelength at 5% transmission.

Ideally the filters should have a square “boxcar” shape for maximum throughput. The specifications on the average transmission, ripple, cut-on and cut-off uncertainty, and roll-off uncertainty are a compromise between the desired sharp edge, practical manufacturing specifications, and cost.
7. Substrate surfaces parallel to $``$5′′.
8. Substrate flatness: $`<`$$`0.0138\lambda /(n1)`$, where $`n`$ is the index of refraction of the substrate (compatible with AO systems).

For example, for $`n=1.5`$, $`\lambda =2200`$ nm, the substrate flatness should be $`<`$61 nm. For $`n=3.4`$, $`\lambda =2200`$ nm, the substrate flatness should be $`<`$13 nm.

This is a flatness specification, not a roughness specification. The Strehl ratio (SR) can be approximated as SR $`1[(2\pi /\lambda )w]^2`$, where $`w`$ is the rms wavefront error (see Schroeder 1987). This expression is valid for a single reflective surface. For two surfaces in transmission, we have SR $`1[(n1)(2\pi /\lambda )w\sqrt{2}]^2`$, where $`n`$ is the index of refraction of the substrate. For SR = 0.985, $`w=0.01378[\lambda /(n1)]`$. Alternatively, one could specify $`\lambda /10`$ peak-to-valley at 0.63 $`\mu `$m, which is approximately $`\lambda /40`$ rms, or $``$15 nm rms. This would be suitable for glass substrates ($`n=1.5`$) at 1.0–2.2 $`\mu `$m as well as silicon substrates ($`n=3.4`$) at 2.2–4.8 $`\mu `$m. It should be noted that after coating, the substrate will deform under the stress of the coatings, taking on a concave or convex shape. However, this does not have any effect on the wavefront error in transmission since the surfaces are parallel.
9. Single-substrate to be used (cemented filters not acceptable).

Cemented filters are not considered acceptable because of potential large surface deformations when cooled. In cases where a high SR is not needed, cemented filters may be acceptable.

Items 7–9 are the main requirements for filters to be used in adaptive optics. The specified flatness conforms to an SR of 0.985.
10. Filter to be designed for a tilt of 5 (to suppress ghost images).

Items 9–10 assume that these filters will be placed at the pupil image within a camera and the filters are tilted to the optical axis by 5 to avoid ghost images in the focal plane. Alternatively, these filters could be designed to be used with no tilt, but with a wedged substrate.
11. Scratch/Dig: 40/20.

This item specifies the maximum size of scratches and digs (pits) that are permitted on the surface of the substrate. The specification 40/20, which is acceptable in most ground-based astronomy applications, requires that scratches be no wider than 4 $`\mu `$m and that the pits be no wider than 200 $`\mu `$m in diameter. A technical description of this specification is given by Bennett & Mattsson (1999).
12. Diameter: 60 mm, or cut to requested size.

The coating is typically not good near the edge and some of the filter is blocked by the filter mount. Therefore the filter diameter is typically specified to be about 10% larger than the actual size required by the optics.
13. Maximum thickness: 5 mm, including blocker for InSb detectors.

Items 12–13 depends on the application. The thickness/diameter ratio must be of a practical value for polishing and mechanical strength and is usually specified to be about 0.1.
14. No radioactive materials such as thorium to be used (to avoid spurious detector noise spikes).

## 3 Filter Fabrication and Filter Profiles

A number of vendors were contacted to solicit bids on the filter production for the consortium of buyers. The vendor selected was OCLI of Santa Rosa, California. This selection was determined both by price of production, ability to meet the specifications, maturity of the production facilities, and willingness to accept individual orders from the consortium. The filters were cut to size by the vendor for each order.

A production run was ordered to accommodate the filter needs of the following observatories and institutions: Anglo-Australian Observatory, California Institute of Technology, Canada-France-Hawaii Telescope, Carnegie Institution, Center for Astrophysics, Cornell University, European Southern Observatory, Gemini Telescopes, Kiso Observatory, Korean Astronomy Observatory, Kyoto University, MPE-Garching, MPI-Heidelberg, NASA Goddard Space Flight Center, NASA Infrared Telescope Facility, National Astronomical Observatory of Japan, National Optical Astronomy Observatories, Nordic Optical Telescope, Ohio State University, Osservatorio Astrofisico di Arcetri, Subaru Telescope, United Kingdom Infrared Telescope, University of California Berkeley, University of California Los Angeles, Univ. of Cambridge, University of Grenoble, Univ. of Hawaii, University of Wyoming, and the William Herschel Telescope.

Figures 1(a)–(g) show the transmission profiles of the filters that were produced. The profiles are shown at 65 K, as obtained by extrapolating the filter shift with temperature from room temperature to 65 K. Figure 2 shows the $`J,H,K_s,L^{}`$, and $`M^{}`$ filter profiles superimposed on the atmospheric transmission at Mauna Kea.

## 4 Photometric properties

Using a modified version of the ATRAN program (Lord 1992), we constructed models of the telluric absorption as a function of air mass to investigate the photometric properties of the new filter set. Absorption spectra covering the range 0.8–5.0 $`\mu `$m were generated for 10 zenith angles between 0 and 70 for each of five values of the precipitable water vapor content (0.5, 1.0, 2.0, 3.0, and 4.0 mm). The zenith angles correspond to a range of air masses between 1.0 and 2.9. The values of the water vapor content span the range typically found on Mauna Kea (see, e.g., Morrison et al. 1973; Warner 1977; Wallace & Livingston 1984). To examine the photometric characteristics down to an air mass of 0.0, we also calculated spectra for five additional air mass values $`<`$1.0 (0.05, 0.10, 0.25, 0.50, and 0.75) for each of the water vapor values. These spectra were generated with ATRAN using a model atmosphere in which both the gas and water vapor content were given by the standard model values (appropriate for Mauna Kea at an air mass of 1) scaled by the air mass value. The water content was then scaled by an additional factor given by the ratio of the desired water vapor value to the standard model value computed by ATRAN for the atmosphere above Mauna Kea. All spectra were computed for an altitude of 4158 m (which is the altitude of the NASA Infrared Telescope Facility) and with a resolution of approximately 100,000.

The telluric absorption spectra were then multiplied by a similar resolution model spectrum of Vega obtained from R. Kurucz (2001). For each filter, we computed synthetic magnitudes by multiplying the resulting spectra by the normalized filter transmission curves and integrating over the wavelength range spanned by the filter. Magnitude differences relative to the original model spectrum of Vega were calculated. In Fig. 3 we show the magnitude differences versus air mass for the case of 2.0 mm precipitable water.

For each filter and water vapor value, the run of magnitude differences as a function of air mass was fitted with a least-squares routine with two functional forms: (1) a linear fit between air masses of 1 and 3 and (2) the rational expression given by Young et al. (1994; see their eq. 7) for the entire range of air masses. The latter expression has the form

$$\mathrm{\Delta }m=(a+bX+cX^2)/(1+dX)$$
(1)
where $`a`$, $`b`$, $`c`$, and $`d`$ are fitting constants and $`X`$ is the air mass. A detailed discussion of this function is given by Young (1989). Both fits are shown in Fig. 3, and the coefficients are given in Tables 2 and 3.

The linear fit simulates the typical photometric reduction technique employed by ground-based observers to obtain the extinction coefficient in magnitude per air mass. The constant term (the intercept of the linear fit) yields the systematic error in the photometric magnitude incurred by adopting this particular functional form of the atmospheric extinction and extrapolating to zero air mass. Inspection of the values given in Table 2 indicates that systematic errors $`>`$0.05 mag result from a linear extrapolation for the $`K^{}`$ and $`M^{}`$ filters for all values of the water vapor. The errors are decreased only slightly if the fit is restricted to the air mass range 1–2.

As can be seen in Fig. 3, equation (1) provides an excellent fit to all the data points. However, the determination of the four coefficients of such a function from actual observations of standard stars is impractical for most observing programs, as a large number of measurements over a range of air masses are required. Even if these observations were carried out, reliable values can be obtained only if measurements in the completely inaccessible air mass range between 0 and 1 were available. The values in Table 3 are presented for comparison to other photometric filters, such as those discussed by Young et al. (1994).

The proposed filters show greatly improved photometric properties over older filter sets. Krisciunas et al. (1987) summarized the average extinction coefficients at Mauna Kea. Comparison of our values given in column 3 of Table 2 with those given in Table I by Krisciunas et al. (1987) reveals that at the typical water vapor of 2–3 mm of precipitable water for Mauna Kea, the reduction in the extinction coefficient is about a factor of 8 at J. In addition, comparison of the values given in column 2 of Table 2 with those given in Table III of Manduca & Bell (1979) shows that the reduction in the error of extrapolation to zero air mass (which they denote as $`\mathrm{\Delta }`$) is about a factor of 9 at J compared to the KPNO winter values.

As was pointed out by Manduca & Bell (1979), the extinction error upon extrapolation to zero air mass does not affect differential photometry, in which one compares observed fluxes of target objects to those of established standard stars. The extinction errors in this case will be absorbed into the photometric zero points. It should be noted, however, that the extinction error may vary with the temperature of the object being observed. This variation leads to a color term in the extinction coefficient. However, as demonstrated by Manduca & Bell, this term is generally small, amounting to an error of $`<`$0.001 mag for most filters. Only the $`J`$ band filter analyzed by Manduca & Bell exhibited a significant color dependence. Given that the new $`J`$ filter is less affected by the atmospheric absorption, we expect that the extinction color term for this filter will be as small as those for the other filters.

The extinction errors introduced by adopting a linear fit may, however, affect absolute photometric measurements as, for example, in the establishment of an exo-atmospheric magnitude system. The latter subject has a number of additional complications beyond the scope of the present paper. We simply wish to point out that such a system requires some way to eliminate the extinction error (for example, with narrowband filters; see Young et al. 1994).

## 5 Summary

A set of 1–5 $`\mu `$m filters designed for improved photometry and for use with adaptive optics is described. Although the filter profiles are optimized for maximum throughput, they avoid most of the detrimental effects of atmospheric absorption bands. A production run of these filters has been completed. It is our hope that widespread use of these filters will help produce more uniform photometric results, reduce the magnitude of the color transformation among observatories, and reduce the uncertainty in reductions to zero or unit air mass.

We thank E. Atad, J. Elias, S. Leggett, K. Matthews, and especially T. Hawarden for discussions and input on the filter specifications. We also thank S. Lord for providing a modified version of ATRAN that allowed the calculation of the atmospheric transmission between air masses of 0.0 and 1.0. Filter profiles in electronic form may be found at http://irtf.ifa.hawaii.edu/Facility/nsfcam/hist /newfilters.html. ATT acknowledges the support of NASA Contract NASW-5062 and NASA Grant No. NCC5-538. DS acknowledges the support of the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the Particle Physics and Astronomy Research Council (United Kingdom), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), CNPq (Brazil), and CONICET (Argentina).

Fig. 1(a)–(g). Filter profiles at a temperature of 65 K. These were measured by the vendor at room temperature, then shifted according to the measured change with temperature of a witness sample.

Fig. 2. $`J,H,K_s,L^{},M^{}`$ filter profiles superimposed on the atmospheric transmission at Mauna Kea kindly provided by G. Milone for 1 mm precipitable water vapor and an air mass of 1.0.

Fig. 3. Extinction plot for 2.0 mm precipitable water. The dashed lines are the linear fits to the magnitudes in the air mass range of 1.0–3.0. The solid lines are the fit to equation (1) in the air mass range of 0.0–3.0.


text:
# 1 The handbag-type diagram for leptoproduction of mesons. The large blob represents a sum over all spectator configuration. 𝑘 and 𝑘' denote the momenta of the active partons. The small blob stands for meson leptoproduction off partons.

Leptoproduction of Polarized Vector Mesons
B. Postler
Fachbereich Physik, Universität Wuppertal, D-42097 Wuppertal, Germany

## Abstract

We present a status report on a study of vector meson leptoproduction for the HERA energy range on the basis of the generalized handbag approach for which the amplitudes factorize in the parton subprocess $`\gamma ^{}gVg`$ and generalized parton distributions (full report in ). In contrast to the leading twist approach transverse degrees of freedom as well as Sudakov suppressions are taken into account in the subprocess. First results for the cross section of the reaction $`\gamma _L^{}p\rho _L^0p`$ are found to be in fair agreement with experiment.

We are going to investigate vector meson leptoproduction off protons at high energies, large virtuality, $`Q^2`$, of the exchanged photon and small invariant momentum transfer, $`t`$, from the initial to the final proton. In the kinematical region we are interested in, the process factorizes into a hard subprocess - meson leptoproduction off partons - and a soft proton matrix element which represents a generalized or skewed parton distribution (GPD) . It can be shown that the cross section is dominated by longitudinally polarized photons for $`Q^2\mathrm{}`$; the cross section for transversely polarized photons is suppressed by $`1/Q`$.

The only application of the GPD approach for the reaction has been performed by Mankiewicz et al. . It turns out however that the cross section comes out to large by order of magnitude. This is similar to the two–gluon exchange approach of Brodsky et al.  which can be understood as a LO $`\alpha _s\mathrm{log}(\mathrm{Q}^2/\mathrm{LQCD})`$ calculation, as shown by Frankfurt et al. . The latter authors argue that transverse momentum effects in the photon wave function lead to a suppression factor multiplying the leading twist amplitude which leads then to agreement with experiment. Motivated by experience with meson form factors  we include transverse degrees of freedom as well as Sudakov suppressions in the subprocess. By these effects contributions from the end-point regions, in which one of the partons entering the meson wave function becomes soft and where factorization breaks down, are suppressed. A similar idea has been advocated by Vanderhaeghen et al.  for the quark subprocess, however both approaches differ in detail.

The kinematics of the process are depicted in figure 1. The lower blob represents the soft hadron part, which is parametrized111Throughout this paper we will use Ji’s notation for GPDs . by a GPD $`H(x,\xi ,t)`$ in case of unpolarized gluons and $`\stackrel{~}{H}(x,\xi ,t)`$ in case of polarized gluons222The GPDs $`H(x,\xi ,t)`$ and $`\stackrel{~}{H}(x,\xi ,t)`$ are universal functions and may also be accessed e.g. in deeply virtual compton scattering., whereas the upper blob contains the hard scattering amplitude and the soft meson production.

Considering the factorization of the process, we will now discuss the different components of the amplitude of leptoproduction of polarized vector mesons. For the GPD we use a double distribution ansatz  where a normal parton distribution is convoluted with a profile function for which we choose $`h(x,y)=6y(1xy)`$. The meson production is described by a convolution of the distribution amplitude, for which we took the asymptotic form, and the hard scattering part. The hard amplitude was calculated considering only gluonic contributions since the ZEUS data — with which we compared our calculations — are at high enough energies to neglect quark contributions.

Furthermore, according to the modified perturbative approach  we retain the quark transverse degrees of freedom and take into account Sudakov suppression333It has been shown in  that factorization is still valid in the modified perturbative approach.. The inclusion of these effects leads to suppression of the contributions from the soft end–point regions which results in an overall smaller cross section for the $`\gamma _L^{}pV_Lp`$ process.

The final expression for the amplitude in this approach is

$`_{0+,0+}^{V(g)}`$ $`=`$ $`{\displaystyle \frac{2e}{\sqrt{N_c}}}\sqrt{1\xi ^2}(1+\xi )Q{\displaystyle 𝑑x𝑑\tau \tau \overline{\tau }H^g(x,\xi ,t)}`$
$`\times `$ $`{\displaystyle d^{\mathrm{\hspace{0.17em}2}}𝐛\widehat{\mathrm{\Psi }}(\tau ,𝐛)\widehat{T}_g(\tau ,Q,𝐛)\alpha _s(\mu _R)\mathrm{exp}[S(\tau ,b,Q)]},`$
with

$`\widehat{\mathrm{\Psi }}(\tau ,𝐛)`$ $`=`$ $`2\pi {\displaystyle \frac{f_V}{\sqrt{2N_c}}}\mathrm{\hspace{0.17em}6}\tau \overline{\tau }\mathrm{exp}\left[\tau \overline{\tau }{\displaystyle \frac{𝐛^2}{4a_V^2}}\right],`$
$`S(\tau ,b,Q)`$ $`=`$ $`s(\tau ,b,Q)+s(\overline{\tau },b,Q){\displaystyle \frac{4}{\beta _0}}\mathrm{ln}{\displaystyle \frac{\mathrm{ln}\left(\mu _R/\mathrm{\Lambda }_{QCD}\right)}{\mathrm{ln}\left(1/b\mathrm{\Lambda }_{QCD}\right)}},`$
$`s(\beta ,b,Q)`$ $`=`$ $`{\displaystyle \frac{8}{3\beta _0}}\left(\widehat{q}\mathrm{ln}\left({\displaystyle \frac{\widehat{q}}{\widehat{b}}}\right)\widehat{q}+\widehat{b}\right)+\mathrm{NLL}\mathrm{terms}\text{[12]},`$
$`\widehat{q}=\mathrm{ln}\left(\beta Q/(\sqrt{2}\mathrm{\Lambda }_{\mathrm{QCD}})\right),\widehat{b}=\mathrm{ln}\left(1/(b\mathrm{\Lambda }_{\mathrm{QCD}})\right)`$
and the hard scattering amplitude is denoted by $`\widehat{T}_g(\tau ,Q,𝐛)`$.

As can be seen in figure 2, our calculation describes experimental data quite well. Without the inclusion of transverse quark momenta and sudakov suppression the curves would be approximately an order of magnitude above the data. The modified perturbative approach therefore is needed to correctly reproduce the data.

Due to the uncertainties, the cross section of vector meson photoproduction ist not well suited to obtain informations on the GPDs, in this context the measurement of asymmetries seems more promising. Especially the measurement of the double–spin asymmetry in vector meson production — a task, for which newer experiments like HERMES are well suited — may provide interesting insights concerning the polarized gluon GPDs.

In the future it is planned to calculate all $`\gamma _{}^{}V_{}`$ amplitudes — $`\gamma _{}^{}V_L`$ is already available — to examine spin–dependent effects by e.g. calculating the spin density matrix elements of the $`\rho `$–meson. In this context it is worthwhile to note that the mechanism which we propose for the suppression of the cross section also serves to regularize the infrared divergences, which appear in the calculation of the $`\gamma _{}^{}V_L`$ amplitudes.

The interest in the $`\gamma _{}^{}V_L`$ transitions bases not only in possible measurements of spin-dependent observables, e.g. at COMPASS, but some information on these amplitudes is already available from the H1  and ZEUS  measurements of the decay density matrix elements in leptoproduction of $`\rho `$-mesons. The analysis of these matrix elements demonstrates that both the $`\gamma _{}^{}V_{}`$ and the $`\gamma _{}^{}V_L`$ transition amplitudes are non-negligible even at photon virtualities as large as $`10\text{GeV}^2`$. It is common to analyse the density matrix elements under the assumption of small relative phases between the amplitudes as is characteristic of Pomeron-type models . In this case the $`\gamma _{}^{}V_L`$ transition amplitudes amount to about $`20\%`$ of the dominant $`\gamma _L^{}V_L`$ amplitude  and, hence, s-channel helicity conservation seems to hold with a fairly high precision.

Also of interest is the application of our method to the calculation of the $`\mathrm{\Phi }`$ and $`J/\mathrm{\Psi }`$ photoproduction cross sections which we intend to do in the near future.


text:
# Constrained Curve Fitting

## 1 INTRODUCTION

With recent developments in lattice QCD, a wide range of high-precision calculations is now possible. To achieve high precision we need tight control of all the systematic errors inherent in the analysis of simulation data. A particularly irritating source of such systematic errors has been curve fitting, which plays a central role in almost all applications of lattice simulations.

Lattice theorists use curve fitting in many different ways: we use it, for example, to extract hadronic masses and matrix elements from Monte Carlo estimates of correlators, to extrapolate light-quark masses to physical values, and to fit short-distance Monte Carlo results to perturbative expansions. In each such case the theory we fit to Monte Carlo data has an infinite number of parameters. In order to fit such a theory using a finite amount of data, we normally truncate both the data set and the theory. Thus we might retain only data for the smallest quark masses, and extrapolate linearly to physical masses. Or we might fit only the large-$`t`$ behavior of a hadronic correlator, ignoring contributions from all but the lowest-mass hadron. Such truncations are usually necessary to obtain good fits with reasonable error estimates for the fit parameters, but they necessarily increase both the statistical and systematic uncertainties in the results. Statistical uncertainties are increased because Monte Carlo data is discarded. Systematic uncertainties are introduced by omitting parts of the theory that could conceivably be significant. Such systematic effects have proven particularly difficult to estimate in past analyses, and these are the principle focus of this article.

In this paper we show how to circumvent the shortcomings of the traditional approach by using constrained curve fitting, a simple modification of standard maximum likelihood techniques . Constrained curve fits, while widely used in other research areas, are not used much by the lattice QCD community. Here we will show how they can be adapted for use in lattice calculations, and we will document their striking advantages over traditional methods. Constrained curve fits provide an elegant procedure for incorporating systematic uncertainties due to underconstrained parts of a theory (high-energy states, for example). Furthermore they allow us to fit much more data — for example, correlators down to $`t=0`$. And they are numerically more robust than unconstrained fits. Finally the formalism for constrained curve fitting can be used to estimate the lattice spacings and quark masses that will minimize the errors in a large-scale simulation.

## 2 CONSTRAINED FITS

### 2.1 The Problem

The central problem is illustrated by the analysis of a meson correlator. Simulations generate a Monte Carlo estimate, $`\overline{G(t)}`$, of the correlator for a finite number of time steps, say $`t=0,1\mathrm{}23`$. Theory tells us that the exact correlator has the form

$$G_{\mathrm{th}}(t;A_n,E_n)=\underset{n=1}{\overset{\mathrm{}}{}}A_n\mathrm{e}^{E_nt},$$
(1)
where we assume that the energies $`E_n`$ are in order of increasing size. The challenge is to fit an infinite number of amplitudes $`A_n`$ and energies $`E_n`$ using only 24 $`\overline{G(t)}`$’s.

Traditional fits minimize $`\chi ^2(A_n,E_n)`$ by varying $`A_n`$ and $`E_n`$, where

$$\chi ^2(A_n,E_n)\underset{t,t^{}}{}\mathrm{\Delta }G(t)\sigma _{t,t^{}}^2\mathrm{\Delta }G(t^{}),$$
(2)
and

$$\mathrm{\Delta }G(t)\overline{G(t)}G_{\mathrm{th}}(t;A_n,E_n).$$
(3)
The correlation matrix is estimated from the Monte Carlo:

$$\sigma _{t,t^{}}^2\overline{G(t)G(t^{})}\overline{G(t)}\overline{G(t^{})}.$$
(4)
Unfortunately this fitting procedure is singular here since there are more fit parameters, $`A_n`$ and $`E_n`$, than data; the final uncertainties in the fit parameters are infinite. Additional information is needed if we are to proceed.

The information we normally add is that the $`A_n`$’s are well behaved, and therefore contributions from high-energy states are suppressed at large $`t`$ by the exponentials in the correlator, Eq. (1). Thus there is a $`t_{\mathrm{min}}`$ above which only the first one or two terms in $`G_{\mathrm{th}}`$ make statistically significant contributions. The standard procedure therefore is to retain, say, only the first two terms in $`G_{\mathrm{th}}`$ and to fit them using only Monte Carlo results from $`tt_{\mathrm{min}}`$. The trick is to find the best $`t_{\mathrm{min}}`$. Choosing $`t_{\mathrm{min}}`$ too small biases the $`E_n`$’s and $`A_n`$’s away from their true values, introducing systematic errors (because two terms is not enough in $`G_{\mathrm{th}}`$). Choosing $`t_{\mathrm{min}}`$ too large gives statistical errors $`\sigma _{E_n}`$ and $`\sigma _{A_n}`$ that are too large, since useful data is discarded. One typically tries to increase $`t_{\mathrm{min}}`$ until the statistical errors mask any possible systematic error. Without a reliable quantitative estimate of the systematic error, however, any procedure for setting $`t_{\mathrm{min}}`$ is necessarily ad hoc.

To illustrate the dependence on $`t_{\mathrm{min}}`$, we plot results for $`E_1`$ and $`E_2`$ from 2-term fits with various $`t_{\mathrm{min}}`$’s in Fig. 1. The Monte Carlo data for these fits was obtained by averaging 840 $`\mathrm{{\rm Y}}`$ correlators evaluated at (quenched) $`\beta =6`$ with local sources and sinks . The competition between large systematic errors for small $`t_{\mathrm{min}}`$ and large statistical errors for large $`t_{\mathrm{min}}`$ is particularly apparent for $`E_2`$ in this plot.

### 2.2 A Solution

Our goal should be to fit all the Monte Carlo data ($`t_{\mathrm{min}}=0`$) using as many terms as we wish in $`G_{\mathrm{th}}`$. As we add more terms to $`G_{\mathrm{th}}`$, however, the errors on the leading parameters grow steadily in a traditional analysis, as is evident in Fig. 2. The reason is easily understood. The large uncertainties in $`E_1`$ and $`E_2`$ for the 8-term fit, for example, result because the parameters for higher-energy states are poorly constrained by the data and therefore wander off to unphysical values. Thus amplitude $`A_4`$ ranges between five and ten times $`A_1`$ in the 8-term fit, while quark models suggest that $`A_4`$ is of order $`A_1`$ or smaller. Since the allowed range for $`A_4`$ affects the error estimates for other parameters, the errors on $`E_1`$ and $`E_2`$ will be unreasonably large so long as the fitting code assumes that $`A_410A_1`$ is plausible. We need some way to teach physics to the fitting code.

To constrain fit parameters to physically reasonable ranges, we augment the $`\chi ^2`$ before minimizing:

$$\chi ^2\chi _{\mathrm{aug}}^2\chi ^2+\chi _{\mathrm{prior}}^2,$$
(5)
where

$$\chi _{\mathrm{prior}}^2\underset{n}{}\frac{(A_n\stackrel{~}{A}_n)^2}{\stackrel{~}{\sigma }_{A_n}^2}+\underset{n}{}\frac{(E_n\stackrel{~}{E}_n)^2}{\stackrel{~}{\sigma }_{E_n}^2}.$$
(6)
The extra terms in $`\chi _{\mathrm{aug}}^2`$ favor $`A_n`$’s in the interval $`\stackrel{~}{A}_n\pm \stackrel{~}{\sigma }_{A_n}`$ and $`E_n`$’s in $`\stackrel{~}{E}_n\pm \stackrel{~}{\sigma }_{E_n}`$. The $`\stackrel{~}{A}_n`$’s, $`\stackrel{~}{\sigma }_{A_n}`$’s … are inputs to the fitting procedure. We choose reasonable values for them on the basis of prior knowledge. This set of input parameters is referred to collectively as the “priors.”

Having chosen the priors, the procedure for a constrained fit is to minimize $`\chi _{\mathrm{aug}}^2`$ fitting all of the Monte Carlo data ($`t_{\mathrm{min}}=0`$). The number of terms in $`G_{\mathrm{th}}`$ is increased until fit results converge for the parameters of interest. Unlike $`t_{\mathrm{min}}`$, the number of terms in $`G_{\mathrm{th}}`$ need not be optimized; it is simply increased until results converge. This is illustrated by fit results for $`E_1`$ and $`E_2`$ from our $`\mathrm{{\rm Y}}`$ data, which are plotted in Fig. 3 for fits with different numbers of terms.

The numerics are greatly improved by the constraints. For example, one can easily fit 100 terms in $`G_{\mathrm{th}}`$ to the $`\mathrm{{\rm Y}}`$ data, even though there are only 24 data points. The fit results for all but the first few parameters simply reproduce the prior information in such a highly overparameterized fit.

The error estimates for the fit parameters in our $`\mathrm{{\rm Y}}`$ fits automatically combine both the statistical errors in the Monte Carlo data, and the systematic errors due to our limited knowledge, as specified by the priors, about the poorly constrained high-energy states. To see how much error is due to each source we refit with all $`\stackrel{~}{\sigma }`$’s doubled. Results for the energies change as follows:

$$\begin{array}{ccc}E_1=0.4526(15)\hfill & & 0.4528(14)\hfill \\ E_2=0.683(49)\hfill & & 0.697(50)\hfill \\ E_3=1.05(12)\hfill & & 1.10(21).\hfill \end{array}$$
(7)
Evidently $`E_1`$ and $`E_2`$ are determined largely by the Monte Carlo data, while $`E_3`$ is strongly affected by the priors. The insensitivity of the leading parameters to the priors is typical for high-quality data. The result for $`E_2`$ is impressive given that it comes from a single local-local correlator. (It also agrees with results obtained from multi-source/sink fits .)

We actually parameterized our $`\mathrm{{\rm Y}}`$ fits in terms of parameters $`a_n\mathrm{log}A_n`$ and $`\epsilon _n\mathrm{log}(E_nE_{n1})`$. This parameterization builds in a priori requirements, $`A_n>0`$ and $`E_n>E_{n1}`$, which improve the fits. Using previous simulations as a guide, we chose priors that favored $`a_n\mathrm{log}0.02\pm \mathrm{log}2`$ and $`\epsilon _n\mathrm{log}0.2\pm \mathrm{log}2`$ or:

$$A_n0.02\begin{array}{c}+0.02\\ 0.01\end{array}E_nE_{n1}0.2\begin{array}{c}+0.2\\ 0.1\end{array}.$$
(8)
Our best 5-term fit to the $`\mathrm{{\rm Y}}`$ data is shown in Fig. 4, together with the data. The fit is excellent all the way down to $`t=0`$: the minimum $`\chi _{\mathrm{aug}}^2`$ divided by the number of data points is 0.8. (This is the correct ratio to examine since $`\chi _{\mathrm{aug}}^2`$ has one extra term beyond those in $`\chi ^2`$ for each fit parameter.) While $`\chi ^2`$, $`Q`$, etc. play a crucial role in traditional fits, where $`t_{\mathrm{min}}`$ must be optimized, they are of secondary importance in constrained curve fitting. The key criterion is convergence as the number of parameters is increased. If $`\chi _{\mathrm{aug}}^2`$ per data point is significantly larger than 1 after convergence, then there is likely a mistake in either the data or the theory.

We have experimented with constrained fits for a wide variety of other correlators, including fits of 30–40 parameters for $`4\times 3`$ matrix $`G`$’s, simultaneous fits of multiple channels (e.g., $`\pi `$ and $`\rho `$), static potentials and glueball masses, and correlators involving staggered quarks. All work well. In some cases fits are greatly simplified. An example is the fit shown in Fig. 5 for a $`B`$ meson correlator made from an NRQCD propagator for the $`b`$ quark and a staggered quark propagator for the $`d`$ quark . Such fits are complicated by contributions from opposite parity states, introduced by the staggered quarks, that oscillate with an overall factor $`(1)^t`$. Traditional fits have difficulty quantifying the importance of these contributions since they are small at large $`t`$. Constrained fitting down to $`t=0`$, however, makes it easy to account and correct for them.

## 3 THEORY AND INTERPRETATION

### 3.1 Bayesian Statistics

Bayesian statistics provides a useful framework for understanding the assumptions that go into constrained curve fitting . In this approach the analysis is recast in terms of probabilities: $`P(A|B)`$ denotes the probability that $`A`$ is true or correct assuming $`B`$ is true. For compactness we denote the set of all parameters by $`\rho =\{A_1,E_1\mathrm{}\}`$, the set of prior parameters by $`\eta =\{\stackrel{~}{A}_1,\stackrel{~}{\sigma }_{A_1}\mathrm{}\}`$, and the Monte Carlo data by $`\overline{G}`$.

We begin with two assumptions. The first is that the Monte Carlo data set is sufficiently large that $`\overline{G}`$ has Gaussian statistics (Central Limit Theorem). Then the probability density for obtaining a particular $`\overline{G}`$ given a particular theory, specified by $`\rho `$, is

$$P(\overline{G}|\rho )\mathrm{e}^{\chi ^2(\rho )/2}$$
(9)
where $`\chi ^2`$ is defined as in Eq. (2).111This is not strictly correct since the correlation matrix $`\sigma _{t,t^{}}^2`$ in $`\chi ^2`$ will in general depend upon $`\rho `$, while in practice we use a fixed Monte Carlo estimate of it.

What we need ultimately is the probability that a particular set of parameters $`\rho `$ is correct given the Monte Carlo data — that is, we need $`P(\rho |\overline{G})`$, not $`P(\overline{G}|\rho )`$. $`P(\rho |\overline{G})`$ is connected to $`P(\overline{G}|\rho )`$ by Bayes Theorem:222This formula follows from the trivial identity for probabilities: $`P(\rho \overline{G})=P(\overline{G}|\rho )P(\rho )=P(\rho |\overline{G})P(\overline{G})`$.

$$P(\rho |\overline{G})=\frac{P(\overline{G}|\rho )P(\rho )}{P(\overline{G})}P(\overline{G}|\rho )P(\rho ).$$
(10)
Here $`P(\overline{G})`$ is the probability of obtaining a particular $`\overline{G}`$ from any theory; it is $`\rho `$ independent. More important is $`P(\rho )`$ which is the probability that a particular set of parameters $`\rho `$ is correct in the absence of any new data. It contains what we know about the parameters before we begin the fit. This is called the “prior” distribution, or simply the prior.

Our second assumption is that the prior distribution can be approximated by the Gaussian

$$P(\rho )=\mathrm{e}^{\chi _{\mathrm{prior}}^2(\rho )/2}$$
(11)
where $`\chi _{\mathrm{prior}}^2`$ is defined as in Eq. (6). The a priori assumption therefore is that $`\rho _i\stackrel{~}{\rho _i}\pm \stackrel{~}{\sigma }_{\rho _i}`$. With this assumption our final probability function is

$$P(\rho |\overline{G})\mathrm{e}^{\chi _{\mathrm{aug}}^2(\rho )/2}$$
(12)
where $`\chi _{\mathrm{aug}}^2`$ is the augmented $`\chi ^2`$ introduced in Section 2.2.

The choice of a Gaussian for the prior distribution is arbitrary; other choices might well be appropriate. There is, however, an argument that suggests Gaussians. Beyond specifying an average value for each parameter and a standard deviation, we want the prior distribution to be as unbiased as possible. In general the least biased choice for a probability density is the one that minimizes the information content, or, equivalently, maximizes the entropy,

$$SP(\rho )\mathrm{log}P(\rho )𝑑\rho .$$
(13)
Given the constraints $`\rho =\stackrel{~}{\rho }`$ and $`\rho ^2\rho ^2=\stackrel{~}{\sigma }^2`$, a simple variational calculation shows that the entropy is maximized by the Gaussian

$$P(\rho )=\frac{1}{\sqrt{2\pi }\stackrel{~}{\sigma }}\mathrm{e}^{(\rho \stackrel{~}{\rho })^2/2\stackrel{~}{\sigma }^2}.$$
(14)
This argument is less compelling than it seems, however, because it relies upon the implicit (and arbitrary) assumption that all $`\rho `$’s are equally likely in the absence of any information about the mean or standard deviation. A different parameterization implies a different assumption.

### 3.2 Error Estimation

Given $`P(\rho |\overline{G})`$ we can compute everything we want to know using integrals; in principle no minimization is needed. For example, we obtain a statistical estimate of an arbitrary function of the parameters using

$$f(\rho )=B^1\mathrm{e}^{\chi _{\mathrm{aug}}^2(\rho )/2}f(\rho )d^n\rho $$
(15)
where

$$B\mathrm{e}^{\chi _{\mathrm{aug}}^2(\rho )/2}d^n\rho ,$$
(16)
and the variance is $`\sigma _f^2f^2f^2`$, as usual. In practice these integrals are quite difficult to evaluate for all but the simplest of fits. This is because $`P(\rho |\overline{G})`$ is typically very sharply peaked about its maximum. For smaller problems, adaptive Monte Carlo integrators, such as `vegas`, are effective. For larger problems Monte Carlo simulation techniques, such as the Metropolis or hybrid Monte Carlo methods, can be effective. Still the cost of evaluating the integrals is often prohibitive, particularly when there are lots of poorly constrained parameters (which lead to long, narrow, high ridges in the probability distribution). Consequently efficient approximations are useful.

One approximation assumes that enough statistics have been accumulated so that $`\chi _{\mathrm{aug}}^2(\rho )`$ is approximately quadratic in the $`\rho `$’s throughout the dominant integration region:

$$\chi _{\mathrm{aug}}^2(\rho )\chi _{\mathrm{aug},\mathrm{min}}^2+(\rho \rho ^{})_i(\rho \rho ^{})_jC_{ij}^1$$
(17)
where parameters $`\rho =\rho ^{}`$ minimize $`\chi _{\mathrm{aug}}^2`$. In this quadratic limit we can approximate

$$f(\rho )f(\rho ^{}),\sigma _f^2C_{ij}_if(\rho ^{})_jf(\rho ^{})$$
(18)
provided $`f(\rho )`$ is sufficiently smooth as well. This is the constrained curve fitting procedure outlined in Section 2.2.

A different, more robust approximation to the Bayesian integrals uses a bootstrap analysis, suitably modified to account for the prior distribution . A bootstrap analysis generates an ensemble of $`\rho `$’s whose distribution approximates the Bayesian distribution $`P(\rho |\overline{G})`$. Each $`\rho `$ in this ensemble is obtained by minimizing $`\chi _{\mathrm{aug}}^2`$, but with different Monte Carlo data and different means for the priors for each $`\rho `$. We replace $`\overline{G}`$ in $`\chi _{\mathrm{aug}}^2`$ by the average of a random selection of $`G`$’s, allowing duplicates, from the original Monte Carlo data ensemble, just as in the standard bootstrap method. The new means, $`\stackrel{~}{\rho }_i^{}`$, for the priors are random numbers drawn from a Gaussian distribution with mean $`\stackrel{~}{\rho }_i`$ and standard deviation $`\stackrel{~}{\sigma }_{\rho _i}`$. The new means incorporate the effects of the priors into the bootstrap $`\rho `$ distribution. Given 100 or 1000 $`\rho `$’s generated in this fashion, we then average over the ensemble to compute estimates of any function $`f(\rho )`$. Although this procedure entails a minimization for each $`\rho `$ in the ensemble, we find that is generally much faster than Monte Carlo evaluation of Bayes integrals (Eq. (15)).

The distribution obtained from this modified bootstrap algorithm is not precisely the Bayes distribution $`P(\rho |\overline{G})`$. It has additional factors such as $`\sqrt{detg_{ij}}`$ where

$$g_{ij}\underset{t,t^{}}{}\sigma _{t,t^{}}^2\frac{G(t;\rho )}{\rho _i}\frac{G(t^{};\rho )}{\rho _j}$$
(19)
is a metric induced on $`\rho `$ space . These factors become constants for sufficiently high statistics and so make no difference in that limit. This particular factor is interesting, however, because it makes the measure in $`\rho `$ space invariant under reparameterizations. This suggests that

$$P^{}(\rho |\overline{G})\sqrt{detg_{ij}}\mathrm{e}^{\chi _{\mathrm{aug}}^2/2}$$
(20)
might be a better choice for our Bayesian probability. Then only $`\chi _{\mathrm{prior}}^2`$ would depend upon the parameterization of $`G_{\mathrm{th}}`$, and $`\chi _{\mathrm{prior}}^2`$ usually has little influence on the best-determined fit parameters.

## 4 OTHER APPLICATIONS

One approach to evaluating high-order QCD perturbation theory is to simulate the quantity of interest at high $`\beta `$’s, where QCD is perturbative, and fit perturbative expansions to the results. An example, from , is the $`5\times 5`$ Wilson loop, $`W(5,5)`$, which was measured in Monte Carlo simulations with 9 different lattice spacings that gave couplings $`\alpha _P(3.4/a)`$ covering the range 0.01–0.07. The measured values of $`W(5,5)`$ were fit using an expansion

$$\frac{1}{20}\mathrm{log}W(5,5)=\underset{n=1}{\overset{M}{}}c_n\alpha _P^n(q^{})$$
(21)
where $`q^{}=2.23/a`$ and $`M5`$. Priors $`\stackrel{~}{c}_n=0`$ and $`\stackrel{~}{\sigma }_{c_n}=5`$ were used in the fits to obtain the results in Table 1 (for $`M=5`$). The agreement with exact results is excellent. Refitting with $`c_1`$ and $`c_2`$ set equal to their exact values gives $`c_3=3.91(21)`$, which agrees well with the result in Table 1 but which is much more precise.

The errors for the $`c_n`$’s again include both statistical errors, reflecting the precision of the data, and systematic errors, reflecting the degree of our ignorance of higher-order coefficients. The low-order coefficients ($`n3`$) are largely determined by the Monte Carlo data; their values and errors are essentially unchanged if the $`\stackrel{~}{\sigma }_{c_n}`$’s are doubled, or if $`M`$ is doubled. The high-order coefficients are controlled mostly by the priors. Were we to increase the statistics, more $`c_n`$’s would be determined by the data. Thus the effective order of the fit increases automatically as warranted by the statistics. This is one of the most attractive aspects of constrained fitting: We don’t have to agonize about whether a fit should be linear or quadratic or…; instead we set $`M=10`$ and allow the data to tell us how much it can determine.

Similar issues arise in extrapolations, for example, to the chiral mass limit or to zero lattice spacing. By using constrained fits, we can escape the limitations of linear or quadratic extrapolations. Instead we can include many orders, with constrained coefficients, and allow the Monte Carlo data to determine the relevant order automatically. (See  for an illustration of an extrapolation in lattice size $`L`$.)

## 5 VARIATIONS

### 5.1 Optimizing Simulations

In our numerical fits of perturbation theory to simulation data (Section 4), the final errors for the various perturbative coefficients, $`c_n`$, are affected by the $`\alpha _P`$’s at which we chose to simulate. Our fitting formalism can be used to determine the choice of $`\alpha _P`$’s that minimizes these errors, thereby optimizing the simulation . To illustrate how this is done, we consider a simpler problem. We measure a quantity $`w`$ whose perturbative expansion is

$$w=1+c_1\alpha _P+c_2\alpha _P^2+\mathrm{}.$$
(22)
For simplicity we ignore all terms beyond second order, and we limit ourselves to a single measurement of $`w`$. The challenge is to determine which value of $`\alpha _P`$ (and therefore which lattice spacing) will lead to the smallest errors when we determine $`c_1`$ from a Monte Carlo measurement, $`\overline{w}\pm \sigma _w`$, of $`w`$. This error can be determined from the second derivatives of

$$\chi _{\mathrm{aug}}^2\frac{(1+c_1\alpha _P+c_2\alpha _P^2\overline{w})^2}{\sigma _w^2}+\frac{c_1^2+c_2^2}{\stackrel{~}{\sigma }_c^2}$$
(23)
using Eqns. (17) and (18). The last term in this equation is our prior; we expect $`c_n0\pm \stackrel{~}{\sigma }_c`$. We obtain the following formula for the fit error on $`c_1`$:

$$\sigma _{c_1}^2=\frac{(\alpha _P^4\stackrel{~}{\sigma }_c^2+\sigma _w^2)\stackrel{~}{\sigma }_c^2}{(\alpha _P^4+\alpha _P^2)\stackrel{~}{\sigma }_c^2+\sigma _w^2},$$
(24)
where we will assume that $`\sigma _w`$ is independent of $`\alpha _P`$ (which is not true for Wilson loops, for example, unless numerical roundoff errors dominate). Note that $`\sigma _{c_1}\stackrel{~}{\sigma }_c`$ in either of the limits $`\alpha _P0`$ or $`\alpha _P\mathrm{}`$; we would learn nothing new from a simulation in either limit. The error is minimum at

$$\alpha =\alpha ^{}\sqrt{\sigma _w/\stackrel{~}{\sigma }_c},$$
(25)
where $`\sigma _{c_1}/\alpha `$ vanishes. At this value

$$\sigma _{c_1}^2\stackrel{~}{\sigma }_c^2\frac{2\sigma _w}{\stackrel{~}{\sigma }_c+2\sigma _w}.$$
(26)
Thus $`\alpha =\alpha ^{}`$ is the preferred value for a simulation; the optimal value depends upon the amount of computing that is available (to reduce $`\sigma _w`$).

The $`c_1`$ error decreases like $`\sqrt{\sigma _w}`$ as the simulation error, $`\sigma _w`$, vanishes. This is rather slow; we would prefer an error that vanishes like $`\sigma _w`$. At some small value of $`\sigma _w`$ it will be better to divide the available computer resources between simulations at two values of $`\alpha _P`$, rather than just one. The analysis above can be repeated for this case, to determine the two optimal $`\alpha _P`$’s, and the resulting $`\sigma _{c_1}`$ can be compared with that above, to determine when to switch from one $`\alpha _P`$ to two.

Our simple analysis illustrates a generic design step that should precede any large-scale simulation in lattice QCD. By examining $`\chi _{\mathrm{aug}}^2`$’s for important quantities, we can estimate the optimal lattice spacings and quark masses that should be used in the simulation before we simulate. Furthermore we can estimate the optimal fraction of the total computing resource to spend on each case. Simulation costs should be significantly reduced by such optimizations.

### 5.2 Tuning Priors: Empirical Bayes

The denominator in the Bayes expression, Eq. (10), is the probability of obtaining the Monte Carlo average $`\overline{G}`$ given the prior information:

$$P(\overline{G})=P(\overline{G}|\eta )\frac{\mathrm{e}^{\chi _{\mathrm{aug}}^2/2}}{\stackrel{~}{\sigma }_{\rho _i}}d^n\rho ,$$
(27)
which, in the Gaussian limit (Eq. (17)), becomes

$$P(\overline{G}|\eta )\frac{\sqrt{detC_{ij}}}{\stackrel{~}{\sigma }_{\rho _i}}\mathrm{e}^{\chi _{\mathrm{aug},\mathrm{min}}^2/2}.$$
(28)
This probability will be small if the priors are inconsistent with the data. This suggests that one might tune the prior parameters $`\eta _i`$ to maximize probability $`P(\overline{G}|\eta )`$, thereby determining the priors from the data. This is referred to as the “empirical Bayes” method. It can be useful if only one or two prior parameters is being optimized. For example, when fitting the perturbation series in Eq. (21), one might vary all the $`\stackrel{~}{\sigma }_{c_n}`$’s together to find the $`\stackrel{~}{\sigma }_c`$ that maximizes $`P(G_{\mathrm{th}}|\eta )`$. (Note that the Gaussian approximation, Eq. (28), is exact in this case). Then one would fit with all $`\stackrel{~}{\sigma }_{c_n}=\stackrel{~}{\sigma }_c`$. Note, however, that separately optimizing all $`\eta _i`$’s leads to complete nonsense.333The “best” fit is then the $`\rho `$ that minimizes $`\chi ^2`$ (not $`\chi _{\mathrm{aug}}^2`$), and the fit errors are zero!

### 5.3 Nonparameteric Fits and Maximum Entropy

Bayesian methods have previously been used for lattice data, and indeed for fitting correlators . The priors and parameterizations used in these earlier papers, however, are quite different from what we have discussed above. These papers express the correlator in terms of its spectral function, $`\rho (\omega )`$:

$$G(t)𝑑\omega \mathrm{e}^{\omega t}\rho (\omega ).$$
(29)
The $`\omega `$ axis is divided into a large number of increments, say 750 or 1000, with centers $`\omega _j`$, and the integral is approximated by a sum. The fit parameters are the values of $`\rho (\omega _j)`$ for all $`j`$. This is an example of a “nonparametric” fit. Peaks in $`\rho (\omega )`$ signal bound states. Since the number of fit parameters is typically (much) larger than the number data points, a prior distribution is essential.

The Maximum Entropy principle provides a framework for constructing priors for this purpose. The spectral function shares several properties with probability densities: it has random noise, $`\rho (\omega )0`$, and $`\rho (\omega )d\omega `$ is physically meaningful. This suggests that we think of the spectral function as a probability density. The most probable density is the one with the largest entropy, which suggests that we choose a monotonic function of the entropy as a prior: for example,

$$P(\rho )\mathrm{e}^{\alpha S(\rho )}$$
(30)
is maximum when the entropy $`S(\rho )`$ is maximum. Using a prior of this form we can fit all of the $`\rho (\omega _j)`$’s to the Monte Carlo data. (Parameter $`\alpha `$ can be optimized using the empirical Bayes method; see the previous section.)

This algorithm has been used to fit the same $`\mathrm{{\rm Y}}`$ data that we used for our constrained fits . The errors that result for the first excited state ($`E_2`$) are 2–3 times larger than we obtained from our constrained fits in Section 2.2. This is not surprising. The maximum entropy prior contains far less information than we included in the priors for the constrained fits, and that is its weakness in this context. One should use all the prior information one has when fitting data. The maximum entropy prior is useful when there is little initial information and no simple parameterization of the theory, as is the case in image processing, for example.

## 6 CONCLUSION

Constrained curve fitting is a major improvement over the traditional fitting procedures used in lattice QCD analyses. Neither the data nor the theory is truncated. Consequently more information, about excited states or higher-order terms, is extracted from the same Monte Carlo simulations. Most importantly these techniques provide a transparent and controlled procedure for including systematic errors due to the infinitely many parameters in the theory that are not constrained by the Monte Carlo data. Constrained curve fitting is useful for any kind of analysis that involves correlators. It is also useful for extrapolating to small quark masses or lattice spacings; the order of extrapolation increases automatically to match the statistical precision of the data. Finally we can use the formalism to estimate the quark masses and lattice spacings that optimize a large-scale simulation given a finite computing resource.

Priors make some people uneasy because the constraints are added information, beyond the Monte Carlo data, that can bias the fit results. Such people might argue that one should “let the data speak for itself.” However, it is precisely our ability to fold extra information into the fit that makes constrained curve fitting so much more powerful and useful than unconstrained fitting. Better a priori knowledge, for example from earlier analyses, results in smaller $`\stackrel{~}{\sigma }`$’s and smaller systematic errors. No a priori knowledge, which corresponds to the unconstrained case, (correctly) results in infinite systematic errors. Traditional fitting is a special case of constrained fitting, with infinitely broad priors for the first few parameters and infinitely narrow priors, centered around zero, for the others. The priors we use in our constrained fits are far more reasonable.

## ACKNOWLEDGMENTS

We thank Persis Drell, Saul Teukolsky, and, especially, Tom Loredo for introducing us to Bayesian techniques. We thank Junko Shigemitsu and Matthew Wingate for the use of their Monte Carlo data and for many useful discussions. This work was supported in part by the National Science Foundation, and by the Department of Energy.