text:
# Glass transition in a simple stochasic model with back-reaction

## 1 Introduction

Glass transition and slow relaxation in systems characterised by weak ergodicity breaking remains still an open area, despite many efforts and numerous significant results in the last decade angell\_95 ; pa\_ste\_ab\_an\_84 ; pa\_me\_dedo\_85 ; ki\_thi\_87 ; ki\_thi\_87a ; ki\_thi\_wo\_89 ; ritort\_95 ; fra\_he\_95 ; cu\_ku\_93 ; fra\_me\_pa\_pe\_98 ; ba\_bu\_me\_96a ; kurchan\_92 ; cu\_ku\_95 ; ma\_pa\_ri\_94b ; nieuwenhuizen\_98 ; ca\_gia\_pa\_99 ; ca\_gia\_pa\_99a ; me\_pa\_99 ; co\_pa\_ve\_00 ; cu\_ku\_pa\_ri\_95 ; par\_sla\_00 ; kaw\_kim\_01 . Among the host of diverse phenomena we are motivated here mostly by the effect of dynamical arrest in colloidal matter fie\_sol\_cat\_99 ; fu\_cat\_03 ; cates\_02 ; pue\_fu\_cat\_02 ; fu\_cat\_02 ; zac\_fof\_daw\_bul\_sci\_tar\_02 ; law\_rea\_mcc\_deg\_tar\_daw\_02 ; sti\_vla\_lop\_roo\_mei\_02 ; fab\_got\_sci\_tar\_thi\_99 , observed experimentally and thoroughly investigated by numerical simulations and Mode-Coupling method. Below the transition point, the dynamics effectively leads to a glassy state with diverging viscosity, however the static thermodynamic transition may not be identifiable. Indeed, dynamical or structural arrest demonstrates the glass transition as a purely dynamic and self-consistent phenomenon, where casual slowdown of certain particles prevents some other particles from moving, which may slow down the others even more etc. The self-consistent nature of the phenomenon is reflected by the analytical approaches available now. One of the most striking phenomena in colloids, suspensions and granular matter is the non-Newtonian response to mechanical perturbation. On one hand, we can have shear-thinning, which amounts to a decrease if viscosity due to applied field, which can be interpreted as restoration of ergodicity due to perturbation cates\_02 . On the other hand, increase of viscosity may result in shear thickening or even jamming, typically observed in particulate or granular matter cat\_wit\_bou\_cla\_99 ; ber\_bib\_sch\_02 .

From the theoretical side, the Mode-Coupling (MC) equations provide us with a well-established framework, capable of explaining a good deal of experimental data gotze\_91 ; go\_sj\_92 ; bouchaud ; sjogren\_03 . The attempts to derive the MC equations starting from the Hamiltonian of the system were successful in the mean-field approximation. It was perhaps the $`p`$-spin spherical model ki\_thi\_87 ; ki\_thi\_87a ; ki\_thi\_wo\_89 ; cu\_ku\_95 , where the machinery reached the farther edge of our current understanding of the phenomenon.

However, the bottom-up approach starting with writing explicit Hamiltonians is far from being complete. The presence of the reparametrisation invariance cu\_ku\_95 ; cu\_ku\_98  leaves the numerical solution of the MC equations as the only means for obtaining the true time-dependence of the correlation and response functions. Also the mean-field approximation generally used now seems to be very difficult to overcome.

The serious difficulties remaining in using the more advanced MC techniques leave the space for more simple phenomenological approaches. We want to follow this path in the present work.

Indeed, the mathematical substance of the Mode Coupling method can be summarised by saying that the time dependence of the correlation (and response) functions depends non-linearly and in time-delayed manner on these functions themselves. Actually, the memory kernel in the MC equations, which is primarily dictated by the properties of the reservoir, depends of the system dynamics.

We may represent the dynamics of the system by a stochastic process and the parameters of the process depend on time through the averaged properties of the process itself. In order to study generic properties of such problems it can be useful to establish a simple idealised model, which would capture the essential mathematic ingredients while avoiding the complications which arise from choosing a specific Hamiltonian at the beginning. The most important ingredient in such an idealised model should be the mechanism of the back-reaction.

We introduced recently chv\_sla\_02  a very simple stochastic process, in which the back-reaction leads to rich dynamic behaviour. The main characteristics was the presence of a phase transition from ergodic to non-ergodic phase. The principal aim of the present work is to investigate analytically some of the properties of the transition and from the numerical solution of the corresponding differential equations infer the non-trivial critical behaviour.

## 2 Langevin equation with back–reaction

### 2.1 System of coupled particles

The model system we will have in mind will be composed of particles, relaxing to their equilibrium positions under the influence of surrounding particles. They can be viewed as colloidal particles immersed in a solvent, but the formulation of our model is generic enough to allow for other interpretations as well, e. g. they can be viewed as micro-domains in a relaxor ferroelectric material.

The time evolution of the model can lead to dynamical arrest, where particles are locked in their positions by surrounding particles, which are also locked in their turn. Therefore, the dynamics can lead to the spatially disordered but time-stable stationary state with glass properties. The indication of the glass transition will be the non-zero value of the Edwards-Anderson order parameter and sensitivity to initial conditions. The interaction between particles will be taken into account on a phenomenological level; if we concentrate on a randomly selected particle (single relaxor), the external field from the rest of the system (reservoir) will change as the states of the other particles (relaxors in reservoir) evolve. The changes in the local external field will be the more rapid the faster is the evolution of the other particles. This leads to the idea of expressing the intensity of the changes in th external field through the velocity of movement of the relaxors in the reservoir. As we suppose all particles to be statistically identical, the movement of our single relaxor should be in probabilistic sense equivalent to the movement of any relaxor within the reservoir. This consideration closes the loop.

We will try to express the intensity of the changes of the local field through the averaged properties of the movement of the single relaxor itself. This introduces the idea of back-reaction: the probabilistic properties of the reservoir dictate the system evolution and the averaged system dynamics tunes the properties of the reservoir itself.

To be more specific, our single relaxor will be described by the continuous real stochastic variable $`𝖷(t)`$. It will evolve under influence of the environmental force, represented by the stochastic variable $`𝖰(t)`$. The force will be modelled by a two–valued random process $`𝖰(t)\{q,+q\}`$ horsthemke , jumping at random instants. Occurrence of the jumps are governed by a self-exciting point process cox ; snyder  with time-dependent intensity $`\frac{1}{2}\lambda (t)`$.

For given (friction-reduced) force $`𝖰(t)`$ the single relaxor is described by the Langevin equation pomeau ; sibani ; morita

$$\frac{d}{dt}𝖷(t)=\gamma 𝖷(t)+𝖰(t)$$
(1)
with initial conditions $`𝖷(0)=𝖷_0`$ and $`𝖰(0)=𝖰_0`$.

We may consider the process $`𝖷(t)`$ as a movement of an over-damped particle which slides in a parabolic potential well; the parabola jumps between two positions at random instants, with time-dependent rate $`\frac{1}{2}\lambda (t)`$.

The intensity of the process, or the frequency of jumps, $`\frac{1}{2}\lambda (t)`$, is related to the movement of the surrounding relaxors, considered as a reservoir. Intuitively the frequency must be smaller if the movement of the relaxors is slower. Therefore, the function $`\lambda (t)`$ should be coupled to the velocity $`𝖵(t)=\frac{d}{dt}𝖷(t)`$.

However, it is not obvious a priori what should be the specific functional dependence. We only require that the dependence is described by a non-negative function analytic at the origin. The simplest choice satisfying this property is

$$\lambda (t)=ϵ\frac{\gamma }{q^2}𝖵^2(t)$$
(2)
where $`ϵ`$ is the dimensionless coupling strength parameter. The latter prescription is the form of the back-reaction we will study in the following. The model described above is sketched schematically in the Fig. 1.

The parameters $`\gamma `$ and $`q`$ can be in principle rescaled to 1 by appropriate choice of the units of time and length. Therefore, the coupling strength $`ϵ`$ remains to be the only physically relevant parameter tuning the behaviour of the system. As we will see, there is a qualitative change in the behaviour of the system at a certain critical value of $`ϵ`$.

### 2.2 Properties of the environmental force $`𝖰(t)`$

The force $`𝖰(t)`$ is a time-inhomogeneous Markov process, therefore its properties are fully described by a master equation. More specifically, let us define the probabilities

$$\pi _\pm (t)=\mathrm{Prob}\left\{𝖰(t)=\pm q\right\}$$
(3)
making the vector $`𝝅(t)=\left(\begin{array}{c}\pi _{}(t)\\ \pi _+(t)\end{array}\right)`$ which satisfies the Pauli master equation

$$\frac{d}{dt}𝝅(t)=\frac{1}{2}\lambda (t)\left(\begin{array}{cc}\hfill 1& \hfill 1\\ \hfill 1& \hfill 1\end{array}\right)𝝅(t).$$
(4)
Solving the equation amounts to calculation of the corresponding time-ordered exponential. The averages and correlation functions can be expressed through the integrated intensity

$$\mathrm{\Lambda }(t)=_0^t\lambda (t^{})dt^{}.$$
(5)
The time dependence of the vector $`𝝅(t)`$ can be written through the semi-group operator

$$R(t,t_0)=\left(\begin{array}{cc}\hfill 1& \hfill 0\\ \hfill 0& \hfill 1\end{array}\right)+\frac{1}{2}\left(\mathrm{e}^{\mathrm{\Lambda }(t)+\mathrm{\Lambda }(t_0)}1\right)\left(\begin{array}{cc}\hfill 1& \hfill 1\\ \hfill 1& \hfill 1\end{array}\right)$$
(6)
as $`𝝅(t)=R(t,t_0)𝝅(t_0)`$. Note that the semi-group operator obeys $`R(t,\tau )R(\tau ,t_0)=R(t,t_0)\tau [t_0,t]`$ which testifies the Markov property of the process.

More explicitly, we find

$$𝖰(t)=𝖰_0\mathrm{exp}(\mathrm{\Lambda }(t))$$
(7)
$$𝖰(t)𝖰(t_1)=q^2\mathrm{exp}(|\mathrm{\Lambda }(t)\mathrm{\Lambda }(t_1)|)$$
(8)
Similarly, also the higher correlation functions can be written as products of exponentials with combinations of $`\mathrm{\Lambda }(t)`$ with appropriate time arguments in the exponents. In fact, higher order correlation functions factorise into product of first and second order correlations, e. g.

$$\underset{l=1}{\overset{2k}{}}𝖰(t_l)=\underset{l=1}{\overset{k}{}}𝖰(t_{2l})𝖰(t_{2l1})$$
(9)
for $`t_{2k}t_{2k1}\mathrm{}t_2t_1`$. Another consequence is, that the cumulants of higher order than two vanish.

Note that the function $`\mathrm{\Lambda }(t)`$ is non-decreasing and can either diverge (if $`lim_t\mathrm{}\lambda (t)>0`$) or assume a finite limit for $`t\mathrm{}`$, if $`\lambda (t)`$ approaches 0 fast enough.

## 3 Glass transition and asymptotic relaxation

### 3.1 Equations for moments

For any given realisation of the process $`𝖰(t)`$, the formal solution of Eq. (1) is

$$𝖷(t)=𝖷_0\mathrm{e}^{\gamma t}+_0^t\mathrm{e}^{\gamma (tt^{})}𝖰(t^{})dt^{}.$$
(10)
If the function $`\lambda (t)`$ were known, various moments (and correlation functions) of the random process $`𝖷(t)`$ could have been computed from (10) using the expressions (7) and (8). However, in our case the function $`\lambda (t)`$ should be computed from the condition (2), relating it to the second moment of the time derivative of $`𝖷(t)`$. This suggests that sufficiently broad set of moments of $`𝖷(t)`$ and $`𝖵(t)`$ may provide a closed set of ordinary differential equations. The solution of this set will yield the closed description of the behaviour of our model.

Indeed, we can define four auxiliary functions

$`s_1(t)`$ $`=`$ $`\mathrm{e}^{\mathrm{\Lambda }(t)}`$ (11)
$`s_2(t)`$ $`=`$ $`\mathrm{e}^{\gamma t}{\displaystyle _0^t}dt^{}\mathrm{e}^{\gamma t^{}\mathrm{\Lambda }(t^{})}`$ (12)
$`s_3(t)`$ $`=`$ $`\mathrm{e}^{\gamma t\mathrm{\Lambda }(t)}{\displaystyle _0^t}dt^{}\mathrm{e}^{\gamma t^{}+\mathrm{\Lambda }(t^{})}`$ (13)
$`s_4(t)`$ $`=`$ $`\mathrm{e}^{2\gamma t}{\displaystyle _0^t}dt^{}\mathrm{e}^{\gamma t^{}\mathrm{\Lambda }(t^{})}{\displaystyle _0^t^{}}dt^{\prime \prime }\mathrm{e}^{\gamma t^{\prime \prime }+\mathrm{\Lambda }(t^{\prime \prime })}`$ (14)
and express the requested quantities through these functions. For example the average coordinate can be written as

$$𝖷(t)=𝖷_0\mathrm{e}^{\gamma t}+𝖰_0s_2(t).$$
(15)
Similarly, the second moment of the coordinate is

$$𝖷^2(t)=𝖷_0^2\mathrm{e}^{2\gamma t}+2𝖷_0𝖰_0\mathrm{e}^{\gamma t}s_2(t)+2q^2s_4(t).$$
(16)
The functions $`s_1(t)`$ to $`s_4(t)`$ can be found by solving the set of non-linear differential equations

$`\dot{s}_1(t)`$ $`=`$ $`\lambda (t)s_1(t)`$ (17)
$`\dot{s}_2(t)`$ $`=`$ $`\gamma s_2(t)+s_1(t)`$ (18)
$`\dot{s}_3(t)`$ $`=`$ $`1\left(\gamma +\lambda (t)\right)s_3(t)`$ (19)
$`\dot{s}_4(t)`$ $`=`$ $`2\gamma s_4(t)+s_3(t)`$ (20)
with initial conditions $`s_1(0)=1`$, $`s_2(0)=s_3(0)=s_4(0)=0`$. The function $`\lambda (t)`$ occurring in the latter equations is itself a combination of the functions $`s_1(t)`$ to $`s_4(t)`$

$$\begin{array}{cc}\hfill \lambda (t)& =ϵ\left[\gamma 2\gamma ^2s_3(t)+2\gamma ^3s_4(t)\right]+\hfill \\ & +\frac{ϵ\gamma ^2}{q^2}[\gamma 𝖷_0^2\mathrm{e}^{2\gamma t}+2\gamma 𝖷_0𝖰_0\mathrm{e}^{\gamma t}s_2(t)\hfill \\ & 2𝖷_0𝖰_0\mathrm{e}^{\gamma t}s_1(t)].\hfill \end{array}$$
(21)
In the following we will assume that the initial condition of the stochastic process is $`𝖷_0=0`$ and $`𝖰_0=0`$, except explicitly mentioned cases. This lead to significant simplification of the mathematical structure of the equations. Indeed, the equations for $`s_3(t)`$ and $`s_4(t)`$ form a closed pair of equations

$$\begin{array}{cc}\hfill \dot{s}_3(t)=& 1(1+ϵ)\gamma s_3(t)+\hfill \\ & +2ϵ\gamma ^2(s_3(t)\gamma s_4(t))s_3(t)\hfill \\ \hfill \dot{s}_4(t)=& s_3(t)2\gamma s_4(t).\hfill \end{array}$$
(22)
Unfortunately, the system (22) cannot be solved analytically. The best one can do is to transform the non-linear Ricatti-type set (22) to one differential equation of Abel type, whose solution, however, is not generally known. Therefore, we will solve the equations (22) numerically. Nevertheless, there is still a significant amount of information which can be extracted analytically.

The typical results of numerical solution are shown in figures 2 and 3. In Fig. 2 we can see the time evolution of the auxiliary functions $`s_1(t)`$ to $`s_4(t)`$. Fig. 3 shows the evolution of the switching rate $`\lambda (t)`$ and average coordinate $`𝖷(t)`$ for non-zero value of the initial condition $`𝖰_0`$. We can observe qualitatively different behaviour for $`ϵ<1`$ and $`ϵ>1`$: first, the switching rate approaches a non-zero limit for $`ϵ>1`$, while for $`ϵ<1`$ it decays to zero. This means that in the latter case the system effectively freezes. This is further confirmed by the observation that for $`ϵ<1`$ the limit value of the average coordinate depends on the initial conditions, while in the opposite case the dependence on initial conditions is lost for large times, the system equilibrates and the average coordinate converges always to zero. The following sections are mainly devoted to the analytical investigation of the above observations.

### 3.2 Fixed points

The first step in investigating the behaviour of the system (22) is the search for the fixed points $`[s_3^{},s_4^{}]`$ of the dynamics. We found that there are only two fixed points, namely

$$[s_3^{},s_4^{}]=[\frac{1}{\gamma ϵ},\frac{1}{2\gamma ^2ϵ}]$$
(23)
and

$$[s_3^{},s_4^{}]=[\frac{1}{\gamma },\frac{1}{2\gamma ^2}].$$
(24)
Let us denote $`\lambda _{\mathrm{}}`$ the value of $`\lambda (t)`$ calculated at the corresponding fixed point. Using (21) the fixed point (23) yields $`\lambda _{\mathrm{}}=(ϵ1)\gamma `$ and the fixed point (24) yields $`\lambda _{\mathrm{}}=0`$.

The linear stability analysis reveals that for $`ϵ>1`$ the fixed point (23) is stable, while (24) is unstable. On the other hand, for $`ϵ<1`$ the fixed point (24) is stable, while (23) is unstable. The case $`ϵ=1`$ is a marginal one, where both fixed points have one of the eigenvalues equal to 0. Therefore, the value $`ϵ=1`$ marks a transition, whose nature will be further pursued in the following.

### 3.3 Ergodic regime $`ϵ>1`$

In this case the relevant fixed point is (23) and inserting its value to the expressions for the moments of $`𝖷(t)`$ we find that both the average coordinate and the average velocity relaxes to zero. On the other hand, the fluctuations of the coordinate reach positive value, so

$$\begin{array}{cc}\hfill \underset{t\mathrm{}}{lim}𝖷(t)& =0\hfill \\ \hfill \underset{t\mathrm{}}{lim}𝖷^2(t)& =\frac{q^2}{\gamma ^2ϵ}.\hfill \end{array}$$
(25)
The external force switching rate converges to positive constant $`lim_t\mathrm{}\lambda (t)=\gamma (ϵ1)`$. In all cases the quantities of interest converge exponentially to their limit values The rate of convergence is determined by the lowest in absolute value eigenvalue, which is

$$\mu _1=\frac{\gamma }{2}\left(ϵ\sqrt{ϵ^28ϵ+8}\right).$$
(26)
In the interval $`ϵ(42\sqrt{2},4+2\sqrt{2})`$ the eigenvalue acquires a non-zero imaginary part, which means that oscillatory behaviour is superimposed over the exponential relaxation.

The overall picture is the following. The back-reaction leads to self-adjustment of the switching rate of the external force exerted by the reservoir. The coordinate fluctuates around the origin and these fluctuations are stationary. Therefore, the stationary regime of the system corresponds to the primitive version with fixed $`\lambda `$, except the fact that the value of $`\lambda `$ is not given from outside, but tuned by the dynamics itself. We call this regime ergodic, because the particles do not freeze at some value of the coordinate $`𝖷(𝗍)`$ but fluctuate forever.

The probability density for the coordinate

$$P(x,t)=\frac{\mathrm{d}}{\mathrm{d}x}\mathrm{Prob}\{𝖷(t)x\}$$
(27)
approaches for $`t\mathrm{}`$ the function chv\_sla\_02 ; morita

$$\underset{t\mathrm{}}{lim}P(x,t)=\frac{\gamma }{q}\frac{(1\stackrel{~}{x}^2)^{(ϵ3)/2}}{B(\frac{ϵ1}{2},\frac{1}{2})}\mathrm{\Theta }(1\stackrel{~}{x}^2),$$
(28)
where $`\stackrel{~}{x}=x\gamma /q`$, $`\mathrm{\Theta }(a)`$ is the Heaviside unit-step function, and $`B(a,b)`$ denotes the Beta-function abramowitz . We can observe a qualitative change at the value $`ϵ=3`$. For $`ϵ>3`$ the limiting distribution (28) has a maximum for $`x=0`$ and approaches $`0`$ at the edges of the support $`[q/\gamma ,q/\gamma ]`$, while for $`ϵ<3`$ it has a minimum at $`x=0`$ and diverges at the edges of the support. The tendency for accumulating the probability close to the points $`\pm q/\gamma `$ when $`ϵ`$ decreases can be regarded as a precursory phenomenon of the transition to the non-ergodic regime, investigated in the next sub-section.

### 3.4 Non-ergodic regime $`ϵ<1`$

In this case we have (24) as stable fixed point. For $`𝖰_0=0`$ the average coordinate converges to 0 again, but the second moment approaches the $`ϵ`$-independent maximum value

$$\underset{t\mathrm{}}{lim}𝖷^2(t)=\frac{q^2}{\gamma ^2}.$$
(29)
As the probability density for the coordinate $`P(x,t)`$ has support limited to the interval $`[q/\gamma ,q/\gamma ]`$, it follows from (29) the the limiting probability density is composed of two $`\delta `$-functions of equal weight $`\frac{1}{2}`$ located at the edges of the latter interval. More generally, for non-symmetric initial condition for the noise, $`𝖰_00`$, the limiting probability density is the sum of two $`\delta `$-functions,

$$\underset{t\mathrm{}}{lim}P(x,t)=\rho _+\delta \left(x\frac{q}{\gamma }\right)+\rho _{}\delta \left(x+\frac{q}{\gamma }\right)$$
(30)
the weights of which depend non-trivially on $`ϵ`$ and the initial condition

$$\rho _\pm =\frac{1}{2}\left(1\pm \frac{𝖰_0}{q}\sigma (ϵ)\right)$$
(31)
where $`\sigma (ϵ)=lim_t\mathrm{}s_1(t)`$.

The function $`\lambda (t)`$ relaxes to zero. Therefore, in this regime, the switching of the external force asymptotically stops and the coordinate $`𝖷(t)`$ approaches either the value $`+q/\gamma `$ or $`q/\gamma `$, where it freezes. So, the coordinate acquires a random but time-independent asymptotic value. More precisely, the mean coordinate $`𝖷(t)`$ approaches a generally non-zero asymptotic value, which depends on the initial condition. This is the manifestation of glassy state in the regime $`ϵ<0`$, characterised by broken ergodicity and non-zero Edwards-Anderson order parameter. This point will be discussed more in detail later in the presentation of correlation functions.

As in the ergodic phase, all quantities relax toward their limit values exponentially for large times. The rate of convergence is governed by the eigenvalue with smallest modulus, which is now

$$\mu _1=2\gamma (1ϵ).$$
(32)
Note that, contrary to the ergodic regime, the eigenvalue is always a real number, so no oscillations occur, at least in the linearised approximation.

### 3.5 Marginal case $`ϵ=1`$

Let us proceed by approaching the marginal case $`ϵ=1`$ from the non-ergodic side, i. e. from below. It might be instructive to cast the equations (22) in terms of the eigenmodes of the linearised approximation. Namely, we can introduce the functions

$$\begin{array}{cc}\hfill \xi (t)& =\frac{1}{(2ϵ1)\gamma }\left(s_3(t)\gamma s_4(t)\frac{1}{2\gamma }\right)\hfill \\ \hfill \eta (t)& =\frac{1}{(2ϵ1)\gamma }\left(s_3(t)+2ϵ\gamma s_4(t)+\frac{1ϵ}{\gamma }\right)\hfill \end{array}$$
(33)
and express the equations (22) in the form

$$\begin{array}{cc}\hfill \dot{\xi }& =2\gamma (1ϵ)\xi +2ϵ\gamma ^3(2ϵ\xi +\eta )\xi \hfill \\ \hfill \dot{\eta }& =\gamma \eta 2ϵ\gamma ^3(2ϵ\xi +\eta )\xi \hfill \end{array}$$
(34)
The function $`\xi (t)`$ has a straightforward physical interpretation: it describes the time evolution of the switching rate of the external force. Indeed, inserting (33) into (21) we get

$$\lambda (t)=2\gamma ^3(12ϵ)\xi (t).$$
(35)
The equations (34) are a convenient starting point for the investigation of the marginal regime. Taking $`ϵ=1`$, the linear term in the equation for $`\dot{\xi }`$ vanishes, while the linear term in the equation for $`\dot{\eta }`$ remains. This suggests that in the long-time regime the value of $`\eta `$ will be negligible compared to $`\xi `$. This consideration will yield the leading term in the relaxation.

Thus, supposing $`|\eta ||\xi |`$ we get the following approximate equation for $`\xi `$

$$\dot{\xi }=4\gamma ^3\xi ^2$$
(36)
which leads to the following asymptotic behaviour

$$\xi (t)\frac{1}{4\gamma ^3}\frac{1}{t},t\mathrm{}.$$
(37)
Now we must check the assumption that $`\eta `$ is negligible compared to $`\xi `$. However, from (34) we can see that the leading term in the relaxation of $`\eta `$ is

$$\eta (t)\frac{1}{\gamma }\frac{1}{t^2},t\mathrm{}$$
(38)
and the assumption is therefore consistent.

The consequence to draw is that in the marginal regime the relaxation becomes power-law with exponent $`1`$. Especially, the relaxation of the switching rate follows the behaviour

$$\lambda (t)\frac{1}{2t},t\mathrm{}.$$
(39)
In Fig. 4 we can compare the numerical solution with the asymptotic behaviour (39). We can see not only that that the function $`\lambda (t)`$ approaches zero according to the power decay (39), but also the corrections to the asymptotic behaviour can be well approximated by a power. Indeed, from the inset in Fig. 4 we can see that

$$\frac{1}{2t\lambda (t)}13t^{0.9},t\mathrm{}.$$
(40)
It is interesting to note that the power in the correction is not an integer, so the naive expansion of the solution in powers of $`t^1`$ cannot be used here. Instead, the behaviour (40) suggests the expression in the form of a continued fraction

$$\lambda (t)=\frac{1}{a_1t^{\alpha _1}+{\displaystyle \frac{1}{a_2t^{\alpha _2}+{\displaystyle \frac{1}{a_3t^{\alpha _3}+\mathrm{}}}}}}$$
(41)
where the values $`a_1=2`$ and $`\alpha _1=1`$ are known exactly and the next pair of parameters is estimated from the numerical solution as $`a_21/6`$ and $`\alpha _20.9`$.

## 4 Correlation functions

Additional information on the properties of the transition from ergodic to non-ergodic behaviour which occurs at the value $`ϵ=1`$ can be gained from the two-time correlation functions. Let us have $`t>t_1>0`$ and define the correlation function

$$C(t,t_1)=𝖷(t)𝖷(t_1).$$
(42)
It can be expressed through the functions $`s_1(t)`$ to $`s_4(t)`$. The most general formula is

$$\begin{array}{cc}\hfill C(t,t_1)& =𝖷_0\mathrm{e}^{\gamma (t+t_1)}+\hfill \\ & +𝖷_0𝖰_0\left[\mathrm{e}^{\gamma t_1}s_2(t)+\mathrm{e}^{\gamma t}s_2(t_1)\right]+\hfill \\ & +q^2[2\mathrm{e}^{\gamma (tt_1)}s_4(t_1)+\hfill \\ & +(s_2(t)\mathrm{e}^{\gamma (tt_1)}s_2(t_1))\frac{s_3(t_1)}{s_1(t_1)}]\hfill \end{array}$$
(43)
although we suppose throughout this section that $`𝖷_0=0`$.

We show in figures 5 (ergodic regime) and 6 (non-ergodic regime) the evolution of correlation functions using the numerical solution for $`s_1`$ to $`s_4`$. We can observe the damping of the correlations in the ergodic regime, while in the non-ergodic regime the correlations converge to a finite limit. Let us now turn to the analytic investigation of the long-time behaviour of the correlation function.

For long enough times we can suppose that we are in the regime of exponential asymptotic relaxation of the functions $`s_1(t)`$ to $`s_4(t)`$ and $`\lambda (t)`$, which is governed by the eigenvalue closest to 0, as given by (26) and (32).

Let us start with the ergodic regime $`ϵ>1`$. We find that for both $`t\mathrm{}`$ and $`t_1\mathrm{}`$ the correlation function behaves like

$$\begin{array}{cc}\hfill C(t,t_1)\frac{q^2}{\gamma ^2ϵ(2ϵ)}& [\mathrm{e}^{\gamma (ϵ1)(tt_1)}\hfill \\ & (ϵ1)\mathrm{e}^{\gamma (tt_1)}].\hfill \end{array}$$
(44)
Important feature of this result is that the correlation in asymptotic regime depends only on the time difference $`tt_1`$ and decays to zero when this difference increases. This supports the picture of the $`ϵ>1`$ phase as a usual ergodic regime without long-time correlations.

The situation is dramatically different for $`ϵ<1`$. Technically speaking, it is important that the function $`\mathrm{\Lambda }(t)`$ has a finite limit for large times, $`lim_t\mathrm{}\mathrm{\Lambda }(t)=\mathrm{\Lambda }_{\mathrm{}}<\mathrm{}`$. For $`ϵ`$ close to $`1`$ (i. e. $`1ϵ1`$) the approach to this limit value can be written in the form

$$\mathrm{\Lambda }(t)\mathrm{\Lambda }_{\mathrm{}}\theta \mathrm{e}^{2\gamma (1ϵ)t}$$
(45)
where $`\theta `$ is a constant depending on $`ϵ`$.

For large times the correction will be small and we can formally write the correlation function as expansion in powers of $`\theta `$. However, we should bear in mind that it is not $`\theta `$ itself, which is small, but the factor $`\mathrm{e}^{2\gamma (1ϵ)t}`$ which appears always together with $`\theta `$.

Finally,

$$\begin{array}{cc}\hfill C(t,t_1)& \frac{q^2}{\gamma ^2}+\frac{\theta q^2}{(\gamma \mu )\gamma }\left(\mathrm{e}^{\mu t}\mathrm{e}^{\mu t_1}\right)\hfill \\ & \frac{2\theta q^2\mu }{(\gamma \mu )(2\gamma \mu )\gamma }\mathrm{e}^{\gamma (tt_1)\mu t_1}+O(\theta ^2)\hfill \end{array}$$
(46)
using $`\mu =2\gamma (1ϵ)`$ for shorter notation.

We can see that the qualitative difference from the ergodic regime consists in the fact that for $`ϵ<1`$ the correlation function converges to a positive $`ϵ`$-independent constant $`q^2/\gamma ^2`$. If we define the Edwards-Anderson order parameter

$$q_{\mathrm{EA}}=\underset{\tau \mathrm{}}{lim}\underset{t_1\mathrm{}}{lim}C(t_1+\tau ,t_1)$$
(47)
we can see that $`q_{\mathrm{EA}}`$ jumps discontinuously from the value $`q_{\mathrm{EA}}=0`$ for $`ϵ>1`$ to $`q_{\mathrm{EA}}=q^2/\gamma ^2`$ for $`ϵ<1`$. This observation represents another evidence that there is a transition from ergodic regime to non-ergodic glassy regime at $`ϵ=1`$. As the Edwards-Anderson parameter is discontinuous at the critical point, the transition should be classified as first-order from this point of view. However, because we do not deal with an equilibrium transition and the phenomenon is of purely dynamical origin, the canonic classification of phase transition as first or second order is of limited relevance here.

## 5 Critical behaviour at $`ϵ1^{}`$

We have already seen that it is possible to characterise the glass transition at $`ϵ=1`$ through the Edwards-Anderson order parameter $`q_{\mathrm{EA}}`$. It has discontinuity at the transition, so the corresponding critical exponent is 0. Here we investigate another quantity, which can play the role of an order parameter, being zero in the ergodic and non-zero in the non-ergodic phase.

The quantity in question will describe how the initial conditions affect the asymptotic vale of the average coordinate. We have already touched this point in Sec. 3.4. Up to now we assumed that the initial condition for the noise is such that $`𝖰_0=0`$. This assumption has the consequence that both in ergodic and non-ergodic regime the average coordinate converges to $`0`$. In this section we investigate the case $`𝖰_00`$.

From (15) and (18)we can see that

$$\underset{t\mathrm{}}{lim}𝖷(t)=\frac{𝖰_0}{\gamma }\sigma (ϵ)$$
(48)
where we defined, as in Sec. 3.4,

$$\sigma (ϵ)=\underset{t\mathrm{}}{lim}s_1(t)$$
(49)
stressing explicitly the dependence on $`ϵ`$. The equations (22) hold also in case $`𝖰_00`$. Therefore, we can proceed without further complications. We only need to plug the solution obtained the same way as in the sections 3 and 4 into the definition (11) of the function $`s_1(t)`$ and find the limit (49).

Let us first present the results of numerical integration. We will turn to the analytic estimate afterwards. The Fig. 7 shows the results of numerical integration, indicating that asymptotically for $`ϵ1`$ the behaviour follows the power law

$$\sigma (ϵ)(1ϵ)^{\frac{1}{2}}.$$
(50)
The exponent $`1/2`$ is still observed only empirically and we do not possess any proof that this is the exact value.

However, we may get some analytical argument in favour of this type of behaviour from the equations (19), (20), and (21), which can be rewritten in slightly different form. Defining new function $`\psi (t)=s_3(t)1/\gamma `$ we can write the set of equations for the pair $`\psi (t)`$ and $`\lambda (t)`$

$$\begin{array}{cc}\hfill \dot{\psi }& =\gamma \psi \frac{1}{\gamma }\lambda \lambda \psi \hfill \\ \hfill \dot{\lambda }& =2\gamma (1ϵ)\lambda +2ϵ\gamma ^2\lambda \psi \hfill \end{array}$$
(51)
with initial conditions $`\psi (0)=1/\gamma `$, $`\lambda (0)=ϵ\gamma `$. If we further define $`\mathrm{\Lambda }_{\mathrm{}}=_0^{\mathrm{}}\lambda (t)dt`$ and $`\mathrm{\Psi }_{\mathrm{}}=_0^{\mathrm{}}\psi (t)dt`$, we can integrate both LHS and RHS of the equations (51) and obtain the exact relation between $`\mathrm{\Lambda }_{\mathrm{}}`$ and $`\mathrm{\Psi }_{\mathrm{}}`$

$$ϵ=2\mathrm{\Lambda }_{\mathrm{}}2ϵ\gamma ^2\mathrm{\Psi }_{\mathrm{}}.$$
(52)
Knowing $`\mathrm{\Lambda }_{\mathrm{}}`$ would solve the problem, because $`\sigma (ϵ)=\mathrm{e}^\mathrm{\Lambda }_{\mathrm{}}`$. However, in addition to (52) we need some other condition. It can be established from the observation that the equation for $`\lambda (t)`$ can be formally solved in the form

$$\lambda (t)=ϵ\gamma \mathrm{exp}\left(2\gamma (1ϵ)t+2ϵ\gamma ^2_0^t\psi (t^{})dt^{}\right).$$
(53)
and because $`\psi (t)<0`$, we have $`_0^t\psi (t^{})dt^{}>\mathrm{\Psi }_{\mathrm{}}`$. Therefore, we can write the following upper bound

$$\sigma (ϵ)<\mathrm{exp}\left(\frac{1}{2}W_\mathrm{L}\left(\frac{ϵ\mathrm{e}^ϵ}{1ϵ}\right)\right)\sigma _{\mathrm{upper}}(ϵ)$$
(54)
where $`W_\mathrm{L}(x)`$ is the Lambert function defined by the equation $`W_\mathrm{L}(x)\mathrm{e}^{W_\mathrm{L}(x)}=x`$.

The leading term in the asymptotic behaviour of the Lambert function for large argument is $`W_\mathrm{L}(x)\mathrm{ln}x`$. If we use it as an approximation for calculating the asymptotic behaviour of $`\sigma (ϵ)`$, starting with (54) we finally get

$$\sigma (ϵ)\sqrt{\mathrm{e}}\sqrt{1ϵ}$$
(55)
which is compatible with the behaviour (50). Actually, we can see in the Fig. 7 that the approximation (55) fits very well the results from numerical integration and lies much closer than the exact upper bound (54). Thus, we conjecture that the formula (55) is in fact the correct asymptotic behaviour for $`ϵ1`$.

To sum up, in the regime $`ϵ>1`$, the initial conditions are irrelevant for long-time dynamics, as expected in the ergodic phase. On the other hand, for $`ϵ<1`$, we observe that the asymptotic value of the average coordinate depends on the initial condition for the noise, which is yet another signature of ergodicity breaking. The factor $`\sigma (ϵ)`$ measures the sensitivity to initial conditions: it vanishes in ergodic phase but remains non-zero in non-ergodic phase. So, it may be considered as a kind of order parameter. Close to the critical point $`ϵ=1`$ it approaches zero continuously as a power with critical exponent $`1/2`$.

## 6 Response to harmonic perturbation

Let us investigate now the response of the particle to the external driving force. Adding the additional term $`F(t)=F_0\mathrm{cos}(\omega t)`$ at the right hand side of Eq. (1), one can repeat, mutatis mutandis, all steps leading to the system of equations for the functions $`s_1(t)`$ to $`s_4(t)`$. We find that the equations (17)-(20) hold unchanged, while the influence of the external force modifies the expression (21) for $`\lambda (t)`$. Actually, in the present case one gets

$$\begin{array}{cc}\hfill \frac{\lambda (t)}{ϵ\gamma }& =12\gamma s_3(t)+2\gamma ^2s_4(t)+\hfill \\ & +\frac{1}{q^2}\left(\gamma \mathrm{e}^{\gamma t}_0^tF(t^{})\mathrm{e}^{\gamma t^{}}dt^{}F(t)\right)^2.\hfill \end{array}$$
(56)
We assumed $`𝖰_0=0`$ and $`𝖷_0=0`$ here.

First quantity to study is the response of the average coordinate. We find

$$\begin{array}{cc}\hfill 𝖷(t)=& 𝖰_0\mathrm{e}^{\gamma t}_0^t\mathrm{e}^{\mathrm{\Lambda }(t^{})+\gamma t^{}}dt^{}+\hfill \\ & +\mathrm{e}^{\gamma t}_0^tF(t^{})\mathrm{e}^{\gamma t^{}}dt^{}.\hfill \end{array}$$
(57)
Obviously enough, in the stationary regime the average coordinate oscillates around $`0`$ with the same frequency as the driving force $`F(t)`$. We arrive at the standard Debye-type dynamic susceptibility

$$\chi (tt^{})=\mathrm{\Theta }(tt^{})\mathrm{e}^{\gamma (tt^{})}.$$
(58)
Thus, the exact response in terms of the coordinate is linear. This behaviour also does not depend on the value of $`ϵ`$.

The situation becomes much more complicated when we turn to quantities, which depend non-linearly on the coordinate, especially the switching rate $`\lambda (t)`$. We solved numerically the set of equations (17) to (20) and (56). We can see in Fig. 8 the evolution of the function $`\lambda (t)`$ within the non-ergodic regime, with $`ϵ=0.8`$. Comparing the behaviour with the evolution in absence of the external harmonic perturbation we can see that the switching rate $`\lambda (t)`$ does not approach zero any more, but it oscillates around some finite value, which we will denote $`A_0`$. This qualitative feature holds for whatever small external force. However, when the amplitude $`F_0`$ of the external field goes to zero, also the value of $`A_0`$ approaches zero according to $`A_0F_{0}^{}{}_{}{}^{2}`$, as can be seen from Fig. 9. Generalising the linear stability analysis of Sec. 3.2 to harmonic oscillations, we conclude that the external field continuously shifts the fixed point with $`\lambda =0`$, (i. e. also $`A_0=0`$) to a position with positive $`A_0`$, but the value of the shift vanishes when the amplitude of the perturbation goes to zero.

Figure 10 exemplifies the evolution of the functions $`s_3(t)`$ and $`s_4(t)`$. We can clearly see that the oscillations are not harmonic. Generally, in the stationary regime these functions are non-harmonic but periodic with the doubled frequency $`2\omega `$. Thus, the same holds also for the function $`\lambda (t)`$. Using (56) the stationary response in terms of the switching rate can be written as Fourier series

$$\lambda _{\mathrm{st}}(t)=A_0+\underset{k=1}{\overset{\mathrm{}}{}}\left(A_k\mathrm{sin}2k\omega t+B_k\mathrm{cos}2k\omega t\right).$$
(59)
The amplitudes of the harmonic components $`A_0`$, $`A_k,B_k`$, $`k=1,2,\mathrm{}`$, satisfy a complicated infinite set of quadratic equations.

To asses the weight of the higher harmonics we performed the fast Fourier transform of the time evolutions obtained by numerical solution, throwing away the initial transient regime. To illustrate the presence of higher harmonics we chose the function $`s_3(t)`$. The modulus of its Fourier transform $`\widehat{s}_3(\nu )=s_3(t)\mathrm{e}^{2\pi \mathrm{i}\nu t}dt`$ is shown in Fig. 11. We can clearly see the peaks at the multiples of the basic frequency. We can also observe that the higher harmonics have quite considerable weight. In the inset of Fig. 11 we show also the Fourier transform of the function $`\lambda (t)`$. Here, the higher harmonics are much less pronounced.

The most important feature of the time evolution under the influence of external harmonic force is the observation that the switching rate remains always positive. This leads to the already mentioned fact that whatever is the coupling strength parameter $`ϵ`$, the mean coordinate in stationary regime oscillates around zero, irrespectively of the initial conditions. However, this is the signature of ergodicity, so the glassy behaviour disappears under the influence of arbitrarily small external perturbation. Such a behaviour was already observed also in the model of sheared colloid cates\_02 .

Both observations can be easily understood. The external force drags the particle back and forth. Once moving, the particle induces through the back-reaction (2) the fluctuations of the environmental force, which prevents the system from freezing in a non-ergodic state. This holds for any positive amplitude of the driving force. However, when the force diminishes, there is still longer transient period, where the system apparently relaxes toward the arrested state, as can be seen qualitatively in Fig. 8. In the limit of infinitesimally small driving, the transient time blows up and the asymptotic state corresponds to the dynamically arrested state. This picture is consistent with the view of glassiness as a purely dynamical phenomenon.

We also observed the response of the system to a signal, which is switched on only after the system relaxed very close to the arrested state. The results can be seen in Fig. 12. Initially, the switching rate relaxes toward zero, but after the perturbation it settles on oscillating behaviour. The average coordinate initially approaches non-zero value (we have chosen initial condition $`𝖰_0>0`$), but the perturbation brings it to oscillations around zero. This can be interpreted as a schematic picture of shear thinning, although the model is too much simplified to account for the shear thinning quantitatively.

## 7 Conclusions

Motivated by the dynamical arrest phenomenon in colloids we formulated and solved a stochastic dynamical model, where the coordinate of a single particle evolves under the influence of stochastic environmental force. The back-reaction couples the switching rate of the force to the average of the square of the velocity of the particle. The strength of the coupling $`ϵ`$ is the crucial parameter which determines the behaviour of the system. The problem reduces to a set of coupled non-linear differential equations, which was investigated both analytically and numerically.

The back-reaction induces a phase transition from the ergodic phase for $`ϵ>1`$ to the non-ergodic glassy state for $`ϵ<1`$. The transition is observed qualitatively in the behaviour of the switching rate, decaying to zero in non-ergodic state, while staying positive in the ergodic state. The Edwards-Anderson parameter, established from the two-time correlation functions, is discontinuous at the transition; it is zero in ergodic phase ($`ϵ>1`$), while for $`ϵ<1`$ it acquires finite value independent of $`ϵ`$. The critical point $`ϵ=1`$ is characterised by power-law decay of the switching rate. The leading term $`t^1`$ in the long-time behaviour was calculated analytically.

We investigated the critical behaviour at $`ϵ1`$ through the dependence of the average coordinate in the long-time limit on the initial condition. We find that the asymptotic value of the average coordinate is proportional to the average initial value of the force, where the proportionality factor $`\sigma (ϵ)`$ is singular at the transition. In the ergodic phase we have $`\sigma (ϵ)=0`$ identically, while in the non-ergodic phase, close to the transition, we found $`\sigma (ϵ)(1ϵ)^{1/2}`$ for $`ϵ1^{}`$.

Therefore, we find a situation quite unusual from the point of view of static equilibrium phase transition. Indeed, we have two variables, which may be considered as order parameter, namely the Edwards-Anderson parameter and the quantity $`\sigma (ϵ)`$. While the former is discontinuous at the critical point, thus indicating first-order transition, the latter is continuous, suggesting second-order transition. the discrepancy is to be attributed to purely dynamical nature of the transition.

Finally, we investigated the response of the system to harmonic external perturbation. We found that the exact response of the coordinate is linear. On the other hand, the response of the variables which are quadratic functions of the coordinate and velocity, like the switching rate, was non-linear and generically contains all higher harmonics, as was seen in the Fourier transform of the signal. We also observed that arbitrarily weak external perturbation is sufficient to “melt” the non-ergodic glassy state and bring it back to ergodic behaviour. We may relate this feature to the notion of stochastic stability parisi\_00 ; in this view our system is not stochastically stable. However, our finding is in accord with previously observed behaviour of sheared colloids cates\_02 . It makes also connection to the rheological properties of thixotropic fluids vol\_nit\_hey\_reh\_02 ; cou\_ngu\_huy\_bon\_02 , although our model is too simplified to give quantitative predictions in this direction.

The back-reaction mechanism described by (2) represent the simplest choice. One may ask what would happen if we tried another prescription. We expect that the methods used here will be as well applicable if we generalise (2) as $`\lambda (t)=F(𝖵^2(t))`$ for an analytic function $`F(x)`$. More complicated situation would appear if the dependence was non-local in time, e. g. of the form $`\lambda (t)=^t𝖵(t)𝖵(t^{}))K(tt^{})\mathrm{d}t^{}`$ with some kernel $`K(t)`$. Such an approach would bring our model closer to the well-studied Mode Coupling equations, but it goes beyond the scope of the present work.

###### Acknowledgements.

We want to dedicate this paper to the memory of Prof. Vladislav Čápek, a passionate theoretical physicist, our dedicated teacher and good friend, who passed away shortly before this work was completed. This work was supported by the project No. 202/00/1187 of the Grant Agency of the Czech Republic.


format: text

text:
# 1 Introduction

## 1 Introduction

A reliable calculation of long-distance contributions to non-leptonic kaon-decay rates, and, in particular, to the CP-violating part parametrized by the quantity $`\epsilon ^{}/\epsilon `$ has been a longstanding challenge. Ideally, one would expect such calculations to be in the domain of lattice QCD, but in practice many theoretical and numerical difficulties have made progress in this direction rather slow. Recently however, two lattice collaborations have reported on numerical results for both the real and imaginary parts of $`\mathrm{\Delta }I=1/2`$ and $`\mathrm{\Delta }I=3/2`$ $`K\pi \pi `$ matrix elements with a rather satisfactory control over statistical errors . These lattice computations were done with the effective weak $`\mathrm{\Delta }S=1`$ hamiltonian with three flavors, i.e. with the charm integrated out, and they were possible because of the use of lattice fermions with good chiral symmetry. Both groups reported values of $`\epsilon ^{}/\epsilon `$ which are non-zero, and thus consistent with the existence of direct CP violation, but of opposite sign (and comparable size) to the experimentally measured value.

While statistical errors for these lattice computations seem to be reasonably under control, this is not the case for a large class of systematic errors, which will need to be studied further in the future. One source of systematic error is the use of the quenched approximation. In a previous paper  we pointed out that, in addition to the fact that quenched QCD is just not the same theory as full QCD, an ambiguity arises in the definition of the quenched version of penguin operators appearing in the $`\mathrm{\Delta }S=1`$ effective weak hamiltonian. The ambiguity originates in the difference of the chiral transformation properties of penguin operators within the quenched and unquenched theories. This implies that not only will quenched lattice results be hampered by the fact that we do not really know whether quenched values of given matrix elements are close to their real-world values, but, in the case of penguins, they also depend on which definition of the operators is chosen, since more than one definition is possible.

In fact, both lattice computations  did not directly compute $`K\pi \pi `$ matrix elements, but $`K\pi `$ (with $`M_\pi =M_K`$) and $`K0`$ transition amplitudes, and used chiral perturbation theory (ChPT) to convert them into the desired $`K\pi \pi `$ matrix elements . In ref.  we explained how the chiral properties of penguin operators change in the transition to the (partially) quenched theory, and how, in principle, more than one definition of a quenched penguin operator is possible. Using ChPT, we traced how this affects $`K\pi `$ and $`K0`$ matrix elements. We restricted ourselves to $`LR`$ penguins (i.e. $`Q_5`$ and $`Q_6`$), because the effects in this case already appear at leading order in ChPT, while they are a next-to-leading order effect for $`LL`$ penguins. Because of the fact that the ambiguity is already present at leading order for matrix elements of $`Q_{5,6}`$, this may be an important issue for $`\epsilon ^{}/\epsilon `$ (while it is expected to be less important for the real parts of $`K\pi \pi `$ amplitudes, and thus the $`\mathrm{\Delta }I=1/2`$ rule).

In this paper we extend our ChPT calculations to the effect of the ambiguity on $`K\pi \pi `$ matrix elements, again to leading chiral order. This is important for two reasons. First, lattice computations may be done directly for $`K\pi \pi `$ matrix elements, and their chiral behavior needs to be known in order to fit lattice results as a function of quark masses. Since lattice computations are typically done with unphysical (i.e. energy/momentum non-conserving) kinematics, we present our results for the most general kinematics possible, both in the quenched and partially quenched cases. Second, once the complete (leading-order) ChPT expressions for $`K0`$, $`K\pi `$ and $`K\pi \pi `$ matrix elements are available, it is possible to give a more detailed discussion of the systematic error introduced by the ambiguity in the definition of quenched penguin operators.

The paper is organized as follows. In section 2, we review the main observation of ref. . We show how a $`LR`$ penguin, which transforms in an irreducible representation (irrep) of $`SU(3)_L\times SU(3)_R`$, splits into two operators in the (partially) quenched theory, with each of them transforming in a different irrep of the (partially) quenched chiral symmetry group. One of these irreps corresponds “naturally” to the single irrep of the unquenched theory, while the other irrep can be considered as “new,” and an artifact of quenching. We give the ChPT realization of all relevant operators at leading and next-to-leading chiral order, introducing new low-energy constants (LECs) which appear in correspondence to the new irrep. In section 3, we present our results for $`K^0\pi ^+\pi ^{}`$ penguin matrix elements to leading order in ChPT, with general kinematics, and specialize these results to physical (i.e. energy-momentum conserving) kinematics. In section 4, we discuss different strategies available for using quenched lattice results to estimate real-world $`K\pi \pi `$ matrix elements, and give some numerical examples. We provide a simple prescription for determining the leading-order LEC representing the new irrep in section 5, and section 6 contains our conclusions. Some of this work has already been presented in ref. .

## 2  Review of $`LR`$ penguins in (partially) quenched QCD and ChPT

A lagrangian definition for partially quenched QCD can be constructed as explained in ref.  (see also ref.  for an alternative realization using the replica method). In addition to the valence quarks $`q_{vi}`$, $`i=u,d,s`$, with masses $`m_{vi}`$, one introduces a separate set of sea quarks $`q_{si}`$, $`i=1,\mathrm{},N`$, with masses $`m_{si}`$, and a set of “ghost” quarks $`q_{gi}`$, $`i=u,d,s`$, with masses equal to those of the valence quarks $`m_{gi}=m_{vi}`$ . Ghost quarks are given bosonic statistics, such that the ghost-quark determinant cancels the valence-quark determinant, thus leaving only the sea-quark determinant present in the path integral. Therefore, only sea quarks propagate in internal loops.

Since partially quenched QCD thus contains more flavors than unquenched QCD, its flavor symmetry group is larger than the QCD one. The full chiral symmetry group relevant for light meson physics is the graded extension of the ordinary chiral group $`SU(3+N|3)_L\times SU(3+N|3)_R`$ . It is graded because part of its elements transform fermions into bosons and vice versa. The quenched theory, which has no sea quarks at all, corresponds to the special case $`N=0`$ .

We consider $`LR`$ penguin operators of the form

$$Q_{penguin}=(\overline{s}d)_L(\overline{u}u+\overline{d}d+\overline{s}s)_R,$$
(2.1)
where

$`(\overline{q}_1q_2)_{L,R}`$ $`=`$ $`\overline{q}_1\gamma _\mu P_{L,R}q_2,`$ (2.2)
$`P_{L,R}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1\gamma _5),`$
and color contractions are not specified, so that $`Q_{penguin}`$ can represent both $`Q_5`$ and $`Q_6`$. As already pointed out in ref. , the $`u`$, $`d`$ and $`s`$ fields in eq. (2.1) represent valence quarks in the partially quenched theory, and the penguin operator can be decomposed as111In a theory with $`K`$ valence quarks, all ratios $`3/N`$ get replaced by $`K/N`$.

$`Q_{penguin}`$ $`=`$ $`{\displaystyle \frac{3}{N}}\text{str}(\mathrm{\Lambda }\psi \overline{\psi }\gamma _\mu P_L)\text{str}(\psi \overline{\psi }\gamma _\mu P_R)+\text{str}(\mathrm{\Lambda }\psi \overline{\psi }\gamma _\mu P_L)\text{str}(A\psi \overline{\psi }\gamma _\mu P_R),`$ (2.3)
$``$ $`{\displaystyle \frac{3}{N}}Q_{penguin}^{PQS}+Q_{penguin}^{PQA},`$
$`A`$ $`=`$ $`\mathrm{diag}(1{\displaystyle \frac{3}{N}},1{\displaystyle \frac{3}{N}},1{\displaystyle \frac{3}{N}},{\displaystyle \frac{3}{N}},\mathrm{},{\displaystyle \frac{3}{N}}),`$ (2.4)
$`\mathrm{\Lambda }_{ij}`$ $`=`$ $`\delta _{is}\delta _{jd},`$
where the first $`3`$ (valence) entries of $`A`$ are equal to $`13/N`$, and the next $`N+3`$ (sea and ghost) entries are equal to $`3/N`$. Here $`\psi `$ collects all quark fields in the theory, $`\psi =(q_{vi},q_{si},q_{gi})`$. The tensor $`\mathrm{\Lambda }`$ projects onto the valence $`(\overline{s}d)_L`$ term in the first factor of $`Q_{penguin}`$. The motivation for splitting $`Q_{penguin}`$ this way is that $`Q_{penguin}^{PQS}`$ and $`Q_{penguin}^{PQA}`$ form different representations of the partially quenched symmetry group: $`Q_{penguin}^{PQS}`$ ($`Q_{penguin}^{PQA}`$) transforms in the trivial (adjoint) irrep of $`SU(3+N|3)_R`$. As a consequence, there are at least two different ways of embedding the QCD penguin operator into the partially quenched theory. One is to choose the partially quenched penguin to be a singlet under $`SU(3+N|3)_R`$, i.e. $`Q_{penguin}^{PQS}`$, as in the unquenched theory, whereas the other choice is to use the original operator, which is seen to be a linear combination of two irreducible operators. This latter choice was made in refs. . In ref.  non-singlet penguin operators such as $`Q_{penguin}^{PQA}`$ were not considered, because singlet factors such as $`(\overline{u}u+\overline{d}d+\overline{s}s)_R`$ in eq. (2.1) had been implicitly extended to singlets under the full (partially-)quenched symmetry group. Therefore, the analysis of ref.  was not complete, and ref.  and this paper remedy this for LR penguins.

To leading order (LO), the operators $`Q_{penguin}^{PQS}`$, $`Q_{penguin}^{PQA}`$ are represented in ChPT by

$`Q_{penguin}^{PQS}`$ $``$ $`\alpha _1^{(8,1)}\text{str}(\mathrm{\Lambda }L_\mu L_\mu )+\alpha _2^{(8,1)}\text{str}(\mathrm{\Lambda }X_+),`$ (2.5)
$`Q_{penguin}^{PQA}`$ $``$ $`f^2\alpha ^{(8,8)}\text{str}(\mathrm{\Lambda }\mathrm{\Sigma }A\mathrm{\Sigma }^{}),`$ (2.6)
where

$$L_\mu =i\mathrm{\Sigma }_\mu \mathrm{\Sigma }^{},X_\pm =2B_0(\mathrm{\Sigma }M^{}\pm M\mathrm{\Sigma }^{}),$$
(2.7)
with $`M`$ the quark-mass matrix, $`B_0`$ the parameter $`B_0`$ of ref. , $`\mathrm{\Sigma }=\text{exp}(2i\mathrm{\Phi }/f)`$ the unitary field describing the partially quenched Goldstone-meson multiplet, and $`f`$ the bare pion-decay constant normalized such that $`f_\pi =132`$ MeV. The $`\alpha `$’s are the corresponding LECs. Notice that $`Q_{penguin}^{PQA}`$, unlike $`Q_{penguin}^{PQS}`$, is of order $`p^0`$, due to the fact that the right-handed current in $`Q_{penguin}^{PQA}`$ is not a partially quenched singlet (cf. electro-magnetic penguins222In fact, $`Q_{penguin}^{PQA}`$ is a component of the same irrep as the electro-magnetic penguin, except for $`N=0`$ .). As already observed in ref. , the new operator $`Q_{penguin}^{PQA}`$ does not contribute at tree level to matrix elements with only valence quarks on external lines, since the matrix $`A`$ is effectively proportional to the unit matrix in the valence sector. Indeed, replacing $`A`$ by the unit matrix in eq. (2.6) would make the operator vanish. This is no longer true at next-to-leading-order (NLO), i.e. at order $`p^2`$, where one-loop contributions from $`Q_{penguin}^{PQA}`$ to valence-quark matrix elements are non-zero.

Since the singlet-operator contributions also start at order $`p^2`$, the NLO contributions from $`Q_{penguin}^{PQA}`$ compete with the LO contributions from $`Q_{penguin}^{PQS}`$, and thus need to be taken into account already in a leading-order analysis of $`K0`$, $`K\pi `$ and $`K\pi \pi `$ matrix elements. This also implies that a renormalization scale dependence already appears at leading order for those partially quenched matrix elements. That scale dependence is absorbed by new $`O(p^2)`$ counterterms for $`Q_{penguin}^{PQA}`$. The complete list of $`CPS`$-even  operators is

$`Q_1^{PQA}`$ $`=`$ $`{\displaystyle \frac{\beta _1^{(8,8)}}{(4\pi )^2}}\text{str}(\mathrm{\Lambda }\{\mathrm{\Sigma }A\mathrm{\Sigma }^{},L_\mu L_\mu \}),`$ (2.8)
$`Q_2^{PQA}`$ $`=`$ $`{\displaystyle \frac{\beta _2^{(8,8)}}{(4\pi )^2}}\text{str}(\mathrm{\Lambda }L_\mu \mathrm{\Sigma }A\mathrm{\Sigma }^{}L_\mu ),`$
$`Q_3^{PQA}`$ $`=`$ $`{\displaystyle \frac{\beta _3^{(8,8)}}{(4\pi )^2}}\text{str}(\mathrm{\Lambda }\{\mathrm{\Sigma }A\mathrm{\Sigma }^{},X_+\}),`$
$`Q_4^{PQA}`$ $`=`$ $`{\displaystyle \frac{\beta _4^{(8,8)}}{(4\pi )^2}}\text{str}(\mathrm{\Lambda }[\mathrm{\Sigma }A\mathrm{\Sigma }^{},X_{}]),`$
$`Q_5^{PQA}`$ $`=`$ $`{\displaystyle \frac{\beta _5^{(8,8)}}{(4\pi )^2}}\text{str}(\mathrm{\Lambda }\mathrm{\Sigma }A\mathrm{\Sigma }^{})\text{str}(L_\mu L_\mu ),`$
$`Q_6^{PQA}`$ $`=`$ $`{\displaystyle \frac{\beta _6^{(8,8)}}{(4\pi )^2}}\text{str}(\mathrm{\Lambda }\mathrm{\Sigma }A\mathrm{\Sigma }^{})\text{str}(X_+),`$
$`Q_7^{PQA}`$ $`=`$ $`{\displaystyle \frac{\beta _7^{(8,8)}}{(4\pi )^2}}i_\mu \text{str}(\mathrm{\Lambda }[\mathrm{\Sigma }A\mathrm{\Sigma }^{},L_\mu ]),`$
where we introduced the $`O(p^2)`$ LECs $`\beta _{1,\mathrm{},7}^{(8,8)}`$.

The partially quenched theory with $`N=3`$ light sea quarks represents a special case. For $`N=3`$, the LECs of the partially quenched theory must be the same as those of the physical, unquenched theory , basically because they represent the coefficients in an expansion in powers of quark masses, and thus only depend on the number of dynamical (sea) quarks, and not on the quark masses themselves. Therefore, in the $`N=3`$ partially quenched theory, what one should do is to omit $`Q_{penguin}^{PQA}`$ altogether, because the aim is to obtain the values of $`\alpha _{1,2}^{(8,1)}`$, and not the $`(8,8)`$ LECs $`\alpha ^{(8,8)}`$ and $`\beta _i^{(8,8)}`$.333As long as one does not consider electro-magnetic penguin contributions.

The quenched case, $`N=0`$, is different. In this case, the decomposition reads

$`Q_{penguin}^{QCD}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\text{str}(\mathrm{\Lambda }\psi \overline{\psi }\gamma _\mu P_L)\text{str}(\psi \overline{\psi }\gamma _\mu P_R)+\text{str}(\mathrm{\Lambda }\psi \overline{\psi }\gamma _\mu P_L)\text{str}(\widehat{N}\psi \overline{\psi }\gamma _\mu P_R),`$ (2.9)
$``$ $`{\displaystyle \frac{1}{2}}Q_{penguin}^{QS}+Q_{penguin}^{QNS},`$
$`\widehat{N}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{diag}(1,1,1,1,1,1),`$ (2.10)
where valence entries of $`\widehat{N}`$ are equal to $`\frac{1}{2}`$, and ghost entries are equal to $`\frac{1}{2}`$. The first operator in the decomposition is a singlet under $`SU(3|3)_R`$, while the second is not ($`NS`$ for non-singlet). However, $`Q_{penguin}^{QNS}`$ can transform into the singlet operator, implying that the non-singlet operators do not form a representation by themselves. In other words, $`Q_{penguin}^{QNS}`$ can mix into $`Q_{penguin}^{QS}`$, which is possible because $`\widehat{N}`$ is not supertrace-less, unlike $`A`$ in the partially quenched case. The ChPT realization of both operators is obtained from the expressions given in eqs. (2.6,2.8) and by replacing $`A\widehat{N}`$. When referring to the quenched theory, we will add a subscript $`q`$ to the LECs, and rename444In the quenched theory there is no relation between $`Q_{penguin}^{QNS}`$ and the electro-magnetic penguins . the LECs for $`Q_{penguin}^{QNS}`$ as $`\alpha ^{(8,8)}\alpha _q^{NS}`$, $`\beta _i^{(8,8)}\beta _{qi}^{NS}`$.

Within the fully quenched approximation, as also in the partially quenched case with $`N3`$, there is no reason that the LECs should have the same values as those of the unquenched theory. In general, the non-analytic terms are modified by quenching, and even the scale dependence of LECs is different between the quenched and unquenched theories. It is thus not a priori clear what choice to make for the embedding of penguin operators into the quenched theory. We will return to this issue in section 4 below.

## 3 Non-singlet $`K^\mathrm{𝟎}\mathbf{}\pi ^\mathbf{+}\pi ^{\mathbf{}}`$ matrix elements with general kinematics

In this section we present the partially quenched and quenched results for the contribution of strong penguin operators to the $`K^0\pi ^+\pi ^{}`$ matrix element to order $`p^2`$ in ChPT. We will restrict ourselves to the isospin limit in the valence sector, $`m_{vu}=m_{vd}`$, not assuming momentum conservation, so that $`qp_1+p_2`$, with $`q`$ the (ingoing) $`K^0`$ momentum and $`p_1`$ ($`p_2`$) the (outgoing) $`\pi ^+`$ ($`\pi ^{}`$) momentum. All momenta are onshell and we work in euclidean space, i.e. $`p_1^2=p_2^2=M_\pi ^2`$, $`q^2=M_K^2`$. $`M_{jvi}`$ ($`M_{jsi}`$) is the mass of a meson made out of the $`j`$th valence quark and the $`i`$th valence (sea) quark. In all results presented below, we have symmetrized the expressions in the pion momenta $`p_1`$ and $`p_2`$.

For the contribution of the singlet operator in eq. (2.5) a simple tree-level calculation yields the order $`p^2`$ result

$`\pi ^+\pi ^{}|Q^{PQS}|K^0`$ $`=`$ $`{\displaystyle \frac{4i}{f^3}}\{\alpha _1^{(8,1)}({\displaystyle \frac{1}{2}}q(p_1+p_2)+p_1p_2)`$
$`+{\displaystyle \frac{2}{3}}\alpha _2^{(8,1)}(M_K^2M_\pi ^2)(1+{\displaystyle \frac{1}{2}}{\displaystyle \frac{2p_1p_2+q(p_1+p_2)2M_\pi ^2}{(qp_1p_2)^2+M_K^2}})\}.`$
Note that the contribution proportional to $`\alpha _2^{(8,1)}`$ vanishes for energy-conserving kinematics, because the corresponding operator can be written as a total derivative by using equations of motion . As mentioned before, the tree-level contribution from the new non-singlet operator in eq. (2.6) vanishes. At one loop, it contributes at order $`p^2`$. For general kinematics, we obtain, in the isospin limit,

$`\pi ^+\pi ^{}|Q^{PQA}|K^0_{1loop}`$ $`=`$ $`{\displaystyle \frac{1}{2}}[𝒜(q;p_1,p_2)+𝒜(q;p_2,p_1)],`$ (3.2)
$`𝒜(q;p_1,p_2)`$ $`=`$ $`{\displaystyle \frac{4i}{3f^3}}\alpha ^{(8,8)}{\displaystyle \underset{isea}{}}\{3(p_1p_2M_\pi ^2)I(M_{2si}^2,M_{2si}^2,(p_1+p_2)^2)`$
$`+{\displaystyle \frac{2M_K^2M_\pi ^2+qp_1}{(qp_1)^2}}\left(L(M_{2si}^2)L(M_{3si}^2)\right)`$
$`6{\displaystyle \frac{(qp_1+M_K^2)(qp_1+M_\pi ^2)}{(qp_1)^2}}I(M_{2si}^2,M_{3si}^2,(qp_1)^2)`$
$`+{\displaystyle \frac{(2p_1p_22qp_1+4qp_22M_\pi ^2)}{(qp_1p_2)^2+M_K^2}}(L(M_{3si}^2)L(M_{2si}^2))\}`$
$`{\displaystyle \underset{ighost}{}}(M_{jsi}M_{jvi},j=2,3),`$
where we used that $`M_{3si}^2M_{1si}^2=M_K^2M_\pi ^2`$ at this order, and

$`L(M^2)`$ $`=`$ $`{\displaystyle \frac{d^D\mathrm{}}{(2\pi )^D}\frac{1}{\mathrm{}^2+M^2}},`$ (3.3)
$`I(M_1^2,M_2^2,p^2)`$ $`=`$ $`{\displaystyle \frac{d^D\mathrm{}}{(2\pi )^D}\frac{1}{(\mathrm{}^2+M_1^2)((\mathrm{}p)^2+M_2^2)}},`$
where $`D`$ is the number of space-time dimensions. Explicit expressions for these integrals in the $`\overline{MS}`$ scheme are given in the appendix. The one-loop diagrams contributing to this result are displayed in fig. (1). In particular, the vacuum-tadpole diagram (last diagram in fig. (1)) gives a non-zero contribution for general kinematics which corresponds to the next-to-last line of eq. (3.2).

In order to obtain the fully quenched result with no sea quarks at all, one simply drops the sum over sea quarks in eq. (3.2), keeping only the sum over ghost quarks, and replaces $`\alpha ^{(8,8)}\alpha _q^{NS}`$. There are no contributions from $`\eta ^{}`$ double-poles.

The new $`O(p^2)`$ counterterms in eq. (2.8) give, again for general kinematics

$`\pi ^+\pi ^{}|Q^{PQA}|K^0_{ct}`$ $`=`$ $`(1{\displaystyle \frac{3}{N}}){\displaystyle \frac{4i}{(4\pi ^2)f^3}}\{(2\beta _1^{(8,8)}+\beta _2^{(8,8)})({\displaystyle \frac{1}{2}}q(p_1+p_2)+p_1p_2)`$ (3.4)
$`{\displaystyle \frac{4}{3}}\beta _3^{(8,8)}(M_K^2M_\pi ^2)(1+{\displaystyle \frac{1}{2}}{\displaystyle \frac{2p_1p_2+q(p_1+p_2)2M_\pi ^2}{(qp_1p_2)^2+M_K^2}})\}.`$
The structure of these contributions is such that the scale dependence contained in eq. (3.2) can be fully compensated by the LECs $`\beta _i^{(8,8)}`$, $`i=1,2,3`$. It is easy to verify that the other LECs $`\beta _i^{(8,8)}`$, $`i=4,\mathrm{},7`$ do not contribute. Since external legs only contain valence quarks, the tree-level $`K\pi \pi `$ matrix element does not involve the diagonal $`A_{ii}`$ elements of the tensor $`A`$ with $`i`$ referring to a sea or ghost quark. This implies that, for this calculation, we may replace $`A`$ with the unit matrix, and the operators $`Q_i^{PQA}`$ in eq. (2.5) for $`i=4,\mathrm{},7`$ vanish, while they become proportional to the operators in eq. (2.5) for $`i=1,2,3`$. It follows that the LECs $`\alpha _{1,2}^{(8,1)}`$ together with $`\beta _{1,2,3}^{(8,8)}`$ will always appear in the specific combinations

$`\alpha _1^{(8,1)}\left(1{\displaystyle \frac{3}{N}}\right){\displaystyle \frac{1}{(4\pi )^2}}(2\beta _1^{(8,8)}+\beta _2^{(8,8)}),`$ (3.5)
$`\alpha _2^{(8,1)}+\left(1{\displaystyle \frac{3}{N}}\right){\displaystyle \frac{2}{(4\pi )^2}}\beta _3^{(8,8)}`$
for all tree-level matrix elements with only valence quarks on the external legs. For the analogue of eq. (3.5) in the quenched case, one replaces the factor $`13/N`$ by $`1/2`$, $`\alpha _i^{(8,1)}\alpha _{qi}^{(8,1)}`$, $`\alpha ^{(8,8)}\alpha _q^{NS}`$ and $`\beta _i^{(8,8)}\beta _{qi}^{NS}`$.

We conclude this section with the expressions for the same matrix elements in the case of “physical” kinematics, i.e. with the choice $`q=p_1+p_2`$. Setting $`q=p_1+p_2`$ and using that $`q^2=M_K^2`$ and $`p_1^2=p_2^2=M_\pi ^2`$, one obtains from eqs. (3), (3.2) and (3.4)

$`\pi ^+\pi ^{}|Q^{PQS}|K^0^{phys}`$ $`=`$ $`{\displaystyle \frac{4i}{f^3}}\alpha _1^{(8,1)}(M_K^2M_\pi ^2),`$ (3.6)
$`\pi ^+\pi ^{}|Q^{PQA}|K^0_{1loop}^{phys}`$ $`=`$ $`{\displaystyle \frac{4i}{3f^3}}\alpha ^{(8,8)}{\displaystyle \underset{isea}{}}\{{\displaystyle \frac{3}{2}}M_K^2I(M_{2si}^2,M_{2si}^2,M_K^2)`$
$`\left({\displaystyle \frac{3}{2}}{\displaystyle \frac{M_K^4}{M_\pi ^2}}3M_K^2\right)I(M_{2si}^2,M_{3si}^2,M_\pi ^2)`$
$`({\displaystyle \frac{3}{2}}{\displaystyle \frac{M_K^2}{M_\pi ^2}}3)[L(M_{2si}^2)L(M_{3si}^2)]\}`$
$`{\displaystyle \underset{ighost}{}}(M_{jsi}M_{jvi},j=2,3),`$
$`\pi ^+\pi ^{}|Q^{PQA}|K^0_{ct}^{phys}`$ $`=`$ $`{\displaystyle \frac{4i}{(4\pi ^2)f^3}}\left(1{\displaystyle \frac{3}{N}}\right)\left(2\beta _1^{(8,8)}+\beta _2^{(8,8)}\right)(M_K^2M_\pi ^2).`$ (3.8)
Dropping the sea-quark terms and carrying out the sum over ghost quarks, we obtain a more explicit expression for the one-loop contribution in the quenched case, in the $`\overline{MS}`$ scheme:

$`\pi ^+\pi ^{}|Q^{QNS}|K^0_{1loop}^{phys}`$ $`=`$
$`{\displaystyle \frac{i}{16\pi ^2f^3}}\alpha _q^{NS}\{12(M_K^2M_\pi ^2)(\mathrm{log}{\displaystyle \frac{M_\pi ^2}{\mathrm{\Lambda }^2}}1)`$
$`+\left({\displaystyle \frac{M_K^6}{M_\pi ^4}}2{\displaystyle \frac{M_K^4}{M_\pi ^2}}+2M_K^2\right)\mathrm{log}{\displaystyle \frac{M_K^2}{M_\pi ^2}}`$
$`+\left({\displaystyle \frac{M_K^6}{M_\pi ^4}}6{\displaystyle \frac{M_K^4}{M_\pi ^2}}+10M_K^24M_\pi ^2\right)\mathrm{log}{\displaystyle \frac{2M_K^2M_\pi ^2}{M_\pi ^2}}`$
$`+2M_K^2\left(F(M_\pi ^2,M_\pi ^2,M_K^2)2i\pi \theta (M_K^24M_\pi ^2)\sqrt{1{\displaystyle \frac{4M_\pi ^2}{M_K^2}}}+{\displaystyle \frac{\pi }{3}}\sqrt{3}\right)`$
$`+({\displaystyle \frac{M_K^4}{M_\pi ^2}}2M_K^2)(2F(M_\pi ^2,M_K^2,M_\pi ^2)+F(M_K^2,2M_K^2M_\pi ^2,M_\pi ^2))\},`$
where the function $`F`$ is given in the appendix.

## 4 Strategies for quenched estimates of real-world penguin matrix elements

Recent numerical estimates of $`K\pi \pi `$ matrix elements reported in refs.  have been obtained via the indirect method, where the simpler $`K\pi `$ and $`K0`$ transition amplitudes are computed on the lattice and then converted into estimates for $`K\pi \pi `$ matrix elements using ChPT. However, the fact that those numerical results are still obtained in the quenched approximation introduces a source of systematic error which is in principle uncontrolled. As already explained, the LECs of the quenched theory do not have to have values equal to those of the unquenched theory. Typically, even the scale dependence of quenched and unquenched LECs is not the same; it depends on the number of light dynamical (sea) quarks in the theory.

In the case of penguin operators, an additional ambiguity arises because it is a priori unclear whether it would be best to take matrix elements of $`Q_{penguin}^{QCD}`$, i.e. a linear combination of $`Q_{penguin}^{QS}`$ and $`Q_{penguin}^{QNS}`$ as in eq. (2.9), or to drop the contribution from $`Q^{QNS}`$ under the assumption that $`\alpha _{q1}^{(8,1)}`$ is the best estimate of $`\alpha _1^{(8,1)}`$.

In order to discuss possible strategies in more detail, we first recall the leading order ChPT expressions for $`K\pi `$ and $`K0`$ matrix elements of strong penguin operators, from ref. . In the quenched approximation ($`N=0`$) one has

$`\pi ^+|Q_{penguin}^{QCD}|K^+`$ $`=`$ $`{\displaystyle \frac{4M^2}{f^2}}\left\{\alpha _{q1}^{(8,1)}\alpha _{q2}^{(8,1)}{\displaystyle \frac{1}{(4\pi )^2}}\left(\beta _{q1}^{NS}+{\displaystyle \frac{1}{2}}\beta _{q2}^{NS}+\beta _{q3}^{NS}\right)\right\},`$ (4.1)
$`0|Q_{penguin}^{QCD}|K^0`$ $`=`$ $`{\displaystyle \frac{4i}{f}}\{(\alpha _{q2}^{(8,1)}+{\displaystyle \frac{1}{(4\pi )^2}}\beta _{q3}^{NS})(M_K^2M_\pi ^2)`$
$`+\alpha _q^{NS}{\displaystyle \underset{ivalence}{}}(L(M_{3vi}^2)L(M_{2vi}^2))\},`$
where $`M_K=M_\pi =M`$ in the case of the $`K\pi `$ matrix element, and contributions of both singlet and non-singlet operators are included. Notice that only the (quenched versions of the) combinations (3.5) of LECs appear in these expressions, as expected. Assuming that one can limit the analysis to leading-order in ChPT, there are at least three different strategies for estimating $`K\pi \pi `$ penguin matrix elements from LECs obtained by fitting eqs. (4.1,4) to quenched numerical results:

* Ignore $`\alpha _q^{NS}`$, but not the other LECs associated with the non-singlet operator, $`\beta _{q1,2,3}^{NS}`$. Both $`\beta _{q3}^{NS}`$ and $`\beta _{q1}^{NS}+\frac{1}{2}\beta _{q2}^{NS}`$ are scale dependent (however, their sum is not, as can be seen from eq. (4.1)), implying that this strategy is scale dependent. However, it still makes sense in case the non-analytic contribution proportional to $`\alpha _q^{NS}`$ (cf. eq. (3)) is numerically small compared to all other contributions at a reasonable scale $`\mathrm{\Lambda }`$ of order 1 GeV. Thus, the linear combination $`\alpha _{q1}^{(8,1)}\frac{1}{(4\pi )^2}(\beta _{q1}^{NS}+\frac{1}{2}\beta _{q2}^{NS})`$ is taken as the best estimate for the unquenched $`\alpha _1^{(8,1)}`$, and eq. (3.6) can then be used to obtain the physical $`K\pi \pi `$ matrix element (at tree level). This is the strategy followed in refs. . In fact, in these works, it was assumed that the contribution proportional to $`\alpha _q^{NS}`$ in eq. (4) is small.
* Drop all the non-singlet operators. It was shown in ref.  that this can be done by dropping, in the fully quenched case, all eye-diagrams in which the right-handed quarks in eq. (2.1) are contracted. This can be easily deduced from eq. (2.9). This strategy was explored for $`Q_6`$ in ref. . (For the partially-quenched case, see below.)
* Perform a complete quenched calculation including all contributions from singlet and non-singlet operators. After extracting all the LECs, singlet and non-singlet, one can use the sum of eqs. (3.6), (3) and the quenched version of eq. (3.8) to determine the quenched $`K\pi \pi `$ matrix element at the physical point.

Strategy 2 isolates $`\alpha _{q1}^{(8,1)}`$, and might thus appear to be the obvious choice, since it is this LEC that is needed for calculating $`K\pi \pi `$ matrix element (to chiral leading-order) in the unquenched theory. However, as we already mentioned, the values of LECs in the quenched and unquenched theories do not have to be equal, and it might happen that (at some scale $`\mathrm{\Lambda }`$) the quenched combination $`\alpha _{q1}^{(8,1)}\frac{1}{(4\pi )^2}(\beta _{q1}^{NS}+\frac{1}{2}\beta _{q2}^{NS})`$, determined from strategy 1, is indeed a better estimate of $`\alpha _1^{(8,1)}`$. Strategy 2 can be viewed as the situation in which the strong interactions are quenched at all scales between the weak and hadronic scales, because in that case only singlet penguin operators would appear in the evolution from the weak to the hadronic scale. So, while, on the one hand, it appears natural to assume only a mild flavor dependence of the LECs, in particular $`\alpha _1^{(8,1)}`$, one might, on the other hand, argue that it is better to calculate the evolution from the weak to the hadronic scale in the unquenched theory, even if the matrix element at the hadronic scale is finally computed in the quenched approximation. The key point is that it is impossible to decide which strategy is best.

The exception to these observations is the case of partially quenched QCD in which the number of light sea quarks is equal to that of the real world, in which $`N=3`$. In the partially quenched theory, the singlet operator is (cf. eq. (2.3))

$$Q_{penguin}^{PQS}=\frac{3}{N}(\overline{s}d)_L\left(\overline{u}_vu_v+\overline{d}_vd_v+\overline{s}_vs_v+\underset{i}{}\overline{q}_{si}q_{si}+\overline{u}_gu_g+\overline{d}_gd_g+\overline{s}_gs_g\right)_R,$$
(4.3)
where the subscripts $`v`$, $`s`$ and $`g`$ denote valence, sea and ghost quarks, respectively. Strategy 2 now corresponds to dropping all diagrams in which the right-handed valence and ghost quarks in the second factor of eq. (4.3) are contracted . If the number of sea quarks $`N=3`$ (but with the sea- and valence-quark masses not necessarily equal), the singlet LECs $`\alpha _{1,2}^{(8,1)}`$ are those of the real world , and therefore strategy 2 is the only correct one in this case.

For any other case, fully quenched or partially quenched with $`N3`$, there is a priori no preferred choice; the spread in results obtained by employing all three strategies should be taken as a (lower bound of the) systematic error due to quenching. The extent to which strategies 1 and 3 lead to numerically different results depends on the size of $`\alpha _q^{NS}`$ contributions (at a given scale $`\mathrm{\Lambda }`$). From eqs. (3.6,3) we find, taking physical values for all parameters, $`M_K=500`$ MeV, $`M_\pi =140`$ MeV, $`f=f_\pi =132`$ MeV, $`M_\rho =770`$ MeV and $`M_\eta =550`$ MeV, that

$`i[K^0\pi ^+\pi ^{}]_q`$ $`=`$ $`400.7(\alpha _{q1}^{(8,1)}\alpha _q^{(27,1)})+(28.27.2i)\alpha _q^{NS}(\mathrm{\Lambda }=1\mathrm{GeV}),`$ (4.4)
$`=`$ $`400.7(\alpha _{q1}^{(8,1)}\alpha _q^{(27,1)})+(32.27.2i)\alpha _q^{NS}(\mathrm{\Lambda }=M_\rho ),`$
$`=`$ $`400.7(\alpha _{q1}^{(8,1)}\alpha _q^{(27,1)})+(37.37.2i)\alpha _q^{NS}(\mathrm{\Lambda }=M_\eta ),`$
where we added in the tree-level ChPT contribution from the $`SU(3)_L`$ 27-plet operator . If $`\alpha _q^{NS}`$ is of the same order as $`\alpha _{q1}^{(8,1)}\alpha _q^{(27,1)}`$, the contribution of the terms proportional to $`\alpha _q^{NS}`$ is indeed small. The smallness of the coefficient of $`\alpha _q^{NS}`$ is due to a $`1/(4\pi )^2`$ suppression factor coming from the loop integral, and one might argue that $`\alpha _q^{NS}/(4\pi )^2`$ is the “natural” parameter to compare with $`\alpha _{q1}^{(8,1)}\alpha _q^{(27,1)}`$, in which case the contribution would not be small. Notice also that a small spurious imaginary part is generated by the non-singlet operator via the ghost-pion one-loop rescattering diagram. It is clear that the value of $`\alpha _q^{NS}`$ will have to be determined from a lattice computation. While this can be done by including the $`\alpha _q^{NS}`$ terms, of e.g. eq. (4), in a fit to lattice data, there exists a much simpler and more reliable way of estimating the size of $`\alpha _q^{NS}`$, as will be explained in the next section.

Under the assumption that $`\alpha _q^{NS}`$ can be neglected without introducing a large uncertainty into the final estimate of strong penguin $`K\pi \pi `$ matrix elements, the question remains whether (to leading order in ChPT) $`\alpha _{q1}^{(8,1)}`$ or $`\alpha _{q1}^{(8,1)}\frac{1}{(4\pi )^2}(\beta _{q1}^{NS}+\frac{1}{2}\beta _{q2}^{NS})`$ would be a better estimate of $`\alpha _1^{(8,1)}`$. The issue was investigated in ref. , where it was found that the difference between the two choices is numerically significant. At the physical kaon mass, the numerical value of the $`B`$ parameter corresponding to $`Q_6`$ turns out to be approximately twice as large when the contribution of the non-singlet operator $`Q_{penguin}^{QNS}`$ is omitted altogether. Translated into estimates for the leading-order LECs, this implies that $`\alpha _{q1}^{(8,1)}`$ is approximately twice as large as $`\alpha _{q1}^{(8,1)}\frac{1}{(4\pi )^2}(\beta _{q1}^{NS}+\frac{1}{2}\beta _{q2}^{NS})`$. This may lead to substantial modifications in quenched estimates of $`\epsilon ^{}/\epsilon `$, as discussed in ref. .

We emphasize that the whole discussion here is based on leading-order ChPT, and that NLO contributions may still lead to a substantial correction. However, it is reasonable to believe that NLO effects will not invalidate the basic content of our observations.

Finally, we give a few more numerical examples of the partially quenched case with $`N=2`$, always keeping the valence quark masses at their physical values (in the isospin limit), and choosing the two sea quarks to be degenerate in mass. Taking $`m_{sea}=m_u=m_d`$, we find

$`i[K^0\pi ^+\pi ^{}]_{N=2}`$ $`=`$ $`400.7(\alpha _1^{(8,1)}\alpha ^{(27,1)})+12.7\alpha ^{(8,8)}(\mathrm{\Lambda }=1\mathrm{GeV}),`$ (4.5)
$`=`$ $`400.7(\alpha _1^{(8,1)}\alpha ^{(27,1)})+14.0\alpha ^{(8,8)}(\mathrm{\Lambda }=M_\rho ),`$
$`=`$ $`400.7(\alpha _1^{(8,1)}\alpha ^{(27,1)})+15.7\alpha ^{(8,8)}(\mathrm{\Lambda }=M_\eta ),`$
whereas taking $`m_{sea}=m_s`$ we obtain

$`i[K^0\pi ^+\pi ^{}]_{N=2}`$ $`=`$ $`400.7(\alpha _1^{(8,1)}\alpha ^{(27,1)})+(2.87.2i)\alpha ^{(8,8)}(\mathrm{\Lambda }=1\mathrm{GeV}),`$ (4.6)
$`=`$ $`400.7(\alpha _1^{(8,1)}\alpha ^{(27,1)})+(4.17.2i)\alpha ^{(8,8)}(\mathrm{\Lambda }=M_\rho ),`$
$`=`$ $`400.7(\alpha _1^{(8,1)}\alpha ^{(27,1)})+(5.87.2i)\alpha ^{(8,8)}(\mathrm{\Lambda }=M_\eta ).`$
Recall that for the $`N=2`$ theory values of the LECs do not have to equal those of the $`N=3`$ theory. However, the partially quenched theory with two light sea quarks is closer to the real-world theory than the quenched ($`N=0`$) theory. This is reflected by the fact that the coefficients of $`\alpha ^{(8,8)}`$ are small compared to those of $`\alpha _q^{NS}`$ in eq. (4.4). Notice in addition that in the $`N=2`$ case, with $`m_{sea}=m_u=m_d`$, the small spurious imaginary part vanishes, since it comes entirely from the pion-rescattering loop diagram where the sea-quark contribution is now fully cancelled by the corresponding ghost-quark contribution.

In the case of $`N=3`$ sea quarks with masses equal to the three valence quarks, ghost- and sea-quark contributions in eqs. (3.2,3.4) cancel,555In eq. (3.4) the cancellation already occurs by just setting $`N=3`$, because in this tree-level expression the sea- and ghost-quark masses do not appear. ¶as they should, because this choice of parameters corresponds precisely to unquenched QCD.

## 5 How to determine $`\alpha _q^{NS}`$ on the lattice

In principle, it is possible to determine $`\alpha _q^{NS}`$ from matrix elements with only physical (valence) particles as external states. For instance, given good enough statistics and a wide enough range of quark masses, it can be determined from a fit to eq. (4). However, as also pointed out in ref. , the logarithmic terms in eq. (4) can look very linear in the typical range of quark masses used in lattice computations, making it hard to disentangle $`\alpha _q^{NS}`$ from $`\alpha _{q2}^{(8,1)}+\frac{1}{(4\pi )^2}\beta _{q3}^{(8,1)}`$. It would therefore be preferable to determine $`\alpha _q^{NS}`$ from a matrix element to which it contributes at order $`p^0`$, because no other operators can “contaminate” the result at that order.

It is very simple to do so, by considering matrix elements with ghost quarks on the external lines instead of valence quarks. Since this corresponds to a flavor rotation on the external lines, one needs to rotate the operator $`Q_{penguin}^{QNS}`$ accordingly. A key point is that, while of course ghost quarks are not explicitly present in a quenched computation, their propagators are identical to those of the valence quarks, which are available in the actual computation.

So, in order to determine $`\alpha _q^{NS}`$, we propose to consider the following matrix element. First, we rotate $`Q_{penguin}^{QNS}`$ by an $`SU(3|3)_L`$ rotation into

$$\stackrel{~}{Q}_{penguin}^{QNS}=(\overline{s}\gamma _\mu P_L\stackrel{~}{d})\text{str}(\widehat{N}\psi \overline{\psi }\gamma _\mu P_R).$$
(5.1)
This operator is in the same irrep of the group $`SU(3|3)_L\times SU(3|3)_R`$, is thus parametrized by the same LECs as $`Q_{penguin}^{QNS}`$, and in particular, to leading order, by $`\alpha _q^{NS}`$. We then consider the matrix element of this operator between a fermionic kaon $`\stackrel{~}{K}\overline{\stackrel{~}{d}}\gamma _5s`$ and the vacuum. To leading order,

$$0|\stackrel{~}{Q}_{penguin}^{QNS}|\stackrel{~}{K}=2if\alpha _q^{NS}+O(p^2),$$
(5.2)
thus isolating $`\alpha _q^{NS}`$. Carrying out all quark Wick contractions, one finds that

$`0|\stackrel{~}{Q}_{penguin}^{QNS}(y)\overline{\stackrel{~}{d}}(x)\gamma _5s(x)|0`$ $`=`$
$`{\displaystyle \frac{1}{2}}\{\text{tr}\left[\gamma _5s(x)\overline{s}(y)\gamma _\mu P_L\stackrel{~}{d}(y)\overline{\stackrel{~}{d}}(x)\right]\text{tr}[\gamma _\mu P_R(u(y)\overline{u}(y)+d(y)\overline{d}(y)+s(y)\overline{s}(y)`$
$`+\stackrel{~}{u}(y)\overline{\stackrel{~}{u}}(y)+\stackrel{~}{d}(y)\overline{\stackrel{~}{d}}(y)+\stackrel{~}{s}(y)\overline{\stackrel{~}{s}}(y))]`$
$`\text{tr}\left[\gamma _5s(x)\overline{s}(y)\gamma _\mu P_Rs(y)\overline{s}(y)\gamma _\mu P_L\stackrel{~}{d}(y)\overline{\stackrel{~}{d}}(x)\right]`$
$`+\text{tr}\left[\gamma _5s(x)\overline{s}(y)\gamma _\mu P_L\stackrel{~}{d}(y)\overline{\stackrel{~}{d}}(y)\gamma _\mu P_R\stackrel{~}{d}(y)\overline{\stackrel{~}{d}}(x)\right]\},`$
where the traces are over spin and color indices only. A key observation is now that ghost propagators and valence propagators are equal flavor by flavor, $`\stackrel{~}{d}(y)\overline{\stackrel{~}{d}}(x)=d(y)\overline{d}(x)`$ , etc.. Using this property, eq. (5) simplifies to

$`0|\stackrel{~}{Q}_{penguin}^{QNS}(y)\overline{\stackrel{~}{d}}(x)\gamma _5s(x)|0`$ $`=`$
$`\text{tr}\left[\gamma _5s(x)\overline{s}(y)\gamma _\mu P_Ld(y)\overline{d}(x)\right]\text{tr}\left[\gamma _\mu P_R(u(y)\overline{u}(y)+d(y)\overline{d}(y)+s(y)\overline{s}(y))\right]`$
$`+{\displaystyle \frac{1}{2}}\text{tr}\left[\gamma _5s(x)\overline{s}(y)\gamma _\mu P_Rs(y)\overline{s}(y)\gamma _\mu P_Ld(y)\overline{d}(x)\right]`$
$`{\displaystyle \frac{1}{2}}\text{tr}\left[\gamma _5s(x)\overline{s}(y)\gamma _\mu P_Ld(y)\overline{d}(y)\gamma _\mu P_Rd(y)\overline{d}(x)\right].`$
We conclude that it is possible to estimate $`\alpha _q^{NS}`$ as a leading-order effect using only combinations of contractions of valence-quark propagators. For the $`K0`$ matrix element of $`Q_{penguin}^{QNS}`$, cf. eqs. (2.9,4), the contractions in terms of valence quarks are of the same form, but the first two terms have the opposite sign, while the last term has the same sign as in eq. (5). Since the $`K0`$ matrix element is of order $`p^2`$, we may combine the two results to obtain

$`\text{tr}\left[\gamma _5s(x)\overline{s}(y)\gamma _\mu P_Ld(y)\overline{d}(y)\gamma _\mu P_Rd(y)\overline{d}(x)\right]_{amputated}`$
$`=\sqrt{Z}\left(0|\stackrel{~}{Q}_{penguin}^{QNS}|\stackrel{~}{K}+0|Q_{penguin}^{QNS}|K\right)`$
$`=\sqrt{Z}\left(2if\alpha _q^{NS}+O(p^2)\right),`$
making it even easier to determine $`\alpha _q^{NS}`$. The wave-function renormalization $`Z`$ is defined from $`\overline{d}\gamma _5s=\sqrt{Z}K`$.

The analysis for a similar determination of $`\alpha ^{(8,8)}`$ in the partially quenched theory is analogous. There, of course, the observation is not new, since $`\alpha ^{(8,8)}`$ is also the leading LEC for the electro-magnetic penguin (which in the partially quenched theory with $`N1`$ is in the same irrep as $`Q_{penguin}^{PQA}`$ of eq. (2.6) ). The main differences between a determination of $`\alpha ^{(8,8)}`$ and $`\alpha _q^{NS}`$ are that, first, $`\alpha _q^{NS}`$ is not related to the electro-magnetic penguin in the quenched case , and second, that in order to determine it using leading-order (in this case $`O(p^0)`$) ChPT, one is forced to consider ghost quarks, as we did above.

## 6 Conclusion

In this paper, we continued our investigation of the ambiguities afflicting strong penguin contributions to $`K\pi \pi `$ weak matrix elements due to the use of the quenched approximation.

The fact that the way of embedding penguin operators of the effective weak hamiltonian in the quenched theory is not unique tells us that, in the enlarged context of electro-weak interactions, the usual definition of the quenched theory is not complete. If only strong interactions are considered, it is sufficient to define quenched QCD as the modified version of QCD in which the quark determinant is set equal to a constant. A field-theoretic definition can be given through the introduction of ghost quarks into the path-integral , giving access to a complete picture of the symmetries of the quenched theory . As soon as one considers operators external to QCD (i.e. the addition of electro-weak interactions), one has to answer the question how these operators should be incorporated into the quenched theory. Usually, this is straightforward. One classifies the operator by its flavor quantum numbers, in other words, one determines the irrep of $`SU(3)_L\times SU(3)_R`$ under which this operator transforms. If there exists a larger irrep of the quenched symmetry group which reduces to the unquenched irrep, the corresponding component of the quenched irrep can be taken as the quenched definition of the operator. However, in the case of strong penguins the operator, while irreducible in the unquenched theory, is a linear combination of components of more than one irrep of the quenched symmetry group. Therefore each LEC of the unquenched theory corresponds to a set of LECs in the quenched theory. The ambiguity arises, because there is a priori no criterium for which linear combination of quenched LECs (if any) would yield the best estimate of the unquenched LEC. In the case of $`LR`$ penguins considered here and in ref. , this phenomenon produces an effect already at leading order in ChPT.

We remark that even in the simplest case when there exists a one-to-one correspondence between unquenched and quenched irreps, there is still the freedom to choose any component of the quenched irrep, and this flexibility can be used to extract LECs in the most convenient way . However, in this case there is no ambiguity in the relation between unquenched and quenched LECs (even though their values may differ). This is in principle not different from the situation within the unquenched theory, where in general any component of an irrep can be used to extract the corresponding LEC. A classic example for weak matrix elements is the relation between $`B_K`$ and the $`K^+\pi ^+\pi ^0`$ decay rate . (At non-leading order, it may not be possible to determine all LECs describing an operator in ChPT from one process, of course.)

The ambiguity affecting penguin operators is fundamental, since there exists no solid theoretical argument that can be used to decide the issue. Therefore, we argue that one should compare all choices that can be reasonably made, and take the resulting spread of estimated values as a lower bound on the systematic error due to quenching. It appears that in the case of $`\epsilon ^{}/\epsilon `$ this systematic error is rather large . A leading-order analysis of currently available lattice data  seems to indicate that quenched lattice computations cannot even confirm that this parameter is non-vanishing in the Standard Model. It could also be that the large numerical difference found between $`\alpha _{q1}^{(8,1)}`$ and $`\alpha _{q1}^{(8,1)}\frac{1}{(4\pi )^2}(\beta _{q1}^{NS}+\frac{1}{2}\beta _{q2}^{NS})`$ would be explained by the fact that higher orders in ChPT have not been taken into account, but we consider this to be unlikely. While it is clear that higher orders are numerically important, there appears to be no reason to assume that $`\beta _{q1}^{NS}+\frac{1}{2}\beta _{q2}^{NS}`$ is small. It could also be that $`\alpha _q^{NS}`$, which appears in $`K^00`$, and needs to be subtracted to obtain $`\alpha _{q1}^{(8,1)}\frac{1}{(4\pi )^2}(\beta _{q1}^{NS}+\frac{1}{2}\beta _{q2}^{NS})`$ from $`K^+\pi ^+`$, is not small. This would affect the determination of $`\alpha _{q2}^{(8,1)}+\frac{1}{(4\pi )^2}\beta _{q3}^{NS}`$, and hence the size of the subtraction. It is therefore important to obtain a reliable estimate of $`\alpha _q^{NS}`$. We suggested a simple method for extracting its value.

The above argument does not imply that lattice computations of $`\epsilon ^{}/\epsilon `$ are doomed to fail. On the contrary, quenched estimates of $`\epsilon ^{}/\epsilon `$ with a particular choice for the strong penguins demonstrate that this computation is feasible, thanks to major advances in both theory and computational power. However, what will be needed in order to eliminate systematic errors due to quenching is a partially quenched study with $`N=3`$ light sea quarks. This is the only approximation to unquenched QCD which is reliable in that it can be extrapolated systematically to the real world . Currently existing quenched results give us invaluable information on what is needed to promote them to the required $`N=3`$ world. For partially quenched QCD with $`N3`$, the situation is essentially the same as for quenched QCD, modulo differences in detail.

Summarizing, we presented the quenched and partially quenched results for the non-singlet contribution to the $`K^0\pi ^+\pi ^{}`$ matrix element. This made it possible to discuss in detail various strategies one might follow to use quenched computations in order to estimate the real-world value of this amplitude. Since $`\alpha _q^{NS}`$ contributes to this matrix element, but not to the $`K^+\pi ^+`$ transition amplitude used in ref. , this introduces an additional ambiguity already at leading chiral order. The importance of this ambiguity depends on the size of $`\alpha _q^{NS}`$ and we have proposed a simple recipe for its determination. Our expressions for $`K^0\pi ^+\pi ^{}`$ with the most general possible kinematics and the inclusion of the non-singlet contributions, are appropriate for the analysis of direct quenched computations of this matrix element at leading order in ChPT. Beyond leading order, new problems arise , which may invalidate current methods for the direct determination of $`K\pi \pi `$ amplitudes with $`\mathrm{\Delta }I=1/2`$ in quenched and partially quenched QCD.

## Appendix

In this appendix we collect explicit expressions for the basic loop integrals appearing in eq. (3.2) etc. Using dimensional regularization, we have

$$L(M^2)=\frac{d^D\mathrm{}}{(2\pi )^D}\frac{1}{\mathrm{}^2+M^2}=\frac{M^2}{16\pi ^2}\left(\left[\frac{2}{ϵ}+\gamma \mathrm{log}4\pi \right]+\mathrm{log}\frac{M^2}{\mathrm{\Lambda }^2}1\right),$$
(A.1)
where $`\mathrm{\Lambda }`$ is the running scale, $`ϵ=4D`$, and

$`\mathrm{Re}I(M_1^2,M_2^2,p^2)`$ $`=`$ $`\mathrm{Re}{\displaystyle \frac{d^D\mathrm{}}{(2\pi )^D}\frac{1}{(\mathrm{}^2+M_1^2)((\mathrm{}p)^2+M_2^2)}}`$
$`={\displaystyle \frac{1}{16\pi ^2}}\{[{\displaystyle \frac{2}{ϵ}}+\gamma \mathrm{log}4\pi ]`$
$`1+\mathrm{log}{\displaystyle \frac{M_2^2}{\mathrm{\Lambda }^2}}+{\displaystyle \frac{1}{2}}(1{\displaystyle \frac{M_1^2}{p^2}}+{\displaystyle \frac{M_2^2}{p^2}})\mathrm{log}{\displaystyle \frac{M_1^2}{M_2^2}}+{\displaystyle \frac{1}{2}}F(M_1^2,M_2^2,p^2)\},`$
in which

$`F(M_1^2,M_2^2,p^2)`$ $`=`$ $`\sqrt{\lambda (1,{\displaystyle \frac{M_1^2}{p^2}},{\displaystyle \frac{M_2^2}{p^2}})}\mathrm{log}{\displaystyle \frac{p^2+M_1^2+M_2^2+p^2\sqrt{\lambda (1,M_1^2/p^2,M_2^2/p^2)}}{p^2+M_1^2+M_2^2p^2\sqrt{\lambda (1,M_1^2/p^2,M_2^2/p^2)}}},`$
$`\lambda (x,y,z)`$ $`=`$ $`(xy+z)^2+4xy.`$ (A.3)
$`\overline{MS}`$ expressions are obtained by dropping the contact terms in square brackets.

For $`p^2>0`$, the argument of the logarithm in eq. (A.3) is positive, and $`I`$ is real. For $`p^20`$, $`F`$ is obtained by analytic continuation. $`\lambda (1,M_1^2/p^2,M_2^2/p^2)`$ turns negative for $`(M_1+M_2)^2<p^2<(M_1M_2)^2`$, and we find that $`I`$ is still real with $`F`$ now given by

$$F(M_1^2,M_2^2,p^2)=2\sqrt{\lambda (1,\frac{M_1^2}{p^2},\frac{M_2^2}{p^2})}\mathrm{arctan}\frac{p^2\sqrt{\lambda (1,M_1^2/p^2,M_2^2/p^2)}}{p^2+M_1^2+M_2^2}.$$
(A.4)
At $`p^2=(M_1^2+M_2^2)`$ the argument of the arctangent has a singularity, across which the branch of the arctangent has to be chosen continuously:

$`\mathrm{arctan}{\displaystyle \frac{p^2\sqrt{\lambda (1,M_1^2/p^2,M_2^2/p^2)}}{p^2+M_1^2+M_2^2}}`$ $`=`$ $`\mathrm{Arctan}{\displaystyle \frac{p^2\sqrt{\lambda (1,M_1^2/p^2,M_2^2/p^2)}}{p^2+M_1^2+M_2^2}}+\pi ,`$ (A.5)
$`p^2<(M_1^2+M_2^2),`$
where Arctan denotes the principal value of the arctangent. Again continuing analytically across $`p^2=(M_1+M_2)^2`$, $`F`$ is again given by eq. (A.3), but $`I(M_1^2,M_2^2,p^2)`$ picks up an imaginary part:

$$\mathrm{Im}I(M_1^2,M_2^2,p^2)=\frac{1}{16\pi ^2}\pi \sqrt{\lambda (1,\frac{M_1^2}{p^2},\frac{M_2^2}{p^2})}.$$
(A.6)
## Acknowledgements

We would like to thank Norman Christ, Bob Mawhinney, Martin Savage and Steve Sharpe for discussions. MG thanks the Institute for Nuclear Theory at the University of Washington for hospitality, and EP thanks the Department of Physics and Astronomy at San Francisco State University for hospitality. MG is supported in part by the US Department of Energy, and EP by the Italian MURST under the program Fenomenologia delle Interazioni Fondamentali.