text:
# 1 Probability density distributions for factors 𝐵̃₁, 𝐵̃₄ and 𝐵̃₅ with 𝐵̃₈=1.

DIRECT CP VIOLATION IN NONLEPTONIC KAON DECAYS
BY AN EFFECTIVE CHIRAL LAGRANGIAN APPROACH AT $`O(p^6)`$111Talk at the XV Int. Seminar on High Energy Physics Problems ”Relativistic Nuclear Physics & Quantum Chromodynamics”, Dubna, September 25–29, 2000.

A.A. Bel’kov1†, G. Bohm2, A.V. Lanyov1 and A.A. Moshkin1

(1) Particle Physics Laboratory, Joint Institute for Nuclear Research,
141980 Dubna, Moscow region, Russia
(2) DESY–Zeuthen, Platanenallee 6, D–15735 Zeuthen, Germany
$``$ E-mail: belkov@cv.jinr.dubna.su

Abstract

A self-consistent analysis of $`K2\pi `$ and $`K3\pi `$ decays within a unique framework of chiral dynamics applied to the QCD-corrected weak nonleptonic quark lagrangian has been performed. The results on $`K2\pi `$ amplitudes at $`O(p^6)`$, including the value for $`\epsilon ^{^{}}/\epsilon `$, are compared with the experiment to fix phenomenological B-factors for mesonic matrix elements of nonpenguin and penguin four-quark operators. The dependence of B-factors on different theoretical uncertainties and experimental errors of various input parameters is investigated. Finally, we present our estimates at $`O(p^6)`$ for the CP-asymmetry of linear slope parameters in the $`K^\pm 3\pi `$ Dalitz plot.

Keywords: direct CP violation, nonleptonic decays, K-mesons, chiral model

The starting point for most calculations of nonleptonic kaon decays is an effective weak lagrangian of the form

$$_w^q\left(|\mathrm{\Delta }S|=1\right)=\sqrt{2}G_FV_{ud}V_{us}^{}\underset{i}{}\stackrel{~}{C}_i𝒪_i,$$
(1)
which can be derived with the help of the Wilson operator product expansion (OPE) from elementary quark processes, with additional gluon exchanges. In the framework of perturbative QCD the coefficients $`\stackrel{~}{C}_i`$ are to be understood as scale and renormalization scheme dependent functions.

In (1), $`𝒪_i`$ are the four-quark operators, defined either by combinations of products of quark currents ($`i=1,2,3,4`$, non-penguin diagrams) or, in case of gluonic ($`i=5,6`$) and electro-weak ($`i=7,8`$) penguin operators, by products of quark densities. The operators $`𝒪_i`$ with $`i=1,2,3,5,6`$ describe weak transitions with isospin change $`\mathrm{\Delta }I=1/2`$ while the operator $`𝒪_4`$ corresponds to $`\mathrm{\Delta }I=3/2`$ transition and operators $`𝒪_{7,8}`$ – to mixture of $`\mathrm{\Delta }I=1/2`$ and $`\mathrm{\Delta }I=3/2`$ amplitudes.

The general scheme of the meson matrix element calculation by using of the weak lagrangian (1) is based on the quark bosonization approach . The bosonization procedure establishes the correspondence between quark and meson currents (densities) and product of currents (densities). Finally, it leads to the effective lagrangian for nonleptonic kaon decays in terms of bosonized (meson) currents and densities:

$$\overline{q}\gamma _\mu \frac{1}{4}(1\gamma ^5)\lambda ^aqJ_{L/R\mu }^{a(mes)},\overline{q}\frac{1}{4}(1\gamma ^5)\lambda ^aqJ_{L/R}^{a(mes)}.$$

The meson currents/densities $`J_{L/R\mu }^a`$ and $`J_{L/R}^a`$ are obtained from the quark determinant by variation over additional external sources associated with the corresponding quark currents and densities . From the momentum expansion of the quark determinant to $`O(p^{2n})`$ one can derive the strong lagrangian for mesons $`_{eff}`$ of the same order and the corresponding currents and densities $`J_{L/R\mu }^a`$ and $`J_{L/R}^a`$ to the order $`O(p^{2n1})`$ and $`O(p^{2n2})`$, respectively. Thus, the bosonization approach gives us the correspondence between power counting for the momentum expansion of the effective chiral lagrangian of strong meson interactions,

$$_s^{(mes)}=_s^{(p^2)}+_s^{(p^4)}+_s^{(p^6)}+...,$$

and power counting for the meson currents and densities:

$`_s^{(p^n)}J_\mu ^{(p^{n1})}\text{(currents)};_s^{(p^n)}J^{(p^{n2})}\text{(densities)}.`$
Some interesting observations on the difference of the momentum behavior of penguin and non-penguin operators can be drawn from power-counting arguments. The leading contributions to the vector currents and scalar densities are of $`O(p^1)`$ and $`O(p^0)`$, respectively. Since in our approach the non-penguin operators are constructed out of the products of currents $`J_{L\mu }^a`$, while the penguin operators are products of densities $`J_L^a`$, the lowest-order contributions of non-penguin and penguin operators are of $`O(p^2)`$ and $`O(p^0)`$, respectively. However, due to the well-known cancelation of the contribution of the gluonic penguin operator $`𝒪_5`$ at the lowest order , the leading gluonic penguin as well as non-penguin contributions start from $`O(p^2)`$ 222There is no cancellation of the contribution of the electromagnetic penguin operator $`𝒪_8`$ at the lowest order and the leading contributions start in this case from $`O(p^0)`$. Consequently, in order to derive the currents which contribute to the non-penguin transition operators at the leading order, it is sufficient to use the terms of the quark determinant to $`O(p^2)`$ only. At the same time the terms of the quark determinant to $`O(p^4)`$ have to be kept for calculating the penguin contribution at $`O(p^2)`$, since it arises from the combination of densities, which are of $`O(p^0)`$ and $`O(p^2)`$, respectively. In this subtle way the difference in momentum behavior is revealed between matrix elements for these two types of weak transition operators; it manifests itself more drastically in higher-order lagrangians and currents.

This fact makes penguins especially sensitive to higher order effects. In particular, the difference in the momentum power counting behavior between penguin and non-penguin contributions to the isotopic amplitudes of $`K3\pi `$ decays, which appears in higher orders of chiral theory, leads to the dynamical enhancement of the charge asymmetry of the Dalitz-plot linear slope parameter .

In our approach the Wilson coefficients $`\stackrel{~}{C}_i`$ in the effective weak lagrangian (1) are treated as the phenomenological parameters which should be fixed from the experiment. They are related with the Wilson coefficients $`\stackrel{~}{C}_i^{QCD}(\mu )`$, calculated in perturbative QCD , via the $`\stackrel{~}{B}_i`$-factors:

$$\stackrel{~}{C}_i^{ph}=\stackrel{~}{C}_i^{QCD}(\mu )\stackrel{~}{B}_i(\mu ).$$

The coefficients $`\stackrel{~}{C}_i^{QCD}(\mu )`$ contain a small imaginary parts and can be presented in the form of a sum of $`z`$ and $`y`$ components

$$\stackrel{~}{C}_i^{QCD}(\mu )=\stackrel{~}{C}_i^{(z)}(\mu )+\tau \stackrel{~}{C}_i^{(y)}(\mu ),\tau =\frac{V_{td}V_{ts}^{}}{V_{ud}V_{us}^{}}.$$
(2)
Respectively, the amplitudes of nonleptonic kaon decays also contain $`z`$\- and $`y`$\- components,

$$𝒜=𝒜^{(z)}+\tau 𝒜^{(y)},$$
(3)
which can be expressed in terms of $`z`$\- and $`y`$-components of Wilson coefficients. The dominating contributions $`𝒜^{(i)}`$ of the four-quark operators $`𝒪_i`$ and $`\stackrel{~}{B}_i`$-factors may be written as

$`𝒜^{(z,y)}`$ $`=`$ $`\left[\stackrel{~}{C}_1^{(z,y)}(\mu )+\stackrel{~}{C}_2^{(z,y)}(\mu )+\stackrel{~}{C}_3^{(z,y)}(\mu )\right]\stackrel{~}{B}_1(\mu )𝒜^{(1)}`$
$`+\stackrel{~}{C}_4^{(z,y)}(\mu )\stackrel{~}{B}_4(\mu )𝒜^{(4)}+\stackrel{~}{C}_5^{(z,y)}(\mu )\stackrel{~}{B}_5(\mu )𝒜^{(5)}+\stackrel{~}{C}_8^{(z,y)}(\mu )\stackrel{~}{B}_8(\mu )𝒜^{(8)}.`$
The observable effects of direct CP-violation in the nonleptonic kaon decays are caused by the $`y`$-components of their amplitudes (3). In particular, the ratio $`\epsilon ^{^{}}/\epsilon `$ can be expressed as

$$\text{Re}\left(\frac{\epsilon ^{^{}}}{\epsilon }\right)=\text{Im}\lambda _t(P_0P_2),P_I=\frac{\omega }{\sqrt{2}|\epsilon ||V_{ud}||V_{us}|}\frac{𝒜_I^{(y)}}{𝒜_I^{(z)}},$$

where $`\text{Im}\lambda _t=\text{Im}V_{ts}^{}V_{td}=|V_{ub}||V_{cb}|\text{sin}\delta `$; $`\omega =𝒜_2^{(z)}/𝒜_0^{(z)}`$, and $`𝒜_I`$ are isotopic amplitudes of $`K2\pi `$ decay corresponding to the $`\pi \pi `$ final states with the isospin $`I=0,2`$.

There are the following theoretical uncertainties which appear both from short distance (Wilson coefficients) and long distance (effective chiral lagrangians and $`\stackrel{~}{B}_i`$-factors) contributions to the kaon decay amplitude:

* Dependence on the parameter $`\text{Im}\tau \text{Im}\lambda _t`$ arising from the imaginary part of Wilson coefficient (2).
* Regularization scheme dependence which arises when calculating Wilson coefficients beyond the leading order (LO) of QCD in the next-to-leading orders: naive dimensional regularization (NDR), t’Hooft-Veltman regularization (HV).
* Dependence of Wilson coefficients on choice of the renormalization point $`\mu `$ and QCD scale $`\mathrm{\Lambda }_{\overline{MS}}^{(4)}=(325\pm 110)`$ MeV.
* Dependence of meson matrix elements on the structure constants $`L_2`$, $`L_3`$, $`L_4`$, $`L_5`$, $`L_8`$ of the general form of the effective chiral lagrangian introduced at $`O(p^4)`$ by Gasser and Leutwyler .
* Dependence of meson matrix elements on the structure constants of the effective chiral lagrangian at $`O(p^6)`$ .
* Factors $`\stackrel{~}{B}_i`$ $`(i=1,4,5,8)`$ for dominating contributions of four-quark operators $`𝒪_i`$ to the meson matrix element (S0.Ex6).
* Dependence on choice of the regularization scheme to fix the UV divergences resulting from meson loops.

In the present paper we have combined a new systematic calculation of mesonic matrix elements for nonleptonic kaon decays from the effective chiral lagrangian approach with Wilson coefficients $`\stackrel{~}{C}_i^{QCD}(\mu )`$, derived by the Munich group  for $`\mu =1`$ GeV and $`m_t=170`$ GeV. For the parameter $`\text{Im}\lambda _t`$ we have used the result obtained in : $`\text{Im}\lambda _t=(1.33\pm 0.14)10^4`$. We performed a complete calculation of $`K2\pi `$ and $`K3\pi `$ amplitudes at $`O(p^6)`$ including the tree level, one- and two-loop diagrams. For the structure coefficients $`L_i`$ of the effective chiral lagrangian at $`O(p^4)`$ we used the values fixed in  from the phenomenological analysis of low-energy meson processes. The structure coefficients of the effective chiral lagrangian at $`O(p^6)`$ have been fixed theoretically from the modulus of the logarithm of the quark determinant of the NJL-type model (see  for more details). The superpropagator regularization has been applied to fix UV divergences in meson loops. The isotopic symmetry breaking ($`\pi ^0`$-$`\eta `$-$`\eta ^{^{}}`$ mixing) was taken into account at the tree level. In  one can find more technical details of the calculation of $`K2\pi `$ amplitudes.

In our phenomenological analysis the results on the $`K2\pi `$ amplitudes, including the value for $`\epsilon ^{^{}}/\epsilon `$, are compared with the experiment to fix the phenomenological $`\stackrel{~}{B}_i`$-factors for the mesonic matrix elements of nonpenguin and penguin four-quark operators. As experimental input we used the experimental values of the isotopic amplitudes $`𝒜_{0,2}^{(exp)}`$ fixed from the widths of $`K2\pi `$ decays and the world average value $`\text{Re}\epsilon ^{^{}}/\epsilon =(19.3\pm 3.8)\times 10^4`$ which includes both old results of NA31  and E731  experiments and recent results from KTeV  and NA48 . The output parameters of the performed $`K2\pi `$ analysis are the factors $`\stackrel{~}{B}_1`$, $`\stackrel{~}{B}_4`$ and $`\stackrel{~}{B}_5`$ for a fixed value of $`\stackrel{~}{B}_8`$.

The dependence of $`\stackrel{~}{B}_i`$-factors on different theoretical uncertainties and experimental errors of various input parameters is investigated by applying the “Gaussian” method. Using Wilson coefficients derived in  in various regularization schemes (LO, NDR, HV) for different values of the QCD scale $`\mathrm{\Lambda }_{\overline{MS}}^{(4)}`$ ,we calculated the probability density distributions for $`\stackrel{~}{B}_i`$-factors obtained by using Gaussian distribution for all input parameters with their errors. As an example, the probability densities for the parameters $`\stackrel{~}{B}_1`$, $`\stackrel{~}{B}_4`$, $`\stackrel{~}{B}_5`$ calculated with $`\stackrel{~}{B}_8=1`$ and $`\mathrm{\Lambda }_{\overline{MS}}^{(4)}=325`$ MeV are shown in Figure 1. We have shown the necessity for a rather large gluonic penguin contribution to describe the recently confirmed large experimental $`\epsilon ^{^{}}`$ value (the factor $`\stackrel{~}{B}_5`$ is found well above 1). Figure 2 shows the correlations between $`\stackrel{~}{B}_5`$ and $`\stackrel{~}{B}_8`$ calculated for central values of all input parameters. From this figure one can see that even for $`\stackrel{~}{B}_8=0`$ values of $`\stackrel{~}{B}_5>2`$ are necessary to explain the large value of $`\epsilon ^{^{}}/\epsilon `$. It should be emphasized, that for even larger values of $`\stackrel{~}{B}_5`$, the contributions of nonpenguin operators to the $`\mathrm{\Delta }I=1/2`$ amplitude are still dominating (see Figure 3). The large $`\stackrel{~}{B}_1`$ and $`\stackrel{~}{B}_5`$ values may be a hint that the long-distance contributions, especially to $`\mathrm{\Delta }I=1/2`$ amplitudes, are still not completely understood. An analogous conclusion has been drawn in , where possible effects from physics beyond the Standard Model are also discussed.

Finally, we present our predictions for the CP-asymmetry of linear slope parameters in the $`K^\pm 3\pi `$ Dalitz plot. These predictions are based on a new calculation of $`K3\pi `$ amplitudes at $`O(p^6)`$ within the same effective lagrangian approach. The obtained $`K3\pi `$ amplitudes include the same theoretical uncertainties as in case of the $`K2\pi `$ analysis. The values of $`\stackrel{~}{B}_1`$, $`\stackrel{~}{B}_4`$, $`\stackrel{~}{B}_5`$ fixed from the $`K2\pi `$ analysis are used as phenomenological input to the $`K3\pi `$ estimates which have been performed in a self-consistent way by the Gaussian method.

The linear slope parameter $`g`$ of the Dalitz plot for $`K3\pi `$ decays is defined through the expansion of the decay probability over the kinematic variables

$$|T(K3\pi )|^2\mathrm{\hspace{0.17em}\hspace{0.17em}1}+gY+...$$

where $`Y`$ is a Dalitz variable:

$$Y=\frac{s_3s_0}{m_\pi ^2},s_3=(p_Kp_{\pi _3})^2,s_0=\frac{m_K^2}{3}+m_\pi ^2,$$

and the notation $`\pi _3`$ belongs to the non-equivalent pion in the decays $`K^\pm \pi ^\pm \pi ^\pm \pi ^{}`$ and $`K^\pm \pi ^0\pi ^0\pi ^\pm `$. Direct CP-violation leads to a charge asymmetry of the linear slope parameter,

$$\mathrm{\Delta }g(K^\pm 3\pi )=\frac{g(K^+3\pi )g(K^{}3\pi )}{g(K^+3\pi )+g(K^{}3\pi )}.$$

In Figure 4 we show the probability density distributions for $`K^\pm \pi ^\pm \pi ^\pm \pi ^{}`$ and $`K^\pm \pi ^0\pi ^0\pi ^\pm `$ decays calculated with $`\stackrel{~}{B}_8=1`$ and $`\mathrm{\Lambda }_{\overline{MS}}^{(4)}=325`$ MeV. Upper and low bounds for $`\mathrm{\Delta }g_{++}`$ and $`\mathrm{\Delta }g_{00+}`$ for different values of $`\mathrm{\Lambda }_{\overline{MS}}^{(4)}`$ in LO, NDR and HV regularization schemes ($`\stackrel{~}{B}_8=1`$) obtained by the Gaussian method are shown in table 1. The limits without brackets correspond to the confidence level of 68% while the limits in brackets – to the confidence level of 95%. Summarizing these results, we have obtained the following upper and lower bounds for the charge symmetries of the linear slope parameter:

$$2.1<\mathrm{\Delta }g_{++}10^4<\mathrm{\hspace{0.17em}\hspace{0.17em}6.2},1.4<\mathrm{\Delta }g_{00+}10^4<\mathrm{\hspace{0.17em}\hspace{0.17em}3.5}\text{with\hspace{0.17em}\hspace{0.17em}CL=68\%};$$

$$1.4<\mathrm{\Delta }g_{++}10^4<\mathrm{\hspace{0.17em}10.4},1.0<\mathrm{\Delta }g_{00+}10^4<\mathrm{\hspace{0.17em}5.1}\text{with\hspace{0.17em}\hspace{0.17em}CL=95\%}.$$

Table 1. Upper and low bounds for $`\mathrm{\Delta }g_{++}`$ and $`\mathrm{\Delta }g_{00+}`$ (in units $`10^4`$).

| $`\mathrm{\Delta }g`$ | $`\mathrm{\Lambda }_{\overline{MS}}^{\left(4\right)}`$ | LO | | NDR | | HV | |
| --- | --- | --- | --- | --- | --- | --- | --- |
|  | (MeV) | min | max | min | max | min | max |
|  | 215 | 2.5 | 5.9 | 2.1 | 4.7 | 2.6 | 6.0 |
|  |  | ( 1.8 | 10.1) | ( 1.5 | 7.3 ) | ( 1.9 | 10.2 ) |
| $`\mathrm{\Delta }g_{++}`$ | 325 | 2.6 | 6.0 | 2.2 | 4.7 | 2.6 | 5.8 |
|  |  | ( 1.9 | 10.2) | ( 1.5 | 7.2 ) | ( 1.8 | 9.7 ) |
|  | 435 | 2.7 | 6.2 | 2.2 | 4.7 | 2.6 | 5.8 |
|  |  | ( 1.9 | 10.4) | ( 1.4 | 7.0) | ( 1.8 | 9.5 ) |
|  | 215 | 1.5 | 3.3 | 1.4 | 2.9 | 1.6 | 3.4 |
|  |  | ( 1.1 | 4.8 ) | ( 1.0 | 4.2 ) | ( 1.2 | 5.1 ) |
| $`\mathrm{\Delta }g_{00+}`$ | 325 | 1.6 | 3.4 | 1.4 | 3.0 | 1.6 | 3.3 |
|  |  | ( 1.2 | 5.0 ) | ( 1.0 | 4.3 ) | ( 1.2 | 4.9 ) |
|  | 435 | 1.7 | 3.5 | 1.5 | 3.1 | 1.7 | 3.4 |
|  |  | ( 1.2 | 5.1 ) | ( 1.0 | 4.3 ) | ( 1.2 | 5.0 ) |

With some experimental updates and theoretical refinements our new estimates confirm the dynamical enhancement mechanism for the charge asymmetry $`\mathrm{\Delta }g`$ by higher order contributions in the effective chiral lagrangian approach which was first observed in . The predicted slope parameter asymmetry, although small, may be in reach of current high statistics experiments . Already with lower statistics, new measurements of quadratic slope parameters of $`K3\pi `$ decays, including the neutral channels, would lead to the improved theoretical understanding of the nonperturbative part of nonleptonic kaon decay dynamics.


format: text

text:
# Edinburgh-2000/20 Leptonic and semi-leptonic B decays

## 1 INTRODUCTION

The determination of decay constants and form factors arising from the decay of B and D mesons is of great importance to phenomenology, in particular in constraining the Cabibbo-Kobayashi-Maskawa (CKM) matrix. These quantities can be calculated on the lattice.

The action used in this calculation is the $`𝒪(a)`$ non-perturbatively improved action , where $`c_{SW}`$ and the coefficients needed to improve the currents, are determined non-perturbatively where possible. The calculation was done at two values of the coupling, $`\beta =6.0`$ and $`\beta =6.2`$. The simulation parameters are shown in Table 1. We simulate with three light quark masses around strange and four propagating heavy quarks around charm. Heavy Quark Symmetry (HQS) is used to motivate the form of the heavy extrapolation to the $`B`$ scale. The light quark mass parameters ($`\kappa _{\mathrm{crit}}`$ etc) have been determined by the UKQCD collaboration .

The vector and axial currents used in this calculation are improved and renormalised as follows:

$`V_\mu ^R`$ $`=`$ $`Z_V(1+b_Vam_q)\left(V_\mu +ac_V\stackrel{~}{}^\nu T_{\mu \nu }\right)`$ (1)
$`A_\mu ^R`$ $`=`$ $`Z_A(1+b_Aam_q)\left(A_\mu +ac_A\stackrel{~}{}_\nu P\right)`$ (2)
## 2 DECAY CONSTANTS

The results for the decay constants can be found in , as well as a discussion of the improvement coefficients, but we present the main results again. The definition of the vector decay constants differs from the pseudoscalar and so we give it here:

$$0|V_\mu ^\mathrm{R}(0)|\mathrm{V}=iϵ_\mu \frac{M_V^2}{f_V}$$
(3)
The decay constants were extracted from appropriate ratios of correlation functions, and then extrapolated linearly in the light quark mass to $`\kappa _n`$ and $`\kappa _s`$. The extrapolation to the $`B`$ meson is motivated by HQS, of the form

$$\mathrm{\Phi }_\mathrm{i}(M_i)C(M_i)f_i\sqrt{M}_i=\gamma _\mathrm{i}\left(1+\frac{\delta _\mathrm{i}}{M_i}+\frac{\eta _\mathrm{i}}{M_i^2}\right)$$
(4)
where $`i`$ is $`P`$ or $`V`$111In this case $`f_i=M_V/f_V`$ and the function $`C(M_i)`$ contains the leading logarithmic $`\alpha _s`$ corrections to the scaling relation

The results for the decay constants are shown in table 2. The large value of the decay constants is in part due to the choice of quantity used to set the scale. $`f_B`$ varies from $`218`$ MeV using $`r_0`$ to set the scale to $`186`$ using $`m_\rho `$. This dominates the systematic uncertainty. The remaining systematic error comes from difference between $`\beta =6.2`$ and $`\beta =6.0`$, and choice of functional form for the heavy extrapolation (quadratic vs linear). According to HQS, the $`\mathrm{\Phi }`$ functions for the $`P`$ and $`V`$ decay constants should agree in the heavy quark limit. However, the $`\beta =6.0`$ data fail to satisfy this constraint, which may be indicative of large discretisation errors on the coarser lattice. To examine this further we considered, the KLM norm 222modified to account for $`𝒪(a)`$ improvement and the dispersion relation. At $`\beta =6.2`$, the KLM factor was close to unity and the dispersion relation showed good continuum behaviour. However, at $`\beta =6.0`$, the KLM factor has a finite effect, and the dispersion relation shows a deviation from the continuum version. These all suggest that the $`\beta =6.0`$ data are being affected by $`𝒪((am_Q)^2)`$ errors.

## 3 RENORMALISATION

The matrix element of a $`PP`$ transition can be parameterised in terms of two form factors,

$`Z_V^{\mathrm{eff}}P(\stackrel{}{p_2})|V_\mu |P(\stackrel{}{p_1})^{\mathrm{latt}}`$ (5)
$`=f_+(q^2)(p_1+p_2\mathrm{\Delta }_{m^2}q)_\mu +f_0(q^2)\mathrm{\Delta }_{M^2}q_\mu `$
where $`q=p_1p_2`$ and $`\mathrm{\Delta }_{m^2}`$ is the difference of the square of the masses normalised by $`q^2`$. Recalling that the vector current is conserved, we know that $`f_+(0)`$ must be one and therefore consider the forward, degenerate matrix element. Here, $`Z_V^{\mathrm{eff}}`$ can then be measured directly by taking the appropriate ratio of two- to three-point correlation functions. This can then be compared to $`Z_V(1+b_Vamq)`$ with the values of $`Z_V`$ and $`b_V`$ determined by , at the values of quark mass for which we measure $`Z_V^{\mathrm{eff}}`$. This is shown in figure 1.

These data have already been shown for $`\beta =6.2`$ in . The ALPHA determinations of $`Z_V`$ and $`b_V`$ used chiral ward identities, and so used zero or very small quark mass. Our results show excellent agreement with the mass-independent scheme even at large quark mass.

## 4 THE SOFT PION RELATION

HQS can be used to relate the scalar form factor, $`f_0`$ to the ratio of the $`B`$ to pion decay constant,

$$f_0(q_{\mathrm{max}}^2)=\frac{f_B}{f_\pi }$$
(6)
A recent review  found that the majority of lattice calculations found this relation violated. However, the key issue for the SPR on the lattice is the chiral extrapolation. This is complicated by the dependence of the form factors on $`q^2`$, itself heavily dependent on the light quark mass. In this calculation the active light quark (A), which occurs in the heavy-light current and the spectator light quark (S) which takes no part in the weak decay have different mass. Hence,

$$f=f(q_{(m_A,m_S)}^2,m_A,m_S)$$
(7)
A Taylor expansion of $`q^2`$ shows $`q^2M_{P_S}^{\mathrm{light}}`$, and we know from PCAC that $`(M_{P_S}^{\mathrm{light}})^2\overline{m}`$ where the average quark mass $`\overline{m}=(m_A+m_S)/2`$. Thus the light quark mass dependence of the form factor is

$$f(q^2,\overline{m},m_S)=\alpha +\beta \overline{m}^{1/2}+\gamma \overline{m}+\delta m_S$$
(8)
The second term in this expression makes the extrapolation difficult to control, as it is non-analytic and $`\overline{m}^{1/2}/\overline{m}`$ blows up as we approach the chiral limit, in the region where we have no data.

Instead of extrapolating at fixed pion three-momentum we separate the $`q^2`$ and chiral behaviour of the form factor. We first interpolate in $`q^2`$ at fixed quark mass so that at each quark mass we have the form factors at the same momentum value, and then extrapolate at fixed $`q^2`$ without the troublesome second term. However, because we extrapolate at fixed $`q^2`$, we no longer have data at $`q_{\mathrm{max}}^2`$ for physical quark masses. This method has already been successfully employed to calculate the differential decay rate for $`B\pi `$ at $`\beta =6.2`$ in . In this analysis, an extra-momentum channel is included ($`1\sqrt{2}_{}`$), and the range of $`q^2`$ altered to allow interpolation at fixed quark mass for both values of $`\beta `$, in the same range of $`q^2`$.

The form of interpolating function is motivated by pole-dominance models.

$$f_i(q^2)=\frac{f_i}{(1q^2/M_i^2)}$$
(9)
with kinematic constraint $`f_0(0)=f_+(0)`$ imposed. We have also tried other functions for the interpolation, including dipole/pole for $`f_+`$ and $`f_0`$ respectively and with or without the kinematic constraint. This is shown in figure 2.

Also shown in figure 2 are other less well motivated functions, such as linear or quadratic in $`q^2`$. As we are interpolating any reasonable, smooth function should produce similar results.

For the extrapolation to the $`B`$ meson mass, we are again guided by HQS,

$$C(M,M_B)f_i(vk)M^{n_i/2}=\epsilon _i\left(1+\frac{\zeta _i}{M}+\frac{\xi _i}{M^2}\right)$$
(10)
where $`n_i=1,+1`$ when $`i=+,0`$ and $`vk`$ is the kinematic variable,

$$vk=\frac{M^2+M_\pi ^2q^2}{2M}$$
(11)
The fixed $`q^2`$ method allows us to choose the $`q^2`$ values at each heavy quark mass such that $`vk`$ remains constant during the heavy quark extrapolation.

We can now consider the momentum dependence of the form factors at the $`B`$ scale. Combining HQS scaling relations and pole dominance models suggest the momentum dependence is a dipole for $`f_+`$ and a pole for $`f_0`$ (DPP). A more physical model has been suggested by Becirevic and Kaidalov (BK)

$`f_+(q^2)`$ $`=`$ $`{\displaystyle \frac{c_B(1\alpha )}{(1q^2/M_B^{}^2)(1\alpha q^2/M_B^{}^2)}}`$
$`f_0(q^2)`$ $`=`$ $`{\displaystyle \frac{c_B(1\alpha )}{(1q^2/\beta M_B^{}^2)}}.`$ (12)
These are shown in figure 3.

The value obtained for $`f_0(q_{\mathrm{max}}^2)`$ can then be compared to the value of $`f_B/f_\pi `$. To try and estimate the systematic errors we considered the following; $`r_0`$ versus $`m_\rho `$ to set the scale, $`\beta =6.2`$ versus $`\beta =6.0`$, different models for the fixed $`q^2`$ interpolation, quadratic vs linear heavy extrapolation. We see the SPR satisfied on the lattice, with large systematic errors, but with a tendency for $`f_0(q_{\mathrm{max}}^2)`$ to lie below $`f_B/f_\pi `$.

We can also check to see if the extrapolations in $`q^2`$ and heavy quark mass commute. In figure 4 the heavy extrapolation of $`f_B/f_\pi `$ and $`f_0(q_{\mathrm{max}}^2)`$ are shown after $`f_0(q^2)`$ has been extrapolated with a pole model to $`q_{\mathrm{max}}^2`$.

Again we see the SPR is satisfied.

A recent calculation  using a NRQCD formulation of heavy quarks confirms the difficulties with the chiral extrapolations of the form factors and also adopts a fixed kinetic variable approach, which reduces their observed violation of the SPR. Another calculation  using the relativistic heavy Wilson formulation, and a fixed kinematic variable approach see the SPR satisfied if a quadratic heavy extrapolation is used.

The SPR is rather difficult to satisfy on the lattice because we have to extrapolate the light quark mass into a region where $`q^2`$ changes rapidly. This can be circumvented by using the fixed $`q^2`$ approach. However, lattice calculations can reliably determine the differential decay rate for $`B\pi `$ away from $`q_{\mathrm{max}}^2`$ where the phase space is non-zero to extract $`V_{\mathrm{ub}}`$. This work was supported by EPSRC grant GR/K41663 and PPARC grants PPA/G/S/1999/00022 and PPA/GS/1997/00655. CMM acknowleges PPARC grant PPA/P/S/1998/00255.


text:
# Attractions between charged colloidal spheres mediated by correlated absorbed ions

## Abstract

We describe arrangements of ions capable of producing short-range attractive interactions between pairs of charged colloidal spheres in the low temperature strongly correlated limit. For particles of radius $`R`$ with bare charge $`Z`$ and comparable absorbed charge $`N`$ ($`NZ`$), the correlations contribution to the spheres self-energy scales as $`N^{3/2}/R`$, and as $`N/R`$ for the interaction energy between two touching spheres. We show that the re-arrangement of charges due to polarization plays an insignificant role in the nature and magnitude of the interaction.

Electrostatic attractions between colloids have been reported in a large number of circumstances . Experiments in colloidal suspension confined into a small (almost two dimensional) region  have shown direct evidence of these attractions. Simulations of colloidal suspensions in bulk  and confined systems  have also reported electrostatic attractions. The simulations have explored, in particular, the strongly correlated limit in which the thermal energy of the system is smaller than the typical interaction energy between the charges absorbed to the colloids.

Recent studies on colloidal suspensions have elucidated the effect of depletion forces  and of Van der Waals interactions , as well as the importance of including correlations  to predict attractions. In this letter we derive analytic expressions for the interaction energy between two colloids in the strongly correlated limit by considering the distribution of absorbed counterions around them as they approach. Though correlations between absorbed counterions are known to induced short range attractions between charged surfaces , and between aligned rods , it has been noted  that the geometry change to colloids is not trivial.

Systems involving large charged objects (polyions) whose charge is compensated by mobile low-valence ions can be study by considering first the electrostatic ground state of a subsystem that consists of one or few polyions with a fixed number of counterions. Once the ground state is constructed, is possible to draw conclusions about the system by introducing thermal fluctuations and setting the subsystem in contact with the rest of the system which then acts as a thermal and particle reservoir. This approach is justified by the fact that there is a small region surrounding the polyion in which the electrostatic interactions are so dominant that thermal fluctuations are only a minor correction. It has been used, for example, to study the interactions between charged surfaces  and between charged linear polyelectrolytes , the conditions of solubility of DNA and other linear polyelectrolytes in multivalent salt solutions , and to describe the differences between charging colloids or surfaces with monovalent or multivalent counterions . Results from these theories were used in Ref.  to obtain simple but accurate predictions that explain many colloidal simulations results.

While in principle it might be possible to solve the ground state of a system consisting of one charged sphere and a set of $`N`$ absorbed point-charges, such a ground state has to be solved almost independently for each different $`N`$. It would be even more difficult to obtain exact results for the case of two spheres, although numerical results can be obtained for $`N`$ not too large. Instead of attempting such solutions, we propose an approximate solution that satisfies a number of desirable physical conditions. We first decompose the charge distribution on the sphere into a multipolar series for each sphere. We then truncate the series and propose values for the dominant terms. Using this ansatz, we obtain an approximation for the interaction energy between two spheres.

Our main results are simple to understand. We recover the expected negative self-energy for a single charged sphere due to the discrete nature of the absorbed counterions, which is of the form $`N^{3/2}/R`$ . Indeed, since $`RN^{1/2}a`$, where $`a`$ is the typical distance between the absorbed charges, this self-energy is roughly proportional to $`N/a`$, as in the case of a flat surface. If the spheres have bare charge $`Z`$, and $`N`$ absorbed (opposite) charges, the effective charge is $`Q=ZN`$, and thus there is an interaction energy between two spheres at a distance $`r`$ of the form $`Q^2/r`$, but there are important contributions arising from the particle nature of the absorbed charge. When the spheres are in contact there is a region in which the charges from both spheres are strongly correlated. In this region the bare and absorbed charge density are doubled and the typical distance between charges is now of the order of $`a/\sqrt{2}`$. The charges in both spheres remain correlated up to when the separation between corresponding points in the spheres are separated by a distance of the order $`a`$. This defines an effective area of contact between the spheres given by $`\pi Ra`$, that contains $`Na/4R`$ charges. The interaction energy is then negative and of order $`N/R`$. Thus, there can be a very large attractive component to the interaction between the spheres capable of producing a net attraction between them for almost neutral spheres ($`NZ`$).

We assume that the total bare charge of the colloidal sphere $`Z`$ is uniformly distributed on its surface or, equivalently, concentrated at the center of the sphere. The surface charge density of the N absorbed particles is then written as $`\rho (𝐫)=\mathrm{\Sigma }_i\delta (𝐫𝐫_i)`$ where $`𝐫_i`$ are vectors pointing from the center of sphere to a position at the surface. The multipolar moments $`M_n^m`$ for this distribution are

$$M_n^m=\sqrt{\frac{4\pi }{2n+1}}𝑑𝐫Y_n^m(\widehat{𝐫})r^n\rho (𝐫)$$
(1)
where $`Y_n^m`$ are spherical harmonics. The charge density can be then re-expressed as the multipolar expansion:

$$\rho (𝐫)=\underset{n,m}{}\sqrt{\frac{2n+1}{4\pi }}M_n^m\frac{\delta (rR)}{R^{n+2}}Y_n^m(\widehat{𝐫}).$$
(2)
Using Maxwell’s representation  of the spherical harmonics it is also possible to replace the surface density by a singular distribution placed at the center of the sphere given by

$$\rho (𝐫)=\underset{n,m}{}\frac{(1)^{n|m|}M_n^{\pm m}_z^{n|m|}(_xi_y)^{|m|}\delta (𝐫)}{((nm)!(n+m)!)^{1/2}}.$$
(3)
We will use the more compact notation $`\rho =M_n^m\{n,m\}`$ where the brace $`\{n,m\}`$ stands in for the suitable functional form in this expression. The potential generated by this charge density outside of the sphere, at $`r>R`$, is

$$\varphi (𝐫)=\mathrm{\Sigma }_{n,m}\sqrt{\frac{4\pi }{2n+1}}M_n^mY_n^m(\widehat{𝐫})\frac{1}{r^{n+1}}.$$
(4)
In the low temperature limit we expect the charges to occupy positions that up to permutations and rotations of the whole sphere are uniquely determined. For large $`N`$, the inhomogeneities arising from the particle nature of the charges in this static ground-state distribution should create only negligible contributions for the first few multipolar moments other than $`n=0`$, with the next non-vanishing moments arising at $`n`$ values that start to resolve distances of order $`a`$, namely $`nN^{1/2}`$. To make analytical calculations feasible we truncate the multipolar expansion and represent the sum of the $`N`$ delta functions as a smooth function with $`N`$ clearly distinguished maxima. Since the interaction between the absorbed charges is repulsive one can expect four nearest neighbors per particle in the ground state and we also require this from our truncated expansion. Our ansatz is then a truncation of the series Eq.(2) that only consider moments with $`n`$ values up to $`n=3l/2`$ with $`l=\sqrt{2N}`$. Further, within this range we take only four moments to be non zero:

$`\rho =N\{0,0\}`$ $`+`$ $`{\displaystyle \frac{(2N)^{3/4}R^l}{\pi ^{1/2}}}\{l,0\}+`$ (6)
$`{\displaystyle \frac{2^{7/8}N^{5/8}R^l}{\pi ^{3/4}}}(\{l,l\}+\{l,l\}).`$
The precise values of the moments given here are discussed below.

The proposed charge density Eq.(6) has $`N`$ maxima, distributed between $`l=\sqrt{N/2}`$ meridians generated by the $`m=l`$ terms. Each meridian has $`l/2`$, non uniformly distributed maxima generated by the $`m=0`$ term. A consistent extension of our ansatz would give the next non-zero moments at $`n`$ and $`\pm m`$ values that are integer multiples of $`l`$. In the limit of infinite radius and near the equator of the sphere the distribution becomes equivalent with the description, by a few Fourier modes, of a square lattice of charges on a plane. This limit is indeed the inspiration for our choice. The major inaccuracies of this distribution are near the poles where the maxima are closer than they should be.

For large $`N`$, and near the equator, the spherical harmonic $`Y_l^0`$ has an oscillatory form $`\pi ^1\mathrm{cos}(l\theta +\psi )`$, with $`\theta `$ the polar coordinate and $`\psi `$ a phase. If we make the position of the charges to coincide with the peaks of this approximate distribution, we obtain, using Eq.(1), $`M_l^0=\pi ^{1/2}(2N)^{3/4}R^l`$. Again, near the equator, we have that $`Y_l^{\pm l}=e^{\pm il\varphi }f(\theta )`$, where $`\varphi `$ is the azimuthal angle, and $`f`$ is a multiple of the Legendre function $`P_l^l`$ that has a maximum in the equator and is important only in a narrow region. We assume that the particles positions coincide with the maxima of the real parts of the waves in the $`\varphi `$ coordinate while the sum over their contributions along the $`\theta `$ is approximated by a multiple of the integral $`f(\theta )𝑑\theta `$ from which $`M_l^l=2^{7/8}(N)^{5/8}R^l/(\pi ^{3/4})`$.

To discuss the interaction between spheres we must consider the redistribution of charges that occurs when they are put in close contact with each other. To first order we can represent this redistribution as a flow generated by the projection, onto the sphere, of a uniform displacement of charges in the direction of the polarization. If the displacement of the charges is $`𝐝`$, the new charge distribution is approximately $`\rho (𝐫)+\mathrm{\Delta }\rho (𝐫)=\rho (𝐫)+𝐝\rho (𝐫)`$. Using the representation Eq. (3) and assuming the displacement to be in the $`\widehat{𝐱}`$ direction, it can be shown that the displacement applied to our ansatz Eq.(6) gives rise to new non-zero extra moments:

$`\mathrm{\Delta }\rho `$ $`=`$ $`{\displaystyle \frac{P}{\sqrt{2}}}\{1,1\}+{\displaystyle \frac{PlM_l^0}{2N}}\{l+1,1\}+{\displaystyle \frac{PlM_l^l}{N}}\{l+1,l+1\}`$ (8)
$`+{\displaystyle \frac{PM_l^l}{\sqrt{2}N}}\{l+1,l1\}+(mm).`$
This equation displays only the positive $`m`$ terms. The magnitude of the generated dipole moment is $`P=dN`$, but as this expression shows, the flow of charges also affects the higher order multipoles. The linear approximation is valid only if the average displacements of the charges are smaller than the typical distance between them $`d<a`$

The energy of interaction between the absorbed charges and the colloid is $`ZN/R`$. If the decomposition of the ground state configuration of the absorbed charges into multipolar moments $`M_n^m`$ is known, their self energy is $`E=_{n,m}M_n^mM_n^m/2R^{2n+1}`$. This expression is infinite since it contains the self-energy of every absorbed charge. The harmonic decomposition of a single point charge shows that its self-energy splits into equal contributions of magnitude $`1/2R`$ for each $`n`$ number giving, symbolically, $`E_{single}=(1/2R)_n1`$. Subtracting $`N`$ times this amount from the energy at each $`n`$ level we obtain:

$$E=\frac{ZN}{R}+\underset{n}{}\left(\frac{N}{2R}+\underset{m}{}\frac{M_n^mM_n^m}{2R^{2n+1}}\right).$$
(9)
Evaluating the energy of our ansatz for the charge distribution and including the possibility of a charge redistribution with dipole moment $`P`$ the main terms in the evaluation of the energy are

$$E=\frac{Z^22ZN+N(N1)}{2R}\alpha \frac{N^{3/2}}{R}+\beta \frac{P^2N^{1/2}}{2R^3},$$
(10)
where $`\alpha =(3/2^{3/2})1/(\pi 2^{1/2})`$, and $`\beta =1/(2^{1/2}\pi )`$. In the second and third terms we have retained only the contributions with highest powers in $`N`$. The monopole contributions have been explicitly separated. In the absence of polarization, this result is clearly consistent with the expected behavior of the self-energy in the $`N\mathrm{}`$ limit that is equivalent to charges in a flat surface in which the self-energy scales as $`N(1/a)N(N/R^2)^{1/2}`$.

To consider two spheres we orient the reference systems for each of them so that their centers lie in the x-axis of both systems. The total distance between centers is $`r`$. The interaction energy can be written as the integral of the potential generated by one sphere $`\varphi `$ times the charge density of the second sphere $`\rho ^{}`$, $`E=𝑑𝐫\varphi \rho ^{}`$. Writing down both $`\varphi `$ and $`\rho ^{}`$ in their multipolar expansions and using the Maxwell representation, the integral becomes a sum of terms that are products of derivatives of $`1/r`$ times derivatives of a delta function at the center of the second sphere. Repeated integration by parts leads to

$$E=\underset{n,m,p,q}{}E_{n,m,p,q}\frac{M_n^mM_p^q}{r^{n+p+1}}$$
(11)
where the coefficient is, using $`s=n+m,t=nm,s^{}=p+q,t^{}=pq,S=s+s^{},T=t+t^{}`$,

$$E_{n,m,p,q}=\left[\frac{S!T!}{s!s^{}!t!t^{}!}\right]^{1/2}\sqrt{\frac{4\pi }{2S+1}}Y_{(S+T)/2}^{(ST)/2}(\frac{\pi }{2},0).$$
(12)
The point of evaluation of the spherical harmonic reflects the direction of the relative position of the spheres in the chosen coordinates, namely the $`\widehat{x}`$ axis. It can be shown that for large moment numbers the magnitude of the coefficient has its maximum values when $`n=p`$ and $`m=\pm q`$, and decays away from this situation as $`\mathrm{exp}((s(S/2))^2(t(T/2))^2)`$. Therefore, the main contributions to the energy from high order multipoles appears when the multipoles in both spheres correspond, that is, when the charge distributions are similar and interlock. In our truncation we are able to pick the relative phases of the multipoles in the spheres to get an attractive energy.

We will consider only the case of two spheres of equal radius. Since we are interested in exploring the possibility of attractive interactions between spheres we consider only the case in which they touch $`r=2R`$. If attractions are present this will be their optimal position. Considering the case of spheres that have equal total charge and magnitude of polarization we obtain the following interaction energy:

$`E_{int}`$ $`=`$ $`{\displaystyle \frac{(ZN)^2}{2R}}\gamma {\displaystyle \frac{N}{2R}}{\displaystyle \frac{(PP^{})(ZN)}{(2R)^2}}+`$ (14)
$`{\displaystyle \frac{PP^{}}{(2R)^3}}(1+2\gamma )`$
where $`\gamma =\pi ^2+2^{3/4}\pi ^{5/2}`$. The total energy is obtained from adding the self-energy and the interaction energy. Minimizing the total energy with respect to the polarization we arrive to the conclusion that the polarization of the spheres should be antiparallel $`𝐏=𝐏^{}`$, and with magnitude $`P=RQ/\beta N^{1/2}`$. The polarization is negligible in the limit of large $`N`$. The reason for this is that a flow of mobile charges towards the point of contact requires a large increase of energy in a large depleted area, while the decrease in energy occurs only in the small region of contact. Thus, we find that the ground state is very rigid, and the interaction between spheres does not alter much its local structure. The polarization corrections can then be ignored and the interaction energy is simply the result given in Eq.(14) without polarization terms. If the effective charge is small, the interaction is dominated by the high multipole correlations and it scales as $`N/R`$. The rigidity of the local structure makes this result robust when applied to considerations of more than two spheres.

Given our result for the evaluation of the interaction energy, we are ready to comment on the problem of attraction in more realistic conditions. The Bjerrum number $`B=e^2/a4\pi ϵ_0k_BT`$ compares the electrostatic interaction between charges at a distance $`a`$ with the typical thermal energy $`k_BT`$. We can consider $`z`$-valent counterions and then, if we have $`z^2B>1`$, the correlation energy can be very strong. The exact number of charges in the spheres is controlled by the equilibrium of the energy between absorbed charges and a chemical potential $`\mu `$ for the extraction of them from the volume not occupied by the colloidal spheres. To solve for the number of absorbed charges we can minimize the free energy $`F=E(N)\mu N`$ with respect to $`N`$. Short range attractions between two charged colloids proportional to $`N/R`$ result from the ”interlocking” correlations of the thin layer of counterions at the surface of the colloids in the area of contact. These correlations easily arise, for high $`Bz^2`$, when the number of absorbed counterions is large and comparable to the bare charge. The relative location of spheres influences the amounts of absorbed charges but their specific values are mostly dependent on the interaction parameter $`Bz^2`$ and on the chemical potential $`\mu `$. This suggests that it is more likely to observe attractions in solutions of high valence counterions or in systems where the accessible volume for the counterions is small as in the case of colloids in confined geometries. Though in principle the confinement of colloidal suspensions by weakly charged surfaces imposes non-trivial boundary conditions that affect the local electrostatic fields, we can always analyze the local configurations of ions under these externally imposed conditions. If the external conditions are successful in creating attractions, the local arrangement of charges are as we have described here.

The explanation for the presence of long-range attractions put forward in Ref.  is that there are attractive metastable states characterized by pairs of spheres with opposite total charges. The overcharging of some spheres is possible due to the negative correlation self-energy in each sphere Eq.(10). What our calculation shows is that even if the effective charges are equal, the ground state is attractive. The metastable states, on the other hand, reduce the potential barrier to achieve close contacts, and at the free energy level might indeed provide a long range attractive interaction. At short distances we expect the charge to redistribute to the type of configurations discussed here.

We thank Cherry Murray for useful conversations. This work was sponsored by the National Science Foundation, grants DMR9807601 and DMR9632472.


text:
# Classical Correlation-Length Exponent in Non-Universal Quantum Phase Transition of Diluted Heisenberg Antiferromagnet

\[

## Abstract

Critical behavior of the quantum phase transition of a site-diluted Heisenberg antiferromagnet on a square lattice is investigated by means of the quantum Monte Carlo simulation with the continuous-imaginary-time loop algorithm. Although the staggered spin correlation function decays in a power law with the exponent definitely depending on the spin size $`S`$, the correlation-length exponent is classical, i.e., $`\nu =4/3`$. This implies that the length scale characterizing the non-universal quantum phase transition is nothing but the mean size of connected spin clusters.

\]

Ground-state phase transitions in two-dimensional (2D) diluted quantum Heisenberg antiferromagnets (HAF’s) have attracted much interests because they are caused by the coexistence of quantum fluctuations and randomness. Numerical works and theoretical works have given various estimates of the critical concentration of the system with $`S=1/2`$. All of them are above the purely-geometrical percolation threshold on a square lattice, $`p_{\mathrm{cl}}`$=0.5927460(5). This suggests that the phase transition could be seriously affected by quantum fluctuations.

Recently, Kato et al. have investigated the diluted 2D HAF with $`S=`$ 1/2, 1, 3/2, and 2 by means of the quantum Monte Carlo (QMC) method with the continuous-imaginary-time loop algorithm. Their conclusion is qualitatively different from the above mentioned results: the critical concentration coincides with $`p_{\mathrm{cl}}`$ and does not depend on $`S`$. The coincidence of the critical concentration with $`p_{\mathrm{cl}}`$ has also been reported on the bond-diluted HAF.

The critical exponents of the phase transition at $`p_{\mathrm{cl}}`$ have also been estimated. Kato et al. have obtained the critical exponent $`\beta `$ of the zero-temperature staggered magnetization. Interestingly, the value of $`\beta `$ is different from that of the classical ($`S=\mathrm{}`$) exponent and depends on $`S`$. They have also estimated other critical exponents by the finite-size scaling (FSS) analysis assuming the following form for the static staggered structure factor at zero temperature:

$$S_\mathrm{s}(L,0,p)L^\mathrm{\Psi }\stackrel{~}{S}_\mathrm{s}(L^{1/\nu }(pp_{\mathrm{cl}})),$$
(1)
where $`L`$ is the system size and

$$S_\mathrm{s}(L,T,p)\frac{1}{L^d}\underset{i,j}{}\mathrm{e}^{\mathrm{i}\stackrel{}{k}(\stackrel{}{r}_i\stackrel{}{r}_j)}S_i^zS_j^z$$
(2)
with $`\stackrel{}{k}=(\pi ,\pi )`$ and $`d=2`$. The bracket $`\mathrm{}`$ in Eq. (2) denotes both the thermal and random averages. The scaling function $`\stackrel{~}{S}_\mathrm{s}(x)`$ in Eq. (1) has the following asymptotic form:

$`\stackrel{~}{S}_\mathrm{s}(x)`$ $``$ $`\{\begin{array}{cc}x^{2\beta }\hfill & \mathrm{for}x1\hfill \\ |x|^{\nu \mathrm{\Psi }}\hfill & \mathrm{for}x1,\hfill \end{array}`$ (5)
where the exponent $`\beta `$ is related to $`\mathrm{\Psi }`$ and $`\nu `$ by the scaling relation

$$2\beta =(d\mathrm{\Psi })\nu .$$
(6)
They have estimated the critical exponent $`\mathrm{\Psi }`$ by the FSS analysis exactly at $`p=p_{\mathrm{cl}}`$ and have found that it depends on $`S`$.

Kato et al. have attributed the $`S`$-dependence in exponent $`\beta `$ or $`\mathrm{\Psi }`$ to quantum fluctuations, while the length exponent $`\nu `$ is assumed to be given by the classical one, $`\nu _{\mathrm{cl}}=4/3`$, which governs a power-law divergence of the mean size of connected spin clusters as $`\lambda (p)|pp_{\mathrm{cl}}|^{\nu _{\mathrm{cl}}}`$. Namely, they have assumed that the staggered spin correlation function between two spins in a cluster is described by the scaling expression

$$C(i,j;p)r_{i,j}^\alpha \stackrel{~}{C}(r_{i,j}/\lambda (p))$$
(7)
with $`\stackrel{~}{C}(x)`$ const. at $`x1`$. The power-law decay of the correlation function is due to quantum fluctuations and its $`S`$-dependent exponent $`\alpha `$ is related to $`\mathrm{\Psi }`$ by

$$\mathrm{\Psi }=2Dd\alpha ,$$
(8)
where $`D`$ is the fractal dimension ($`91/48`$ for $`d=2`$). Kato et al. have checked this scenario simply by evaluating $`\nu `$ through Eq. (6) using $`\beta `$ and $`\mathrm{\Psi }`$ obtained by their simulation.

In the present paper, we perform the FSS analysis on $`S_\mathrm{s}(L,0,p)`$ of systems with $`S=1/2`$ and 1 more systematically by carrying out the QMC simulation at various concentrations not only at $`p_{\mathrm{cl}}`$. It is found, as we see below, that our QMC data are well described by the FSS form of Eq. (1) in a whole range of $`p`$ studied, which includes the asymptotic ranges $`x1`$ and $`x1`$ in Eq. (5). Within numerical accuracy of our analysis, exponent $`\mathrm{\Psi }`$ turns out to definitely depend on the spin size $`S`$, while exponent $`\nu `$ does not. The latter $`S`$-independent value of $`\nu `$ coincides with the classical one (=4/3) within the error bar. These results further support the above mentioned scenario of the quantum phase transition in the diluted 2D HAF, particularly, the ansatz that the length scale characterizing the transition is nothing but the mean size of connected spin clusters.

The system we study is the site-diluted HAF on a square lattice described by the Hamiltonian

$$=J\underset{i,j}{}ϵ_iϵ_j𝐒_i𝐒_j,$$
(9)
where $`J(>0)`$ is the antiferromagnetic coupling constant, $`_{i,j}`$ denotes the summation over all nearest-neighbor pairs and $`𝐒_i`$ is the quantum spin operator at site $`i`$. The quenched magnetic occupation factors $`\{ϵ_i\}`$ independently take 1 or 0 with probability $`p`$ and $`1p`$, respectively. We simulate $`L\times L`$ square lattices with the periodic boundary condition by means of the same QMC method with the continuous-imaginary-time loop algorithm as that adopted by Kato et al. An improved estimator is used to calculate the static staggered structure factor. At each parameter set ($`L,T,p`$) physical quantities of interest are averaged over 1000 – 3000 samples. At each sample, 103 – 104 Monte Carlo steps (MCS) are spent for measurement after 500 – 103 MCS for thermalization.

In a finite system, $`S_\mathrm{s}(L,T,p)`$ converges rapidly to its zero-temperature values at temperatures lower than the gap due to the finiteness of the system. The saturation temperature turns out to be smaller for smaller $`|pp_{\mathrm{cl}}|`$. In the present work $`S_\mathrm{s}(L,0,p)`$ close to $`p_{\mathrm{cl}}`$ is approximated by $`S_\mathrm{s}(L,T,p)`$ at the following temperatures: for the $`S=1/2`$ system $`T=0.002J`$ for $`L=24`$ and $`T=0.001J`$ for $`L=32`$, 40, and 48, and for the $`S=1`$ system $`T=0.01J`$ for $`L=24`$ and $`T=0.005J`$ for $`L=32`$, 40, and 48.

Let us first examine the FSS analysis of $`S_\mathrm{s}(L,0,p)`$ at $`p`$ close to $`p_{\mathrm{cl}}`$. The results of the $`S=1/2(S=1)`$ system at $`0.580p0.605`$ ($`0.585p0.600`$) are shown in Fig. 1(a) (1(b)). The error bars in the figures represent the standard deviation. In the FSS fit the critical concentration is set to be $`p_{\mathrm{cl}}`$ ($`0.5927460`$), and the scaling function $`\stackrel{~}{S}_\mathrm{s}(x)`$ is approximated by a polynomial of order 2. As seen in Fig. 1(a), the QMC data for $`S=1/2`$ are well scaled with $`\mathrm{\Psi }=1.192`$ and $`\nu =1.23`$. The statistical accuracy of the fit is shown in Fig. 2, where we draw the confidence region within which the true values of $`\mathrm{\Psi }`$ and $`\nu `$ fall with probability 68.3% (1-$`\sigma `$), 95.3% (2-$`\sigma `$), or 99.7% (3-$`\sigma `$). Similarly, the data for $`S=1`$ are well scaled with $`\mathrm{\Psi }=1.555`$ and $`\nu =1.09`$ as seen in Figs. 1(b) and 2. From these results we can conclude that exponent $`\mathrm{\Psi }`$ definitely depends on the spin size $`S`$. The exponents $`\nu `$ of $`S=1/2`$ and 1, on the other hand, coincide with each other and with its classical value ($`=4/3`$) within the numerical accuracy of the present analysis. The obtained values of $`\mathrm{\Psi }`$, $`\nu `$, $`\alpha `$, and $`\beta `$ are summarized in Table I. The exponents $`\alpha `$ and $`\beta `$ are calculated by the scaling relations, Eqs. (6) and (8).

The value of $`\mathrm{\Psi }`$ for $`S=1/2`$ obtained above is consistent with the one ($`=1.17(6)`$) that has been estimated by Kato et al. from $`S_\mathrm{s}(L,0,p)`$ at $`p=p_{\mathrm{cl}}`$. In their analysis $`S_\mathrm{s}(L,0,p_{\mathrm{cl}})`$ is approximated by $`S_\mathrm{s}(L,T,p_{\mathrm{cl}})`$ at low temperatures where its $`T`$-dependence becomes not discernible within the error bars. Kato et al. have also performed the FSS analysis making use of the data at all temperatures they have simulated. This analysis yields $`\mathrm{\Psi }=1.27(2)`$ which differs distinctly from the present result listed in Table I. This discrepancy may be due to the systematic error in our estimates where $`S_\mathrm{s}(L,0,p)`$ is approximated by $`S_\mathrm{s}(L,T,p)`$ at a small but finite temperature as described before. It should be remedied when the finite-temperature FSS analysis is carried out as done by Kato et al.

From Fig. 2 one sees that the accuracy of $`\nu `$ is significantly poorer than that of $`\mathrm{\Psi }`$, which is interpreted as follows. The value of $`\mathrm{\Psi }`$ can be essentially extracted solely from $`S_\mathrm{s}(L,0,p_{\mathrm{cl}})`$, while to evaluate $`\nu `$ we need $`S_\mathrm{s}(L,0,p)`$ at $`p`$ other than $`p_{\mathrm{cl}}`$, or $`\delta S_\mathrm{s}(L,0,p)/\delta p`$ at $`p=p_{\mathrm{cl}}`$. Naturally, the statistical error of $`\nu `$ is expected to be larger than that of $`\mathrm{\Psi }`$. Concerned with the systematic error in $`\nu `$ arising from finite-temperature corrections, on the other hand, the saturation temperature of $`S_\mathrm{s}(L,T,p)`$ becomes smaller for smaller $`|pp_{\mathrm{cl}}|`$, i.e., the effect of finite-temperature corrections to $`S_\mathrm{s}(L,T,p)`$ becomes maximal exactly at $`p=p_{\mathrm{cl}}`$. This implies that the leading order finite-temperature correction to $`\delta S_\mathrm{s}(L,0,p)/\delta p`$ at $`p=p_{\mathrm{cl}}`$ disappears, and so the systematic error in $`\nu `$ is much smaller than that in $`\mathrm{\Psi }`$.

Next let us discuss $`S_\mathrm{s}(L,0,p)`$ in the full range of $`p`$ we have simulated, i.e., $`0.2p1`$, which is much wider than that in Fig. 1. As shown in Figs. 3(a) and (b) respectively for $`S=1/2`$ and 1, all QMC data turn out to lie on a universal curve when $`L^\mathrm{\Psi }S_\mathrm{s}(L,0,p)`$ are plotted against $`L^{1/\nu }|pp_{\mathrm{cl}}|`$ by using the exponents $`\mathrm{\Psi }`$ and $`\nu `$ listed in Table I. For $`L^{1/\nu }|pp_{\mathrm{cl}}|>1`$, the data points merge to the dashed line, which represents $`ax^{2\beta }`$ for $`p>p_{\mathrm{cl}}`$ and $`b|x|^{\nu \mathrm{\Psi }}`$ for $`p<p_{\mathrm{cl}}`$, where $`\beta `$, $`\nu `$, and $`\mathrm{\Psi }`$ are those listed in Table I, and $`a`$ and $`b`$ are arbitrary constants adjusted to fix the position of the asymptotic lines. At much larger $`L^{1/\nu }|pp_{\mathrm{cl}}|`$ the scaling fit becomes deteriorated, indicating that the corresponding $`p`$ is out of the scaling region. The results shown in Fig. 3 also support the scenario due to Kato et al. on the non-universal quantum phase transition.

In summary, in order to establish nature of the non-universal quantum phase transition of 2D site-diluted HAF’s for $`S=1/2`$ and $`1`$, we have performed the QMC simulation in a relatively wider region of concentration than that of Kato et al. and have estimated the critical exponent $`\nu `$ more systematically. We have observed that the static staggered structure factor is well described by the scaling form of Eq. (1) with Eq. (5). In particular, the exponent $`\nu `$ is confirmed to coincide with the classical one. These results support the following arguments by Kato et al.: 1) there exists no other macroscopic characteristic lengths than $`\lambda (p)`$, the mean size of connected spin clusters at concentration $`p`$ and 2) the staggered spin correlation function between two sites on a fractal cluster decays in a power law as $`C(i,j;p)r_{i,j}^\alpha `$, where $`\alpha `$ depends on the strength of quantum fluctuations specified by the spin size $`S`$.

Most of numerical calculations for the present work have been performed on the CP-PACS at University of Tsukuba, Hitachi SR-2201 at Supercomputer Center, University of Tokyo, and SGI2800 at Institute for Solid State Physics, University of Tokyo. The present work is supported by the “Large-scale Numerical Simulation Program” of Center for Computational Physics, University of Tsukuba, and also by the “Research for the Future Program” (JSPS-RFTF97P01103) of Japan Society for the Promotion of Science.


text:
# Shock structures in time averaged patterns for the Kuramoto-Sivashinsky equation

## Abstract

The Kuramoto-Sivashinsky equation with fixed boundary conditions is numerically studied. Shocklike structures appear in the time-averaged patterns for some parameter range of the boundary values. Effective diffusion constant is estimated from the relation of the width and the height of the shock structures.
PACS numbers:05.45.-a, 47.54.+r

An ordered structure can be obtained by time averages of disordered spatiotemporal patterns. Such ordered structures are experimentally found in Faraday waves, rotating thermal convection and electroconvection \[1-3\]. These ordered structures seem to come from the boundary effect. The Kuramoto-Sivashinsky equation is one of the simplest partial differential equations exhibiting spatiotemporal chaos. Average patterns of spatiotemporal chaos in the Kuramoto-Sivashinsky equation were studied in . We study average patterns of the Kuramoto-Sivashinsky equation with different boundary conditions. The equation in one dimension has the form

$$h_t=h_{xx}h_{xxxx}+\frac{1}{2}(h_x)^2,$$
(1)
where $`h=h(x,t)`$ is a real function, $`x[0,L]`$, and the subscripts stand for derivatives. An equivalent equation is obtained for $`u=h_x`$ as

$$u_t=u_{xx}u_{xxxx}+uu_x.$$
(2)
Equation (2) possesses the Galilean symmetry$`(xx+vt,uu+v)`$. The Kuramoto-Sivashinsky equation exhibits spatiotemporal chaos and it is conjectured that the large-scale properties of the equation are described by the stochastic Burgers equation \[5-8\].

$$u_t=Du_{xx}+\lambda uu_x+\xi ,$$
(3)
where $`D>0`$, $`\xi (x,t)`$ denotes the noise term, and the coefficient $`\lambda `$ of the nonlinear term is expected to be 1 owing to the same Galilean symmetry as the Kuramoto-Sivashinsky equation. The Fourier transform $`\xi _k(t)`$ of $`\xi (x,t)`$ satisfies

$$<\xi _k(t)\xi _k^{}(t^{})>=Tk^2\delta _{k,k^{}}\delta (tt^{}).$$

Characteristic solutions to the Burgers equation without the last noise term are the shock solutions:

$$u(x)=A\mathrm{tanh}(\kappa x),$$

where $`\kappa =\lambda A/(2D)`$. We have performed numerical simulations to seek for some shocklike structures in large-scale averaged patterns for the Kuramoto-Sivashinsky equation. We have used a simple real-space scheme obtained by replacing spatial and temporal derivatives by finite differences. The spatial discretization step is $`\delta x=0.25`$ and the temporal discretization step $`\delta t=2.5\times 10^4`$, and the system size $`L=400`$. We have considered fixed boundary conditions, that is,

$$u(0,t)=U,u(L,t)=U,u_x(0,t)=u_x(L,t)=0,$$

where $`\pm U`$ are the fixed values at the boundaries and $`U`$ is changed as a parameter. The initial condition is

$`u(x,0)`$ $`=`$ $`U+Ux/(2l),\mathrm{for}\mathrm{\hspace{0.33em}0}<x<l`$
$`=`$ $`U(xL/2)/(L2l),\mathrm{for}l<x<Ll`$
$`=`$ $`U+U(xL)/(2l),\mathrm{for}Ll<x<L`$
where $`l=4`$ is assumed. The simulation is performed until $`t=t_2=1500`$ and the temporal average is calculated as

$$\overline{u}(x)=(1/N)\underset{n=1}{\overset{N}{}}u(x,t_1+n\mathrm{\Delta }t)$$

where $`t_1=500`$, $`\mathrm{\Delta }t=50\delta t`$ and $`N=(t_2t_1)/\mathrm{\Delta }t`$.

Figure 1 displays snapshot patterns at $`t=1500`$ and the averaged pattern. In the central region of the averaged pattern, $`\overline{u}`$ takes nearly 0 and the mean slope of $`\overline{u}`$ is very small for $`U<1.3`$. This type of averaged patterns is obtained for $`U=0`$ in .

A stationary shocklike pattern appears at $`U=1.65`$. Since the solution is stationary in time, the averaged pattern and the snapshot pattern overlap in Fig. 1(c). This stationary solution can be understood from the absolute stability of a constant solution as follows. A constant solution $`u(x,t)=U`$ is convectively unstable, but it becomes absolutely stable for $`U>1.622`$. The stability of the constant solution $`u(x,t)=U`$ is determined from the linear dispersion relation for the perturbation $`\delta u(x,t)u_ke^{ikxi\omega t}`$ around $`u(x,t)=U`$ . The linear dispersion relation is written as

$$\omega (k)=Uk+i(k^2k^4).$$
(5)
The solution is absolutely stable under the condition that $`d\omega /dk=U+2i(kk^3)=0`$ and $`\mathrm{Im}\omega (k)=U\mathrm{Im}k+\mathrm{Re}(k^2k^4)<0`$ where $`k`$ is generally a complex number. The transition occurs at $`U=U_c`$ where $`U_c`$ satisfies

$$27U_c^468U_c^28=0,$$

that is, $`U_c=\{(34+\sqrt{34^2+827})/27\}^{1/2}1.622`$. The convectively unstable but absolutely stable solution is observed in our system, since the fixed boundary conditions are assumed and no disturbance is created at the boundaries. A solution which connects the two absolutely stable solutions $`u(x,t)=U`$ and $`u(x,t)=U`$ appears as a shocklike structure with oscillating tails in Fig. 1(c).

Another type of average patterns appears for $`1.3<U<1.6`$. Two nearly flat regions are observed and a shocklike structure connecting the two flat regions appears at $`U=1.45`$ as shown in Fig.1(b). The averaged value of $`u(x,t)`$ in the flat region is not $`\pm U`$ and $`u(x,t)`$ exhibits spatiotemporal chaos around the averaged value. We have assumed the form of the averaged pattern as $`A\mathrm{tanh}(\kappa x)`$ and estimated the value of $`A`$ and $`\kappa `$. Figure 2 displays the averaged patterns and the corresponding profile $`\overline{u}(x)=A\mathrm{tanh}(\kappa x)`$ at $`U=1.33`$ and $`1.45`$. As $`U`$ is increased, $`A`$ is increased and $`\kappa `$ is decreased. We have observed this shocklike structure in the time-averaged pattern even if the initial value of $`u(x,0)`$ is a random number between -0.02 and 0.02. We have checked up that the shocklike structure appears in the time-averaged pattern for a larger system with $`L=3200`$. It is fairly robust, but we are not sure that the shocklike structure is maintained in an infinite size system.

Figure 3 displays the height $`A`$ of the shocklike structure as a function of $`U`$. The diamonds denote simulation results for $`L=400`$ and the crosses denote simulation results for $`L=1600`$. The shock amplitude $`A`$ seems to be independent of the system size for sufficiently large $`L`$. The height $`A`$ seems to increases from zero continuously, however, the transition is not clear, since the spatial fluctuation of the averaged pattern is large near $`U1.3`$. Figure 4 displays the effective diffusion constant $`D=\lambda A/(2\kappa )`$ at $`L=400`$ (denoted by diamonds) and $`L=1600`$ (denoted by crosses), where $`\lambda =1`$ is assumed. In our simulations $`A,\kappa `$ and $`D`$ depends on $`U`$. We have plotted a relation of $`D(U)`$ v.s. $`\kappa (U)`$ in Fig.4. As $`\kappa (U)`$ is decreased, $`D(U)`$ seems to increase. It might be interpreted that the effective diffusion constant tends to increase, as the length scale is longer. Sneppen et al. estimated the effective diffusion constant $`D`$ and $`\lambda `$ as $`D10.5`$ and $`\lambda =2.32`$ from the analysis of the interface width $`<(h(x,t)<h>)^2>`$ . Our effective diffusion constant seems to approach about 10 as $`\kappa 0`$. This result may be consistent with the result of Sneppen et al., although they consider a very large system and the effective coefficient $`\lambda `$ of the nonlinear term is not 1 in their analysis.

In summary, we have numerically studied the Kuramoto-Sivashinsky equation with fixed boundary conditions. At both the boundaries, $`u(x,t)`$ are fixed to $`\pm U`$. If the $`U`$ is smaller than a critical value, the averaged pattern takes nearly zero. If $`U`$ is larger than another critical value, a stationary shock pattern with oscillating tails appears owing to the absolute stability. We have found a shocklike structure of the spatiotemporal chaos between the two critical values. We do not understand well the fairly stable structure in the spatiotemporal chaos, but it may be partly related to the weakness of the absolute instability. We are not sure whether the shocklike structures are due to finite-size effects.

In our finite-size simulations, the time averaged patterns take fairly definite nonzero amplitudes $`\pm A`$. The time averaged pattern is approximated with the shock solution of the Burgers equation and the effective diffusion constant is estimated. The effective diffusion constant depends on the length scale of the shock width and it seems to approach about 10 for a very long scale.

This work was supported in part by the Grant-in-Aid for Scientific Research No.11837014 from the Ministry of Education, Science, Sports and Culture.

Figure captions
Fig. 1. Snapshot patterns (solid lines) $`u(x,t)`$ at $`t=1500`$ and the time averaged patterns (dashed lines) $`\overline{u}(x)`$ at $`U=1.25`$ (a), 1.45 (b) and 1.65 (c).
Fig. 2. Time averaged patterns (solid lines) and their fitting with $`\overline{u}=A\mathrm{tanh}(\kappa x)`$ (dashed lines) at $`U=1.33`$ (a) and $`1.45`$ (b).
Fig. 3. Relation of the height $`A`$ v.s. $`U`$ for the time averaged patterns at $`L=400`$ (denoted by diamonds) and $`L=1600`$ (denoted by crosses).
Fig. 4. Relation of the effective diffusion constant $`D`$ v.s. $`\kappa `$ for the time averaged patterns at $`L=400`$ (denoted by diamonds) and $`L=1600`$ (denoted by crosses).


text:
# Damping and the Hartree Ensemble Approximation

## 1 Introduction

Dynamics of quantum fields in real time is a lot more complicated than statics in imaginary time and one has to make approximations before giving the problem to the computer . As an improvement on classical dynamics we are studying the Hartree approximation. The Hartree equations of motion for the $`\phi ^4`$ model in 1+1 dimensions are given by

$`\ddot{\phi }_x`$ $`=`$ $`\mathrm{}\phi _x(\mu ^2+\lambda \phi _x^2+3\lambda C_x)\phi _x,`$
$`\ddot{f}_x^\alpha `$ $`=`$ $`\mathrm{}f_x^\alpha (\mu ^2+3\lambda \phi _x^2+3\lambda C_x)f_x^\alpha ,`$
$`C_x`$ $`=`$ $`{\displaystyle \underset{\alpha }{}}[(1+2n_0^\alpha )|f_x^\alpha |^2],`$
where $`\phi _x=\widehat{\phi }_x`$ is the mean field and the $`f_x^\alpha `$ are a complete set of mode functions used for the parametrisation of the operator field $`\widehat{\phi }_x`$ in the Hartree approximation,

$$\widehat{\phi }_x=\phi _x+\underset{\alpha }{}(\widehat{a}_\alpha f_x^\alpha +\widehat{a}_\alpha ^{}f_x^\alpha ).$$

The initial conditions are specified by a suitable choice for $`\phi _x`$, $`\dot{\phi }_x`$ and the mode functions, as well as the initial occupation numbers $`n_0^\alpha =\widehat{a}_\alpha ^{}\widehat{a}_\alpha `$.

Hartree-like approximations have been widely used, but when the mean field is homogeneous they do not lead to proper thermalization. This may be ascribed to lack of sufficient scattering in the interaction term involving a spatially constant $`C_x`$. We try to improve on this  by writing a homogeneous initial density matrix as a superposition of gaussian pure states. The Hartree approximation is then applied to realizations of this ensemble, which are typically inhomogeneous. This allows for scattering of the modes (‘particles’) via the fluctuating mean field and redistribution over various momenta. To test for thermalization we start out of equilibrium with ‘the modes in their vacuum’ ($`n_0^\alpha =0`$) and compute equal-time- and auto-correlation functions. For the resulting particle distributions, see . Here we concentrate on damping phenomena.

We found funny modulations in auto-correlation functions. Fig. 1 shows such modulations on top of a roughly exponential decay. The correlation function corresponds to the zero momentum mode of the mean field,

$$F_{\mathrm{mf}}(t)=𝑑x\overline{\phi (x,t)\phi (0,0)}_{\text{ }\mathrm{conn}},$$

where the over-bar denotes a time average, taken after waiting a long time for the system to be in approximate equilibrium. As discussed in , this equilibrium is approximately thermal. The function $`F_{\mathrm{mf}}(t)`$ is analogous to the symmetric correlation function of the quantum field theory:

$$F(t)=𝑑x\frac{1}{2}\{\widehat{\phi }(x,t),\widehat{\phi }(0,0)\}_{\text{ }\mathrm{conn}}.$$
(1)
Does $`F(t)`$ also have the modulations?

## 2 Calculation of $`F(t)`$

The function $`F(t)`$ can be expressed in terms of the (zero momentum) spectral function $`\rho (p^0)`$,

$$F(t)=_{\mathrm{}}^{\mathrm{}}\frac{dp^0}{2\pi }e^{ip^0t}\left(\frac{1}{e^{p^0/T}1}+\frac{1}{2}\right)\rho (p^0)$$
(2)
which is determined by the retarded selfenergy $`\mathrm{\Sigma }(p^0)`$,

$$\rho (p^0)=\frac{2\mathrm{I}\mathrm{m}\mathrm{\Sigma }(p^0)}{[m^2p_0^2+\mathrm{Re}\mathrm{\Sigma }(p^0)]^2+[\mathrm{Im}\mathrm{\Sigma }(p^0)]^2}.$$
(3)
The selfenergy can be calculated in perturbation theory. The relevant diagrams are shown in Fig. 2.

The one-loop diagram is present only in the ‘phase of broken symmetry’. It leads to damping for frequencies $`p_0^2>4m^2`$, which are irrelevant for the quasiparticle damping at $`p_0^2=m^2`$. (There is really only a symmetric phase in 1+1 dimensions, but this is due to symmetry restoration by nonperturbative effects which will not obliterate the damping mechanism described by the diagrams.) So we concentrate on the two-loop diagram. This is given by the sum of two terms, $`\mathrm{\Sigma }_1+\mathrm{\Sigma }_2`$. The first has an imaginary part corresponding to $`13`$ processes, requiring $`p_0^2>9m^2`$, so it does not contribute to plasmon damping. The second is given by

$`\mathrm{\Sigma }_2`$ $`=`$ $`{\displaystyle \frac{9\lambda ^2}{16\pi ^2}}{\displaystyle \frac{dp_2dp_3}{E_1E_2E_3}}`$ (4)
$`{\displaystyle \frac{(1+n_1)n_2n_3n_1(1+n_2)(1+n_3)}{p^0+iϵ+E_1E_2E_3}}`$
$`+\left[(p^0+iϵ)(p^0+iϵ)\right],`$
where $`\lambda `$ is the coupling constant (introduced as $`_1=\lambda \phi ^4/4`$ ), and

$`E_1`$ $`=`$ $`\sqrt{m^2+(p_2+p_3)^2},E_{2,3}=\sqrt{m^2+p_{2,3}^2},`$
$`n_i`$ $`=`$ $`[\mathrm{exp}(E_i/T)1]^1,i=1,2,3.`$
Its imaginary part corresponds to $`22`$ processes, which contribute to plasmon damping (the regions near $`p_0=\pm m`$ ).

Now the formula for the thermal plasmon damping rate (at zero momentum) in terms of the retarded selfenergy,

$$\gamma =\mathrm{Im}\mathrm{\Sigma }(m)/2m,$$
(5)
leads to a logarithmically divergent answer. This is a collinear divergence which is absent in more than one space dimension. Inspection shows that the singular part of $`\mathrm{\Sigma }_2`$ is given by the non-relativistic region of the integral in (4). Using polar coordinates $`p_2=p\mathrm{cos}\varphi `$, $`p_3=p\mathrm{sin}\varphi `$ this non-relativistic ($`p<\kappa m`$) contribution is proportional to

$$_0^\kappa p𝑑p_0^{2\pi }\frac{d\varphi }{2\pi }\frac{1}{p^0+iϵm+(p^2/2m)\mathrm{sin}2\varphi }$$

$$=m\left[i\mathrm{ln}\left(\frac{i}{z}+\sqrt{1\frac{1}{z^2}}\right)\right],$$

$$z=\frac{p^0+iϵm}{\kappa ^2/2m},$$

where $`p^0m`$. This function is plotted in Fig. 3.

A natural way out of the difficulty posed by the divergent $`\mathrm{\Sigma }(m)`$ may be to continue the selfenergy analytically into the lower half of its second Riemann sheet, $`p^0mi\gamma `$, and replace (5) by the improved definition

$$m^2(mi\gamma )^2+\mathrm{\Sigma }(mi\gamma )=0.$$
(6)
For weak coupling $`\lambda /m^21`$ we then get the equation

$$\frac{\gamma }{m}=\frac{9\lambda ^2}{16\pi m^4}\frac{e^{m/T}}{\left(e^{m/T}1\right)^2}\left[\mathrm{ln}\frac{m}{\gamma }+c(T)\right],$$
(7)
where the constant $`c`$ has to be determined by matching a numerical evaluation of $`\mathrm{\Sigma }`$ to the form (2) for $`p^0m`$.

We evaluated $`\mathrm{\Sigma }_2`$ in (4) for $`T=m`$ by numerical integration with $`ϵ/m=0.02`$, $`0.01`$ and linear extrapolation $`ϵ0`$. The result is shown in Figs. 4 and 5, together with a matching to the logarithmic singularity, giving $`c0.51`$. For example, Eq. (7) now gives $`\gamma /m=0.061`$, for $`\lambda /m^2=0.4`$.

To see how well this $`\gamma `$ describes the decay of the correlator $`F(t)`$ we evaluated this function directly from (2) and (3). The divergence in $`\mathrm{Im}\mathrm{\Sigma }(p^0)`$ at $`p^0=m`$ leads to a zero in the spectral function $`\rho (p^0)`$. So is there a peak at all in $`\rho (p^0)`$? Fig. 6 shows what happens: the ‘usual’ peak has separated into two twins!

Fig. 7 shows the resulting $`F(t)`$. The effect of the double peak is indeed an oscillating modulation on top of the roughly exponential decay. The decay corresponding to $`\mathrm{exp}(\gamma t)`$, with $`\gamma `$ given by (7), is also indicated in the plot: it does not do a good job in describing the average decay beyond the first interference minimum. The ‘Twin Peaks’ phenomenon implies that the usual definition of damping rate (6) is unreliable in 1+1 dimensions.

Fig. 8 shows the result of a calculation of $`F(t)`$ with parameters taken from the numerical simulation in Fig. 1 with the larger volume. In this case $`ϵ`$ was kept finite, $`ϵ/m=0.005`$, which may be more realistic since one expects anyway the infinity in $`\mathrm{Im}\mathrm{\Sigma }`$ to be smeared out by damping effects in the propagators in Fig. 2. Figs. 1 and 8 are reasonably similar, but how to compare the the average slopes is somewhat ambiguous.

## 3 Outlook

Summarizing, we are encouraged by the similarities in the qualitative features of the numerical and analytical auto-correlation functions. Quantitatively, the damping times are also of the same order of magnitude (105 vs 67 $`m^1`$), but we did not really accurately compute the relevant auto-correlation function yet in our simulations, which will require a lot more numerical effort.

## Acknowledgements

We thank Gert Aarts for useful conversations. This work is supported by FOM/NWO.


text:
# XMM-Newton Detection of a Comptonized Accretion Disc in the Quasar PKS 0558-504

## 1 Introduction

The most prominent feature in the broad-band spectral energy distribution of many active galactic nuclei (AGN) is the ‘big blue bump’ (BBB), which contains a large fraction of the bolometric luminosity. In most broad-line AGN the BBB begins around 1 $`\mu `$m, rises to a peak in the ultraviolet and then declines again into the unobservable EUV (e.g. Elvis et al. el (1994)). The soft X-ray excesses seen in many AGN may be the high-energy tail of the BBB (e.g. Turner & Pounds tutwo (1989); Walter et al. wa (1994); Laor et al. la (1997)). The harder 2–10 keV X-ray spectrum can usually be quite well described by a power-law with photon-index $`\mathrm{\Gamma }1.9`$. Most AGN also show an iron K$`\alpha `$ fluorescence feature around 6–7 keV and/or a hump of emission above 10 keV (e.g. Pounds et al. pou (1990)), widely interpreted in terms of Compton reflection of hard X-rays by optically-thick matter probably in the form of an accretion disc (George & Fabian ge (1991)).

Three basic classes of model have been proposed to explain the BBB. The first class involves thermal emission from an optically thick accretion disc surrounding the central black hole in which the primary heating is by viscous dissipation within the disc itself (e.g. Shields sh (1979)). The disc thermal spectrum may be comptonized by a hot disc corona; a process which changes the observed shape of the continuum, helping to explain soft X-ray emission from discs which are otherwise too cool (e.g. Czerny & Elvis cz (1987)). A second class of models also involves thermal emission from an accretion disc, but with irradiation by a hard X-ray power-law (perhaps arising from magnetic flares) as the primary energy source (e.g. Ross & Fabian ro (1993)). In these models the disc reprocesses the hard X-ray emission into thermal emission which is re-radiated at the local disc temperature. Irradiated discs may have flatter radial temperature profiles, and hence different spectra, compared to discs dominated by viscous dissipation. The final class of model involves optically thin, thermal free-free emission from a hot plasma (Barvainis ba (1993)), although single temperature models of this kind do not fit the observed continuum very well (Siemiginowska et al. si (1995)).

While the spectral form of the optical and UV continuum in AGN are well constrained by existing data, the soft X-ray continuum shape is often quite poorly defined for individual QSOs and therefore places only limited constraints on BBB models. These limitations have been due to the low spectral resolution, signal-to-noise ratio and/or narrow spectral bandpass provided by previous X-ray instruments, thereby making it difficult to quantify absorption/emission features which may contaminate the X-ray spectrum and hence distort the derived continuum shape. The advent of XMM-Newton with its broad bandpass, high throughput and good spectral resolution now provides the opportunity to obtain high-quality X-ray spectra and hence make progress in defining the intrinsic shape of the X-ray continuum in AGN.

In this paper we present XMM-Newton observations of PKS 0558$``$504, a low-redshift quasar (z $`=0.137`$; $`M_\mathrm{V}24.5`$) identified during the HEAO-1 survey. It has quite narrow Balmer lines (FWHM $`=1500`$ km s-1, Remillard et al. re (1986)), which has led to it being routinely described as a Narrow-Line Seyfert 1 (NLS1), although we note that the UV lines, such as Ly$`\alpha `$, are broader (FWHM $`3800`$ km s-1; Kuraszkiewicz et al. ku (2000)). PKS 0558$``$504 is unusual in being one of the few radio-loud NLS1-type galaxies (Schartel et al. sch (1996)). We find that the X-ray spectrum of PKS 0558$``$504 is dominated by a strong soft X-ray excess, which may be the high-energy tail of the BBB. The X-ray spectrum is not contaminated by either cold or warm (ionised) absorption features. Soft X-ray narrow emission-lines are also excluded at a high significance level. Thus, we claim to have detected the featureless X-ray continuum emission from a luminous quasar.

## 2 XMM-Newton Observations

We present observations of PKS 0558$``$504 taken with the European Photon Imaging Camera (EPIC) during orbits 42 and 84 (1–2 March 2000 and 24–25 May 2000) of the calibration and performance verification phase of XMM-Newton. The EPIC instrument is described in Turner et al. (tuone (2001)) and Strueder et al. (str (2001)). During both orbits the EPIC MOS cameras were operated in several different observing modes. Here we include only those observations taken with the medium filter, for which the total EPIC MOS exposure times for orbits 42 and 84 are 85.5 ksecs and 53.7 ksecs respectively. Only the full-frame mode EPIC PN data for orbit 84 were used, for which the total exposure time is 10.2 ksecs.

For operational reasons, the EPIC observations in each orbit were divided into several separate exposures of varying lengths. These exposures were reduced using the XMM-Newton SAS (Science Analysis Software) package. The event lists output from the standard emchain and epchain scripts were further filtered for each EPIC exposure using the SAS xmmselect task. For the MOS exposures only X-ray events corresponding to patterns 0–12 (similar to ASCA event grades 0–4) were selected. Only pattern 0 (single pixel) events were selected for the PN. Known hot or flickering pixels and electronic noise were rejected using the SAS and the low-energy cutoff was set to 200 eV. One of the EPIC MOS-1 images from orbit 42 is shown in Fig. 1.

Source spectra and light curves were extracted from each EPIC exposure using a circular region of diameter 1 arcminute centred on the source position. PKS 0558$``$504 is by far the brightest X-ray source in the 30 arcminute EPIC field-of-view. The EPIC cameras are subject to non X-ray background events due to particles which can vary significantly during an exposure. Background spectra and light curves for use in later data analysis were derived in a similar way to the source data using an offset region close to the source. We note that the observed count-rate for PKS 0558$``$504 implies little effect due to pile-up and no significant spectral differences were found between observations taken in different MOS window modes. The Xspec v11.0 software package was used to calibrate the background-subtracted EPIC spectra using the most recent camera response matrices derived from ground-based and in-orbit data. Before further analysis the EPIC spectra were binned to give a minimum of 20 counts per bin.

Reflection Grating Spectrograph (RGS; den Herder et al. den (2001)) exposures from orbit 84 were also reduced to search for weak/narrow soft X-ray absorption or emission-line features. Three exposures were analysed, two from RGS-1 (each of 25 ksecs) and one from RGS-2 (31 ksecs). The RGS data were first processed using the SAS rgsproc  script. The merged RGS event lists were then fed to the SAS evselect and xmmselect tasks, and both source and background spectra for each exposure were extracted interactively. These spectra were then calibrated using response matrices derived from the individual RGS exposures using the SAS rmfgen task.

The Optical Monitor (OM; Mason et al. ma (2001)) imaging data for orbits 42 and 84 were reduced using standard SAS tasks and calibrated with the instrument response information for each filter.

## 3 Variability

No substantial variability was detected within individual EPIC exposures on timescales less than a few minutes, so the data were rebinned and assembled into orbit light curves by combining the individual light curves for each exposure in order to search for longer duration events. The PKS 0558$``$504 broad-band EPIC MOS light curve for orbit 84, binned into 125 seconds bins, is shown in Fig. 2.

Strong broad-band X-ray flux variability is seen during both orbits, with similar rates of increase and decrease. The largest change in count rate occurs during orbit 84 – an increase of 35% in 3000 seconds. Based on the spectral fits discussed below and using $`H_0=50`$ km s-1 Mpc-1 and $`q_0=0.5`$, this corresponds to a rate of increase in the 0.2–10 keV luminosity of $`\mathrm{\Delta }\mathrm{L}/\mathrm{\Delta }\mathrm{t}=5\times 10^{41}`$ erg s-2. Assuming photon diffusion through a spherical mass of accreting matter in which the opacity is dominated by Thomson scattering, the observed luminosity change implies a minimum efficiency of converting matter into energy of $`\eta =(\mathrm{\Delta }\mathrm{L}/\mathrm{\Delta }\mathrm{t})/(2\times 10^{42})`$ (Guilbert, Fabian & Rees gu (1983)). In the simplest interpretation, the change in luminosity observed during orbit 84 therefore implies $`\eta =0.25`$, an efficiency which exceeds the theoretical maximum for a non-rotating (Schwarzschild) black hole but is consistent with that for a maximally rotating Kerr black hole. Some contribution of beamed X-rays cannot be ruled out, but the shape of the continuum (discussed below) is not consistent with that of a typical beamed source. A tendency for the X-ray continuum to harden slightly when brighter is seen when comparing the data from orbits 42 and 84. The X-ray variability of PKS 0558$``$504 is discussed in more detail in a separate paper (Gliozzi et al. gltwo (2001)).

During a Ginga observation, the PKS 0558$``$504 continuum appeared to vary by 67% in 3 minutes (implying $`\eta 2`$; Remillard et al. rem (1991)). Given the high count-rates and well determined background of the XMM-Newton observations, we could have detected such a flare on timescales as short as a few tens of seconds, but none were observed. ROSAT HRI observations also failed to detect such rapid flaring (Gliozzi et al. gl (2000)), but did show luminosity variability rates (allowing for their use of $`H_0=70`$ km s-1 Mpc-1) which imply accretion efficiencies similar to that derived from the XMM-Newton data. Without confirmation using an imaging detector it is impossible to confirm or deny the association of the Ginga flare with intrinsic variability in PKS 0558$``$504, so we do not consider the Ginga flare any further in this paper.

## 4 Spectral Analysis

Previous analyses of ASCA observations of PKS 0558$``$504 imply an X-ray spectrum which is softer at lower energies. Fixing the absorbing column density at the Galactic value of $`4.4\times 10^{20}`$ cm-2, Gliozzi et al. (gl (2000)) find $`\mathrm{\Gamma }=2.99\pm 0.09`$ over 0.1–2.4 keV from ROSAT data. ASCA data give $`\mathrm{\Gamma }=2.26\pm 0.03`$ over 0.6–10 keV (Vaughan et al. va (1999)), whereas the Ginga data reveal a typical AGN power-law spectrum with $`\mathrm{\Gamma }=1.89_{0.09}^{+0.04}`$ over 2–18 keV (Lawson & Turner law (1997)). Similarly, Comastri (com (2000)) finds $`\mathrm{\Gamma }=1.97\pm 0.04`$ over 2–10 keV from BeppoSAX  data.

As noted above, no large ($`10`$%) changes in spectral shape were observed during each orbit. Thus, to improve the signal-to-noise ratio the MOS data were combined to form a mean spectrum for each orbit before detailed spectral analysis.

The background-subtracted, summed MOS observations provide the highest S/N data for PKS 0558$``$504 to date. These spectra are very poorly fitted by a single power-law ($`\chi _\nu ^2=2.75`$) due to the presence of a large, broad soft X-ray excess which dominates the X-ray emission from PKS 0558$``$504 in the XMM-Newton bandpass. To illustrate this feature, Fig. 3 shows the MOS spectra for each orbit with a power-law fitted to the 4–10 keV band, allowing for Galactic absorption (fit 1 in Table 1 – all the fits are summarised in Table 1). The power-law is fitted above 4 keV as the soft excess appears to extend to $`3`$ keV. The PN spectrum for orbit 84 is consistent with the MOS allowing for the current calibration uncertainties.

No evidence is found for a significant change in X-ray spectral shape between the two orbits, but there is a small decrease in the mean X-ray flux of $`18`$% from orbit 42 to 84. The OM data show that the UV flux also decreased by $`10`$% from orbit 42 to 84. For orbit 42 the 0.2–2.0 keV and 2–10 keV fluxes are $`2.44\times 10^{11}`$ and $`1.01\times 10^{11}`$ erg cm-2 s-1 respectively, corresponding to luminosities of $`3.7\times 10^{45}`$ and $`0.98\times 10^{45}`$ erg s-1.

The soft X-ray excess in PKS 0558$``$504 is far too broad to be fitted by a single temperature blackbody. To parameterise the spectra with a relatively simple model, we fitted a combination of blackbodies and a power-law. At least three blackbodies are required of 80, 175 and 475 eV, to give an acceptable fit ($`\chi _\nu ^2=1.2`$; fit 2 in Table 1). Other than the relative normalisation (consistent with the small overall change in flux), the derived power-law indices and blackbody temperatures required to fit the mean spectra from orbits 42 and 84 are statistically consistent. Therefore, we combined the entire MOS data before proceeding. The multiple blackbody and power-law fit to the combined MOS data is shown in Fig. 4. The high-energy power-law slope in the multi-component model is consistent with the $`\mathrm{\Gamma }=1.9`$ commonly seen in radio-quiet AGN.

No iron fluorescence line is detected from either neutral or ionised material. For the combined MOS data, the upper limits to the EW of a 6.4 keV line are $`<10`$ eV for a narrow line ($`\sigma =10`$ eV) and $`<90`$ eV for a broad line ($`\sigma =300`$ eV). For a narrow ionised line at 6.7 keV, the upper limit is 42 eV. In addition, no evidence is found in the EPIC or RGS spectra for any intrinsic cold or warm absorption. The O vii or O viii edge optical depth limits are $`\tau <0.1`$. The upper limit for an ionised iron edge is $`\tau <0.2`$. We note that an apparent feature in Figs 3, 4 and 6 at $`7`$ keV (which is the wrong energy for an iron line in the QSO rest-frame) is due to noise in some of the MOS spectra that went to form the mean. This feature is not present in the PN spectra. Overall, PKS 0558$``$504 is a powerful quasar in the X-ray band, with a 0.25–10 keV (unabsorbed) luminosity of $`4.5\times 10^{45}`$ erg s-1 mostly emitted below 3 keV. The lack of spectral features strongly implies that we are seeing the bare quasar continuum emission.

## 5 Discussion

The multiple-blackbody plus power-law model discussed above provides an acceptable fit to the EPIC spectra, but has no physical basis. To examine the broad-band spectral energy distribution (SED) of PKS 0558$``$504 the dereddened optical and UV data from the OM (assuming Galactic E(B$``$V) $`=0.044`$) together with the X-ray continuum derived from the combined MOS data (using the multiple blackbody model) are shown in Fig. 5.

The SED of PKS 0558$``$504 clearly shows an optical/UV BBB, and the XMM-Newton data suggest this feature extends into the X-ray band. To power an object as luminous as PKS 0558$``$504 ($`L_{\mathrm{bol}}10^{46}`$ erg s-1) would require a black hole of mass $`M_{\mathrm{BH}}10^8`$ $`M_{\mathrm{}}`$ accreting at the Eddington rate. Based on accretion disc models (e.g. Czerny & Elvis cz (1987)), the thermal spectrum of such a system could explain the optical/UV portion of the BBB but would be unable to produce the broad, soft X-ray excess or the high-energy power-law-like emission.

One model suggested to explain the steep X-ray spectra of extreme soft X-ray sources, such as NLS1s, is that these source accrete at an unusually high rate leading to a large soft photon flux. These photons then Compton cool the hard photon source, generated in a hot accretion disc corona, resulting in a steeper X-ray spectrum (Pounds, Done & Osborne po (1995)). The temperature, geometry and covering factor of disc coronae are poorly defined from theory, but multiple-temperatures are likely to be required due to vertical temperature stratification.

To explore Comptonization models, the results of fitting the Xspec model CompTT to the combined orbit 42 and 84 data are given in Table 1 (fit 3) and shown in Fig. 6. A model involving two electron populations, both scattering photons from a single temperature thermal source, provides a good fit to the mean X-ray spectrum. Allowing the thermal source temperature to vary independently does not produce a significantly better fit. The temperature and luminosity of the soft photon source are consistent with that of a thermal accretion disc spectrum for a $`10^8M_{\mathrm{}}`$ black hole accreting at close to the Eddington rate. The Comptonizing electron temperatures of 4 and 55 keV are at the low end of the expected temperature range (Haardt & Maraschi ha (1993)), possibly due to the large soft photon flux Compton cooling the corona. We note that the exact values of the parameters for the two electron populations in the comptt model are quite poorly constrained. This is because increasing the optical depth slightly can harden the spectrum in a similar way to adding a higher temperature electron population, thus giving a variety of two-component models that adequately fit the data.

We note that, as suggested by the power-law and multiple blackbody fits, the high-energy tail of the Comptonized spectrum has a shape consistent with that of ‘normal’ broad-line AGN power-laws ($`\mathrm{\Gamma }1.9`$), and is significantly harder than that derived from the ASCA data. These differences are presumably due to the inability of the ASCA data to delineate the tail of the soft X-ray excess which extends to around 3 keV in this source. The observed mean slope of the hard power-law in AGN ($`\mathrm{\Gamma }=1.9`$) is predicted from models in which feedback occurs between soft photons produced by hard X-ray heating Compton cool a moderate optical-depth corona (Haardt & Maraschi ha (1993)). Such a model cannot explain the bulk of the soft X-ray luminosity in PKS 0558$``$504 as it dominates the luminosity, but may produce the residual hard power-law emission.

The reduced chi-squared value for the Comptonized model ($`\chi _\nu ^2=1.3`$) shows that the Comptonization model fit leaves some spectral curvature, mainly above a few keV. Some of this imperfection is due to remaining calibration uncertainties around features such as detector Gold edge, but most of it is due to the Comptonization model simply being too simplistic to adequately fit such high-quality data. Some of the remaining spectral curvature can be reduced by including an additional Comptonizing component (i.e. hotter electrons) but this component does not formally improve the fit in terms of the reduced chi-squared presumably because the XMM-Newton high-energy cutoff is too low to adequately constrain a very-hot electron population.

An alternative explanation for a high-energy ‘power-law-like’ component is a small amount of reflection from an ionised disc. The disc is likely to be ionised given the luminosity of the source. Indeed, the absence of any strong iron K$`\alpha `$ emission line or iron edge at high energies and the lack of absorption features due to moderately ionised material at low energies suggests that any X-ray heated disc will have a substantial ‘skin’ which is fully ionised and effectively acting as a ‘mirror’. This has been suggested for high accretion rate sources in the context of a magnetic-flare model for X-ray variability (Nayakshin na (2000)). The small iron K$`\alpha `$ emission line and iron edge strengths predicted by the magnetic-flare model, assuming $`\mathrm{\Gamma }2`$, are consistent with the upper limits from the XMM-Newton spectra (Nayakshin, private communication). What is less clear is if a fully ionised corona can be maintained in a source like PKS 0558$``$504 with a large soft photon flux if those photons are generated close to the magnetic-flares. Clearly the XMM-Newton data for PKS 0558$``$504 demand comparison with more detailed accretion disc models.

## 6 Conclusions

The XMM-Newton X-ray spectrum of PKS 0558$``$504 is dominated by a large, variable soft excess. The fastest observed variability implies a minimum efficiency of converting matter to energy of $`\eta =0.25`$, suggesting accretion onto a Kerr black hole. Combining the X-ray data with simultaneous OM optical and UV data reveal a large big blue bump; too broad to be due purely to a thermal accretion disc spectrum. The X-ray spectrum can be fitted by a Comptonization model in which thermal photons, presumably from the disc, are scattered by a hot corona. At least two electron temperatures are required. Overall, PKS 0558$``$504 can be explained by an accretion disc running at a high accretion rate producing copious amounts of soft EUV photons which are then Compton cooling a hot disc corona.

## Acknowledgements

Based on observations obtained with the XMM-Newton, and ESA science mission with instruments and contributions diectly funded by ESA member states and the USA (NASA). EPIC was developed by the EPIC Consortium led by the Principal Investigator, Dr. M. J. L. Turner. The consortium comprises the following Institutes: University of Leicester, University of Birmingham, (UK); CEA/Saclay, IAS Orsay, CESR Toulouse, (France); IAAP Tuebingen, MPE Garching, (Germany); IFC Milan, ITESRE Bologna, IAUP Palermo, (Italy). EPIC is funded by: PPARC, CEA, CNES, DLR and ASI. We thank the EPIC and SSC teams for much help and advice, particularly Steve Sembay, Gareth Griffiths, Richard West and Ian Stewart. We also thank Kim Page for help with data reduction.