format: text

text:
# 1 Introduction

## 1 Introduction

Precise measurements of the top quark properties and interactions are planned at TESLA  and will most certainly belong to the research program of any future $`e^+e^{}`$ collider . The measurements should be confronted with theoretical predictions matching the same high precision level of a few per mille. It is obvious that in order to reach that high precision it is mandatory to include radiative corrections and, as the measurements of some top quark physical properties will be carried out at high energies, much above the $`t\overline{t}`$ threshold, it is crucial to know off-resonance background contributions to any specific 6 fermion decay channel and to estimate the effects related to the off-shellness of the $`t\overline{t}`$-pair.

In this article, a technical documentation of a numerical program, eett6f, is presented which allows for computer simulation of the 6 fermion reactions which are relevant for the top quark pair production and its decay into a specific 6 fermion final state

$$e^+(p_1)+e^{}(p_2)b(p_3)+f_1(p_4)+\overline{f_1^{}}(p_5)+\overline{b}(p_6)+f_2(p_7)+\overline{f_2^{}}(p_8),$$
(1)
where $`f_1=\nu _\mu ,\nu _\tau ,u,c`$, $`f_2=\mu ^{},\tau ^{},d,s`$, and $`f_1^{}`$, $`f_2^{}`$ are the corresponding weak isospin partners of $`f_1`$, $`f_2`$, $`f_1^{}=\mu ^{},\tau ^{},d,s`$, $`f_2^{}=\nu _\mu ,\nu _\tau ,u,c`$, and the particle four momenta have been indicated in the parentheses. In the present version, v. 1.0, of the program, reactions (1) are treated in the lowest order of SM. Moreover, for the sake of simplicity, it is assumed that the actual values of $`f_1`$ and $`f_2^{}`$ in (1) are different from each other, and that neither $`f_1^{}`$ nor $`f_2`$ is an electron. The physics contents of the program is described in , where the details on the physical model, a method of calculation and physical results obtained with eett6f v1.0 have been discussed.

A similar analysis of the 6 fermion processes relevant for a $`t\overline{t}`$ production in $`e^+e^{}`$ annihilation have been performed in , , where semileptonic channels of reaction (1) have been studied, and in , where purely hadronic channels of (1) have been analysed. Moreover, irreducible QCD background to top searches in semileptonic channels of (1) has been discussed in . After this work has already been completed an exstensive study of six fermion reactions, including (1) among others, in the massless fermion limit has appeared . Using the same input parameters and separation cuts as in , the present program gives 5.8160(32) fb and 17.223(15) fb for the cross sections of $`e^+e^{}b\nu _\mu \mu ^+\overline{b}\tau ^{}\overline{\nu }_\tau `$ and $`e^+e^{}bu\overline{d}\overline{b}\mu ^{}\overline{\nu }_\mu `$, respectively, at CMS energy of 500 GeV. These results are 2–3 standard deviations off the corresponding results of . This difference can most probably be traced back to somewhat different implementation of the top quark width in the top quark propagator. In the present work, the top quark width is introduced both in the denominator and in the numerator of the top propagator while, in , it seems to be introduced in the denominator only. The invariant mass and angular distributions of the $`bu\overline{d}`$ quark triple of $`e^+e^{}bu\overline{d}\overline{b}\mu ^{}\overline{\nu }_\mu `$ at $`s^{1/2}=500`$ GeV shown in  are nicely reproduced by eett6f within accuracy of the plots.

The main advantage of the present program is that it allows for taking into account both the electroweak and QCD lowest order contributions. As light fermion masses are not neglected, the cross sections can be calculated without any kinematical cuts. Moreover, a number of options have been implemented in the program which make possible calculation of the cross sections while switching on and off different subsets of the Feynman diagrams. It is possible to calculate cross sections in two different narrow width approximations discussed below, too.

Besides of a possibility of taking into account solely the electroweak contributions, or switching off the Higgs boson exchange diagrams, the program allows also for a simplified treatment of reaction (1) by utilizing a few different approximations: the double resonance approximation for $`W`$ bosons

$`e^+(p_1)+e^{}(p_2)`$ $``$ $`b(p_3)+W_{}^{+}{}_{}{}^{}(p_{45})+\overline{b}(p_6)+W_{}^{}{}_{}{}^{}(p_{78})`$ (2)
$``$ $`b(p_3)+f_1(p_4)+\overline{f_1^{}}(p_5)+\overline{b}(p_6)+f_2(p_7)+\overline{f_2^{}}(p_8),`$ (3)
where only those 61 Feynman diagrams are taken into account which contribute to $`e^+e^{}bW^+\overline{b}W^{}`$ and the W bosons are considered as being off-mass-shell and the double resonance approximation for a $`t`$\- and $`\overline{t}`$-quark

$`e^+(p_1)+e^{}(p_2)`$ $``$ $`t^{}(p_{345})+\overline{t}^{}(p_{678})`$ (4)
$``$ $`b(p_3)+f_1(p_4)+\overline{f_1^{}}(p_5)+\overline{b}(p_6)+f_2(p_7)+\overline{f_2^{}}(p_8),`$ (5)
with only two ‘signal’ diagrams contributing. The intermediate state momenta are $`p_{45}=p_4+p_5`$, $`p_{78}=p_7+p_8`$ in Eq. (2) and $`p_{345}=p_3+p_4+p_5`$, $`p_{678}=p_6+p_7+p_8`$ in Eq. (4). The narrow width approximation for the $`W`$ bosons (top quarks) is obtained by replacing the intermediate state $`W`$ bosons (top quarks) in Eq. (2) (in Eq. (4)) with on mass shell $`W`$’s (top quarks) and by multiplying the cross section of the on mass shell reaction (2) $`\left((\text{4})\right)`$ with the branching ratios corresponding to final state of reaction (1). The necessary partial decay width of the $`W`$ and $`t`$ are calculated according to the lowest order SM, and the fermion masses are neglected in the three body top decay with.

The necessary matrix elements of reactions (1), (2) and (4) are calculated with the helicity amplitude method described in  and  and phase space integrations are performed with the Monte Carlo (MC) method. More details on the multi-channel MC algorithm (see e.g. ), on which the MC integration and event generation is based, are given in the next section.

## 2 Monte Carlo integration and event generation

In order to improve the convergence of the MC integration the most relevant peaks of the matrix element squared related to the Breit-Wigner shape of the $`W,Z`$, Higgs and top quark resonances as well as to the exchange of a massless photon or gluon have to be mapped away. As it is not possible to find out a single parametrization of the $`n`$-dimensional phase space which would allow to cover the whole resonance structure of the integrand, it is necessary to utilize a multi-channel MC approach. In the multi-channel MC approach, the random numbers $`0x_j1`$, $`j=1,\mathrm{},n`$, and hence four momenta of the final state particles, are generated according to one of $`N`$ different probability densities $`f_i(x)`$, $`i=1,\mathrm{},N`$, $`x=(x_1,\mathrm{},x_n)`$. The distribution $`f_i(x)`$, which accounts for several different peaks of the integrand, is selected with probability $`a_i`$. All the distributions $`f_i(x)`$ are then combined in a single probability distribution $`f(x)`$ which should cover possibly all the peaks of the integrand, or at least the most relevant ones

$`f(x)={\displaystyle \underset{i=1}{\overset{N}{}}}a_if_i(x),\mathrm{with}{\displaystyle \underset{i=1}{\overset{N}{}}}a_i=1.`$ (6)
The second condition of Eq. (6) guaranties that, if every distribution $`f_i(x)`$ is normalized to unity, then the combined distribution $`f(x)`$ is normalized to unity, too. Obviously, all the distributions $`f_i(x)`$ and weights $`a_i`$, $`i=1,\mathrm{},N`$, are to be non negative.

The basic phase space parametrizations which are used in the program are listed below. The 6 particle phase space of reaction (1) is parametrized in 3 different ways:

$`\mathrm{d}^{14}Lips`$ $`=`$ $`(2\pi )^{14}\mathrm{d}PS_2(s,s_{345},s_{678})\mathrm{d}PS_2(s_{345},m_3^2,s_{45})\mathrm{d}PS_2(s_{678},m_6^2,s_{78})`$ (7)
$`\times `$ $`\mathrm{d}PS_2(s_{45},m_4^2,m_5^2)\mathrm{d}PS_2(s_{78},m_7^2,m_8^2)\mathrm{d}s_{345}\mathrm{d}s_{678}\mathrm{d}s_{45}\mathrm{d}s_{78},`$
$`\mathrm{d}^{14}Lips`$ $`=`$ $`(2\pi )^{14}\mathrm{d}PS_2(s,s_{34},s_{5678})\mathrm{d}PS_2(s_{5678},s_{56},s_{78})\mathrm{d}PS_2(s_{34},m_3^2,m_4^2)`$ (8)
$`\times `$ $`\mathrm{d}PS_2(s_{56},m_5^2,m_6^2)\mathrm{d}PS_2(s_{78},m_7^2,m_8^2)\mathrm{d}s_{34}\mathrm{d}s_{5678}\mathrm{d}s_{56}\mathrm{d}s_{78}`$
and

$`\mathrm{d}^{14}Lips`$ $`=`$ $`(2\pi )^{14}\mathrm{d}PS_2(s,m_3^2,s_{45678})\mathrm{d}PS_2(s_{45678},s_{45},s_{678})\mathrm{d}PS_2(s_{678},m_6^2,s_{78})`$ (9)
$`\times `$ $`\mathrm{d}PS_2(s_{45},m_4^2,m_5^2)\mathrm{d}PS_2(s_{78},m_7^2,m_8^2)\mathrm{d}s_{45678}\mathrm{d}s_{45}\mathrm{d}s_{678}\mathrm{d}s_{78}.`$
The 4 particle phase space of reaction (2) in the narrow $`W`$ width approximation is parametrized in 2 different ways:

$$\mathrm{d}^8Lips=(2\pi )^8\mathrm{d}PS_2(s,m_3^2,s_{456})\mathrm{d}PS_2(s_{456},m_4^2,s_{56})\mathrm{d}PS_2(s_{56},m_5^2,m_6^2)\mathrm{d}s_{456}\mathrm{d}s_{56}$$
(10)
and

$$\mathrm{d}^8Lips=(2\pi )^8\mathrm{d}PS_2(s,s_{34},s_{56})\mathrm{d}PS_2(s_{34},m_3^2,m_4^2,)\mathrm{d}PS_2(s_{56},m_5^2,m_6^2)\mathrm{d}s_{34}\mathrm{d}s_{56}.$$
(11)
In Eqs. (711), $`s_{ijk\mathrm{}}=(p_i+p_j+p_k+\mathrm{})^2,i,j,k=3,\mathrm{},8`$, and $`\mathrm{d}PS_2(s,s^{},s^{\prime \prime })`$ is a two particle (subsystem) phase space element defined by

$$\mathrm{d}PS_2(s,s^{},s^{\prime \prime })=\delta ^4\left(pp^{}p^{\prime \prime }\right)\frac{\mathrm{d}^3p^{}}{2E^{}}\frac{\mathrm{d}^3p^{\prime \prime }}{2E^{\prime \prime }}=\frac{|\stackrel{}{p}^{}|}{4\sqrt{s}}\mathrm{d}\mathrm{\Omega }^{},$$
(12)
where $`\stackrel{}{p}^{}`$ is the momentum and $`\mathrm{\Omega }^{}`$ is the solid angle of one of the particles (subsystems) in the relative centre of mass system, $`\stackrel{}{p}^{}+\stackrel{}{p}^{\prime \prime }=0`$. Making use of the rotational symmetry with respect to the $`e^+e^{}`$ beam axis reduces the dimension of the phase space elements to 13 in Eqs. (7)–(9) and to 7 in Eqs. (10) and (11).

Parametrizations (711) are used with different permutations of external particle momenta such that invariants $`s_{ijk\mathrm{}}`$ possibly correspond to the virtuality of propagators of the gauge bosons, Higgs boson and/or top quarks. The invariants are then transformed to the interval $`[0,1]`$ by performing mappings smoothing out peaks related to the propagators. All other integration variables of Eqs. (711), not related to peaks of the matrix element squared, are also transformed to the interval $`[0,1]`$ with a simple mapping defined in the following way. Let $`y_j[c_j,b_j]`$ be a variable not related to a peak in the matrix element squared and $`x_j`$ be random variable uniformly distributed in the interval $`[0,1]`$, then the relation

$$y_j=\left(b_jc_j\right)x_j+c_j$$
(13)
defines the necessary mapping. After having all the integration variables transformed to the interval $`[0,1]`$, any of phase space parametrizations (79) for reaction (1), or (10), (11) for reaction (2) in the narrow width approximation, for a given permutation of particle momenta, can be directly associated with the probability density $`f_i(x)`$, which is referred to as a kinematical channel. Altogether 59 (21) different channels are used in order to integrate the matrix element squared of reaction (1) (reaction (2) in the narrow width approximation).

The weights $`a_i`$ are calculated in the initial scanning run that starts with all the weights equal to each other. The weights $`a_i`$, $`i=1,\mathrm{},N`$, for the actual run are calculated as ratios

$$a_i=\sigma _i/\underset{j=1}{\overset{N}{}}\sigma _j,$$
(14)
where $`\sigma _j`$ denotes the cross section obtained with the $`i`$-th kinematical channel in the initial scan.

Calculation of the cross section of reaction (4) in the narrow top width approximation does not require the multi-channel MC approach and can be performed with the single phase space parametrization (12).

## 3 Description of the program

eett6f is a package written in FORTRAN 90. It consists of 50 files including a makefile. They are stored in one working directory. The user should specify the physical input parameters in module inprms.f and select a number of options in the main program csee6f.f.

### 3.1 Program input

The default values of the input parameters and options used in the program are those specified below.

#### 3.1.1 Physical parameters

The initial physical parameters to be specified in module inprms.f are the following.

Gauge boson masses and widths (GeV):
mw=80.419 GeV, the W mass,
gamw=2.12 GeV, the W width,
mz=91.1882 GeV, the Z mass,
gamz=2.4952 GeV, the Z width.
The electroweak (EW) mixing parameter sw2 is then calculated from

sw2 = 1-mw2/mz2

utilizing mw2 = mw\**2, mz2 = mz\**2 in the fixed width scheme, and
mw2 = mw\**2 - i\*mw\*gamw, mz2 = mz\**2 - i\*mz\*gamz in the complex mass scheme.
ralp0 = 137.03599976, an inverse of the fine structure constant in the Thomson limit,
gmu=1.16639$`\times \mathrm{𝟷𝟶}^\mathrm{𝟻}`$ GeV<sup>-2</sup>, the Fermi coupling constant,
alphas=0.1185, the strong coupling constant at mz.
An inverse of the fine structure constant at mw, ralpw, is then calculated from

ralpw=4.44288293815837/(2\*sw2\*gmu\*mw\**2).

mh=115.0 GeV, the Higgs boson mass,
gamh=0 GeV, the Higgs boson width.
If gamh = 0, then the Higgs width is calculated according to the lowest order of SM.
Fermion masses and widths:
me = 0.510998902 MeV, game = 0 MeV, for an electron,
mmu = 105.658357 MeV, gammu = 0 MeV, for a muon,
mtau=1.77703 GeV, gamtau = 0 GeV, for a lepton tau,
mu = 5 MeV, gamu = 0 MeV, for an up quark,
md=9 MeV, gamd = 0 MeV, for a down quark,
mc=1.3 GeV, gamc = 0 GeV, for a charm quark,
ms=150 MeV, gams = 0 MeV, for a strange quark,
mt=174.3 GeV, gamt = 1.5 GeV, for a top quark,
mb=4.4 GeV, gamb = 0 GeV, for a bottom quark.
ncol=3, the number of colours,
convc=0.389379292$`\times 10^{12}`$ fb $`\mathrm{GeV}^2`$, a conversion constant.

#### 3.1.2 Main options

The following main options should be selected in the main program csee6f.f.
The number of different centre of mass (CMS) energies ne
ne = 1. Recommended if the unweighted event are to be generated.
The actual values of the CMS energies in the array aecm of size ne:
aecm=(/500.d0/).
The final state of (1) by selecting a decay mode of the intermediate $`W^+`$, $`W^{}`$ bosons, iwp, iwm = 1 or 2, where 1 corresponds to a leptonic, and 2 to a hadronic decay mode, and by specifying the family index ifp,ifm = 1,2,3 of the fermion pair resulting from the $`W^+`$ or $`W^{}`$ decay. For example, the reaction $`e^+e^{}b\nu _\mu \mu ^+\overline{b}d\overline{u}`$ corresponds to $`W^+`$ decaying leptonically and $`W^{}`$ decaying hadronically and the corresponding flags are:
iwp = 1
ifp = 2
iwm = 2
ifm = 1.
Whether or not to calculate the Born cross section utilizing the full 6 particle kinematics, iborn = 1 (yes) / else (no), with ncall0 calls to the integrand
iborn = 1
ncall0 = 20000. Recommended No of calls is a few millions.
Generate the unweighted events or not, imc = 1(yes)/else(no)?
imc = 0.
No standard event record is used. If imc = 1, then the final state particle momenta of the accepted unweighted events are printed in the output.
Scan the Born cross section with nscan0 calls, iscan0 = 1(yes)/else(no), in order to find the dominant kinematical channels, adjust integration weights and find out the maximum value of the cross section
iscan0 = 1. This option is strongly recommended.
nscan0 = 200. About 1 thousand is recommended.
Calculate an approximate cross section of (2), iwwa = 1 (yes) / else (no)?
iwwa = 1
ncallww = 20000. Recommended value is several hundred thousands.
Scan the Born cross section in the narrow W-width approximation with nscanww calls, iscanww = 1(yes)/else(no)?
iscanww = 1. Recommended.
nscanww = 200. Recommended value is about 1 thousand.
Calculate an approximate cross section of (4), itopa = 1(yes)/else(no)?
itopa = 1
ncalla = 20000.
Include the Higgs boson exchange, ihiggs = 1(yes)/else(no)?
ihiggs = 1
Calculate the electroweak contributions only, iew = 1(yes)/else(electroweak + QCD
Feynman diagrams)?
iew = 0.
Calculate the $`t\overline{t}`$ resonant diagrams only, itt = 1(yes)/else(all the diagrams)?
itt = 0.
Calculate the double W resonant diagrams only, iww = 1(yes)/else(all the diagrams)?
iww = 0.
Choose the scheme: ischeme = 1(complex mass scheme)/else(fixed width scheme)
ischeme = 1.
If ischeme=1, then should alpha\_W be complex (iaplw=1) or real (ialpw=0)?
ialpw=1.

#### 3.1.3 Auxiliary options

Impose cuts, icuts = 1(yes)/else(no)?
icuts = 0
If icuts = 1, then specify the kinematical cuts, e.g.:
ctlb = cos(10\*pi/180) – cosine of the charged lepton–beam angle,
ctqb = cos(5\*pi/180) – cosine of the quark–beam angle,
ctll = 1 – cosine of the charged lepton–charged lepton angle,
ctlq = cos(5\*pi/180) – cosine of the charged lepton–quark angle,
ecutl = 1 GeV – minimum charged lepton energy,
ecutq = 3 GeV – minimum quark energy,
mqq = 10 GeV – minimum invariant mass of a quark pair,
mll = 0 GeV – minimum invariant mass of a charged lepton pair,
Calculate distributions, idist = 1(yes)/0(no)?
idis=0.
If idis = 1, then specify parameters of the distributions:
xmin = (/ xl1, xl2, xl3, xl4/) – lower bounds,
xmax = (/ xu1, xu2, xu3, xu4/) – upper bounds,
nbs = (/ n1, n2, n3, n4/) – the corresponding numbers of bins in each distribution.
Constants xli, xui should be of type real(8) and ni of type integer. The number of desired distributions and the maximum number of bins, nbmax = max{n1, n2, n3, n4} should be specified in module distribs.f.
The maximum of the fully differential cross section crmax, relevant only if iscan0 = 0 or iscanww = 0,
crmax=1000.

### 3.2 Routines of eett6f

The main program csee6f.f, each subroutine, function or module are located in a file named the same way as the routine itself, except for 3 functions: srr, src and scc located in a single file dotprod.f. The program consists of the following routines.

* SUBROUTINE boost – returns a four vector boosted to the Lorentz frame of velocity –v.
* SUBROUTINE carlos – the MC integration routine.
* SUBROUTINE couplsmb – returns the SM couplings.
* FUNCTION cross – calculates the cross section of (1), approximated cross sections of (25) and the cross section of (2) in the narrow $`W`$-width approximation.
* FUNCTION crosstopa – calculates the cross section of (4) in the narrow top width approximation.
* The MAIN PROGRAM csee6f.
* MODULE distribs – contains parameters of distributions.
* MODULE drivec – contains driving flags and some kinematical variables.
* SUBROUTINE eebwbw – returns the squared matrix element averaged over initial spins and summed over final spins and colours of $`e^+e^{}bW^+\overline{b}W^{}`$.
* SUBROUTINE eee – returns a contraction of a triple gauge boson coupling with three complex four vectors.
* SUBROUTINE eeee – returns a contraction of a quartic gauge boson coupling with four complex four vectors.
* SUBROUTINE eeff1 – returns the squared matrix element averaged over initial spins and summed over final spins and colours of $`e^+e^{}f\overline{f}`$.
* SUBROUTINE eeh – returns a contraction of the Higgs-gauge boson coupling with two complex four vectors.
* SUBROUTINE eett6f – returns the squared matrix element averaged over initial spins and summed over final spins and colours of reaction (1) and approximations (25).
* SUBROUTINES eev, eve and vee – return contractions of a triple gauge boson coupling with two polarization vectors leaving, respectively, the third, second and first Lorentz index uncontracted.
* SUBROUTINE fhfb – returns a set complex scalars: $`scal(\lambda _1,\lambda _2)=g_{S12}\overline{u_1}u_2\mathrm{\Delta }_F`$, where $`u_1`$, $`u_2`$ are fermion spinors, $`\mathrm{\Delta }_F`$ is a scalar boson propagator and $`g_{S12}`$ is a coupling of the scalar boson to fermion 1 and 2.
* SUBROUTINE fvfa – returns a set of four vectors $`\epsilon ^\mu (\lambda _1,\lambda _2)=\overline{u}_1(\lambda _1)\gamma _\nu (g_V^{()}P_{}+g_V^{(+)}P_+)`$ $`\times u_2(\lambda _2)D_V^{\nu \mu }`$, where $`u_1`$, $`u_2`$ are fermion spinors, $`P_\pm =(1\pm \gamma _5)/2`$, are chirality projectors, $`g_V^{(\pm )}`$ are the fermion-gauge boson couplings of definite chirality and $`D_V^{\nu \mu }`$ is the propagator of a gauge boson $`V`$ in arbitrary linear gauge.
* SUBROUTINE fef – returns a set of matrix elements $`mat(\lambda _1,\lambda ,\lambda _2)=\overline{u_1}(\lambda _1)/\epsilon (\lambda )(g_V^{()}P_{}+g_V^{(+)}P_+)u_2(\lambda _2)`$, where notation is the same as in fvfa.
* SUBROUTINE impcuts – imposes cuts on the kinematical variables of process (1).
* MODULE inprms – contains initial input parameters.
* MODULE kincuts – contains kinematical cuts.
* SUBROUTINE kinee1212 – returns the four momenta, phase space normalization and flux factor for a 2 $``$ 6 process in the CMS. The phase space parametrization used is that of Eq. (9).
* SUBROUTINE kinee6f – returns the four momenta, phase space normalization and flux factor for a 2 $``$ 6 process in the CMS. The phase space parametrization used is that of Eq. (8).
* SUBROUTINE kineeff1 – returns the four momenta, phase space normalization and flux factor for a $`22`$ process in the CMS. The phase space is parametrized according to Eq. (12).
* SUBROUTINE kineeffff – returns the four momenta, phase space normalization and flux factor for a 2 $``$ 4 process in the CMS. The phase space is parametrized according to Eq. (11).
* SUBROUTINE kineetbff – returns the four momenta, phase space normalization and flux factor for a 2 $``$ 4 process in the CMS. The phase space is parametrized according to Eq. (10).
* SUBROUTINE kineett6f – returns the four momenta, phase space normalization and flux factor for a 2 $``$ 6 process in the CMS. The phase space parametrization used is that of Eq. (7).
* SUBROUTINE kinff – returns the final state four momenta and phase space normalization for a $`22`$ process in the CMS.
* FUNCTION lamsq – the kinematic lambda function, $`\lambda (\sqrt{x},\sqrt{y},\sqrt{z})`$.
* MODULE mathprms – contains necessary arithmetical constants.
* SUBROUTINE parfix – returns parameters of a specific process (1) and the branching ratios for the narrow width approximations; initializes weights for the Monte Carlo integration.
* MODULE parproc – contains parameters for a specific process.
* SUBROUTINE peu – returns a set of generalized spinors: $`u(\lambda ,\lambda _1)=S_F/\epsilon (\lambda )(g_V^{()}P_{}+g_V^{(+)}P_+)u_1(\lambda _1)`$, where $`S_F`$ is a Feynman propagator of an internal fermion and the remaining notation is the same as in fvfa.
* SUBROUTINE psnee1212 – returns a phase space normalization as the one of kinee1212 for a given set of external particle momenta.
* SUBROUTINE psnee6f – returns a phase space normalization as the one of kinee6f for a given set of external particle momenta.
* SUBROUTINE psneeffff – returns a phase space normalization as the one of kineeffff for a given set of external particle momenta.
* SUBROUTINE psneetbff – returns a phase space normalization as the one of kineetbff for a given set of external particle momenta.
* SUBROUTINE psneett6f – returns a phase space normalization as the one kineett6f for a given set of external particle momenta.
* SUBROUTINE recpol – returns real polarization vectors of a vector boson in the rectangular basis.
* FUNCTION scc (contained in dotprod.f) – returns the Minkowski dot product of two complex four vectors.
* SUBROUTINE spheric – returns spherical components ps of a four vector $`p^\mu `$, ps$`=(|𝐩|,\mathrm{cos}\theta ,\mathrm{sin}\theta ,\mathrm{cos}\varphi ,\mathrm{sin}\varphi )`$, with $`\theta `$ and $`\varphi `$ being a polar and azimuthal angles of momentum p.
* SUBROUTINE spinc – returns the contractions: $`p^0I𝐩\sigma `$ and $`p^0I+𝐩\sigma `$, where $`p^\mu =(p^0,𝐩)`$ is a complex four vector, $`I`$ is the $`2\times 2`$ unit matrix and $`\sigma `$ are the Pauli matrices.
* SUBROUTINE spinornew – returns helicity spinors in the Weyl representation.
* SUBROUTINE spinr – returns the same contractions as in spinc for a real four vector $`p^\mu `$.
* FUNCTION src (contained in dotprod.f) – returns the Minkowski dot product of a real and a complex four vector.
* FUNCTION srr (contained in dotprod.f) – returns the Minkowski dot product of two real four vectors.
* SUBROUTINE uep – returns a set of generalized spinors: $`u(\lambda _1,\lambda )=u_1(\lambda _1)/\epsilon (\lambda )(g_V^{()}P_{}+g_V^{(+)}P_+)S_F`$, where notation is the same as in peu.
* SUBROUTINE wwidth – returns the partial width of the $`W`$-boson averaged over initial spins and summed over final spins, and colours for hadronic decay modes.

### 3.3 Run output

A sample of the listing of the test run output is given in Appendix A. It contains a specification of the process, information on the scheme choice, and values of the relevant physical parameters in the very beginning. Then, for each value of the CMS energy and for each cross section calculated, results of the initial scan and resulting weights for the actual run are printed, if iscan = 1. Sometimes the weights calculated according to Eq. (14) do not add up exactly to 1. Then, the last weight is changed slightly in order the second condition of Eq. (6) to be fulfilled. This change is supposed to be completely irrelevant numerically, but the message about it is printed. The final result for the total cross section is called Integral. Its value is printed together with the standard deviation and the actual number of calls used in the calculation. Finally, information on events acceptance efficiency is given, that means a fraction of accepted weight 1 events. If the event generation option is switched on, i.e. imc = 1, then the corresponding number of unwighted events will be printed out as collections of final state particle four momenta.

Whenever a maximum value of the cross section initially assigned in csee6f.f or found in the result of the initial scan is overflown by more than a factor 1.5, a corresponding message informing about it is printed. The integration is still valid, however, if the program is run as an event generator, it should be rerun. How to proceed in this situation is described in the next section. In the very end, all the calculated total cross sections together with the corresponding standard deviations are collected in the tabular form.

## 4 Use of the program

Up to now the program have been run only on Unix/Linux platforms. In order to run the program, the user should select a specific name of a FORTRAN 90/95 compiler, desired options and the name of the output file in the makefile. The output file is called test at present. The program can be then compiled, linked and run by executing a single command
make test.
The results of the test run should reproduce those contained in file testrun.

The program can be run as the MC event generator of unweighted events by selecting
imc = 1
in csee6f.f. It is then recommended to run the program for a single CMS energy and to perform a scan with a relatively large number of calls, nscan0 or nscanww. Attention should be paid to possible messages informing about updates of the maximum weight. In this case the program should be rerun, however, this time without the initial scan. The initial integration weights aw0 or awa in parfix.f and the maximum weight crmax in csee6f.f should be updated according to the results of the prior scan. At present the efficiency of events acceptance is relatively low. However, in view of a relatively fast performance of the program, this should not be a serious limitation for the user.

## 5 Outlook

The following improvements of the program are envisaged in the near future. First of all, a new version is being prepared which will allow for calculations of reactions involving the same external fermion flavours. Another step will be an inclusion of the most relevant factorizable higher order effects: electroweak and QCD corrections to the top pair production, $`e^+e^{}t\overline{t}`$, to the $`t`$ and $`\overline{t}`$ quark decays, as well as to the $`W^+`$ and $`W^{}`$ decays.

Appendix

## Appendix A Listing of test run output

------------------------------------------------------------------------------

Approaching the process:

e+ e- -\> b vm mu+ b d u

with the program ’eett6f v. 1.0’ by

Karol Kolodziej

University of Silesia

Katowice, Poland

e-mail: kolodzie@us.edu.pl

------------------------------------------------------------------------------

Complex mass scheme
with sw2 = ( 0.222271237327414,-0.000778147496385)
and complex 1/alpha\_W = (132.386398893114400, 3.953796222427256)
Gauge boson masses and widths: W, Gm\_W, Z, Gm\_Z
80.41900 2.12000 91.18820 2.49520
Higgs boson mass and width: H, Gm\_H
0.11500E+03 0.49657E-02
1/alpha0 = 137.0359998 alphas = 0.1185000 G\_mu = 0.116639E-04
t-quark mass and width:
174.30 1.50000
External fermion masses: e- b vm mu+ d u :
0.51099890E-03 0.44000000E+01 0.00000000E+00 0.10565836E+00
0.90000000E-02 0.50000000E-02
------------------------------------------------------------------------------

E\_cm = 500.000 GeV

------------------------------------------------------------------------------
All cross sections in femtobarns (fb)
------------------------------------------------------------------------------
The Born cross section with full 6 particle kinematics:
------------------------------------------------------------------------------
Results of the initial scan:
Kinematics Weight Integral St.dev. No of calls
1 0.17241379E-01 0.77943026E+01 0.53027331E+00 200
2 0.17241379E-01 0.70433588E+00 0.14557680E+00 200
3 0.17241379E-01 0.46787757E+00 0.14942639E+00 200
4 0.17241379E-01 0.40042444E-01 0.17572985E-01 200
5 0.17241379E-01 0.61761884E-01 0.17401289E-01 200
6 0.17241379E-01 0.70179588E-02 0.50890185E-02 200
7 0.17241379E-01 0.18918195E-04 0.14788334E-04 200
8 0.17241379E-01 0.59135660E-02 0.26123591E-02 200
9 0.17241379E-01 0.65946999E-05 0.37561359E-05 200
10 0.17241379E-01 0.71059323E-02 0.32477486E-02 200
11 0.17241379E-01 0.12331601E+00 0.52550633E-01 200
12 0.17241379E-01 0.50746197E-01 0.85455398E-02 200
13 0.17241379E-01 0.52341207E-01 0.13094625E-01 200
14 0.17241379E-01 0.16697241E-02 0.91218523E-03 200
15 0.17241379E-01 0.21727387E-02 0.96316906E-03 200
16 0.17241379E-01 0.10069801E+00 0.44469111E-01 200
17 0.17241379E-01 0.55418681E-01 0.10461875E-01 200
18 0.17241379E-01 0.47690194E-03 0.25400283E-03 200
19 0.17241379E-01 0.12270501E-05 0.68571740E-06 200
20 0.17241379E-01 0.57639038E-02 0.20367189E-02 200
21 0.17241379E-01 0.10948697E-03 0.78871872E-04 200
22 0.17241379E-01 0.38725698E-03 0.32739976E-03 200
23 0.17241379E-01 0.31913662E-03 0.31771241E-03 200
24 0.17241379E-01 0.18057468E-04 0.15825867E-04 200
25 0.17241379E-01 0.75726379E-02 0.10353862E-02 200
26 0.17241379E-01 0.40524491E-01 0.40417673E-01 200
27 0.17241379E-01 0.48346637E-03 0.25434560E-03 200
28 0.17241379E-01 0.17689713E-05 0.10354945E-05 200
29 0.17241379E-01 0.24943550E-03 0.85703152E-04 200
30 0.17241379E-01 0.68130329E+01 0.55624830E+00 200
31 0.17241379E-01 0.80133941E+00 0.14035567E+00 200
32 0.17241379E-01 0.37030192E+00 0.13116529E+00 200
33 0.17241379E-01 0.27961633E-01 0.79658869E-02 200
34 0.17241379E-01 0.12381885E-02 0.65639060E-03 200
35 0.17241379E-01 0.38608372E-05 0.23584221E-05 200
36 0.17241379E-01 0.22584220E-03 0.16055018E-03 200
37 0.17241379E-01 0.10262033E-05 0.42675278E-06 200
38 0.17241379E-01 0.59892263E-04 0.17034981E-04 200
39 0.17241379E-01 0.11805922E-04 0.96111792E-05 200
40 0.17241379E-01 0.37482281E-02 0.13933930E-02 200
41 0.17241379E-01 0.69961565E-02 0.29473064E-02 200
42 0.17241379E-01 0.26317460E-03 0.96401737E-04 200
43 0.17241379E-01 0.14981430E-04 0.14211747E-04 200
44 0.17241379E-01 0.36580810E-02 0.31164792E-02 200
45 0.17241379E-01 0.65694693E-05 0.59639961E-05 200
46 0.17241379E-01 0.59283100E-01 0.19915949E-01 200
47 0.17241379E-01 0.34966532E-02 0.32846811E-02 200
48 0.17241379E-01 0.39467871E-01 0.24414577E-01 200
49 0.17241379E-01 0.33496371E-02 0.76289358E-03 200
50 0.17241379E-01 0.65454536E+00 0.11507542E+00 200
51 0.17241379E-01 0.11775709E-01 0.58766647E-02 200
52 0.17241379E-01 0.12389848E-01 0.28548786E-02 200
53 0.17241379E-01 0.19047822E-01 0.13103839E-01 200
54 0.17241379E-01 0.29502206E-01 0.16858165E-01 200
55 0.17241379E-01 0.62510122E+00 0.13584265E+00 200
56 0.17241379E-01 0.82405870E-02 0.22527034E-02 200
57 0.17241379E-01 0.55814360E-02 0.23231953E-02 200
58 0.17241379E-01 0.78977140E-02 0.36720050E-02 200
------------------------------------------------------------------------------
Born scanned with 11565 calls = 0.19039197E+02
The maximum value of the cross section= 0.27872147E+04
shifted to 0.41808220E+04
------------------------------------------------------------------------------
Weight 58 = 0.4148134056762879E-03 rescaled to 0.4148134056756181E-03
Calculating the Born cross section with the weights:
1 0.40938191E+00 0.40938191E+00
2 0.36993992E-01 0.44637590E+00
3 0.24574439E-01 0.47095034E+00
4 0.21031583E-02 0.47305350E+00
5 0.32439333E-02 0.47629743E+00
6 0.36860583E-03 0.47666604E+00
7 0.99364459E-06 0.47666703E+00
8 0.31059956E-03 0.47697763E+00
9 0.34637490E-06 0.47697798E+00
10 0.37322648E-03 0.47735120E+00
11 0.64769547E-02 0.48382816E+00
12 0.26653539E-02 0.48649351E+00
13 0.27491290E-02 0.48924264E+00
14 0.87699297E-04 0.48933034E+00
15 0.11411925E-03 0.48944446E+00
16 0.52889841E-02 0.49473344E+00
17 0.29107678E-02 0.49764421E+00
18 0.25048428E-04 0.49766926E+00
19 0.64448631E-07 0.49766932E+00
20 0.30273882E-03 0.49797206E+00
21 0.57506089E-05 0.49797781E+00
22 0.20339986E-04 0.49799815E+00
23 0.16762084E-04 0.49801491E+00
24 0.94843646E-06 0.49801586E+00
25 0.39773936E-03 0.49841360E+00
26 0.21284769E-02 0.50054208E+00
27 0.25393213E-04 0.50056747E+00
28 0.92912076E-07 0.50056757E+00
29 0.13101157E-04 0.50058067E+00
30 0.35784246E+00 0.85842313E+00
31 0.42088930E-01 0.90051206E+00
32 0.19449451E-01 0.91996151E+00
33 0.14686352E-02 0.92143014E+00
34 0.65033652E-04 0.92149518E+00
35 0.20278362E-06 0.92149538E+00
36 0.11861960E-04 0.92150724E+00
37 0.53899508E-07 0.92150729E+00
38 0.31457348E-05 0.92151044E+00
39 0.62008511E-06 0.92151106E+00
40 0.19686903E-03 0.92170793E+00
41 0.36746070E-03 0.92207539E+00
42 0.13822779E-04 0.92208921E+00
43 0.78687301E-06 0.92209000E+00
44 0.19213421E-03 0.92228213E+00
45 0.34504971E-06 0.92228248E+00
46 0.31137396E-02 0.92539622E+00
47 0.18365550E-03 0.92557987E+00
48 0.20729799E-02 0.92765285E+00
49 0.17593374E-03 0.92782879E+00
50 0.34378833E-01 0.96220762E+00
51 0.61849819E-03 0.96282612E+00
52 0.65075477E-03 0.96347687E+00
53 0.10004530E-02 0.96447733E+00
54 0.15495510E-02 0.96602688E+00
55 0.32832332E-01 0.99885921E+00
56 0.43282221E-03 0.99929203E+00
57 0.29315502E-03 0.99958519E+00
58 0.41481341E-03 0.10000000E+01
59 0.00000000E+00 0.00000000E+00
60 0.00000000E+00 0.00000000E+00
------------------------------------------------------------------------------
Final result:
Integral = 0.201003E+02 Standard dev. = 0.27E+00 with 19980 calls
------------------------------------------------------------------------------
Efficiency of the events acceptance = 0.410E-02
------------------------------------------------------------------------------
The narrow W-width approximation, e+ e- -\> b W+ b W- -\> 6f:
------------------------------------------------------------------------------
Results of the initial scan:
Kinematics Weight Integral St.dev. No of calls
1 0.47619048E-01 0.23618560E+00 0.10375292E+00 200
2 0.47619048E-01 0.25574899E+00 0.11393103E+00 200
3 0.47619048E-01 0.90318927E-02 0.27624478E-02 200
4 0.47619048E-01 0.14178256E+00 0.63376575E-01 200
5 0.47619048E-01 0.29188135E-02 0.15778093E-02 200
6 0.47619048E-01 0.75109770E-02 0.28446236E-02 200
7 0.47619048E-01 0.87010486E-01 0.36285359E-01 200
8 0.47619048E-01 0.20008583E-01 0.10505322E-01 200
9 0.47619048E-01 0.25432033E+00 0.12331680E+00 200
10 0.47619048E-01 0.14796519E+00 0.51432260E-01 200
11 0.47619048E-01 0.14443649E+00 0.47481858E-01 200
12 0.47619048E-01 0.13431004E+00 0.57817377E-01 200
13 0.47619048E-01 0.62963498E-01 0.18297859E-01 200
14 0.47619048E-01 0.14057286E+00 0.83384476E-01 200
15 0.47619048E-01 0.16790199E+02 0.67749184E+00 200
16 0.47619048E-01 0.61926860E-01 0.12581699E-01 200
17 0.47619048E-01 0.12722393E+00 0.37044922E-01 200
18 0.47619048E-01 0.15741819E+00 0.64646409E-01 200
19 0.47619048E-01 0.16399987E+00 0.54156119E-01 200
20 0.47619048E-01 0.21743058E+00 0.13237257E+00 200
21 0.47619048E-01 0.13765610E+00 0.84254619E-01 200
------------------------------------------------------------------------------
Born scanned with 4200 calls = 0.19300621E+02
The maximum value of the cross section= 0.98788719E+03
shifted to 0.14818308E+04
------------------------------------------------------------------------------
Weight 21 = 0.7132210911535855E-02 rescaled to 0.7132210911536152E-02
Calculating the Born cross section with the weights:
1 0.12237203E-01 0.12237203E-01
2 0.13250817E-01 0.25488019E-01
3 0.46795867E-03 0.25955978E-01
4 0.73460104E-02 0.33301989E-01
5 0.15122900E-03 0.33453218E-01
6 0.38915728E-03 0.33842375E-01
7 0.45081704E-02 0.38350545E-01
8 0.10366808E-02 0.39387226E-01
9 0.13176795E-01 0.52564021E-01
10 0.76663435E-02 0.60230365E-01
11 0.74835153E-02 0.67713880E-01
12 0.69588457E-02 0.74672726E-01
13 0.32622525E-02 0.77934978E-01
14 0.72833337E-02 0.85218312E-01
15 0.86993051E+00 0.95514882E+00
16 0.32085424E-02 0.95835737E+00
17 0.65917015E-02 0.96494907E+00
18 0.81561205E-02 0.97310519E+00
19 0.84971292E-02 0.98160232E+00
20 0.11265471E-01 0.99286779E+00
21 0.71322109E-02 0.10000000E+01
The maximum weight updated from 0.14818308E+04 to 0.16809386E+04
The maximum weight updated from 0.16809386E+04 to 0.17190927E+04
------------------------------------------------------------------------------
Final result:
Integral = 0.206457E+02 Standard dev. = 0.18E+00 with 20000 calls
------------------------------------------------------------------------------
Efficiency of the events acceptance = 0.129E-01
------------------------------------------------------------------------------
Approximated cross section e+ e- -\> t t -\> 6f:
------------------------------------------------------------------------------
Final result:
Integral = 0.193058E+02 Standard dev. = 0.66E-01 with 20000 calls
------------------------------------------------------------------------------
Born cross sections (fb):
------------------------------------------------------------------------------
ecm cs sd cstta sda csbbww sdbbww
0.5000E+03 0.20100E+02 0.27E+00 0.19306E+02 0.66E-01 0.20646E+02 0.18E+00


text:
# Dynamics of a two-component atomic Bose-Einstein condensate

## Abstract

The dynamical population oscillations between two internal states of a Bose-Einstein condensate are investigated within the rotating wave approximation. Analytical expressions for the population imbalance in the number states and coherent states have been derived, which predict different revival periods. Thus the true quantum state of the condensates may be unambiguously determined by detecting the atom intensity evolution for one internal state.
PACS: 03.75 Fi, 05.30 -d, 32.80.Pj, 74.20.D
Keywords: Bose-Einstein condensates, collapse and revival, Rabi oscillation.

The experimental observation of the Bose-Einstein condensation of a trapped, dilute gas of alkali atoms, and the high accuracy of the engineering, are opening a new avenue to investigate the interplay between macroscopics and quantum coherence. In particular, recent experiments on two-component Bose-Einstein condensates (BECs) in 87Rb have stimulated considerable theoretical work on the dynamics of the phase and number fluctuations of Bose condensates. Vigorous efforts have concentrated on the issue of temporal phase coherence between two coupled Bose condensates. Along semiclassical (mean) analysis much attention has focused on the coherent phenomena such as Josephson effect and macroscopic quantum self-trapping (MQST). Studies involving small quantum corrections show that due to the nonlinearities arising from the atom-atom interactions and the discrete spectrum of the many-body system, the dynamics of the population oscillations arising in the mean-field approximation is modulated by the collapses and revivals.

In the present work, we study quantum dynamics of a two-component atomic BEC such as 87Rb. Our purpose is to provide an analysis on how the coherent atom oscillations between two components depend on the initial states chosen for the condensates. We chose two different initial states of the condensate. In the first example we consider Fock state for each component of the condensate. In this case, we show that the revival period depends on whether the total atom number $`N`$ is even or odd. In the other example we invoke Bose-broken symmetry and consider the product of two coherent states. In this case we show that the revival period is always definite. In experiments, such as 87Rb, initial Fock states can be created by condensing independently two thermal atomic clouds trapped in different Zeeman level. Whereas, the initial coherent state may be created by splitting a condensate (trapped, say, in the Zeeman state $`|f=2,m_f=2`$) in two condensates trapped in the Zeeman states $`|f=2,m_f=2`$ and $`|f=2,m_f=1`$, respectively. Another aspect we are interested in this paper is the dependence of the coherent atom oscillations on the relative phase between two internal states and nonlinear atomic collisions.

We consider here Bose-Einstein condensation of a trapped gas of atoms that have two internal states $`|1`$ and $`|2`$ with energies $`\mathrm{}\omega _o/2`$ and $`\mathrm{}\omega _o/2`$, respectively. There is a spatially uniform radiation field with frequency $`\omega _e`$ that couples the two internal states with a Rabi frequency $`\mathrm{\Omega }_0`$. The atoms in states $`|1`$ and $`|2`$ are subject to isotropic harmonic trapping potentials $`V_i(𝐫)=\frac{1}{2}m\omega _ir^2`$ for $`i=\{1,2\}`$, respectively. Furthermore, the atoms interact via elastic two-body collisions through the interaction potentials $`V_{ij}(𝐫𝐫^{})=\frac{4\pi \mathrm{}^2a_{ij}}{m}\delta (𝐫𝐫^{})`$, Where $`a_{ij}`$ is the $`s`$-wave scattering length between atoms in states $`i`$ and $`j`$. It is assumed that $`a_{ij}>0`$ corresponding to repulsive interactions. In the formalism of second quantization, the system is described by the Hamiltonian

$$\widehat{H}=\widehat{H}_{atom}+\widehat{H}_{coll},$$
(1)
$`\widehat{H}_{coll}`$ $`=`$ $`{\displaystyle }d^3r\{{\displaystyle \frac{4\pi \mathrm{}^2a_{11}}{2m}}\widehat{\mathrm{\Psi }}_1^+(𝐫)\widehat{\mathrm{\Psi }}_1^+(𝐫)\widehat{\mathrm{\Psi }}_1(𝐫)\widehat{\mathrm{\Psi }}_1(𝐫)`$ ()
$`+{\displaystyle \frac{4\pi \mathrm{}^2a_{22}}{2m}}\widehat{\mathrm{\Psi }}_2^+(𝐫)\widehat{\mathrm{\Psi }}_2^+(𝐫)\widehat{\mathrm{\Psi }}_2(𝐫)\widehat{\mathrm{\Psi }}_2(𝐫)`$
$`+{\displaystyle \frac{4\pi \mathrm{}^2a_{12}}{m}}\widehat{\mathrm{\Psi }}_1^+(𝐫)\widehat{\mathrm{\Psi }}_2^+(𝐫)\widehat{\mathrm{\Psi }}_1(𝐫)\widehat{\mathrm{\Psi }}_2(𝐫)\},`$
$`\widehat{H}_{atom}`$ $`=`$ $`{\displaystyle }d^3r\{\widehat{\mathrm{\Psi }}_1^+(𝐫)[{\displaystyle \frac{\mathrm{}^2}{2m}}^2+V_1(𝐫)+{\displaystyle \frac{\mathrm{}\delta }{2}}]\widehat{\mathrm{\Psi }}_1(𝐫)`$ ()
$`+\widehat{\mathrm{\Psi }}_2^+(𝐫)\left[{\displaystyle \frac{\mathrm{}^2}{2m}}^2+V_2(𝐫){\displaystyle \frac{\mathrm{}\delta }{2}}\right]\widehat{\mathrm{\Psi }}_2(𝐫)`$
$`+{\displaystyle \frac{\mathrm{}\mathrm{\Omega }_0}{2}}(\widehat{\mathrm{\Psi }}_1^+(𝐫)\widehat{\mathrm{\Psi }}_2(𝐫)+\widehat{\mathrm{\Psi }}_2^+(𝐫)\widehat{\mathrm{\Psi }}_1(𝐫))\}.`$
Here, the atomic field operators $`\widehat{\mathrm{\Psi }}_i(𝐫)`$ and $`\widehat{\mathrm{\Psi }}_i^+(𝐫)`$ have been written in the field interaction representation which is rotating at the frequency of the external field $`\omega _e`$, $`\delta =\omega _0\omega _e`$ denotes the field detuning from resonance excitation.

Now we can approximate the field operators in the zero-temperature limit by a two-mode model such that $`\widehat{\mathrm{\Psi }}_i(𝐫)\widehat{c}_i\varphi _i(𝐫)`$, where $`\widehat{c}_i`$ is the annihilation operator which obeys the usual boson commutation relations. In the two-mode approximation the Hamiltonian becomes (with $`\mathrm{}=1`$),

$`\widehat{H}`$ $`=`$ $`({\displaystyle \frac{1}{2}}\delta +E_1)\widehat{c}_1^+\widehat{c}_1+({\displaystyle \frac{1}{2}}\delta +E_2)\widehat{c}_2^+\widehat{c}_2+{\displaystyle \frac{\mathrm{\Omega }}{2}}(\widehat{c}_1^+\widehat{c}_2+\widehat{c}_2^+\widehat{c}_1)`$ ()
$`+{\displaystyle \frac{\kappa _{11}}{2}}\widehat{c}_1^+{}_{}{}^{2}\widehat{c}_{1}^{}{}_{}{}^{2}+{\displaystyle \frac{\kappa _{22}}{2}}\widehat{c}_2^+{}_{}{}^{2}\widehat{c}_{2}^{}{}_{}{}^{2}+\kappa _{12}\widehat{c}_1^+\widehat{c}_1\widehat{c}_2^+\widehat{c}_2,`$
where

$$E_i=d^3r\varphi _i^{}(𝐫)[^2/(2m)+V_i(𝐫)]\varphi _i(𝐫),$$
(7)
$$\mathrm{\Omega }=\mathrm{\Omega }_0d^3r\varphi _2^{}(𝐫)\varphi _1(𝐫),$$
(8)
$$\kappa _{ij}=\frac{4\pi a_{ij}}{m}d^3r|\varphi _i(𝐫)|^2|\varphi _j(𝐫)|^2.$$
(9)
The analysis of Eq. (2) is greatly simplified by the introduction of the angular momentum operators

$$\widehat{J}_x=\frac{1}{2}(\widehat{c}_2^+\widehat{c}_2\widehat{c}_1^+\widehat{c}_1),$$
(10)
$$\widehat{J}_y=\frac{i}{2}(\widehat{c}_2^+\widehat{c}_1\widehat{c}_1^+\widehat{c}_2),$$
(11)
$$\widehat{J}_z=\frac{1}{2}(\widehat{c}_1^+\widehat{c}_2+\widehat{c}_2^+\widehat{c}_1).$$
(12)
The Casimir invariant is therefore

$$\widehat{J}^2=\frac{\widehat{N}}{2}(\frac{\widehat{N}}{2}+1),$$
(13)
where $`\widehat{N}=\widehat{c}_1^+\widehat{c}_1+\widehat{c}_2^+\widehat{c}_2`$ is the total number operator, which is a constant of motion and thus can be set equal to the total number of atoms $`N`$. Therefore the present two-mode approximation is transformed into an angular momentum model with total angular momentum given by $`j=N/2`$. In terms of angular momentum operators the Hamiltonian (2) may be rewritten as

$$\widehat{H}=\mathrm{\Omega }\widehat{J}_z+2\kappa \widehat{J}_x^2+f(j),$$
(14)
where

$$\kappa =\frac{1}{2}(\frac{\kappa _{11}}{2}+\frac{\kappa _{22}}{2}\kappa _{12}),$$
(15)
$$\kappa ^{}=\frac{1}{2}(\frac{\kappa _{11}}{2}+\frac{\kappa _{22}}{2}+\kappa _{12}),$$
(16)
$$f(j)=2\kappa ^{}j^2+(E_1+E_2\kappa \kappa ^{})j.$$
(17)
In deriving Eq. (6) we have assumed that

$$E_1+\delta E_2+(\frac{\kappa _{11}}{2}\frac{\kappa _{22}}{2})(2j1)=0,$$
(18)
a condition that can always be achieved by adjusting field detuning $`\delta `$ and shifting the energy levels of the condensate components. Note that some control of the nonlinear parameter $`\kappa `$ can be achieved through the proper engineering of the trapping potential or tuning the scattering lengths via Feshbach resonances.

In general the Hamiltonian (6) is difficult to treat in an exact way because of the presence of nonlinear atomic interactions. In order to gain physical insight into the dynamics of such a two-component BEC problem, some approximations are necessary: a common assumption is the so-called rotating wave approximation (RWA) $`\widehat{U}=e^{i\mathrm{\Omega }t\widehat{J}_z}`$. Using the operator identity

$$\mathrm{exp}(i\lambda \widehat{J}_z)\widehat{J}_x\mathrm{exp}(i\lambda \widehat{J}_z)=\widehat{J}_x\mathrm{cos}\lambda \widehat{J}_y\mathrm{sin}\lambda ,$$
(19)
the Hamiltonian in the rotating frame (RF) becomes

$`\widehat{H}_{RF}`$ $`=`$ $`\widehat{U}^+\widehat{H}\widehat{U}+i{\displaystyle \frac{d\widehat{U}^+}{dt}}\widehat{U}`$ ()
$`=`$ $`{\displaystyle \frac{1}{2}}\kappa (e^{2i\mathrm{\Omega }t}\widehat{J}_+^2+e^{2i\mathrm{\Omega }t}\widehat{J}_{}^2+\widehat{J}_+\widehat{J}_{}+\widehat{J}_{}\widehat{J}_+)+f(j),`$ (20)
where $`\widehat{J}_\pm =\widehat{J}_x\pm i\widehat{J}_y`$ are raising and lowing operators. Now we invoke RWA by suppressing rapidly oscillating time-dependent terms in Eq. (10). Thus we obtain the RWA Hamiltonian in the rotating frame as follows

$$\widehat{H}_{RF}^{(r)}=\kappa \widehat{J}_z^2+g(j),$$
(21)
where $`g(j)=\kappa j^2+f(j)`$. Instead of discussing dynamics in the rotating frame, we return to the Hamiltonian in the previous field interaction representation

$$\widehat{H}_{eff}=\widehat{U}\widehat{H}_{RF}^{(r)}\widehat{U}^+i\frac{d\widehat{U}^+}{dt}\widehat{U}=\kappa \widehat{J}_z^2+\mathrm{\Omega }\widehat{J}_z+g(j),$$
(21)
which is our starting Hamiltonian for the following discussions. It reveals in Eq. (12) that $`\widehat{J}_z`$ is now a constant of motion and solely determined by the initial condition. The frozen $`\widehat{J}_z`$ behavior can be understood when $`\mathrm{\Omega }`$ is sufficiently larger than $`\kappa `$, i.e., $`\mathrm{\Omega }\kappa `$. In this case, the external field forces the total spin to remain polarized in the $`z`$-direction because it costs energy to change the spin direction. To see the validity of our RWA Hamiltonian, we numerically integrate the Schödinger equation with the original Hamiltonian (6) and present in Fig. 1(a) the expectation values $`\widehat{J}_z`$ against time for various values of $`\mathrm{\Omega }`$ with the initial angular momentum state $`|j,j_x`$ where the subscript denotes the quantization axis along the $`x`$ direction. We remark that for the system starts from $`|j,j_x`$, the only nonvanishing spin component is $`\widehat{J}_z`$ because $``$ $`\widehat{J}_x=`$ $`\widehat{J}_y=0`$ at all times. Clearly it reveals in Fig. 1(a) that when $`\mathrm{\Omega }=0`$, $`\widehat{J}_z`$ vanishes after some time, which means that the initial coherence and polarization of $`\widehat{J}_z`$ is completely destroyed by the phase diffusion due to nonlinear term $`2\kappa \widehat{J}_x^2`$. However, when $`\mathrm{\Omega }`$ is sufficiently larger than $`\kappa `$, the value of $`\widehat{J}_z`$ remains almost unchanged during time evolution, consistent with Eq. (12) in which $`\widehat{J}_z`$ is a constant of motion. As a further illustration we calculate the energy spectrums of the two Hamiltonians and present the results in Fig. 1(b). It is found in calculations that the difference between the energies of two Hamiltonians approaches to be negligible when increasing the value of $`\mathrm{\Omega }`$. Thus we expect that the RWA Hamiltonian is valid in the parameter regime $`\mathrm{\Omega }j\kappa `$. In the following we will neglect the constant term $`g(j)`$ because this term denotes a simple constant energy shift and has no contribution to the dynamics of the system.

Now the Heisenberg equations of motion for the raising and lowing operators ($`\widehat{J}_\pm =\widehat{J}_x\pm i\widehat{J}_y`$ ) can be readily solved as follows

$$\widehat{J}_+(t)=e^{it(2\kappa \widehat{J}_z+\kappa +\mathrm{\Omega })}\widehat{J}_+,$$
(22)
$$\widehat{J}_{}(t)=e^{it(2\kappa \widehat{J}_z+\kappa \mathrm{\Omega })}\widehat{J}_{}.$$
(23)
Given the initial state $`\psi (0)`$ of the system, a quantity tailored to the dynamics is the population imbalance between the two components

$`N_{}(t)`$ $`=`$ $`2\psi (0)\widehat{J}_x(t)\psi (0)`$ ()
$`=`$ $`e^{i(\kappa +\mathrm{\Omega })t}\psi (0)e^{2i\kappa t\widehat{J}_z}\widehat{J}_+\psi (0)+c.c..`$ (24)
Number state description. We shall now proceed to demonstrate how the population imbalance evolves for two different kinds of initial states of the condensates. First we consider the case that the condensate is in the Fock (number) states, equally the angular momentum states $`j,m_x`$ where the value of $`2m`$ gives the atom number difference between the two BECs. We suppose that the condensates are initially localized in one well, corresponding to the state $`j,j_x`$. Substituting this initial state into Eq. (14) gives

$`N_{}(t)`$ $`=`$ $`e^{i(\kappa +\mathrm{\Omega })t}{\displaystyle \underset{m=j}{\overset{j}{}}}\sqrt{(j+m)(jm+1)}`$ ()
$`\times d_{j,m}^j({\displaystyle \frac{\pi }{2}})d_{m1,j}^j({\displaystyle \frac{\pi }{2}})e^{i2\kappa mt}+c.c.,`$
where the matrix elements $`d_{m,m^{}}^j(\theta )=j,m\mathrm{exp}(i\theta \widehat{J}_y)j,m^{}`$. Using the identity

$$d_{j,m}^j(\theta )=()^{jm}\left(\begin{array}{c}2j\\ j+m\end{array}\right)^{1/2}\left(\mathrm{cos}\frac{\theta }{2}\right)^{j+m}\left(\mathrm{sin}\frac{\theta }{2}\right)^{jm},$$
(26)
and recursive relation

$$d_{j,m}^j(\theta )=\frac{\mathrm{sin}(\theta /2)}{\mathrm{cos}(\theta /2)}\left(\frac{j+m+1}{jm}\right)^{1/2}d_{j,m+1}^j(\theta ),$$
(27)
we obtain

$`N_{}(t)`$ $`=`$ $`{\displaystyle \frac{1}{2^N}}e^{i(\kappa +\mathrm{\Omega })t}e^{i2j\kappa t}{\displaystyle \underset{m=j}{\overset{j}{}}}(j+m)\left(\begin{array}{c}2j\\ j+m\end{array}\right)e^{i2\kappa (j+m)t}+c.c.`$ (30)
$`=`$ $`N[\mathrm{cos}(\kappa t)]^{N1}\mathrm{cos}(\mathrm{\Omega }t).`$ ()
From Eq. (18) one can see that the population imbalance involves a rapidly oscillating part and an envelope function that is responsible for the collapse and revival of the population imbalance. The most prominent feature in Eq. (18) is that the parity of total atom number $`N`$, whether it is odd or even, plays an essential role in determining revival periods. When $`N`$ is odd, the revival period is $`\pi /\kappa `$, and the population imbalance is always larger than zero. That is, the population undergoes the macroscopic quantum self-trapping. However, when the number of atoms is even, $`N1`$ is odd and the condensate will undergo antirevivals when $`t=(2n+1)\pi /\kappa `$ so that $`N_{}(t)=N`$ at these times. In this case, the revival period becomes $`2\pi /\kappa `$ and the time average of $`N_{}(t)`$ is zero. Two examples are plotted in Figs. 2(a)-(b). To show the validity of the analytical results of Eq. (18) which is derived from the RWA Hamiltonian, we numerically calculate the Schrödinger equation with the original Hamiltonian (6). The results are shown in Fig. 2(c)-(d), corresponding to Fig. 2(a)-(c), respectively. Clearly, the RWA formula, equation (18), describes the system’s evolution very well when compared with the exact numerical results shown in Fig. 2(c)-(d), implying the different periods of quantum revivals for different parities of $`N`$. We notice that the reduction in revival amplitude shown in Fig. 2(c)-(d) is due to the state mixing induced by angular momentum operator $`\widehat{J}_x`$, which is eliminated in the RWA treatment. If the Fock states with a fixed number of atoms are appropriate for the description of the condensates, then by monitoring the initially unoccupied internal state using off-resonant light scattering, this change for the magnitude of the revival period of the population imbalance between the two components may be observed. A similar phenomenon has been discussed in Ref. in the context of the interference pattern, without considering coupling and nonlinear atomic collision between the two components. However, because the visibility of the interference pattern is never negative, so the revival period considered in Ref. is always $`\pi /\kappa `$ with a $`\pi `$ phase shift for different parities of $`N`$.

Coherent state description. We turn now to the coherent state description of the problem. It is customary to consider the condensates to be in two coherent states, associated with a macroscopic wave function with both magnitude and phase, the presence of which is due to the spontaneous breaking of gauge symmetry. This is the view we will adopt here, to assume the condensate to be such coherent states with definite relative phase

$$\psi (0)=\alpha ,\beta .$$
(31)
The amplitude of the two components may be presented by $`\left|\alpha \right|=\sqrt{N}\mathrm{sin}(\theta )`$ and $`\left|\beta \right|=\sqrt{N}\mathrm{cos}(\theta )`$, conserving total atom number. The relative phase between the two components is defined to be $`\varphi `$. Thus the two complex amplitudes are related by $`\beta \alpha ^{}=\frac{1}{2}N\mathrm{sin}(2\theta )\mathrm{exp}(i\varphi )`$. To give a closed expression for Eq. (14) with the initial state Eq. (19), we need the normal-ordering arrangement for the operator function $`\mathrm{exp}(\lambda \widehat{J}_z)\widehat{J}_+`$. With the aid of the following identity

$$e^{\lambda \widehat{c}_i^+\widehat{c}_i}=\underset{r=0}{\overset{\mathrm{}}{}}\frac{x^r}{r!}\widehat{c}_i^{+r}\widehat{c}_i^r,$$
(32)
where $`x=e^\lambda 1`$, and after a straightforward calculation, we obtain

$`e^{2\lambda \widehat{J}_z}\widehat{J}_+`$ $`=`$ $`{\displaystyle \frac{1}{2}}e^\lambda (\widehat{c}_2^++\widehat{c}_1^+){\displaystyle \underset{m,n,p,q=\mathrm{}}{\overset{\mathrm{}}{}}}\{{\displaystyle \frac{x^mx^ny^qz^p}{m!n!p!q!}}`$ ()
$`\times \widehat{c}_1^+{}_{}{}^{m}\widehat{c}_{1}^{+n}\widehat{c}_2^{+p}\widehat{c}_2^{+q}[{\displaystyle \underset{l=\mathrm{}}{}}{\displaystyle \frac{z^l}{l!}}(\widehat{c}_2^++\widehat{c}_1^+)^l(\widehat{c}_2+\widehat{c}_1)^l]`$
$`\times \widehat{c}_1^m\widehat{c}_1^p\widehat{c}_2{}_{}{}^{n}\widehat{c}_{2}^{q}\}(\widehat{c}_2\widehat{c}_1),`$
where $`y=(e^{2\lambda }+1)/2`$, and $`z=(e^{2\lambda }1)/2`$. Substituting Eqs. (20)-(21) into Eq. (14), we have the final result for the population imbalance at time $`t`$ in the coherent state representation

$$N_{}(t)=A\mathrm{exp}(2N\mathrm{sin}^2\frac{\kappa t}{2})\mathrm{cos}[W\mathrm{sin}(\kappa t)+\mathrm{\Omega }t+\phi ],$$
(35)
where we have defined $`A=N\sqrt{1\mathrm{cos}^2\varphi \mathrm{sin}^2(2\theta )}`$, $`W=N\mathrm{sin}(2\theta )\mathrm{cos}\varphi `$, and a constant phase shift $`\phi =\mathrm{arctan}[\mathrm{sin}\varphi \mathrm{tan}(2\theta )]`$. We can see from Eq. (22) that the dynamics of the population imbalance is mainly determined by the envelope function $`\mathrm{exp}(2N\mathrm{sin}^2\frac{\kappa t}{2})`$, with modulated by a rapidly varying cosine term. The characteristics of collapses and revivals can be clearly seen from Eq. (22). For comparison with Eq. (18) derived in the number state space, let us chose the value of $`\theta =0`$, corresponding the case that the condensate is initially localized in one internal state. Then $`N_{}(t)`$ in Eq. (22) reduces to

$$N_{}(t)=N\mathrm{exp}(2N\mathrm{sin}^2\frac{\kappa t}{2})\mathrm{cos}(\mathrm{\Omega }t).$$
(36)
Although no relative phase information occurs in Eq. (23), quite different from number states, revivals for the condensate in the coherent states occur at definite times $`t=2\pi n/\kappa `$ ($`n`$ is an integer ), independent of the parities of the average particle number $`N`$. As illustrated in Ref., this difference of revival periods between the Fock state and coherent state descriptions originates from the fact that the atom number difference operator $`\widehat{c}_2^+\widehat{c}_2\widehat{c}_1^+\widehat{c}_1`$ is quantized in units of $`2`$ when the condensates are in Fock states and units of $`1`$ for coherent states.

We emphasize here that for any realistic experimental set up, the exact atom number is unknown and maybe has a large number fluctuation over many different experimental runs. As a result, what we obtain in Fock states should be wieghted by a factor $`p(N)`$ describing the probability the state $`|j,j_x`$ occurs at the begaining. According to quantum measurement theory, however, the atom number in any given experiment would be well defined. Thus in the Fock-state description, the occurence of antirivaval phenomenon, which feaures the parity of the atom number, is robust for any single experiment. On the other hand, the antirevival phenomenon never occurs in the coherent-state description, which can be seen from Eq. (23). This is because the coherent state is a superposition of the Fock states with the weighting factor satisfying the Poisson number distrbution. Consequently, in each single experiment quantum fluctuations will be present in the coherent state description, which smears out the antirevivals. This fact allows us the possibility to detect the true quantum state of the two-component BEC system in realistic experiment.

In the case of the values of nonlinear interactions satisfying the following condition

$$\kappa =\frac{1}{2}(\frac{\kappa _{11}}{2}+\frac{\kappa _{22}}{2}\kappa _{12})0,$$
(37)
our RWA treatment of Hamiltonian (6) is exact and Eq. (22) reduces to description of a simple Rabi oscillation with period $`2\pi /\mathrm{\Omega }`$. This form of atomic collective excitation has been observed in Ref.. The existence of Rabi oscillation also depends on the initial state conditions, which is reflected by the fact that the oscillation amplitude $`A`$ in Eq. (22) is determined by the initial population distribution and relative phase between the two components. In particular when the value of $`\theta =\pi /4`$, corresponding to equal atom distribution between the two internal states, and $`\varphi =n\pi `$, we arrive at the conclusion that the oscillation is completely suppressed and no quantum tunneling occurs during time evolution. Note that the present conditions of zero oscillation is identical to the fixed points discussed in Ref., in which the semiclassical approximation has been adopted. In general the value of effective interaction $`\kappa `$ is not zero and it can be adjusted through the proper engineering of the trapping potential, and hence of the condensate wave functions $`\varphi _i(𝐫)`$. To illustrate the interplay between the weak nonlinearity and strong Josephson-like coupling, we present in Fig. (3) $`N_{}(t)`$ \[Eq. (23)\] for several values of $`\kappa `$. It reveals in Fig. (3) that the presence of the effective nonlinear interaction modulates the fringe visibility of the internal Rabi oscillations as a manner of collapses and revivals, as has been verified by JILA group. When $`\kappa `$ increases to large values, the rapid phase diffusion caused by the nonlinear interaction destroys the collective excitation of atoms, corresponding to the complete smearing of the Rabi oscillations, as illustrated in Fig. 3(d).

To see the effect of the initial relative phase on the dynamics of the system, we present in Fig. 4 the time evolution of $`N_{}(t)`$ for four different values of $`\varphi `$. The value of $`\theta `$ in these figures is chosen to be $`\pi /4`$, corresponding to the two condensates with equal atoms. We can see in Fig. 4 that the oscillation amplitude is maximal for the values $`\varphi =\pi /2`$. Experimental observation of this change of oscillation amplitudes may help us determine the relative phase between the two condensates.

We have restricted our analysis to zero temperature. Our study of the effect of initial state selection on the consequent dynamics of a two-component condensate should be extended to finite temperatures, where dissipative effects associated with the thermal cloud of noncondensate atoms must be included. Further study of damping in this system is needed.

In summary, with the help of rotating wave approximation, we have derived analytical expressions for the population imbalance between the two Bose condensates in the number states and coherent states. The time evolution of the coherent atom oscillations is shown to be quite different between these two state descriptions. For the condensates in number states, the revival periods of the oscillations are either $`\pi /\kappa `$ or $`2\pi /\kappa `$, depending on whether the total atom number is odd or even. Whereas, for the condensates in coherent states, the revival periods are always $`2\pi /\kappa `$. So the true quantum state of the condensates may be unambiguously determined by detecting the atom intensity evolution for one trap.

Figure captions

Fig. 1. (a) Time evolution of $`\widehat{J}_z`$ for the value of total atom number $`N=100`$. The initial state is $`|j,j_x`$; (b) energy spectrums of the exact Hamiltonian (closed triangles) and the RWA Hamiltonian (closed squares) for the values of $`N=100,\mathrm{\Omega }=30\kappa `$.

Fig. 2. Analytical results of population imbalance $`N_{}(t)`$ for the value of total atom number (a) $`N=100`$ and (b) $`N=101`$. The exact numerical results corresponding to (a) and (b) are shown in (c) and (d), respectively. Other parameters are $`\kappa =0.001\mathrm{\Omega }`$.

Fig. 3. Population imbalance $`N_{}(t)`$ for the value of effective nonlinear coupling (a) $`\kappa =0`$, (b) $`\kappa =0.001\mathrm{\Omega }`$, (c) $`\kappa =0.005\mathrm{\Omega }`$, and (d) $`\kappa =0.025\mathrm{\Omega }`$. Other parameters are $`N=100`$ and $`\theta =0`$.

Fig. 4. Population imbalance $`N_{}(t)`$ for the value of relative phase (a) $`\varphi =0`$, (b) $`\varphi =\pi /6`$, (c) $`\varphi =\pi /4`$, and (d) $`\varphi =\pi /2`$. Other parameters are $`N=100`$ and $`\kappa =0.01\mathrm{\Omega }`$.


text:
# 1 Introduction

## 1 Introduction

The understanding of the origin of fermion families and the observed pattern of fermion masses and mixings, together with the origin of CP violation, are among the major outstanding problems in particle physics.

Most extensions of the standard model (SM) naturally include new sources of CP violation. In supersymmetric (SUSY) extensions of the SM we find additional sources of CP violation, due to the presence of new CP violating phases. However, these new phases are severely constrained to be small by experimental bounds on the electric dipole moments (EDM’s) .

Since the question of CP-violation is closely related to the general flavour problem, one may wonder whether it is possible, within a supersymmetric extension of the SM, to establish a connection between the need for small CP violating phases and the observed pattern of quark masses and mixings. Such a connection might be possible if one imposes universality of strength for Yukawa couplings (USY) on a supersymmetric extension of the SM. In Ref. , we introduce the USY ansatz in the framework of SUSY, analysing the compatibility of this scenario with the EDM bounds, and discussing the contributions to $`K`$ and $`B`$ system CP violation observables.

## 2 Universal strength of Yukawa couplings

One of the most simple, and yet very attractive suggestions for the structure of the Yukawa couplings consists in assuming that they have a universal strength $`g_y`$, so that the mass matrices can be written as $`m_{ij}=g_yv/\sqrt{2}\mathrm{exp}(i\varphi _{ij})`$, where $`i,j`$ denote family indices. In the limit of small phases, the USY matrices can be viewed as a perturbation of the democratic type matrices . Within the framework of universal strength of Yukawa couplings, the quark Yukawa matrices can be parametrized as $`U_{ij}=\lambda _u/3\mathrm{exp}[i\mathrm{\Phi }_{ij}^u]`$ and $`D_{ij}=\lambda _d/3\mathrm{exp}[i\mathrm{\Phi }_{ij}^d]`$, where $`\lambda _{u,d}`$ are overall real constants, and $`\mathrm{\Phi }^{u,d}`$ are pure phase matrices. By performing appropriate weak-basis transformations the matrices $`U`$ and $`D`$ can, without loss of generality, be put in the form

$$U=\frac{\lambda _u}{3}\left(\begin{array}{ccc}e^{ip_u}& e^{ir_u}& 1\\ e^{iq_u}& 1& e^{it_u}\\ 1& 1& 1\end{array}\right);D=\frac{\lambda _d}{3}K^{}\left(\begin{array}{ccc}e^{ip_d}& e^{ir_d}& 1\\ e^{iq_d}& 1& e^{it_d}\\ 1& 1& 1\end{array}\right)K,$$
(1)
where $`K=\mathrm{diag}(1,e^{i\kappa _1},e^{i\kappa _2})`$. The phases $`(p,q,r,t)`$ affect both the spectrum and the $`V_{\mathrm{CKM}}`$, while $`\kappa _i`$ only influence the $`V_{\mathrm{CKM}}`$. The quark mass hierarchy forces the modulus of each $`(p,q,r,t)`$ phase to be small, at most of order $`m_2/m_3`$. Also, the fact that $`(|V_{cb}|^2+|V_{ub}|^2)`$ is small, constrains each one of the phases $`\kappa _1,\kappa _2`$ to be at most of order $`|V_{cb}|`$ .

In Ref. , it was shown that USY can easily account for the quark mass spectrum at low energy scales. This class of Yukawa couplings can be generated at some high energy scale through the breaking of a flavour symmetry, as discussed in , and it has been also recently argued that USY textures may naturally arise within the framework of a theory with two “large extra dimensions” . Furthermore, RG evolution preserves the USY texture that is responsible for the non-degenerate quark spectrum .

By scanning the parameter space, one can find values of $`(p,q,r,t,\kappa _1,\kappa _2)`$, such that the quark spectrum and the $`V_{\mathrm{CKM}}`$ can be correctly reproduced. To illustrate some of the features of this ansatz, let us consider, as an example, the following choice of USY phases:

The overall factors have been defined as $`\lambda _u/3=m_t/v\mathrm{sin}\beta `$ and $`\lambda _d/3=m_b/v\mathrm{cos}\beta `$. The Yukawa couplings can be diagonalized as $`S_R^{d(u)}Y^{d(u)^T}S_L^{d(u)^{}}=Y_{\mathrm{diag}}^{d(u)}`$, where $`S_R`$ and $`S_L`$ are unitary matrices. It is worth stressing that within the context of the SM, the matrices $`S_{L,R}^{u,d}`$ do not have, by themselves, any physical meaning, only the combination $`S_L^u.S_L^d^{}V_{\mathrm{CKM}}`$ is physically meaningful. However the matrices $`S_{L,R}^{u,d}`$ do play a significant rôle in some extensions of the SM, as for example the MSSM with non-universal soft-breaking terms. In the present USY model the $`V_{\mathrm{CKM}}`$ is given by

$$V_{\mathrm{CKM}}=\left(\begin{array}{ccc}0.9756.88\times 10^4i& 0.221+3.96\times 10^3i& 0.00274\times 10^5i\\ 0.2209.71\times 10^3i& 0.973+0.061i& 0.041+1.92\times 10^3i\\ 0.0063+1.10\times 10^4i& 0.0411.75\times 10^3i& 0.999+0.027i\end{array}\right).$$
(2)
Regarding CP violation, the above $`V_{\mathrm{CKM}}`$ leads to a strength of CP breaking (measured by $`J\mathrm{Im}(V_{ud}V_{cs}V_{us}^{}V_{cd}^{})`$) that is too small to account for the observed value of $`\epsilon _K`$. This is a generic feature of the USY ansatz, and it provides further motivation to embed USY in a larger framework, where new sources of CP violation naturally arise, like the MSSM.

## 3 Supersymmetric USY models and CP violation

In Ref. , we considered the minimal supersymmetric standard model (MSSM), where a minimal number of superfields is introduced and $`R`$ parity is conserved, with the following soft SUSY breaking terms

$`V_{SB}`$ $`=`$ $`m_{0\alpha }^2\varphi _\alpha ^{}\varphi _\alpha +ϵ_{ab}(A_{ij}^uY_{ij}^uH_2^b\stackrel{~}{q}_{L_i}^a\stackrel{~}{u}_{R_j}^{}+A_{ij}^dY_{ij}^dH_1^a\stackrel{~}{q}_{L_i}^b\stackrel{~}{d}_{R_j}^{}+A_{ij}^lY_{ij}^lH_1^a\stackrel{~}{l}_{L_i}^b\stackrel{~}{e}_{R_j}^{}`$ (3)
$``$ $`B\mu H_1^aH_2^b+\mathrm{H}.\mathrm{c}.){\displaystyle \frac{1}{2}}(m_3\overline{\stackrel{~}{g}}\stackrel{~}{g}+m_2\overline{\stackrel{~}{W^a}}\stackrel{~}{W}^a+m_1\overline{\stackrel{~}{B}}\stackrel{~}{B}),`$
where the Yukawa couplings are of the USY form, and $`\varphi _\alpha `$ denotes all the scalar fields of the theory. The $`\mu `$ term and the gaugino soft terms are assumed to be real, so that the structure and phases of the trilinear terms will play a key rôle, regarding both the stringent bounds coming from the EDM’s, and the enhancement of the CP observables. In this work we will also take as a guideline the assumption that all supersymmetric phases should be no greater than the largest of the USY phases in Table 1.

We begin by investigating how the EDM bounds constrain the parameter space for SUSY models with universal strength of Yukawa couplings. The current experimental bound on the EDM of the neutron and mercury atom are $`d_n<6.3\times 10^{26}e\text{cm}`$ and $`d_{Hg}<2.1\times 10^{28}e\text{cm}`$, respectively . These bounds can be translated into constraints for the imaginary parts of the flavour conserving $`LR`$ mass insertions. Requiring that $`d_n`$ does not exceed the experimental limit compels $`\mathrm{Im}(\delta _{11}^{d(u)})_{LR}`$ to be less than $`10^6`$, while compatibility with the mercury EDM corresponds to having $`\mathrm{Im}(\delta _{11}^{d(u)})_{LR}<10^710^8`$ and $`\mathrm{Im}(\delta _{22}^d)_{LR}<10^510^6`$ .

We showed that, even in the limit of small (or vanishing) supersymmetric phases, the large mixing inherent to USY couplings (displayed in $`S_L`$ and $`S_R`$), together with the EDM experimental limits, severely constrain any non-universality of the trilinear terms. In Fig. 1 we present the constraints from the EDM’s on the off-diagonal entries of the $`A`$-terms, in particular $`A_{12}`$ and $`A_{13}`$. In our analysis we assumed $`\mathrm{tan}\beta =5`$, $`m_0=m_{1/2}=250\mathrm{GeV}`$ and $`A_{ij}=m_0`$ for all elements except $`A_{12,13}`$, which we set in the range $`[3m_0,3m_0]`$.

From these figures, it is clear that any significant non-universality among the $`A`$ terms leads to unacceptably large contributions to the EDM’s. Similar constraints hold for the other off-diagonal elements $`A_{21}`$, $`A_{31}`$, $`A_{23}`$ and $`A_{32}`$, and this situation is far more severe than the case of hierarchical Yukawa couplings. Since Hermitian USY couplings lead to an unrealistic sum rule for the quark masses, the associated EDM problem cannot be overcome by taking simultaneously Hermitian Yukawa and trilinear couplings . An interesting possibility is to have trilinear terms that can be factorized as $`\widehat{A}=A_{ij}Y_{ij}=A.Y`$ or $`Y.A`$ . This factorization implies that the mass insertion $`(\delta _{11}^d)_{LR}`$ is suppressed by the ratio $`m_d/m_{\stackrel{~}{q}}`$. To be specific, let us consider the following structure

$$A=m_0\left(\begin{array}{ccc}a& a& a\\ b& b& b\\ c& c& c\end{array}\right).$$
(4)
In this case the trilinear couplings $`\widehat{A}`$ can be written as $`\widehat{A}=\mathrm{diag}(a,b,c).Y`$, with complex $`(a,b,c)`$. For $`|A_{ij}|3m_0`$, the EDM’s constraint on these phases is quite stringent:

$$\phi _a<0.02\mathrm{rad};\phi _b<0.35\mathrm{rad};\phi _c<0.02\mathrm{rad}.$$
(5)
Thus, we can conclude that in USY models the maximal allowed non-universality for the $`A`$-terms is the structure presented in Eq. (4), with the associated SUSY phases ($`\phi _{ij}`$) within the limits given in Eq. (5). In view of this, we analyse the contributions to the kaon system CP violating observables, $`\epsilon _K`$ and $`\epsilon ^{}/\epsilon `$.

Let us start our analysis by considering the indirect CP violating parameter of the $`K`$ sector, $`\epsilon _K`$. In the presence of supersymmetric ($`\stackrel{~}{g}`$ and $`\stackrel{~}{\chi }^\pm `$) contributions, the off-diagonal entry in the kaon mass matrix can be decomposed as $`M_{12}=M_{12}^{\mathrm{SM}}+M_{12}^{\stackrel{~}{g}}+M_{12}^{\stackrel{~}{\chi }^\pm }`$. Within the present scenario, with USY phases as in Table 1, we find that the SM contribution to $`\epsilon _K`$ is $`𝒪(10^5)`$. Regarding the supersymmetric contributions, and taking the $`A`$-terms as in Eq. (4), we have studied the behaviour of the $`(\delta _{12}^d)_{LL}`$, $`(\delta _{12}^d)_{LR}`$ and $`(\delta _{12}^u)_{LL}`$, mass insertions, focusing on the correlation between the $`(\delta _{12}^d)_{LL}`$ and $`(\delta _{12}^d)_{LR}`$. We have verified that gluino contributions clearly dominate over those of the charginos, and that the $`LL`$ mass insertions provide the dominant gluino contribution. We notice that in our model, the large mixing displayed by the rotation matrices enhances any non-diagonal contribution to $`M_{\stackrel{~}{Q}_L}^2`$ induced by the trilinear terms from RG evolution, hence $`LL`$ mass insertions can account for the experimental value of $`\epsilon _K`$.

We then consider supersymmetric contributions to the direct CP violation parameter, $`\epsilon ^{}/\epsilon `$, which can be symbolically decomposed as a sum of the SM, gluino, and chargino terms, the latter receiving contributions from diagrams with one and two mass insertions in the squark internal line.

In USY scenarios, the SM prediction for $`\epsilon ^{}/\epsilon `$ is found to be $`𝒪(10^6)`$. The supersymmetric contributions were analysed with detail in . In order to discuss the rôle of each contribution within our model, in Fig. 2(a) we plot the gluino and chargino contributions to $`\epsilon ^{}/\epsilon `$ as function of $`\delta A`$, which is defined as

$$\delta A=(A_{3i}A_{1i,2i})/m_0=ca(\mathrm{assuming}a=b).$$
(6)
The SUSY CP violating phases have been fixed as: $`\phi _a=\phi _c=0`$ and $`\phi _b=0.1`$. As in the previous figures, we have assumed $`\mathrm{tan}\beta =5`$ and $`m_0=m_{1/2}=250`$ GeV.

As can be seen from this figure, the chargino contributions with one mass insertion give the dominant contribution to $`\epsilon ^{}/\epsilon `$. The chargino contributions with two mass insertions that arise from the SUSY effective $`\overline{s}dZ`$ vertex are also relevant and considerably larger than that of the gluino.

Finally, we address the issue of CP violation in the $`B_d`$ meson system. The present world average, $`a_{J/\psi K_S}=0.79\pm 0.12`$, is dominated by the BaBar and Belle results . In the case of supersymmetric contributions to the $`\mathrm{\Delta }B=2`$ transition, the ratio of the total and partial SM contributions can be parametrized as $`r_d^2e^{2i\theta _d}=M_{12}(B_d)/M_{12}^{\mathrm{SM}}(B_d)`$, hence the measurement of $`a_{J/\psi K_S}`$ does not determine $`\mathrm{sin}2\beta `$ but rather $`\mathrm{sin}(2\beta +2\theta _d)`$. The SM contributions to $`\mathrm{sin}2\beta `$ associated with the phases of Table 1 are negligible, $`𝒪(10^3)`$.

Regarding SUSY contributions, we took into account charged-Higgs, chargino and gluino contributions to $`M_{12}^{\mathrm{SUSY}}(B_d)`$. Since in this class of models flavour mixing is not suppressed, one finds that the $`(\delta _{13}^{d,u})_{LL}`$ and $`(\delta _{13}^{d,u})_{LR}`$ mass insertions can be significantly large and saturate the experimental result. We have found that within our scenario, the leading supersymmetric contribution to $`\mathrm{sin}(2\beta +2\theta _d)`$ was associated with the $`\mathrm{\Delta }B=2`$ chargino mediated box diagrams. In Fig. 2(b) we present the correlation between the values of $`\epsilon _K`$ and $`\mathrm{sin}2\theta _d`$ for $`\mathrm{tan}\beta =5`$ and $`m_0=m_{1/2}=250`$ GeV. The values of the parameters $`a,b,`$ and $`c`$ are randomly selected in the range $`[3,3]`$ and the phases fixed as $`\phi _a=\phi _c=0`$ and $`\phi _b=0.1`$. From Fig. 2(b), it is clear that one can have $`\mathrm{sin}2\theta _d`$ within the experimental range while having a prediction for $`\epsilon _K`$ and $`\epsilon ^{}/\epsilon `$ compatible with the measured value. This result appears as characteristic of the kind of models under consideration. Despite the smallness of the phases introduced in these models via the Yukawa and trilinear couplings, they are sufficient to account for CP violation in the $`K`$-system, as well as the observed value of $`a_{J/\psi K_S}`$.

## 4 Discussion and Conclusions

In this work we have studied the implications of having universal strength of Yukawa couplings within the unconstrained MSSM.

We have argued that the trilinear soft terms play a key rôle in embedding USY into SUSY. In fact, due to the large mixing and associated phases, the constraints from the EDM’s on the SUSY parameter space are far more stringent than in the case of a standard Yukawa parametrization. We found that in order to satisfy the bound of the mercury EDM, the $`A`$-terms should be matrix factorizable, with phases constrained to be in the range $`10^210^1`$.

We have investigated the new contributions to both $`K`$ and $`B`$ system CP observables, finding that gluino mediated boxes with $`LL`$ mass insertions provide the leading contributions to $`\epsilon _K`$, while $`\epsilon ^{}/\epsilon `$ is dominated by chargino loops, through $`LL`$ flavour mixing. Regarding the $`B`$ system, we argued that within this model supersymmetric chargino exchanges, provide the leading contributions, which are in agreement with the recent measurements at BaBar and Belle.

In conclusion, we have presented an alternative scenario for CP violation, where the strength of CP violation stemming from the SM is naturally small, so that new SUSY contributions are essential to generate the correct value of $`\epsilon _K`$ and $`\epsilon ^{}/\epsilon `$, as well as the recently observed large value of $`a_{J/\psi K_S}`$.

Acknowledgments

This work was supported in part by the Portuguese Ministry of Science through project CERN/P/Fis/40134/2000, CERN/P/Fis/43793/2001, and by the E.E. through project HPRN-CT-2000-001499. M.G. and A.T. acknowledge support from ’Fundacão para a Ciência e Tecnologia’, under grants SFRH/BPD/5711/2001 and PRAXIS XXI BD/11030/97, respectively. The work of S.K. was supported by PPARC.