Approximate Actions for Lattice QCD Simulation

TL;DR

This paper introduces a systematic method for creating approximate actions in lattice QCD simulations, demonstrating their effectiveness in simpler models and exploring their potential for high-statistics studies with minimal full simulation tuning.

Contribution

It develops a set of tuning conditions for approximate actions and evaluates their performance in lattice models, providing a framework for efficient QCD simulations.

Findings

Approximate actions can match true actions under certain conditions.

Delicate tuning is required for accurate physical predictions.

The approach is promising for high-statistics studies after initial full simulations.

Abstract

We describe a systematic approach to generating approximate actions for the lattice simulation of QCD. Three different tuning conditions are defined to match approximate with true actions, and it is shown that these three conditions become equivalent when the approximate and true actions are sufficiently close. We present a detailed study of approximate actions in the lattice Schwinger model together with an exploratory study of full QCD at unphysical parameter values. We find that the technicalities of the approximate action approach work quite well. However, very delicate tuning is necessary to find an approximate action which gives good predictions for all physical observables. Our best view of the immediate applicability of the methods we describe is to allow high statistics studies of particular physical observables after a low statistics full fermion simulation has been used to prepare the stage.

Tables

Table 1. TABLE I.: Sequence of approximations to T 𝑇 T used to test the effects of adding loops of increasing size to generate an approximate Schwinger model action. The loops included in these approximations are defined in Figure 3 . The superscripts ( n ) 𝑛 (n) on the coefficients a i ( n ) superscript subscript 𝑎 𝑖 𝑛 a_{i}^{(n)} indicate the approximation used to calculate those coefficients. Note that, in all cases, approximation k 𝑘 k uses coefficients determined by simulating with approximation k − 1 𝑘 1 k-1 .

Approximation	Form of $T$ used
0	$T=0$
1	$T=a_1^{\left(0\right)}{\widehat{𝒲}}_1$
2	$T=a_1^{\left(1\right)}{\widehat{𝒲}}_1+a_2^{\left(1\right)}{\widehat{𝒲}}_2$
3	$T=a_1^{\left(2\right)}{\widehat{𝒲}}_1+a_2^{\left(2\right)}{\widehat{𝒲}}_2+a_3^{\left(2\right)}{\widehat{𝒲}}_3+a_{3^{\prime }}^{\left(2\right)}{\widehat{𝒲}}_{3^{\prime }}$

Table 2. TABLE II.: Numerical values for the coefficients a i subscript 𝑎 𝑖 a_{i} in the expansion of T = ∑ i a i ​ 𝒲 ^ i 𝑇 subscript 𝑖 subscript 𝑎 𝑖 subscript ^ 𝒲 𝑖 T=\sum_{i}a_{i}\hat{\cal W}_{i} and their relative contributions (see text).

	$a_i^{\left(0\right)}$	$a_i^{\left(1\right)}$	$a_i^{\left(2\right)}$	$a_i^{\left(3\right)}$	${\sigma }^{\left(3\right)}\left({\widehat{𝒲}}_i\right)$	$a_i^{\left(3\right)}{\sigma }_t^{\left(3\right)}\left({\widehat{𝒲}}_i\right)$
$\widehat{{𝒲}_1}$	$0.212\left(10\right)$	$0.208\left(10\right)$	$0.197\left(11\right)$	$0.222\left(10\right)$	$10.2\left(2\right)$	$2.26\left(11\right)$
$\widehat{{𝒲}_2}$	$0.085\left(7\right)$	$0.083\left(7\right)$	$0.073\left(7\right)$	$0.083\left(7\right)$	$13.8\left(3\right)$	$1.15\left(10\right)$
$\widehat{{𝒲}_3}$	$0.041\left(6\right)$	$0.018\left(11\right)$	$0.034\left(7\right)$	$0.048\left(8\right)$	$12.7\left(2\right)$	$0.61\left(10\right)$
$\widehat{{𝒲}_{3^{\prime }}}$	$0.025\left(4\right)$	$0.026\left(5\right)$	$0.034\left(5\right)$	$0.029\left(4\right)$	$23.9\left(4\right)$	$0.69\left(10\right)$

Table 3. TABLE III.: Rectangular Wilson loop measurements made on configurations generated by full dynamical fermion simulation (true), by an approximate action with 1 × 1 1 1 1\times 1 loop added (column 2), with 1 × 1 1 1 1\times 1 and 3 × 3 3 3 3\times 3 loops added (column 3), and with 1 × 1 1 1 1\times 1 , 3 × 3 3 3 3\times 3 , and 5 × 5 5 5 5\times 5 loops added (column 4). In a given column, expectations of loops included in the corresponding approximate action are marked with an asterisk.

	true	$1\times 1$	$+3\times 3$	$+5\times 5$
$1\times 1$	805.1(3)	805.5(2)*	804.9(2)*	804.4(2)*
$2\times 2$	415.7(9)	391.8(7)	410.4(8)	408.2(8)
$3\times 3$	158.9(11)	118.9(9)	162.4(10)*	160.6(10)*
$4\times 4$	51.2(10)	23.7(7)	47.0(8)	46.7(9)
$5\times 5$	15.1(8)	3.8(7)	8.9(7)	13.8(6)*
$1\times 2$	640.2(6)	633.5(4)	637.6(4)	636.4(4)
$1\times 3$	512.2(8)	498.2(5)	507.5(5)	506.3(6)
$1\times 4$	411.2(10)	392.0(6)	403.9(6)	402.9(7)
$1\times 5$	330.9(11)	308.5(7)	321.5(7)	320.8(8)
$1\times 8$	172.9(11)	150.4(7)	162.0(8)	162.3(9)
$2\times 6$	88.8(9)	58.4(6)	77.9(7)	78.0(8)
$3\times 8$	14.4(6)	3.7(5)	10.0(5)	10.5(6)

Table 4. TABLE IV.: Comparison of Wilson loop expectations calculated with a true full dynamical simulation, and with approximate actions containing 1 × 1 1 1 1\times 1 , 3 × 3 3 3 3\times 3 and 5 × 5 5 5 5\times 5 loops. Column 2 lists the loop expectations determined with approximate action coefficients tuned on a quenched ensemble. The approximate action coefficients in column 4 are tuned on an ensemble generated from the approximate action of column 2, and represent therefore a retuning of the coefficients generating column 2. Columns 3 and 5 list the corresponding deficit corrected loop expectations for each approximate action.

loop	$trueaction$	$quenchedtuning$	$quenchedtuning$	$retuned$	$retuned$
			$\left(corrected\right)$		$\left(corrected\right)$
$1\times 1$	$201.5\left(02\right)$	$201.2\left(02\right)$	$201.2\left(02\right)$	$201.5\left(01\right)$	$201.5\left(01\right)$
$2\times 2$	$105.2\left(06\right)$	$101.8\left(08\right)$	$104.4\left(11\right)$	$103.3\left(05\right)$	$104.3\left(08\right)$
$3\times 3$	$40.5\left(07\right)$	$39.8\left(07\right)$	$39.8\left(07\right)$	$42.6\left(08\right)$	$42.6\left(08\right)$
$4\times 4$	$12.9\left(06\right)$	$12.4\left(05\right)$	$13.1\left(10\right)$	$13.7\left(03\right)$	$14.1\left(08\right)$
$5\times 5$	$3.6\left(04\right)$	$3.7\left(04\right)$	$3.7\left(04\right)$	$5.1\left(04\right)$	$5.1\left(04\right)$
$1\times 2$	$160.7\left(04\right)$	$159.2\left(04\right)$	$160.2\left(06\right)$	$159.8\left(02\right)$	$160.2\left(03\right)$
$1\times 3$	$128.9\left(05\right)$	$126.5\left(06\right)$	$128.5\left(07\right)$	$127.5\left(02\right)$	$128.5\left(04\right)$
$1\times 4$	$103.4\left(06\right)$	$100.5\left(07\right)$	$103.2\left(09\right)$	$101.6\left(03\right)$	$102.9\left(06\right)$
$1\times 5$	$82.9\left(06\right)$	$79.8\left(08\right)$	$83.2\left(12\right)$	$81.0\left(03\right)$	$82.4\left(08\right)$
$1\times 8$	$42.5\left(07\right)$	$39.5\left(08\right)$	$43.2\left(14\right)$	$41.1\left(05\right)$	$42.9\left(12\right)$
$2\times 6$	$22.2\left(06\right)$	$19.2\left(07\right)$	$21.0\left(13\right)$	$20.5\left(03\right)$	$21.6\left(06\right)$
$3\times 8$	$3.8\left(03\right)$	$2.8\left(04\right)$	$2.0\left(08\right)$	$3.6\left(04\right)$	$3.3\left(08\right)$
$5\times 7$	$1.4\left(03\right)$	$0.8\left(02\right)$	$2.0\left(06\right)$	$0.9\left(02\right)$	$1.2\left(06\right)$
$5\times 8$	$0.4\left(02\right)$	$0.4\left(02\right)$	$0.4\left(05\right)$	$-0.1\left(03\right)$	$1.2\left(05\right)$

Table 5. TABLE V.: The residual variance corresponding to an approximate action incorporating the Wilson loop terms shown. The variance is a measure of the quality of the approximate action.

Terms	Loops	${\sigma }^2\left(T-{\sum }_i{a}_i{𝒲}_i\right)$
0		$2.69\left(27\right)$
1	$3\times 3$	$1.02\left(09\right)$
2	$+1\times 6$	$0.61\left(07\right)$
3	$+5\times 5$	$0.59\left(07\right)$
4	$+1\times 2$	$0.31\left(05\right)$

Table 6. TABLE VI.: Effect of tuning condition choice on the loop coefficients in a tuned approximate action. Loop coefficient values as determined by minimizing σ 2 superscript 𝜎 2 \sigma^{2} are given in the second column. The shifts in these values when higher order tuning conditions are incorporated are given in the remaining columns.

coefficient	${\sigma }^2$	$Operators$	$Acceptance$
	$Minimized$	$Matched$	$Minimized$
$a_{3\times 3}$	$0.2150\left(12\right)$	$-0.013\left(45\right)$	$-0.026\left(94\right)$
$a_{6\times 1}$	$0.0325\left(18\right)$	$-0.001\left(5\right)$	$-0.001\left(10\right)$
$a_{5\times 5}$	$0.0105\left(20\right)$	$-0.004\left(5\right)$	$-0.006\left(8\right)$
$a_{1\times 2}$	$0.0109\left(22\right)$	$-0.003\left(9\right)$	$-0.004\left(18\right)$

Table 7. TABLE VII.: Approximate action optimization data for the lattice QCD test described in the text. Wilson loop operators are defined in column 2. Column 3 lists the value of σ 2 ​ ( T − a i ​ 𝒲 i ) superscript 𝜎 2 𝑇 subscript 𝑎 𝑖 subscript 𝒲 𝑖 \sigma^{2}(T-a_{i}{\cal W}_{i}) (no sum on i 𝑖 i ) for the loop operator defined in that row. Column 4 lists the value of σ 2 ​ ( T − ∑ k = 1 i a i ​ 𝒲 i ) superscript 𝜎 2 𝑇 superscript subscript 𝑘 1 𝑖 subscript 𝑎 𝑖 subscript 𝒲 𝑖 \sigma^{2}(T-\sum_{k=1}^{i}a_{i}{\cal W}_{i}) where the sum includes all the loop operators at that row or higher in the table. The order in which loop operators are listed is the order which minimizes σ 2 ​ ( T − ∑ k = 1 i a i ​ 𝒲 i ) superscript 𝜎 2 𝑇 superscript subscript 𝑘 1 𝑖 subscript 𝑎 𝑖 subscript 𝒲 𝑖 \sigma^{2}(T-\sum_{k=1}^{i}a_{i}{\cal W}_{i}) when at most i 𝑖 i loops are included. Note that we present information only on those loops which generate better approximate actions than the single plaquette loop, which appears at position 35.

Order	Specification	${\sigma }^2\left[i\right]$	${\sigma }^2\left[1+\cdots +i\right]$
1	$\left\{11234\overline{3}\overline{1}\overline{1}\overline{4}\overline{2}\right\}$	$8.7\left(\pm 4.4\right)$	$8.7\left(\pm 4.4\right)$
2	$\left\{111222\overline{1}\overline{1}\overline{1}\overline{2}\overline{2}\overline{2}\right\}$	$82.7\left(\pm 9.2\right)$	$7.2\left(\pm 4.1\right)$
3	$\left\{123413\overline{1}\overline{4}\overline{3}\overline{2}\overline{1}\overline{3}\right\}$	$12.6\left(\pm 4.4\right)$	$6.0\left(\pm 4.0\right)$
4	$\left\{123\overline{1}\overline{2}\overline{3}\right\}$	$14.4\left(\pm 5.0\right)$	$5.8\left(\pm 3.9\right)$
5	$\left\{1213\overline{1}\overline{2}\overline{1}\overline{3}\right\}$	$9.3\left(\pm 4.5\right)$	$5.6\left(\pm 3.9\right)$
6	$\left\{12\overline{1}\overline{2}34\overline{3}\overline{4}\right\}$	$22.5\left(\pm 5.6\right)$	$5.3\left(\pm 3.8\right)$
7	$\left\{1112\overline{1}\overline{1}\overline{1}\overline{2}\right\}$	$19.0\left(\pm 4.6\right)$	$5.1\left(\pm 3.7\right)$
8	$\left\{112\overline{1}\overline{1}\overline{2}334\overline{3}\overline{3}\overline{4}\right\}$	$17.1\left(\pm 4.8\right)$	$5.0\left(\pm 3.8\right)$
9	$\left\{1233\overline{2}\overline{1}\overline{3}\overline{3}\right\}$	$9.5\left(\pm 4.6\right)$	$4.8\left(\pm 3.7\right)$
10	$\left\{111122\overline{1}\overline{1}\overline{1}\overline{1}\overline{2}\overline{2}\right\}$	$60.1\left(\pm 8.4\right)$	$4.6\left(\pm 3.6\right)$
11	$\left\{112334\overline{3}\overline{3}\overline{1}\overline{1}\overline{4}\overline{2}\right\}$	$9.6\left(\pm 4.6\right)$	$4.4\left(\pm 3.7\right)$
12	$\left\{11223\overline{2}\overline{2}\overline{1}\overline{1}\overline{3}\right\}$	$14.4\left(\pm 4.6\right)$	$4.0\left(\pm 3.9\right)$
13	$\left\{112\overline{1}\overline{1}\overline{2}\right\}$	$15.2\left(\pm 4.5\right)$	$3.6\left(\pm 3.7\right)$
14	$\left\{12\overline{1}\overline{2}13\overline{1}\overline{3}\right\}$	$50.1\left(\pm 7.5\right)$	$3.5\left(\pm 3.8\right)$
15	$\left\{12\overline{1}31\overline{2}\overline{1}\overline{3}\right\}$	$20.0\left(\pm 5.2\right)$	$3.2\left(\pm 3.9\right)$
16	$\left\{122\overline{1}\overline{2}\overline{2}13\overline{1}\overline{3}\right\}$	$33.1\left(\pm 5.6\right)$	$3.1\left(\pm 3.8\right)$
17	$\left\{11233\overline{1}\overline{1}\overline{2}\overline{3}\overline{3}\right\}$	$17.7\left(\pm 4.6\right)$	$3.0\left(\pm 3.8\right)$
18	$\left\{112233\overline{2}\overline{2}\overline{1}\overline{1}\overline{3}\overline{3}\right\}$	$23.0\left(\pm 4.7\right)$	$2.7\left(\pm 3.8\right)$
19	$\left\{11233\overline{2}\overline{1}\overline{1}\overline{3}\overline{3}\right\}$	$14.9\left(\pm 4.6\right)$	$2.6\left(\pm 3.7\right)$
20	$\left\{122\overline{1}331\overline{2}\overline{2}\overline{1}\overline{3}\overline{3}\right\}$	$24.3\left(\pm 5.2\right)$	$2.6\left(\pm 3.7\right)$
21	$\left\{123\overline{2}\overline{1}\overline{3}\right\}$	$14.2\left(\pm 5.1\right)$	$2.5\left(\pm 3.6\right)$
22	$\left\{12133\overline{1}\overline{2}\overline{1}\overline{3}\overline{3}\right\}$	$14.5\left(\pm 4.4\right)$	$2.5\left(\pm 3.7\right)$
23	$\left\{1234\overline{3}\overline{1}\overline{4}\overline{2}\right\}$	$12.6\left(\pm 4.8\right)$	$2.4\left(\pm 3.8\right)$
24	$\left\{12\overline{1}\overline{2}12\overline{1}\overline{2}\right\}$	$52.8\left(\pm 7.8\right)$	$2.4\left(\pm 3.8\right)$
25	$\left\{1234\overline{3}\overline{2}\overline{1}\overline{4}\right\}$	$12.6\left(\pm 4.3\right)$	$2.4\left(\pm 3.8\right)$
26	$\left\{12213\overline{1}\overline{2}\overline{2}\overline{1}\overline{3}\right\}$	$10.7\left(\pm 4.4\right)$	$2.3\left(\pm 3.7\right)$
27	$\left\{11122\overline{1}\overline{1}\overline{1}\overline{2}\overline{2}\right\}$	$32.6\left(\pm 5.6\right)$	$2.3\left(\pm 3.8\right)$
28	$\left\{1123\overline{1}\overline{1}\overline{2}\overline{3}\right\}$	$10.5\left(\pm 4.7\right)$	$2.3\left(\pm 3.7\right)$
29	$\left\{112334\overline{3}\overline{3}\overline{2}\overline{1}\overline{1}\overline{4}\right\}$	$12.3\left(\pm 4.5\right)$	$2.3\left(\pm 3.8\right)$
30	$\left\{122\overline{1}31\overline{2}\overline{2}\overline{1}\overline{3}\right\}$	$10.9\left(\pm 4.5\right)$	$2.3\left(\pm 3.8\right)$
31	$\left\{1123\overline{2}\overline{1}\overline{1}\overline{3}\right\}$	$10.9\left(\pm 4.4\right)$	$2.3\left(\pm 4.0\right)$
32	$\left\{11234\overline{3}\overline{2}\overline{1}\overline{1}\overline{4}\right\}$	$10.4\left(\pm 4.4\right)$	$2.3\left(\pm 4.0\right)$
33	$\left\{1122\overline{1}\overline{1}\overline{2}\overline{2}\right\}$	$26.6\left(\pm 4.9\right)$	$2.2\left(\pm 4.0\right)$
34	$\left\{12\overline{1}\overline{2}1\overline{2}\overline{1}2\right\}$	$71.4\left(\pm 9.4\right)$	$2.2\left(\pm 3.9\right)$
35	$\left\{12\overline{1}\overline{2}\right\}$	$22.6\left(\pm 5.6\right)$	$2.2\left(\pm 4.0\right)$
	$T$	$295.0\left(\pm 20.0\right)$

Equations

$$S_t=-\beta \sum _{\square }\frac{1}{3}\mathrm{Re}\mathrm{Tr}{U}_{\square }-\mathrm{Tr}\log M_{\kappa }^{†}{M}_{\kappa }$$

$$S_a=-{\beta }^{\prime }\sum _{\square }\frac{1}{3}\mathrm{Re}\mathrm{Tr}{U}_{\square }.$$

$$⟨F{⟩}_t=\frac{\int [dU]F{e}^{-S_t}}{\int [dU]e^{-S_t}}\text{and}⟨F{⟩}_a=\frac{\int [dU]F{e}^{-S_a}}{\int [dU]e^{-S_a}}.$$

$$⟨F{⟩}_j=\frac{⟨F{e}^{S_i-S_j}{⟩}_i}{⟨e^{S_i-S_j}{⟩}_i}=⟨F{⟩}_i+⟨\left(F-⟨F{⟩}_i\right)C_i({\Delta }_{ij}){⟩}_i,$$

$${\Delta }_{ij}=S_i-S_j=-{\Delta }_{ji}$$

$$C_i(\Delta )=\frac{1}{⟨e^{\Delta }{⟩}_i}{e}^{\Delta }.$$

$$C_i(\Delta )$$

$$⟨F{⟩}_j-⟨F{⟩}_i$$

$$S_a\equiv S_a(a_1,\dots ,a_n).$$

$$(F,G)=⟨F^{\ast }G⟩=\frac{\int [dU]F^{\ast }G{e}^{-S}}{\int [dU]e^{-S}},$$

$$d(F,G)=((F-G),(F-G))=⟨(F-G{)}^{\ast }(F-G)⟩.$$

$$d(S_a-⟨S_a⟩,S_t-⟨S_t⟩)=⟨({\Delta }_{at}-⟨{\Delta }_{at}⟩{)}^2⟩$$

$${\Delta }_{at}=S_a-S_t=-{\Delta }_{ta}$$

$${\sigma }^2({\Delta }_{at})=⟨({\Delta }_{at}-⟨{\Delta }_{at}⟩{)}^2⟩$$

$$⟨({\mathcal{S}}_i-⟨{\mathcal{S}}_i⟩)({\Delta }_{at}-⟨{\Delta }_{at}⟩)⟩=0\text{ for i = 1, …, n,}$$

$${\mathcal{S}}_i=\frac{\partial }{\partial a_i}{S}_a(a_1,\dots ,a_n).$$

$$⟨{\mathcal{S}}_i{⟩}_t=⟨{\mathcal{S}}_i{⟩}_a.$$

$$⟨{\mathcal{S}}_i{⟩}_t=⟨{\mathcal{S}}_i{⟩}_t+⟨({\mathcal{S}}_i-⟨{\mathcal{S}}_i{⟩}_t)C_t(-{\Delta }_{at}){⟩}_t,$$

$$⟨({\mathcal{S}}_i-⟨{\mathcal{S}}_i{⟩}_t)C_t(-{\Delta }_{at}){⟩}_t=0.$$

$${e^{-S_a}|}_U{P}_a(U\to U^{\prime })={e^{-S_a}|}_{U^{\prime }}{P}_a(U^{\prime }\to U),$$

$${e^{{\Delta }_{at}}|}_U{A}_{at}(U\to U^{\prime })={e^{{\Delta }_{at}}|}_{U^{\prime }}{A}_{at}(U^{\prime }\to U)$$

$$A_{at}(U\to U^{\prime })=\min \left(1,e^{{\delta }_{U\to U^{\prime }}({\Delta }_{at})}\right)$$

$${\delta }_{U\to U^{\prime }}({\Delta }_{at})={{\Delta }_{at}|}_{U^{\prime }}-{{\Delta }_{at}|}_U.$$

$$⟨⟨A_{at}⟩⟩=\int d[U^{\prime }]d[U]A_{at}{\frac{e^{-S_t}}{\int d[U]e^{-S_t}}|}_U{P}_a(U\to U^{\prime }).$$

$$e^{⟨⟨{\delta }_{U\to U^{\prime }}({\Delta }_{at})⟩⟩}\le ⟨⟨e^{{\delta }_{U\to U^{\prime }}({\Delta }_{at})}⟩⟩=1$$

$$⟨⟨{\delta }_{U\to U^{\prime }}({\Delta }_{at})⟩⟩\le 0.$$

$$⟨⟨A_{at}⟩⟩\approx erfc\left(\frac{1}{2}\sqrt{-⟨⟨\delta ({\Delta }_{at})⟩⟩}\right).$$

$$P_a(U\to U^{\prime })\to {\frac{e^{-S_a}}{\int d[U]e^{-S_a}}|}_{U^{\prime }}.$$

$$-⟨⟨\delta ({\Delta }_{at})⟩⟩$$

$$=$$

Full text

††preprint: Liverpool Preprint: LTH 365, hep-lat/9608145

I Introduction

The design of improved algorithms for simulation with dynamical fermions has provided an important challenge for lattice QCD. Perturbative arguments and numerical experiments suggest that a significant component of the effect of dynamical quarks can be accounted for by a shift in the effective lattice spacing. The quenched approximation takes this view to its extreme, and ignores all dynamical quark effects except for this shift in scale. In most cases, the systematic errors introduced by quenching have so far turned out to be of the same order as the best presently achievable statistical errors in full dynamical simulations. As available computer power grows, and as we study larger volumes and smaller lattice spacings, we can expect this situation to change, and to begin to see significant discrepancies between quenched and full theory results. Systematic errors due to quenching will then turn out be larger than full theory statistical errors, and the quenched approximation will loose its usefulness.

In a recent paper [1] it has been shown that the quenched algorithm is a single member of a large class of approximate QCD simulation algorithms. The members of this class interpolate in a smooth fashion between the extreme cases of quenched QCD simulations on the one hand, and full dynamical QCD simulations on the other. Individual members can be considered as different approximations to the full dynamical QCD action. So long as the quenched approximation represents the only significant example of an approximate simulation algorithm, it is possible to consider it in isolation and to accept its historically important contributions to lattice QCD while remaining wary of its mathematical correctness. However, once one begins to consider variants which seek to improve on the quenched approximation without themselves being exact, a number of important issues arise, and a more critical approach is needed.

First, of course, is the question of what we mean mathematically by the term ‘approximate action’. The qualitative view of quenched QCD which we normally adopt is that, at fixed lattice spacing and coupling, it produces ensembles which are reasonably close to the ensembles which would be generated for full QCD at some shifted coupling. Reasonable closeness is however a vague term. What it tends to mean in practice is that ensemble expectations of a set of operators in quenched and full QCD seem to agree within statistical errors, so long as one allows for an effective lattice spacing shift.

In this paper we consider the general principles which can be applied to define what is meant by a simulation algorithm based on an approximate action in the sense defined by the quenched approximation. We will describe a systematic approach to generating approximate actions and will discuss conditions which allow the tuning of any candidate approximate action to optimize its matching to a true action of interest. The methods we develop, while motivated by the approximate action view, are also applicable whenever we wish to match results from different lattice simulations. In particular, the tuning we propose defines the optimum matching conditions between Wilson and Kogut-Susskind fermion simulations, and between naive and improved action simulations.

A second issue of concern, when considering simulations with approximate actions, is the question of their computational performance. The results reported in [1] were obtained on very small systems, and for these systems it was found that approximate action simulations were very competitive with exact simulations including full fermions. However, a question left unanswered was whether this competitiveness would continue when the method was applied to more realistic situations. The basic proposal of [1] was to attempt an approximation of the fermion determinant in QCD using Wilson loop operators of increasing size. On a small system there are relatively few such operators, but the number of operators grows very fast as the size of the system increases and the computational work required would also grow. Some preliminary data which we reported in [2] for the case of the Schwinger model seems to bear out this view. There we found that only very slow improvements in the quality of the approximation were achieved as Wilson loops of systematically increasing size were added to the approximate action. An independent study of QCD at realistic lattice sizes [3] was also somewhat negative. This study found that the plaquette operator was actually a very poor approximation to the fermion determinant in typical large simulation. Both of these results suggest that very many Wilson loops would be needed to get a numerically useful approximation in realistic situations.

We report in this paper, however, that the pessimistic view just described is an artifact of the particular approximation prescription defined in [1]. The systematic approach to tuning approximate actions which we develop here allows us to generate much more aggressive approximations than are otherwise possible, and we find that it is not necessary to include Wilson loops of all shapes and sizes to get useful approximate actions. Indeed, the best approximate actions in the Schwinger model and QCD cases which we have studied are found to include only a few Wilson loop operators of quite surprising shape and dimension.

An issue which we do not address in this paper is that of the behaviour of approximate actions in the continuum limit. The quenched case is clearly problematic in this limit [4]. However, the strategy we propose here is a numerical one in which one may choose the most appropriate approximate action for any given simulation at fixed lattice spacing and coupling. If one varies the lattice spacing and coupling, then the corresponding approximate action would also vary. It seems conceivable to us that the variations in approximate action which result would be sufficient to allow a mathematically correct continuum limit of the approximate action to be taken.

The plan of the remainder of this paper is as follows. In Section II, we define a simple cumulant expansion method which allows direct comparison of the expectations of operators calculated in two different path integral measures. We argue that this expansion method is the most efficient way of determining the shifts introduced in expectations of physical observables due to changes in the underlying path integral measure.

The concept of an approximate action is defined in Section III. The quenched QCD action is the prototypical example here. Approximate actions allow us to define corresponding approximate simulation algorithms, and we develop in this section three different ways in which approximate actions might be defined. We then argue that all three definitions become equivalent when the approximate action is sufficiently close to the action which it attempts to approximate.

As an initial test, we consider the Schwinger model. Section IV gives the explicit form of this two-dimensional model which we have used, and describes some necessary technical features of the calculations. Simulation results are presented in the following two sections. In Section V, we consider the quenched approximation in some detail. Physical quantities calculated here include the static potential and the mass of the mesonic bound state of the model. We present comparisons of these observables calculated in the quenched approximate simulation and full dynamical simulation. Section VI describes a variety of different extensions to the quenched approximation. The analysis begins with the scenario proposed in [1] in which the fermion determinant is expanded in terms of Wilson loops of increasing size. We show first that this expansion technique is quite straightforward to execute and appears to be numerically stable. We argue, however, that this ordered expansion is very slow to converge. Instead, we propose an aggressive approach to developing approximate actions which directly targets the interesting physics. This approach is seen to improve considerably the slow convergence apparent in [1].

In Section VII, we describe a trial application of the approximate action approach to QCD. The purpose of the test is to gain some understanding of the potential of the method for QCD simulation and to identify practical problems which remain to be solved. The test described is for a small lattice with heavy quarks.

Finally, our conclusions are presented in Section VIII.

II The Comparison of Path Integral Measures

In all lattice simulation, the goal is to calculate expectation values of interesting operators in some path integral measure. The measure is given by some action which we think of as the true action and denote by $S_{t}$. Our goal is to find an action $S_{a}$ which allows us to approximate the measure generated by $S_{t}$. The prototypical example is given by lattice QCD. The fundamental degrees of freedom here are the $SU(3)$ link matrices $U_{\mu}(x)$ assigned to edges which join nearest neighbouring sites on a four dimensional hypercube. A possible ‘true’ action for this system is the lattice QCD action which includes a Wilson fermion determinant term,

[TABLE]

while a possible approximate action is the quenched lattice QCD action,

[TABLE]

$U_{\Box}$ here represents the product of gauge links around a fundamental plaquette, while $M_{\kappa}$ is the Wilson fermion matrix for hopping parameter $\kappa$. We have also implicitly allowed for a bare coupling shift between exact and approximate cases by including different $\beta$’s in the different actions.

Expectation values of operators $F$ in the path integral measures generated by $S_{t}$ and $S_{a}$ are denoted respectively by $\langle F\rangle_{t}$ and $\langle F\rangle_{a}$, and are given by

[TABLE]

Before addressing the problem of how to choose an approximate action, we first need a technique to compare these expectations. The simplest numerical comparison technique is to generate independent sequences of configurations for each measure, to evaluate expectations on these sequences, and then compare. For our purposes, this simple approach is not computationally feasible. We are attempting to find approximate actions for a given true action because the computational work needed to generate independent configurations in the true measure is very large. Our eventual goal is to find a tuning technique which allows us to generate an approximate action without ever having to simulate the measure of the true action. Thus we are lead to consider methods by which we can use ensembles of configurations generated in one measure to estimate expectations of operators in a second measure.

Provided the fundamental degrees of freedom of two action functionals, $S_{1}$ and $S_{2}$, are the same (*i.e. *provided the configuration space upon which they are defined is the same), we can directly relate expectations in the corresponding measures,

[TABLE]

where

[TABLE]

and

[TABLE]

This comparison formula is exact. However, the calculation of the overlap expectations of $C_{i}(\Delta_{ij})$ with $F-\langle F\rangle_{i}$ is very difficult numerically, since $\Delta_{ij}$ is typically a quantity which is of $O(V)$, where $V$ is the volume of our system. If we use this formula, therefore, to calculate expectations for $S_{j}$ by correcting an ensemble of configurations generated for $S_{i}$, we find that we must re-weight with factors which are $O(e^{V})$. The fluctuations in expectations so calculated are also of $O(e^{V})$, so the work needed to get reasonable statistics is exponential in $V$. As a result, a direct application of this formula to compare $\langle F\rangle_{i}$ and $\langle F\rangle_{j}$ is not feasible numerically.

We adopt instead a cumulant expansion technique which is generated by expanding $C_{i}(\Delta)$ in a Taylor series,

[TABLE]

Substituting this expansion in (4), we find

[TABLE]

The leading term on the right hand side here is the connected part of the correlation of the operator $F$ with the difference in the actions $\Delta_{ij}$. As before, $\Delta_{ij}$ is a quantity of $O(V)$; but now we don’t have to exponentiate, and so we expect fluctuations only of $O(V)$ rather than $O(e^{V})$. Since any numerical simulation must be carried out at finite volume, we can expect that at least this leading correction between the expectation values in the different measures will be calculable with finite work. Of course, higher order terms will require calculation of connected correlations involving higher powers of $\Delta_{ij}$. These correlations will have correspondingly larger fluctuations. The connected correlations involving $\Delta_{ij}^{n}$ will have fluctuations of $O(V^{n})$, for example, and will be correspondingly harder to calculate. One should view this cumulant expansion technique, therefore, as the numerical version of an asymptotic expansion. If the difference $\Delta_{ij}$ is sufficiently small, one or perhaps two terms in this expansion may be calculable with reasonable work.

III Tuning the Approximate Action

In our prototypical example introduced in the last section, the true action $S_{t}$ is the lattice QCD action which combines a single plaquette gauge term and a Wilson fermion determinant term expressed as the trace of the log of the Wilson hopping matrix. Our prototypical approximate action $S_{a}$ includes just the single plaquette gauge term but with a possibly different coupling from that in the true action. More generally, we imagine that the approximate action depends on a finite number of couplings $a_{1},\ldots,a_{n}$,

[TABLE]

One possibility, for example, is that $S_{a}$ includes a sum of different translationally invariant Wilson loop operators [1]. The $a_{i}$ are then the couplings multiplying each different Wilson loop operator in the sum. Another possibility is to replace the fermion determinant term with an approximate bosonic term proposed by Lüscher [5]. The $a_{i}$ then include the plaquette coupling for the gauge term together with all the roots of the polynomial used to approximate the inverse of the fermion matrix in the Lüscher approach. The question we now address is how to choose the couplings $a_{i}$ so that $S_{a}$ becomes a good approximation to the true action $S_{t}$. We define three different tuning schemes for the coefficients, and show that these three different schemes are closely related.

A Minimizing the Distance Norm

The first tuning procedure is described in [1]. There it is observed that the measure defined by an action functional $S$ allows us to define an inner product $(\,,\,)$ on the space of gauge invariant functionals of the link matrices, which takes the form,

[TABLE]

where $F$, $G$, and $S$ are all functions of the link matrices $U_{\mu}(x)$. This inner product turns the space of gauge invariant functions into a Hilbert space, and allows us to define the distance $d(F,G)$ between any two functions as

[TABLE]

Functions of link matrices can now be considered as vectors in this Hilbert space and the natural tuning condition in this language is to minimize the distance, or equivalently to maximize the overlap, between the two vectors which represent the true and approximate actions. Since an overall constant shift in an action functional has no effect on the measure it generates, the tuning should minimize the connected part only. Thus, the first prescription is to minimize

[TABLE]

over the space of couplings $a_{1},\dots,a_{n}$. $\Delta_{at}$ is the difference in the approximate and true actions,

[TABLE]

and the object to be minimized is simply the variance $\sigma^{2}(\Delta_{at})$ of this difference,

[TABLE]

evaluated in the measure $\langle\ \rangle$. Differentiating with respect to each of these couplings in turn, we generate a set of $n$ simultaneous equations to determine the minimum,

[TABLE]

where

[TABLE]

In deriving this particular tuning condition, we have been purposefully vague about the measure in which expectations are to be evaluated. There are two obvious measures available, that generated by $S_{t}$, and that generated by $S_{a}$. Either measure will allow a calculation of the tuned values for the couplings $a_{i}$, so this first tuning condition actually has two possible variants. Our preferred choice is to minimize in the true measure. This minimization is an explicit problem, but in principle requires generating a set of configurations in the full measure. In practice, we could also use the cumulative expansion methods of the last section to estimate expectations in the true measure using our best guess as to the approximate measure. Minimizing in the approximate measure, on the other hand, is an implicit problem since we only learn the correct values for the coefficients $a_{i}$ after the minimization is completed. Minimization in this case requires an iterative method which first guesses the correct couplings, then evaluates the corrections to those guesses.

For notational simplicity in what follows, we will assume that equations of the form (18) apply for all values of $i$, and we will drop the phrase ‘for $i=1,\dots,n$’ in such equations.

B Operator Matching

A different approach to tuning the approximate action is suggested by the observation that, in some cases at least, both true and approximate actions represent different regularizations of the same continuum theory. For the prototypical case of full Wilson fermion action versus quenched action this is, of course, not correct since the quenched approximate action does not have a unitary continuum limit. However, it is possible to imagine cases where this is indeed true. For example, one could consider the case where one action was the full Wilson fermion action, while the second action was an $O(a)$ improved action. In this view, the couplings in both actions are bare parameters in different regularizations of the same continuum theory.

Normally, one fixes the mapping from bare to renormalized couplings by specifying renormalization conditions for each bare parameter whose finite part needs to be fixed. These conditions directly relate regularized theories to the renormalized continuum theories which they regulate. In our case, imposing renormalization conditions which relate the regularized theories to continuum renormalized theories is not feasible, since this requires that we separately generate a full continuum limit extrapolation of each regularized theory. What one can do, however, is to impose renormalization conditions directly between the two lattice regularized theories. If these conditions are maintained in the passage to the continuum limit, then both regularized theories will proceed to the same renormalized continuum theory as we require.

Consider, therefore, matching the approximate action $S_{a}(a_{1},\ldots,a_{n})$ with the true action $S_{t}$. The approximate action has $n$ couplings which must be determined, so we must impose $n$ independent matching conditions between the approximate and true measures to fix these $n$ couplings. The obvious possibilities include matching operators such as $m_{\pi}$, $m_{\rho}$, $m_{N}$ [6]. These particular operators however are not ideal for our purposes since they are difficult to measure accurately on a small number of configurations, and they are not particularly sensitive to small changes in the couplings $a_{i}$ which we need to determine. The choice of matching conditions which we make, however, is not really very important since any set of conditions will correctly fix the cutoff dependent divergences of the couplings being matched, and will differ only in the way they fix finite parts. Different matchings will differ by finite renormalizations only, so we are free to choose some simply implemented matching conditions. An alternative set which is simple to calculate, and which should be sensitive to changes in the $a_{i}$, is to match the operators ${\cal S}_{i}=\partial S/\partial a_{i}$ between the approximate and true actions.

Thus, the second tuning condition we propose is to match the expectations of the $n$ operators ${\cal S}_{i}$ ($i=1,\dots,n$) in the approximate and true measures,

[TABLE]

Applying the cumulant expansion to $\langle S_{i}\rangle_{a}$, we find

[TABLE]

which directly gives

[TABLE]

C Maximizing Acceptance

Consider using an approximate action as part of an exact Markov process which generates configurations distributed according to the measure of the true action. The particular exact algorithm we propose involves a standard Markov transition step which satisfies detailed balance for the approximate action. This step generates new trial configurations which do not have the correct true action distribution. One can then correct to the true distribution by executing a Metropolis accept/reject on the trial configurations. Our third tuning condition for the approximate action is to maximize acceptance in this exact algorithm.

Let $P_{a}(U\rightarrow U^{\prime})$ be a Markov transition probability for link configuration $U$ to go (reversibly) to link configuration $U^{\prime}$. If this transition probability satisfies detailed balance for the approximate action $S_{a}$, then

[TABLE]

where the notation $\left.\right|_{U,U^{\prime}}$ indicates that the relevant action terms are evaluated on the link configurations labeled by $U$ and $U^{\prime}$. To get an exact algorithm for the true action $S_{t}$ we must now execute an accept/reject step. The probability, $A_{at}$, to accept a trial step generated by $P_{a}$ must satisfy a modified detailed balance,

[TABLE]

and the optimum choice for this acceptance probability is

[TABLE]

where

[TABLE]

Maximizing acceptance for this compound algorithm requires maximizing the weighted average of $A_{at}$ over all configurations $U$ which can occur at the start of a trial step, and all configurations $U^{\prime}$ which can occur at the end of such a step. If we denote this weighted average as $\langle\langle A_{at}\rangle\rangle$, we have

[TABLE]

Since $P_{a}$ satisfies detailed balance for action $S_{a}$, it is straightforward to show that

[TABLE]

which implies

[TABLE]

Further, if $\langle\langle\delta(\Delta_{at})\rangle\rangle$ is small, and if only the leading terms in the Taylor expansion of $e^{\langle\langle\delta(\Delta_{at})\rangle\rangle}$ are significant, the acceptance probability is given approximately by [13]

[TABLE]

Acceptance, therefore, is maximized when $-\langle\langle\delta(\Delta_{at})\rangle\rangle$ is minimized. The worst case for acceptance occurs when the transition probability $P_{a}(U\rightarrow U^{\prime})$ defines a perfect heat bath transition for the approximate action,

[TABLE]

Trial configurations $U^{\prime}$ in this case are completely independent of the configurations $U$ which they attempt to replace. In this worst case, we have:

[TABLE]

Maximizing acceptance requires maximizing this last term over the couplings $a_{1},\dots,a_{n}$ in the approximate action $S_{a}$. If we differentiate with respect to the couplings $a_{i}$, we find the following condition:

[TABLE]

An important question to ask once an approximate action has been tuned to maximize acceptance is whether the acceptance achieved by this maximization is large enough to generate a practical exact algorithm. The answer, of course, is given by the value of $\langle\langle A_{at}\rangle\rangle$ which results. $\langle\langle A_{at}\rangle\rangle$ is given in (30). We can rearrange this equation by using the worst case formula for $-\langle\langle\delta(\Delta_{at})\rangle\rangle$ given in (32). When expanded, this latter equation takes the very simple form,

[TABLE]

(The subscript $t$ on $\sigma_{t}^{2}$ denotes that the expectations involved are to be evaluated in the true action measure). Thus, we find

[TABLE]

The actual minimum value for $\langle\langle A_{at}\rangle\rangle$ which defines a practical algorithm is a function of the ratio of work involved in generating uncorrelated trial configurations with approximate and true actions. In QCD, the generation of uncorrelated configurations with, for example, an exact hybrid algorithm is very expensive. This suggests that even very small acceptances might suffice to generate practical exact algorithms based on approximate actions. A 50% acceptance is achieved when $\sigma^{2}_{t}(\Delta_{at})\approx 1$, while a 10% acceptance is achieved when $\sigma^{2}_{t}(\Delta_{at})\approx 5$.

D Tuning Procedure

We have, at this point, proposed three different tuning conditions, (18, 22, 35), to define the adjustable parameters $a_{i}$ in an approximate action. Two of these equations involve the cumulant expansion factor $C_{t}(\pm\Delta_{at})$. As we have argued in the last section, this factor is not directly calculable in a single numerical experiment with a single ensemble generated by either the approximate or true action. Instead, to apply these tuning conditions, we must expand $C_{t}(\pm\Delta_{at})$ in powers of $\Delta_{at}$ using (8). The tuning conditions which result when we keep just the leading and next to leading terms of this expansion take the following forms. First, if vector overlap in the true measure is maximized using (18), we have unchanged

[TABLE]

Second, if operators are matched in approximate and true measures (22), then we have

[TABLE]

Finally, if acceptance is maximized using (35), we have

[TABLE]

These three conditions all agree to leading order in $\Delta_{at}$, and only begin to disagree at next to leading order in the cumulant expansion. If we have managed to find approximate and true actions which agree sufficiently closely that the cumulant expansion is a numerically useful tool for comparing expectations in the two measures, then we would expect the terms at next to leading order in these equations to be much smaller than the leading order terms.

These observations suggests the following self-consistent approach to generating approximate actions for lattice simulation. First, apply one of the three conditions listed here to tune the effective action. This tuning can be carried out for an ensemble of configurations generated either for the approximate or for the full measure, so long as the cumulant expansion is used appropriately. Once the approximate action is determined, a variety of operators should be calculated in this approximate action, and further, the shift between expectations in approximate and true actions should be determined for these operators using the cumulant expansion. If the shift in operators so calculated is small, then we expect that the approximate action provides a good representation of the true action. Further, we expect that all three possible tuning conditions will be approximately satisfied, and that the cumulant expansion will work well for most operators. On the other hand, if the shift in the operators so calculated is large, then the approximate action is not a good representation of the true action. In this case, we find that the different matching conditions will produce significantly different values for the approximate action couplings. The situation which will then prevail is that we will be able to adjust certain operators to match between the approximate and true measures. However, this adjustment will not cause other different operators to match, and simulation with the approximate action will not provide a useful approximate ensemble for the study of properties of the true action ensemble.

At leading order in the cumulant expansion, the common tuning condition from all three methods takes the form

[TABLE]

As written, this condition requires that expectations are to be evaluated in the true measure. As we suggested above, we can use the cumulant expansion to apply this condition using operators evaluated in the approximate measure also, so long as we use the cumulant expansion to correct appropriately. We have

[TABLE]

These two forms differ only by terms which we have already neglected when we work to leading order. Thus we can quite consistently choose to drop the next to leading terms on the right hand side of (43) when we wish to use an approximate action ensemble rather than a true action ensemble to perform tuning.

Our leading order approximate action tuning condition now takes the simple form

[TABLE]

where expectations are evaluated in either true or approximate measures. This condition represents a set of coupled equations for the tunable coefficients $a_{i}$ in the approximate action $S_{a}$, and are precisely the equations which minimize $\sigma^{2}(\Delta_{at})$ over these coefficients. $\sigma^{2}(\Delta_{at})$ therefore provides a key measure of the quality of an approximate action. It also, according to (36), determines whether a practical exact algorithm can be built using a global accept/reject step to correct from approximate to true action. In what follows, we use this variance to discriminate between candidate approximate actions and to guide our search for an exact update algorithm based on the approximate action.

E Wilson Loop Approximations

To close this section, let us consider the technicalities of minimizing $\sigma^{2}$. The approximate actions which we study in the remainder of this paper are simple sums of Wilson loop operators of different sizes and shapes. We will use the notation ${\cal W}_{i}$ to denote a particular Wilson loop operator which we always take to be translationally and rotationally invariant. A typical example of such an operator in QCD is the $1\times 1$ plaquette summed over all sites and orientations,

[TABLE]

Wilson loop operators are always summed over all distinct orientations and locations on the lattice. Thus the operators we consider are rotationally and translationally invariant. This is sufficient for our purposes since our major effort is to find an approximation for the fermion determinant which is also rotationally and translationally invariant.

The class of approximate actions which we will consider now takes the simple form

[TABLE]

and our goal will be to find the best approximation within this class to the true action

[TABLE]

where $T$, defined as,

[TABLE]

is the trace of the log of the fermion determinant for $n_{f}$ equal mass fermions. We have immediately that

[TABLE]

The difference between approximate and true actions for this case is

[TABLE]

where, for convenience, we have absorbed the pure gauge term by redefining the coefficients $a^{\prime}_{i}$,

[TABLE]

and the work of finding an approximate action reduces to that of finding an approximation for $T$ of the form $T\approx\sum_{i}a_{i}{\cal W}_{i}$. For convenience in what follows, we drop the prime superscript on the parameters $a_{i}$ which define this approximation. Our lowest order minimization condition (45) is now simply expressed as

[TABLE]

To determine the parameters $a_{i}$ which define the best approximation to $T$ we must calculate connected correlations of each of the Wilson loop operators included in the approximate action both with themselves (i.e. $\langle({\cal W}_{i}\!-\!\langle{\cal W}_{i}\rangle)({\cal W}_{j}\!-\!\langle{\cal W}_{j}\rangle)\rangle$), and with $T$ (i.e. $\langle({\cal W}_{i}\!-\!\langle{\cal W}_{i}\rangle)(T\!-\!\langle T\rangle)\rangle$). Once these correlators are available, then the evaluation of the $a_{i}$’s is achieved by solving the linear system (56).

It is also useful, in the process of generating the coefficients $a_{i}$ of the approximation to $T$, to define a Gram-Schmidt orthogonalized basis of Wilson loop operators. If the original loops are ${\cal W}_{i}$ for $i=1,\dots,n$, then the orthogonalized loops are $\hat{{\cal W}_{i}}$ and are defined recursively by

[TABLE]

and

[TABLE]

where

[TABLE]

The operators $\hat{\cal W}_{i}$ therefore satisfy

[TABLE]

and, in terms of these orthogonalized operators, the solution of the linear system (56) is given by

[TABLE]

IV The Schwinger Model

In order to test the ideas developed in the proceeding section, we need a model which exhibits the important features of the QCD lattice model, but is computationally accessible. The Schwinger model offers a suitable testing ground. This theory is a two-dimensional gauge theory with gauge group $U(1)$ rather than $SU(3)$. With no fermions present the theory is (almost) trivial. All the interesting dynamics are generated by fermions, and we expect very significant dynamical fermion effects to be present as a result. Computationally, it is possible to generate full dynamical fermion simulations in just a few days on a workstation, so the model allows considerable experimentation with multiple different actions over reasonable time scales. The fact that much is known exactly and perturbatively about the corresponding continuum theory provides an important bonus.

The explicit form of Schwinger Model action which we adopt is

[TABLE]

$n_{f}$ is the number of equal mass fermions included, which we take to be 2 in all our simulations. Gauge elements $U_{\mu}(x)$ are defined on the links of a two dimensional square lattice and take values in the group $U(1)$. The hopping matrix $M_{\kappa}$ is given by the usual Wilson form

[TABLE]

$M$ can be decomposed into its red/black (even/odd) components,

[TABLE]

Using this decomposition, one may construct a Hermitian, positive definite, red/black preconditioned matrix $H^{(e)}$ such that

[TABLE]

For two flavours of equal mass fermions, the fermion contribution to the Schwinger model action is then given by

[TABLE]

In what is to follow, we will use preconditioned iterative methods to calculate inverses and logs of the fermion matrix. We have found that the red-black form here provides the best preconditioning.

In the language of the last section, our true action is the Schwinger model action with dynamical fermions included and is given by (62). We used the hybrid molecular dynamics algorithm [7] with a five step integration scheme [8] to generate ensembles of configurations distributed according to this true action. All the approximate actions which we considered involved varying sets of Wilson loop operators. To generate ensembles distributed according to these actions, we used a simple Metropolis link update algorithm.

The only other simulation issue worth commenting upon is that of the calculation of $T$ and $T^{2}$ for individual configurations in an ensemble. These calculations are required whenever we need to find correlations between action terms. $T$ is defined in (67) as the trace of the log of the matrix $H^{(e)}$. To calculate $T$ we adopt a stochastic approach [1]. A number $N_{\phi}$ of Gaussian distributed vectors $\phi_{a}(x)$ are generated. The index $x$ here runs over the (even) sites of the lattice, and $a=1,\dots,N_{\phi}$ labels the different vectors generated. We use the notation $\langle\ \rangle_{\phi}$ to denote expectations over the vectors $\phi$ and the probability distribution of the components of these vectors are normalized so

[TABLE]

The stochastic estimators, $E_{T}$ and $E_{T^{2}}$, which we adopted for $T$ and $T^{2}$ respectively, are given by

[TABLE]

where, for convenience, we define

[TABLE]

We have immediately that

[TABLE]

Note that $E_{T^{2}}$ as defined here is an unbiased estimator of $T^{2}$. The other obvious choice of estimator for $T^{2}$ is $(E_{T})^{2}$. We see from these equations that $E_{T}^{2}$ is a biased estimator for $T^{2}$ which only becomes unbiased in the limit $N_{\phi}\rightarrow\infty$. The variances of the estimators $E_{T}$ and $E_{T^{2}}$ are also easy to calculate. We have

[TABLE]

The variance in $E_{T^{2}}$ has a more complicated form, but a similar $1/N_{\phi}$ overall factor. When using these estimators to calculate $T$ and $T^{2}$, we must therefore choose $N_{\phi}$ sufficiently large that the fluctuations due to $\phi$ are small.

The remaining technical problem is to calculate $L\phi=\log H^{(e)}\phi$. We use a Chebychev polynomial approximation (of order $N_{C}$) to perform this calculation. To optimize the polynomial approximation used, it is necessary to know the maximum and minimum eigenvalues of the matrix $H^{(e)}$. These values are obtained with a Lanczos method [9]. The calculation of the Chebychev polynomial approximation itself then requires $N_{C}$ matrix-times-vector products for each vector $\phi$. In all our analyses, we have found that the time needed to estimate $T$ and $T^{2}$ is a significant fraction (5%-50%) of the time needed to perform a full dynamical simulation.

Arbitrary numbers of equal mass fermions can, of course, be dealt with by incorporating a factor $n_{f}/2$ along with $T$ in the accompanying formalism.

V Lowest Order Results

A The Quenched Approximation

Our preliminary application of the methods introduced in Section III to the Schwinger model was performed on lattices of size $16^{2}$ and $32^{2}$, at $\beta$ values from 2.0 to 3.0, and $\kappa$ from $0.25$ to $0.265$. Typically, we generated 1000 or 2000 configurations using either the true full fermion action or an approximate action. Maximum and minimum eigenvalues of $H^{(e)}$ were calculated on each of these configurations using a Lanczos algorithm. The major computational effort is then to evaluate the stochastic estimators $E_{T}$ and $E_{T^{2}}$ of $T$ and $T^{2}$ respectively on each configuration. We used 30-50 different Gaussian $\phi$’s per configuration to accumulate these estimates, and we used Chebychev polynomials of order 100-200 to calculate $L\phi\equiv\log H^{(e)}\phi$ for each $\phi$. We normally worked with configurations which were first fixed to Landau gauge. This was expected to reduce the influence of gauge artifacts in the Lanczos procedure and allowed the use of Fourier acceleration methods. In practice, we found that preconditioning additional to that provided by the red/black construction gave little extra benefit.

We also accumulated measurements of a variety of different Wilson loop operators on each configuration. As described in Section III D, we will use the notation ${\cal W}_{i}$ for such a loop operator labeled by the generic index $i$. Where appropriate, the generic index will be replaced by a descriptive index. For example, ${\cal W}_{1\times 1}$ denotes the simple $1\times 1$ plaquette Wilson loop operator which, for the Schwinger model, is defined by

[TABLE]

Typical loop operators considered include $1\times 1$, $2\times 1$, $2\times 2$, $3\times 1$ etc.

When the raw data for the estimators $E_{T}$, $E_{T^{2}}$ and the loop operators ${\cal W}_{i}$ have been calculated on a configuration by configuration basis, the final analysis required is to evaluate correlations of different loop operators with $E_{T}$, and to derive the values of the approximate action coefficients according to the conditions of Section III D.

The first and simplest case to consider is the quenched approximation. The true action (for two flavours of equal mass fermions, (62)) is

[TABLE]

The approximate action which defines the quenched approximation includes just the single plaquette Wilson loop operator, ${\cal W}_{1\times 1}$. The parameters of this approximate action must however be adjusted according to the conditions in Section III D and so we write it with a different coupling, $\beta^{\prime}$. Thus,

[TABLE]

The difference between these two actions is

[TABLE]

where

[TABLE]

$\delta\beta$ is the shift in plaquette coupling induced by approximating $T$ with a plaquette term only. Our lowest order tuning condition is that the variance, $\sigma^{2}(\Delta)$, is to be minimized in either the true or approximate gauge measure. Explicitly we have

[TABLE]

The condition on $\delta\beta$ to minimize this expectation is

[TABLE]

In the initial test, we used the quenched algorithm to generate configurations. Thus the ensemble which we have generated is the appropriate one for the approximate action, and we choose to minimize the variance of $\Delta$ in this measure. For $\beta^{\prime}=2.50$ and $\kappa=0.26$, we found $\delta\beta=-0.212\pm 0.010$. These and all other errors were evaluated using the bootstrap method, with subensembling to detect and remove autocorrelations. The interpretation of this result is that, within the space of quenched actions defined by (78), the best approximation to the true action (77) with $\beta=2.500-0.212(\pm 0.010)$ and $\kappa=.26$ is given by the quenched action with $\beta^{\prime}=2.500$.

Having defined an approximate action, it is immediately important to check the quality of the approximation which that action generates. In particular, we ask whether we have been justified in using just the leading term in the cumulant expansion when tuning the approximate action. To investigate this, we generated a new ensemble of 2000 configurations using the hybrid molecular dynamics algorithm to incorporate dynamical fermions ($\kappa=0.26$) at $\beta=\beta^{\prime}+\delta\beta=2.29$. As a simple preliminary test, we compared $\langle{\cal W}_{1\times 1}\rangle$ as measured using the approximate and true actions. The quenched, approximate action value was

[TABLE]

while the true dynamical fermion action value was

[TABLE]

The error $\pm 1.4$ in this last equation is due to the uncertainty with which we have determined $\delta\beta$. The agreement between the two measurements is a very encouraging result. The expectations of ${\cal W}_{1\times 1}$ match in the true and approximate measures. We are therefore in the situation where tuning using the common lowest order tuning conditions (45) produces an approximate action which matches the expectation of an operator in the approximate and true measures. This is our first evidence that, in this instance, we have been justified in taking only the leading term in the cumulant expansion.

B Bound state mass

We now investigate how well the quenched approximation performs in calculating expectations of some interesting physical operators in the Schwinger model. The continuum model with a single species of massless fermions has a ‘vector’ bound state with mass $e/\sqrt{\pi}$ where $e$ is the fermion charge. For $n_{f}$ flavours of massless ‘quarks’, there is an $SU(n_{f})$ symmetry, leading to an iso-singlet ‘vector’ particle with mass

[TABLE]

and a massless $(n_{f}^{2}-1)$-plet [10]. On the lattice we have two flavours and look for a bound state by studying the large Euclidean time dependence of zero momentum correlators involving $\overline{\psi}\gamma_{1}\psi$. The latter were constructed from ‘quark’ propagators which were obtained by a Conjugate Gradient solver applied to the red/black preconditioned Hermitian matrix $H^{(e)}$ or by the BiCGStab algorithm [11] used with $M^{(e)}$. The latter turned out to be the more efficient method.

Figure 1 shows the vector state effective mass $m_{v}$ obtained from a full dynamical fermion simulation with the true action at $\beta=2.29$, $\kappa=0.26$ compared with that obtained from a quenched simulation at $\beta^{\prime}=2.50$. The quenched values plotted for $m_{v}$ include the leading ‘deficit correction’ which is defined by the first term on the right hand side of the cumulant expansion formula (LABEL:Eqn:FCumExpans). The dashed line in this figure serves both to indicate the (uncorrected) value for $m_{v}$ obtained by a quenched simulation with $\beta^{\prime}=2.29$, and the Euclidean time range used to extract mass estimates. From the true action simulation and from the deficit corrected estimate from the quenched simulation respectively, we find $m_{v}=0.369(3)$ and $m_{v}=0.367(7)$ where the bracketed digit is the error in the last significant digit. These vales are in complete agreement. In the deficit corrected value, the contribution due to the $\beta^{\prime}$ shift from $\beta^{\prime}=2.29$ to $\beta^{\prime}=2.5$ is $-.041$ and that due to the deficit correction term is $-.037$.

This is pleasing but not yet surprising. The approximate action simulation is again generating predictions which are in agreement with those from the correct true action simulation. It must, however, be pointed out that the statistical fluctuations in the approximate action estimates are, for the same number of configurations, somewhat larger ($\times 2.5$). So any comparison of algorithm efficiency must take this into account. We will return to this issue later.

The relatively small difference between quenched and full theory estimates in this particular observable is in accord with experience gained in lattice studies of the hadron spectrum. It also means that there is little point in using $m_{v}$ when studying more complicated approximate actions. To really stretch the method we now study an observable which is explicitly dependent on the long distance dynamics.

C Static potential

Figure 2 shows a comparison of the static potential using quenched simulation data at $\beta^{\prime}=2.50$, full dynamical simulation data at $\beta=2.29$, $\kappa=0.26$ and the deficit corrected estimate from the quenched data. The results have been expressed in dimensionless units using our result $m_{v}=0.37$ in lattice units. Also shown is the (exact) infinite volume quenched result (dashed line) and the continuum massless Schwinger result (full line) obtained by evaluating Wilson loops due to a vector field of mass $m_{v}$. The continuum result is

[TABLE]

Note that, expressed this way, the continuum result is independent of the number of flavours. The comparison between the lattice estimates is satisfactory although the errors associated with the deficit estimate (LABEL:Eqn:FCumExpans) get quite large at large distances. This could signal a systematic under-estimate of fermion effects at increasing $R$, but it is not yet statistically significant. The same behaviour is seen in a direct comparison using Wilson loops with increasing area (see below). It seems therefore, that the potential screening effects due to light dynamical fermions are capable of being reproduced by these first order estimates based on quenched configurations.

VI Extensions

A Naive Improvement

The results of the last section are quite encouraging. Approximating $T$ with a $1\times 1$ plaquette is seen to shift significantly the quenched results towards the full theory results. The absolute measure of how well we have actually done was defined in the discussion following (36) and is given by the corresponding values of $\sigma^{2}(\Delta)$. For the $32^{2}$ lattice at $\beta=2.50$, $\kappa=0.26$, we find $\sigma^{2}(T)=99.0(2)$, while $\sigma^{2}(T+\delta\beta{\cal W}_{1\times 1})=3.1(7)$. These values allow us to predict acceptances according to (38) for algorithms which first generate trial configurations, then accept/reject those configurations according to the change in $T$ or $T+\delta\beta{\cal W}_{1\times 1}$ respectively. The value of $\sigma^{2}(T)=99.0(2)$ implies that an algorithm which uses the pure gauge term of the Schwinger model at unshifted $\beta$ to generate trial configurations will have an acceptance of $2\times 10^{-12}$ (${\rm erfc}(0.5\sqrt{99.0})$). The value $\sigma^{2}(T+\delta\beta{\cal W}_{1\times 1})=3.1(7)$, on the other hand, implies that an algorithm which uses the pure gauge term at a properly tuned shifted $\beta$ will have an acceptance of 0.2. For the Schwinger model, at the parameters we have analyzed, the $1\times 1$ loop is therefore seen to be a reasonable approximation to the fermion determinant, since it produces a finite (non zero) value of acceptance. It is not perfect, however, since acceptance is still quite small.

We wish to consider now whether, and to what extent, the addition of further loops to the approximate action can improve the situation. For this and all remaining analysis in this paper, we will concentrate on the effects of approximate actions on loop operators only. The single loop analysis just described shows that the loop expectation values are the most sensitive to changes in action. The effective mass is rather insensitive, and in any case, is already quite well described by the addition of just the $1\times 1$ plaquette term in the approximate action.

The first addition which we considered involves an approximate action where we included loops containing up to three plaquettes. The loops involved are the $1\times 1$, $1\times 2$, and $1\times 3$ rectangles (denoted by 1, 2, and 3 for short), together with a ‘chair’ shaped loop of length 8 (denoted by $3^{\prime}$). These four loops are shown in Figure 3 and are some, but not all, of the simple loops of perimeter length 8 or less. An immediate problem with this loop expansion is that the loop possibilities expand very rapidly as the loop length increases. Even at this simple expansion level, we are already ignoring some non-trivial length 8 loops (for example, the $1\times 1$ squared loop operator).

The target full fermion action of our analysis has $\beta=2.50$ and $\kappa=0.26$. Expressed in terms of Gram-Schmidt orthogonalized loop operators, the generic approximation of $T$ is now

[TABLE]

Initially we have no information about the values of the coefficients $a_{i}$, so we set them to zero and simulate with an approximate action which contains only the pure gauge part of the Schwinger action with $\beta=2.50$ (*i.e. *we simulate with the unshifted quenched version of the Schwinger action). The ensemble so generated allows us to make new estimates of the coefficients $a_{i}$, which can be used to generate a new approximation. We repeated this process four times in all. Table I summarizes the sequence followed.

Ensembles of 2000 configurations were generated for each approximation, and the approximate action coefficients for the next approximation were calculated on these ensembles using the formulae of Section III E. The results achieved are given in Table II. Approximation 0 in a quenched simulation with $\beta=2.50$. The value which results for $a_{1}^{(0)}$ from this approximation is $0.212(10)$ which is exactly the value of $-\delta\beta$ determined in the last section. Here however, we have specified that the true action has $\beta=2.50$. Thus we interpret $a_{1}^{(0)}$ as defining the best quenched approximation coupling shift to use to simulate the full theory with $\beta=2.50$. This best quenched theory has $\beta^{\prime}=2.50+0.212(10)$ and defines our second approximation (approximation 1 in Table I). Each approximation of course generates a different ensemble, and the values for the approximate action coefficients $a_{i}$ therefore change from approximation to approximation. This is clear also in the table. But a very nice feature to note is the fundamental stability of the procedure we have adopted. The coefficients which result after approximation 3 are, within errors, the same as the coefficients after approximation 0.

For the highest approximation attempted (approximation 3), we also list the corresponding root variance $\sigma(\hat{{\cal W}}_{i})$ and the product $a_{i}\sigma(\hat{\cal W}_{i})$. We have, when the coefficients $a_{i}$ are tuned according to 61, that

[TABLE]

Thus, $a_{i}\sigma(\hat{\cal W}_{i})$ is a measure of the importance of a particular loop operator in approximating $T$. If the addition of loops of increasing size generated a rapidily converging approximation to the true action, we would expect $a_{i}\sigma(\hat{\cal W}_{i})$ to decrease rapidly as the loop size increases. It is clear from the table, however, that this is not the case for the Schwinger model. The contributions of two and three plaquette loops are not significantly smaller than the contribution of the single plaquette loop.

In Figure 4, we show the static potential, at fixed distance in lattice units (for $R=a$ and $R=5a$), versus the approximation number, $n$. We show both the raw measurements made using each set of configurations, and the deficit-corrected estimates given by (LABEL:Eqn:FCumExpans) as described above. The corresponding dynamical fermion result (2,000 configurations) is also shown for each inter-quark distance. As expected, we see that adding loops of increasing size has very little effect on the short distance potential at $R=a$. However, at $R=5a$ there is an effect, and a decreasing trend is evident in the (uncorrected) data obtained as larger loops are added to approximate action. This decrease is quite slow unfortunately, and confirms that the systematic addition of loops of increasing size produces only a slowly converging approximation to the true action.

One significant positive observation can made from Figure 4. This observation, anticipated in the previous section, is that at each approximation order, $n>0$, the deficit-corrected estimate given by (LABEL:Eqn:FCumExpans) is compatible with the dynamical fermion result. At first sight, it may therefore seem that there is no need to go beyond $n=1$, *i.e. *to use other than quenched data to make properly corrected estimates of the full theory. However, as emphasized in Section III, one can only be confident of such cumulant expansion-based estimates when the shifts involved are known to be small. This will only be so when the corrections are based on a suitably accurate approximate action.

Note finally that the corrections are not simply shifts in $\beta$, but are non-trivial effects of dynamical fermions. One may think of the estimates as a form of re-weighting in which the approximate vacuum used becomes closer to the true vacuum as more loops are added. However, what is very clear is that the addition of small loops only in the approximate action improves physics only at the scales appropriate for small loops. Thus the static potential at short range is almost perfectly determined by the approximate action. At longer range, however, the approximate action is not very successful in improving the quenched result.

B Selective Improvement

This analysis suggests that an approximate action based on a systematic expansion of $T$ in terms of loops of increasing size is going to require the inclusion of a very large number before convergence is achieved. However, the impression already gained is that, once a loop of a particular size is added to the approximate action, the addition of more loops of a similar size is not particularly productive. This suggests an alternative approach of adding loops selectively rather than systematically. The aim would be to include typical loops of different sizes rather than all loops up to a given size. As an experiment, we chose to consider approximate actions with $1\times 1$, $3\times 3$ and $5\times 5$ loops included in turn. For simplicity, tuned coefficients defining all these actions were determined on a single quenched, $\beta=2.50$ ensemble.

The results achieved by these new approximate actions are given in Table III. The leftmost column shows the values for various loop expectations from a full dynamical fermion simulation. The remaining columns gives the values achieved by the various different approximate actions. What we see is that the expectation value achieved for a loop in a simulation which includes that loop in the approximate action correctly matches its full dynamical simulation value. This is in good agreement with our analysis of Section III where we argued that our approximate action tuning conditions are such as to match loops included in the approximate action to their correct values in the true action. Loops not explicitly included in the approximate action are not so well tracked. There is however considerable improvement as we add more loops to the approximate action. For example, the $3\times 8$ loop changes from value 3.7 when only the $1\times 1$ term is included in the approximate action to the value 10.5 when $1\times 1$, $3\times 3$, and $5\times 5$ loops are included. This is to be compared with its exact full theory value of 14.4.

A check on the results achieved is to execute an iteration of the tuning algorithm to see if our results are sensitive to the choice of coefficients in the $1\times 1+3\times 3+5\times 5$ approximate action used. To execute this iteration, we used our $1\times 1+3\times 3+5\times 5$ approximate action to generate a new ensemble on which we recalculated the coefficients $a_{1\times 1},a_{3\times 3}$ and $a_{5\times 5}$. Table IV shows the results of this retuning. Again we find that the results are quite stable. In the Schwinger model, a set of quenched configurations seems sufficient to calculate the approximate action coefficients.

C Optimized Improvement of The Effective Action

As a final test of the method, we attempt a very direct approach to defining the effective action by using an ensemble of configurations generated using the full fermion algorithm. The ensemble is for a $16^{2}$ lattice at $\beta=2.50$ and $\kappa=0.26$. In all we generated 1,000 configurations, and systematically considered a large number of different possible approximations for $T$ of the form

[TABLE]

where the index $i$ here runs over some selected subset of the Wilson loop operators which we have calculated. For example, one possible approximation to consider is $T\approx a_{2\times 2}{\cal W}_{2\times 2}+a_{1\times 3}{\cal W}_{1\times 3}$. For each possible approximation, we choose values of $a_{i}$ which minimize $\sigma^{2}(T-\sum_{i}a_{i}{\cal W}_{i})$ for that approximation, and then determine the approximate action which generates the minimum value for $\sigma^{2}$. The results of this minimization are listed in Table V, and the optimum loop choice is shown in Figure 5. The first line of the table gives the value of $\sigma^{2}(T)$ with no approximation attempted. The second line indicates that the best single loop approximation to $T$ is obtained when the loop is the $3\times 3$ loop. The best two loop approximation is obtained by using the $3\times 3$ and $1\times 6$ loops, etc. The measure of the quality of a given approximation is the value of $\sigma^{2}(T-\sum_{i}a_{i}{\cal W}_{i})$ listed in the final column of this table. The best one loop approximation reduces $\sigma^{2}$ from $2.69$ to $1.02$. Adding more loops reduces $\sigma^{2}$ further, but not at the same rate as adding only a single loop.

D Stability

In all cases considered so far, we have tuned the effective action by minimizing $\sigma^{2}$. In Section III we proposed three different tuning conditions. These all agree to lowest order in the cumulant expansion. They differ at higher order, however. To test whether the differences between the three tuning conditions are significant, we compared the coefficients for the effective actions tuned according to (18), (22) and (35). The results of this comparison are shown in table VI. The values determined by minimizing $\sigma^{2}$ for the coefficients of the $3\times 3$, $1\times 6$, $5\times 5$ and $1\times 2$ loops in the best four loop approximation for $T$ are given in the second column of the table. If instead, we tune to match the expectations of the four loops in the approximate action between the true and approximate measures, then the values of the coefficients are shifted by the amounts listed in the third column of the table. If we tune to maximize acceptance, then the shifts necessary are given in the final column of the table. In all cases the shifts are very small, compatible with zero at this level of statistics. This is indicative that we have found an approximate action which is sufficiently close to the true action that all the tuning conditions which we have considered are self-consistent.

Having arrived at a tuned approximate action (*e.g. *coefficients given in Table VI), we have made preliminary tests of an exact algorithm as described in Section III C. A similar exact algorithm was introduced in the context of the bosonic algorithm [5] by Peardon [12] and further developed by others [13]. These tests have shown that correct results can indeed be obtained with an experimental acceptance as predicted in Section III C. However, careful tuning is required when applying the algorithm to large systems. The approximate action must yield $\sigma^{2}(\Delta_{at}$ not too different from $1$ so that reasonable acceptances are obtained.

VII Application to QCD

So far, the main numerical results in this paper have concerned the application of approximate actions to the Schwinger model. We chose this model because it allowed useful comparisons between the ensembles generated by an exact full fermion algorithm and those generated by a variety of different approximate actions in reasonable computational time. However, our main interest in the application of approximate actions is in lattice QCD. A complete repetition of the analysis so far described for lattice QCD is a considerable computational undertaking, and beyond the scope of this exploratory paper. Nevertheless, it has been possible to gain some indication of both the potential of the method and the practical problems which still remain for lattice QCD.

The exploratory simulation which we have carried out is for full QCD on a small lattice. The true action which we considered was the standard single plaquette pure gauge action with coupling $\beta$, combined with a two-flavour Wilson fermion action with hopping parameter $\kappa$. The lattice size adopted was $6^{4}$ with $\beta=5.7$ and $\kappa=0.16$. This is obviously a very small lattice with very heavy quarks. The basic goal was to determine the extent to which it is possible to approximate the true action for this theory with an action made up from a finite set of Wilson loops.

The direct approach to analyzing this question follows exactly the strategy adopted at the end of the last section, and is simply to generate an ensemble of configurations using the full QCD action, to evaluate $T$ and various different loops on these configurations, and then to determine the set of loops which best approximates $T$. This is the approach we have adopted. We used hybrid molecular dynamics to generate configurations. The ensemble analyzed contained 516 configurations. Each configuration was separated by 10 trajectories of length 0.5. $T$ was calculated as for the Schwinger model on each configuration. We also calculated a variety of different Wilson loop operators on each configuration. These loops were chosen to be representative rather than complete. Four dimensions allow many more possible loop configurations than two dimensions, and an exhaustive study of all the loop possibilities is not justified in the present instance where we are working with a very unphysical set of lattice parameters. Our goal is simply to explore the possibilities rather than to present a definitive calculation.

Following our experience with the Schwinger model, we looked at a variety of different loop shapes and sizes. These included two, three, and four dimensional loops, with total length varying from 4 to 12 lattice spacings. The only non-trivial loop of length 4, of course, is the pure gauge $1\times 1$ plaquette action, while the length 12 loops considered included $5\times 1$, $4\times 2$, and $3\times 3$ 2-dimensional loops together with more complicated three and four dimensional constructions. A particular loop construction is defined by specifying the steps in the various possible directions which a link in that loop can represent. Thus the specification $\{12\bar{1}\bar{2}\}$ denotes a loop generated by first stepping in the direction $1$ with positive sense, then in direction $2$ with positive sense, then in direction $1$ with negative sense, and finally in direction $2$ with negative sense. This set of four steps brings the loop operator back to its starting site. We can then take the trace in the usual way. The digits $1$ and $2$ can represent any two distinct directions in the four dimensional lattice. An unbarred digit defines a step in the positive direction along the corresponding axis, a barred digit defines a step in the negative direction. The specification implicitly assumes a sum over all independent orientations and starting locations of the loop. For example, the specification we have just written represents the simple $1\times 1$ plaquette summed over all orientations and lattice sites, and is therefore the standard one plaquette pure gauge action term.

To test how well the true action can be approximated, we then seek to minimize $\sigma^{2}(T-\sum_{i}a_{i}{\cal W}_{i})$ over the different loops considered. In all, we calculated over 40 loops. This allows us to consider a very large number of different approximate actions. These possibilitys include actions with a single loop include, actions with two loops included etc.

Table VII shows some of the results we achieved. The first column lists the variety of loops which we considered. The order in which these loops are listed represents their relative importance in reducing fluctuations in the quantity $\sigma^{2}(\Delta)$ and is determined only after all analysis is complete. The second column lists the minimum achieved for $\sigma^{2}(\Delta)$ when just a single loop is included in the approximate action. This minimum, of course, varies depending on which particular loop is considered as the single term approximate action, and the second column shows this variation. The raw value of the variance in $T$ is listed at the end of the table for comparison purposes.

It is immediately noteworthy that there is a large variation in the values achieved for $\sigma^{2}(\Delta)$, indicating that there is a large variation in the effectiveness with which different loops can individually approximate $T$. The effective loops tend to be loops with four dimensional structure, and it is quite striking that the simplest loop of all, the single plaquette loop (defined as $\{12\bar{1}\bar{2}\}$) is actually quite unimportant. The best match between true and approximate actions is given by the loop which is the first entry in the table. This loop, on its own, reduces the fluctuations in $S_{t}$ from a value of $295(20)$ to a value $8.7(4.4)$. This is a remarkable 97% reduction in the raw fluctuations of $S_{t}$. The loop which does this is not a simple loop but rather an oddly shaped four dimensional construction which has total length 12. We have not been able to identify any remarkable features about this particular loop relative to all the other loops which we considered. Other similarly shaped loops are seen to reduce the fluctuations in $S_{t}$ almost as much, so it is reasonable to imagine that the important feature is the four dimensional nature of the loop and its approximate size rather than the exact details of how the loop is generated.

The determination of the best single loop approximate action is only the first analysis step. The next is to determine the best two loop approximate action. Our approach is exactly as before. We minimize $\sigma^{2}(\Delta)$ for all two loop approximate actions. In fact, we only need to minimize $\sigma^{2}(\Delta)$ over loops with $1$ and $n$ where $n$ is any of the other available loops. This minimization determines that the approximate action which uses the loops placed 1 and 2 in the table is the best two loop approximate action. The final column of Table VII indicates the value which results for $\sigma^{2}(\Delta)$ from this minimization. The minimization for higher numbers of terms in the effective action can now proceed by iteration. The final result we find is that the best $n$ loop approximate action is the action which contains loops $1,2,\ldots,n$. The values which result for $\sigma^{2}(\Delta)$ with this best $n$ term approximate action are all given by the final column of Table VII.

The second interesting feature of this exploratory study is that the addition of extra loops beyond the first loop in the approximate action is not particularly helpful. The approximate action with a single loop reduces the raw fluctuations in $S_{t}$ from $295(20)$ to $8.7(4.4)$. Adding a second and third term reduces these fluctuations only to $7.2(4.1)$ and $5.95(4.00)$ respectively. These shifts are very much smaller than that achieved by a single term alone.

As even more loops are added, the error in $\sigma^{2}(\Delta)$ becomes greater than its actual value, and we conclude that we have insufficient data to calculate reliably the best approximate action for more than a small number of terms, so the actual ordering after the first three or four terms should be considered suspect. However, a significant observation which can be made is that the $1\times 1$ plaquette loop contribution comes at position 35, which is very far down the list. We conclude therefore that the pure gauge plaquette action is actually a rather poor approximant of the the full fermion action.

VIII Conclusions

We have proposed a general framework for the construction and evaluation of approximate actions for lattice QCD. We have demonstrated how the coefficients in the simplest form of such actions can be evaluated directly from a true action simulation or by using an iterative procedure which begins with an approximate action simulation, and never requires simulation with the true action. We have also shown that the results of the iterative procedure are quite stable, and converge very quickly.

As a first application, we have tested approximate actions based on the systematic expansion of the trace of the log of the fermion matrix in Wilson loops of increasing size as originally proposed in [1]. This expansion works well for observables sensitive to short range physics. Using the Schwinger model, we were able to test longer range observables where it becomes clear that an impractically large number of loops will be required to study physics over a range of scales using small lattices spacings.

The systematic error correction estimate proposed by [1] is shown to give reliable estimates of observables based on approximate action simulations. The better the approximate action, the more reliable these ‘deficit’ estimates become. However, since these estimated corrections require a good, unbiased, stochastic estimate of $T$, the trace of the log of the fermion matrix, their practical value is dependent on further progress in developing efficient calculational methods for this stochastic estimation. The computational methods which we have used here to evaluate $T$ are quite demanding. Typically they require work equivalent to that needed to evaluate all meson and hadron propagators at a number of different hopping parameters on a given configuration.

Having established the above behaviour, we have gone on to show that the greatest determinant of the quality of an effective action is not the number of loops included, but the choice of scales represented by the loops which are included. We have found dramatic evidence of this both for the two-dimensional test model and for QCD itself. The systematic matching and tuning methods have allowed us to propose approximate actions significantly closer to the true QCD action than the quenched approximation or, indeed, of an approximation based on the ordered inclusion of loops of increasing size.

We presented two methods for capitalising on approximate actions in lattice simulation. The first, as originally proposed in [1], is to simulate with the approximate action and to make (presumably) reliable estimates of the small corrections to desired observables. The challenge to progress in this strategy is to improve stochastic evaluation methods for $T$. The second method described is to use the approximate action to generate candidate update configurations within an exact algorithm, and we proved by construction that such an exact simulation was feasible in the Schwinger model example. The challenge to the construction of such an exact algorithm in QCD is to implement the fine tuning required to achieve viable acceptance rates. It is perhaps too optimistic to imagine that this fine tuning can be achieved by considering actions made of Wilson loop operators only. The QCD test described here suggests that the return achieved in adding more Wilson loops to the approximate action becomes very poor after the first few loops. On the other hand, the tuning methods we have defined are completely general, and can be applied to other candidates for approximate actions. In particular, a combination approximate action including a small number of Wilson loop operators and a multiboson Lüscher [5] term seems to us to be an interesting possibility as a tunable candidate approximate action.

While awaiting these developments, an immediate and exciting application of the methods we have described is to perform relatively high statistics measurements using the approximate action directly. The strategy here is to make low statistics measurements of the quantity of interest using the full theory. These can then be used to help tune an approximate action to reproduce this measurement approximately (and other more accessible ones to higher levels of accuracy). The approximate action can then be used to give high statistics estimates of the quantity of interest (*e.g. *glueball mass).

IX Acknowledgements

One of us (J.S.) would like to thank D. Weingarten for much discussion on the subject matter of this paper. We would also like to thank The Hitachi Dublin Laboratory for providing computer resources which we used to perform some of the calculations described here.