format: text

text:
# Dynamics of bright matter wave solitons in a quasi 1D Bose-Einstein condensate with a rapidly varying trap

## I Introduction

Recently bright matter wave solitons has been observed in a Bose-Einstein condensate. In the experiment described in Ref. propagation in anisotropic BEC with cigar type geometry was considered. In Ref.  the soliton is monitored by projecting the bound state of $``$5000 atoms into expulsive harmonic potential. The soliton was observed propagating without changing of the form for distances is order of $`1`$mm. The expulsive harmonic trap used, in the experiment, corresponds to an unstable potential for the GP equation and a soliton can exists only as a metastable state in BEC . It is of interest to investigate the dynamics of a bright matter wave soliton in inhomogeneous time dependent systems, in particular possible stabilization of unstable dynamics or complicated dynamics of solitons in BEC.

An interesting example of such a trap is an optical trap which is indeed rapidly varied in time. It is realized for instance for a condensate in dipole trap formed by a strong off-resonant laser field . A typical model is a trap formed by a harmonic potential with different frequencies which is cutting at an energy $`V_c`$. Numerical simulations of the condensate dynamics and splits in 1D BEC with a positive scattering length under periodically shaken trap has been performed in Ref.. The same form of wave equation leads to the propagation of spatial solitons in periodically modulated parabolic waveguides. This problem appears in investigation of intense light beams in a nonlinear waveguide with inhomogeneous distribution of transverse refractive index. A typical distribution can be approximated by a quadratic profile. The rapid variation of this profile along the longitudinal direction leads to a rapidly varying quadratic potential in the NLS equation.

The general mathematical problem is to investigate localized states for the nonlinear wave equation with cubic nonlinearity and with rapidly varying (not small) potential of the form $`V_0=f(x)\alpha (t/ϵ),ϵ1`$. Also it is of interest to investigate possibilities of stabilization of a soliton by rapid perturbations. For multi dimensional NLS equation under rapidly varying in space periodic potential such a problem is discussed in .

The paper is organized as follows:

In section 2 we describe the procedure of obtaining the averaged equation for the GP equation with an external potential, that is rapidly varying in time and inhomogeneous in space, with the form $`f(x)\mathrm{cos}(\mathrm{\Omega }t)`$,where $`\mathrm{\Omega }=1/ϵ`$. In Section 3, we considered some important potentials for applications, which have the forms given by $`f(x)=x^2,1/\mathrm{cosh}(bx)`$, and $`\mathrm{cos}(kx)`$.

## II Averaged Gross-Pitaevskii equation

The dynamics of BEC is described by the time dependent GP equation

$$i\psi _\tau =\frac{\mathrm{}^2}{2m}\mathrm{\Delta }\psi +V(r,t)\psi +g_0|\psi |^2\psi ,$$
(1)
where $`m`$ is the atom mass, $`g=4\pi \mathrm{}^2a_s/m`$,with $`a_s`$ the atomic scattering length. $`a_s>0`$ corresponds to the BEC with repulsive interaction between atoms and $`a_s<0`$ to the attractive interaction. The trap potential is given by $`V=m\omega ^2(y^2+z^2)/2+\alpha (t)(m\omega _1^2x^2/2+V_1(x,t))`$, where $`V_1(x,t)`$ is bounded potential or the optical lattice potential, $`\alpha (t)`$ describes the time dependence of the potential. We will specify the form of $`V_1`$ later. Below we will consider the cigar type condensate with $`\omega ^2\omega _1^2`$. Within such restrictions we can look for the solution of (1) with form $`\psi (x,y,z,t)=R(y,z)U(x,t)`$, where $`R`$ satisfies the equation

$$\frac{\mathrm{}^2}{2m}\mathrm{\Delta }R+\frac{m\omega ^2}{2}(y^2+z^2)R=\lambda _\rho R,\lambda =\mathrm{}\omega .$$

Averaging over the transverse mode $`R`$(i.e. multiplying by $`R^{}`$,

$$|R_0|^2=\frac{m\omega }{\pi \mathrm{}}\mathrm{exp}(\frac{m\omega }{\mathrm{}}\rho ^2),$$

and integrating over $`\rho `$) we obtain the quasi 1D GP for $`U`$

$$i\mathrm{}U_\tau =\frac{\mathrm{}^2}{2m}U_{xx}+(\frac{m\omega _1^2x^2}{2}+V_1(x,t))U+G|U|^2U,$$
(2)
where $`G=g|R|^4𝑑x𝑑y/|R|^2𝑑x𝑑y=(2\mathrm{}|a_s|\omega )`$. In the dimensionless variables

$$t=\omega \tau /2,x=x/l,l=\sqrt{\frac{\mathrm{}}{m\omega }},u=\sqrt{2|a_s|}U,$$

we have the governing equation

$$iu_t+\beta u_{xx}+2\sigma |u|^2u=\alpha (t)f(x)u.$$
(3)
where we will assume a periodic modulation of the trap

$$\alpha (t)=\alpha _0+\alpha _1\mathrm{sin}(\mathrm{\Omega }t).$$
(4)
For example for $`V=m\omega _1^2x^2/2,\alpha _0f(x)=(\omega _1/\omega )^2x^2.`$ Also $`\sigma =\pm 1`$ are the attractive and repulsive two body interactions respectively. We introduce the parameter $`\beta =\pm 1`$ to have possibility to use the results for optical beam propagation. For BEC system $`\beta =1`$.

The field $`u(x,t)`$ can be represented in the form of sum of slowly and rapidly varying parts $`U(x,t)`$ and $`\xi (x,t)`$

$$u(x,t)=U(x,t)+\xi (x,t).$$
(5)
For obtaining the equation for an averaged field we will apply the asymptotic procedure suggested in , namely we will present the rapidly varying part of the field as expansion on Fourier series

$$\xi =A\mathrm{sin}(\mathrm{\Omega }t)+B\mathrm{cos}(\mathrm{\Omega }t)+C\mathrm{sin}(2\mathrm{\Omega }t)+D\mathrm{cos}(2\mathrm{\Omega }t)+\mathrm{}$$
(6)
where $`A,B,C,D`$ are functions of (x,t) that are slowly varying in the scale $`O(1)`$ functions. By substituting the equations (5),(6) into (3) we obtain the next set of equations for the slowly varying field and the coefficients of the expansion for the rapidly varying component

$`iU_t+\beta U_{xx}+2|U|^2U+U^{}(A^2+B^2+C^2+D^2)+`$ (7)
$`+2U(|A|^2+|B|^2+|C|^2+|D|^2)+(BCA^{}|A|^2D+ACB^{}+|B|^2D+`$ (8)
$`+ABC^{}+{\displaystyle \frac{A^2D^{}}{2}}+{\displaystyle \frac{B^2D^{}}{2}})=\alpha f(x)U+{\displaystyle \frac{\alpha _1}{2}}f(x)A`$ (9)
$`iA_ti\mathrm{\Omega }B+\beta A_{xx}+4|U|^2A+2U^{}(BCAD)+2U^2A^{}+`$ (10)
$`2U(DA^{}+CB^{}+BC^{}AD^{})+({\displaystyle \frac{3}{2}}|A|^2A+\mathrm{}.)=\alpha _1f(x)U{\displaystyle \frac{\alpha _1}{2}}f(x)D,`$ (11)
$`i\mathrm{\Omega }A+iB_t+\beta B_{xx}+4|U|^2B+2U^2B^{}+2U^{}(AC+BD)+\mathrm{}`$ (12)
$`={\displaystyle \frac{\alpha _1}{2}}f(x)C,`$ (13)
$$iC_t2i\mathrm{\Omega }D+\beta C_{xx}+4|U|^2C+2U^{}AB+\mathrm{}=\frac{\alpha _1}{2}f(x)B,$$
(14)
$$2i\mathrm{\Omega }C+iD_t+\beta D_{xx}+4|U|^2D+U^{}(B^2A^2)+\mathrm{}=\frac{\alpha _1}{2}f(x)A.$$
(15)
The parameters $`\alpha _1,\mathrm{\Omega }`$ are assumed $`1`$. As we can find from this system, the coefficients can be represented in the form of expansion

$`A={\displaystyle \frac{a_1}{\mathrm{\Omega }^2}}+{\displaystyle \frac{a_2}{\mathrm{\Omega }^4}},B={\displaystyle \frac{b_1}{\mathrm{\Omega }}}+{\displaystyle \frac{b_2}{\mathrm{\Omega }^3}}`$ (16)
$`C={\displaystyle \frac{c_1}{\mathrm{\Omega }^3}}+{\displaystyle \frac{c_2}{\mathrm{\Omega }^5}},D={\displaystyle \frac{d_1}{\mathrm{\Omega }^2}}+{\displaystyle \frac{d_2}{\mathrm{\Omega }^4}}.`$ (17)
For the coefficients of expansion we have the expressions

$`a_1=i\alpha _1f(x)U_t\alpha _1(f(x)U)_{xx})2\alpha _1f(x)|U|^2U+\alpha _0\alpha _1f^2(x)U,`$ (18)
$`b_1=i\alpha _1f(x)U,d_1={\displaystyle \frac{\alpha _1^2f(x)^2}{4}}U,c_1=\mathrm{}.`$ (19)
From (7),(16),(18) we obtain the averaged equation for $`U`$

$`iU_t+\beta U_{xx}+2|U|^2U=\alpha _0f(x)Ui{\displaystyle \frac{ϵ^2}{2}}f^2(x)U_t`$ (20)
$`\beta {\displaystyle \frac{ϵ^2}{2}}f(x)(f(x)U)_{xx}2ϵ^2f^2(x)|U|^2U+{\displaystyle \frac{ϵ^2\alpha _0}{2}}f^2(x)U.`$ (21)
Here $`ϵ=\alpha _1/\mathrm{\Omega }`$.

This equation has the conserved quantity - the number of atoms

$$_{\mathrm{}}^{\mathrm{}}𝑑x(1+\frac{ϵ^2f^2(x)}{2})|U|^2=const.$$
(22)
So, it is useful to introduce the new field $`V`$ by

$$V=(1+\frac{ϵ^2f^2(x)}{2})^{1/2}U.$$
(23)
Substituting (23) into (20) and keeping the terms in order $`ϵ^2`$ we obtain the equation

$$iV_t+\beta V_{xx}+2|V|^2V=(\alpha _0f(x)+\frac{ϵ^2}{2}\beta [f_x(x)]^2)V+O(ϵ^4).$$
(24)
The averaged equation has the form of a modified NLS equation with a slowly varying potential

$$W(x)=\alpha _0f(x)+\frac{ϵ^2}{2}\beta [f_x(x)]^2.$$
(25)
This result shows that the soliton dynamics can have more complicated character than in the case of BEC with slowly varying parameters. We can expect here the stabilization of the unstable dynamics of a soliton by rapidly varying perturbation. It is a direct analogy with the stabilization of systems with a few degrees of freedom under rapidly varying in time parameters.

## III DYNAMICS OF DIFFERENT MODELS OF TRAP POTENTIALS

### A A quadratic potential

As an example, we will consider the attractive BEC with a quadratic trap potential. In this case we have a perturbed trap potential $`[\alpha _0+\alpha _1\mathrm{sin}(\mathrm{\Omega }t)]f(x)`$ with $`f(x)=x^2`$ and $`\beta =1`$ . An averaged effective potential is given by

$$W=(\alpha _0+2ϵ^2)x^2,$$
(26)
where $`ϵ=\alpha _1/\mathrm{\Omega }`$. The case when $`\alpha _0>0`$ corresponds to the inverted (repulsive) harmonic potential case which is well known to be unstable . As seen that when $`2ϵ^2>|\alpha _0|`$ the sign of the effective potential is reversed and the solitary wave dynamics becomes stable. The frequency of oscillations of the soliton center now is

$$\mathrm{\Omega }_0=2\sqrt{(2ϵ^2\alpha _0)}.$$
(27)
This result can be also obtained by using the moments method. It is useful to introduce two new variables: a position of a solitary solution of the governing equation (3)

$$X(t)=_{\mathrm{}}^{\mathrm{}}𝑑xx|u(x,t)|^2,$$
(28)
and a field momentum

$$P(t)=_{\mathrm{}}^{\mathrm{}}𝑑x(u^{}(x,t)u_x(x,t)u(x,t)u_x^{}(x,t)).$$
(29)
By using governing equation (3) we come to the following set of integro-differential equations :

$`iP_t+2{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x^{}{\displaystyle \frac{W(x^{})}{x^{}}}|u(x^{},t)|^2=0,`$ (30)
$`iX_t\beta P=0.`$ (31)
In the case of a quadratic effective potential ($`W(x)=\alpha x^2`$) the set of equations (30) is closed and we get the equation of motion for the soliton solution center, $`X(t)`$, in the form

$$X_{tt}+4\beta \alpha X=0.$$
(32)
We point out that the set (30) is exact. No additional assumption was made in obtaining this set of equations, and it is valid for any solitary solution of the governing equation (3). Also, it should be noted two remarkable facts connected with the dynamics of a solitary solution of NLS equation in the presence of a quadratic potential (or dynamics of BEC in a quadratic trap): i).In this case the equations of motion (30) for the soliton center, are exact, ii) after averaging, the initially quadratic potential remains to be quadratic but its sign may be reversed.

Fig. 1. Numerical simulations of the full GP equation for $`f(x)=x^2`$ with $`\omega =10`$, $`\alpha _0=0.0493`$ and values of the perturbation parameter $`\alpha _1`$ providing reversing of the initially inverted harmonic potential. Two cases are presented : a) $`\alpha _1=4.7124`$ and b) $`\alpha _1=3.1415`$.

The results of numerical simulation of the full Gross-Pitaevskii equation with the modulated in time expulsive (inverted) harmonic potential are shown in Fig. 1. Two cases with $`\alpha _1=4.7124`$ and $`\alpha _1=3.1415`$ were considered. For both cases $`\alpha _0=0.0493`$ and $`\mathrm{\Omega }=10`$. The results of computation are the following : a) $`\alpha _1=4.7124`$, $`\mathrm{\Omega }=1.309,(T=4.8)`$ Calculation by Eq. (32) gives $`\mathrm{\Omega }_{cal}=1.257,(T_{cal}=5.0)`$. b) $`\alpha _1=3.1415`$, $`\mathrm{\Omega }=0.7705,(T=8.155)`$ Calculation by Eq. (32) gives $`\mathrm{\Omega }_{cal}=0.7695,(T_{cal}=8.165)`$.

In both cases we observe the appearance of stable oscillations of the soliton position. We obtain good agreement between the full simulations and the averaged equation. As seen, the greater is the value of $`ϵ`$, the discrepancy we found between results of calculation of Eq.(32) and the ones of full numerical simulations of NLS equation. We can expect this since the ratio $`\alpha _1/\mathrm{\Omega }`$ is not small for these values of parameters. Let us estimate the values of parameters for the experimental situation. In the experiment  the longitudinal frequency was $`\omega _1=4`$Hz, so, at the temporal modulations of this potential, with the frequency $`\mathrm{\Omega }=40`$Hz we should observe a stable bright matter wave soliton.

### B The case of a bounded potential

Let us consider a bounded in the space trap potential for the atomic Bose-Einstein condensate  under fast temporal perturbations $`(\alpha _0+\alpha _1\mathrm{sin}(\mathrm{\Omega }t))f(x)`$ with $`f(x)=1/\mathrm{cosh}(bx)`$ and $`\beta =1`$.

After averaging over the period of fast perturbations, the effective potential $`W(x)`$ has the form

$$W(x)=\frac{\alpha _0}{\mathrm{cosh}(bx)}+\frac{ϵ^2b^2}{2}\frac{\mathrm{sinh}(bx)^2}{\mathrm{cosh}(bx)^4}.$$
(33)
Provided that

$$\frac{ϵ^2b^2}{2}>\frac{27}{2}\alpha _0$$
(34)
this effective potential takes a triple well structure having a central minimum and two lateral minima. In Fig. 2, a typical form of the potential is presented, where $`\alpha _0=0.05,\alpha _1=9,\omega =10`$ and $`b=0.6`$. (As seen, the condition (34) makes one to take large ratio of parameters $`\alpha _1`$ and $`\alpha _0`$.)

Fig. 2. Effective potential for $`f(x)=1/\mathrm{cosh}(bx)`$ after temporal averaging. The parameters of perturbation $`\omega =10`$, $`\alpha _0=0.05`$ and $`\alpha _1=9`$. Dashed line is for the case with $`\alpha _1=0`$ when temporal perturbation is turned off.

The two lateral minima are positioned at $`x=\pm 3.96`$ and the central one at the point $`x=0`$.

The results of numerical simulations of the Gross-Pitaevskii equation are shown in Fig. 3.

Fig. 3. Numerical simulations of the full GP equation for the bounded potential with $`f(x)=1/\mathrm{cosh}(bx)`$. The case a) corresponds to the finite motion of the soliton center in the central minimum with $`X_0=0`$. The case b) corresponds to the motion in one of the lateral minima with $`X_0=3.96`$. For both cases $`\omega =10`$, $`\alpha _0=0.05`$ and $`\alpha _1=9`$.

To describe small oscillations of the soliton solution center and to obtain equations of motion in general case, we will have to make some approximations. As the ansatz, we suppose the solution to take the form $`u(x,t)=u(xX(t),t)`$. Here and after, for the sake of simplicity, we also suppose the position of the minimum at $`X_0=0`$ and a symmetrical form for the soliton solution (dynamics of which we are investigating)

$$|u(x,t)|^2=|u(x,t)|^2.$$
(35)
Now we can consider the second term of the equation (30). At first we rewrite the effective potential as a sum of antisymmetric and symmetric parts $`W(x)=W^A(x)+W^S(x)`$, where $`W^A(x)=\frac{1}{2}(W(x)W(x))`$ and $`W^S(x)=\frac{1}{2}(W(x)+W(x))`$. Then, substituting $`W(x)`$ into the second term of Eq. (30) expanding $`|u(xX(t),t)|^2`$ in terms of $`X(t)`$ and by holding the first two terms,

$$|u(x^{}+X,t)|^2=|u(x,t)|^2+X\frac{|u(x^{},t)|^2}{x^{}},$$
(36)
we have

$`2{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x^{}{\displaystyle \frac{W(x^{})}{x^{}}}|u(x^{}+X,t)|^2=2{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x^{}{\displaystyle \frac{W^A(x^{})}{x^{}}}|u(x^{},t)|^2+`$ (37)
$`2X{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x^{}{\displaystyle \frac{W^S(x^{})}{x^{}}}{\displaystyle \frac{|u(x^{},t)|^2}{x^{}}}.`$ (38)
Introducing new parameters

$$\omega ^2=2\beta _{\mathrm{}}^{\mathrm{}}𝑑x^{}\frac{W^S(x^{})}{x^{}}\frac{|u(x^{},t)|^2}{x^{}},\mathrm{\Delta }X=\frac{_{\mathrm{}}^{\mathrm{}}𝑑x^{}\frac{W^A(x^{})}{x^{}}|u(x^{},t)|^2}{\omega ^2},$$

we get the following equation of motion for the center of soliton

$$X_{tt}+\omega ^2(X+\mathrm{\Delta }X)=0.$$
(39)
Calculation of the parameters $`\mathrm{\Delta }X`$ and $`\omega ^2`$ for our case ($`\alpha _0=0.05,\alpha _1=9,\omega =10`$ and $`b=0.6`$ with $`f(x)=1/\mathrm{cosh}(bx)`$) gives

1) central minimum: $`X_{0}^{}{}_{}{}^{W}=0`$, $`\omega =0.2539`$ ($`T=24.74`$), $`\mathrm{\Delta }X=0`$, $`X^{st}=0`$. Results of numerical simulation: $`\omega =0.251`$ (or $`T=25.0`$), $`\mathrm{\Delta }X=0`$.

2) a lateral minimum: $`X_{0}^{}{}_{}{}^{W}=3.96`$, $`\omega =0.00599`$ ($`T=104.73`$), $`\mathrm{\Delta }X=0.321`$, $`X^{st}=X_{0}^{}{}_{}{}^{W}+\mathrm{\Delta }X=4.18`$. Results of numerical simulation: $`\omega =0.0539`$ (or $`T=116.8`$), $`X^{st}=4.3`$. Here $`X_{0}^{}{}_{}{}^{W}`$ is the position of a given minimum of the potential W(x); $`X_{0}^{}{}_{}{}^{st}`$ is the stationary point (about which oscillations occur); $`\mathrm{\Delta }X`$ is the shift of a stationary point; $`\omega `$ and $`T`$ are the frequency and period of oscillations.

From the present analysis we conclude that the theory of averaging of GP equation over rapidly varying perturbations works well in the case of bounded perturbation potential even when the expansion parameter $`ϵ`$ is sufficient large (in our case $`ϵ=0.9`$). The effective trap potential bifurcates from the one well form to the triple well structure. This can lead to the splitting of single attractive BEC into three parts.

### C A periodic potential

Let $`f(x)=\mathrm{cos}(kx)`$. Then the averaged equation coincide with the unperturbed one but with the another potential:

$$F(x)=\alpha _0\mathrm{cos}(kx)+\frac{ϵ^2k^2}{4}(1\mathrm{cos}(2kx)).$$
(40)
This potential has a more complicated structure and, as a result, the motion of the soliton center has new properties. Let us consider the motion of one soliton in such a potential. The single soliton solution is

$$V(x,t)=2i\eta \text{sech}[2\eta (x\zeta )]\mathrm{exp}[i\frac{\xi z}{\eta }+4i(\eta ^2+\xi ^2)t+i\delta ],$$
(41)
where $`z=2\eta (x\zeta )`$, and $`\zeta =4\xi t`$ for the unperturbed problem. The effective potential acting on the soliton mass center is

$$W_{sol}=\frac{\pi \alpha k}{2\eta \mathrm{sinh}(\frac{\pi k}{2\eta })}\mathrm{cos}(k\zeta )\frac{\pi ϵ^2k^3}{2\eta \mathrm{sinh}(\frac{\pi k}{4\eta })}\mathrm{cos}(2k\zeta ).$$
(42)
The equation of motion is obtained via

$$\zeta _{tt}=\frac{W_s}{\zeta }.$$
(43)
So, for the soliton center $`\zeta `$, we have the following equation of motion

$$\zeta _{tt}=\frac{\pi ϵ^2k^4}{2\eta \mathrm{sinh}(\frac{\pi k}{2\eta })}\mathrm{sin}(2k\zeta )+\frac{\pi k^2\alpha }{2\eta \mathrm{sinh}(\frac{\pi k}{4\eta })}\mathrm{sin}(k\zeta ).$$
(44)
As seen the motion of the soliton center occurs under the action of two harmonic effective perturbation. It represents a close analogy with the stabilization of a pendulum with the rapidly oscillating support point (the Kapitza pendulum). As it is well known , the oscillations of the pendulum at the point $`k\zeta =\pi `$ are unstable. The rapid perturbation stabilize this fixed point. The stabilization condition follows from the effective potential for the soliton in the periodic field (42) and it is given by

$$ϵ^2k^2>\alpha \mathrm{cosh}(\frac{k\pi }{4\eta }).$$
(45)
Fig. 4. Numerical simulations of the full GP equation for a spatial periodic potential $`f(x)=\mathrm{cos}(bx)`$ with the initial soliton position at an unstable point $`X_0=\pi /b`$. The case a) (dashed line) corresponds to infinite motion of the soliton center when $`\alpha _1=0`$ with the fast perturbation turned off. The case b) corresponds to $`\alpha _1=9`$ when the fast perturbation stabilizes the motion of the soliton center. The parameters of perturbation are $`\omega =10`$, $`\alpha _0=0.05`$, as in the case of a bounded potential.

The results of numerical simulation of the full GP equation are presented in Fig.4. The parameters $`\alpha _0`$ and $`\omega `$ are the same as for the above considered bounded potential. Fig. 4a depicts the dynamics of a soliton center in a periodic potential $`\alpha _0\mathrm{cos}(kx)`$ with its initial position $`\zeta =\pi /k`$ corresponding to the unstable point when $`\alpha _1=0`$ (the fast perturbation is turned off). As seen the motion is infinite. Fig. 4b shows the case when $`\alpha _1=9`$ (the fast perturbation is turned on). In this case the motion of the soliton becomes finite and it oscillates with the frequency $`\mathrm{\Omega }_1=0.303`$ whereas the predicted value by Eq. (44) is $`0.308`$. One can see that the agreement with the theory is remarkable.

## IV Conclusion

We have investigated the propagation of bright matter wave soliton in the Bose-Einstein condensate with trap potential rapidly varying in time.

The cases of modulated in time quadratic, bounded and periodic trap potentials have been analyzed. For the repulsive (unstable) trap potential it is shown that there exists a critical value of modulation parameters, when the matter wave soliton is stabilized. Analogous phenomenon of stabilization of unstable fixed points is found for periodic modulations (optical lattices). For the bounded potential it is shown that the effective trap potential can bifurcate from the one well to the triple-well structure; and so, may give rise to the splitting of single attractive BEC into three parts.

## V Acknowledgments

The work is partially supported by Fund of support of the fundamental research of the Uzbek Academy of Sciences(Award 15-02). F.A. is grateful to FAPESP for the support of his work. We are grateful to J.G. Caputo, for the participation at the initial stage of this work and to R. Kraenkel and L. Tomio for useful discussions.


text:
# An Exact Algorithm for Any-flavor Lattice QCD with Kogut-Susskind Fermion

## 1 Introduction

In numerical studies of lattice QCD, advancing simulations including dynamical quarks is the most pressing issue to confirm the validity of QCD and to extract low energy hadronic properties from it. While including dynamical quark effects is still a difficult task, recent developments of computational power and algorithms have enabled dynamical simulations of reasonable scale. Much efforts are thus being spent to accurately compute physical quantities in full QCD simulations .

A large number of these simulations are being made with Wilson-type fermion action using the Hybrid Monte Carlo (HMC) algorithm , which is one of the exact dynamical fermion algorithms. A limitation of the HMC algorithm is that the number of flavors has to be even to express the fermion determinant in terms of pseudo-fermion fields. Recently, however, the polynomial Hybrid Monte Carlo (PHMC) algorithm  has been proposed to simulate the odd number of flavors with Wilson-type fermion as an exact algorithm . The combination of the HMC and PHMC algorithms enables us to simulate the realistic world of 2+1-flavors of quarks with lattice QCD .

The Kogut-Susskind (KS) fermion action has also been widely used in full QCD simulations. For the KS action, the application of the HMC algorithm is restricted to a multiple of four flavors due to the four-fold Dirac fermion content of a single KS fermion. Even if one adopts the usual assumption that the $`1/4`$ power of the KS fermion determinant provides a lattice discretization of a single Dirac fermion determinant, efficient exact algorithms have not been known for two- or single-flavor quark with the KS formalism. For this reason, dynamical KS fermion simulations for two- or three-flavor QCD are made with the $`R`$-algorithm  even today, which is an approximate algorithm.

The $`R`$-algorithm involves a systematic error from a finite step size of the molecular dynamics integrator. Strictly speaking, a careful extrapolation of physical quantities to the zero step size is required. Since this is too computer time consuming, numerical simulations are usually carried out at a finite step size which is chosen so that the systematic error is considered invisible compared to the statistical error. While checks on small lattices are usually made to ensure smallness of the systematic error at least for several quantities, the possibility that other physical quantities are spoiled by the finite step size effects is difficult to exclude. Even such checks become progressively more difficult as smaller quark masses require vastly increasingly computer time. Clearly, an exact and efficient algorithm is desired for dynamical QCD simulations with the KS fermion action not only for two flavors but also for a single flavor of quarks.

In this article, we propose an exact algorithm for KS fermions which is capable of simulating an arbitrary number of flavors. Our algorithm is an application of the PHMC algorithm. It is an extension of the idea of the Rational Hybrid Monte Carlo (RHMC) algorithm  put forward by Horváth, Kennedy, and Sint a few years ago. They briefly described their idea and tested their algorithm together with the PHMC algorithm for two- and four-flavor QCD with the KS fermions on a small size lattices. We shall comment on the difference between their algorithm and ours below. Recently Hasenfratz and Knechtli  also proposed an exact algorithm for KS fermions with hyper-cubic smeared links, which makes use of the polynomial approximation and global update algorithm. The algorithm is considered to be effective for the action for which the HMC type algorithm cannot be applied.

The outline of our algorithms goes as follows. Applying a polynomial approximation to the fractional power of the KS fermion matrix, we rewrite the original partition function to a form suitable for the PHMC algorithm. The resulting effective partition function has two parts; one part is described by an effective action for the polynomial approximation of the fermion action, which can be treated by the HMC algorithm. The other part is the correction term which removes the systematic errors from the polynomial approximation. The correction term can be evaluated by the Kennedy-Kuti  noisy Metropolis test, which has been successfully used in the multi-boson algorithm . With this combination, we can make an exact algorithm for KS fermions. In this work we describe the details of the algorithm, and report results of our numerical test on the applicability of the algorithm to realistic simulations.

Since the polynomial approximation to rewrite the partition function is not unique, there can be several realizations of the PHMC algorithm for the KS fermion. We construct two types of realizations of the polynomial approximation to the PHMC algorithm with $`N_f`$ quark flavors:

Use a polynomial which approximates the $`N_f/8`$ power of the KS fermion matrix which corresponds to $`N_f/2`$ quark flavors. Introducing a single pseudo-fermion field with squaring the fermion matrix, we obtain $`N_f`$ flavors of quarks.

Use an even-order polynomial which approximates the $`N_f/4`$ power of the KS fermion matrix, which corresponds to $`N_f`$ flavors of quarks. To express $`N_f`$ quark flavor with a single pseudo-fermion field, the even-order polynomial is factored into a product of two polynomials using the method by Alexandrou et al. .

We estimate the computational cost of the two algorithms in terms of the number of multiplication by the KS fermion matrix, and find that the algorithm B is roughly a factor two more efficient.

We investigate the property and efficiency of the Chebyshev polynomial approximation as an example of the polynomial approximation. The method to split the even-order polynomial into two polynomials, which is used in case B, is also described. The limitation of our method for splitting the even-order polynomial is discussed.

The Kennedy-Kuti noisy Metropolis test involves a non-trivial factor. Since we take the fractional power of the KS operator, the correction term also includes fractional powers of fermion matrices. In order to evaluate the measure for the noisy Metropolis test acceptance rate, we need to evaluate the fractional power of the correction matrix. To do this, we make use of the Lanczos-based Krylov subspace method by Boriçi , originally proposed for the inverses square root of the squared Hermitian Wilson-Dirac operator in the Neuberger’s overlap fermion. We modify the Boriçi’s algorithm to our purpose.

Finally we investigate the validity and property of the PHMC algorithm (case B) in the case of two flavors of quarks numerically. We first confirm that our algorithm works correctly using a small lattice of a size $`8^3\times 4`$, where a comparison to the $`R`$-algorithm is performed. On an $`8^3\times 16`$ lattice we investigate the mass dependence of the algorithm. We find that for quark masses lighter than $`m_{\mathrm{PS}}/m_\mathrm{V}0.60`$ the polynomial with order larger than $`O(600)`$ is required. Finally, we apply our algorithm to a moderately large lattice of $`16^4`$ with a rather heavy quark mass $`m_{\mathrm{PS}}/m_\mathrm{V}0.69`$ as a prototype for future realistic simulations. On this lattice violation of reversibility and convergence of the Lanczos-based Krylov subspace method for the noisy Metropolis test are investigated. We find satisfactory results; there is no visible reversibility violation, and the Krylov subspace method converges within the limit of double precision arithmetic. The test run on the large size lattice shows reasonable efficiency on the computational time.

The rest of the paper is organized as follows. In section 2, we introduce the lattice QCD partition function with the KS fermions, and rewrite it to a form suitable for the PHMC algorithm. We give two forms for the partition function as described above. The outline of the PHMC algorithm is also shown in this section. In section 3, we describe the Chebyshev polynomial as a specific choice for the polynomial approximation. The error of the polynomial approximation is investigated. The molecular dynamics (MD) step with the polynomial approximation in the PHMC algorithm is explained in section 4. The details of the noisy Metropolis test is given in section 5. We estimate the computational cost of the algorithms in section 6, where we find that the case B algorithm is better. In section 7, we show the test results with the PHMC algorithm. Conclusions are given in the last section.

## 2 Effective Action for PHMC algorithm

The QCD partition function for $`N_f`$ flavors of quarks using the KS fermion is defined by

$$Z=𝒟Udet[D]^{N_f/4}e^{S_g[U]},$$
(2.1)
where $`S_g[U]`$ is the gauge action, $`U_\mu (n)`$ is gauge links, and $`det[D]`$ is the determinant of the KS fermion operator $`D`$. The KS fermion operator $`D`$ is defined by

$$D(n,m)=am\delta _{n,m}+\frac{1}{2}\underset{\mu }{}\eta _\mu (n)[U_\mu (n)\delta _{n+\widehat{\mu },m}U_\mu ^{}(n\widehat{\mu })\delta _{n\widehat{\mu },m},].$$
(2.2)
where $`am`$ is the bare quark mass with lattice unit $`a`$, and $`\eta _\mu (n)`$ is the usual KS fermion phase. In our work we adopt the usual assumption that taking the fourth root of the KS fermion operator represents a lattice discretization of a single Dirac fermion operator in the continuum.

The determinant of the KS fermion operator can be rewritten as

$$det[D]=det\left[\begin{array}{cc}am\mathrm{𝟏}_{ee}& M_{eo}\\ M_{oe}& am\mathrm{𝟏}_{oo}\end{array}\right]=det\left[\begin{array}{cc}\mathrm{𝟏}_{ee}& 0\\ 0& \widehat{D}_{oo}\end{array}\right],$$
(2.3)
where $`\widehat{D}_{oo}=(am)^2\mathrm{𝟏}_{oo}M_{oe}M_{eo}`$ with $`M_{eo}(M_{oe})`$ the hopping matrix from odd site to even site (even site to odd site) defined in Eq. (2.2). Since $`\widehat{D}_{oo}`$ is nothing but the odd part of $`D^{}D`$, the eigenvalues are real and positive semi-definite, which enable us to take the fractional power of the KS fermion operator except for vanishing quark masses. Thus the QCD partition function is reduced to

$$Z=𝒟Udet[\widehat{D}_{oo}]^{N_f/4}e^{S_g[U]}.$$
(2.4)
To apply the PHMC algorithm, we approximate the fractional power of $`\widehat{D}_{oo}`$ by a polynomial of $`\widehat{D}_{oo}`$. We consider two methods for the polynomial approximation. We restrict ourselves to the case that the number of flavors is one or two. Generalization to any integer flavors is achieved by combining the single-flavor and two-flavor cases.

#### Case A

We introduce a polynomial $`P_{N_{\mathrm{poly}}}[x]`$ of order $`N_{\mathrm{poly}}`$, which approximates $`x^{N_f/8}`$ for real and positive (non-zero) $`x`$ as

$$x^{N_f/8}P_{N_{\mathrm{poly}}}[x]=\underset{i=0}{\overset{N_{\mathrm{poly}}}{}}c_ix^i,$$
(2.5)
where $`c_i`$’s are real coefficients. The property that the eigenvalues of the KS fermion operator $`\widehat{D}_{oo}`$ are bounded below by $`(am)^2`$ enables us to substitute $`\widehat{D}_{oo}`$ into Eq. (2.5) to approximate the fractional power of the KS fermion operator. Using this polynomial, we can rewrite the determinant as

$`det\left[\widehat{D}_{oo}\right]^{N_f/4}`$ $`=`$ $`\left({\displaystyle \frac{det\left[\widehat{D}_{oo}\left(P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\right)^{8/N_f}\right]}{det\left[\left(P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\right)^{8/N_f}\right]}}\right)^{N_f/4}`$ (2.6)
$`=`$ $`{\displaystyle \frac{det\left[W[\widehat{D}_{oo}]\right]^{N_f/4}}{det\left[P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\right]^2}},`$
where we introduced

$$W[\widehat{D}_{oo}]=\widehat{D}_{oo}\left(P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\right)^{8/N_f},$$
(2.7)
whose deviation from the identity matrix indicates the residual of the polynomial approximation. We refer to $`W[\widehat{D}_{oo}]`$ as the correction matrix. Note that the exponent of $`P_{N_{\mathrm{poly}}}`$ becomes an integer when $`N_f=2`$ or $`N_f=1`$ as we assumed. Introducing pseudo-fermion fields to the denominator of Eq. (2.6), we obtain

$`det\left[\widehat{D}_{oo}\right]^{N_f/4}`$ $`=`$ $`det\left[W[\widehat{D}_{oo}]\right]^{N_f/4}{\displaystyle 𝒟\varphi _o^{}𝒟\varphi _oe^{S_q[U,\varphi _o^{},\varphi _o]}},`$ (2.8)
$`S_q[U,\varphi _o^{},\varphi _o]`$ $`=`$ $`\left|P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\varphi _o\right|^2.`$ (2.9)
Thus the QCD partition function Eq. (2.4) becomes

$$Z=𝒟U𝒟\varphi _o^{}𝒟\varphi _odet[W[\widehat{D}_{oo}]]^{N_f/4}e^{S_g[U]S_q[U,\varphi _o^{},\varphi _o]}.$$
(2.10)
#### Case B

We define a polynomial $`P_{N_{\mathrm{poly}}}[x]`$ with $`N_{\mathrm{poly}}`$ even to approximate $`x^{N_f/4}`$ for real and positive (non-zero) $`x`$ by

$$x^{N_f/4}P_{N_{\mathrm{poly}}}[x]=\underset{i=0}{\overset{N_{\mathrm{poly}}}{}}c_ix^i.$$
(2.11)
Similarly to the case A, we can rewrite the determinant as

$`det\left[\widehat{D}_{oo}\right]^{N_f/4}`$ $`=`$ $`\left({\displaystyle \frac{det\left[\widehat{D}_{oo}\left(P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\right)^{4/N_f}\right]}{det\left[\left(P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\right)^{4/N_f}\right]}}\right)^{N_f/4}`$ (2.12)
$`=`$ $`{\displaystyle \frac{det\left[W[\widehat{D}_{oo}]\right]^{N_f/4}}{det\left[P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\right]}},`$
where

$$W[\widehat{D}_{oo}]=\widehat{D}_{oo}\left(P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\right)^{4/N_f}.$$
(2.13)
For $`N_{\mathrm{poly}}`$ even, $`P_{N_{\mathrm{poly}}}[x]`$ can be factored into two polynomials as employed in the multi-boson algorithm for single-flavor QCD . Making use of this property, we obtain

$$det\left[\widehat{D}_{oo}\right]^{\frac{N_f}{4}}=\frac{det\left[W[\widehat{D}_{oo}]\right]^{N_f/4}}{\left|det\left[Q_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\right]\right|^2},$$
(2.14)
where $`Q_{N_{\mathrm{poly}}}[x]`$ is defined by

$$\left|Q_{N_{\mathrm{poly}}}[x]\right|^2=P_{N_{\mathrm{poly}}}[x],Q_{N_{\mathrm{poly}}}[x]=\underset{i=0}{\overset{N_{\mathrm{poly}}/2}{}}d_ix^i,$$
(2.15)
with complex coefficients $`d_i`$. The factoring is not unique. In the next section we describe the method for dividing the polynomial. Introducing pseudo-fermion fields to the denominator of Eq. (2.14), we obtain the QCD partition function for the PHMC algorithm in the case B as

$`Z`$ $`=`$ $`{\displaystyle 𝒟U𝒟\varphi _o^{}𝒟\varphi _odet[W[\widehat{D}_{oo}]]^{N_f/4}e^{S_g[U]S_q[U,\varphi _o^{},\varphi _o]}},`$ (2.16)
$`S_q[U,\varphi _o^{},\varphi _o]`$ $`=`$ $`\left|Q_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\varphi _o\right|^2.`$ (2.17)
The PHMC algorithm follows two steps:

1. The HMC step with the effective action Eq. (2.10) for the case A (Eq. (2.16) for the case B).
2. The noisy Metropolis test to remove the systematic error represented by $`W[\widehat{D}_{oo}]`$ in Eq. (2.6) for the case A (Eq. (2.14) for the case B).

The noisy Metropolis test in step (2) has been used in the multi-boson algorithm . The reweighting technique  and stochastic noisy estimator  method can be applied to remove the systematic error. In this paper we employ the noisy Metropolis test to make our algorithm exact .

The original idea by Horváth et al. differs from ours. We separate the action into two pieces; the effective action with polynomial approximation and the correction factor. They do not separate the full action and apply the Kennedy-Kuti noisy estimator for the energy conservation violation $`dH`$ itself to make the algorithm exact, while we apply the method to the correction factor. It is not known which approach is more efficient. We leave a study of this issue to future work as a comparison of the two algorithms is beyond the scope of the present paper.

## 3 Polynomial approximation

We employ the Chebyshev polynomial approximation to $`x^s`$ with a real positive $`x`$, where $`s`$ takes the following values: $`1/2`$, $`1/4`$, $`1/8`$, depending on the choice of case A or B, and on $`N_f`$. We explain the application to the KS fermion operator and investigate the relation among the order of polynomial, the residual of the approximation, and the choice of case A or B.

Several polynomial approximations for $`1/x`$ have been proposed in studies of the multi-boson algorithm. They include Chebyshev , adapted , and least square  polynomials. The choice of the polynomial affects the efficiency, i.e., how much one can decrease the polynomial order so as to make the correction matrix as close to the unit matrix as possible. Since this is not a problem of the simulation algorithm itself, we employ the simple choice of the Chebyshev polynomial approximation for the fractional power of the KS fermion operator.

### 3.1 Chebyshev approximation

The Chebyshev polynomial expansion of $`x^s`$ is

$$x^s=\left[1+(1ϵ)y\right]^s=\underset{k=0}{\overset{\mathrm{}}{}}c_kT_k[y],$$
(3.1)
where $`y=(x1)/(1ϵ)`$, $`ϵ`$ is optimized so as to satisfy $`y[1,1]`$ depending on the support of $`x`$ and $`ϵ(0,1)`$. We restrict the support of exponent $`s`$ to $`s(0,1]`$. 111The fractional power inverse of a matrix is also given by Gegenbauer polynomial expansion . $`T_k[x]`$ is the $`k`$-th order Chebyshev polynomial defined by

$$T_k[x]=\mathrm{cos}(k\mathrm{arccos}(x)).$$
(3.2)
The polynomial has the following recurrence formula:

$$T_k[x]=\alpha _{k1}xT_{k1}[x]+\beta _{k2}T_{k2}[x](k1),$$
(3.3)
where $`\alpha _k=2/(1+\delta _{k,0})`$ ($`k0`$), $`\beta _1=0`$, $`\beta _k=1`$ ($`k0`$), and $`T_0[x]=1`$. The Chebyshev polynomial expansion coefficients $`c_k`$ in Eq. (3.1) are calculated as

$`c_k`$ $`=`$ $`{\displaystyle _{1}^{}{}_{}{}^{1}}[1+(1ϵ)y]^s{\displaystyle \frac{T_k[y]}{\sqrt{1y^2}}}𝑑y`$ (3.4)
$`=`$ $`{\displaystyle \frac{2r^k}{1+\delta _{k,0}}}(1+r^2)^sF(s,s+k;1+k;r^2){\displaystyle \frac{\mathrm{\Gamma }(s+k)}{\mathrm{\Gamma }(s)\mathrm{\Gamma }(1+k)}},`$
where $`r=\left(1+\sqrt{ϵ(2ϵ)}\right)/(1ϵ)`$, $`F(\alpha ,\beta ;\gamma ;z)`$ Gaussian hyper-geometric function, and $`\mathrm{\Gamma }(z)`$ Gamma function. Truncating Eq. (3.1) at order $`N_{\mathrm{poly}}`$, we approximate $`x^s`$ by

$$x^s=\left[1+(1ϵ)y\right]^sP_{N_{\mathrm{poly}}}[x]=\underset{k=0}{\overset{N_{\mathrm{poly}}}{}}c_kT_k[y].$$
(3.5)
The truncation error is bounded by

$$\left|x^sP_{N_{\mathrm{poly}}}[x]\right|\frac{2}{\mathrm{\Gamma }(s)}\left(\frac{1+r^2}{1r^2}\right)^s\frac{(r)^{N_{\mathrm{poly}}+1}}{1+r},$$
(3.6)
where $`r`$ takes a value in the region $`1r0`$. This inequality is not optimal. It demonstrates, however, an exponential decrease of the residual with increasing $`N_{\mathrm{poly}}`$. When $`ϵ1`$, it is expected that $`N_{\mathrm{poly}}\sqrt{ϵ}`$ at a constant residual.

The operator corresponding to $`y`$ in the above relation is given by shifting and changing the normalization of the KS fermion operator:

$$\widehat{D^{}}_{oo}=\mathrm{𝟏}_{oo}(a\lambda )^2\widehat{M}_{oo}=\frac{2}{2(am)^2+(a\mathrm{\Lambda }_{\mathrm{max}})^2}\widehat{D}_{oo},$$
(3.7)
where $`(a\lambda )^2=2(a\mathrm{\Lambda }_{\mathrm{max}})^2/(4(am)^2+2(a\mathrm{\Lambda }_{\mathrm{max}})^2)`$, and $`\widehat{M}_{oo}=\mathrm{𝟏}_{oo}+2M_{oe}M_{eo}/(a\mathrm{\Lambda }_{\mathrm{max}})^2`$. $`a\mathrm{\Lambda }_{\mathrm{max}}`$ is introduced to keep all eigenvalues of $`\widehat{M}_{oo}`$ in the domain $`[1,1]`$. Since the normalization in front of $`\widehat{D}_{oo}`$ in Eq. (3.7) can be absorbed into the normalization of the partition function, we use $`\widehat{D^{}}_{oo}`$ instead of $`\widehat{D}_{oo}`$ and omit the prime symbol from $`\widehat{D^{}}_{oo}`$ in the rest of the paper. For the free case, $`a\mathrm{\Lambda }_{\mathrm{max}}=2`$ is sufficient to satisfy the condition for the eigenvalues of $`\widehat{M}_{oo}`$. In the interacting case, it becomes larger than two. This is seen, for example, in a study for SU(2) lattice gauge theory where the complete eigenvalue distribution has been investigated . With this expression for the KS fermion operator, the polynomial approximation of $`\widehat{D}_{oo}^s`$ becomes

$$\widehat{D}_{oo}^sP_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]=\underset{i=0}{\overset{N_{\mathrm{poly}}}{}}c_kT_k[\widehat{M}_{oo}],$$
(3.8)
where $`c_i`$ is obtained from Eq. (3.4) with $`ϵ=1(a\lambda )^2`$. The exponent $`s`$ is chosen to be $`s=N_f/8`$ for the case A and $`s=N_f/4`$ for the case B.

For the case B, we have to solve Eq. (2.15) to obtain the half-order polynomial $`Q_{N_{\mathrm{poly}}}`$. Here we choose to construct $`Q_{N_{\mathrm{poly}}}`$ so as to have the form of the Chebyshev polynomial expansion as Eq. (3.8). Since this problem is rather complicated, we postpone the discussion to the subsection at the end of this section, and proceed assuming that the polynomial $`Q_{N_{\mathrm{poly}}}`$ and its coefficients $`d_i`$ are already given.

Given the Chebyshev polynomial expansion coefficients $`c_m`$ ($`d_m`$), we can evaluate the polynomial $`P_{N_{\mathrm{poly}}}[x]`$ ($`Q_{N_{\mathrm{poly}}}[x]`$) using the Clenshaw’s recurrence formula (for example, see ). When $`x`$ is the KS fermion operator $`\widehat{D}_{oo}`$, the multiplication of the operator $`P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]`$ on a vector $`v_o`$ is carried out by the following three step recurrence formula:

$$y_o^{(k)}=\alpha _k(\widehat{M}_{oo})y_o^{(k+1)}+\beta _ky_o^{(k+2)}+c_kv_o,$$
(3.9)
where $`y^{(k)}`$ is a working vector labeled $`k`$, $`\alpha _k`$ and $`\beta _k`$ are given in Eq. (3.3), and $`c_k`$ is the Chebyshev polynomial expansion coefficient. Solving for $`y_o^{(k)}`$ with this equation from $`k=N_{\mathrm{poly}}`$ to $`k=0`$ with the initial condition $`y_o^{(N_{\mathrm{poly}})}=y_o^{(N_{\mathrm{poly}}+1)}=0`$, we obtain

$$P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]v_o=y_o^{(0)}.$$
(3.10)
For $`Q_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]v_o`$, $`d_k`$ is used instead of $`c_k`$ and $`k`$ runs from $`N_{\mathrm{poly}}/2`$ to $`0`$. Note that we need not store all working vectors $`y_o^{(k)}`$; only two working vectors are required for the computation. This method can be applied to any matrix polynomial which has the same structure for the recurrence relation as Eq. (3.3). The computational cost to calculate $`P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]v_o`$ ($`Q_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]v_o`$) is $`N_{\mathrm{poly}}`$ ($`N_{\mathrm{poly}}/2`$) by means of the number of matrix-vector multiplication.

Now we discuss the quality of the polynomial approximation. We evaluate the polynomial approximation residual using

$$R(x)=\left|x(P_{N_{\mathrm{poly}}}[x])^{1/s}1\right|.$$
(3.11)
Figure 1 shows the residual of the approximation as a function of $`y=(x1)/(1ϵ)`$. The calculation is performed using Clenshaw’s recurrence in double precision arithmetic. In order to compare the approximation at the fixed computational cost, the number of matrix-vector multiplication, in calculating the correction matrix irrespective of the choice of $`s`$ (case A or B and $`N_f`$), we keep $`N_{\mathrm{poly}}/s`$ constant ($`N_{\mathrm{poly}}/s=800`$) as an example. We observe that the approximation becomes worse for smaller $`s`$ ($`s=1`$ case reaches the limit of double precision arithmetic). We also investigate the $`N_{\mathrm{poly}}/s`$ dependence of the polynomial approximation. In figure 2 we plot the $`N_{\mathrm{poly}}/s`$ dependence of the residual defined by

$$\sqrt{R^2}=\sqrt{_{1}^{}{}_{}{}^{1}𝑑yR(x(y))^2}.$$
(3.12)
Clear exponential decay is observed. The dependence of the slope on the choice of $`s`$ indicates that the computational cost increases with decreasing $`s`$. In our case, defining the polynomial order by $`N_{\mathrm{poly}}^A`$ and $`N_{\mathrm{poly}}^B`$ for the case A and B, respectively, it is expected that $`N_{\mathrm{poly}}^A/(N_f/8)2N_{\mathrm{poly}}^B/(N_f/4)`$ holds when we require that the integrated residual takes a similar value for the two cases. We thus find $`N_{\mathrm{poly}}^AN_{\mathrm{poly}}^B`$ at the same value of the polynomial residual.

### 3.2 Determination of $`Q_{N_{\mathrm{poly}}}`$ in case B

Here we discuss how to solve Eq. (2.15) to obtain the half-order polynomial $`Q_{N_{\mathrm{poly}}}`$. In our work $`Q_{N_{\mathrm{poly}}}`$ should have the Chebyshev polynomial expansion form of Eq. (3.8). A simple procedure to obtain $`Q_{N_{\mathrm{poly}}}`$ is as follows: (i) express Eq. (3.5) as a polynomial of $`y`$ instead of the expansion of the Chebyshev polynomial, (ii) split the polynomial into the product of two polynomials like Eq. (2.15) by means of usual polynomial expansion, and (iii) recover the Chebyshev polynomial expansion for the half-order polynomial. However, this method has a numerical problem in step (iii). In order to make this point clear, and to present an alternative procedure, let us elaborate on the procedure above.

In the step (i), we expand the Chebyshev polynomial to write

$$P_{N_{\mathrm{poly}}}[x]=\underset{k=0}{\overset{N_{\mathrm{poly}}}{}}c_kT_k[y]=\underset{k=0}{\overset{N_{\mathrm{poly}}}{}}c_k^{}y^k,$$
(3.13)
where $`c_k^{}`$ can be obtained from the original $`c_k`$ using the appropriate recurrence formula (see ex. ). In the step (ii), finding the roots of the polynomial $`_{k=0}^{N_{\mathrm{poly}}}c_k^{}y^k=0`$, we obtain the product representation:

$$\underset{k=0}{\overset{N_{\mathrm{poly}}}{}}c_k^{}y^k=c_{N_{\mathrm{poly}}}^{}\underset{k=1}{\overset{N_{\mathrm{poly}}}{}}(yz_k),$$
(3.14)
where $`z_k`$’s are the roots of the polynomial. Because $`N_{\mathrm{poly}}`$ is even and the coefficients $`c_k^{}`$’s are real, the root $`z_k`$ pairs with its complex conjugate $`z_k^{}=z_k^{}`$. Thus we can split the polynomial to the product of two polynomials :

$`{\displaystyle \underset{k=0}{\overset{N_{\mathrm{poly}}}{}}}c_k^{}y^k`$ $`=`$ $`c_{N_{\mathrm{poly}}}^{}{\displaystyle \underset{k=1}{\overset{N_{\mathrm{poly}}/2}{}}}(yz_{j(k)})(yz_{j(k)}^{}),`$ (3.15)
$`=`$ $`\left|\left(c_{N_{\mathrm{poly}}}^{}\right)^{1/2}{\displaystyle \underset{k=1}{\overset{N_{\mathrm{poly}}/2}{}}}(yz_{j(k)})\right|^2,`$
where $`j(k)`$ is a reordering index to represent the pairing condition above. The explicit form of the reordering is not unique depending on how to distribute the complex pairs in the monomials, and several methods have been proposed .

In the third step (iii), we calculate the half-order polynomial according to

$`\left(c_{N_{\mathrm{poly}}}^{}\right)^{1/2}{\displaystyle \underset{k=1}{\overset{N_{\mathrm{poly}}/2}{}}}(yz_{j(k)})`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{N_{\mathrm{poly}}/2}{}}}d_k^{}y^k,`$ (3.16)
$`=`$ $`{\displaystyle \underset{k=0}{\overset{N_{\mathrm{poly}}/2}{}}}d_kT_k[y],`$
where $`d_k^{}`$’s are obtained from $`z_k`$ and $`c_{N_{\mathrm{poly}}}^{}`$ by expanding the product representation, and $`d_k`$’s are extracted from $`d_k^{}`$’s using an appropriate reverse recurrence formula  (i.e. the relation opposite to Eq. (3.13)). In this way, we could derive $`d_k`$’s from $`c_k`$’s so as to satisfy

$$P_{N_{\mathrm{poly}}}[x]=\underset{k=0}{\overset{N_{\mathrm{poly}}}{}}c_kT_k[y]=\left|\underset{k=0}{\overset{N_{\mathrm{poly}}/2}{}}d_kT_k[y]\right|^2=\left|Q_{N_{\mathrm{poly}}}[x]\right|^2.$$
(3.17)
Unfortunately, we find that the reverse recurrence to extract the Chebyshev polynomial expansion coefficients $`d_k`$ from the usual polynomial coefficients $`d_k^{}`$’s is numerically unstable . Therefore we decided to directly solve Eq. (3.17) with respect to $`d_k`$.

Using the addition relation, $`T_k[y]T_l[y]=(T_{k+l}[y]+T_{|kl|}[y])/2`$, to Eq. (3.17), we extract the following second-order simultaneous equations,

$$f_k(\{d\},\{d^{}\})=c_k,(k=0,\mathrm{},N_{\mathrm{poly}}),$$
(3.18)
where $`f_k(\{d\},\{d^{}\})`$ depends on the sets $`\{d\}=\{d_1,d_2,\mathrm{},d_{N_{\mathrm{poly}}/2}\}`$ and $`\{d^{}\}=\{d_1^{},d_2^{},\mathrm{},d_{N_{\mathrm{poly}}/2}^{}\}`$. We do not write down the explicit form of $`f_k(\{d\},\{d^{}\})`$ here because of its length and complicated form. The solution to these equations is not unique. This corresponds to the reordering ambiguity in the previous case of splitting the product representation.

We solve numerically the simultaneous equation Eq. (3.18) using Mathematica with desired accuracy starting from an initial choice of $`\{d\}`$ and $`\{d^{}\}`$. The accuracy of the solution is examined by numerically evaluating the relation Eq. (3.17). Figure 3 shows the polynomial coefficients $`d_i`$ for $`Q_{N_{\mathrm{poly}}}[x]`$ derived by the direct method using Eq. (3.18). The accuracy of the solution stays at satisfactory level within double precision arithmetic as observed in Figure 4, where the polynomials are evaluated with Clenshaw’s recurrence formula in double precision arithmetic.

A practical limitation with our direct method is that it is rather slow. On a Linux PC with 1 GHz Pentium III CPU, we could find the coefficients only for $`N_{\mathrm{poly}}<800`$ within a tolerable computational time. A more efficient methods to solve Eq. (2.15) is desired.

## 4 Calculation of the polynomial pseudo-fermion force in the Hybrid Monte Carlo algorithm

We apply the usual HMC algorithm to the partition function Eq. (2.10) (case A) or Eq. (2.16) (case B), introducing a fictitious time and canonical momenta $`P_\mu (n)`$ to the link variables $`U_\mu (n)`$. A nontrivial task is the calculation of the molecular dynamics (MD) force from the quark action expressed by the polynomial approximation. We utilize the Clenshaw’s recurrence formula for this purpose. Here we first discuss the variation of the polynomial of a matrix $`A`$ for the general case, and then describe the force calculation in the HMC algorithm.

Let $`P_N[A]`$ be a matrix polynomial of A with order $`N`$ described by

$$P_N[A]=\underset{i=0}{\overset{N}{}}c_i\mathrm{\Phi }_i[A],$$
(4.1)
where $`\mathrm{\Phi }_i[A]`$ is defined to have the following recurrence relation:

$$\mathrm{\Phi }_i[A]=\alpha _{i1}A\mathrm{\Phi }_{i1}[A]+\beta _{i2}\mathrm{\Phi }_{i2}[A].$$
(4.2)
In most cases $`\mathrm{\Phi }_0[A]`$ is a constant and set to be unity. As in Eq. (3.9), $`P_N[A]`$ is evaluated by the Clenshaw’s recurrence formula:

$$Y^{(k)}=\alpha _kAY^{(k+1)}+\beta _kY^{(k+2)}+c_k\mathrm{𝟏},$$
(4.3)
where $`Y^{(k)}`$’s are working matrices, $`Y^{(N+1)}=Y^{(N)}=0`$, $`k`$ runs from $`N`$ to $`0`$, and $`P_N[A]=Y^{(0)}`$. We take the variation of Eq. (4.1) with respect to $`A`$, and denote it by $`\delta P_N[A]`$:

$$\delta P_N[A]=\underset{i=0}{\overset{N}{}}\delta \mathrm{\Phi }_i[A]c_i.$$
(4.4)
The variation $`\delta \mathrm{\Phi }_i[A]`$ also has the recurrence formula obtained by differentiating Eq. (4.2):

$$\delta \mathrm{\Phi }_i[A]=\alpha _{i1}\delta \mathrm{\Phi }_{i1}[A]A+\beta _{i2}\delta \mathrm{\Phi }_{i2}[A]+\alpha _{i1}\mathrm{\Phi }_{i1}[A]\delta A.$$
(4.5)
Substituting Eq. (4.3) and Eq. (4.5) into Eq. (4.4), we have

$$\delta P_N[A]=\underset{i=1}{\overset{N}{}}\alpha _{i1}\mathrm{\Phi }_{i1}[A]\delta AY^{(i)},$$
(4.6)
where $`\delta \mathrm{\Phi }_0[A]=0`$ and a similar technique to the Clenshaw’s formula is used.

Applying Eq. (4.6) to our problem, we can evaluate the variation of the pseudo-fermion action Eq. (2.9) with respect to the infinitesimal change of link variables as

$`\delta S_q`$ $`=`$ $`\delta \left|P_{N_{\mathrm{poly}}}[\widehat{M}_{oo}]\varphi _o\right|^2`$ (4.7)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N_{\mathrm{poly}}}{}}}\left(\alpha _{i1}T_{i1}[\widehat{M}_{oo}]y_o^{(0)}\right)^{}\left(\delta \widehat{M}_{oo}\right)y_o^{(i)}+\text{h.c.},`$
where h.c. means the Hermitian conjugate of the preceding term, and $`y_o^{(i)}`$ satisfies Eq. (3.9) with $`v_o=\varphi _o`$. We used the Hermiticity of $`\alpha _k`$ and $`P_{N_{\mathrm{poly}}}[\widehat{M}_{oo}]`$, and $`y_o^{(0)}=P_{N_{\mathrm{poly}}}[\widehat{M}_{oo}]\varphi _o`$ was applied in the last line. A more convenient form of $`\delta S_q`$ for practical calculations is given by

$`\delta S_q={\displaystyle \frac{2}{(a\mathrm{\Lambda }_{\mathrm{max}})^2}}{\displaystyle \underset{i=1}{\overset{N_{\mathrm{poly}}}{}}}\alpha _{i1}X_{}^{(i)}{}_{}{}^{}\delta MZ^{(i)}+\text{h.c.},`$ (4.8)
where $`X^{(i)}`$, $`Y^{(i)}`$, and $`M`$ are defined as

$`X^{(i)}`$ $`=`$ $`\left(\begin{array}{c}M_{eo}x_o^{(i)}\\ x_o^{(i)}\end{array}\right),`$ (4.11)
$`Z^{(i)}`$ $`=`$ $`\left(\begin{array}{c}M_{eo}y_o^{(i)}\\ y_o^{(i)}\end{array}\right),`$ (4.14)
$`\delta M`$ $`=`$ $`\left(\begin{array}{cc}0& \delta M_{eo}\\ \delta M_{oe}& 0\end{array}\right),`$ (4.17)
with

$`y_o^{(i)}`$ $`=`$ $`\alpha _i(\widehat{M}_{oo})y_o^{(i+1)}+\beta _iy_o^{(i+2)}+c_i\varphi _o(i=N_{\mathrm{poly}},\mathrm{},0),`$
$`x_o^{(i)}`$ $`=`$ $`\alpha _{i2}(\widehat{M}_{oo})x_o^{(i1)}+\beta _{i3}x_o^{(i2)}(i=2,\mathrm{},N_{\mathrm{poly}}),`$
$`x_o^{(1)}`$ $`=`$ $`y_o^{(0)}.`$ (4.18)
In the actual calculation, we first calculate $`y_o^{(i)}`$ from $`i=N_{\mathrm{poly}}`$ to $`i=0`$ and store $`Z^{(i)}`$ on memory. We then sum up each of the force contribution $`X_{}^{(i)}{}_{}{}^{}\delta MZ^{(i)}`$ evaluating $`X^{(i)}`$ from $`i=1`$ to $`i=N_{\mathrm{poly}}`$. A similar method, which has a more complicated form, has been obtained for the HMC algorithm with the overlap fermions by C. Liu . Since the even component of $`X^{(i)}`$ and $`Z^{(i)}`$ appear as a byproduct of the multiplication of $`\widehat{M}_{oo}`$ in the recurrence equation, we do not need extra calculations for the even components. The explicit form of the force contribution is obtained from these equations as usual.

We have described the force calculation for the case A as an example. The force calculation for the case B is almost identical except for the replacement $`N_{\mathrm{poly}}N_{\mathrm{poly}}/2`$ and $`c_id_i`$.

Eqs. (4.8), (4.17), and (4.18) do not contain an iterative procedure. This leads us to expect that the reversibility violation of the MD evolution is smaller than in the usual HMC algorithm with four-flavor KS fermions in which an iterative solver such as the Conjugate Gradient (CG) method is used to invert the KS fermion operator. Our implementation, however, still involves the possibility that round-off errors grow to violate the reversibility in the summation of $`X_{}^{(i)}{}_{}{}^{}\delta MZ^{(i)}`$ from $`i=1`$ to $`i=N_{\mathrm{poly}}`$. We investigate this issue in section 7 on a realistic size lattice.

The external field $`\varphi _o`$ is generated at the beginning of every MD trajectory according to the pseudo-fermion action Eq. (2.9) (case A) (Eq. (2.17) (case B)) using the global heat-bath method. This is achieved by solving the following equations with respect to $`\varphi _o`$:

$`P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\varphi _o`$ $`=`$ $`\chi _o\text{ (case A)},`$ (4.19)
$`Q_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]\varphi _o`$ $`=`$ $`\chi _o\text{ (case B)},`$ (4.20)
where $`\chi _o`$ is a Gaussian noise vector with unit variance. Using the identity

$`\varphi _o`$ $`=`$ $`(P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}])^1\chi _o`$ (4.21)
$`=`$ $`\widehat{D}_{oo}(P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}])^{(8/N_f1)}W[\widehat{D}_{oo}]^1\chi _o\text{ (case A)},`$
$`\varphi _o`$ $`=`$ $`(Q_{N_{\mathrm{poly}}}[\widehat{D}_{oo}])^1\chi _o`$ (4.22)
$`=`$ $`\widehat{D}_{oo}(Q_{N_{\mathrm{poly}}}[\widehat{D}_{oo}])^{}(P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}])^{(4/N_f1)}W[\widehat{D}_{oo}]^1\chi _o\text{ (case B)},`$
we invert the coefficient matrix $`P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]`$ (case A) ($`Q_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]`$ (case B)) using the CG solver to $`W[\widehat{D}_{oo}]`$.

## 5 Noisy Metropolis test

In order to make algorithm exact, we have to take into account the correction term $`det[W[\widehat{D}_{oo}]]^{N_f/4}`$. This is achieved by a noisy Kennedy-Kuti  Metropolis test. This method has been used in the multi-boson algorithm . In this section we explain the details of the noisy Metropolis test.

### 5.1 Definition

When a trial configuration $`U^{}`$ which is generated by the molecular dynamics step in the HMC part of the algorithm is accepted at the HMC Metropolis test, we make the noisy Metropolis test for the correction factor $`det[W[\widehat{D}_{oo}]]^{N_f/4}`$. The acceptance probability of the trial configuration $`U^{}`$ from an initial configuration $`U`$ is defined by

$$P_{\mathrm{𝑐𝑜𝑟𝑟}}[UU^{}]=\mathrm{min}[1,e^{dS[U,U^{}]}],$$
(5.1)
with

$$dS[U,U^{}]=\left|W[\widehat{D}_{oo}^{}]^{N_f/8}W[\widehat{D}_{oo}]^{N_f/8}\eta _o\right|^2|\eta _o|^2,$$
(5.2)
where $`\eta _o`$ is a Gaussian random vector with unit variance. This probability satisfies the detailed balance relation . The factor $`W[\widehat{D}_{oo}^{}]`$ is calculated on the trial configuration $`U^{}`$, while $`W[\widehat{D}_{oo}]`$ on the initial configuration $`U`$.

Eq. (5.2) can be modified to

$`dS[U,U^{}]`$ $`=`$ $`\zeta _o^{}W[\widehat{D}_{oo}^{}]^{N_f/4}\zeta _o|\eta _o|^2,`$
$`\zeta _o`$ $`=`$ $`W[\widehat{D}_{oo}]^{N_f/8}\eta _o.`$ (5.3)
We employ Eq. (5.3) for $`dS`$ in the present work. It is not known at present which of the expressions, Eq. (5.2) or (5.3), is more useful for calculating $`dS`$ in respect of computational efficiency and accuracy, which is left for future studies.

In order to evaluate Eq. (5.3), we need to calculate the fractional power of the matrix $`W[\widehat{D}_{oo}]`$ (or $`W[\widehat{D}_{oo}^{}]`$). In Ref.  we used a Taylor expansion method for the (inverse-)square root of the correction matrix with the $`O(a)`$-improved Wilson fermions. In order to suppress the truncation error of the Taylor expansion we explicitly calculate and monitor the residual for the Taylor expansion for which we need an extra computational cost. The direct application of this method to Eq. (5.2) requires more computational overhead compared to the Wilson case, since the residuals contains several application of $`W[\widehat{D}_{oo}]^{N_f/8}`$ and $`W[\widehat{D}_{oo}]^{N_f/4}`$, (e.g. for $`N_f=2`$ case, 8 times, and 4 times to define the residuals, respectively). Instead of this method we employ the Krylov subspace method to avoid the explicit calculation of the residuals as described in the following.

Here we roughly estimate the order of polynomial $`N_{\mathrm{poly}}`$ required to achieve a given acceptance rate at a volume $`V`$ and $`a\mathrm{\Lambda }_{\mathrm{max}}`$. The acceptance rate Eq. (5.3) behaves as $`dS=\text{[constant]}\times V\times (r)^{N_{\mathrm{poly}}+1}`$, where $`r`$ is defined in Eq. (3.4). To keep $`dS<\delta `$ with a small constant $`\delta `$ that maintains sufficient acceptance rate, we need

$$N_{\mathrm{poly}}>\frac{a\mathrm{\Lambda }_{\mathrm{max}}}{2am}\left(\mathrm{ln}V\mathrm{ln}\delta +\text{[constant]}\right),$$
(5.4)
where we used $`r1+2(am)/(a\mathrm{\Lambda }_{\mathrm{max}})`$ and $`ama\mathrm{\Lambda }_{\mathrm{max}}`$ from the definition of $`r`$ in Eq. (3.7). Therefore we need to increase $`N_{\mathrm{poly}}`$ linearly as increasing the condition number $`(a\mathrm{\Lambda }_{\mathrm{max}})/(am)`$, and logarithmically with volume $`V`$. A similar discussion on the required $`N_{\mathrm{poly}}`$ can be found in Refs.  in the literature of the multi-boson algorithm.

### 5.2 The Krylov subspace method

Since the matrix is Hermitian and positive definite with the KS fermion, and already well preconditioned, we can take the fractional power with the Krylov subspace method with better efficiency . A Lanczos-based Krylov subspace method was developed by Boriçi , which was utilized for the calculation of inverse square root of squared Hermitian Wilson-Dirac operator in the Neuberger’s overlap operator. The method is an application of the Krylov subspace method to calculate functions of large sparse matrices . The Lanczos-based Krylov subspace method enables us to define a kind of residual without explicit residual calculation. We apply this method to evaluate the fractional power of the correction matrix.

Consider a matrix function multiplied on a vector, $`f(A)b`$, with an $`n\times n`$ Hermitian matrix $`A`$ and an $`n`$ component vector $`b`$ in general. The $`k`$-dimensional orthogonal basis $`Q_k=(q_1,\mathrm{},q_k)`$, which spans the Krylov subspace with respect to the matrix $`A`$, is obtained by the Lanczos algorithm. This basis satisfies the following condition:

$$AQ_k=Q_kT_k+\beta _kq_{k+1}(e_k^{(k)})^T,q_1=\rho _1b,\rho _1=1/|b|,$$
(5.5)
where $`T_k`$ is a $`k\times k`$ tridiagonal matrix whose diagonal and sub-diagonal parts are $`\alpha _1,\alpha _2,\mathrm{},\alpha _k`$ and $`\beta _1,\beta _2,\mathrm{},\beta _{k1}`$, respectively. $`e_m^{(n)}`$ is the unit vector in the $`m`$-th direction in $`n`$-dimension. The superscript $`T`$ means transpose.

If $`\beta _k`$ is sufficiently small after $`k`$ iterations of the Lanczos method, the matrix function $`f(A)`$ acting on the vector $`b`$ is approximated by

$$f(A)bQ_kf(T_k)e_1^{(k)}|b|.$$
(5.6)
In our case, $`A=W[\widehat{D}_{oo}]`$, $`f(x)=x^{N_f/4}`$, and $`b=\eta _o`$, or $`A=W_{oo}^{}`$, $`f(x)=x^{N_f/8}`$, and $`b=\zeta _o`$. The dimension $`k`$ is much smaller than that of $`A`$ when $`A`$ is close to unity. Thus, the calculation of $`f(A)`$ is reduced to that of $`f(T_k)`$ with smaller computational cost.

Our algorithm is almost identical to that of Ref. . We show the algorithm to calculate the matrix function $`f(A)b`$ with the Lanczos-based method in Algorithm 1. We compute $`f(T_k)`$ by diagonalizing $`T_k`$ with LAPACK subroutines. If the algorithm stops after $`k`$ iterations, we have an approximate solution to the matrix function $`f(A)`$ given by $`x_kf(A)b`$.

In our algorithm, a large cancellation error can occur in the Gram-Schmidt orthogonalization step because our correction matrix is well conditioned $`W[\widehat{D}_{oo}]1`$. We therefore implement a full reorthogonalization step in our algorithm to keep the orthogonality among the Lanczos vectors $`q_i`$.

Our algorithm requires $`k`$ working vectors to store the Lanczos vectors. However, we already employ $`N_{\mathrm{poly}}`$ working vectors to calculate the quark force in the HMC step, which we can reuse for the Lanczos vectors. A possible problem that the dimension of Krylov subspace, $`k`$, exceeds the number of working vectors, $`N_{\mathrm{poly}}`$, do not arise in practice for large $`N_{\mathrm{poly}}`$ because large $`N_{\mathrm{poly}}`$ means that $`W`$ is very close to unity so that the Lanczos algorithm stops at earlier steps.

We can consider various types of stopping condition. For example, $`\mathrm{𝑒𝑟𝑟}_1`$ is based on the residual in CG algorithm for the calculation of $`A^1b`$ , and $`\mathrm{𝑒𝑟𝑟}_2`$ on a comparison of the successive approximation $`x_i`$ as described in Algorithm 1. In order to avoid explicit residual calculation, these stopping criteria should ensure that the residual decreases to sufficient level during the iteration. For the CG-based stopping condition, which was originally introduced by Boriçi , the analytical relation between the CG-based stopping condition and the truncation error of the Lanczos iteration is discussed by van den Eshof et al. , where it is proved that the CG-based stopping condition for the inverse square root in the overlap operator is a safe stopping condition.

We cannot directly apply their proof to our case because the exponent of the correction matrix is not limited to $`1/2`$. We employ the CG-based stopping condition $`\mathrm{𝑒𝑟𝑟}_1`$ in the noisy Metropolis test, and numerically verify the validity of this choice by observing the convergence behavior of the residuals and $`dS`$. This analysis is described in section 7. We leave the mathematical proof whether the CG stopping condition is safe or not for our case for future studies.

## 6 Cost estimate

The computational cost of our algorithm is measured by the number of multiplication with $`\widehat{M}_{oo}`$ for traversing a single unit of trajectory. The number of multiplication is separately estimated for the HMC part and the noisy Metropolis part. We define the algorithmic parameters as follows:

* the number of MD step : $`N_{\mathrm{MD}}^A`$
* the number of iteration of the CG algorithm : $`N_{\mathrm{CG}}^A`$
* the order of polynomial : $`N_{\mathrm{poly}}^A`$

where the superscript $`A`$ refers to the case A algorithm, which should be replaced with $`B`$ for the case B algorithm. We use the single leapfrog scheme for the MD integrator.

#### Cost of HMC part

The computational cost of the HMC part of the algorithm is divided into three pieces; calculation of the MD force with the polynomial pseudo-fermion, the generation of pseudo-fermion field according to the polynomial action, and the calculation of the total Hamiltonian for the HMC Metropolis test.

From Eqs. (4.8), (4.17), and (4.18) in section 4, the cost of the force calculation in the HMC algorithm is estimated as

$$N_{\mathrm{MD}}^A\times (2N_{\mathrm{poly}}^A1)\text{ ( case A)},$$

$$N_{\mathrm{MD}}^B\times (N_{\mathrm{poly}}^B1)\text{ ( case B)}.$$

From Eqs. (4.21) and (4.22), the computational cost to generate the pseudo-fermion field is estimated as

$$((8/N_f)\times N_{\mathrm{poly}}^A+1)\times N_{\mathrm{CG}}^A+(8/N_f1)\times N_{\mathrm{poly}}^A+1\text{ (case A)},$$

$`((4/N_f)\times N_{\mathrm{poly}}^B+1)\times N_{\mathrm{CG}}^B+(4/N_f1)\times N_{\mathrm{poly}}^B+N_{\mathrm{poly}}^B/2+1`$
$`\text{(case B)}.`$
The computational cost of the Hamiltonian comes from the calculation of the pseudo-fermion action at the end of the MD step. The number of multiplication of $`\widehat{M}_{oo}`$ is estimated as $`N_{\mathrm{poly}}^A`$ for the case A and $`N_{\mathrm{poly}}^B/2`$ for the case B.

Summarizing, the computational cost of the HMC part of our algorithm is given by

$`N_{\mathrm{HMC}\mathrm{cost}}^A`$ $`=`$ $`(2N_{\mathrm{poly}}^A1)\times N_{\mathrm{MD}}^A+((8/N_f)\times N_{\mathrm{poly}}^A+1)\times N_{\mathrm{CG}}^A`$ (6.1)
$`+(8/N_f1)\times N_{\mathrm{poly}}^A+1,`$
for the case A, and

$`N_{\mathrm{HMC}\mathrm{cost}}^B`$ $`=`$ $`(N_{\mathrm{poly}}^B1)\times N_{\mathrm{MD}}^B+((4/N_f)\times N_{\mathrm{poly}}^B+1)\times N_{\mathrm{CG}}^B`$ (6.2)
$`+(4/N_f1)\times N_{\mathrm{poly}}^B+N_{\mathrm{poly}}^B/2+1,`$
for the case B.

We observed in section 3 that $`N_{\mathrm{poly}}^AN_{\mathrm{poly}}^B`$ at a comparable value of the polynomial approximation residual Eq. (3.12). Assuming $`N_{\mathrm{MD}}^AN_{\mathrm{MD}}^B`$ and $`N_{\mathrm{CG}}^AN_{\mathrm{CG}}^B`$, we find $`N_{\mathrm{HMC}\mathrm{cost}}^A2N_{\mathrm{HMC}\mathrm{cost}}^B`$. We conclude that the cost of HMC part of the algorithm is twice better for the case B than for the case A.

#### Cost of noisy Metropolis part

Here we estimate the cost of the noisy Metropolis test, $`N_{\mathrm{NMP}\mathrm{cost}}^A`$ and $`N_{\mathrm{NMP}\mathrm{cost}}^B`$, for both cases. The cost arises from Eq. (5.3), and is estimated as

$$N_{\mathrm{NMP}\mathrm{cost}}^A=((8/N_f)\times N_{\mathrm{poly}}^A+1)\times N_{\mathrm{CG}}^A\times 2,$$
(6.3)
for the case A, and

$$N_{\mathrm{NMP}\mathrm{cost}}^B=((4/N_f)\times N_{\mathrm{poly}}^B+1)\times N_{\mathrm{CG}}^B\times 2,$$
(6.4)
for the case B. The factor $`2`$ arises since we need to call twice the Lanczos-based algorithm to calculate $`W[\widehat{D}_{oo}]^{N_f/8}\eta _o`$ and $`W[\widehat{D}_{oo}^{}]^{N_f/4}\zeta _o`$. Here we used $`N_{\mathrm{CG}}^A`$ ($`N_{\mathrm{CG}}^B`$) as the number of Lanczos iteration. This is because the number of Lanczos iteration is expected to be almost identical to that of CG iteration to generate the pseudo-fermion field, when we employ the CG based stopping criterion. We expect $`N_{\mathrm{CG}}^AN_{\mathrm{CG}}^B`$ and $`N_{\mathrm{poly}}^AN_{\mathrm{poly}}^B`$ as before. Then we find $`N_{\mathrm{cost}\mathrm{A}}^{\mathrm{NMP}}2N_{\mathrm{cost}\mathrm{B}}^{\mathrm{NMP}}`$.

#### Total cost

Combining the result on the computational cost for the HMC step and that for the noisy Metropolis test described above, we find that the cost for the case A algorithm is larger than that for the case B algorithm by a factor two when the approximated polynomials for the two algorithms have the same residual on the correction matrix. In our numerical test described in the next section, we employ the case B algorithm.

## 7 Numerical tests

We test our algorithm (case B) with the Chebyshev polynomial on three lattices in the two-flavor case. The simulation program is written in the optimized FORTRAN90 on SR8000 model F1 at KEK. Double precision arithmetic is applied to the whole numerical operations. The lattice volume, gauge coupling, and quark masses we employed are shown in Table 1. The “small size” lattice is used to investigate the basic property of the algorithm. We show the $`N_{\mathrm{poly}}`$ dependence of $`dS`$ and the averaged acceptance rate of the noisy Metropolis test $`P_{\mathrm{𝑐𝑜𝑟𝑟}}`$. The $`N_{\mathrm{poly}}`$ dependence of averaged plaquette $`P`$ is also presented on this lattice, and a comparison to results from the $`R`$-algorithm is also shown. Using the “middle size” lattice, we compare the $`N_{\mathrm{poly}}`$ dependence of $`dS`$ for different quark masses. We employ the “large size” lattice parameter to see the applicability of our algorithm to realistic simulations, where we check the reversibility of the MD evolution and the validity of the CG-based stopping criterion for the Lanczos algorithm in the noisy Metropolis test. A comparison of $`P`$ to the $`R`$-algorithm is also made.

The unit trajectory length is chosen to be $`1`$ in the following. We employ the single leapfrog scheme for the MD evolution, and call the number of MD step as $`N_{\mathrm{MD}}`$. The parameter $`a\mathrm{\Lambda }_{\mathrm{max}}`$ is roughly optimized during the thermalization period for each lattice parameter.

### 7.1 Results on the small size lattice

On the small size lattice, $`a\mathrm{\Lambda }_{\mathrm{max}}`$ is chosen to be $`2.37`$. Figure 5 shows the $`N_{\mathrm{poly}}`$ dependence of $`dS`$ for this lattice, where the number of MD step is $`N_{\mathrm{MD}}=25`$. The dotted line shows a two-parameter fit to $`a\mathrm{exp}(bN_{\mathrm{poly}})`$. We observe a clear exponential decay as expected.

Figure 6 shows the $`N_{\mathrm{poly}}`$ dependence of the averaged noisy Metropolis acceptance rate $`P_{\mathrm{𝑐𝑜𝑟𝑟}}`$. We observe consistent results to the theoretical ansatz $`\mathrm{erfc}(\left(a\mathrm{exp}(bN_{\mathrm{poly}})\right)^{1/2}/2)`$, where the dotted curve shows the ansatz with $`a`$ and $`b`$ obtained in Figure 5.

We show the $`N_{\mathrm{poly}}`$ dependence of the averaged plaquette $`P`$ in Figure 7. The horizontal line shows the constant fit to the results. The fit error is presented by the dashed lines. The results does not depend on $`N_{\mathrm{poly}}`$ as it should be. In Figure 8, we plot the MD step size, $`dt=1/N_{\mathrm{MD}}`$, dependence of the averaged plaquette together with the results with the $`R`$-algorithm. We employ $`N_{\mathrm{poly}}=200`$ for the PHMC algorithm in the figure. Since the $`R`$-algorithm has $`O(dt^2)`$ errors, we fit the results with the $`R`$-algorithm with $`f(dt)=adt^2+bdt^3+c`$ as shown by the dotted curve. The horizontal lines show the constant fit to the PHMC results again. The result with the PHMC algorithm does not depend on $`dt`$ and is consistent with that of the zero step size limit of the $`R`$-algorithm. With these observations we conclude that the PHMC algorithm works perfectly on the small size lattice.

### 7.2 Results on the middle size lattice

We plot the $`N_{\mathrm{poly}}`$ dependence of $`dS`$ for $`am=0.01`$ and $`0.02`$ in Figs. 9 and 10, respectively. Here $`N_{\mathrm{MD}}=50`$ and $`a\mathrm{\Lambda }_{\mathrm{max}}=2.28`$ are used for both masses. A clear exponential decay is observed in both figures. These behaviors are similar to those seen for the small size lattice.

For the $`am=0.01`$ case, which corresponds to the ratio of pseudo-scalar and vector meson masses $`m_{\mathrm{PS}}/m_\mathrm{V}0.61`$ , we need $`N_{\mathrm{poly}}`$ 500-600 for reasonable acceptance rate for the noisy Metropolis test. Moreover, assuming $`a\mathrm{\Lambda }_{\mathrm{max}}/am1`$ and independence of $`a\mathrm{\Lambda }_{\mathrm{max}}`$ from $`am`$, we expect that $`N_{\mathrm{poly}}`$ behaves as $`a\mathrm{\Lambda }_{\mathrm{max}}/am`$ in order to keep the residual at a constant level (see below Eq. (3.6)). This leads us to suspect that for simulations with much lighter quark masses a polynomial of order $`O(1000)`$ is required.

Our PHMC algorithm has several problems with such a large order polynomial. One is the memory cost in the calculation of the MD force from the polynomial pseudo-fermion. Fortunately this would not be a serious hindrance in nowadays high performance computing since memory cost is relatively low compared to the Wilson fermions (the KS fermions do not have spin indices). Another problem is the extraction of the polynomial coefficients of $`Q_{N_{\mathrm{poly}}}`$ from the original polynomial $`P_{N_{\mathrm{poly}}}`$ as described in section 2.

### 7.3 Results on the large size lattice

We show the results on the large size lattice in Table 2. We quote the averaged plaquette value with the $`R`$-algorithm from Ref. . We observe a roughly $`2\sigma `$ deviation, which may be ascribed to a finite step size error inherent in the $`R`$-algorithm. The number of multiplication of hopping matrices $`M_{eo}`$ and $`M_{oe}`$ measured in the program is roughly proportional to $`N_{\mathrm{poly}}`$, which is expected from the discussion in section 6. We will discuss the finite step size error of the $`R`$-algorithm and compare the computational cost of the PHMC algorithm to that of the $`R`$-algorithm after describing the numerical property of the reversibility violation, the Lanczos algorithm, and $`dS`$ for the PHMC algorithm.

As described in section 4, we investigate the reversibility violation on the large size lattice using the following observables:

$`\mathrm{\Delta }H/H`$ $``$ $`|H(t_rt_r)H(0)|/H(0),`$ (7.1)
$`\mathrm{\Delta }U`$ $``$ $`\sqrt{{\displaystyle \underset{n,\mu ,a,b}{}}\left|(U_\mu )_{a,b}(n)(t_rt_r)(U_\mu )_{a,b}(n)(0)\right|},`$ (7.2)
$`\mathrm{\Delta }P`$ $``$ $`\sqrt{{\displaystyle \underset{n,\mu ,a,b}{}}\left|(P_\mu )_{a,b}(n)(t_rt_r)(P_\mu )_{a,b}(n)(0)\right|},`$ (7.3)
where $`X(0)`$ means the observable $`X`$ calculated at the initial configuration at $`t=0`$ in the MD evolution and $`X(t_rt_r)`$ the observable calculated at the reversed configuration which is obtained from the initial configuration at $`t=0`$ by integrating the MD equation to $`t=t_r`$ and then integrating back to $`t=0`$. The trajectory length is $`t_r=1`$. We measured these quantities using 20 configurations which are separated by 5 trajectories. We observe

$`\mathrm{\Delta }H/H`$ $`=`$ $`0.26(6)\times 10^{15},`$
$`\mathrm{\Delta }P/\sqrt{9\times 4\times N_{\mathrm{vol}}}`$ $`=`$ $`0.4162(7)\times 10^{14},`$
$`\mathrm{\Delta }U/\sqrt{9\times 4\times N_{\mathrm{vol}}}`$ $`=`$ $`0.1484(2)\times 10^{14},`$ (7.4)
with the $`N_{\mathrm{poly}}=300`$ PHMC algorithm. The errors are estimated with the binned jackknife method. These values are at a completely satisfactory level with the double precision arithmetic. Consequently it is concluded that the method we employed for the force calculation of the polynomial pseudo-fermion is stable in the MD evolution and does not cause violation of reversibility.

The validity of the CG based stopping criterion for the Lanczos method to calculate $`dS`$ in the noisy Metropolis test is also investigated on the same 20 configurations. We measured $`dS`$, and two residuals defined by

$`r_1`$ $`=`$ $`\left|\left((W[\widehat{D}_{oo}])^{N_f/8}\right)^{8/N_f}\eta _oW[\widehat{D}_{oo}]\eta _o\right|/|W[\widehat{D}_{oo}]\eta _o|,`$ (7.5)
$`r_2`$ $`=`$ $`\left|W[\widehat{D}_{oo}^{}]\left((W[\widehat{D}_{oo}^{}])^{N_f/4}\right)^{4/N_f}\zeta _o\zeta _o\right|/|\zeta _o|,`$ (7.6)
where $`\zeta _o`$ is defined in Eq. (5.3), and the number of iteration of the Lanczos iteration by varying the tolerance level $`\mathrm{𝑡𝑜𝑙}`$ for $`err_1`$.

Figure 11 shows the convergence behavior of the above quantities. The number of iteration increases step by step with the decreasing stopping condition (top of the figure). The residuals Eqs. (7.5) and (7.6) stagnate around a $`O(10^{12})`$ level (middle of the figure). As the residuals are defined as relative error, one may suspect that the number of $`O(10^{12})`$ is not sufficient with the double precision arithmetic. We think this is due to the accumulation of round-off errors in applying the Lanczos iteration several times to calculate the residuals. Namely, for the explicit residual calculation, we need four times application of the Lanczos iteration for $`r_1`$ and twice for $`r_2`$ in the $`N_f=2`$ case, and the Lanczos iteration does not span the completely same Krylov subspace in each application.

For the correctness of the algorithm, $`dS`$ itself is more important. To see the stopping condition dependence of $`dS`$, we measured the convergence of $`|dSdS_{end}|`$ as the metric and show the result in the bottom of Figure 11, where $`dS_{end}`$ is $`dS`$ at the most stringent stopping condition $`\mathrm{𝑡𝑜𝑙}=10^{14}`$. Since the change is too rapid against the change of the number of iteration, we could not see the exponential decay. The metric stagnates around $`O(10^{10})`$. The reason is understood as follows. $`dS`$ is defined as the difference of $`\zeta _o^{}W[\widehat{D}_{oo}^{}]^{N_f/4}\zeta _o`$ and $`|\eta _o|^2`$ in Eq. (5.3). Numerically it is observed that $`\zeta _o^{}W[\widehat{D}_{oo}^{}]^{N_f/4}\zeta _o`$ and $`|\eta _o|^2`$. have almost the same values of about $`3\times 16^4/2O(10^5)`$, and the resulting $`dS`$ is of $`O(10^1)`$. Within double precision arithmetic, the subtraction of $`O(10^5)`$ from $`O(10^5)`$ yielding $`O(10^1)`$ for $`dS`$ means that $`dS`$ only has 9-10 significant figures. The stagnation around $`O(10^{10})`$ would occur in such a case. We stress that the 9-10 significant figures for $`dS`$ is sufficient in the realistic simulations with $`O(10^4)`$ trajectories. No visible effect from the choice of the CG-based stopping criterion is observed.

The plaquette value of the $`R`$-algorithm from Ref.  is larger than that of the PHMC algorithm about $`2\sigma `$ as shown in Table 2. This may contradict the previous observation on the small size lattice where the plaquette value of the $`R`$-algorithm approaches the value of the PHMC algorithm from smaller value. To make clear the $`dt`$ dependence of the plaquette with the $`R`$-algorithm, we imported the program of the $`R`$-algorithm into SR8000-F1 model and produced several trajectories with the $`R`$-algorithm on the large size lattice.

Table 3 shows the results with the $`R`$-algorithm on the large size lattice. The definition of the norm $`\mathrm{𝑟𝑒𝑠}`$ for the stopping criterion of the CG algorithm required in the $`R`$-algorithm is the same as that of Ref.  (where the symbol $`r`$ is used instead of $`\mathrm{𝑟𝑒𝑠}`$). We do not observe clear stopping criterion dependence on $`P`$. Figure 12 shows the $`dt`$ dependence of $`P`$ on the large size lattice. Open circles are the results with the $`R`$-algorithm produced for the comparison. Filled square is the previous result with the $`R`$-algorithm . Open triangles are the results with the PHMC algorithm scatted around $`dt=0.02`$ for clarity of presentation. We observe that the $`dt`$ behavior of the $`R`$-algorithm in small $`dt`$ region is slightly different from that of the small size lattice (Fig. 8), while the behavior at large $`dt`$ region is similar to each other. We also observe that the value of $`dt`$ where $`P`$ largely deviate from the value at the limit $`dt0`$ is different (it is $`dt0.01`$ in Fig. 8 and $`dt0.040.05`$ in Fig. 12).

Although we cannot make clear statement on the discrepancy of $`dt`$ behavior between two lattice sizes, we can say that the $`dt`$ dependence is affected by the physical situation and parameters. The reason is as follows. The error analysis of the $`R`$-algorithm by $`dt`$ perturbation fails at large $`dt`$. More precisely it is said that the point where the perturbative analysis fails is governed by $`dt/m`$ with $`m`$ the lightest fermion mode in the $`R`$-algorithm . We do not tune the input parameters for these two simulation sets. It is natural that the physical lightest fermion mode is different between the small and large size lattices. Thus we consider that this discrepancy comes from the different physical parameter and situation between the small and large size lattices and $`P`$ with the $`R`$-algorithm at $`dt=0.02`$ is accidentally larger than that with the PHMC algorithm in the large size lattice. We stress that the result from the PHMC algorithm is consistent with the limit $`dt0`$ of the results with the $`R`$-algorithm.

The computational cost of the PHMC algorithm of $`N_{\mathrm{poly}}=400`$ is comparable to that of the $`R`$-algorithm with $`\mathrm{𝑟𝑒𝑠}=10^8`$ at $`dt=0.02`$ in terms of #Mult/traj on the large size lattice. We need to take the limit $`dt0`$ and use sufficiently small $`\mathrm{𝑟𝑒𝑠}`$ in the $`R`$-algorithm theoretically, which requires significantly large amount of computational cost for the $`R`$-algorithm. The actual computational time for the PHMC algorithm with $`N_{\mathrm{poly}}=400`$ on the large size lattice was measured as $`136`$ sec for unit trajectory with SR8000-F1 4-nodes (peak speed : 12$`\times `$4 GFlops, sustained speed : 3.5$`\times `$ 4 GFlops) and it is $`120`$ sec for the $`R`$-algorithm with $`\mathrm{𝑟𝑒𝑠}=10^8`$ (sustained speed : 3.4$`\times `$ 4 GFlops) at $`dt=0.02`$. The same computational speed is achieved because both programs make use of the common subroutine for the hopping matrix multiplication. Thus we conclude that the PHMC algorithm works on the lattice size of $`16^4`$ with a quark mass $`am=0.02`$ corresponding to $`m_{\mathrm{PS}}/m_\mathrm{V}0.69`$  within reasonable computational time and is applicable to realistic simulations. The advantage of exact algorithm is very clear in this situation.

## 8 Conclusion

In this paper, we have presented an exact algorithm for dynamical lattice QCD simulation with the KS fermions where single flavor quark is defined as the $`1/4`$ power of the KS fermion matrix. The algorithm is an extension of the polynomial Hybrid Monte Carlo (PHMC) algorithm in which the Hermitian polynomial approximation is applied to the fractional power of the KS fermion matrix. The Kennedy-Kuti noisy Metropolis test is incorporated to make the algorithm exact.

We introduced two types of polynomial approximations and corresponding PHMC algorithms. One algorithm approximates $`\widehat{D}_{oo}^{N_f/8}`$ with a single polynomial $`P_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]`$ (case A), and the other $`\widehat{D}_{oo}^{N_f/4}`$ with a squared polynomial $`|Q_{N_{\mathrm{poly}}}[\widehat{D}_{oo}]|^2`$ (case B). For the noisy Metropolis test, we made use of a Lanczos-based Krylov subspace method to calculate the fractional power of the correction matrix. The efficiency of the two algorithms was estimated, and it was found that the latter (case B) had better performance than that of the former by a factor two.

We tested our algorithm (case B) using the Hermitian Chebyshev polynomial for the case of two-flavor QCD on three lattice sizes. Results on a small lattice of $`8^3\times 4`$ demonstrated that the algorithm works correctly, e.g., the averaged plaquette value agrees with that of the $`R`$-algorithm after extrapolation to the zero step size limit. We have also shown that our algorithm works on a moderately large lattice of size $`16^4`$, albeit for a rather heavy quark mass of $`m_{\mathrm{VS}}/m_\mathrm{V}0.69`$, within reasonable simulation costs compared to that of the $`R`$-algorithm.

There are several points that require improvements with our work. One of the points concerns the fact that the calculation of the polynomial coefficients in case B for splitting the original polynomial becomes progressively difficult toward lighter quark masses. While this is not a limitation of the PHMC algorithm itself, solutions should be found to solve this problem for future realistic simulations since the case A algorithm, which has no such problem, is expected to be twice slower than than the case B algorithm. Another point is that further improvement of the algorithm may be achieved by optimizing the polynomial approximation, and by combining the preconditioning technique and the polynomial approximation.

Anticipating progress on these fronts, we conclude that our algorithm provides an attractive method for dynamical KS fermion simulations for $`2+1`$-flavor QCD without systematic errors originating from the simulation algorithm.

## Acknowledgments

This work is supported by the Supercomputer Project No. 79 (FY2002) of High Energy Accelerator Research Organization (KEK), and also in part by the Grant-in-Aid of the Ministry of Education (Nos. 11640294, 12304011, 12640253, 12740133, 13135204, 13640259, 13640260, 14046202, 14740173). N.Y. is supported by the JSPS Research Fellowship.


text:
# High-t Meson Photo- and Electroproduction: a Window on the Partonic Structure of Hadrons

## 1 Introduction

Exclusive photoproduction of mesons at high momentum transfer $`t`$ offers us a fantastic tool to investigate the partonic structure of hadrons. The high momentum transfer $`t`$ implies that the impact parameter is small enough to allow a short distance interaction between, at least, one constituent of the probe and one constituent of the target. In addition, the exclusive nature of the reaction implies that all the constituents of the probe and the target be in the small interaction volume, in order to be able to recombine into the well defined particles emitted in the final state.

This enables us to identify, to determine the role and to access the interactions between each constituent of hadrons.

The relevant constituents depend on the scale of observation. At low momentum transfer $`t`$, a comprehensive description  of available data is achieved by the exchange of a few Regge trajectories, between the probe and the target. At very high (asymptotic) momentum transfer $`t`$, the interaction should reduce to the exchange of the minimal number of gluons, in order to share the momentum transfer between all the current quarks which recombine into the initial and final hadrons . Here, dimensional counting rules lead to the famous power law behavior of the various cross sections, but it is unlikely that the present generation of high luminosity facilities gives access to this regime.

To date, CEBAF at Jefferson Laboratory is the only facility which allows to reach momentum transfers $`t`$ up to 6 GeV2 in exclusive reactions, at reasonably high energy ($`s`$ up to 12 GeV2). This range will be considerably enlarged when the 12 GeV CEBAF upgrade becomes a reality.

The range of momentum transfers accessible at JLab corresponds to a resolving power of the order of $`0.1`$ to $`0.2`$ fm. It is significantly smaller than the size of a nucleon, but comparable to the correlation lengths of partons (the distance beyond which a quark or a gluon cannot propagate and hadronizes). At this scale, the relevant degrees of freedom are the constituent partons, whose lifetime is short enough to prevent them to interact and form the mesons which may be exchanged in the $`t`$-channel, but is too long to allow to treat them as current quarks or gluons. This is the non perturbative partonic regime, where the amplitude can be computed as a set of few dominant Feynman diagrams which involve dressed quarks and gluons, effective coupling constants and quark distributions . This resembles the treatment of meson exchange mechanisms at low energies .

The first results , recently released by JLab, support such a picture and provide a link with Lattice QCD predictions on the gluon propagator , modelisations of quark wave functions in hadrons  and Regge saturating trajectories .

## 2 Gluon exchange

Measurement of phi photoproduction selects gluon exchanges. Since the $`\varphi `$ meson is predominantly made of a strange $`s\overline{s}`$ quark pair, and to the extent that the strangeness content of the nucleon is small, quark interchange mechanisms are suppressed.

The destructive interference between the two graphs in the bottom of Fig. 1 (where each gluon couples to the same quark in the nucleon, but may couple to a different quark in the vector meson) leads to a node in the cross-section (dashed line) depicted in Fig. 2. This node is filled when the two gluons are allowed to couple also to two different quarks in the proton (faint solid line), giving access to their correlation. At the largest momentum transfer $`t`$, nucleon exchange in the $`u`$-channel takes over, but the corresponding peak move towards higher $`t`$ when the incoming photon energy increases (solid curves, marked 3.5 and 4.5 GeV), leaving more room to access the two gluon exchange contribution.

The combined use of the high luminosity of CEBAF and of the large acceptance of the CLAS set up at JLab made possible the measurement of cross sections, more than two orders of magnitude below previous ones, in a virgin domain up to $`t`$= 5.5 GeV2. The data at $`E_\gamma `$= 3.5 GeV have been published , while the data at $`E_\gamma `$= 4.5 GeV are still preliminary and their errors will soon decrease. This is however good enough to exhibit the trend of the cross section, to confirm the move of the $`u`$-channel peak toward higher momentum transfers and to establish the relevance of the two gluon exchange description.

The key to the success of such a good agreement is the use  of the correlated quark wave function of Ref.  and the Lattice gluon propagator of Ref. . For instance, if a perturbative gluon propagator had been used, one would have obtained a much more steep $`t`$ dependency of the cross-section (see Ref.  for more details). In the gluon loop (Fig. 1), the virtuality of each gluon is on average $`t/4`$, i.e. about 1 GeV2 at $`t=`$ 4 GeV2 where, according to Fig. 2 in Ref. , the lattice gluon propagator exhibits strong non perturbative corrections. On the contrary, it reaches it asymptotic behavior around a virtuality of 4 GeV2: this requires a momentum transfer $`t`$ of about 16 GeV2, which will be achievable when CEBAF is upgraded to 12 GeV.

I refer to the talk of K. McCormick  for the analysis of the tensor polarisation of the $`\varphi `$ which confirms the dominance of the $`u`$-channel contribution at the highest momentum transfer.

## 3 Quark exchange

Quark exchange mechanisms are not supressed in the photoproduction of $`\rho `$ and $`\omega `$ mesons, which are mostly made of light quarks.

The most stricking example is the photoproduction of $`\omega `$ meson (Fig. 3). At low momentum transfer, $`\pi `$ exchange (dashed line) takes over two gluon exchange (dash-dotted line). At the largest momentum transfer, $`u`$-channel proton exchange dominates and accounts for the backward angle cross section. Here the experimental node is well reproduced by the use of a non-degenerate Regge trajectory for the nucleon. At intermediate transfer, the use of a linear Regge trajectory for the $`\pi `$ leads to a vanishing cross section, while the use the saturated Regge trajectory (full line), which already led to a good accounting of the cross section of the $`\gamma p\pi ^+n`$ reaction , enhances the cross section by two orders of magnitude.

This is really a parameter free prediction, which has been beautifully confirmed by the recent CLAS data: I refer to the talk of M. Battaglieri , for a more detailed discussion in the $`\omega `$ sector as well as the $`\rho `$ sector.

## 4 Compton scattering

A few GeV real photon has a significant hadronic component. Due to the uncertainty principle, it fluctuates into vector mesons (or quark anti-quark pairs of various flavors) over a distance which exceeds the nucleon size. For instance, a 4 GeV real photon fluctuates into a $`\rho `$ meson over about 2.7 fm. The real Compton scattering cross section is therefore related to the $`\rho `$ meson photoproduction cross section by a simple multiplicative factor: $`4\pi \alpha _{em}/f_V^2`$, where $`f_V`$ is the radiative decay constant of the vector meson. As shown in Fig. 4, the comparison with this model  and the old Cornell data confirm this conjecture.

More interesting, the model predicts also spin observables. As an example, Fig. 5 compares the predicted longitudinal spin (helicity) tranfer (solid line in the upper part) to the hard scattering models  (solid line in the bottom part) and the soft ”handbag” model  (dashed line). The strong variation at backward angles comes from the $`u`$-channel baryon exchange. The preliminary data from Hall A at JLab confirms this prediction. I refer to the talk of A. Nathan , for a more detailed account and discussion of future experimental prospects.

## 5 Conclusion

A consistent picture is emerging from the study of exclusive photoproduction of vector meson and Compton scattering. At low momentum transfer, it relies on diffractive scattering of the hadronic contents of the photon, in a wide energy range from threshold up to the HERA energy domain. At higher momentum transfer, it relies on a partonic description of hard scattering mechanisms which provides us with a bridge with Lattice Gauge calculations. The dressed gluon and quark propagators have already been estimated on lattice. One may expect that correlated constituent quark wave functions (at least their first moments) will soon be available from lattice. An estimate, on lattice, of the saturating part of the Regge trajectories is definitely called for.

To day, JLab is the only laboratory which allows to explore this regime, thanks to its high luminosity. Its current operation, at 4–6 GeV, has already revealed a few jewels. Its energy upgrade to 12 GeV will expand to higher values the accessible range in momentum transfer and allow to extend these studies to virtual photon induced reactions. It will permit, among others, a more comprehensive study of the onset of asymptotic hard scattering regime and of the role of correlations between quarks.

## Acknowledgments

This work has been partly supported by European Commission under Contract HPRN-CT-2000-00130. Part of the results presented in this note arose from collaboration and discussions with F. Cano, over the past two years.