text:
# QUARKONIUM FORMATION FROM UNCORRELATED QUARK-ANTIQUARK PAIRS

## 1 INTRODUCTION

The utility of heavy quarkonium production rates in nuclear collisions as a signature of color deconfinement was proposed more than 15 years ago . Since one expects that the long-range color confining potential will be screened in a deconfined medium, the quark and antiquark constituents of bound states will be liberated. As the system expands and cools, these constituents will, in general, diffuse away from each other to separations larger than typical hadronic dimensions. When the confining potential reappears, a given heavy quark will not be able to “find” its heavy antiquark partner and form heavy quarkonium. It must then bind with one of the antiquarks within range at hadronization. Since these antiquarks are predominantly the lighter $`u`$, $`d`$, and $`s`$ flavors, the final hadronic states will preferentially be those with “open” heavy flavor. The result will be a decreased population of heavy quarkonium relative to that which would have formed if a region of deconfinement had not been present. This scenario as applied to the charm sector is known as $`J/\psi `$ suppression.

At LHC energy, perturbative QCD estimates predict that hundreds of pairs of charm-anticharm quarks will be produced in a central lead-lead collision. This situation provides a “loophole” in the Matsui-Satz argument since there will be copious numbers of heavy antiquarks in the interaction region with which any given heavy quark may combine. In order for this to happen, however, one must invoke a physical situation in which quarkonium states can be formed from all combinations of heavy quarks and antiquarks. This of course would be expected to be valid in the case that a space-time region of color deconfinement is present but it is not necessarily limited to this possibility.

One can make a model-independent estimate of how such a “recombination” mechanism would depend on nuclear collision observables. For a given charm quark, the probability $`𝒫`$ to form a $`J/\psi `$ is proportional to the number of available anticharm quarks relative to the number of light antiquarks,

$$𝒫\frac{N_{\overline{c}}}{N_{\overline{u},\overline{d},\overline{s}}}\frac{N_{c\overline{c}}}{N_{\mathrm{ch}}}.$$
(1)
In the second step, we have replaced the number of available anticharm quarks by the total number of pairs initially produced, assuming that the total number of bound states formed remains a small fraction of the total $`c\overline{c}`$ production. We normalize the number of light antiquarks by the number of produced charged hadrons. Since this probability is generally very small, one can simply multiply by the total number of charm quarks, $`N_c`$, to obtain the number of $`J/\psi `$ expected in a given event,

$$N_{J/\psi }\frac{N_{c\overline{c}}^{}{}_{}{}^{2}}{N_{\mathrm{ch}}},$$
(2)
where the use of the initial values $`N_{c\overline{c}}=N_c=N_{\overline{c}}`$ is again justified by the relatively small number of bound states formed.

The essential property of this result is that the growth of $`N_{J/\psi }`$, quadratic in total charm, with energy  is expected to be much faster than the growth of total particle production in heavy ion collisions . Without this quadratic mechanism, $`J/\psi `$ production is typically some small energy-independent fraction of total initial charm production . We thus anticipate that the quadratic formation will become dominant at sufficiently high energy. Generic estimates of the significance of this type of formation process can be made . Here we look at specific predictions of two models which share the above properties and update the expectations to LHC energies.

## 2 Statistical Hadronization

The statistical hadronization model is motivated by the successful fits of relative abundances of light hadrons produced in high energy heavy ion interactions according to a hadron gas in chemical and thermal equilibrium . Extension of the model to hadrons containing heavy quarks underpredicts the observed abundances. This effect may be attributed to the long time scales associated with thermal production and annihilation of heavy quarks. The statistical hadronization model as first formulated for charm quarks  assumes that the $`c\overline{c}`$ pairs produced in the initial hadronic interactions survive until their subsequent hadronization, at which time they are distributed into hadrons according to the same thermal equilibrium parameters that fit the light hadron abundances. Chemical equilibrium abundances are adjusted by a factor $`\gamma _c`$ which accounts for the non-thermal heavy quark density. One power of this factor multiplies a given thermal hadron population for each heavy quark or antiquark contained in the hadron. Thus the relative abundance of the $`J/\psi `$ to that of $`D`$ mesons, for example, may be enhanced in this model.

The value of $`\gamma _c`$ is determined by conservation of the heavy quark flavor. For the charm sector, the conservation constraint relates the number of initially-produced c-cbar pairs $`N_{c\overline{c}}`$ to their distribution into open and hidden charm hadrons,

$$N_{c\overline{c}}=\frac{1}{2}\gamma _cN_{\mathrm{open}}+\gamma _{c}^{}{}_{}{}^{2}N_{\mathrm{hidden}},$$
(3)
where $`N_{\mathrm{open}}`$ is the number of hadrons containing one $`c`$ or $`\overline{c}`$ quark and $`N_{\mathrm{hidden}}`$ is the number of hadrons containing a $`c\overline{c}`$ pair. For most applications, $`N_{\mathrm{hidden}}`$ (and also multi-charm hadrons) can be neglected compared with $`N_{\mathrm{open}}`$ due to the mass differences. Thus the charm enhancement factor is simply

$$\gamma _c=\frac{2N_{c\overline{c}}}{N_{\mathrm{open}}},$$
(4)
leading directly to the quadratic dependence of the hidden charm hadron population on $`N_{c\overline{c}}`$. One can then express the total number of $`J/\psi `$ in terms of the various thermal densities, $`n_i`$, and the total number of $`c\overline{c}`$ pairs, $`N_{c\overline{c}}`$. One factor of system volume $`V`$ remains implicit here. It is generally replaced by the ratio of number to density for total charged hadrons, $`n_{\mathrm{ch}}/N_{\mathrm{ch}}`$. Then the number of $`J/\psi `$ produced obeys the generic form anticipated in Eq. (2).

$$N_{J/\psi }=4\frac{n_{\mathrm{ch}}n_{J/\psi }}{n_{\mathrm{open}}^{}{}_{}{}^{2}}\frac{N_{c\overline{c}}^2}{N_{\mathrm{ch}}}$$
(5)
For collider experiments such as those at the LHC and RHIC, relating the corresponding central rapidity densities will be more relevant. Since Eq. (5) is homogeneous in the total particle and quark pair numbers, it will also be valid if these are replaced by their rapidity densities. To get an order of magnitude estimate, we choose a “standard” set of thermal parameters, $`T=170`$ MeV and $`\mu _B0`$, for which the thermal density ratio is approximately $`0.5`$. For a specific normalization, we assume $`dN_{\mathrm{ch}}/dy`$ = 2000 for a central collision at the LHC and take the initial charm rapidity density to be $`dN_{c\overline{c}}/dy=25`$, roughly corresponding to $`N_{c\overline{c}}=200`$ for central collisions $`(b=0)`$. Using these inputs, one predicts $`dN_{J/\psi }/dy=0.625`$, indicating that several $`J/\psi `$ will form through statistical hadronization in a central collision. To put this number in perspective, it is revealing to form the $`J/\psi `$ to $`N_{c\overline{c}}`$ rapidity density ratio, 0.025 with the same assumptions. For comparison, one expects the corresponding hadronic production ratio to be of $`0.01`$. This number would then be significantly reduced if placed in a region of color deconfinement. Thus the efficiency of $`J/\psi `$ formation via statistical hadronization at the LHC is expected to be substantial.

These numbers can be easily adjusted to other charm and charged particle densities using Eq. (5). Variations of the thermal parameters can also be investigated. For example, if the hadronization temperature is decreased to 150 MeV, the prefactor combination of thermal densities increases by approximately a factor of two.

The centrality dependence is controlled by the behavior of $`N_{c\overline{c}}`$ and $`N_{\mathrm{ch}}`$. The former should be proportional to the nuclear overlap function $`T_{AA}(b)`$ but is generally recast in terms of the dependence on the number of nucleon participants, $`N_{\mathrm{part}}`$. The calculation of $`N_{\mathrm{part}}`$ requires a model calculation dependent on the total inelastic cross section, $`\sigma _{\mathrm{in}}`$, as well as the nuclear geometry. We parameterize the expected behavior as a power-law $`N_{c\overline{c}}N_{\mathrm{part}}^{}{}_{}{}^{4/3}`$. However, there will be deviations from this behavior for the larger values of $`\sigma _{\mathrm{in}}`$ expected at the LHC . The centrality dependence of $`N_{\mathrm{ch}}`$ at RHIC is also consistent with a power-law with exponent $``$ 1.2 . We will use the same dependence for our estimates at the LHC. It is clear that for sufficiently peripheral collisions one will encounter situations in which the average number of initially produced $`c\overline{c}`$ pairs is of order unity or less. At this point, one must revisit the assumptions of the original statistical hadronization model which assumed a grand canonical ensemble. The grand canonical approach is valid only when $`N_{c\overline{c}}`$ is large enough for the fluctuations about the average value to be neglibile. Thus for peripheral collisions, one must recalculate the statistical results in the canonical approach where the charm number is exactly conserved, as noted in . Charm conservation can be implemented via a correction factor ,

$$N_{c\overline{c}}=\frac{1}{2}\gamma _cN_{\mathrm{open}}\frac{I_1(\gamma _cN_{\mathrm{open}})}{I_0(\gamma _cN_{\mathrm{open}})}+\gamma _{c}^{}{}_{}{}^{2}N_{\mathrm{hidden}}.$$
(6)
In the limit of large $`\gamma _cN_{\mathrm{open}}`$, the ratio of Bessel functions approaches unity and the grand canonical result is recovered. In the opposite limit when $`\gamma _cN_{\mathrm{open}}0`$, the ratio of Bessel functions approaches $`\frac{1}{2}\gamma _cN_{\mathrm{open}}`$. In this limit, the dependence on $`N_{c\overline{c}}`$ in Eq. (5) changes from quadratic to linear. At the LHC this effect will not be relevant until one reaches very peripheral events, but at lower energies it can be significant over a much larger range of centralities .

The results for $`dN_{J/\psi }/dy`$ as a function of $`N_{\mathrm{part}}`$ at the LHC are shown in Fig. 1. The results are shown for three different values of $`dN_{c\overline{c}}/dy(b=0)`$, corresponding to $`N_{c\overline{c}}(0)200`$, 150, and 100. There is a rapid increase with centrality due to the quadratic dependence of $`N_{J/\psi }`$ on $`N_{c\overline{c}}`$.

It is also interesting to look at these results normalized by $`N_{\mathrm{part}}`$, shown in Fig. 2. This ratio also increases with centrality, providing a signature for the statistical hadronization process that is less dependent on $`N_{c\overline{c}}`$ for the overall normalization. The corresponding results when normalized by $`dN_{c\overline{c}}/dy`$ are shown in Fig. 3. The same general behavior is seen but the increase with centrality is less pronounced since the $`dN_{c\overline{c}}/dy`$ is assumed to vary with a larger power, $`N_{\mathrm{part}}^{}{}_{}{}^{4/3}`$. All of these ratios are at the percent level for central collisions and hence are larger than expected if the total $`J/\psi `$ population were due to initial production followed by suppression in a deconfined medium.

The region of very peripheral collisions deserves some separate comments. First, there is a rise at low $`N_{\mathrm{part}}`$ in both Figs. 2 and 3 due to the onset of corrections from the canonical ensemble treatment. However, the extremely large values of the ratios as $`N_{\mathrm{part}}0`$ are an artifact of the decreasing interaction volume, $`V0`$. This calculation must be cut off before $`N_{\mathrm{part}}=2`$, i.e. only one interacting pair. It is also in this region where one must take into account the remaining $`J/\psi `$ from initial production. Since the survival probability is maximum for very peripheral collisions and the statistical hadronization process is least effective in this same region, there will be a crossover in the relative importance of the two mechanisms. Some studies have already been performed for this situation at SPS and RHIC energies .

Finally, there is another lower cutoff in centrality for the statistical hadronization results, needed to avoid a contradiction with the $`\psi /(J/\psi )`$ ratio at the SPS. Since both of these states receive identical factors of $`\gamma _c`$, their ratio must be that predicted for chemical equilibrium in the absence of any charm enhancement or suppression. Although the measured ratio appears to be consistent for more central collisions , there is an indication that it rises sharply for more peripheral collisions. Most treatments have thus inserted a cutoff of $`N_{\mathrm{part}}`$ = 100, below which model predictions become inconsistent .

The numerical values for $`dN_{J/\psi }/dy`$ are tabulated as a function of impact parameter in Table 1 for the three choices of initial charm multiplicity density and the default values of all other quantities.

## 3 Kinetic Formation in a Deconfied Region

The kinetic model has been developed  to investigate the possibility that $`J/\psi `$ may form directly in a deconfined medium. This formation takes advantage of the mobility of the initially-produced charm quarks in a deconfined region. In order to motivate this view, consider the “standard” physical picture of deconfinement in which quarkonium is suppressed by collisions with free gluons in the medium . Then the formation process, in which a $`c`$ and $`\overline{c}`$ in a relative color octet state are captured into a color-singlet quarkonium bound state and emit a color octet gluon, is simply the inverse of the the breakup reaction responsible for the suppression. This is an inevitable consequence of the suppression picture.

The proper time evolution of the $`J/\psi `$ population is given by the rate equation

$$\frac{dN_{J/\psi }}{d\tau }=\lambda _\mathrm{F}\frac{N_cN_{\overline{c}}}{V(\tau )}\lambda _\mathrm{D}N_{J/\psi }\rho _g,$$
(7)
where $`\rho _g`$ is the gluon number density and $`V(\tau )`$ is the time-dependent volume of the deconfined spatial region. The reactivities $`\lambda _{\mathrm{F},\mathrm{D}}`$ are the reaction rates, $`\sigma v_{\mathrm{rel}}`$, averaged over the momentum distributions of the initial participants, i.e. $`c`$ and $`\overline{c}`$ for $`\lambda _\mathrm{F}`$ and $`J/\psi `$ and $`g`$ for $`\lambda _\mathrm{D}`$.

The solution of Eq. (7) grows quadratically with $`N_{c\overline{c}}`$, as long as $`N_{J/\psi }N_{c\overline{c}}`$. In this case, we have

$$N_{J/\psi }(\tau _f)=ϵ(\tau _f)\left[N_{J/\psi }(\tau _0)+N_{c\overline{c}}^2_{\tau _0}^{\tau _f}𝑑\tau \lambda _\mathrm{F}[V(\tau )ϵ(\tau )]^1\right].$$
(8)
The function $`ϵ(\tau _f)=\mathrm{exp}(_{\tau _0}^{\tau _f}𝑑\tau \lambda _\mathrm{D}\rho _g)`$ would be the suppression factor if formation were neglected.

The quadratic factor $`N_{c\overline{c}}^2`$ is present, as expected, for the additional formation process. The normalization factor of $`N_{\mathrm{ch}}`$ is not immediately evident, but is implicit in the system volume factor. This volume is now time-dependent, accounting for the decreasing charm quark density during expansion. Here the transverse area of the deconfined region is determined not just by the nuclear geometry but by the dynamics which determine the extent of the deconfined region. This area is modeled by the energy density in terms of the local participant density in the transverse plane, $`n_{\mathrm{part}}(b,s=0)`$. The transverse area is defined by the ratio of the participant number to the local participant density. Note that the maximum local density is at $`s=0`$. Thus,

$$A_T(b)=A_T(0)[N_{\mathrm{part}}(b)n_{\mathrm{part}}(0,s=0)/N_{\mathrm{part}}(0)n_{\mathrm{part}}(b,s=0)]$$
(9)
These area effects will be more explicit when the centrality dependence is considered.

The numerical results depend on a number of parameters, including the initial volume and temperature, the time expansion profile, the reaction cross sections, the behavior of the quarkonium masses and binding energies in the deconfined region, and the charm quark momentum distributions. For specifics, see Ref. . Our previous results have used initial values $`N_{c\overline{c}}=200`$, 150, and 100, spanning a reasonable range of expectations . The results are very sensitive to the initial charm quark momentum distributions, as may be expected. We assume the charm $`p_T`$ distributions are gaussian and the charm rapidity distributions are flat over a plateau of variable width, $`\mathrm{\Delta }y`$. The range $`1<\mathrm{\Delta }y<7`$ spans the range between an approximate thermal momentum distribution, $`\mathrm{\Delta }y1`$, to a distribution similar to that of the initial pQCD production, $`\mathrm{\Delta }y7`$.

The results as a function of the initial number of $`c\overline{c}`$ pairs produced in central collisions are shown in Fig. 4. There is a rapid decrease in formation with increasing $`\mathrm{\Delta }y`$. The quadratic dependence on $`N_{c\overline{c}}`$ is evident, but there is also a substantial linear component in some of the curves. This linear contribution arises because the final $`J/\psi `$ formation by this mechanism is large enough for exact charm conservation to reduce the number of $`c`$ and $`\overline{c}`$ quarks available to participate in the formation process. The curve labeled “Quadratic Extrapolation” uses a quadratic dependence derived from a fit valid only for low $`N_{c\overline{c}}`$. Note that the result for a thermal distribution is very similar to the assumption $`\mathrm{\Delta }y=1`$.

The corresponding centrality dependence is presented in Fig. 5, where we give $`N_{J/\psi }`$ at hadronization for three different initial charm quark momentum distributions, thermal, $`\mathrm{\Delta }y=4`$ and $`\mathrm{\Delta }y=7`$, as well as for our three choices of $`N_{c\overline{c}}(b=0)`$.

Finally, the ratio of final $`J/\psi `$ to initial charm production is shown in Fig. 6 using the same parameters in Fig. 4. These ratios are most easily compared to either initial production or suppression. There is a substantial variation in the predictions and it is evident that a simultaneous measurement of open charm will be required for an interpretation. However, the centrality dependence is opposite to that expected in any pure suppression scenario.

We have updated the calculations to include the charm quark momentum distribution from a leading order pQCD calculation . The rapidity distribution has a somewhat larger effective $`\mathrm{\Delta }y`$ and the $`p_T`$ distribution does not fall as fast as a simple gaussian. As a result, the formation efficiency is further reduced. Such distributions may be most relevant, given preliminary results from RHIC .

The numerical values for $`N_{J/\psi }`$ are compared with the statistical hadronization model results for $`dN_{J/\psi }/dy`$ in Table 1. The overall magnitudes are comparable, although the centrality dependences differ somewhat. Thus details such as the resulting $`J/\psi `$ momentum distributions will be required to differentiate between these two models . For completeness, $`N_{J/\psi }`$ for the thermal distributions and the assumption $`\mathrm{\Delta }y=4`$ are presented in Table 2.

## 4 Conclusions

The ”smoking gun” signature of the quarkonium formation mechanism is the quadratic dependence on total charm. For central collisions at the LHC one expects that this feature will lead to a total $`J/\psi `$ rate greater than that produced by an incoherent superposition of the initial nucleon-nucleon collisions, even without any subsequent suppression due to deconfinement effects. In addition, the centrality dependence can be used to identify the quadratic dependence on charm assuming that the initial charm production scales with the number of binary collisions. Binary scaling leads to an increase of the ratio of $`J/\psi `$ to initial charm as the collision centrality increases, independent of specific parameters which control the overall magnitudes. A simultaneous measurement of total charm will be essential for such conclusions to be drawn.

Uncertainties in the absolute magnitude of the formation process are inherent in the model parameters. For statistical hadronization, one can constrain the thermal parameters to within a factor of two using the observed hadron populations. There is some additional uncertainty related to the lower cutoff on centrality needed to ensure the quarkonium ratios are consistent with an overall thermal picture. There is also the possibility that the correction for canonical ensemble effects will involve a thermal volume parameter not necessarily equal to the total system volume . In addition, the formation mechanism could be limited to those charm quarks whose phase space separation is within some maximum value, introducing another as yet unconstrained parameter . With kinetic formation, a similar set of uncertainties exist. There are uncertainties in the space-time properties of the deconfinement region. In addition, possible variations of charmonium binding energies and reaction cross sections in a deconfined region are at present not well understood. There are indications that the efficiency of the formation mechanism is considerably reduced when included in a partonic transport calculation .

The primary uncertainty in both models is still the initial number of charm quarks and their momentum distributions. The tabulated $`J/\psi `$ results should be regarded in this light. Thus numbers may be only an order of magnitude estimate. However, the variation with centrality and total initial charm should provide experimental signatures which are largely independent of the overall magnitudes.

## 5 Acknowledgments

My thanks to Anton Andronic for discussions on the Statistical Model and Martin Schroedter for updates on the Kinetic Model calculations. This work was supported in part by U.S. Department of Energy Grant DE-FG03-95ER40937.


format: text

text:
# Jet Tomography in Heavy Ion Collisions

## 1 Introduction

Experiments at the Relativistic Heavy Ion Collider RHIC in Brookhaven show for the first time a significant quenching of high-$`p_{}`$ leading hadron spectra. In particular, transverse momentum spectra for neutral pion  and charged hadron  are suppressed if compared to spectra in p+p collisions rescaled by the number of binary collisions. This suppression is most pronounced (up to a factor $`5`$) in central Au+Au collisions and smoothly approaches the binary scaling case with decreasing centrality. The azimuthal anisotropy $`v_2(p_{})`$ of hadroproduction stays close to maximal up to the highest transverse momentum . Moreover, the disappearance of back-to-back high-$`p_{}`$ hadron correlations  provides an additional indication that final state medium effects play a decisive role in hadroproduction up to $`p_{}10`$ GeV.

Parton energy loss has been proposed to account for the small nuclear modification factor , the azimuthal anisotropy  and the disappearance of dijets . Several studies (see references given in ) indicate, however, that in the kinematical regime relevant for RHIC ($`p_{}<12`$ GeV), $`p_{}`$-broadening, shadowing, formation time and possibly other effects contribute significantly to the high-$`p_{}`$ nuclear modification as well.

Here, we discuss the status of parton energy loss calculations and their comparison to data in the kinematical $`p_{}`$-range probed at RHIC.

## 2 Medium-induced gluon radiation from a static medium

Several groups  calculated recently the modification of the elementary splitting processes $`qqg`$ and $`ggg`$ due to multiple scattering. The inclusive energy distribution of gluon radiation off an in-medium produced parton takes the form

$`\omega {\displaystyle \frac{dI}{d\omega }}`$ $`=`$ $`{\displaystyle \frac{\alpha _sC_R}{(2\pi )^2\omega ^2}}\mathrm{\hspace{0.17em}2}\mathrm{Re}{\displaystyle _{\xi _0}^{\mathrm{}}}𝑑y_l{\displaystyle _{y_l}^{\mathrm{}}}𝑑\overline{y}_l{\displaystyle d^2𝐮_0^{\chi \omega }d^2𝐤e^{i𝐤_{}𝐮}e^{\frac{1}{2}_{\overline{y}_l}^{\mathrm{}}𝑑\xi n(\xi )\sigma (𝐮)}}`$ (1)
$`\times {\displaystyle \frac{}{𝐲}}{\displaystyle \frac{}{𝐮}}{\displaystyle _{𝐲=0}^{𝐮=𝐫(\overline{y}_l)}}D𝐫\mathrm{exp}\left[i{\displaystyle _{y_l}^{\overline{y}_l}}𝑑\xi {\displaystyle \frac{\omega }{2}}\left(\dot{𝐫}^2{\displaystyle \frac{n(\xi )\sigma \left(𝐫\right)}{i\omega }}\right)\right],`$
where gluons of transverse momentum $`k_{}<\chi \omega `$ are included. The radiation of hard quarks or gluons differs by the Casimir factor $`C_R=C_F`$ or $`C_A`$, respectively. The properties of the medium enter eq. (1) by the product of the time-dependent density $`n(\xi )`$ of scattering centers times the strength of a single elastic scattering $`\sigma (𝐫)=2\frac{d𝐪}{(2\pi )^2}|a(𝐪)|^2\left(1e^{i𝐪𝐫}\right)`$, where $`|a(𝐪)|^2`$ denotes the elastic high-energy cross section of a single scatterer.

For explicit calculations, one has to approximate the path integral in (1). This is done either by a saddle point approximation  obtained in the medium-induced soft multiple scattering approximation,

$`n(\xi )\sigma (𝐫){\displaystyle \frac{1}{2}}\widehat{q}(\xi )𝐫^2.`$ (2)
Here, the only medium-dependent quantity is the transport coefficient  $`\widehat{q}(\xi )`$ which characterizes the transverse momentum squared $`\mu ^2`$ transferred to the projectile per mean free path $`\lambda `$. Alternatively, one can proceed in the opacity expansion which amounts to expanding the path integral in (1) in powers of the elastic scattering center

$`K(𝐫,y_l;\overline{𝐫},\overline{y}_l)=D𝐫\mathrm{exp}\left[i{\displaystyle _{y_l}^{\overline{y}_l}}𝑑\xi {\displaystyle \frac{\omega }{2}}\left(\dot{𝐫}^2{\displaystyle \frac{n(\xi )\sigma \left(𝐫\right)}{i\omega }}\right)\right]`$ (3)
$`=K_0(𝐫,y_l;\overline{𝐫},\overline{y}_l){\displaystyle \underset{z}{\overset{z^{}}{}}}𝑑\xi n(\xi ){\displaystyle 𝑑\rho K_0(𝐫,y_l;\rho ,\xi )\frac{\sigma (\rho )}{2}K_0(\rho ,\xi ;\overline{𝐫},\overline{y}_l)}+\mathrm{}.`$
Fig.1 compares numerical results for the gluon energy distribution obtained in these two approximations. Qualitatively, this figure can be understood by estimating to what degree initial state gluons decohere from the partonic projectile. Decoherence depends on the relative phase $`\phi `$ accumulated by the gluon due to scattering. This phase grows with the transverse momentum accumulated by the emitted gluon.

For multiple soft scattering, we have

$$\phi =\frac{k_{}^2}{2\omega }\mathrm{\Delta }z\frac{\widehat{q}L}{2\omega }L=\frac{\omega _c}{\omega },$$
(4)
which defines the “characteristic gluon frequency” $`\omega _c=\frac{1}{2}\widehat{q}L^2`$. Given the number $`N_{\mathrm{coh}}`$ of scattering centers which add coherently in the gluon phase (4), one find $`k_{}^2N_{\mathrm{coh}}\mu ^2`$. With the coherence time of the emitted gluon, $`t_{\mathrm{coh}}\frac{\omega }{k_{}^2}\sqrt{\frac{\omega }{\widehat{q}}}`$ and $`N_{\mathrm{coh}}=\frac{t_{\mathrm{coh}}}{\lambda }=\sqrt{\frac{\omega }{\mu ^2\lambda }}`$, one estimates for the gluon energy spectrum per unit pathlength

$$\omega \frac{dI^{(mult)}}{d\omega dz}\frac{1}{N_{\mathrm{coh}}}\omega \frac{dI^{1\mathrm{scatt}}}{d\omega dz}\frac{\alpha _s}{t_{\mathrm{coh}}}\alpha _s\sqrt{\frac{\widehat{q}}{\omega }}.$$
(5)
This $`1/\sqrt{\omega }`$-energy dependence agrees with the small-$`\omega `$ behaviour in Fig.1. It is cut off for $`\omega >\omega _c`$ where the phase (4) is smaller than unity and the reduced decoherence suppresses gluon emission.

For single hard scattering ($`N=1`$ opacity expansion) with momentum transfer $`\mu `$, decoherence occurs if the typical gluon formation time $`\overline{t}_{\mathrm{coh}}=\frac{2\omega }{\mu ^2}`$ is smaller than the typical distance $`L`$ between the production point of the parton and the position of the scatterer. The relevant phase is

$$\gamma =\frac{L}{\overline{t}_{\mathrm{coh}}}\frac{\overline{\omega }_c}{\omega },$$
(6)
which defines the characteristic gluon energy $`\overline{\omega }_c=\frac{1}{2}\mu ^2L`$. The gluon energy spectrum per unit pathlength can be estimated in terms of the coherence time $`\overline{t}_{\mathrm{coh}}`$,

$$\omega \frac{dI^{N=1}}{d\omega dz}\frac{\alpha _s}{\overline{t}_{\mathrm{coh}}}\alpha _s\frac{\mu ^2}{\omega }.$$
(7)
The full calculation in Fig.1 agrees with this $`1/\omega `$-dependence in the range $`\omega >\overline{\omega }_c`$.

In QCD, collinear gluons are hard and softer gluons tend to be emitted under larger angles. As a consequence, the limitations on transverse momentum phase space translate into a depletion of the infrared region of the gluon energy distribution, seen in Fig.1. The estimates (5) and (7) determine the dependence on the characteristic gluon energies $`\omega _c`$ ($`\overline{\omega }_c`$) only. They do not include the constraint on transverse phase space. The latter is determined by the parameters

$$R_\chi =\frac{1}{2}\widehat{q}\chi ^2L^3=\chi ^2\omega _cL,\overline{R}_\chi =\frac{1}{2}\chi ^2\mu ^2L^2=\chi ^2\overline{\omega }_cL.$$
(8)
Clearly, keeping $`\omega _c`$ ($`\overline{\omega }_c`$) fixed and taking $`R`$ ($`\overline{R}`$) to infinity amounts to the limit of infinite pathlength. In this limit, the projectile can accumulate an arbitrarily large medium-induced transverse momentum and the limit on transverse momentum phase space is removed. For a medium of finite size, however, gluon emission at angles $`\mathrm{\Theta }_c^2\frac{k_{}^2_{\mathrm{med}}}{\omega ^2}`$ is suppressed since the emitted gluons are sensitive to the kinematical constraint $`k_{}O(\omega )`$. In the multiple soft scattering approximation, this translates into

$$\mathrm{\Theta }_c^2\frac{\sqrt{\omega \widehat{q}}}{\omega ^2}\left(\frac{\omega }{\omega _c}\right)^{3/2}\frac{1}{R}1\frac{\widehat{\omega }}{\omega _c}\left(\frac{1}{R}\right)^{2/3}.$$
(9)
The position of the maximum of $`\omega \frac{dI^{(mult)}}{d\omega }`$ as a function of $`R`$ is consistent with this dependence on $`\widehat{\omega }`$, see Fig. 1.

For the single hard scattering approximation, the corresponding estimate for the infrared cut-off $`\widehat{\omega }`$ due to transverse momentum phase space constraints reads

$$\mathrm{\Theta }_c^2\frac{\mu ^2}{\widehat{\omega }^2}\left(\frac{\overline{\omega }_c}{\widehat{\omega }}\right)^2\frac{1}{\overline{R}}1\frac{\widehat{\omega }}{\overline{\omega }_c}\frac{1}{\sqrt{\overline{R}}}.$$
(10)
The position of the maximum of $`\omega \frac{dI^{N=1}}{d\omega }`$ in Fig. 1 changes $`\frac{1}{\sqrt{\overline{R}}}`$, in accordance with this estimate. We thus have a semi-quantitative understanding of how phase space constraints deplete the non-perturbative soft region of the medium-induced gluon energy distribution. This suppression of the non-perturbative small-$`\omega `$ contributions helps to make the calculation of medium-induced energy loss perturbatively stable.

In Ref., we compare the single hard and multiple soft scattering approximations of (1) in detail. In general, one finds that the gluon energy distribution is significantly harder in the single hard scattering approximation. For example, the average energy loss $`\mathrm{\Delta }E=𝑑\omega \omega \frac{dI}{d\omega }`$ receives in the multiple soft scattering approximation a dominant contribution from the region $`\omega <\omega _c`$. In contrast, the dominant contribution comes from the hard region $`\omega >\overline{\omega }_c`$ in the single hard scattering approximation.

Remarkably, although both approximations emphasize different kinematical regions, the observable results are quantitatively comparable if comparable sets of model parameters are used. To relate the model parameters in both approximations, one observes that up to logarithmic accuracy $`\mu ^2n_0L\widehat{q}L`$, where $`n_0L`$ determines the average number of scattering centers within the in-medium pathlength $`L`$. This implies

$$R(n_0L)\overline{R},\omega _c(n_0L)\overline{\omega }_c.$$
(11)
In the multiple soft scattering approximation, one uses the transport coefficient $`\widehat{q}`$ to characterize the average transverse momentum squared per unit pathlength. In the opacity expansion, one specifies not only this information, but additionally the average number of scattering centers in which this momentum transfer takes place. As a consequence, the opacity expansion has one additional model parameter $`n_0L`$. Numerically, we find that for the choice $`n_0L=3`$ in (11), the energy density distribution, quenching weight and even the angular dependence of the average energy loss are quantitatively comparable . We found deviations from this relation for those quantities for which due to kinematical constraints only the region $`\omega <\omega _c`$ contributes in both approximations, see below. This is in particular the case for the comparison to RHIC data in Fig.4, where $`\omega <E<\omega _c`$.

## 3 Quenching Weights

If gluons are emitted independently by a hard parton, then the probability $`P(\mathrm{\Delta }E)`$ that this parton loses an additional energy fraction $`\mathrm{\Delta }E`$ can be calcualted from the normalized sum of the emission probabilities for an arbitrary number of $`n`$ gluons which carry away a total energy $`\mathrm{\Delta }E`$:

$`P(\mathrm{\Delta }E)={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}\left[{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle 𝑑\omega _i\frac{dI(\omega _i)}{d\omega }}\right]\delta \left(\mathrm{\Delta }E{\displaystyle \underset{i=1}{\overset{n}{}}}\omega _i\right)\mathrm{exp}\left[{\displaystyle 𝑑\omega \frac{dI}{d\omega }}\right].`$ (12)
The summation over arbitrarily many gluon emissions in (12) can be performed analytically by Laplace transformation . In general, one finds a discrete and a continuous part,

$$P(\mathrm{\Delta }E)=p_0\delta (\mathrm{\Delta }E)+p(\mathrm{\Delta }E).$$
(13)
The discrete weight $`p_0`$ may be viewed as the probability that no additional gluon is emitted due to in-medium scattering and hence no medium-induced energy loss occurs. For finite in-medium pathlength, there is always a finite probability $`p_0`$ that the projectile is not affected by the medium. For infinite in-medium pathlength, however, one finds $`lim_R\mathrm{}p_0=0`$.

In Fig. 2, we show the continuous part of the quenching weight (13) in the multiple soft scattering approximation. The parton energy loss $`\mathrm{\Delta }E`$ is generally seen to increase with increasing characteristic gluon energy $`\omega _c`$ or increasing factor $`R`$. This implies that $`\mathrm{\Delta }E`$ grows with increasing momentum transfer $`\widehat{q}`$ from the medium, and with increasing in-medium pathlength $`L`$, as naively expected.

To understand the negative contributions of the quenching weight for small values $`R`$, one recalls that the gluon energy distribution $`\omega \frac{dI}{d\omega }`$ calculated in (1) is only the medium-induced modification of a radiation pattern $`\frac{dI^{(vac)}}{d\omega }`$ which occurs in the absence of a medium,

$$\omega \frac{dI^{(tot)}}{d\omega }=\omega \frac{dI^{(vac)}}{d\omega }+\omega \frac{dI}{d\omega }.$$
(14)
By writing the corresponding probabilities (11) in Mellin space, one shows that

$$P^{(tot)}(\mathrm{\Delta }E)=_0^{\mathrm{}}𝑑\overline{E}P(\mathrm{\Delta }E\overline{E})P^{(vac)}(\overline{E}).$$
(15)
The probability $`P^{(tot)}(\mathrm{\Delta }E)`$ is normalized to unity and it is positive definite. In contrast, the medium-induced modification of this probability, $`P(\mathrm{\Delta }E)`$, is a generalized probability. It can take negative values for some range in $`\mathrm{\Delta }E`$, as long as its normalization is unity,

$$_0^{\mathrm{}}𝑑\overline{E}P(\overline{E})=p_0+_0^{\mathrm{}}𝑑\overline{E}p(\overline{E})=1.$$
(16)
The discrete weight $`p_0`$ is found to coincide with this normalization condition. A qualitatively similar behaviour of the quenching weights is found in the single hard scattering approximation.

## 4 Angular dependence of radiation probability

The quenching weights discussed above allow to calculate the average energy loss outside an opening angle $`\mathrm{\Theta }`$,

$`\mathrm{\Delta }E(\mathrm{\Theta })`$ $`=`$ $`{\displaystyle 𝑑\omega \omega \frac{dI^{>\mathrm{\Theta }}}{d\omega }(\omega _c,R=\omega _cL)}`$ (17)
$`=`$ $`{\displaystyle 𝑑\overline{E}\overline{E}\left[P(\overline{E},\omega _c,R=\omega _cL)P(\overline{E},\omega _c,R_\chi =\chi ^2\omega _cL)\right]}.`$
This calculation is straightforward since a finite emission angle $`\mathrm{\Theta }`$ can be taken into account simply by reducing the kinematical constraint $`k_{}<\omega `$ in (1) to $`k_{}<\chi \omega `$, $`\chi =\mathrm{sin}\mathrm{\Theta }`$. In Fig.3, we compare the resulting angular radiation patterns obtained in the single hard and multiple soft scattering approximation. In agreement with the discussion following eq. (11) above, we find that for comparable sets of model parameters $`\omega _c`$, $`R`$ and $`\overline{\omega }_c`$, $`\overline{R}`$, $`n_0L`$ respectively, the multiple soft and single hard scattering approximations lead to a comparable angular dependence of $`\mathrm{\Delta }E(\mathrm{\Theta })`$ for $`\mathrm{\Theta }>10^{}`$. For smaller angles, deviations persist which one can understand quantitatively  and which should be regarded as an intrinsic uncertainty of this calculation.

## 5 Application

The quenching weight $`P(\mathrm{\Delta }E)`$ determines the quenching factor

$`Q(p_{})`$ $`=`$ $`{\displaystyle 𝑑\mathrm{\Delta }EP(\mathrm{\Delta }E)\left(\frac{d\sigma ^{\mathrm{vac}}(p_{}+\mathrm{\Delta }E)/dp_{}^2}{d\sigma ^{\mathrm{vac}}(p_{})/dp_{}^2}\right)},`$ (18)
which determines the reduction of transverse momentum leading hadron spectra due to medium-induced energy loss. Alternatively, the quenching weight can be used to determine medium-modified fragmentation functions

$`D_{h/q}^{(\mathrm{med})}(x,Q^2)={\displaystyle _0^1}𝑑ϵP(ϵ){\displaystyle \frac{1}{1ϵ}}D_{h/q}({\displaystyle \frac{x}{1ϵ}},Q^2).`$ (19)
From these, the reduction of leading transverse momentum spectra can be calculated  by convoluting with the perturbative hard matrix elements which are approximately proportional to $`x^6`$,

$$R_{ff}(p_{})=\frac{x_{\mathrm{max}}^6D_{h/q}^{(\mathrm{med})}(x_{\mathrm{max}},p_{}^2)}{x_{\mathrm{max}}^6D_{h/q}(x_{\mathrm{max}},p_{}^2)}|_{p_{}=x_{\mathrm{max}}E_q}.$$
(20)
In Fig.4, we compare the data on the suppression of the $`\pi ^0`$-spectra to suppression factors calculated from (18) and (20). We find that both definitions of the suppression factors lead to quantitatively comparable results. Moreover, the single hard and multiple soft scattering approximation lead to quantitatively comparable results for suitable choices of the model parameters. Remarkably, we find that in the presence of kinematical constraints depleting the infrared region of $`\omega \frac{dI}{d\omega }`$ in Fig.1, the $`p_{}`$-dependence of the quenching factor flattens considerably. This feature is necessary to find in Fig.4 a shallow $`p_{}`$-dependence which is consistent with the data.

The application of the current calculations of parton energy loss to data below $`p_{}<10`$ GeV entails significant theoretical uncertainties. On the one hand, there are plausible competing physics effects which may be relevant at $`p_{}<10`$ GeV (see Ref. for further discussion). On the other hand, the high-energy “eikonal” approximation used in the derivation of (1) becomes questionable if $`\mathrm{\Delta }EE`$ which is the case in Fig.4. This should motivate improved calculations of (1) for which the finite energy effects of the partonic projectile are taken into account.


text:
# 1 Introduction

## 1 Introduction

In a previous note , an exhaustive study of the channel simulations:

$$BV_1V_2,\gamma V,PV,PP,$$

($`V=1^{},P=0^{})`$ has been performed by stressing the helicity formalism and its consequences. General formulas have been established, notably those giving the final angular distributions in the case of the production of two vector mesons decaying into pseudoscalar mesons.

The squared modulus of the decay amplitude has the following form:

$`|A|^2h_{\lambda ,\lambda ^{}}F_{\lambda ,\lambda ^{}}(\theta _1)G_{\lambda ,\lambda ^{}}(\theta _2,\varphi ),`$ (1)
where (summation over $`\lambda (\lambda ^{})`$ is omitted):

* $`h_{\lambda ,\lambda ^{}}`$ is the matrix density element constructed from the weak effective hamiltonian $`H_w^{eff}`$ taken between the initial state ($`B_0`$) and the final state $`f`$.
* $`F_{\lambda ,\lambda ^{}}(\theta _1)`$ and $`G_{\lambda ,\lambda ^{}}(\theta _2,\varphi )`$ are the matrix elements related to the decays $`V_1a_1+b_1`$ and $`V_2a_2+b_2`$ respectively.
* $`\theta _j`$ is the polar angle of particle $`a_j`$ in the rest frame of the resonance $`V_j`$ while $`\varphi `$ is the angular difference $`\varphi _2\varphi _1`$ , where $`\varphi _j`$ is the polar angle of $`a_j`$ in $`V_j`$ rest frame.

$`\lambda (\lambda ^{})`$ being the helicity state of the vector mesons; $`\lambda =1,0,+1`$.

As it can be noticed, the essential parameters for the determination of the decay dynamics are the unknown matrix elements $`h_{\lambda ,\lambda ^{}}`$; while the two other ones, $`F_{\lambda ,\lambda ^{}}(\theta _1)`$ and $`G_{\lambda ,\lambda ^{}}(\theta _2,\varphi )`$, are kinematic (or geometric) parameters because they are completely determined from the Wigner rotation matrices. The reader is referred to the note $`99051`$ for a full kinematic description of the $`B^0`$ decay and the physical significance of the angles $`\theta _{1,2}`$ and $`\varphi `$.

Before dealing with the mathematical determinations of the $`h_{\lambda ,\lambda ^{}}`$ elements, a simple justification of the two vector meson channel is given below.

## 2 Quantum numbers of the $`V_{1}^{}{}_{}{}^{0}V_{2}^{}{}_{}{}^{0}`$ system

In the case of two vector meson $`B^0`$ decay, the most interesting case is the one related to neutral mesons supplemented by the condition $`C|V_{}^{0}{}_{i}{}^{}=|V_{}^{0}{}_{i}{}^{}`$, where $`C`$ is the charge conjugation operator and $`V_{}^{0}{}_{i}{}^{}`$ is a neutral vector meson eigenstate of $`C`$. Some examples of these channels are:

$$\rho ^0\rho ^0,J/\mathrm{\Psi }\rho ^0,J/\mathrm{\Psi }\mathrm{\Phi },\mathrm{\Phi }\mathrm{\Phi }\mathrm{}$$

These vector mesons have, in addition, the parity quantum number equal to $`1`$. Noticing that the total angular momentum of the $`V_1^0V_2^0`$ system: $`\stackrel{}{J}=\stackrel{}{\mathrm{}}+\stackrel{}{S}=\stackrel{}{s_B}`$ is equal to zero and because the total spin $`\stackrel{}{S}=\stackrel{}{s_1}+\stackrel{}{s_2}`$, with $`s_1=s_2=1`$, the orbital angular momentum can have three different values: $`\mathrm{}=S=0,1,2`$.

Thus, parity, charge conjugation and CP quantum numbers of the $`V_1^0V_2^0`$ system can be computed:

$$P(V_1^0V_2^0)=(1)^2(1)^{\mathrm{}},C(V_1^0V_2^0)=(1)^2,$$

$$$$

$$CP(V_1^0V_2^0)=(1)^{\mathrm{}}.$$

We are led to the important result that the CP value of $`V_1^0V_2^0`$ is a mixing of two different eigenvalues $`+1`$ and $`1`$ whatever the initial state ($`B^0`$ or $`\overline{B^0}`$) is. A direct consequence of this result is that CP symmetry is not an exact one.

The above relation does not hold for reactions involving a neutral $`K^{}`$ like:

$$B_d^0K^0\rho ^0,J/\mathrm{\Psi }K^0\mathrm{}$$

because $`K^0`$ and $`\overline{K}^0`$ are two distinct particles; $`C|K^0=|\overline{K^0}|K^0`$.

However, it is worth noticing two interesting features for channels with an intermediate resonance like $`K^0(\overline{K^0})`$:

$$K^0K^+\pi ^{},K^0\pi ^0,$$

$$\overline{K^0}K^{}\pi ^+,\overline{K^0}\pi ^0.$$

The decay channels are in the ratio $`2/3`$ and $`1/3`$ respectively. On one hand, the sign of the charged kaon shows clearly the nature of the neutral $`K^{}`$ from which it comes and consequently the flavour of the original $`B^0(\overline{B^0})`$. So, a neutral $`K^{}`$ decay is a direct way for $`B^0`$ flavour tagging.

On the other hand, when a neutral $`K^0(\overline{K^0})`$ decays into $`K^0(\overline{K^0})\pi ^0`$, the neutral kaon $`K^0(\overline{K^0})`$ is not the true physical particle, because approximately $`50\%`$ of the $`K^0(\overline{K^0})`$ go into $`K_S^0`$ and $`50\%`$ into $`K_L^0`$ respectively and the true detectable particle is $`K_S^0`$ which goes to $`\pi ^+\pi ^{}`$. Thus, in the special channel:

$$B^0(\overline{B^0})K^0(\overline{K^0})\rho ^0,$$

$$K_S^0\pi ^0,$$

tagging the original $`B^0`$ is no longer possible but, the $`K_S^0\pi ^0`$ being a common final state to both $`B^0`$ and $`\overline{B^0}`$, the above relation $`CP=(1)^{\mathrm{}}`$ is still available .

In the following, emphasis will be put on the channels $`K^{0(\pm )}\rho ^0(\omega )`$ and the physical importance of the $`\rho ^0(\omega )`$ mixing for the determination of CP violation.

## 3 $`\rho ^0(\omega )`$ mixing and its consequence

It is well known from hadronic physics that the neutral isovector $`\rho _8`$ and the isosinglet $`\omega _8`$ mix together, leading to the ”true” physical resonances $`\rho ^0`$ and $`\omega `$. On the phenomenological level, this mixing is made possible because of the existence of a common final state to both $`\rho ^0`$ and $`\omega `$ decays :

$$\rho ^0\pi ^+\pi ^{},(BR100\%),$$

$$\omega \pi ^+\pi ^{},(BR2.2\%).$$

In the same framework, it has been established that the $`\pi \pi `$ final state interaction provides a phase shift $`\delta `$ which reaches $`90^{}`$ when the $`\pi \pi `$ invariant mass is at the $`\omega `$ pole ($`M_\omega =782MeV`$) .

This interesting physical property has important consequences in the case where a $`\rho ^0`$ resonance is produced in some $`B^{0\pm }`$ decays like:

$$B^0K^0\rho ^0,(Fig.1)$$

$$B^+K^+\rho ^0,B^{}K^{}\rho ^0,(Fig.2)$$

These decays require both tree (T) and penguin (P) diagrams. As it is emphasized in reference , the amplitude $`A`$ and $`\overline{A}`$ respectively for $`B^+`$ and $`B^{}`$ decays can be set in the following form:

$`A=A^T+A^P=A^T\left(1+r\mathrm{exp}(i\delta )\mathrm{exp}(i\varphi )\right),`$ (2)
$`\overline{A}=\overline{A}^T+\overline{A}^P=A^T\left(1+r\mathrm{exp}(i\delta )\mathrm{exp}(i\varphi )\right),`$ (3)
where:

$`r=\left|{\displaystyle \frac{A^P}{A^T}}\right|,`$ (4)
$`\overline{A}^T=A^T,\overline{A}^P=|A^P|\mathrm{exp}(i\delta )\mathrm{exp}(i\varphi ).`$ (5)
Expressions of $`A`$ and $`\overline{A}`$ displayed above suppose that final state interactions (FSI) arise essentially from the penguin diagrams; this hypothesis is supported by the fact that, to order $`G_F\alpha _s`$ ($`G_F`$ and $`\alpha _s`$ are respectively the Fermi constant and the QCD fine-structure constant), the absorptive part of the transition amplitude is obtained from the penguin diagrams .

In the special case of $`\rho ^0\omega `$ mixing, another hypothesis is made using more intuitive arguments: the phase shift due to the mixing is included in the FSI and it is predominating at the $`\omega `$ pole, justifying the above expressions of $`A`$ and $`\overline{A}`$ that the phase shift $`\delta `$ is principally the one generated by the $`\rho ^0\omega `$ mixing.

By $`CP`$ transformation, the strong phase $`\delta `$ remains unchanged while the weak phase $`\varphi `$, which is related to the CKM matrix elements, changes sign. Thus, the asymmetry parameter $`a_{CP}^{dir}`$ which can reveal direct CP violation can be deduced in the following way:

$`a_{CP}^{dir}={\displaystyle \frac{A^2\overline{A}^2}{A^2+\overline{A}^2}}={\displaystyle \frac{2\mathrm{sin}\delta \mathrm{sin}\varphi }{1+r^2+2r\mathrm{cos}\delta \mathrm{cos}\varphi }}.`$ (6)
It is straightforward to notice that the parameter $`a_{CP}^{dir}`$ depends both on the strong phase and the weak phase and, consequently, the maximum value of $`a_{CP}^{dir}`$ can be reached if $`\mathrm{sin}\delta =1`$, which allows us to state that the strong final state interaction (FSI) among pions coming from the $`\rho ^0\omega `$ decays enhances the direct $`CP`$ violation in the vicinity of the resonance $`\omega `$ mass.

### Simulation of the $`\rho ^0\omega `$ mixing

A simple and phenomenological relation describing the amplitude of the $`\rho ^0\omega `$ mixing is used for the Monte-carlo simulations . In the $`\rho ^0`$ Breit-Wigner, the ($`\rho ^0`$) propagator is replaced by the following one:

$`𝒜={\displaystyle \frac{1}{s_\rho }}+{\displaystyle \frac{T_\omega }{T_\rho }}{\displaystyle \frac{\mathrm{\Pi }_{\rho \omega }}{s_\rho s_\omega }},`$ (7)
where

* $`1/s_V=1/(sM_{V}^{}{}_{}{}^{2}+i\mathrm{\Gamma }_VM_V)`$ is the $`V`$ resonance propagator, $`M_V`$ and $`\mathrm{\Gamma }_V`$ being respectively the mass and the width of the resonance $`V`$.
* $`T_\omega `$ and $`T_\rho `$ are respectively the $`\omega `$ and $`\rho `$ production amplitudes.
* $`\mathrm{\Pi }_{\rho \omega }`$ is the mixing parameter for which recent values come from $`e^+e^{}`$ annihilations:

$`\mathrm{}e(\mathrm{\Pi }_{\rho \omega })=3500\pm 300MeV^2`$ and $`\mathrm{}m(\mathrm{\Pi }_{\rho \omega })=300\pm 300MeV^2.`$

Due to the same physical processes which enter into the production of the $`\rho ^0`$ and $`\omega `$ resonances (they are both made out from $`u\overline{u}`$ and $`d\overline{d}`$ quark pairs with the same weight $`1/2`$), it seems natural to choose $`T_\omega /T_\rho =1`$. So, the squared mass distribution of the $`\pi \pi `$ system becomes simplified and it is given by:

$$d\sigma /dm^2|𝒜(\rho ^0(\omega ))|^2,$$
(8)
where $`𝒜`$ is the amplitude of the two Breit-Wigner given above and $`m`$ is the $`\pi \pi `$ invariant mass.

In Figure $`3`$, are displayed the $`\pi \pi `$ invariant mass spectra for the $`\rho ^0`$ Breit-Wigner and the $`\rho ^0\omega `$ mixing respectively. Because of the very narrow $`\omega `$ width ($`\mathrm{\Gamma }_\omega =8.44MeV`$), we notice a high and narrow peak at the $`\omega `$ pole ($`782MeV`$).

## 4 Dynamics of the $`BV_1V_2`$ decay

The formalism describing the $`B^{0(\pm )}`$ decay into two vector mesons is derived from the general formalism related to the hadronic weak decay of a heavy meson (or heavy quark). It is based on the new concepts introduced by the Heavy Quark Effective Theory (HQET) which involves additionnal symmetry due to the high mass of the heavy quark ($`b`$ or $`c`$ quark) . Technical calculations require a weak effective hamiltonian, $`H_w^{eff}`$, by using the ”Operator Product Expansions” (OPE) pioneered by Wilson and which involve field operators describing both tree and penguin diagrams, the last ones include both QCD and electroweak penguins (Figures 1 and 2).

The general form of $`H_w^{eff}`$ is given by:

$$H_{w}^{}{}_{}{}^{eff}=\frac{G_F}{\sqrt{2}}\underset{q=d,s}{}\left(V_{ub}V_{uq}^{}(c_1O_1+c_2O_2)V_{tb}V_{tq}^{}\underset{i=3}{\overset{10}{}}c_iO_i\right),$$
(9)
where $`c_i`$ are the Wilson coefficients and $`O_i`$ are field operators with dimension $`d4`$; they are computed at an energy scale $`\mu `$ which is identified, here, with the $`b`$ quark mass $`m_b`$.

In the case of charmless $`B`$ decays, Wilson coefficients have been calculated by Buchalla et al . These coefficients represent the perturbative part of the weak hamiltonian, they are estimated by the Renormalization Group techniques and their values depend on the renormalization scheme which is used. Their physical significance is the weight of each field operator $`O_i(\mu )`$ entering in the weak hamiltonian $`H_w^{eff}`$. From reference , the values of $`c_i`$ which have been computed at the energy scale $`\mu =m_b`$ are:

$`\begin{array}{cc}c_1=0.3125,c_2=1.1502,\hfill & \\ c_3=0.0174,c_4=0.0373,\hfill & \\ c_5=0.0104,c_6=0.0459,\hfill & \\ c_7=1.050\times 10^5,c_8=3.839\times 10^4,\hfill & \end{array}`$
$`c_9=0.0101,c_{10}=1.959\times 10^3.`$ (11)
The first two coefficients, $`c_1`$ and $`c_2`$, are related to the tree diagrams and they show clearly their dominance with respect to the penguin ones. Coefficients $`c_3c_6`$ correspond to QCD penguin operators while $`c_7c_{10}`$ are related to the EW ones.

However, those values of $`c_i`$ must be modified when renormalization of operator $`O_i`$ at one-loop order is taken into account.

Detailed expressions of operators $`O_i(\mu )`$ and their physical interpretation are given in reference .

Thus, a general form for the weak decay amplitude into a final state $`f`$ can be expressed like:

$$A(B^0f)=f|H_w^{eff}|B^0=\frac{G_F}{\sqrt{2}}\underset{i=1}{\overset{10}{}}\underset{q=d,s}{}\lambda _{q}^{}{}_{}{}^{i}c_i(\mu )f|O_i(\mu )|B^0,$$
(12)
where $`\lambda _{q}^{}{}_{}{}^{i}`$ is the product of two CKM matrix elements: $`V_{ub}V_{uq}^{}`$ (for $`i=1,2`$) or $`(V_{tb}V_{tq}^{})`$ (for $`i=3,\mathrm{},10`$).

The hadronic matrix elements $`f|O_i(\mu )|B^0`$ represent the non-perturbative contribution to the amplitude $`A(B^0f)`$. Usually, they are estimated according to some specific models: Non Relativistic Quark Model (NRQM), Form Factor models (BSW) and especially the Lattice QCD calculations.

In the following, calculation of the hadronic matrix elements is performed in the framework of the BSW model  from which form factors are derived by the knowledge of the hadronic wave functions for both initial and final states.

## 5 Determination of the density-matrix elements

The $`B^0`$ decay into two vector mesons requires the helicity formalism which has been intensively used in the previous paper . To each vector meson (spin $`1`$) is assigned a set of three polarization 4-vectors defined in this way:

$$ϵ_1=(0,\stackrel{}{ϵ_1}),ϵ_2=(0,\stackrel{}{ϵ_2}),ϵ_3=(|\stackrel{}{k}|/m,E\widehat{k}/m),$$
(13)
and verifying the following relations:

$$ϵ_{i}^{}{}_{}{}^{2}=1,ϵ_iϵ_j=0,\mathrm{with}ij,$$
(14)
where $`m,E,\stackrel{}{k}`$ are respectively the mass, the energy and the momentum of the vector meson; $`\widehat{k}`$ is defined as the unit vector along the vector momentum, $`\widehat{k}=\stackrel{}{k}/|\stackrel{}{k}|`$.

The three vectors $`\stackrel{}{ϵ_1},\stackrel{}{ϵ_2}`$ and $`\stackrel{}{ϵ_3}=E\widehat{k}/m`$ form an orthogonal basis; $`ϵ_1`$ and $`ϵ_2`$ are called the transverse polarization vectors while $`\stackrel{}{ϵ_3}`$ is the longitudinal polarization one.

From that basis, an helicity basis is defined according to:

$$ϵ(+)=\frac{\left(ϵ_1+iϵ_2\right)}{\sqrt{2}},ϵ()=\frac{\left(ϵ_1iϵ_2\right)}{\sqrt{2}},ϵ(0)=ϵ_3.$$
(15)
These 4-vectors are eigenvectors of the helicity operator $``$ with the eigenvalues $`\lambda =+1,1`$ and $`0`$ respectively. For a clear account of the helicity basis for a spin $`1`$ particle, the reader can consult the book of Dewitt-Smith .

In the case of two vector mesons coming from the $`B`$ decay, their 4-momenta are defined in the $`B`$ rest frame and their corresponding polarization vectors are correlated because $`\widehat{k_1}=\widehat{k_2}`$. For an explicit calculation of their spatial components, see the appendix A.

The weak hadronic amplitude is then decomposed on the helicity basis according to the general formalism developed by the authors BSW . This method allows one to obtain two interesting results:

$``$ the contribution of the tree and penguin operators to the global amplitude via the helicity states.

$``$ the total contribution of each helicity state.

A way of illustrating this method is to study the channel: $`B^0(\overline{B}^0)K^0(\overline{K}^0)\rho ^0`$.

$`(i)`$ First of all, the mass of each resonance ($`K^0`$ and $`\rho ^0`$ ) is generated according to a relativistic Breit-Wigner:

$$\frac{d\sigma }{dM^2}=C\frac{\mathrm{\Gamma }_RM_R}{(M^2M_R^2)^2+(\mathrm{\Gamma }_RM_R)^2},$$
(16)
$`C`$ being a normalization constant.

$`(ii)`$ The weak hadronic matrix element is expressed as the sum of three helicity matrix elements; each one of the form, $`H_\lambda =V_1V_2|H_{w}^{}{}_{}{}^{eff}|B`$, is defined by gathering all the Wilson coefficients of both tree and penguin operators. Linear combinations of those coefficients arise like: $`c_{t1}^\rho ,c_{p1}^\rho ,`$ and $`c_{p2}^\rho `$ (see Appendix B) and the helicity amplitude $`H_\lambda `$ gets the following expression:

$$H_\lambda =(V_{ub}V_{us}^{}c_{t_1}^\rho V_{tb}V_{ts}^{}c_{p_2}^\rho )\{\beta _1\epsilon _{\alpha \beta \gamma \delta }ϵ_K^\alpha (\lambda )ϵ_\rho ^\beta (\lambda )P_B^\gamma P_K^\delta $$

$$+i(\beta _2ϵ_K^{}(\lambda )ϵ_\rho ^{}(\lambda )\beta _3(ϵ_K^{}(\lambda ).P_B)(ϵ_\rho ^{}(\lambda ).P_B))\}$$

$$+(V_{tb}V_{ts}^{}c_{p_1}^\rho )\{\beta _4\epsilon _{\alpha \beta \gamma \delta }ϵ_\rho ^\alpha (\lambda )ϵ_K^\beta (\lambda )P_B^\gamma P_\rho ^\delta ,$$

$$+i(\beta _5ϵ_\rho ^{}(\lambda )ϵ_K^{}(\lambda )\beta _6(ϵ_\rho ^{}(\lambda ).P_B)(ϵ_K^{}(\lambda ).P_B))\}$$
(17)
with:

$``$ $`\epsilon _{\alpha \beta \gamma \delta }`$: antisymmetric tensor in the Minkowski space.

$``$ $`\beta _{1,4}=\frac{G_F}{2}f_{\rho ,K}m_{\rho ,K^{}}\frac{2}{m_B+m_{K^{},\rho }}V^{BK^{},\rho }(m_{\rho ,K^{}}^2)`$.

$``$ $`\beta _{2,5}=\frac{G_F}{2}f_{\rho ,K}m_{\rho ,K^{}}(m_B+m_{K^{},\rho })A_1^{BK^{},\rho }(m_{\rho ,K^{}}^2)`$.

$``$ $`\beta _{3,6}=\frac{G_F}{2}f_{\rho ,K}m_{\rho ,K^{}}\frac{2}{m_B+m_{K^{},\rho }}A_2^{BK^{},\rho }(m_{\rho ,K^{}}^2)`$.

$``$ $`f_K`$, $`f_\rho `$: respectively $`K^0`$ and $`\rho ^0`$ decay constants.

$``$ $`V^{BK^{},\rho }`$, $`A_i^{BK^{},\rho }`$: respectively Vector and Axial form factors (see Appendix C).

$``$ $`ϵ_{K,\rho }(\lambda )`$: $`K^0`$, $`\rho ^0`$ polarization vectors expressed in the $`B`$ rest frame.

It is worth noticing that the tensorial terms which enter $`H_\lambda `$ become simplified in the $`B`$ rest frame because the $`B`$ 4-momentum is given by $`P_B=(m_b,\stackrel{}{0})`$. Then, using the orthogonality properties of $`ϵ_j(\lambda )`$, the helicity amplitude $`H_\lambda `$ acquires a much simpler expression than above:

$$H(\lambda )=iB(\lambda )(V_{ub}V_{us}^{}c_{t_1}^\rho V_{tb}V_{ts}^{}c_{p_2}^\rho )+iC(\lambda )(V_{tb}V_{ts}^{}c_{p_1}^\rho ),$$
(18)
with:

$$B(0)=\beta _2\frac{m_B^2(m_K^2+m_\rho ^2)}{2m_Km_\rho }\beta _3\frac{|\stackrel{}{p}|^2m_B^2}{m_Km_\rho },$$

$$C(0)=\beta _5\frac{m_B^2(m_K^2+m_\rho ^2)}{2m_Km_\rho }\beta _6\frac{|\stackrel{}{p}|^2m_B^2}{m_Km_\rho },$$

$$B(\pm 1)=\beta _1m_B|\stackrel{}{p}|\beta _2,$$

$$C(\pm 1)=\beta _4m_B|\stackrel{}{p}|\beta _5,$$
(19)
$`|\stackrel{}{p}|`$ being the common momentum to $`V_1`$ and $`V_2`$ particles in the $`B`$ rest frame.

$`(iii)`$ Expressing the CKM matrix elements according to Wolfenstein parametrization :

$$V_{CKM}=\left[\begin{array}{ccc}1\frac{\lambda ^2}{2}& \lambda & A\lambda ^3(\rho i\eta )\\ \lambda & 1\frac{\lambda ^2}{2}& A\lambda ^2\\ A\lambda ^3(1\rho i\eta )& A\lambda ^2& 1\end{array}\right]+O(\lambda ^4),$$
(20)
where we use:

$`A=0.815`$, $`\lambda =0.2205`$: well known

$`0.09<\rho <0.254`$, $`0.323<\eta <0.442`$.

Taking into account the preceding relations, we arrive at the final form for the amplitudes $`H_\lambda `$:

$$H\left(\genfrac{}{}{0pt}{}{0}{\pm 1}\right)=A\lambda ^2\{[(\eta \lambda ^2c_{t_1}^\rho \mathrm{}m(c_{p_2}^\rho ))B\left(\genfrac{}{}{0pt}{}{0}{\pm 1}\right)\mathrm{}m(c_{p_1}^\rho )C\left(\genfrac{}{}{0pt}{}{0}{\pm 1}\right)]$$

$$+i[(\rho \lambda ^2c_{t_1}^\rho +\mathrm{}e(c_{p_2}^\rho ))B\left(\genfrac{}{}{0pt}{}{0}{\pm 1}\right)+\mathrm{}e(c_{p_1}^\rho )C\left(\genfrac{}{}{0pt}{}{0}{\pm 1}\right)]\},$$
(21)
from which the density-matrix elements $`h_{\lambda ,\lambda ^{}}`$ can be derived automatically;

$`h_{\lambda ,\lambda ^{}}=H_\lambda H_{\lambda ^{}}^{}{}_{}{}^{}`$.

Due to the hermiticity of the matrix ($`h_{\lambda ,\lambda ^{}}`$), only six elements must be calculated and, furthermore, a normalization condition is applied:

$$N(h_{++}+h_{00}+h_{})=1,$$
(22)
($`N`$ being the normalization constant) which makes easier the comparison of the modulus of the different matrix elements.

In the next histograms (Fig.4 - Fig.9) are displayed the spectra of $`h_{\lambda ,\lambda ^{}}`$ for different values of the Wolfenstein parameters $`\rho `$ and $`\eta `$. In our study, these spectra are obtained for the four couples of values: $`(0.09,0.323)`$; $`(0.09,0.442)`$; $`(0.254,0.323)`$ and $`(0.254,0.442)`$. But, due to the fact that some density matrix elements do not vary too much with $`\rho `$ and $`\eta `$, in most cases only the spectra corresponding to the first couple of values are shown. All the histograms correspond to a sample of $`20000`$ generated events.

It is important to notice that large spectrum of values for $`h_{\lambda ,\lambda ^{}}`$ are obtained and not single ones because of the broad range of both the $`\rho ^0`$ resonance mass and the common momentum $`|\stackrel{}{p}|`$ (see the analytical expressions of $`B(\lambda )`$ and $`C(\lambda )`$ given above).

* Whatever the values of $`\rho `$ and $`\eta `$ are, the dominant value of $`h_{++}=|H_{+1}|^2`$ is $`10^2`$, numerical result which is proved too by complete analytical calculations. Thus the dominant polarization state is the longitudinal one because $`h_{00}=|H_0|^260\%`$, its mean value being around $`85\%`$ (Fig.4).
* Due to the tiny value of $`|H_{+1}|`$, the modulus of the non-diagonal elements $`h_+=H_+H_{}^{}{}_{}{}^{}`$ and $`h_{+0}=H_+H_{0}^{}{}_{}{}^{}`$ are usually smaller than $`0.2`$; while the modulus of $`h_0=H_{}H_{0}^{}{}_{}{}^{}`$ can reach $`0.5`$ (Fig.5).

Fig.6 and Fig.7 display the variations of the diagonal matrix elements $`h_{},h_{00}`$ and $`h_{++}`$ with respect to the four sets of $`\rho `$ and $`\eta `$ values: it can be seen that $`h_{++}`$ has always a tiny value and $`h_{00}`$ is always dominant. Other physical features appear: $`h_{00}`$ is very sensitive to the parameter $`\eta `$; its spectrum is rather wide for $`\eta =0.323`$, while it is bounded between $`0.8`$ and $`1.0`$ for $`\eta =0.442`$. For a fixed value of $`\eta `$, no noticeable variation with the parameter $`\rho `$ is seen. Fig.8 shows the real and imaginary parts of the non-diagonal elements $`h_0,h_+,h_{+0}`$ respectively for $`\rho =0.09`$ and $`\eta =0.323`$. It is worth noticing that both real and imaginary parts of $`h_+`$ and $`h_{+0}`$ are too small and close to zero. Due to the importance of $`h_+`$ matrix element in the $`\varphi `$ angle distribution (see Section 6), a full study of both real and imaginary parts of $`h_+`$ has been done. Fig.9 shows the corresponding spectra according to the values of $`\rho `$ and $`\eta `$. It can be deduced that the real and imaginary parts have very similar distributions and both are dominated by small values ($`0.05`$).

## 6 Decays of vector mesons $`V_1V_2`$ into two pseudoscalar mesons

The matrix elements derived above allow us to compute the degrees of polarization of each resonance like $`K^{}`$ or $`\rho ^0`$. The angular distributions of the pseudoscalar mesons in each $`V_i`$ rest frame depend on:
$`(i)`$ the spin $`1`$ of the vector meson $`V_i`$.
$`(ii)`$ the weight of each helicity state.
$`(iii)`$ the correlations among the helicity states of the two vector mesons.

Complete analytical expression of the final angular distributions is the following one:

$`{\displaystyle \frac{d^3\mathrm{\Gamma }}{d\mathrm{cos}\theta _1d\mathrm{cos}\theta _2d\varphi }}`$ $``$ $`(h_{++}+h_{})\mathrm{sin}^2\theta _1\mathrm{sin}^2\theta _2/4+h_{00}\mathrm{cos}^2\theta _1\mathrm{cos}^2\theta _2`$ (23)
$`+\left(\mathrm{}e(h_{+0})\mathrm{cos}\varphi \mathrm{}m(h_{+0})\mathrm{sin}\varphi +\mathrm{}e(h_0)\mathrm{cos}\varphi \mathrm{}m(h_0)\mathrm{sin}\varphi \right)\mathrm{sin}2\theta _1\mathrm{sin}2\theta _2/4`$
$`+\left(\mathrm{}e(h_+)\mathrm{cos}2\varphi \mathrm{}m(h_+)\mathrm{sin}2\varphi \right)\mathrm{sin}^2\theta _1\mathrm{sin}^2\theta _2/2.`$
Angles $`\theta _1`$ , $`\theta _2`$ and $`\varphi `$ have been defined in Section $`1`$.

Explicit angular distributions for polar and azimuthal angles can be derived from the relation above. It is interesting to notice that, due to the pseudoscalar nature of the final particles, angles $`\theta _1`$ and $`\theta _2`$ have the same distributions:

* $`d\sigma /d\mathrm{cos}\theta _{1,2}(3h_{00}1)\mathrm{cos}^2\theta _{1,2}+(1h_{00}).`$
* $`d\sigma /d\varphi (1+2(\mathrm{}e(h_+)\mathrm{cos}2\varphi \mathrm{}m(h_+)\mathrm{sin}2\varphi )).`$

In Fig.10 are displayed respectively the $`\mathrm{cos}\theta `$ distribution and the azimuthal angle $`\varphi `$ one. Some comments on these curves are necessary:

* The $`\mathrm{cos}\theta `$ distribution is practically the same whatever the values of $`\rho `$ and $`\eta `$ are; no sensitivity to particular values of $`\rho `$ and $`\eta `$ is seen.
* As far as angle $`\varphi `$ is concerned, its distribution depends on the matrix element $`h_+`$. Despite the fact that $`\mathrm{}e(h_+)`$ and $`\mathrm{}m(h_+)`$ do not exhibit sensitive differences (see Fig.9), those parameters present some dependence upon $`\rho `$ and $`\eta `$: full curve corresponds to $`\rho =0.09,\eta =0.323`$; while dashed one is related to $`\rho =0.254,\eta =0.442`$. A visible discrepancy among these two curves is seen.

## 7 Perspectives and conclusion

$``$ Thanks to the HQET approach and the OPE formalism which is used, we have at our disposal rigorous and complete calculations of the dynamics of the $`B^{0\pm }`$ decays into two vector mesons. This formalism is available for all charmless $`B`$ decays provided the spin of the intermediate resonance(s) is less or equal $`1`$; the only changes which must be taken into account are the $`V_{CKM}`$ matrix elements, the masses and the widths of the resonances involved in each decay.

$``$ In the case of leptonic decay of one resonance, like $`J/\mathrm{\Psi }e^+e^{},\mu ^+\mu ^{}`$, the angular distributions are modified because of the spin $`1/2`$ final leptons; which require the use of other Wigner rotation matrices. Those calculations have been already done in our first paper .

$``$ In the case where a $`(c\overline{c})`$ bound state or a charmed meson is produced like:

$$B^0J/\mathrm{\Psi }\rho ^0,D^{}X(X=\rho ^0,\omega ,K^0),$$

the Wilson coefficients involved in the effective hamiltonian have to be modified, but we do not expect big change with respect to the $`c_i(c_i^{})`$ coefficients used in the present paper.

$``$ Other interesting consequences arise from this formalism: it can be easily extended to the numerous channels like: $`BVP,PP`$ where one or two pseudoscalar mesons ($`P=0^+`$) are produced directly from the $`B`$ decay. Because of the simple equality $`\lambda (P)=\lambda (V)=0`$, the number of helicity states is reduced from $`3`$ to $`1`$.

$``$ An important point which has been mentionned in the present note is the role of the $`\rho ^0\omega `$ mixing and its consequence for the determination of the direct CPV parameter (Section $`3`$ and reference ). Tagging of $`B^+`$ and $`B^{}`$ is made easy thanks to the $`K^+`$ and $`K^{}`$ mesons coming from the cascade decays. In our opinion, we can also exploit all the angular distributions of the final particles (and their correlations) in order to detect an eventual discrepancy which can arise between the $`B^+`$ and $`B^{}`$ decays respectively. However, a complete study of those channels and their simulations require the knowledge of the strong phase shift $`\delta `$ (mentionned in Section 3) according to the $`\pi \pi `$ invariant mass. Work is in progress.

$``$ Those calculations and simulations can be implemented into SICBMC, the Monte-Carlo generator of the LHCb experiment, in order to perform afterwards a full analysis of the simulated channels.

### Acknowledgements

The authors are very grateful to Dr P.Perret, leader of the LHCb Clermont-Ferrand team, for his advices and his suggestions.

One of us (Z.J.A.) is very indebted to Professor A.W.Thomas, Director of the Special Research Centre for the Subatomic Structure of Matter, for the very exciting and illuminating discussions he got with him about the QCD penguin diagrams and their importance in the evaluation of the $`B^0`$ decay width.

This work was supported in part by the Australian Research Council and the University of Adelaide.

## Appendix

## Appendix A Polarizations in $`B_d^0`$ rest frame

Momentum:

$$\stackrel{}{k}_K=\stackrel{}{k}_\rho =\stackrel{}{k}=\left(\begin{array}{c}k\mathrm{sin}\theta \mathrm{cos}\varphi \\ k\mathrm{sin}\theta \mathrm{sin}\varphi \\ k\mathrm{cos}\theta \end{array}\right),$$

where $`\theta `$ and $`\varphi `$ are respectively polar and azimuthal angles of the produced $`K^0`$.

Longitudinal polarization:

$$ϵ_K(0)=(\frac{|\stackrel{}{k}|}{m_K},\frac{E_K}{m_K}\widehat{k}),ϵ_\rho (0)=(\frac{|\stackrel{}{k}|}{m_\rho },\frac{E_\rho }{m_\rho }(\widehat{k})).$$

Tranversal polarizations :

$$\stackrel{}{ϵ}_K(1)=\left(\begin{array}{c}\mathrm{cos}\theta \mathrm{cos}\varphi \\ \mathrm{cos}\theta \mathrm{sin}\varphi \\ \mathrm{sin}\theta \end{array}\right)=\stackrel{}{ϵ}_\rho (1),$$

$$\stackrel{}{ϵ}_K(2)=\left(\begin{array}{c}\mathrm{sin}\varphi \\ \mathrm{cos}\varphi \\ 0\end{array}\right)=\stackrel{}{ϵ}_\rho (2).$$

Helicity frame :

$$ϵ_K(+)=\left(ϵ(1)+iϵ(2)\right)/\sqrt{2},ϵ_K()=\left(ϵ(1)iϵ(2)\right)/\sqrt{2},$$

$$\stackrel{}{ϵ}_K(+)=\left(\begin{array}{c}\mathrm{cos}\theta \mathrm{cos}\varphi i\mathrm{sin}\varphi \\ \mathrm{cos}\theta \mathrm{sin}\varphi +i\mathrm{cos}\varphi \\ \mathrm{sin}\theta \end{array}\right)/\sqrt{2}=\stackrel{}{ϵ}_K^{}()=\stackrel{}{ϵ}_\rho (),$$

$$\stackrel{}{ϵ}_K()=\left(\begin{array}{c}\mathrm{cos}\theta \mathrm{cos}\varphi +i\mathrm{sin}\varphi \\ \mathrm{cos}\theta \mathrm{sin}\varphi i\mathrm{cos}\varphi \\ \mathrm{sin}\theta \end{array}\right)/\sqrt{2}=\stackrel{}{ϵ}_K^{}(+)=\stackrel{}{ϵ}_\rho (+).$$

## Appendix B Wilson’s coefficients

We use, in the case of the $`\rho ^0`$ production, the following linear combinations of the effective Wilson coefficients:

$`c_{t1}^\rho =c_1^{}+{\displaystyle \frac{c_2^{}}{N_c}},`$
$`c_{p1}^\rho =(c_4^{}+{\displaystyle \frac{c_3^{}}{N_c}})+{\displaystyle \frac{1}{2}}(c_{10}^{}+{\displaystyle \frac{c_9^{}}{N_c}}),`$
$`c_{p2}^\rho ={\displaystyle \frac{3}{2}}(c_7^{}+{\displaystyle \frac{c_8^{}}{N_c}}+c_9^{}+{\displaystyle \frac{c_{10}^{}}{N_c}}),`$
where $`c_{t1}^\rho `$ relative to tree diagram, $`c_{pi}^\rho `$ relative to penguin diagram and $`0.98<N_c<2.01`$.

When $`q^2/m_b^2=0.3`$:
$`c_1^{}=0.3125,c_2^{}=1.1502,`$

$`c_3^{}=2.443\times 10^2+1.543\times 10^3i,c_4^{}=5.808\times 10^24.628\times 10^3i,`$

$`c_5^{}=1.733\times 10^2+1.543\times 10^3i,c_6^{}=6.668\times 10^24.628\times 10^3i,`$

$`c_7^{}=1.435\times 10^42.963\times 10^5i,c_8^{}=3.839\times 10^4,`$

$`c_9^{}=1.023\times 10^22.963\times 10^5i,c_{10}^{}=1.959\times 10^3.`$

When $`q^2/m_b^2=0.5`$:
$`c_1^{}=0.3125,c_2^{}=1.1502,`$

$`c_3^{}=2.120\times 10^2+2.174\times 10^3i,c_4^{}=4.869\times 10^21.552\times 10^2i,`$

$`c_5^{}=1.420\times 10^2+5.174\times 10^3i,c_6^{}=5.729\times 10^21.552\times 10^2i,`$

$`c_7^{}=8.340\times 10^59.938\times 10^5i,c_8^{}=3.839\times 10^4,`$

$`c_9^{}=1.017\times 10^29.938\times 10^5i,c_{10}^{}=1.959\times 10^3.`$

## Appendix C Form factors (BSW model)

| \[.3cm\] | $`V`$ | $`A_1`$ | $`A_2`$ |
| --- | --- | --- | --- |
| \[.3cm\] $`BK^{}`$ | $`\frac{0.369}{1m_\rho ^2(GeV^2)/5.43^2(GeV^2)}`$ | $`\frac{0.328}{1m_\rho ^2(GeV^2)/5.43^2(GeV^2)}`$ | $`\frac{0.331}{1m_\rho ^2(GeV^2)/5.43^2(GeV^2)}`$ |
| \[.3cm\] $`B\rho `$ | $`\frac{0.329}{1m_K^{}^2(GeV^2)/5.32^2(GeV^2)}`$ | $`\frac{0.283}{1m_K^{}^2(GeV^2)/5.32^2(GeV^2)}`$ | $`\frac{0.283}{1m_K^{}^2(GeV^2)/5.32^2(GeV^2)}`$ |
| \[.3cm\] |  |  |  |

For further details, see reference  and literature quoted therein.


text:
# Development and operation of the twin radio frequency single electron transistor for cross-correlated charge detection

## I Introduction

The radio frequency single electron transistor (rf-SET) is a nano-scale device capable of detecting fractions of an electron charge on sub-microsecond time-scales Schoelkopf et al. (1998). In addition to being an exquisite electrometer Devoret and Schoelkopf (2000), the rf-SET may also find application as a single photon Komiyama et al. (2000); Schoelkopf et al. (1999) and particle detector Bouchiat et al. (1997), micro-mechanical displacement sensor Blencowe and Wybourne (2000), ultra-high impedance voltage amplifier Segall et al. (2002); Stevenson et al. (2001, 2002) and ultimately as a fast single-shot read out device in solid-state quantum computers Aassime et al. (2001a). At present the sensitivity of rf-SETs is primarily limited by the post-amplifier, which is typically a GaAs high-electron mobility transistor (HEMT) operating at $`T`$ 4K. In the near future, this limitation is likely to be overcome via the use of rf amplifiers based on a micro-strip dc superconducting quantum interference device (SQUID) which have demonstrated noise temperatures near 100mK at frequencies suitable for rf-SET operation Bradley et al. (2003). At low frequencies (below 100kHz), rf-SETs suffer from 1/f-type charge noise associated with the motion of charge in the SET dielectric tunnel barriers and substrate. In particular telegraph noise, that results from charge fluctuators that are strongly coupled to the SET, poses a challenge to the widespread use of rf-SETs as both detectors and readout devices.

In an effort to address the issue of telegraph noise we have utilized two superconducting rf-SETs to form a cross-correlation charge detector, termed a twin rf-SET. The cross-correlation technique is advantageous in suppressing charge noise originating from fluctuating traps in the surrounding material system and constitutes a means of discriminating read out signals from spurious charge fluctuations on micro-second time-scales Buehler et al. (2003). Independent control of the source-drain bias across each superconducting SET permits both devices to simultaneously operate at the point of maximum sensitivity and measurement efficiency Clerk et al. (2002). By making use of wavelength division multiplexing (WDM) the twin rf-SET provides two (or more) outputs for cross-correlation measurements using only a single cryogenic following amplifier and a single transmission line. Such multiplexed arrangements are highly desirable in instances where arrays of rf-SETs are used as detectors Stevenson et al. (2002).

The article is organized as follows: in Sec. II we briefly review the rf-SET. Sec. III introduces the twin rf-SET and the multiplexing technique. We also discuss our demodulation method that makes use of a mixer-based detection scheme. In Sec. IV we discuss the operation of the resonant circuits and focus on the degree of cross-talk between the circuits and the issues that affect independent operation of the SETs. Additionally we present electromagnetic (EM) simulation results of our circuit and show good agreement with our data. Sec. V presents frequency and time domain data for the case of simultaneous operation. We demonstrate charge sensitivities of $`\delta q_{SET1}=7.5\mu e/\sqrt{Hz}`$ and $`\delta q_{SET2}=4.4\mu e/\sqrt{Hz}`$ and show how the technique of cross-correlation can be used to suppress spurious charge noise on microsecond time-scales.

## II rf-SET

Normal state SET charge detectors are based on the Coulomb blockade Grabert and Devoret (1992) of electrons tunneling across two low capacitance junctions $`J_{1,2}`$ ($`C_1,R_1`$ and $`C_2,R_2`$). An ‘island’ with a characteristic charging energy $`E_C=e^2/2C_\mathrm{\Sigma }`$, where $`C_\mathrm{\Sigma }`$ is the capacitive coupling to the environment, is formed between the tunnel junctions. The sequential electron tunneling current across the island can be controlled by the voltage on a gate which is capacitively coupled to the island. Coulomb blockade is effective for temperatures $`T<E_C/k_B`$ and for tunneling resistances $`R_{1,2}`$ of order of the resistance quantum $`h/e^2`$ 26k$`\mathrm{\Omega }`$. These restrictions generally limit the operation of SETs to cryogenic temperatures. Traditionally, the bandwidth of SETs has been limited to a few kHz by the large $`RC`$ time constant associated with the SET output resistance and the capacitance of the wiring from the SET to the room temperature amplifiers. In contrast to the conventional SET, the rf-SET Schoelkopf et al. (1998) makes use of a $`LC`$ impedance transformer to match the impedance of the SET to the characteristic impedance of a transmission line. In this regime, either the reflected Schoelkopf et al. (1998) or transmitted Fujisawa and Hirayama (2000) power of an incident rf carrier wave is a function of the SET resistance. Mapping the device resistance to changes in the reflected or transmitted power thereby allows fast charge detection. In the case of reflection measurements, the incident rf carrier is coupled to the $`LCR`$ matching network via a directional coupler and the reflected power is coupled to a low noise cryogenic HEMT amplifier. The $`LCR`$ circuit is described by the loaded quality factor $`Q=(Z_{LC}/R_{SET}+Z_0/Z_{LC})^1`$, where $`Z_0`$ is the characteristic impedance of the transmission line and $`R_{SET}`$ is the differential resistance of the SET and $`Z_{LC}=\sqrt{L/C}`$. The (loaded) $`Q`$ factor is inversely proportional to the bandwidth of the impedance matching network.

SETs fabricated from metals that exhibit superconductivity at these temperatures generally also display rich current - voltage characteristics in association with the interplay between charging effects and Cooper pair transport across junctions $`J_1`$ and $`J_2`$. Notably, in the vicinity of the threshold voltage for current transport (defined by the superconducting gap and the charging energy), resonant tunneling processes can be observed which are governed by the dynamics of coherent Cooper pair and incoherent quasi-particle transport Fulton et al. (1989). Superconducting SETs biased to a Josephson quasi-particle resonance (JQP) have recently been explored in conjunction with reading-out a quantum computer and shown to approach quantum limits for measurement efficiency Clerk et al. (2002).

## III Twin rf-SET setup

The technique of wavelength division multiplexing has long been applied in communications engineering as a means of utilizing channel resources efficiently. In line with earlier work on SQUIDs Muck (1991), recent work by Stevenson et al., Stevenson et al. (2002) has explored the possibility of using multiplexed rf-SETs to amplify signals from high impedance photon detectors for application in astronomy. Here we extend these results and note that for our case where the SETs are independently biased using bias-tees, a lumped-element circuit analysis is generally inadequate in describing our data. Consequently we present circuit simulation results that account for the length of transmission lines used in our setup Fulton and Dolan (1987).

Turning now to the specifics of our setup, we achieve multiplexing of signals from the two rf-SETs via two tuned $`LC`$ impedance transformers (Figure 1). Each tank circuit can be approximated by a parasitic capacitance $`C`$ to ground and a chip inductor with inductance $`L_{1,2}`$. The circuits are tuned by choosing appropriate inductances to transform the impedance of the SET (typically 40 - 100k$`\mathrm{\Omega }`$) downward towards a characteristic impedance of 50$`\mathrm{\Omega }`$ at the respective resonant frequency $`\omega _{1,2}=1/\sqrt{CL_{1,2}}`$. Figure 1 shows a schematic of our cryogenic setup, including two separate bias-tees for independent dc-biasing of each SET. In this setup the length of coaxial transmission line between the ‘T’ section and the SETs (see Figure 1, $`d_1`$ dashed line and $`d_2`$ solid line) is a critical parameter determining the cross-talk between the SET tank circuits. The ‘T’ arrangement can be considered a single shunt ‘stub’, where the length of the stub and transmission line between the load and the stub position tune the admittance seen looking into the line Pozar (1998). In the present experimental arrangement we place our cryogenic HEMT amplifier at the $`T`$ = 4K stage of our dilution refrigerator to avoid the use of cryogenic feed-throughs and to reduce losses between the SET and the amplifier. The twin rf-SET data presented here was taken using tank circuits with resonance frequencies $`f_1335`$MHz and $`f_2360`$MHz, using inductances of 780nH and 660nH, with a total parasitic capacitance estimated to be $`0.3`$pF for each SET. The choice of frequencies corresponds to the bandwidth where our cryogenic rf-amplifier Fulton and Dolan (1987) maintains optimum noise performance and in addition, where chip inductors are readily available. Variations in the parasitic capacitance $`C`$ (and resonance frequency $`\omega `$) for different samples can be attributed to changes in the length of bond-wires from the inductors to the sample contact pads, changing $`C`$ by tens of $`fF`$.

The incident carrier wave is produced by combining the output of two independent rf signal generators (Figure 1). Each generator output is fed to a directional coupler in order to tap off power (-16dB) for the SETs. The remaining signal (-0.11dB) is directed to the mixers as a reference signal. A subsequent splitter/combiner is used to feed the two frequencies to one semi-rigid coaxial waveguide into the cryostat. After the signal enters the cryostat the rf power is coupled to the tank circuits using another directional coupler. The signals $`f_{1,2}`$ are directed to SET1,2 by the tuned impedance transformers, tank1,2. The reflected signal is then coupled to the cryogenic amplifier (40dB gain) and a room temperature amplifier (40dB gain). A power splitter feeds the reflected rf signal to two mixers for demodulation. High roll-off, narrow bandpass filters Fulton and Dolan (1987) are used before both the local oscillator (LO) and RF mixer inputs to suppress inter-modulation distortion Pozar (1998) and higher order components. Active phase-shifters are used to ensure constructive interference of RF and LO inputs at the mixer. Finally, the intermediate frequency (IF) mixer product is low-pass filtered before being output into either a multichannel oscilloscope or spectrum analyzer for time and frequency domain measurements respectively. To further improve signal quality we also use various attenuators and filters at room and cryogenic temperatures.

The lower part of Figure 1 shows a SEM image of an $`Al/Al_2O_3`$ twin-SET device used in these experiments. The device was fabricated using electron beam lithography and the standard shadow mask evaporation technique Fulton and Dolan (1987). The tunnel junctions are formed where the leads (horizontal) overlap with the ‘islands’ (vertical). In this arrangement the SETs are $``$500nm apart and control gates used in these experiments couple to both SETs inducing charges of similar magnitude on both devices. In addition there is a central double-dot structure, developed for read out simulations but not relevant to the measurements described here Buehler et al. (2003).

## IV Tank circuit characterization

In order to maximize the sensitivity of the rf-SET, the device should be operated in a regime where the amount of reflected power depends most strongly on the SET resistance. Figures 2a and 2b show the calculated power reflection coefficient $`\mathrm{\Gamma }^2`$ ($`\mathrm{\Gamma }=(Z_0Z_{LCR})/(Z_0+Z_{LCR})`$) in a single rf-SET arrangement as a function of the SET resistance for two different matching networks. For Figure 2a the tank circuit $`LC`$ values were chosen to transform a 17k$`\mathrm{\Omega }`$ impedance to $`Z_0`$ (50$`\mathrm{\Omega }`$) at resonance and for Figure 2b to transform 50k$`\mathrm{\Omega }`$ to $`Z_0`$. For perfect matching the power reflection coefficient then approaches zero when the SET resistance is 17k$`\mathrm{\Omega }`$ for figure 2a and 50k$`\mathrm{\Omega }`$ for figure 2b. We discuss two modes of operation, the under-matched (or under-coupled) regime, where the device resistance is greater than the matching resistance and the over-matched (or over-coupled) regime where the device resistance is smaller than the matching resistance. Over-matching yields an enhanced depth of modulation in association with the large slope of $`\mathrm{\Gamma }^2`$ as a function of the SET resistance, but reduces the $`Q`$ factor Roschier et al. (2004).

Figures 2c and 2d show measurement results for a single rf-SET. The reflected rf power (top) is measured together with the dc conductance (bottom) as a function of gate voltage. The data shown in Figure 2c is limited to the under-matched regime, where the resistance of the SET is always greater than 17k$`\mathrm{\Omega }`$. The corresponding behavior of the reflection coefficient is indicated in Figure 2a by the darker shaded region. In this regime the reflected power increases monotonically with increasing SET resistance, so that the $`S11`$ data shown in Figure 2c is a mirror image of the SET conductance. In contrast, the data shown in Figure 2d was taken with the SET biased to a point where the differential resistance is less than 50k$`\mathrm{\Omega }`$ at the top of a Coulomb blockade peak. Moving from the blockaded state to a point where the differential resistance is less than 50k$`\mathrm{\Omega }`$ results in a non-monotonic response in the reflected rf power as shown in the upper section of Figure 2d. This behavior is due to the impedance transformer operating initially in the under-matched regime and moving to over-matching as the differential resistance moves through 50k$`\mathrm{\Omega }`$, as indicated by the dark shaded region in Figure 2b. The non-monotonic behavior of the transfer function in the over-matched regime complicates operation of the device as a charge detector, as seen in the upper panel of Figure 2d. Over-matching also reduces the $`Q`$-factor so that in a multiplexed arrangement the resonances need to be separated further in frequency to minimize cross-talk.

Turning now to simultaneous operation of both SETs we note that if there was no cross-coupling between the two circuits, then the output signal for tank2 should be independent of the resistance of SET1. Clearly cross-talk is a function of the overlap between the two resonance traces defined by their ‘single circuit’ $`Q`$ factors, so that the resonance frequency of the two circuits should be well separated to minimize cross-talk. In spite of well separated resonances however, additional effects can strongly modify the cross-talk between the two SETs. We now discuss the dominant mode of cross-talk in our setup. In the multiplexed arrangement shown in Figure 1, the branching ‘T’ section behaves as a stub tuning network since at the resonant frequency of SET1, tank2 appears as a load at the end of a length of transmission line in parallel with the transmission feed-line to SET1. The impedance seen looking into the ‘T’ section is Pozar (1998), $`1/Z_{in}=1/Z_1+1/Z_2`$, where

$$Z_{1,2}(\omega _{1,2},d_{1,2})=Z_0\frac{Z_{tank1,2}+jZ_0tan\beta d_{1,2}}{Z_0+jZ_{tank1,2}tan\beta d_{1,2}}$$
(1)
and $`\beta =2\pi /\lambda _{1,2}`$. For the case where the tank circuits are well separated in frequency, with tank2 close to an open circuit for $`\omega _1`$,

$$Z_2=jZ_0cot\beta d_2$$
(2)
so that the impedance $`Z_{1,2}`$ at $`\omega _1`$ looking into the ‘T’ section is just determined by $`d_1`$, $`d_2`$ and $`Z_{tank1}`$. In an ideal arrangement $`d_2`$ should be of zero length, or adjusted so that the impedance $`Z_2`$ appears as an open line. If this condition is achieved the amount of reflected power is solely determined by $`Z_{tank1}`$. In the case where this is not achieved, the stub geometry and length of the transmission lines between the ‘T’ and the SET tank circuits strongly affect both the impedance matching of the $`LC`$ network to the SETs and the degree of cross-talk between the two circuits. Figure 3a and 3b compare the two cases of ‘length matching’ using an EM simulation software package Fulton and Dolan (1987). Optimum length matching is illustrated in Figure 3a where $`S11`$ is plotted as a function of frequency for different values of SET1 resistance, including both under-matching (SET $`>`$ 50k$`\mathrm{\Omega }`$) and over-matching (SET $`<`$ 50k$`\mathrm{\Omega }`$). The resistance of SET2 is held constant at 100k$`\mathrm{\Omega }`$ and inductances $`L_1`$=780nH, $`L_2`$=660nH define the resonances at $`\omega _1`$=339MHz and $`\omega _2`$=370MHz for a parasitic capacitance $`C`$=0.3pF. With transmission line lengths set at $`d_1`$ = 285mm and $`d_2`$ = 275mm, near independent operation of the SETs can be achieved. The inset shows the variation in S11 at $`\omega _2`$ as a function of SET1 resistance. Only in the over-matched case ($`R_{SET1}=30k\mathrm{\Omega }`$) does the resistance of SET1 weakly affect the reflected power at $`\omega _2`$. Turning to the unoptimized case illustrated in Figure 3b ($`d_1`$ = 137mm and $`d_2`$ = 112mm) we firstly note that the variation in line length has altered the resonant frequencies ($`\omega _{1,2}`$) and $`Q`$-factors of the circuits in comparison with the results shown in Figure 3a. Additionally the degree of cross-talk is now greatly increased as shown in the inset to Figure 3b. In particular, for the case where SET1 is operated in the over-matched regime, a poorly tuned stub network can create cross-talk signals as large as 50%.

The above discussion, in which the length of transmission line between the ‘T’ section and the SET tank circuits strongly affects both matching and cross-talk, is valid in the regime where line lengths approach or exceed $`\lambda /4`$, which at these frequencies is the order of centimeters. Although the use of micro-strip or strip-line techniques can be used to avoid this issue, we note that the constraint of independent source-drain biasing for each SET requires the incorporation of a bias-tee in the transmission line in between the ‘T’ junction and SETs. Bias-tees typically contain several inductors of appreciable size and are not easily miniaturized with low insertion loss and good standing wave ratios. Of further interest, a stub network (or double-stub) could be used as an impedance matching technique without the requirement of chip inductors commonly used for impedance matching in the rf-SET. With the addition of a varactor diode terminating the stub, the value of the reactance provided by the stub can be varied with an applied voltage. Such a technique would have several advantages over the present lumped-element matching technique, including in-situ control of the matching and $`Q`$-factor.

We now compare our simulation results with experimental data. Figure 4 shows $`S11`$ data taken with a network analyzer for varying SET2 resistances, $``$43k$`\mathrm{\Omega }`$ (minimum S11 at 364MHz) to $`>`$10M$`\mathrm{\Omega }`$ (maximum S11 at 364MHz). The resistance of SET1 is kept constant throughout (R$`{}_{SET1}{}^{}`$10M$`\mathrm{\Omega }`$, biased in the superconducting gap). Using the microwave simulation package we are able to produce $`S11`$ traces that approximate our experimental results. As fitting parameters we adjust the length of transmission lines $`d_1`$ and $`d_2`$ in the simulation, since these are difficult to measure experimentally at mK temperatures. The values used here ($`d_1`$=313mm, $`d_2`$=304mm) in the simulation results however, are in agreement with what we expect for our setup. In this arrangement the cross-talk between the tank circuits is limited to $`<`$ 5%, which is close to its minimum value for this frequency separation and the SET resistance range. For the simulation results we also tune the length of the bond wire between our circuit board and chip, and the input impedance of our cryogenic microwave amplifier Fulton and Dolan (1987). The amplifier impedance deviates from the ideal case of 50$`\mathrm{\Omega }`$ since it contains an input network to achieve noise matching to the first-stage transistor. Simultaneous noise and impedance matching are generally not achieved, with the result being a reduction in signal level at the output in order to maximize the signal-to-noise ratio Hagen (1996). Further weak deviations between our experimental data and the simulation results are most likely due to an inaccurate frequency dependence of the simulated amplifier input impedance and the variation in parameters characterizing our circuit elements with temperature.

## V Sensitivity and charge noise rejection

We now present measurement results for simultaneous operation of two rf-SETs, beginning with the charge sensitivities. The data shown in Figure 5 was taken with the SETs in the superconducting regime. Operating the twin rf-SET in the superconducting mode maximizes the device sensitivity, since the transconductance and small signal differential conductance are maximized. Figure 5a shows frequency domain data for both rf-SETs in response to a 2.5 MHz sine wave signal applied to a gate with a rms amplitude of $`0.1e`$ on each SET island. Clear amplitude modulation (AM) of both rf carriers is observed. The sensitivity of the rf-SET can now be determined by measuring the signal to noise ratio (SNR) of either sideband for a rms induced charge signal $`q_{rms}`$ on the SET island. The spectral density of the noise is calculated as follows Aassime et al. (2001b):

$$\delta q=\frac{q_{rms}}{\sqrt{B}\times 10^{SNR/20}}(e/\sqrt{Hz})$$
(3)
where $`B`$ is the chosen resolution bandwidth of the spectrum analyzer. For the twin rf-SET device studied here, we measure an optimum sensitivity of $`\delta q=7.5\mu e/\sqrt{Hz}`$ and $`\delta q=4.4\mu e/\sqrt{Hz}`$ for SET1 and SET2 with resonances at 335MHz and 360MHz respectively. The sensitivities were measured for $`q=0.005e`$ signals at 1.3 MHz in the dc+rf (superconducting) mode Aassime et al. (2001b). The frequency dependence of the noise floor shown in Figure 5a is primarily linked to the frequency dependence of the noise temperature of the cryogenic amplifier, which reaches a minimum of $`T_N1.5`$K at $`\omega `$350MHz. Additional (weak) contributions include shot noise associated with the SET current, phase noise of our rf sources, $`1/f`$ noise in the SETs and electronics and the frequency dependent noise of the room temperature amplifiers and mixers.

Focusing on the response of the SETs in the time domain, we now demonstrate the suppression of charge noise via cross-correlated measurements. Figure 5b shows the response of both SETs (measured simultaneously) when a 1 kHz square wave is applied to a nearby gate via a 40 dB attenuator. This signal induces charges of $`0.1e`$ and $`0.05e`$ on the islands of SET1 and SET2 respectively. Although the output of the two SETs faithfully follow the gate signal, additional unwanted charge noise is also picked up by SET2 (e.g. in the 0.7 - 1 ms time interval). The observed random telegraph signals (RTSs) are probably associated with charge noise arising from two-level charge traps in the SET oxide tunnel barrier or substrate. Using a multi-channel digital oscilloscope, we are able to multiply the time derivative signals from SET1 and SET2 in real-time. The result is shown in the lower portion of Figure 5b and constitutes a real-time cross-correlation of the SET signals. Close inspection of the correlation trace shows that this technique produces sharp spikes when there is a true signal that affects both SETs (such as the rising or falling edges of the square wave applied to the gate) and a clear suppression of the random charge noise by up to 95% of the peak height correlation. This ability to reject events associated with charge noise is likely to be important for the application of SETs as detectors and in the readout of solid state quantum computers.

## VI Conclusion

In conclusion we have described the development of a twin rf-SET device that makes use of two tuned impedance transformers to perform wavelength division multiplexing (WDM) using a single cryogenic following amplifier. In conjunction with a mixer-based demodulation technique, we have explored the issues affecting cross-talk and sensitivity. The twin rf-SET has been shown to perform real time cross-correlation of SET signals. Such a correlation technique suppresses charge noise and enables true signals associated with readout to be distinguished from spurious artifacts on time-scales required for solid state quantum computation. Future work will investigate the use of the twin rf-SET for charge motion detection on microsecond time-scales, and its application to readout of both semiconducting and superconducting charge qubits.

We thank D. Barber for technical support and S. Kenyon, K. Lehnert and R. Schoelkopf for fruitful discussions and insights. This work was supported by the Australian Research Council, the Australian Government and by the US National Security Agency (NSA), Advanced Research and Development Activity (ARDA) and the Army Research Office (ARO) under contract number DAAD19-01-1-0653. DJR acknowledges a Hewlett-Packard Fellowship.