text:
# Symmetries and Triplet Dispersion in a Modified Shastry-Sutherland Model for SrCu2(BO3)2

## I Introduction

The study of quantum spin systems exhibiting a finite spin gap has advanced significantly through the recent synthesis of novel magnetic materials. While quasi-one dimensional systems have been studied for a long time now, new interest arose from the discovery of new quasi-two dimensional materials like CaV4O9  and (VO2)P2O7 . Particularly interesting is SrCu2(BO3)2  since it is an experimental realization of the Shastry-Sutherland model . In Ref.  we presented an extended version of that model (Fig. 1) and deduced its $`T=0`$ phase-diagram in the model parameter space.

The Hamiltonian is given by

$$H=J_1\underset{H_0}{\underset{}{\underset{}{}\stackrel{}{S}_i\stackrel{}{S}_j}}+J_2\underset{H_1}{\underset{}{\underset{}{}\stackrel{}{S}_i\stackrel{}{S}_k}}+J_3\underset{H_2}{\underset{}{\underset{}{}\stackrel{}{S}_i\stackrel{}{S}_l}},$$
(1)
where the bonds corresponding to interactions $`J_1`$, $`J_2`$ and $`J_3`$ are shown in Fig. 1. For $`J_3=0`$ the model reduces to the original Shastry-Sutherland model.
It can be easily verified that the singlet-dimer state (singlets on all strong bonds $`J_1`$, henceforth called dimers) is an exact eigen state of our model: a single spin interacts via $`J_2`$ and $`J_3`$ with $`S=0`$ objects only. Thus the corresponding terms in Hamiltonian (1) do not contribute. The remaining expression simply gives $`(3/8)J_1`$ per spin. This article focuses on the region, where the singlet-dimer state is the ground state, called the dimer phase. In this phase a single excitation (magnon) is introduced by breaking up one singlet and substituting one triplet instead (the triplet’s z-component is irrelevant). By hopping from dimer to dimer this triplet acquires a dispersion, which we intend to calculate. It suffices to do the calculations on an effective square lattice $`\mathrm{\Gamma }_{\mathrm{eff}}`$, where one site represents one dimer. Note that $`\mathrm{\Gamma }_{\mathrm{eff}}`$ has a sublattice structure A/B, where A (B) corresponds to vertical (horizontal) dimers.
We derive a power series expansion about the limit of strong dimerization ($`J_2,J_30`$) for the one-magnon dispersion. By fitting the expansion to INS-Data obtained by Kageyama et al.  we deduce model parameters which show good agreement with parameters previously determined.
The next section gives a short introduction to the method used. Before we calculate and discuss the one-magnon dispersion (Sec. IV) we use some of the model’s symmetries to derive interesting and useful relations between various hopping amplitudes (Sec. III).

## II Method and qualitative pictures

The one-magnon dispersion is calculated perturbatively about the limit of isolated dimers using the flow equation method introduced previously . Given a perturbation problem that can be formulated in the standard way

$$=_0+x_\mathrm{S}$$
(2)
this method in its present formulation relies only on two further prerequisites:

1. The unperturbed Hamiltonian $`_0`$ must have an equidistant spectrum bounded from below. The difference between two successive levels is called an energy quantum.
2. There is a number $`N>0`$ such that the perturbing Hamiltonian $`H_\mathrm{S}`$ can be written as $`H_\mathrm{S}=_{n=N}^NT_n`$ where $`T_n`$ increments (or decrements, if $`n<0`$) the number of energy quanta by $`n`$.

The flow equation method maps the perturbed Hamiltonian $``$ by a continuous unitary transformation to an effective Hamiltonian $`_{\mathrm{eff}}`$, which conserves the number of energy quanta, i.e. $`[_{\mathrm{eff}},_0]=0`$. Thus the effective Hamiltonian is block diagonal and has the form

$$_{\mathrm{eff}}=_0+\underset{k=1}{\overset{\mathrm{}}{}}x^k\underset{|\underset{¯}{m}|=k,M(\underset{¯}{m})=0}{}C(\underset{¯}{m})T(\underset{¯}{m}),$$
(3)
where $`\underset{¯}{m}`$ is a vector of dimension $`k`$ of which the components are in $`\{\pm N,\pm (N1),\mathrm{}\pm 1,0\}`$; $`M(\underset{¯}{m})=0`$ signifies that the sum of the components vanishes which reflects the conservation of the number of energy quanta. The operator product $`T(\underset{¯}{m})`$ is defined by $`T_{\underset{¯}{m}}=T_{m_1}T_{m_2}\mathrm{}T_{m_k}`$, where $`k`$ is the order of the process. The coefficients $`C(\underset{¯}{m})`$ are generally valid fractions, which we computed up to order $`k=15`$ for $`N=1`$ and up to order $`k=10`$ for $`N=2`$ (cf. Ref. ).
We want to emphasize that the effective Hamiltonian $`_{\mathrm{eff}}`$ with known coefficients $`C(\underset{¯}{m})`$ can be used straightforwardly in all perturbative problems that meet conditions (A) and (B). The interested reader can find the $`C(\underset{¯}{m})`$ and additional information on our home-pages. For checking purposes we tested the $`C(\underset{¯}{m})`$ by applying the method to toy models. In the $`N=1`$ case, for instance, we considered the one dimensional harmonic oscillator perturbed by itself and an additional linear potential. Thus, with $`\mathrm{}\omega =1`$

$`_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}(P^2+X^2)=a^{}a+{\displaystyle \frac{1}{2}}`$ (4)
$`_\mathrm{S}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(P^2+X^2)+X=\underset{T_0}{\underset{}{a^{}a+{\displaystyle \frac{1}{2}}}}+\underset{T_1}{\underset{}{{\displaystyle \frac{1}{\sqrt{2}}}a^{}}}+\underset{T_1}{\underset{}{{\displaystyle \frac{1}{\sqrt{2}}}a}}.`$ (5)
In this case Eq. (2) can be solved exactly to give

$$E_n=(1+x)\left(n+\frac{1}{2}\right)\frac{1}{2}\frac{x^2}{1+x}.$$
(6)
We expand this last equation about $`x=0`$ and obtain

$$E_n=n+\frac{1}{2}+\frac{1}{2}(2n+1)x\frac{1}{2}\underset{i=2}{\overset{\mathrm{}}{}}(1)^ix^i.$$
(7)
By inserting the $`T_i`$ defined in Eq. (5) in the effective Hamiltonian (3) and calculating $`n|_{\mathrm{eff}}|n`$ with $`n`$ we retain Eq. (7) exactly up to 15th order.

We now show that Hamiltonian (1) meets conditions (A) and (B) for $`N=1`$. To this end we rewrite Eq. (1)

$$\frac{H}{J_1}=H_0+xH_\mathrm{S},\text{with}$$
(8)
$$H_\mathrm{S}=H_1+\frac{y}{x}H_2,x=\frac{J_2}{J_1}\text{ and }y=\frac{J_3}{J_1}.$$
(9)
In the limit of isolated dimers ($`x=0`$, with $`y/x`$ finite) $`H`$ is bounded from below and has an equidistant energy spectrum since up to a trivial constant $`H_0`$ simply counts the number of excited dimers, i.e. energy quanta.
To decompose $`H_\mathrm{S}`$ we follow the same proceeding as in Ref.  and state the result

$$H_\mathrm{S}=T_1+T_0+T_1$$
(10)
with

$`T_{\pm 1}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(1{\displaystyle \frac{y}{x}}\right){\displaystyle \underset{\nu }{}}𝒯_{\pm 1}(\nu )`$ (11)
$`T_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(1+{\displaystyle \frac{y}{x}}\right){\displaystyle \underset{\nu }{}}𝒯_0(\nu ),`$ (12)
where $`\nu `$ denotes the pairs of adjacent dimers. For some fixed $`\nu `$ we encounter the pair depicted in Fig. 2, the state of which is determined by $`|x_1,x_2`$, where $`x_1,x_2\{s,t^1,t^0,t^1\}`$ are singlets or one of the triplets occupying the vertical (horizontal) dimer, respectively. The superscript $`n\{0,\pm 1\}`$ in $`t^n`$ stands for the $`S^z`$ component. For the pair $`|x_1,x_2`$ we give the action of the local operators $`𝒯_i`$ in Tab. I.

The remaining matrix elements can be constructed by using $`𝒯_n^{}=𝒯_n`$. Note that we need to fix the orientation for singlets, say spin up with positive sign always at the right (upper) site of the dimers. Hence $`𝒯_1`$ and $`𝒯_1`$ acquire a global minus for oppositely oriented dimer pairs (reflection of the dimer pair in Fig. 2 about the vertical dimer).

Instead of Hamiltonian (1) we can now use the effective Hamiltonian $`_{\mathrm{eff}}`$ (Eq. (3)) with the $`T_i`$ defined in Eqs. (11,12). The effective Hamiltonian (3) simplifies the computations considerably. Since $`[_{\mathrm{eff}},_0]=0`$, $`_{\mathrm{eff}}`$ is block diagonal allowing triplet conserving processes only. Thus the effective Hamiltonian acts in a much smaller Hilbert space than the original problem, which is a great advantage in the numerical implementation.
In addition to this simplification the explicit form of $`_{\mathrm{eff}}`$ provides a simple and comprehensive picture of the physics involved. Imagine we were to put a triplet on one of the dimers in the lattice. In the real substance this local excitation would polarize its environment due to the exchange couplings and must be viewed as a dressed quasi particle surrounded by a cloud of virtual excitations fluctuating in space and time. The effective Hamiltonian includes these fluctuations as virtual processes $`T(\underset{¯}{m})`$, each weighted by the factor $`C(\underset{¯}{m})`$. Each process ends with a state having the same number of triplets as the initial state. Other quantum numbers such as the total spin are also conserved. As the order $`k`$ is increased more and longer processes are allowed for and the accuracy of the results will be enhanced. Inspecting the weight factors $`C(\underset{¯}{m})`$ shows that longer processes have less influence.

Let us follow one of the possible virtual processes and understand why the observed triplet dispersion of SrCu2(BO3)2 is rather flat (cf. Fig. 8). Suppose we begin with one triplet in the lattice as depicted in the upper left corner of Fig. 3.

By applying $`T_1`$ once we can create another triplet only on one of the two $`horizontally`$ adjacent dimers as is clear from Tab. I and Fig. 2. From there we might create another one and so on till we can close a circle (bottom left state in Fig. 3). We now start destructing triplets by $`T_1`$ processes and end up with the shifted triplet. The amplitude for this hopping is $`x^6`$. It is the largest amplitude one can find (see also Ref. ). Therefore the triplets are rather localized leading to a flat dispersion.

Note that the $`T_1`$ ($`T_1`$) processes are proportional to $`(xy)`$ (see Eq. (11)). So, the leading triplet motion is $`𝒪((xy)^6)`$ in leading order. Even local processes (without hopping) include at least two $`T_1`$ ($`T_1`$) processes. Hence all perturbative amplitudes are at least of order $`(xy)^2`$. We will make use of this fact later on.

Examining the two-triplet sector  we showed that $`correlated`$ hopping processes occur in second order already. The actual dispersion, however, sets in only in third order. This much lower order ($`x^3`$ instead of $`x^6`$) explains the much stronger two-magnon dispersion .

To quantify the picture constructed let $`|𝐫=|r_1,r_2`$ denote the state of the system with one triplet at $`𝐫\mathrm{\Gamma }_{\mathrm{eff}}`$ and singlets on all other sites. The amplitude $`t_{𝐫^{}𝐫}^{o(𝐫)}`$ for a triplet-hopping from site $`𝐫`$ to site $`𝐫^{}`$ is given by

$$t_{𝐫^{}𝐫}^{o(𝐫)}=𝐫^{}|_{\mathrm{eff}}|𝐫,$$
(13)
where the upper index $`o(𝐫)\{v,h\}`$ allows to distinguish whether the hopping started on a vertically oriented ($`v`$) or a horizontally oriented dimer ($`h`$). Furthermore we choose to split the hopping amplitudes into an average part $`\overline{t}_𝐬`$ and an alternating part $`dt_𝐬`$ ($`𝐬=𝐫^{}𝐫`$)

$$t_𝐬^{o(𝐫)}=\overline{t}_𝐬+e^{i\mathrm{𝐐𝐫}}dt_𝐬,$$
(14)
with $`𝐐=(\pi ,\pi )`$.
The right hand side of Eq. (13) can be easily implemented on a computer. For details see Ref. . We want to point out, however, that all computations are done symbolically, i.e. we obtain results as functions (polynomials) of all model parameters.

## III Symmetries

Before we calculate the one-triplet dispersion quantitatively in the next section it is worthwhile to look at the symmetries the model in Fig. 1 displays. The two-dimensional space group of the model can be identified to be p4mm with the underlying point group 4mm as can be verified in Fig. 4.

We will not treat all symmetry aspects but concentrate on those that will be of use later on. We choose the coordinate system parallel to the dimers such that one dimer (horizontal or vertical) lies in the origin and introduce two diagonals u and v crossing the origin with slope -1 and 1 , respectively. The distance between the centers of two adjacent dimers is set to unity.

Several relations between different hopping amplitudes can be deduced. To this end we define six symmetry operations which map the lattice onto itself. Note that the fixed singlet orientation can lead to negative phase factors.
$`m_{x/y}`$: Reflection about x/y-axis
$`m_{x/y}|r_1,r_2=(1)^{r_1+r_2}|\pm r_1,r_2`$.
$`I`$: Inversion about the origin
$`I|r_1,r_2=|r_1,r_2`$.
$`\sigma _{v/u}`$: Reflection about the diagonals v/u plus translation by $`(0,1)`$
$`\sigma _{v/u}|r_1,r_2=|\pm r_2,\pm r_11`$.
$`R`$: Rotation of $`\pi /2`$ about the origin plus translation by $`(0,1)`$
$`R|r_1,r_2=(1)^{r_1+r_2}|r_2,r_11`$.

Since these operations leave $``$ unchanged they all commute with $`_{\mathrm{eff}}`$. Applying $`m_x`$ to Eq. (13) we find

$`t_{𝐫^{}𝐫}^{o(𝐫)}`$ $`=`$ $`(1)^{r_1^{}+r_2^{}}r_1^{},r_2^{}|m_x_{\mathrm{eff}}|r_1,r_2`$ (15)
$`=`$ $`(1)^{r_1+r_2+r_1^{}+r_2^{}}r_1^{},r_2^{}|_{\mathrm{eff}}|r_1,r_2`$ (16)
$`=`$ $`(1)^{r_1+r_2+r_1^{}+r_2^{}}t_{(r_1^{},r_2^{})(r_1,r_2)}^{o(𝐫)},`$ (17)
or simply

$$t_𝐬^{o(𝐫)}=(1)^{s_1+s_2}t_{s_1,s_2}^{o(𝐫)}.$$
(18)
Analogously using $`m_y`$, $`I`$, $`\sigma _u`$, $`\sigma _v`$ and $`R`$ we find

$`t_𝐬^{o(𝐫)}`$ $`=`$ $`(1)^{s_1+s_2}t_{s_1,s_2}^{o(𝐫)}`$ (19)
$`=`$ $`t_𝐬^{o(𝐫)}`$ (20)
$`=`$ $`t_{s_2,s_1}^{o(𝐫(0,1))}`$ (21)
$`=`$ $`t_{s_2,s_1}^{o(𝐫(0,1))}`$ (22)
$`=`$ $`(1)^{s_1+s_2}t_{s_2,s_1}^{o(𝐫(0,1))}`$ (23)
respectively. In particular Eqs. (19) and (20) yield together

$`t_{(s_1,0)}^{o(𝐫)}`$ $`=`$ $`(1)^{s1}t_{(s1,0)}^{o(𝐫)}=t_{(s1,0)}^{o(𝐫)}`$ (24)
$`t_{(s_1,0)}^{o(𝐫)}`$ $`=`$ $`0\text{, if }s_1\text{ odd, analogously}`$ (25)
$`t_{(0,s_2)}^{o(𝐫)}`$ $`=`$ $`0\text{, if }s_2\text{ odd,}`$ (26)
describing the interesting fact, that hopping along the axis has non zero amplitude only if this hopping moves the triplet an even number of sites ($`\mathrm{\Gamma }_{\mathrm{eff}}`$).

## IV Dispersion

Since $`_{\mathrm{eff}}`$ conserves the number of triplets the one-triplet dispersion is particularly easy to calculate. Starting with one triplet this excitation can only be shifted, i.e. the triplet hops on the effective lattice $`\mathrm{\Gamma }_{\mathrm{eff}}`$. Additional care has to be taken to account for the sublattice structure of $`\mathrm{\Gamma }_{\mathrm{eff}}`$. From Eq. (13) we get

$$_{\mathrm{eff}}|𝐫=\underset{𝐫^{}}{}t_𝐫^{}^{o(𝐫)}|𝐫+𝐫^{}.$$
(27)
We introduce Fourier transformed states

$$|\sigma ,𝐤=\frac{1}{\sqrt{L}}\underset{𝐫}{}|𝐫e^{i(𝐤+\sigma Q)𝐫}$$
(28)
with the number of dimers $`L`$, the new quantum number $`\sigma \{0,1\}`$ reflecting the sublattice structure and $`𝐤`$ a vector in the magnetic Brillouin zone (MBZ). Calculating the action of $`_{\mathrm{eff}}|\sigma ,𝐤`$ on these states yields

$`_{\mathrm{eff}}|\sigma ,𝐤`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{L}}}{\displaystyle \underset{𝐫,𝐫^{}}{}}(\overline{t}_𝐫^{}+e^{iQ𝐫}dt_𝐫^{})|𝐫+𝐫^{}e^{i(𝐤+\sigma Q)𝐫}`$ (29)
$`=`$ $`{\displaystyle \underset{𝐫^{}}{}}\overline{t}_𝐫^{}e^{i(𝐤+\sigma Q)𝐫^{}}\underset{|\sigma ,𝐤}{\underset{}{{\displaystyle \frac{1}{\sqrt{L}}}{\displaystyle \underset{𝐫}{}}|𝐫e^{i(𝐤+\sigma Q)𝐫}}}`$ (31)
$`+{\displaystyle \underset{𝐫^{}}{}}dt_𝐫^{}e^{i(𝐤+\overline{\sigma }Q)𝐫^{}}\underset{|\overline{\sigma },𝐤}{\underset{}{{\displaystyle \frac{1}{\sqrt{L}}}{\displaystyle \underset{𝐫}{}}|𝐫+𝐫^{}e^{i(𝐤+\overline{\sigma }Q)(𝐫+𝐫^{})}}},`$
with $`\overline{\sigma }=1\sigma `$. We used the definitions (13) and (14). Further, since Eq. (20) holds for even and odd $`r_1+r_2`$ one has $`\overline{t}_𝐬=\overline{t}_𝐬\text{ and }dt_𝐬=dt_𝐬`$. Thus we can simplify the sums over $`𝐫^{}`$ in Eq. (31) to give

$`_{\mathrm{eff}}|\sigma ,𝐤`$ $`=`$ $`\underset{a_\sigma }{\underset{}{\left[\overline{t}_\mathrm{𝟎}+2{\displaystyle \underset{𝐫>0}{}}\overline{t}_𝐫\mathrm{cos}((𝐤+\sigma Q)𝐫)\right]}}|\sigma ,𝐤`$ (32)
$`+`$ $`\underset{b}{\underset{}{\left[dt_\mathrm{𝟎}+2{\displaystyle \underset{𝐫>0}{}}dt_𝐫\mathrm{cos}((𝐤+\overline{\sigma }Q)𝐫)\right]}}|\overline{\sigma },𝐤,`$ (33)
with $`𝐫>0`$ if and only if ($`r_1>0`$ or $`r_1=0`$ but $`r_2>0`$). In Appendix A we show that $`dt_\mathrm{𝟎}=0`$ and $`dt_𝐫=0`$ for $`r_1+r_2`$ odd. Hence $`b`$ does not depend on $`\overline{\sigma }`$

$$b=2\underset{\genfrac{}{}{0pt}{}{𝐫>0}{r_1+r_2\text{ even}}}{}dt_𝐫\mathrm{cos}(\mathrm{𝐤𝐫}),$$
(34)
and $`_{\mathrm{eff}}`$ is symmetric in the new states. The remaining $`2\times 2`$ problem can be solved easily to give the dispersion

$$\omega _{1/2}(𝐤)=\underset{\omega _0(𝐤)}{\underset{}{\frac{a_0+a_1}{2}}}\pm \frac{1}{2}\sqrt{(a_0a_1)^2+4b^2}.$$
(35)
Thus the one-triplet dispersion splits into two branches. We want to point out, however, that at $`𝐤=\mathrm{𝟎}`$ and at the borders of the MBZ (i.e. $`|k_x+k_y|=\pi `$ or $`|k_yk_x|=\pi `$) the two branches fall onto each other leading to a two-fold degenerate dispersion. An analogous degeneracy is noticed in the two-triplet sector . In Appendix B we demonstrate that the degeneracy is due to the glide line symmetries $`R`$ and $`\sigma _{u/v}`$ and show that $`(a_0a_1)`$ and $`b`$ both vanish. Moreover, $`(a_0+a_1)`$ is a sum over $`𝐫`$ with $`r_1+r_2`$ even and thus contains “even” hopping amplitudes only. At $`𝐤=\mathrm{𝟎}`$ and at the border of the MBZ we therefore find that one-triplet hopping takes place on one species of the two sublattices A/B in $`\mathrm{\Gamma }_{\mathrm{eff}}`$ only. In other words, the triplets live either on the horizontal or on the vertical dimers.

We calculated the amplitudes $`\overline{t}_𝐫`$ and $`dt_𝐫`$ (and therefore the dispersion) as exact polynomials in $`x`$ and $`y`$ up to and including 15th order. On appearance of this article these polynomials will be published in electronic form on our home pages.
Expanding the square root in Eq. (35) about the limit of vanishing $`x`$ and $`y`$ produces terms $`x^\alpha y^\beta `$ with $`\alpha +\beta 10`$. Hence the energy splitting starts in 10th order and is negligible for all reasonable values of $`x`$ and $`y`$. The fact that the splitting starts four orders later than the dispersion may be understood by observing that ($`a_0a_1`$) is a sum over $`𝐫`$ with $`r_1+r_2`$ odd. From the discussion at the end of section III it is clear that $`t_{\pm 1,0}^{o(𝐫)}`$ and $`t_{0,\pm 1}^{o(𝐫)}`$ vanish so that the leading “odd” process is $`t_{\pm 2,\pm 1}^{o(r)}`$ or $`t_{\pm 1,\pm 2}^{o(r)}`$, which start in 10th order only. One may object that $`b`$ contains $`dt_{1,1}`$, which could be larger, but from relation (A5) in Appendix A follows that $`dt_{1,1}`$=0. The amplitudes $`dt_{\pm 2,0}`$ and $`dt_{0,\pm 2}`$ in $`b`$ also start in 10th order only. Hence the almost degeneracy in the one triplet sector can be understood on the basis of the symmetries of the lattice.

By substituting $`y=0`$ in $`\omega _0(𝐤)`$ we obtain the decimal numbers computed by Zheng et al. . Our series expansion for the dispersion (35) converges nicely but it can be improved further by the use of D-log Padé approximants .

Let us first consider the energy gap $`\mathrm{\Delta }:=\omega (0,0)`$ as function of $`x`$ and $`y`$. We fix the ratio of $`x`$ and $`y`$ and use D-log Padé approximants to extrapolate the remaining expression. For all $`x/y`$ ratios we tested, we find the approximants to be very stable, i.e. most of the possible approximants at a fixed ratio coincide very well. (Some are defective, i.e they show spurious singularities.) Fig. 5 shows D-log Padé extrapolations for $`y=x\mathrm{tan}(\pi /6)`$ and $`y=x\mathrm{tan}(\pi /8)`$. In Ref.  we used this technique to determine the line in the ($`x,y`$)-plane where the gap $`\mathrm{\Delta }`$ vanishes. The vanishing of $`\mathrm{\Delta }`$ indicates definitively the break down of the dimer phase. But it may happen that another excitation becomes soft before the elementary triplet vanishes (cf. Ref. ) or that a first order transition takes place . In the present work we choose to examine the elementary triplet only.

For $`y=x\mathrm{tan}(\pi /6)`$ the gap vanishes at $`x=0.558(1)`$ ($`y=0.323(1)`$). At this point we get the lowest curve in Fig. 6 where we choose to show the dispersion along a triangle in the MBZ.

In Sec. II we saw that all hopping amplitudes are proportional to $`(xy)^6`$. Since the dispersion is a sum over these amplitudes it is clear that $`(\omega (\pi ,0)\omega (0,0))(xy)^6`$. Indeed, for increasing $`y`$ we see that the dispersion decreases (Fig. 6). At $`x=y`$ we find the dispersion to be absolutely flat. This is a signiture of the fact that at $`x=y`$ the total spin on each $`J_1`$ bond is a conserved quantity . Thus there will be no triplet motion.

We turn to the behaviour of the dispersion $`\omega (𝐤)`$ for small $`|𝐤|`$ on the line in the ($`x,y`$)-plane where the gap vanishes. After fixing the ratio $`x/y`$ in Eq. (35) and applying the D-log Padé technique we end up with an expression

$$\omega (x;𝐤)=\mathrm{exp}\left[_0^x\underset{f(x;𝐤)}{\underset{}{A(𝐤)\frac{P(x;𝐤)}{Q(x;𝐤)}}}\right]dx,$$
(36)
where $`A`$ is a function of $`𝐤`$ and $`P`$ and $`Q`$ are polynomials in $`x`$ (where the leading coefficient is unity) of order $`M`$ and $`N`$, respectively (shorthand: \[M,N\]-approximant). At $`𝐤=0`$ the smallest positive zero of $`Q`$, say $`x_0`$, is the point where the exponent diverges logarithmically to $`\mathrm{}`$, hence $`\omega (x_0;0)=0`$. We set $`k_y=0`$ and will show that for small $`k_x`$

$$\omega (x_0;k_x)k_x^\alpha \text{ + higher orders}.$$
(37)
Differentiating $`\mathrm{ln}\omega `$ yields for $`k_x0`$

$$\alpha =\underset{k_x0}{lim}k_x_0^{x_0}_{k_x}f(x;k_x)dx.$$
(38)
We decompose $`P`$ and $`Q`$ in linear factors

$`P(x;k_x)`$ $`=`$ $`(xp_1(k_x))(xp_2(k_x))\mathrm{}(xp_M(k_x))`$ (39)
$`Q(x;k_x)`$ $`=`$ $`(xq_1(k_x))(xq_2(k_x))\mathrm{}(xq_N(k_x)),`$ (40)
such that $`q_1(k_x)|_{k_x=0}=x_0`$. In this way the factor $`(xq_1(k_x))`$ dominates the behaviour of $`\alpha `$ for small $`k_x`$. Further, since the dispersion is invariant under the substitution $`k_xk_x`$ we have $`A,p_i,q_ik_x^2`$ for small $`k_x`$. Thus we are led to write $`q_1(k_x)=x_0+\beta k_x^2`$, with positive $`\beta `$. With these preparations we rewrite Eq. (38)

$`\alpha `$ $`=`$ $`\underset{k_x0}{lim}k_x{\displaystyle _0^{x_0}}{\displaystyle \frac{2\beta k_x}{(x\underset{q_1}{\underset{}{x_0\beta k_x^2}})^2}}g(x;k_x)𝑑x`$ (41)
$`=`$ $`\underset{k_x0}{lim}2{\displaystyle _0^{\frac{x_0}{\beta k_x^2}}}{\displaystyle \frac{g(x_0\beta k_x^2;k_x)}{(1+y)^2}}𝑑y,`$ (42)
where we substituted $`x=x_0\beta k_x^2y`$ and $`g(x;𝐤)`$ for $`A(k_x)P(x;k_x)(xq_1)/Q(x;k_x)`$. Taking the limit yields

$$\alpha =2g(x_0;0)_0^{\mathrm{}}\frac{1}{(1+y)^2}𝑑y=2g(x_0;0).$$
(43)
A straightforward calculation shows that

$$g(x_0;0)=A(0)\frac{P(x_0;0)}{_xQ(x_0,0)},$$
(44)
which we can easily calculate. Fig. 7 shows the exponent $`\alpha `$ as function of the angle $`\varphi `$ measured from the positive $`x`$-axis to the negative $`y`$-axis.

There is only one ratio $`x/y`$ for which $`\alpha =1`$. At this point the dispersion vanishes with finite spin wave velocity ($`\varphi 0.081(1)`$ or $`x0.679(1)`$, $`y0.055(1)`$).

Let us turn to a comparison of the theoretical dispersion to experimental data for SrCu$`2`$(BO3)2. To fit the dispersion to experimental data we make use of the parameter dependance of our results by requiring the curves to go through certain points and solve the resulting set of equations. At $`𝐤=(0,0)`$ ESR , FIR  and INS  data suggest a value of $`\omega (0,0)=2.98`$meV. At finite $`𝐤`$ we have to rely on the INS measurement, which contain rather large errors. In Fig. 8 we show the INS data (bullets and error bars) and two of our fitted curves. For the parameter values given we get an excellent agreement.

Because of the flatness of the dispersion and the comparably large error bars it is not possible to fix the model parameters unambiguously . As sketched in Fig. 8 one can lower $`x`$ and enlarge $`J_1`$ and $`|y|`$ without loosing reasonable agreement. On the basis of the one-triplet dispersion it is not possible to fix the model parameters quantitatively. Our investigation in the two-triplet sector  (where $`y`$ is disregarded, since it is probably not very important) enables us to fix the parameters to the intervals ($`y=0`$,) $`x=0.603(3)`$ and $`J_1=6.16(10)`$meV. The corresponding one-triplet dispersion is the solid curve in Fig. 8.

## V Conclusion

In summary, we have demonstrated the utility of our perturbation method introduced earlier . By the means of a continuous unitary transformation the original Hamiltonian is mapped onto a block diagonal, energy quanta conserving, effective Hamiltonian $`_{\mathrm{eff}}`$. This allows to do the calculations in different energy sectors separately.
Here we used this method to develop a picture of excitation processes in a model for SrCu2(BO3)2 and calculated the series expansion of the one-magnon dispersion $`\omega `$ about the limit of isolated dimers as a 15th order polynomial in the model parameters.

We showed that the dispersion decomposes into two nearly degenerate branches with a splitting proportional to $`x^{10}`$. At $`𝐤=\mathrm{𝟎}`$ and at the border of the magnetic Brillouin zone the two branches fall onto each other. By making use of our detailed analysis of the models symmetries we showed that point group operations alone cannot explain these degenaracies. In fact we showed that the model’s space group symmetries have to be taken into account.
Moreover we analysed the critical behaviour of $`\omega `$ for small $`|𝐤|`$ at various ratios of $`x`$ and $`y`$. It is found that the dispersion vanishes only for $`x0.679(1)`$ and $`y0.055(1)`$ with a finite spin velocity.
Finally, we fitted our dispersion to INS data obtained by Kageyama et al.  and found that it is not possible to fix the model parameters unambiguously from the one-triplet data alone due to very large experimental error bars. The parameter values determined by our investigations in the two-triplet sector (: $`y=0`$, $`x=0.603`$ , $`J=6.67`$meV), however, lead to a one-magnon dispersion agreeing nicely with the experimental findings.
The authors would like to thank Alexander Bühler for helpful discussions and the Regional Computing Centre of the University of Cologne (RRZK) and Dr. Kim Yong Taik in particular for generous assistance in our large scale computations. This work was supported by the DFG-Schwerpunkt 1073 and by the SFB 341.

## A

First we show that $`dt_𝐫`$=0 if $`r_1=r_2`$.
With Eq. (22) (the $`\sigma _v`$ symmetry) we have

$$t_𝐬^{o(𝐫)}=t_{s_2,s_1}^{o(𝐫(0,1))}.$$
(A1)
Splitting both sides according to Eq. (14) we get

$$\overline{t_𝐬}+e^{iQ𝐫}dt_𝐬=\overline{t}_{s_2,s_1}+e^{iQ(𝐫(0,1))}dt_{s_2,s_1},$$
(A2)
leading to

$`\overline{t_𝐬}+dt_𝐬`$ $`=`$ $`\overline{t}_{s_2,s_1}dt_{s_2,s_1},\text{ for }r\text{ even and}`$ (A3)
$`\overline{t_𝐬}dt_𝐬`$ $`=`$ $`\overline{t}_{s_2,s_1}+dt_{s_2,s_1},\text{ for }r\text{ odd.}`$ (A4)
Taking the difference of both equations yields

$$dt_{s_1,s_2}=dt_{s_2,s_1},$$
(A5)
which proves the assertion. In particular $`dt_\mathrm{𝟎}=0`$.
We now show that $`dt_𝐫=0`$ if $`r_1+r_2`$ is an odd number.
According to Eq. (13) we have

$`t_𝐫^{}^{o(𝐫)}`$ $`=`$ $`𝐫+𝐫^{}|_{\mathrm{eff}}|𝐫`$ (A6)
$`=`$ $`𝐫|_{\mathrm{eff}}|𝐫+𝐫^{}=t_𝐫^{}^{o(𝐫+𝐫^{})}`$ (A7)
$`=`$ $`t_𝐫^{}^{o(𝐫+𝐫^{})}`$ (A8)
since $`_{\mathrm{eff}}`$ has only real matrix elements in this basis (cf. Eq. (3)). The last equality follows from Eq. (20). Splitting both sides according to Eq. (14) yields

$`\overline{t_𝐫^{}}+e^{iQ𝐫}dt_𝐫^{}`$ $`=`$ $`\overline{t_𝐫^{}}+e^{iQ𝐫}e^{iQ𝐫^{}}dt_𝐫^{}`$ (A9)
$`dt_𝐫^{}`$ $`=`$ $`e^{iQ𝐫^{}}dt_𝐫^{}`$ (A10)
$`dt_𝐫^{}`$ $`=`$ $`0\text{ for }r_1^{}+r_2^{}\text{ odd.}`$ (A11)
## B Degeneracy of $`\omega `$

The dispersion (35) will be two-fold degenerate if the square root vanishes. We will use the shorthand $`𝐫`$ even/odd for $`r_1+r_2`$ even/odd.
I) $`𝐤=0`$

A) $`b\stackrel{!}{=}0`$
For $`𝐤=0`$ Eq. (34) gives

$`{\displaystyle \frac{b}{2}}`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{𝐫>0}{𝐫\text{ even}}}{}}dt_𝐫={\displaystyle \underset{\genfrac{}{}{0pt}{}{r_1=0}{\genfrac{}{}{0pt}{}{r_2>0}{𝐫\text{ even}}}}{}}dt_𝐫+{\displaystyle \underset{\genfrac{}{}{0pt}{}{r_1>0}{\genfrac{}{}{0pt}{}{r_2=0}{𝐫\text{ even}}}}{}}dt_𝐫+{\displaystyle \underset{\genfrac{}{}{0pt}{}{r_1>0}{\genfrac{}{}{0pt}{}{r_2>0}{𝐫\text{ even}}}}{}}dt_𝐫+{\displaystyle \underset{\genfrac{}{}{0pt}{}{r_1>0}{\genfrac{}{}{0pt}{}{r_2<0}{𝐫\text{ even}}}}{}}dt_𝐫`$ (B1)
$`=`$ $`I_1+I_2+I_3+I_4.`$ (B2)
We rewrite the first and the third sum on the right hand side of Eq. (B2)

$$I_1=\underset{\genfrac{}{}{0pt}{}{r_1>0}{\genfrac{}{}{0pt}{}{r_2=0}{𝐫\text{ even}}}}{}dt_{r_2,r_1}\text{ and }I_3=\underset{\genfrac{}{}{0pt}{}{r_1>0}{\genfrac{}{}{0pt}{}{r_2<0}{𝐫\text{ even}}}}{}dt_{r_2,r1}.$$
(B3)
From Eq. (23) (the $`R`$ symmetry) we deduce

$`dt_{r_1,r_2}`$ $`=`$ $`dt_{r_2,r_1},`$ (B4)
and see that $`I_1=I_2`$ and $`I_3=I_4`$. Hence we find $`b=0`$.

B) $`(a_0a_1)\stackrel{!}{=}0`$
For $`𝐤=0`$ we have (cf. Eq. (33))

$$\frac{a_0a_1}{2}=\underset{𝐫>0}{}\overline{t}_𝐫\left[1\mathrm{cos}(\pi (r_1+r_2))\right]=2\underset{\genfrac{}{}{0pt}{}{𝐫>0}{𝐫\text{ odd}}}{}\overline{t}_𝐫,$$
(B5)
and by recalling that $`dt_𝐫=0`$ for $`𝐫`$ odd Eq. (23) yields

$$\overline{t}_{r_1,r_2}=\overline{t}_{r_2,r_1},$$
(B6)
so that we can use the same splitting as in Eq. (B2). Thus the energy-degeneracy at $`𝐤=0`$ is due to the rotational symmetry $`R`$ as defined in Sec III.
II) $`k_x+k_y=\pi `$

A) $`b\stackrel{!}{=}0`$
Making use of $`dt_𝐫=dt_𝐫`$, which follows from Eq. (20), we have for $`k_y=\pi k_x`$

$`{\displaystyle \frac{b}{4}}`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{𝐫}{𝐫\text{ even}}}{}}dt_𝐫\mathrm{cos}(k_x(r_1r_2)+\pi r_2)`$ (B7)
$`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{r_1,r_2}{𝐫\text{ even}}}{}}dt_{r_2,r_1}\mathrm{cos}(k_x(r_2r_1))\mathrm{cos}(\pi r_1).`$ (B8)
In the last step we choose to rearrange the sum and observe that if $`𝐫`$ is even we have $`r_1`$ and $`r_2`$ both odd or both even. In both cases the identity

$$\mathrm{cos}(\pi r_1)=\mathrm{cos}(\pi r_2)$$
(B9)
holds. Inserting relation (A5) in the last row of Eq. (B8) we end up with

$$\frac{b}{4}=\underset{\genfrac{}{}{0pt}{}{𝐫}{𝐫\text{ even}}}{}dt_𝐫\mathrm{cos}(k_x(r_1r_2))\mathrm{cos}(\pi r_2),$$
(B10)
resulting in $`b=b`$ and thus $`b=0`$.

B) $`(a_0a_1)\stackrel{!}{=}0`$
Analogously to I) B) we have here

$`{\displaystyle \frac{a_0a_1}{4}}`$ $`=`$ $`{\displaystyle \underset{𝐫}{}}\overline{t}_𝐫\left[\mathrm{cos}(\mathrm{𝐤𝐫})\mathrm{cos}(\mathrm{𝐤𝐫}+\pi (r_1+r_2))\right]`$ (B11)
$`=`$ $`2{\displaystyle \underset{𝐫\text{ odd}}{}}\overline{t}_𝐫\mathrm{cos}(k_x(r_1r_2)+\pi r_2)`$ (B12)
$`=`$ $`2{\displaystyle \underset{𝐫\text{ odd}}{}}\overline{t}_{r_2,r_1}\mathrm{cos}(k_x(r_2r_1))\mathrm{cos}(\pi r_1)`$ (B13)
$`=`$ $`2{\displaystyle \underset{𝐫\text{ odd}}{}}\overline{t}_{r_1,r_2}\mathrm{cos}(k_x(r_1r_2))\mathrm{cos}(\pi r_2),`$ (B14)
where the last but one equality follows from Eq. (22) (the $`\sigma _v`$ symmetry), i.e $`\overline{t}_{r_1,r_2}=\overline{t}_{r_2,r_1}`$, and from the fact that if $`𝐫`$ is odd we have that $`r_1`$ odd and $`r_2`$ even or $`r_1`$ even and $`r_2`$ odd. From that we see $`\mathrm{cos}(\pi r_1)=\mathrm{cos}(\pi r_2)`$.
The degeneracy over the remaining three borders of the magnetic Brillouin zone can be shown analogously. It is interesting to note, that the calculations necessarily involved glide line operations. The degeneracies can thus not be explained by considering point group symmetries only.


format: text

text:
# Droplet Formation in Quark-Gluon Plasma at Low Temperatures and High Densities

## I Introduction

Our understanding of the QCD matter at low temperatures ($`T`$) and high baryon number densities ($`n_B`$) has developed rapidly over the past few years; for recent reviews, see . A likely scenario is that in this region, quarks form Cooper pairs and new condensates develop. One of the most interesting questions today concerns the detailed mapping of the QCD phase diagram, especially in the above region where QCD may exhibit a colour superconducting phase. In this context, the light quark masses ($`m_{u,d,s}`$) are expected to play an important role. Figure 1 shows two likely scenarios corresponding to $`m_sm_{u,d}0`$ and $`m_s=m_{u,d}0`$. The transition between the quark-gluon plasma (QGP) phase and the colour superconducting phase is expected to be first order, in QCD .

In this paper, we study the first-order phase transition between the QGP phase and the colour superconducting phase (Fig. 1). Our approach is phenomenological and is similar in spirit to that adopted by Côté and Kharchenko . They, however, studied the Bose-Einstein condensation in a gas of trapped atomic hydrogen at temperatures $`50\mu K`$ and densities $`10^{14}`$ cm-3. Temperatures and densities that we shall consider here are, of course, vastly different, and QGP has very little in common with the atomic hydrogen gas. Interestingly, however, as we shall see, the same physics of activation energy barrier governs the formation of droplets in both cases.

It is important to get a handle on the cold and dense quark matter because of its possible relevance to the neutron-star core. While the region of asymptotic densities is amenable to rigorous calculations, the same cannot be said about the region of intermediate densities studied here. Both numerical (lattice) and direct experimental studies of this region of the QCD phase diagram are out of reach at the moment. Hence a phenomenological study such as the present one is appropriate and worthwhile. It has allowed us to make quantitative statements about the schematic phase diagram in Fig. 1.

Alford et al.  have discussed droplet formation in the transition between the hadronic matter phase and the 2SC phase (Fig. 1), at $`T=0`$, as a function of density. The chiral symmetry is restored inside their droplets (which they call nucleons) but is broken outside. Neergaard and Madsen  also have studied the droplet formation but in the QGP to hadronic matter transition occurring at low baryon number densities. In this paper, we have considered the QGP to colour-superconducting phase transition occurring at high $`n_B`$. To our knowledge, droplet formation in this transition has not been discussed previously in the literature. Our droplets are the opposite of those in  in the sense that the chiral symmetry is broken inside our droplet, but is restored outside.

## II Model

Consider a gas of weakly interacting quarks, antiquarks and gluons in a volume $`V`$. We take $`q`$ and $`\overline{q}`$ to be the massless current quarks (Fig. 1). We are interested in low temperatures ($`T`$) and high number densities ($`n`$) such that $`n\lambda ^31`$ where $`\lambda `$ is the thermal wave length (see Appendix). We also need the baryon number density $`n_B=(n_qn_{\overline{q}})/3`$ to be large compared to that for the normal nuclear matter, namely 0.17 fm-3.

If the interactions among the particles are treated to the lowest order in $`\alpha _s`$, the strong coupling constant (Fig. 2), the baryon number density is given by ($`c=\mathrm{}=1`$)

$$n_B=\left(1\frac{2\alpha _s}{\pi }\right)\frac{d_q}{3}\left[\frac{\mu _q(kT)^2}{6}+\frac{\mu _q^3}{6\pi ^2}\right],$$
(1)
where $`d_q=`$ 2 (spin) $`\times `$ 3 (colour) $`\times `$ 3 (flavour) = 18 is the quark degeneracy factor, $`k`$ is the Boltzmann constant and $`\mu _q`$ is the quark chemical potential . For a given baryon number density $`n_B`$ and temperature $`T`$, Eq. (1) uniquely fixes $`\mu _q`$. This in turn determines the total energy and pressure of the gas as

$`E_{gas}`$ $`=`$ $`\left(1{\displaystyle \frac{15\alpha _s}{4\pi }}\right)d_gV{\displaystyle \frac{\pi ^2(kT)^4}{30}}`$ (2)
$`+`$ $`d_qV\left[\left(1{\displaystyle \frac{50\alpha _s}{21\pi }}\right){\displaystyle \frac{7\pi ^2(kT)^4}{120}}+\left(1{\displaystyle \frac{2\alpha _s}{\pi }}\right)\left({\displaystyle \frac{\mu _q^2(kT)^2}{4}}+{\displaystyle \frac{\mu _q^4}{8\pi ^2}}\right)\right],`$ (3)
$`P_{gas}`$ $`=`$ $`E_{gas}/3V.`$ (4)
Here $`d_g=`$ 2 (spin) $`\times `$ 8 (colour) = 16 is the gluon degeneracy factor. It was possible to obtain the above analytic expressions because we have considered the gas of quarks as well as antiquarks, rather than a gas of quarks alone . At low $`T`$ and high $`n_B`$, the $`\mu _q`$ terms in the pressure dominate and the gluonic pressure arising from the first term in Eq. (2), is quite negligible.

We now consider the first-order phase transition from QGP to the superconducting phase as the temperature is lowered (Fig. 1). In the latter phase, as a result of the attractive one-gluon-exchange interaction between quarks in the colour $`\overline{\mathrm{𝟑}}`$ channel, quarks near the Fermi surface tend to pair up as bosons. These are the Cooper pairs which we assume to be massless. Density fluctuations in the metastable QGP gas may provide suitable conditions for nucleation and growth of a droplet of the superconducting phase. The high-temperature phase (QGP) is present around the droplet and the low-temperature superconducting phase exists inside the droplet. We will presently explore the energetics for the formation of such a droplet by calculating the energy barrier between the two phases.

We have given above the expressions for number density, energy and pressure of the QGP; see Eqs. (1)-(3). We now derive corresponding expressions for the superconducting phase inside the droplet. Before we begin, it is necessary to recall important differences between the BCS superconductor and the colour superconductor. (a) In the BCS theory of superconductivity the ratio of the energy gap ($`\mathrm{\Delta }`$) to the Fermi kinetic energy for a typical metal is of the order of $`10^4`$. In contrast, in the problem under consideration, the gap $`\mathrm{\Delta }`$ is estimated to vary between several tens of MeV to about 100 MeV , and the quark chemical potential $`\mu _q`$ appearing in Eq. (1) takes values between $`250`$ and $`600`$ MeV for the $`n_B`$ and $`T`$ considered here. (b) In the BCS superconductor, the Debye energy $`\omega _D`$ is much smaller than the chemical potential $`\mu `$, while in the colour superconductor, the two may be of the same order of magnitude. (c) Finally, if $`n`$ is the density of electrons and $`\xi `$ the Pippard coherence length which measures the spatial extension of the Cooper pair wave function, then in the BCS theory, $`n\xi ^31`$. In the present case, on the other hand, $`\xi 1/(\pi \mathrm{\Delta })0.6`$ fm for $`\mathrm{\Delta }=100`$ MeV, and $`n\xi ^3<1`$ for the densities considered in this paper. This suggests that the $`qq`$ pairs are rather compact objects unlike Cooper pairs of electrons in a metal. For a good discussion of these and related issues, see a review article by Kerbikov . We shall make use of the above points while deriving expressions for thermodynamics quantities of the droplet. For the sake of simplicity, we shall assume that all quarks in the energy interval $`(\mu _q\mathrm{\Delta },\mu _q)`$ form Cooper pairs and all other quarks remain unpaired, and shall ignore other interactions.

Consider a spherical region of radius $`R_d`$, volume $`V_d`$, containing $`N_1`$ quarks out of which $`N_2`$ remain unpaired and $`N_3`$ undergo Cooper pairing, when the pairing interaction is “switched on”. We have

$`N_1`$ $`=`$ $`{\displaystyle \frac{d_qV_d}{(2\pi )^3}}{\displaystyle _0^{\mathrm{}}}4\pi p^2f𝑑p,`$ (5)
$`N_2`$ $`=`$ $`N_1{\displaystyle \frac{d_qV_d}{(2\pi )^3}}{\displaystyle _{\mu _q\mathrm{\Delta }}^{\mu _q}}4\pi p^2f𝑑p,`$ (6)
$`N_3`$ $`=`$ $`N_1N_2,`$ (7)
where $`f=\left[\mathrm{exp}\beta (p\mu _q)+1\right]^1`$ with $`\beta =1/(kT)`$. The corresponding energies $`E_{1,2,3}`$ are

$`E_1`$ $`=`$ $`{\displaystyle \frac{d_qV_d}{(2\pi )^3}}{\displaystyle _0^{\mathrm{}}}4\pi p^3f𝑑p,`$ (8)
$`E_2`$ $`=`$ $`E_1{\displaystyle \frac{d_qV_d}{(2\pi )^3}}{\displaystyle _{\mu _q\mathrm{\Delta }}^{\mu _q}}4\pi p^3f𝑑p,`$ (9)
$`E_3`$ $`=`$ $`\left({\displaystyle \frac{N_3}{2}}\right)\left({\displaystyle \frac{\pi }{R_d}}\mathrm{\Delta }\right).`$ (10)
In the last equation $`(N_3/2)`$ is the number of Cooper pairs (bosons) in the droplet and $`(\pi /R_d)`$ is the energy of the ground state of the single-particle spectrum of a spherical cavity of radius $`R_d`$. Due to the small size of the droplet, even the lowest excited state has energy which is very large compared to the temperature. Therefore all the $`qq`$ pairs tend to occupy the ground state . Next, entropies $`S_{1,2,3}`$ are given by

$`S_1`$ $`=`$ $`{\displaystyle \frac{d_qV_d}{(2\pi )^3}}{\displaystyle _0^{\mathrm{}}}4\pi p^2𝑑p\left[f\mathrm{ln}f+(1f)\mathrm{ln}(1f)\right],`$ (11)
$`S_2`$ $`=`$ $`S_1+{\displaystyle \frac{d_qV_d}{(2\pi )^3}}{\displaystyle _{\mu _q\mathrm{\Delta }}^{\mu _q}}4\pi p^2𝑑p\left[f\mathrm{ln}f+(1f)\mathrm{ln}(1f)\right],`$ (12)
$`S_3`$ $`=`$ $`0.`$ (13)
$`S_3`$ can be taken to be zero for the reasons given above. Finally, the pressures $`P_{1,2,3}`$ are calculated as follows:

$$P_1=E_1/(3V_d).$$
(14)
The thermodynamic potential $`\mathrm{\Omega }`$ for massless quarks is given by

$`\mathrm{\Omega }=kT{\displaystyle \mathrm{ln}\left[1+\mathrm{exp}\beta (\mu _qp)\right]}.`$

Hence

$`P_2`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }}{V_d}}={\displaystyle \frac{kTd_q}{(2\pi )^3}}\left({\displaystyle _0^{\mathrm{}}}{\displaystyle _{\mu _q\mathrm{\Delta }}^{\mu _q}}\right)4\pi p^2dp\mathrm{ln}\left[1+\mathrm{exp}\beta (\mu _qp)\right]`$ (15)
$`=`$ $`{\displaystyle \frac{E_1}{3V_d}}{\displaystyle \frac{kTd_q}{2\pi ^2}}{\displaystyle _{\mu _q\mathrm{\Delta }}^{\mu _q}}p^2𝑑p\mathrm{ln}\left[1+\mathrm{exp}\beta (\mu _qp)\right].`$ (16)
Performing the last integration by parts and rearranging the terms, one finally gets

$$P_2=\frac{E_2}{3V_d}\frac{kTd_q}{6\pi ^2}\{\mu _q^3\mathrm{ln}2(\mu _q\mathrm{\Delta })^3\mathrm{ln}\left[1+\mathrm{exp}(\beta \mathrm{\Delta })\right]\}.$$
(17)
The quantal pressure $`P_3`$ is given by

$$P_3=(E_3/V_d)_{N_3}=N_3/8R_d^4.$$
(18)
It is easy to verify that if $`\mathrm{\Delta }`$ is set equal to zero, i.e. if no pairing is allowed to take place, $`N_3`$, $`E_3`$, $`P_3`$ vanish and $`N_2`$, $`E_2`$, $`S_2`$, $`P_2`$ become equal to $`N_1`$, $`E_1`$, $`S_1`$, $`P_1`$, respectively. Contribution of antiquarks to the above thermodynamic quantities is smaller by many orders of magnitude. Although we have included it in our calculations its effect is negligible.

The droplet pressure $`P_d`$ is given by

$$P_d=P_2+P_3+\overline{P}_2+\overline{P}_3+P_{gluons},$$
(19)
where $`P_2`$ and $`P_3`$ are as in Eqs. (14)-(15), $`\overline{P}_2`$ and $`\overline{P}_3`$ are the corresponding pressures due to antiquarks and $`P_{gluons}`$ is the pressure due to gluons in the droplet.

Consider once again the spherical region of radius $`R_d`$. Initially, i.e. before the pairing interaction is switched on, there are $`N_1`$ quarks, $`\overline{N}_1`$ antiquarks and a certain number of gluons in it. In the final state, i.e. when the pairing interaction is on, there are $`N_2`$ quarks, $`(N_3/2)`$ quark pairs, $`\overline{N}_2`$ antiquarks, $`(\overline{N}_3/2)`$ antiquark pairs and the same number of gluons as before. The activation energy ($`A`$) is defined as the difference between the Helmholtz free energies ($`F`$) in the final and initial states. Since the gluonic contribution cancels out, one gets

$`A`$ $`=`$ $`F(N_2)+F(N_3/2)F(N_1)+F(\overline{N}_2)+F(\overline{N}_3/2)F(\overline{N}_1)`$ (20)
$`=`$ $`E_2TS_2+E_3E_1+TS_1+\overline{E}_2T\overline{S}_2+\overline{E}_3\overline{E}_1+T\overline{S}_1.`$ (21)
In the limit $`\mathrm{\Delta }0`$, this $`A`$ vanishes.

The reader may wonder why surface tension of the droplet makes no appearance in our formalism. In the classical theory of nucleation, there are two competing contributions to the energy of a droplet: droplet formation reduces bulk free energy but only at the cost of surface energy. The former contribution wins if the droplet radius exceeds a certain critical value, while the surface tension ensures that a droplet with a smaller radius shrinks and disappears. In the formalism presented here, energy of the droplet is calculated quantum mechanically. It is clear from Eq. (9) that small droplets have large quantal energy and tend to dissolve back into the QGP gas. This is similar to the effect of the surface tension in the classical theory. Moreover, the thermal wavelength ($`\lambda `$) of the quarks is so large that a classification of particles as those localized on the surface and those in the interior breaks down. In short, we did not neglect the surface energy; it is included in the total free energy of the quantal droplet.

## III Results and Discussion

The precise boundaries of the various phases in Fig. 1 are not known. Since the baryon number density ($`n_B`$) of the normal nuclear matter is $`0.17`$ fm-3, we have performed numerical calculations for $`n_B=`$ 0.5 to 3 fm-3 and temperature $`T=`$ 0 to 100 MeV. The resulting quark chemical potential $`\mu _q`$ varies over $`250600`$ MeV and the baryon chemical potential***This definition is applicable only for a gas of quarks, and not for the normal nuclear matter. $`\mu _B3\mu _q`$ varies over $`7501800`$ MeV. This fully covers the densities likely to be seen in the neutron-star cores. We have considered four different values of the superconducting gap parameter $`\mathrm{\Delta }`$, namely, 1, 10, 50 and 100 MeV.

Since we are interested in a gas of weakly interacting quarks, antiquarks and gluons, we take $`\alpha _s`$ to be zero. Small values of $`\alpha _s`$ do not affect the results in Eqs. (1)-(3) significantly. We first choose a set of values for $`n_B`$ and $`T`$, and solve Eq. (1) numerically to determine $`\mu _q`$. This $`\mu _q`$ is then used to evaluate $`P_{gas}`$ as given by Eq. (3). Also, once $`\mu _q`$ is determined, $`n_q`$ and $`n_{\overline{q}}`$ can be evaluated separately by using

$`n_q`$ $`=`$ $`{\displaystyle \frac{d_q}{(2\pi )^3}}{\displaystyle _0^{\mathrm{}}}4\pi p^2f𝑑p,`$ (22)
$`n_{\overline{q}}`$ $`=`$ $`n_q(\mu _q\mu _q)`$ (23)
(Above integrations cannot be done analytically, and so it is not possible to get closed expressions like Eq. (1), for $`n_q`$ and $`n_{\overline{q}}`$.) Next, we choose a value for $`\mathrm{\Delta }`$, and treating the droplet radius $`R_d`$ as a free parameter, calculate various other quantities; see Eqs. (4)-(17). Some of these results are shown in Figs. 3-4.

In our formalism, as in the classical theory, the activation energy when plotted as a function of the radius (Fig. 3), initially increases from zero, attains a maximum, then starts decreasing monotonically eventually becoming negative for a sufficiently large value of the radius. In other words, the activation energy exhibits a barrier. The barrier originates as follows. Consider Eq. (17). As already stated, the antiquark terms are extremely small. The terms $`TS_1`$ and $`TS_2`$ are also negligible and in fact they tend to cancel each other. That leaves three terms, namely $`E_1,E_2`$ and $`E_3`$. For small radii, $`E_2`$ and the quantal energy term $`E_3`$ dominate pushing the activation energy curve upwards. For large radii, the $`E_1`$ term dominates pulling the curve to the negative domain. As expected, the barrier becomes wider if $`T`$ is raised or if $`n_B`$ is lowered. The critical radii can be read off from Fig. 3. Droplets with larger radii readily grow because that tends to reduce the activation energy. This in turn easily brings about the transition from QGP to the superconducting phase.

Figure 4 shows the ratio of the droplet pressure to the pressure of the surrounding gas, for $`\mathrm{\Delta }=100`$ MeV. This ratio is significantly larger than unity for small droplets, which will tend to push the radii to larger values, thereby helping to bring about the QGP to superconducting phase transition. In the limit $`\mathrm{\Delta }0`$, this ratio tends to unity and is thus incapable of pushing the droplet radii to larger values.

It is interesting to know the properties of the droplet when the pressure inside the droplet (Eq. 16) equals the pressure of the surrounding QGP gas (Eq. 3). These results are shown in Figs. 5-8. The condition $`P_d=P_{gas}`$ fixes $`R_d`$. We solve the equation $`P_d=P_{gas}`$ numerically to determine $`R_d`$, for various values of $`n_B`$, $`T`$ and $`\mathrm{\Delta }`$. $`n_B`$ was varied from 0.5 to 3 fm-3 in the steps of 0.5 fm-3. Figure 5 shows the results for $`R_d`$. As $`T`$ and/or $`n_B`$ are increased, $`R_d`$ decreases. However, $`R_d`$ is found to be rather insensitive to the variation in $`n_B`$. Although we have varied $`T`$ between 0 and 100 MeV, temperatures relevant for neutron stars are close to zero. Because of the small value of $`R_d`$, the energy of even the lowest excited state inside the droplet turns out to be several hundreds of MeV which is large compared to the temperature, effectively forcing all Cooper pairs (bosons) to occupy the ground state. Note also that at low temperatures, the droplet radius $`R_d`$ is larger than the Pippard coherence length $`\xi 1/(\pi \mathrm{\Delta })`$.

Figure 6 shows the activation energy $`A`$. First see the upper half of Fig. 6. For a fixed $`n_B`$, as $`T`$ increases, $`A`$ changes its sign from negative to positive. This is expected because, at higher $`T`$, droplet formation is energetically unfavourable. Secondly, as $`n_B`$ increases (see curves labelled a to f ), $`A`$ remains negative up to higher and higher $`T`$. This is also expected because at higher $`n_B`$, the transition to the superconducting phase is facilitated. Now see the lower half of Fig. 6. It shows the $`A=0`$ contour, in the $`Tn_B`$ plane. This is also the phase transition curve separating the QGP phase at high $`T`$ and the superconducting phase at low $`T`$: Recall $`AF_{SC}F_{QGP}`$ where $`F`$ are the free energies of the two phases. Both $`F_{SC}`$ and $`F_{QGP}`$ decrease as $`T`$ increases and cross each other at some temperature, say $`T_c`$ (Fig. 7). Since the system follows the path of minimum $`F`$, there is a kink in $`F(T)`$. In other words, the slope $`dF/dT=S`$ has a discontinuity, indicating a first-order phase transition at $`T_c`$. Thus the $`A=0`$ contour serves to demarcate the region of the superconducting phase. Note that as $`n_B`$ increases, the $`A<0`$ region in Fig. 6 survives up to higher and higher $`T`$.

We now present an estimate of the latent heat ($`L`$) associated with the first-order phase transition considered in this paper. The transition temperature $`T_c`$ for any given $`n_B`$ can be determined from Fig. 6. The latent heat is then estimated as $`L=T_c\times (S_1(T_c)S_2(T_c))`$. For $`\mathrm{\Delta }=100`$ MeV, this turns out to be $`150`$ to $`250`$ MeV and it increases with $`n_B`$. In the limit $`\mathrm{\Delta }0`$, the two free-energy curves in Fig. 7 coincide, $`A`$ vanishes, $`S1=S2`$, there is no discontinuity in entropy, and the latent heat is zero. There is only one phase namely QGP, and there is no question of any phase transition.

Finally, Fig. 8 shows 3-dimensional surface plots of the (normalized) activation energy $`A^{}`$ defined as $`A^{}(F_{SC}F_{QGP})/(F_{SC}+F_{QGP})`$ as a function of $`n_B`$ and $`T`$, for $`\mathrm{\Delta }=1`$ and 100 MeV. It is clear that in the limit $`\mathrm{\Delta }0`$, the activation energy surface is flat and there is no energy advantage for the QGP to superconducting phase transition. Indeed there would be only one phase, namely QGP, which would extend up to the $`T=0`$ axis. On the other hand, for $`\mathrm{\Delta }=100`$ MeV, the surface is such that the system would like to “roll down” the slope thereby bringing about the phase transition.

It is interesting to compare the activation energy results for the QCD phase transition, obtained here with those for the Bose-Einstein condensation in a gas of trapped atomic hydrogen, presented in . In  two gas densities were considered: $`10^{14}`$ cm-3 and $`7\times 10^{13}`$ cm-3 both at the same temperature $`50\mu K`$. We have repeated their calculations for $`T=4555\mu K`$. Their activation energy exhibits a barrier, and they denote the maximum value of the activation energy by $`A_0`$. Consider the ratio $`A_0/kT`$ at the peak position and the corresponding number $`\stackrel{~}{N}_0`$ of particles in the droplet. In  these two quantities are highly sensitive to the changes in the gas density and temperature.There is a misprint in : in Fig. 1, the y-axis label should read $`A(10^{10})`$ a.u. For the inset, it is $`A(10^8)`$ a.u. For instance, for the above tiny change in the gas density, the ratio $`A_0/kT`$ changes from 4.8 to 104.1, and $`\stackrel{~}{N}_0`$ changes from 18 to 3862. Similar sensitivity was seen when $`T`$ was changed by a few $`\mu K`$. In contrast, the results presented here are relatively stable with respect to changes in $`n_B`$ and $`T`$.

In conclusion, we have studied hitherto unexplored first-order phase transition from QGP to a colour-superconducting phase, and have presented a simple model of the transition with very few assumptions. It provides numerical support to the prevalent ideas of phase structure of QCD in the low-temperature, high-baryon-number-density regime. We have studied properties of the droplets of the superconducting phase, formed at various points in the $`Tn_B`$ plane, for several values of the gap parameter. We have been able to estimate the extent of the superconducting phase in this plane and the latent heat of the transition. At low temperatures, droplet radii are found to be larger than the Pippard coherence length which measures the spatial extent of the $`qq`$ Cooper pairs.

RSB acknowledges the hospitality of the McMaster University where this work was initiated. We thank Mark Alford and Mustansir Barma for their comments on this work. This research was partially supported by a grant from NSERC (Canada).

Appendix

If $`f(p)`$ is the distribution function of a gas in equilibrium, the number density $`n`$ is given by $`n=f(p)d^3p/(2\pi \mathrm{})^3`$. For a nonrelativistic classical gas at temperature $`T`$, one has

$`f=\mathrm{exp}\left[\beta \left({\displaystyle \frac{p^2}{2m}}\mu \right)\right],`$

where $`\beta =1/kT`$ and $`\mu `$ is the chemical potential. Hence

$`n=\mathrm{exp}(\beta \mu )(2\pi \mathrm{}^2/mkT)^{3/2}\mathrm{exp}(\beta \mu )\lambda ^3,`$

where $`\lambda `$ is the thermal wave length.

For an extreme relativistic gas of massless classical particles,

$`f=\mathrm{exp}\left[\beta (pc\mu )\right],`$

and hence

$`n=\mathrm{exp}(\beta \mu )(c\mathrm{}\pi ^{2/3}/kT)^3.`$

We define the thermal wave length in this case as $`\lambda =c\mathrm{}\pi ^{2/3}/kT,`$ so that the fugacity $`\mathrm{exp}(\beta \mu )`$ once again equals $`n\lambda ^3`$.

FIGURE CAPTIONS

FIG. 1. Schematic phase diagrams of QCD: $`T`$ is the temperature and $`\mu _B`$ is the baryon chemical potential. $`E`$ is the critical point at which the line of first-order phase transitions ends. $`2SC`$ and $`CFL`$ (Colour-Flavour-Locked) are 2-flavour-like and 3-flavour-like colour superconducting phases, respectively. The dashed line denotes the critical temperature at which quark-quark pairing vanishes. The thick arrow indicates the first-order phase transition considered in this paper. For the normal nuclear matter at $`T=0`$, $`\mu _B`$ should be defined as $`m_N+\mathrm{}^2k_F^2/2m_N`$ where $`k_F1.4`$ fm-1  and $`m_N`$ is the nucleon mass, which gives $`\mu _B980`$ MeV. (Figure adapted from .)

FIG. 2. Feynman diagrams contributing to the equation of state of quark-gluon plasma in order $`\alpha _s`$. Curly lines: gluons, dashed lines: ghost particles, solid lines: quarks.

FIG. 3. Activation energy $`A`$ vs droplet radius $`R_d`$ for various gap values $`\mathrm{\Delta }`$. Baryon number density $`n_B`$ (fm-3) and temperature $`T`$ (MeV) are: (a) 0.5 and 5 (b) 0.5 and 100 (c) 3 and 5 (d) 3 and 100, respectively.

FIG. 4. Ratio of the pressure of the droplet to the pressure of the surrounding QGP gas vs droplet radius $`R_d`$, at gap $`\mathrm{\Delta }=100`$ MeV. Curve labels are as in Fig. 3.

FIG. 5. Droplet radius $`R_d`$ vs temperature $`T`$ for various gap values $`\mathrm{\Delta }`$. Results are insensitive to $`n_B`$, so not all curves are labelled.

FIG. 6. Upper panel: Activation energy $`A`$ vs temperature $`T`$. Baryon number densities $`n_B`$ (fm-3) are (a) 0.5 (b) 1 (c) 1.5 (d) 2 (e) 2.5 (f) 3. Lower panel: The $`A=0`$ contour or the phase transition line in the $`Tn_B`$ plane. SC: superconducting phase.

FIG. 7. Free energies $`F_{QGP}`$ and $`F_{SC}`$ in the QGP and superconducting phases, respectively, vs temperature $`T`$.

FIG. 8. Three-dimensional surface plots of the activation energy $`A^{}`$ for $`\mathrm{\Delta }=1`$ MeV (upper part) and $`\mathrm{\Delta }=100`$ MeV (lower part), as a function of $`n_B`$ (fm-3) and $`T`$ (MeV).


text:
# 𝜀'/𝜀 in the standard model

## 1 Introduction

The $`CP`$-violating ratio $`\epsilon ^{}/\epsilon `$ (for a review see, e.g., ) is computed as

$$\text{Re}\epsilon ^{}/\epsilon =e^{i\varphi }\frac{G_F\omega }{2\left|ϵ\right|\text{Re}A_0}\text{Im}\lambda _t\left[\mathrm{\Pi }_0\frac{1}{\omega }\mathrm{\Pi }_2\right],$$
(1)
where, $`\lambda _t=V_{td}V_{ts}^{}`$ is a combination of Cabibbo–Kobayashi-Maskawa (CKM) matrix elements,

$`\mathrm{\Pi }_0`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{cos}\delta _0}}{\displaystyle \underset{i}{}}y_i\text{Re}Q_i_0(1\mathrm{\Omega }_{\eta +\eta ^{}}),`$ (2)
$`\mathrm{\Pi }_2`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{cos}\delta _2}}{\displaystyle \underset{i}{}}y_i\text{Re}Q_i_2,`$ (3)
and $`Q_i_I=2\pi ,I|Q_i|K`$. The $`\mathrm{\Delta }S=1`$ effective lagrangian $`_W`$ is given by

$$_W=\underset{i}{}C_i(\mu )Q_i(\mu ),$$
(4)
where

$$C_i(\mu )=\frac{G_F}{\sqrt{2}}V_{ud}V_{us}^{}\left[z_i(\mu )+\tau y_i(\mu )\right].$$
(5)
In (5), $`G_F`$ is the Fermi coupling, the functions $`z_i(\mu )`$ and $`y_i(\mu )`$ are the Wilson coefficients and $`\tau =V_{td}V_{ts}^{}/V_{ud}V_{us}^{}`$. According to the standard parameterization of the CKM matrix, in order to determine $`\epsilon ^{}/\epsilon `$ , we only need to consider the $`y_i(\mu )`$ components, which control the $`CP`$-violating part of the lagrangian. The coefficients $`y_i(\mu )`$, and $`z_i(\mu )`$ contain all the dependence of short-distance physics, and depend on the $`t,W,b,c`$ masses, the intrinsic QCD scale $`\mathrm{\Lambda }_{\mathrm{QCD}}`$, the $`\gamma _5`$-scheme used in the regularization and the renormalization scale $`\mu `$. They are known to the next-to-leading order (NLO) order in $`\alpha _s`$ and $`\alpha _e`$ .

The $`Q_i`$ in eq. (4) are the effective four-quark operators obtained in the standard model by integrating out the vector bosons and the heavy quarks $`t,b`$ and $`c`$. A convenient and by now standard basis includes the following ten operators:

$$\begin{array}{ccc}\hfill Q_1& =& \left(\overline{s}_\alpha u_\beta \right)_{\mathrm{V}\mathrm{A}}\left(\overline{u}_\beta d_\alpha \right)_{\mathrm{V}\mathrm{A}},\hfill \\ \hfill Q_2& =& \left(\overline{s}u\right)_{\mathrm{V}\mathrm{A}}\left(\overline{u}d\right)_{\mathrm{V}\mathrm{A}},\hfill \\ \hfill Q_{3,5}& =& \left(\overline{s}d\right)_{\mathrm{V}\mathrm{A}}_q\left(\overline{q}q\right)_{\mathrm{V}\mathrm{A}},\hfill \\ \hfill Q_{4,6}& =& \left(\overline{s}_\alpha d_\beta \right)_{\mathrm{V}\mathrm{A}}_q(\overline{q}_\beta q_\alpha )_{\mathrm{V}\mathrm{A}},\hfill \\ \hfill Q_{7,9}& =& \frac{3}{2}\left(\overline{s}d\right)_{\mathrm{V}\mathrm{A}}_q\widehat{e}_q\left(\overline{q}q\right)_{\mathrm{V}\pm \mathrm{A}},\hfill \\ \hfill Q_{8,10}& =& \frac{3}{2}\left(\overline{s}_\alpha d_\beta \right)_{\mathrm{V}\mathrm{A}}_q\widehat{e}_q(\overline{q}_\beta q_\alpha )_{\mathrm{V}\pm \mathrm{A}},\hfill \end{array}$$
(6)
where $`\alpha `$, $`\beta `$ denote color indices ($`\alpha ,\beta =1,\mathrm{},N_c`$) and $`\widehat{e}_q`$ are the quark charges ($`\widehat{e}_u=2/3`$, $`\widehat{e}_d=\widehat{e}_s=1/3`$). Color indices for the color singlet operators are omitted. The labels $`(V\pm A)`$ refer to the Dirac structure $`\gamma _\mu (1\pm \gamma _5)`$.

Notice the explicit presence of the final-state-interaction (FSI) phases $`\delta _I`$ in eqs. (2) and (3). Their presence is a consequence of writing the absolute values of the amplitudes in term of their dispersive parts.

Finally, $`\mathrm{\Omega }_{\eta +\eta ^{}}`$ is the isospin breaking (for $`m_um_d`$) contribution of the mixing of $`\pi `$ with $`\eta `$ and $`\eta ^{}`$.

Experimentally the ratio $`\epsilon ^{}/\epsilon `$ is extracted, by collecting $`K_L`$ and $`K_S`$ decays into pairs of $`\pi ^0`$ and $`\pi ^\pm `$, from the relation

$$\left|\frac{\eta _+}{\eta _{00}}\right|^21+6\text{Re}\frac{\epsilon ^{}}{\epsilon },$$
(7)
where

$$\eta _{00}\frac{\pi ^0\pi ^0|_W|K_L}{\pi ^0\pi ^0|_W|K_S}$$
(8)
and

$$\eta _+\frac{\pi ^+\pi ^{}|_W|K_L}{\pi ^+\pi ^{}|_W|K_S}.$$
(9)
The announcement last year of the preliminary result from the KTeV collaboration (FNAL)  based on data collected in 1996-97, and the present result from the NA48 collaboration (CERN), based on data collected in 1997 and 1998 , settle the long-standing issue of the presence of direct $`CP`$ violation in kaon decays. However, a clearcut determination the actual value of $`\epsilon ^{}`$ at the precision of a few parts in $`10^4`$ must wait for further statistics and scrutiny of the experimental systematics. By computing the average among the two 1992 experiments (NA31  and E731 ) and the KTeV and NA48 data we obtain

$$\text{Re}\epsilon /\epsilon ^{}=(1.9\pm 0.46)\times 10^3,$$
(10)
where the error has been inflated according to the Particle Data Group procedure ($`\sigma \sigma \times \sqrt{\chi ^2/3}`$), to be used when averaging over experimental data—in our case four sets—with substantially different central values.

The value in eq. (10) can be considered the current experimental result. Such a result will be tested within the next year by the full data analysis from KTeV and NA48 and (hopefully) the first data from KLOE at DA$`\mathrm{\Phi }`$NE (Frascati); at that time, the experimental uncertainty will be reduced to a few parts in $`10^4`$.

## 2 Hadronic matrix elements

If there is no sizable cancellation between the relevant effective operators, the order of magnitude of $`\epsilon ^{}/\epsilon `$ is bound to be of the order of $`10^3`$. A simple argument for this is presented in . The problem is that any cancellation, or the lack thereof, among the operators heavily depends on the size of the hadronic matrix elements and there is no available estimate of them that is free of hard-to-control assumptions. To estimate the hadronic matrix elements in a systematic manner without having first to solve QCD (not even by lattice simulation) we needed a model that is simple enough to understand its dynamics and, at the same time, not too simple so as to still include the relevant physics. We chose the chiral quark model ($`\chi `$QM)  in which all coefficients of the relevant chiral lagrangian are parameterized in terms of just three parameters: the quark and gluon condensates, and the quark constituent mass. The model makes possible a complete estimate of all matrix elements, it includes non-factorizable effects, chiral corrections and final-state interaction, all of which we thought to be relevant.

The $`\chi `$QM  is an effective quark model of QCD which can be derived in the framework of the extended Nambu-Jona-Lasinio model of chiral symmetry breaking (for a review, see, e.g., ). In the $`\chi `$QM an effective interaction between the $`u,d,s`$ quarks and the meson octet is introduced via the term

$$_{\chi \text{QM}}=M\left(\overline{q}_R\mathrm{\Sigma }q_L+\overline{q}_L\mathrm{\Sigma }^{}q_R\right),$$
(11)
which is added to an effective low-energy QCD lagrangian whose dynamical degrees of freedom are the $`u,d,s`$ quarks propagating in a soft gluon background. The matrix $`\mathrm{\Sigma }`$ in (11) is the same as that used in chiral perturbation theory and it contains the pseudo-scalar meson multiplet. The quantity $`M`$ is interpreted as the constituent quark mass in mesons (current quark masses are also included in the effective lagrangian).

The hadronic matrix elements are matched with the NLO Wilson coefficients at the scale $`\mathrm{\Lambda }_\chi 0.8`$ ($`m_\rho `$) as the best compromise between the range of validity of chiral perturbation and that of strong coupling expansion.

## 3 The fit of the $`\mathrm{\Delta }I=1/2`$ rule

In order to assign the values of the model-dependent parameters $`M`$, $`\overline{q}q`$ and $`\alpha _sGG/\pi `$ , we consider the $`CP`$-conserving amplitudes in the $`\mathrm{\Delta }I=1/2`$ selection rule of $`K\pi \pi `$ decays. In practice, we compute the amplitudes

$`A_0`$ $`=`$ $`{\displaystyle \frac{G_F}{\sqrt{2}}}V_{ud}V_{us}^{}{\displaystyle \frac{1}{\mathrm{cos}\delta _0}}{\displaystyle \underset{i}{}}z_i(\mu )\text{Re}Q_i(\mu )_0,`$
$`A_2`$ $`=`$ $`{\displaystyle \frac{G_F}{\sqrt{2}}}V_{ud}V_{us}^{}{\displaystyle \frac{1}{\mathrm{cos}\delta _2}}{\displaystyle \underset{i}{}}z_i(\mu )\text{Re}Q_i(\mu )_2`$ (12)
$`+\omega A_0\mathrm{\Omega }_{\eta +\eta ^{}},`$
within the $`\chi `$QM approach and vary the parameters in order to reproduce their experimental values.

This procedure combines a model for low-energy QCD—which allow us to compute all hadronic matrix elements in terms of a few basic parameters—with the phenomenological determination of such parameters. In this way, some shortcomings of such a naive model (in particular, the matching between long- and short-distance components) are absorbed in the phenomenological fit. As a check, we eventually verify the stability of the computed observables against renormalization scale and scheme dependence.

The fit of the $`CP`$-conserving involves the determination of the FSI phases. The absorptive components of the hadronic matrix elements appear when chiral loops are included. In our approach the direct determination of the rescattering phases gives at $`O(p^4)`$ $`\delta _020^0`$ and $`\delta _212^0`$. Although these results show features which are in qualitative agreement with the phases extracted from pion-nucleon scattering the deviation from the experimental data is sizeable, especially in the $`I=0`$ component. On the other hand, at $`O(p^4)`$ the absorptive parts of the amplitudes are determined only at $`O(p^2)`$ and disagreement with the measured phases should be expected. As a matter of fact, the authors of ref.  find that at $`O(p^6)`$ the absorptive part of the hadronic matrix elements are substantially modified to give values of the rescattering phases quite close to those extracted from pion-nucleon scattering. At the same time the $`O(p^6)`$ corrections to the dispersive part of the hadronic matrix elements are very small. This result corroborates our ansatz  of trusting the real parts of the $`O(p^4)`$ matrix elements while inputting the experimental values of the rescattering phases in all parts of our analysis, which amounts to taking $`\mathrm{cos}\delta _00.8`$ and $`\mathrm{cos}\delta _21`$.

Hadronic matrix elements in the $`\chi `$QM depend on the $`\gamma _5`$-scheme utilized . Their dependence partially cancels that of the short-distance NLO Wilson coefficients. Because this compensation is only numerical, and not analytical, we take it as part of our phenomenological approach. A formal $`\gamma _5`$-scheme matching can only come from a model more complete than the $`\chi `$QM. Nevertheless, the result, as shown in Fig. 2, is rather convincing.

By taking

$$\mathrm{\Lambda }_{\mathrm{QCD}}^{(4)}=340\pm 40\text{MeV}$$

and fitting at the scale $`\mu =0.8`$ GeV the amplitudes in eqs. (3)–(12) to their experimental values, allowing for a $`\pm `$ 30% systematic uncertainty, we find (see Fig. 1)

$`M`$ $`=`$ $`195_{15}^{+25}\text{MeV},`$
$`\alpha _sGG/\pi `$ $`=`$ $`\left(330\pm 5\text{MeV}\right)^4,`$
$`\overline{q}q`$ $`=`$ $`\left(235\pm 25\text{MeV}\right)^3.`$ (13)
in the HV-scheme, and

$`M`$ $`=`$ $`195_{10}^{+15}\text{MeV},`$
$`\alpha _sGG/\pi `$ $`=`$ $`\left(333_6^{+7}\text{MeV}\right)^4,`$
$`\overline{q}q`$ $`=`$ $`\left(245\pm 15\text{MeV}\right)^3.`$ (14)
in the NDR-scheme.

As shown by the light (NDR) and dark (HV) curves in Fig. 2, the $`\gamma _5`$-scheme dependence is controlled by the value of $`M`$, the range of which is fixed thereby. The $`\gamma _5`$ scheme dependence of both amplitudes is minimized for $`M190200`$ MeV. The good $`\gamma _5`$-scheme stability is also enjoyed by $`\epsilon ^{}/\epsilon `$ and $`\widehat{B}_K`$.

The fit of the amplitude $`A_0`$ and $`A_2`$ is obtained for values of the quark and gluon condensates which are in agreement with those found in other approaches, i.e. QCD sum rules and lattice, although the relation between the gluon condensate of QCD sum rules and lattice and that of the $`\chi `$QM is far from obvious. The value of the constituent quark mass $`M`$ is in good agreement with that found by fitting radiative kaon decays .

In Fig. 3 we present the anatomy of the relevant operator contributions to the $`CP`$ conserving amplitudes. It is worth noticing that, because of the NLO enhancement of the $`I=0`$ matrix elements (mainly due to the chiral loops), the gluon penguin contribution to $`A_0`$ amounts to about 20% of the amplitude.

Turning now to the $`\mathrm{\Delta }S=2`$ Lagrangian, by using the input values found by fitting the $`\mathrm{\Delta }I=1/2`$ rule we obtain in both $`\gamma _5`$-schemes

$$\widehat{B}_K=1.1\pm 0.2.$$
(15)
The result (15) includes chiral corrections up to $`O(p^4)`$ and it agrees with what we found in .

Notice that no estimate of $`\epsilon ^{}/\epsilon `$ can be considered complete unless it also gives a value for $`\widehat{B}_K`$. The case of the $`\chi `$QM, for instance, is telling insofar as the enhancement of $`B_6`$ is partially compensated for by a large $`\widehat{B}_K`$ (and accordingly a smaller $`\text{Im}\lambda _t`$).

### 3.1 The bag parameters $`B_i`$

The bag parameters $`B_i`$, defined as

$$B_iQ_i^{\mathrm{model}}/Q_i^{\mathrm{VSA}},$$
(16)
have become a standard way of displaying the values of the hadronic matrix elements in order to compare them among various approaches. However they must be used with care because of their dependence on the renormalization scheme and scale, as well as on the choice of the VSA parameters.

They are given in the $`\chi `$QM in table 1 in the HV and NDR schemes, at $`\mu =0.8`$ GeV, for the central value of $`\mathrm{\Lambda }_{\mathrm{QCD}}^{(4)}`$. The uncertainty in the matrix elements of the penguin operators $`Q_{58}`$ arises from the variation of $`\overline{q}q`$. This affects mostly the $`B_{5,6}`$ parameters because of the leading linear dependence on $`\overline{q}q`$ of the $`Q_{5,6}`$ matrix elements in the $`\chi `$QM, contrasted to the quadratic dependence of the corresponding VSA matrix elements. Accordingly, $`B_{5,6}`$ scale as $`\overline{q}q^1`$, or via PCAC as $`m_q`$, and therefore are sensitive to the value chosen for these parameters. For this reason, we have reported the corresponding values of $`B_{5,6}`$ when the quark condensate in the VSA is fixed to its PCAC value. It should however be stressed that such a dependence is not physical and is introduced by the arbitrary normalization on the VSA result. The estimate of $`\epsilon ^{}`$ is therefore almost independent of $`m_q`$, which only enters the NLO corrections and the determination of $`\widehat{B}_K`$.

The enhancement of the $`Q_{5,6}`$ matrix elements with respect to the VSA values—the conventional normalization of the VSA matrix elements corresponds to taking $`\overline{q}q`$$`(0.8\mathrm{GeV})(220\mathrm{MeV})^3`$—is mainly due to the NLO chiral loop contributions. Such an enhancement, due to final state interactions, has been found in $`1/N_c`$ analyses beyond LO , as well as old and recent dispersive studies .

## 4 Bounds on $`\text{Im}\lambda _t`$

The updated measurements for the CKM elements $`|V_{ub}/V_{cb}|`$ implies a change in the determination of the Wolfenstein parameter $`\eta `$ that enter in $`\text{Im}\lambda _t`$. This is of particular relevance because it affects proportionally the value of $`\epsilon ^{}/\epsilon `$.

The allowed values for $`\text{Im}\lambda _t\eta |V_{us}||V_{cb}|^2`$ are found by imposing the experimental constraints for $`\epsilon `$, $`|V_{ub}/V_{cb}|`$, $`\mathrm{\Delta }m_d`$ and $`\mathrm{\Delta }m_s`$. By using the method of Parodi, Roudeau and Stocchi , who have run their program starting from our inputs listed in , it is found that

$$\text{Im}\lambda _t=(1.14\pm 0.11)\times 10^4,$$
(17)
where the error is determined by the Gaussian distribution in Fig. 4.

## 5 Estimating $`\epsilon ^{}/\epsilon `$

The value of $`\epsilon ^{}`$ computed by taking all input parameters at their central values is shown in Fig. 5. The figure shows the contribution to $`\epsilon ^{}`$ of the various operators in two $`\gamma _5`$ renormalization schemes at $`\mu =0.8`$ GeV and $`1.0`$ GeV. The advantage of such a histogram is that, contrary to the $`B_i`$, the size of the individual contributions does not depend on some conventional normalization.

As the histogram makes it clear, the gluon-penguin operator $`Q_6`$ dominates and there is very little cancellation with the electroweak penguin operator $`Q_8`$. The dominance of the $`I=0`$ components in the $`\chi `$QM originates from the $`O(p^4)`$ chiral corrections, the detailed size of which is determined by the fit of the $`\mathrm{\Delta }I=1/2`$ rule. It is a nice feature of the approach that the renormalization scheme stability imposed on the $`CP`$ conserving amplitudes is numerically preserved in $`\epsilon ^{}/\epsilon `$. The comparison of the two figures shows also the remarkable renormalization scale stability of the central value once the perturbative running of the quark masses and the quark condensate is taken into account.

In what follows, the model-dependent parameters $`M`$, $`\alpha _sGG/\pi `$ and $`\overline{q}q`$ are uniformly varied in their given ranges (flat scanning), while the others are sampled according to their normal distributions. Values of $`\epsilon ^{}/\epsilon `$ found in the HV and NDR schemes are included with equal weight.

For a given set, a distribution is obtained by collecting the values of $`\epsilon ^{}/\epsilon `$ in bins of a given range. This is shown in Fig. 6 for a particular choice of bins. Because the skewness of the distribution is less than one, the mean and the standard deviation are a good estimate of the central value and the dispersion of values around it.

This statistical analysis yields

$$\epsilon ^{}/\epsilon =(2.2\pm 0.8)\times 10^3.$$
(18)
A more conservative estimate of the uncertainties is obtained via the flat scanning of the input parameters, which gives

$$0.9\times 10^3<\text{Re}\epsilon ^{}/\epsilon <4.8\times 10^3.$$
(19)
In both estimates a theoretical systematic error of $`\pm 30\%`$ is included in the fit of the $`CP`$ conserving amplitudes $`A_0`$ and $`A_2`$.

The stability of the numerical outcomes is only marginally affected by shifts in the value of $`\mathrm{\Omega }_{\eta +\eta ^{}}`$. Any effective variation of $`\mathrm{\Omega }_{\eta +\eta ^{}}`$ is anti-correlated to the value of $`\alpha _sGG/\pi `$ obtained in the fit of $`A_2`$. We have verified that this affects the calculation of $`\widehat{B}_K`$ and the consequent determination of $`\text{Im}\lambda _t`$ in such a way to compensate numerically in $`\epsilon ^{}/\epsilon `$ the change of $`\mathrm{\Omega }_{\eta +\eta ^{}}`$. Waiting for a confident assessment of the NLO isospin violating effects in the $`K\pi \pi `$ amplitudes, we have used for $`\mathrm{\Omega }_{\eta +\eta ^{}}`$ the ‘LO’ value $`0.25\pm 0.10`$.

The weak dependence on some poorly-known parameters is a welcome outcome of the correlation among hadronic matrix elements enforced in our semi-phenomenological approach by the fit of the $`\mathrm{\Delta }I=1/2`$ rule.

While the $`\chi QM`$ approach to the hadronic matrix elements relevant in the computation of $`\epsilon ^{}/\epsilon `$ has many advantages over other techniques and has proved its value in the prediction of what has been then found in the experiments, it has a severe short-coming insofar as the matching scale has to be kept low, around 1 GeV and therefore the Wilson coefficients have to be run down at scales that are at the limit of the applicability of the renormalization-group equations. Moreover, the matching itself suffers of ambiguities that have not be completely solved. For these reasons we have insisted all along that the approach is semi-phenomenological and that the above shortcomings are to be absorbed in the values of the input parameters on which the fit to the $`CP`$ conserving amplitudes is based. Because of its simplicity, the $`\chi `$QM is clearly not the final word and it can now been abandoned—as a ladder used to climb a wall after we are on the other side—as we work for better estimates, in particular, those from the lattice simulations.

## 6 Other estimates

Figure 7 summarizes the present status of theory versus experiment. In addition to the current $`\chi `$QM calculation (and an independent estimate similarly based on the $`\chi `$QM), there are eight estimates of $`\epsilon ^{}/\epsilon `$ published in the last year. Trieste’s, München’s and Roma’s ranges are updates of their respective older estimates, while the others are altogether new.

The estimates reported are discussed in . The effects first pointed out by the $`\chi `$QM have been taken up and refined by some of the other approaches, in particular: non-factorizable corrections  and FSI . Figure 7 makes clear that, except for estimates heavily based on VSA for the crucial operator $`Q_6`$, all other estimates are in agreement with the experiments.

## Acknowledgments

It is a pleasure to thank my collaborators S. Bertolini and J. O. Eeg and my former students V. Antonelli and E. I. Lashin for the work done together.