text:
# 𝐸1-𝐸2 interference in the Coulomb dissociation of 8B

## I INTRODUCTION

The Coulomb dissociation (CD) method provides an alternate indirect way to determine the cross sections for the radiative capture reaction 7Be (p, $`\gamma `$) 8B at low relative energies , which is a key reaction of the $`pp`$ chain through which energy is generated in the sun. The rate of this reaction is the most uncertain nuclear input to the standard solar model calculations , which affects the high energy solar neutrino flux and bears significantly on the solar neutrino problem. The CD method reverses the radiative capture by the dissociation of a projectile (the fused system) in the Coulomb field of a target, by making the assumption that nuclei do not interact strongly and the electromagnetic excitation process is dominated by a single multipolarity . Therefore, one has to pin down a kinematical domain where breakup cross sections due to different multipolarities ($`E1`$ and $`E2`$) can be clearly separated from each other. At the same time, it is also necessary to investigate the importance of the $`E1E2`$ interference effects which has not been considered so far in most of the analyses of the Coulomb dissociation data of 8B.

The magnitudes of the $`E1E2`$ interference cross sections could be appreciable for certain observables in the Coulomb dissociation of 8B. It was shown in Refs.  that the dissociation of 8B around 45 MeV/nucleon beam energy is dominated by $`E1`$ transitions but the interference with $`E2`$ amplitudes produces large asymmetries in the angular and momentum distributions of the breakup fragments. This fact has been used in a recent measurement of the parallel momentum distribution of the 7Be fragment resulting from the breakup of 8B on a Pb target at 44 MeV/nucleon beam energy to put constraint on the contributions of the $`E2`$ component in the breakup of 8B .

In this paper, our aim is to investigate the role of the $`E1E2`$ interference effects in the analysis of the Coulomb dissociation data at beam energies 25.8 MeV at Notre Dame  and around 50 MeV/nucleon measured at RIKEN . We would like to examine the importance of the interference effects in different cases and try to see whether they could be really absent in the above data. To this end, we have developed a first order semiclassical theory for the electromagnetic excitation of a composite nucleus leading to continuum final states which allows to estimate the $`E1`$ and $`E2`$ breakup contributions separately as well as their interference. Our formulation is different from that of Refs.  in the sense that we do not assume straight line trajectories for the motion of the projectile in the Coulomb field of the target nucleus. Our theory is closer to that of Baur and Weber  where triple differential cross sections with respect to the relative and centre of mass (c.m.) angles and energies of the fragments are determined. However, we have used a more general coupling scheme of angular momenta and spins than those in Ref. .

We organize the paper as follows. The theoretical formalism has been described in section 2. We describe the structure model adopted for 8B in section 3. The results and discussions on them are presented in section 4. Finally, the summary and conclusions are given in section 5.

## II FORMALISM

We consider the following reaction

$`a+Ac+v+A,`$ (1)
in which the composite projectile $`a(=c+v)`$ breaks up into a core $`c`$ and a valence particle $`v`$ in the Coulomb field of the target nucleus $`A`$ which remains in the ground state (elastic breakup). We consider electric transitions of dipole and quadrupole types only. It is to be noted that the $`M1`$ transition contribution is small for the reactions considered here at beam energies $`<`$52 MeV/nucleon .

In the initial channel, for given spins $`s_{c_i}`$ and $`s_{v_i}`$ of the core and valence particles, the projectile wave function is

$`\psi _{I_iM_i}(\stackrel{}{r})`$ $`=`$ $`{\displaystyle \underset{l_iL_i}{}}\left[[Y_{l_i}(\widehat{\stackrel{}{r}})\chi _{s_{v_i}}]_{L_i}\chi _{s_{c_i}}\right]_{I_iM_i}`$ (2)
$`\times `$ $`{\displaystyle \frac{f_{(L_is_{c_i})I_i}(r)}{r}}.`$ (3)
In the above, $`l_i`$ is the orbital angular momentum of the valence particle relative to the core and $`I_i`$ is the total spin of the projectile in its ground state with projection $`M_i`$. $`\chi `$’s are the spin wave functions of the core and the valence particle.

More explicitly

$`\psi _{I_iM_i}(\stackrel{}{r})`$ $`=`$ $`[{\displaystyle \underset{\genfrac{}{}{0pt}{}{\genfrac{}{}{0pt}{}{l_im_{l_i}}{L_iM_{L_i}}}{m_{v_i}m_{c_i}}}{}}(1)^{l_i+s_{v_i}M_{L_i}}\sqrt{2L_i+1}`$ (4)
$`\times `$ $`\left(\begin{array}{ccc}l_i& s_{v_i}& L_i\\ m_{l_i}& m_{v_i}& M_{L_i}\end{array}\right)Y_{l_im_{l_i}}(\widehat{\stackrel{}{r}})\chi _{s_{v_i}m_{v_i}}]`$ (7)
$`\times `$ $`(1)^{L_i+s_{c_i}M_i}\sqrt{2I_i+1}\left(\begin{array}{ccc}L_i& s_{c_i}& I_i\\ M_{L_i}& m_{c_i}& M_i\end{array}\right)`$ (10)
$`\times `$ $`\chi _{s_{c_i}m_{c_i}}{\displaystyle \frac{1}{r}}f_{(L_is_{c_i})I_i}(r).`$ (11)
In the above, $`m_{l_i}`$ is the projection of $`l_i`$.

For the final channel, the wave function is

$`\psi _{\stackrel{}{k}s_{v_f}m_{v_f}s_{c_f}m_{c_f}}(\stackrel{}{r})={\displaystyle \underset{\genfrac{}{}{0pt}{}{\genfrac{}{}{0pt}{}{I_fM_fm_{v_f}^{}m_{c_f}^{}}{l_fm_{l_f}m_{l_f}^{}}}{L_fM_{L_f}M_{L_f}^{}}}{}}[[(1)^{l_f+s_{v_f}M_{L_f}^{}}`$ (12)
$`\times `$ $`\sqrt{2L_f+1}\left(\begin{array}{ccc}l_f& s_{v_f}& L_f\\ m_{l_f}^{}& m_{v_f}^{}& M_{L_f}^{}\end{array}\right)Y_{l_fm_{l_f}^{}}(\widehat{\stackrel{}{r}})\chi _{s_{v_f}m_{v_f}^{}}]`$ (15)
$`\times `$ $`(1)^{L_f+s_{c_f}M_f}\sqrt{2I_f+1}\left(\begin{array}{ccc}L_f& s_{c_f}& I_f\\ M_{L_f}^{}& m_{c_f}^{}& M_f\end{array}\right)`$ (18)
$`\times `$ $`\chi _{s_{c_f}m_{c_f}^{}}](1)^{l_f+s_{v_f}M_{L_f}}\sqrt{2L_f+1}`$ (19)
$`\times `$ $`\left(\begin{array}{ccc}l_f& s_{v_f}& L_f\\ m_{l_f}& m_{v_f}& M_{L_f}\end{array}\right)(1)^{L_f+s_{c_f}M_f}\sqrt{2I_f+1}`$ (22)
$`\times `$ $`\left(\begin{array}{ccc}L_f& s_{c_f}& I_f\\ M_{L_f}& m_{c_f}& M_f\end{array}\right)Y_{l_fm_{l_f}}^{}(\widehat{\stackrel{}{k}}){\displaystyle \frac{1}{r}}g_{(L_fs_{c_f})I_f}(k,r),`$ (25)
where $`\stackrel{}{k}`$ is the wave vector associated with the relative motion between the two fragments in the continuum. $`l_f`$ is the orbital angular momentum for the core-valence relative motion and $`I_f`$ is the total spin of the $`c+v`$ system in the final channel.

We assume that the electromagnetic multipole operators act only on the relative motion variables of the system. This leads to the following conditions: $`s_{c_i}=s_{c_f}=s_c`$, $`s_{v_i}=s_{v_f}=s_v`$ and $`m_{c_i}=m_{c_f}^{}`$, $`m_{v_i}=m_{v_f}^{}`$.

Now we define a reduced matrix element of the multipole operator $`M(El,m)(=[Z_ce(\frac{A_v}{A_a})^l+(1)^lZ_ve(\frac{A_c}{A_a})^l]r^lY_{lm}(\widehat{\stackrel{}{r}})`$ with multipolarity $`l`$ having projection $`m`$) between the final and the initial state:

$`(L_is_c)I_iM(El)(L_fs_c)I_f=(1)^{l_i}`$ (26)
$`\times `$ $`\sqrt{{\displaystyle \frac{(2l_i+1)(2l_f+1)(2l+1)}{4\pi }}}\left(\begin{array}{ccc}l_i& l& l_f\\ 0& 0& 0\end{array}\right)`$ (29)
$`\times `$ $`[Z_ce({\displaystyle \frac{A_v}{A_a}})^l+(1)^lZ_ve({\displaystyle \frac{A_c}{A_a}})^l]`$ (30)
$`\times `$ $`{\displaystyle _0^{\mathrm{}}}𝑑rf_{(L_is_c)I_i}^{}(r).r^l.g_{(L_fs_c)I_f}(k,r).`$ (31)
In Eq.(5), $`Z_i(i=c,v)`$ is the charge of the $`i`$th fragment and $`A_j(j=c,v,A)`$ is the mass number of the $`j`$th partcle.

With this definition, the matrix element of the multipole operator is given by

$`\psi _{I_iM_i}M(El,m)\psi _{\stackrel{}{k}s_{v_f}m_{v_f}s_{c_f}m_{c_f}}`$ (32)
$`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{\genfrac{}{}{0pt}{}{l_im_{l_i}}{L_im_{L_i}}}{m_{v_i}m_{c_i}}}{}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{\genfrac{}{}{0pt}{}{I_fM_f}{l_fm_{l_f}m_{l_f}^{}}}{L_fM_{L_f}M_{L_f}^{}}}{}}(1)^{m_{v_i}+L_i+2L_f+3s_v+3s_cM_iM_{L_f}}`$ (33)
$`\times `$ $`(1)^{M_{L_f}^{}2M_f}\sqrt{(2L_i+1)(2I_i+1)}(2L_f+1)`$ (34)
$`\times `$ $`(2I_f+1)\left(\begin{array}{ccc}l_i& s_v& L_i\\ m_{l_i}& m_{v_i}& M_{L_i}\end{array}\right)\left(\begin{array}{ccc}L_i& s_c& I_i\\ M_{L_i}& m_{c_i}& M_i\end{array}\right)`$ (39)
$`\times `$ $`\left(\begin{array}{ccc}l_f& s_v& L_f\\ m_{l_f}^{}& m_{v_i}& M_{L_f}^{}\end{array}\right)\left(\begin{array}{ccc}L_f& s_c& I_f\\ M_{L_f}^{}& m_{c_i}& M_f\end{array}\right)`$ (44)
$`\times `$ $`\left(\begin{array}{ccc}l_f& s_v& L_f\\ m_{l_f}& m_{v_f}& M_{L_f}\end{array}\right)\left(\begin{array}{ccc}L_f& s_c& I_f\\ M_{L_f}& m_{c_f}& M_f\end{array}\right)`$ (49)
$`\times `$ $`\left(\begin{array}{ccc}l_i& l& l_f\\ m_{l_i}& m& m_{l_f}^{}\end{array}\right)(L_is_c)I_iM(El)(L_fs_c)I_f`$ (52)
$`\times `$ $`Y_{l_fm_{l_f}}^{}(\widehat{\stackrel{}{k}}).`$ (53)
The expression for the triple differential cross section for the above electromagnetic breakup process in the incident beam coordinate system (the spins are not observed) is given by

$`{\displaystyle \frac{d^3\sigma }{d\mathrm{\Omega }_{cv}d\mathrm{\Omega }_{cm}dE_{cv}}}={\displaystyle \frac{d\sigma ^{el}}{d\mathrm{\Omega }_{cm}}}\rho (E_{cv})P(\widehat{𝐤}),`$ (54)
where $`\frac{d\sigma ^{el}}{d\mathrm{\Omega }_{cm}}`$ is the Rutherford scattering cross section and $`\rho (E_{cv})`$ is the density of final states at relative energy $`E_{cv}`$. $`P(\widehat{𝐤})`$, obtained from the matrix element in Eq.(6) above, is associated with the probability that the two fragments are scattered in the final channel with relative motion wave vector $`\stackrel{}{k}`$. We expand it in terms of the spherical harmonic $`Y_{LM}(\widehat{𝐤})`$. It is given by

$`P(\widehat{𝐤})={\displaystyle \underset{LM}{}}A_{LM}Y_{LM}(\widehat{𝐤}),`$ (55)
where

$`A_{LM}`$ $`=`$ $`\left({\displaystyle \frac{4\pi Z_Ae}{\mathrm{}v}}\right)^2{\displaystyle \underset{\genfrac{}{}{0pt}{}{lml^{}m^{}}{L_iL_i^{}L_fL_f^{}}}{}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{l_il_i^{}}{\genfrac{}{}{0pt}{}{l_fl_f^{}}{I_fI_f^{}}}}{}}(1)^{s_c+s_vl_il_i^{}L_iL_i^{}}`$ (56)
$`\times `$ $`(1)^{I_f+I_f^{}I_i+m}^{(l+l^{})}{\displaystyle \frac{(2I_f+1)(2I_f^{}+1)}{(2l+1)(2l^{}+1)}}`$ (57)
$`\times `$ $`\sqrt{(2L_i+1)(2L_i^{}+1)}(2L_f+1)(2L_f^{}+1)`$ (58)
$`\times `$ $`\sqrt{{\displaystyle \frac{(2l_f+1)(2l_f^{}+1)(2L+1)}{4\pi }}}`$ (59)
$`\times `$ $`[{\displaystyle \underset{nn^{}}{}}D_{mn}^l({\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi \theta _{cm}}{2}})D_{m^{}n^{}}^l^{}({\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi \theta _{cm}}{2}})`$ (60)
$`\times `$ $`Y_{ln}({\displaystyle \frac{\pi }{2}},0)Y_{l^{}n^{}}({\displaystyle \frac{\pi }{2}},0)I_{ln}(\theta _{cm},\xi )I_{l^{}n^{}}(\theta _{cm},\xi )]`$ (61)
$`\times `$ $`\left(\begin{array}{ccc}l_f& l_f^{}& L\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccc}l& l^{}& L\\ m& m^{}& M\end{array}\right)\left\{\begin{array}{ccc}l& l^{}& L\\ I_f^{}& I_f& I_i\end{array}\right\}`$ (68)
$`\times `$ $`\left\{\begin{array}{ccc}l^{}& L_i^{}& L_f^{}\\ s_v& l_f^{}& l_i^{}\end{array}\right\}\left\{\begin{array}{ccc}l& L_i& L_f\\ s_v& l_f& l_i\end{array}\right\}\left\{\begin{array}{ccc}L_f& L_f^{}& L\\ l_f^{}& l_f& s_v\end{array}\right\}`$ (75)
$`\times `$ $`\left\{\begin{array}{ccc}l& I_i& I_f\\ s_c& L_f& L_i\end{array}\right\}\left\{\begin{array}{ccc}l^{}& I_i& I_f^{}\\ s_c& L_f^{}& L_i^{}\end{array}\right\}\left\{\begin{array}{ccc}I_f& I_f^{}& L\\ L_f^{}& L_f& s_c\end{array}\right\}`$ (82)
$`\times `$ $`(L_is_c)I_iM(El)(L_fs_c)I_f`$ (83)
$`\times `$ $`(L_i^{}s_c)I_iM(El^{})(L_f^{}s_c)I_f^{}^{}.`$ (84)
In Eq.(9), $``$ is half the distance of closest approach in a head-on collision and $`\xi `$ is the adiabaticity parameter, given as the ratio of the collision time and the excitation time. $`D_{mn}^l`$’s are the Wigner $`D`$-functions and $`\theta _{cm}`$ is the scattering angle of the c.m. of the projectile. $`Z_A`$ is the charge number of the target and $`v`$ is the projectile-target relative velocity in the entrance channel. $`I_{ln}(\theta _{cm},\xi )`$ is the classical orbital integral in the focal system of the hyperbolic orbit of the projectile .

Through our formalism, we can account for pure dipole ($`l=l^{}=1`$) and pure quadrupole ($`l=l^{}=2`$) transitions as well as mixed transitions ($`l=1,l^{}=2`$ or $`l=2,l^{}=1`$) or dipole-quadrupole interference. Also, it is possible to calculate very exclusive observables up to the level of the triple differential cross section. Previous calculations on the breakup of 8B were done by assuming that the angular distribution of fragments is isotropic in the projectile rest frame . This approximation gives

$`{\displaystyle \frac{d^3\sigma }{d\mathrm{\Omega }_{cv}d\mathrm{\Omega }_{cm}dE_{cv}}}={\displaystyle \frac{1}{4\pi }}{\displaystyle \frac{d^2\sigma }{d\mathrm{\Omega }_{cm}dE_{cv}}}.`$ (85)
It should be noted that expression for similar triple differential cross section as in Eq.(7) was given by Baur and Weber in Ref. , through which it is possible to account for $`E1`$ and $`E2`$ contributions separately as well as their interference. However, in contrast to our work, the coupling scheme of angular momenta and spins followed by them does not allow the use of the more detailed coupling scheme adopted here. In fact, the coupling scheme in Ref.  is more restricted in the sense that the orbital angular momentum of the valence particle with respect to the core is coupled to the sum of the spins of $`c`$ and $`v`$ both in the initial and the final channel.

We note that integration over the solid angle associated with the relative motion of the fragments gives

$`{\displaystyle \frac{d^2\sigma }{d\mathrm{\Omega }_{cm}dE_{cv}}}=\sqrt{4\pi }{\displaystyle \frac{d\sigma ^{el}}{d\mathrm{\Omega }_{cm}}}\rho (E_{cv})A_{00},`$ (86)
which is free from the $`E1E2`$ interference term, since $`A_{00}`$ (with $`L,M`$ = 0) does not involve these terms. This is evident from the 3-j symbol $`\left(\begin{array}{ccc}l& l^{}& L\\ m& m^{}& M\end{array}\right)`$ occurring in Eq.(9), because for $`L`$ = 0, $`l`$ and $`l^{}`$ must be the same. Thus, any cross section obtained from the above triple differential cross section (Eq.(7)) by integration with respect to the solid angle $`\mathrm{\Omega }_{cv}`$ is free from $`E1E2`$ interference. This result is independent of the structure models of 8B. Therefore, the analyses presented in Refs.  are indeed free from the $`E1E2`$ interference effects.

## III STRUCTURE MODEL

The calculations of the reduced matrix elements involved in Eq.(9) require the detailed knowledge of ground state as well as continuum structure of 8B, which is not yet known with certainty. In our calculations, we adopt the single particle potential model (SPPM) for 8B . It is to be noted that the matrix elements of the multipole operators enter directly into these calculations. Therefore, they are quite sensitive to the structure model of 8B.

Within the SPPM, the valence proton (with spin $`\frac{1}{2}`$) in 8B (with spin-parity 2+) is assumed to move relative to an inert 7Be core (with intrinsic spin-parity 3/2-) in a Coulomb field and a Woods-Saxon plus spin-orbit potential, with an adjustable depth $`V_0(l(Ls_c)I)`$ for the initial and each final channel:

$`V(r)=V_0(l(Ls_c)I)(1F_{s.o.}(\stackrel{}{l}.\stackrel{}{s}){\displaystyle \frac{r_0}{r}}{\displaystyle \frac{d}{dr}})f(r),`$ (87)
where

$`f(r)=(1+\mathrm{exp}((rR)/a))^1.`$ (88)
Adjusting the depth allows one to reproduce the energy of the known states. We use $`a`$=0.52 fm, $`r_0`$=1.25 fm and $`R`$=2.391 fm . The spin-orbit strength is set to $`F_{s.o.}`$ = 0.351 fm . The rms distance of the core-proton relative motion and the rms size of 8B come out to be 4.24 fm and 2.64 fm respectively.

The well depth for the ground state channel, $`V_0(l_i(L_is_c)I_i)`$ = $`V_0((p_{3/2},3/2^{})2^+)`$, was adjusted to reproduce the one-proton separation energy of 0.137 MeV. It came out to be $``$44.658 MeV). Similarly, the observed $`I_f`$ = 1+ and 3+ resonances in 8B are described as $`p_{3/2}`$ waves coupled to the ground state of the core, and the well depths for these channels, $``$42.14 and $``$36.80 MeV, respectively, have been adjusted to reproduce the known resonance energies (0.637 and 2.183 MeV respectively). A $`p_{3/2}`$ wave and the spin of the core can also couple to the total spin 0+. But we ignore this channel since it appears to be very weak in the low-lying excitation spectrum of 8B. For all other partial waves ($`s_{1/2},p_{1/2},d_{3/2}`$ etc.) we choose identical well depths and set them equal to the value $``$42.14 MeV obtained for the $`(p_{3/2},3/2^{})1^+`$ channel, as suggested by Robertson .

## IV RESULTS and DISCUSSIONS

To check the accuracy of our formulation, we present in Fig. 1 a comparison of our calculations with the data  for the triple differential cross sections for the reaction 6Li + Pb $``$ $`\alpha `$ \+ $`d`$ \+ Pb at 156 MeV beam energy as a function of the relative energy between the alpha particle and the deuteron for $`\theta _\alpha =\theta _d=`$3. Only $`E2`$ excitation contributes in this case. Since the multipole charge for the dipole case is zero for 6Li (Eq.(5)), the $`E1`$ contribution is zero. We have assumed 6Li to be a cluster of $`\alpha `$ particle and deuteron, for which the structure model

of Ref.  has been used. We note that our calculations are in good agreement with the data. It may be noted that for the results presented in Ref.  for the same reaction, the structure part has been obtained from a constant astrophysical $`S`$-factor of 1.7$`\times `$10-5 MeV.mb, and not by using proper wave functions for the ground state and excited states of 6Li as has been done by us. In these data the contributions from nuclear excitation effects are negligible as the freagments have been detected at very forward angles .

In Fig. 2, we present the results of our calculations for the triple differential cross section (Eq.(7)) for Coulomb breakup of 8B on a Pb target at 46.5 MeV/nucleon beam energy as a function of $`\theta _{cm}`$ for relative energy $`E_{cv}`$ = 0.2 MeV and $`\theta _{cv}`$ = 1 (top half) and as a function of $`E_{cv}`$ for $`\theta _{cm}`$ = 3 and $`\theta _{cv}`$ = 1 (bottom half). We see that the $`E1E2`$ interference (long dashed line) contribution is quite important in both cases. In fact, it is larger than the $`E2`$ contribution (dashed line) and modifies the coherent sum (solid line) of $`E1`$ (dotted line) and $`E2`$ cross sections significantly. It may be noted that this type of triple differential cross sections have not yet been measured for reactions involving unstable radioactive nuclei (e.g. 8B).

However, the differential cross section $`\frac{d\sigma }{d\theta _{cm}}`$ has been measured as a function of $`\theta _{cm}`$ in the dissociation of 8B on Pb at beam energies $``$50 MeV/nucleon at RIKEN . At larger scattering angles the cross sections are more sensitive to the $`E2`$ component. A detailed investigation of this reaction was carried out in Ref.  by assuming the fragment emission to be isotropic in the projectile rest frame (Eq.(10)).

In Fig. 3, we present our calculations for the different components in the Coulomb breakup cross section for this reaction at the beam energy of 51.9 MeV/nucleon, by using our formalism (i.e. without making the approximation of isotropic emission). For this observable the coherent sum (solid line) of $`E1`$ (dotted line) and $`E2`$ (dashed line) contributions is simply the sum of these two separate cross sections, as there is no contribution from the $`E1E2`$ interference component as discussed above (see Eq.(11)). We note that the results presented in Fig. 3 are the same as the pure Coulomb dissociation (with semiclassical theory) results presented in Ref. , which have been obtained with the isotropic emission assumption. These calculations, however, use a different structure model for 8B, namely, the shell model embedded in the continuum. Therefore, the present calculations give credence to the calculations reported in Ref. . As discussed earlier in Refs. , the disagreement between the data and the calculations beyond 4 in this figure can be attributed to the point like projectile approximation of the semiclassical theory, which is no longer valid for larger angles. Inclusion of finite size effects in the calculations leads to better agreement between data and the theoretical results . It should also be mentioned that the nuclear breakup effects are important only at larger angles ($`>4^{}`$) in all the three energy bins .

The triple differential cross section $`\frac{d^3\sigma }{d\mathrm{\Omega }_{cv}d\mathrm{\Omega }_{cm}dE_{cv}}`$ can be related to that of the individual fragments $`c`$ and $`v`$ ($`\frac{d^3\sigma }{dE_cd\mathrm{\Omega }_cd\mathrm{\Omega }_v}`$) . Since interference contributions are significant in $`\frac{d^3\sigma }{d\mathrm{\Omega }_{cv}d\mathrm{\Omega }_{cm}dE_{cv}}`$, it is expected that they would also be important in $`\frac{d^3\sigma }{dE_cd\mathrm{\Omega }_cd\mathrm{\Omega }_v}`$ in general. This is, indeed, the case as can be seen in Fig. 4. In this figure, we have plotted $`\frac{d^3\sigma }{dE_cd\mathrm{\Omega }_cd\mathrm{\Omega }_v}`$ (with $`c{}_{}{}^{7}Be`$ and $`vp`$) as a function of $`\theta _{7_{Be}}`$ for Coulomb breakup of 8B on 58Ni at the beam energy of 25.8 MeV, for $`E_{7_{Be}}`$ = 22.575 MeV (the beam velocity energy), $`\theta _p`$ = 20 and $`\varphi _p`$ = 0. We see that the $`E1E2`$ interference term is quite significant and is even larger than that of the dipole component.

However, for less exclusive observables (like the angular distribution of individual fragments) the interference term is not so strong. This can be seen in Fig. 5, where we show the 7Be angular distribution resulting from the Coulomb dissociation of 8B on a Ni target at 25.8 MeV. We note that the $`E1E2`$ interference (long dashed line) is small in magnitude compared to the dipole (dotted line) and quadrupole (dashed line) breakup cross sections, excepting at larger angles (around 75) where it is comparable in magnitude to the $`E1`$ component. The interference component oscillates between positive and negative values which is also reflected in the coherent sum of the $`E1`$ and $`E2`$ cross sections (solid line).

To understand why the interference is not large in the 7Be angular distribution shown in Fig. 5, we have plotted in Fig. 6 the $`\theta _p`$ (top half, with typical value of $`\varphi _p`$ = 10) and $`\varphi _p`$ (bottom half, with typical value of $`\theta _p`$ = 20) variations of the interference term in the triple differential cross section shown in Fig. 4 for typical energy and angle of the 7Be fragment $`E_{7_{Be}}`$ = 22.575 MeV and $`\theta _{7_{Be}}`$ = 10 respectively. The variation patterns are oscillatory and show positive and negative values of almost equal magnitudes. For some angles the cross sections are very small. Therefore, it is no surprise that integrations with respect to the polar and azimuthal angles of the proton will lead to cancellation and hence to small magnitudes of the interference term in the overall 7Be angular distribution.

However, there is still a remaining question of the effects of strong (nuclear) interaction between the colliding nuclei. They can be almost eliminated by choosing the proper kinematical conditions. For example, the 8B breakup data taken by the RIKEN collaboration  have been shown  to be almost free from the nuclear effects for 7Be-$`p`$ relative energies $``$ 0.75 MeV at very forward angles ($``$ 4). Similarly the data taken at GSI  at higher beam energy of $``$ 250 MeV/A are also free from these effects (see $`e.g.`$ ). Therefore, the results presented in Fig. 2 are unlikely to be affected by nuclear excitation. However, these effects can be quite strong  for the breakup studies  at lower beam energies ($``$ 25 MeV). In any case, it is straightforward to calculate  the amplitudes for pure nuclear and Coulomb-nuclear interference terms for the dipole and quadrupole excitations as well as of their interference.

## V SUMMARY and CONCLUSIONS

In summary, we have performed first order semiclassical calculations of the electromagnetic excitation of 8B, where the electric dipole and quadrupole excitation components as well as their interference effects are included. The theory permits consideration of coupling of angular momenta and spins of the projectile fragments in detail. It is possible to calculate exclusive observables to the level of the triple differential cross section within this formalism.

We find that the magnitude of the $`E1E2`$ interference term could be appreciable in the triple differential cross sections. However, observables which are not functions of the solid angle associated with the relative motion of the breakup fragments in the final channel, are free from this term. This result is independent of the structure models of 8B. We have also shown that for double differential cross sections involving angle of scattering of the projectile with respect to the target and energy of relative motion between the projectile fragments, the approximation of the isotropic emission in the rest frame of the projectile is quite good. Therefore, analysis of the RIKEN data presented earlier using this approximation is quite accurate and this data, in the proper kinematical regime ($`\theta _{{}_{}{}^{8}B_{}^{}}<`$4 and $`E_{{}_{}{}^{7}Bep}`$500 keV) can be used to extract rather reliable astrophysical $`S`$-factor $`S_{17}`$(0).

The interference terms are also significant in the triple differential cross sections of the individual fragments. However, in the angular distribution of the individual fragments these terms are not significant. Therefore, the experimental data present in Ref.  are almost free from the $`E1E2`$ interference terms. However, the effects of the three body kinematics (not considered in the analysis of these data) may still be important.


format: text

text:
# CRITICAL BEHAVIOUR IN THE DENSE PLANAR NJL MODEL

\[

## Abstract

We present results of a Monte Carlo simulation of a 2+1 dimensional Nambu – Jona-Lasinio model including diquark source terms. A diquark condensate $`qq`$ is measured as a function of source strength $`j`$. In the vacuum phase $`qq`$ vanishes linearly with $`j`$ as expected, but simulations in a region with non-zero baryon density suggest a power-law scaling and hence a critical system for all $`\mu >\mu _c`$. There is no diquark condensation signalling superfluidity. Comparisons are drawn with known results in two dimensional theories, and with the pseudogap phase in cuprate superconductors. We also measure the dispersion relation $`E(k)`$ for fermionic excitations, and find results consistent with a sharp Fermi surface. Any superfluid gap $`\mathrm{\Delta }`$ is constrained to be much less than the constituent quark mass scale $`\mathrm{\Sigma }_0`$.

preprint: SWAT/00/261

\]

In the absence of a reliable non-perturbative method for calculating the properties of QCD at high densities, models of strongly-interacting matter continue to rely on results derived from effective theories such as the Nambu – Jona-Lasinio (NJL) model . Recently such approaches have suggested the intriguing possibility of superconducting behaviour in quark matter at high density as a result of BCS condensation between diquark pairs at the Fermi surface . Even in simplified models, however, systematic methods of calculation are hard to find, and up to now mostly self-consistent methods have been applied . Scenarios such as color-flavor locking , however, imply that in the high density phase the global U(1)B symmetry in the QCD Lagrangian corresponding to conservation of baryon number is spontaneously broken, leading to superfluid behaviour and a massless scalar excitation via Goldstone’s theorem. A description of the fluctuations of these modes and the resulting long-range interactions may need a more systematic non-perturbative approach; in this Letter we will use numerical results from lattice field theory to argue that the situation is more complex than hitherto thought.

The Lagrangian for the NJL model, incorporating a baryon chemical potential $`\mu `$, is

$$_{NJL}=\overline{\psi }(/+m+\mu \gamma _0)\psi \frac{g^2}{2}[(\overline{\psi }\psi )^2(\overline{\psi }\gamma _5\stackrel{}{\tau }\psi )^2],$$
(1)
where $`\psi ,\overline{\psi }`$ carry a global isopsin index acted on by Pauli matrices $`\tau _i`$. The model in 2+1 dimensions has been studied with $`\mu 0`$ using staggered lattice fermions $`\chi ,\overline{\chi }`$ : the reader is referred to this paper for further details of the formulation and symmetries of the lattice model. Apart from the obvious numerical advantages of working in a reduced dimensionality, the model in this case also has an interacting continuum limit for a critical value of the coupling $`1/g_c^21`$ , making controlled predictions possible in principle. In the simulations of this Letter, as in , $`1/g^2=0.5`$, rather far from the continuum limit. For $`\mu =0`$, this results in a ground state with broken chiral symmetry and a dynamically generated fermion mass $`\mathrm{\Sigma }_00.71`$, together with a triplet of Goldstone pions. The “constituent quark mass” $`\mathrm{\Sigma }_0`$ defines the physical scale in cutoff units. As $`\mu `$ is raised from zero, the system remains essentially unaltered until $`\mu =\mu _c0.65`$, whereupon the constituent quark mass $`\mathrm{\Sigma }`$ falls to zero as chiral symmetry is restored in a strong first-order transition, and baryon density $`n=\overline{\psi }\gamma _0\psi `$ jumps discontinuously from zero to a non-zero value; by $`\mu =0.8`$, the density $`n0.25`$ quarks of each isopsin component per lattice site.

The tendency for diquark pairing at high density was investigated via measurements of the diquark propagator $`G(t)=_𝐱qq(0,\mathrm{𝟎})\overline{q}\overline{q}(t,𝐱)`$ . In particular it was found that for $`\mu >\mu _c`$ and the choice of isoscalar diquark operator $`qq=\chi ^{tr}(x)\tau _2\chi (x)`$ (the “spectral scalar” channel of ), $`G(t)`$ plateaus at a non-zero value as $`t\mathrm{}`$, implying $`qq\overline{q}\overline{q}0`$ by the cluster property. Studies on varying spatial volumes, however, reveal that the plateau height is not an extensive quantity as naively expected, and that diquark condensation signalling a superfluid state is not unambiguously indicated.

The situation can be clarified by the introduction of source terms for diquark and anti-diquark pairs:

$$[j,\overline{ȷ}]=_{NJL}+j\chi ^{tr}(x)\tau _2\chi (x)+\overline{ȷ}\overline{\chi }(x)\tau _2\overline{\chi }^{tr}(x).$$
(2)
The functional integral may now be written using the Gor’kov representation  as

$$Z[j,\overline{ȷ}]=\mathrm{Pf}(𝒜[j,\overline{ȷ}]),$$
(3)
where $`\mathrm{}`$ denotes averaging with respect to bosonic auxiliary fields $`\sigma `$ and $`\stackrel{}{\pi }`$ introduced to make $`_{NJL}`$ bilinear in $`\chi `$ and $`\overline{\chi }`$, and the antisymmetric matrix $`𝒜`$ is

$$𝒜=\left(\begin{array}{cc}\overline{ȷ}\tau _2& \frac{1}{2}M\\ \frac{1}{2}M^{tr}& j\tau _2\end{array}\right),$$
(4)
where $`M=M(\mu ;\sigma ,\stackrel{}{\pi })`$ is the conventional kinetic operator of $`_{NJL}`$ . For convenience we now define diquark operators $`qq_\pm `$ via

$$qq_\pm (x)=\frac{1}{2}\left[\chi ^{tr}(x)\tau _2\chi (x)\pm \overline{\chi }(x)\tau _2\overline{\chi }^{tr}(x)\right],$$
(5)
with corresponding sources $`j_\pm =j\pm \overline{ȷ}`$. The diquark condensate is readily calculated to be

$$qq_+=\frac{1}{V}\frac{\mathrm{ln}Z}{j_+}=\frac{1}{4V}\mathrm{tr}\tau _2𝒜^1.$$
(6)
Under U(1)B transformations $`\chi e^{i\beta }\chi `$, $`\overline{\chi }\overline{\chi }e^{i\beta }`$, the operators $`qq_\pm `$ rotate into each other; condensation of $`qq_+0`$ implies a massless Goldstone mode in the $`qq_{}`$ channel as $`j_+0`$, illustrated by a Ward identity for the susceptibility $`\chi _{}`$:

$$\chi _{}\underset{x}{}qq_{}(0)qq_{}(x)|_{j_{}=0}=\frac{qq_+}{j_+}.$$
(7)
Numerical simulation using the functional measure (3) requires us to address the issue of calculating the Pfaffian . We begin by observing that $`\mathrm{Pf}𝒜=\pm \mathrm{det}^{\frac{1}{2}}𝒜`$, the sign depending on the details of the ordering of the Grassmann integration. Now, using the property of a square block matrix

$$\mathrm{det}\left(\begin{array}{cc}X& Y\\ W& Z\end{array}\right)=\mathrm{det}X\mathrm{det}(ZWX^1Y),$$
(8)
together with the property $`\tau _2M\tau _2=M^{}`$, we find $`\mathrm{det}𝒜=\mathrm{det}(j\overline{ȷ}+M^{}M/4)`$, and hence is real and positive if $`j\overline{ȷ}`$ is chosen real and positive. It follows that $`\mathrm{Pf}𝒜`$ is real. To avoid the sign problem, however, we choose to simulate with measure $`\mathrm{det}^{\frac{1}{2}}(𝒜^{}𝒜)=\mathrm{Pf}^2𝒜`$, corresponding to a doubling in the number of fermion species. This is usual practice in simulations of four-fermi models; “crosstalk” between the two species, which might conceivably lead to problems for $`\mu 0`$ due to light unphysical baryonic states in the spectrum, is subleading in this case and causes no problems . It is straightforward to implement the simulation using a hybrid molecular dynamics ‘R’ algorithm  — the resulting dynamics are those of 2 flavors of staggered lattice fermion, corresponding to $`N_f=4`$ flavors of four-component physical fermion.

We have performed simulations on lattice sizes $`16^3`$, $`16^2\times 24`$, $`24^3`$, $`32^3`$ and $`48^3`$, with chemical potentials $`\mu =0,0.2`$ corresponding to the low density chirally-broken phase, and $`\mu =0.8`$ in the high density chirally-restored phase . Diquark source values ranging from $`j=0.3`$ down to as low as $`j=0.025`$ were used; throughout we took $`\overline{ȷ}=j`$ and kept a bare Dirac mass $`m=0.01`$ to assist identification of chirally broken and restored phases. We plot $`\mathrm{ln}qq_+`$ vs. $`\mathrm{ln}j`$ in Fig. 1. At low density the results for $`\mu =0`$ and $`\mu =0.2`$ are very similar and support a linear relation $`qq_+j`$, suggesting that U(1)B is unbroken in the zero-source limit. For $`\mu =0.8`$, when chiral symmetry is restored, the situation is different: whilst there is still no evidence for a non-zero condensate as $`j0`$, the log-log plot now displays a marked curvature, together with evidence for significant finite volume effects. Empirically we find that on a $`L_s^2\times L_t`$ system the dominant correction scales as $`1/L_t`$ – this enables an extrapolation of the data with $`L_t24`$ to the thermodynamic limit, shown in Fig. 2.

The extrapolated data, denoted by filled squares in Fig. 1, suggest the best description in the high-density phase is a power law:

$$qq_+j^\alpha ,$$
(9)
with $`\alpha 0.31`$. There is thus no superfluid condensate in this region of parameter space, and no consequent Goldstone pole. The Ward identity (7), however, implies that the transverse susceptibility $`\chi _{}`$ diverges as $`j^{\alpha 1}`$; thus long-wavelength fluctuations remain important, and there are strongly-interacting massless states present. This is the first main result of this Letter.

Now, the value $`\mu =0.8`$ was originally chosen merely as characteristic of the high-density phase; moreover data from partially quenched simulations  suggest that the power-law scaling (9) may be generic for $`\mu >\mu _c`$. This is analogous to the behaviour of the $`2d`$ X-Y model , which in accordance with well-known theorems  does not have spontaneous magnetisation, but instead exhibits a line of critical points for all $`T<T_c`$ . In this low-temperature phase the modulus of the spin variables is constant, but the phase fluctuates; the divergence of $`\chi _{}`$ signals correlation of the phase over macroscopic distance scales without long-range phase coherence. This critical behaviour without an order parameter can be characterised by two exponents $`\delta `$ and $`\eta `$, defined in our notation by $`qq_+j^{\frac{1}{\delta }}`$, and $`lim_{j_+0}qq_{}(0)qq_{}(x)1/x^{d2+\eta }`$. For the X-Y model  $`\delta 15`$ and $`\eta \frac{1}{4}`$ (the inequalites being saturated up to logarithmic corrections at $`T=T_c`$), and satisfy the hyperscaling relation

$$\delta =\frac{d+2\eta }{d2+\eta }.$$
(10)
For the dense planar NJL model, by contrast, our results suggest $`\delta 3`$.

The phenomenon of long-range phase coherence being wiped out by soft transverse fluctuations is particularly interesting for a composite order parameter field, suggesting as it does that some dynamical pairing mechanism is in play even in the absence of symmetry breaking. It has been suggested that such a state may describe a pseudogap phase, characterised by a suppression of spectral weight in the vicinity of the Fermi surface, and observed as a precursor to the superconducting state for underdoped cuprate superconductors . It is intriguing in the current context that cuprate superconductivity is a planar phenomenon ***An important distinction is that in cuprates the pairing operator is $`d`$-wave, whereas for the NJL model it is $`s`$-wave..

A field theoretic example of possible relevance is studied in , where results from an exactly soluble $`2d`$ fermionic model are generalised to the Gross-Neveu model with U(1) chiral symmetry . Once again due to the low dimensionality no long-range ordering via $`\overline{\psi }\psi 0`$ is possible; however for sufficiently large number of flavors $`N_f`$ the model has a phase analogous to the low temperature phase of the X-Y model, with $`\eta O(1/N_f)`$. The physical fermion is a superposition of positive and negative chirality states and has zero net chirality; despite the absence of chiral symmetry breaking it propagates as $`C(x)x^{(\frac{1}{2}+\eta )}e^{\mathrm{\Sigma }x}`$, implying massive propagation as $`x\mathrm{}`$, but also a non-trivial spectral function, signalled by the departure from canonical scaling of the pre-exponential factor.

These considerations motivate a study of the spectrum of spin-$`\frac{1}{2}`$ excitations in the NJL model by analysing the fermion propagator $`C(𝐤,t)=_𝐱\mathrm{Re}\chi (0,\mathrm{𝟎})\overline{\chi }(t,𝐱)e^{i𝐤.𝐱}`$, where for staggered fermions the sum only includes x which are an even number of lattice spacings from the origin in each direction. Note that for $`j0`$, the most general definition of $`C(t)`$ contains both “normal” $`\chi \overline{\chi }`$ and “anomalous” $`\chi \chi `$ components. In this Letter we analyse only the isosinglet normal component, with k oriented along a lattice axis. Fig. 3 shows possible spectra $`E(k)`$ for $`\mu >\mu _c`$, $`n>0`$ . Solid lines denote states which are occupied in the ground state; following common usage we label these “Dirac” and “Fermi” seas (note there is no quantum number which distinguishes between these branches). Dashed lines denote empty states, and hence describe the spectrum of allowed excitations. Note that for $`k<k_F`$ the lowest excitations vacate states in the Fermi sea, and hence are ‘hole-like’: in contrast for $`k>k_F`$, excitations add quarks to the system, and are ‘particle-like’. On the right hand side the curve is continuous as it crosses the axis; there is a sharp Fermi surface with well-defined Fermi momentum $`k=k_F`$, where excitations cost zero energy. On the left hand side we have sketched the effect of a discontinuity at the Fermi surface arising from a BCS instability; the lowest excitation energy now corresponds to the gap $`\mathrm{\Delta }`$.

We have measured $`C(𝐤,t)`$ on $`32^3`$ lattices at $`\mu =0.8`$. Empirically we find that for $`t`$ even it is consistent with zero, and that for $`t`$ odd it is well described by a 3-parameter fit $`C(𝐤,t)=Ae^{Et}+Be^{E(L_tt)}`$. For small values of $`k`$ we find $`AB`$, ie. the signal predominantly forwards-moving, and vice-versa for $`k`$ as approaches its maximum value $`\pi /2`$. As $`j`$ is decreased the switch between these two behaviours becomes increasingly pronounced. On comparison with the analytic form for free fermions, we label these cases as ‘hole’ and ‘particle’ respectively. We also find that $`E(k,j)`$ decreases monotonically and roughly linearly with $`j`$, and a linear extrapolation to $`j=0`$ can be made based on $`j[0.1,0.025]`$. The results, shown in Fig. 4, are in good qualitative agreement with the expectations of Fig. 3. Conservatively, we can bound any possible gap from above, and find $`\mathrm{\Delta }<0.05\mathrm{\Sigma }_00.7`$, in contrast with expectations based on the self-consistent gap equation approach of . More boldly, we plot a curve of the form $`E(k)=|P+Q(\frac{\pi }{2}k)^2|`$, which suggests it is plausible the dispersion relation is continuous near $`k=k_F`$, and hence $`\mathrm{\Delta }=0`$. The fit yields a $`k_F0.90`$ slightly larger than the Fermi energy $`\mu `$, which is characteristic of a system with attractive inter-particle forces. We conclude that there is a sharp Fermi surface in this region of parameter space: this is the second main result of the Letter.

In the future, it will be helpful to explore other values of $`\mu >\mu _c`$, to see whether $`\delta =\delta (\mu )`$. Independent estimates of $`\eta `$ from simulations at $`j=0`$ would also be valuable; one could then use hyperscaling (10) to probe the effective dimensionality of the system. It has been suggested that absence of diquark condensation in the planar NJL model is an artifact due to the low dimensionality of the Fermi surface in 2+1 dimensions ; it should be mentioned, however, that partially quenched data taken on small 3+1 dimensional systems also appear to scale according to (9) . Finally, by analogy with cuprate superconductivity, it may be the case that this critical phase we have found is a precursor to a superfluid phase setting in at higher densities. Simulations with the current parameters become plagued by lattice artifacts for $`\mu \stackrel{>}{}1.0`$ . To observe superfluid behaviour simulations closer to the continuum limit may be necessary.

###### Acknowledgements.

This work is supported by the TMR-network “Finite temperature phase transitions in particle physics” EU-contract ERBFMRX-CT97-0122; in addition SJH thanks the Institute for Nuclear Theory at the University of Washington for its hospitality and the Department of Energy and the Leverhulme Trust for partial support. We were helped by discussions with Mark Alford, Nick Dorey, Helen Fretwell, Tim Hollowood, Krishna Rajagopal, Dirk Rischke, Misha Stephanov, Costas Strouthos and Jac Verbaarschot.


text:
# A Model for Addition Spectra in Quantum Dots

## I Introduction

The combined effect of disorder and electron-electron interactions is an old theme, to which Michael Pollak made important early contributions. The topic is still very much alive. In particular, with the advance in fabrication of mesoscopic electronic systems, the new feature of the confinement, in one or more directions, has been added to the problem.

This paper deals with disordered quantum dots, i.e., small islands of interacting electrons, in the presense of the random potential of impurities. Let us denote by $`E_0(N)`$ the ground state energy of a dot with $`N`$ electrons and define the addition spectrum as $`E_0(N+1)E_0(N)\mu (N+1)`$. This quantitiy is a subject of extensive theoretical and experimental study. It is determined experimentally by measuring the conductance of the dot, as a function of the gate voltage $`V_G`$. The latter determines the electrostatic potential on the dot. For a dot weakly coupled to the leads, and at low temperatures, the conductance is negligible (Coulomb blockade), unless $`V_G`$ is egual to $`\mu (N)`$, at which point the blockade is lifted and an additional electron can enter the dot. Thus, the conductance exibits a sequence of peaks. The distance between two consecutive peaks is

$`\mathrm{\Delta }_2(N)=\mu (N+1)\mu (N)=E_0(N+1)2E_0(N)+E_0(N1).`$ (1)
The simplest possible description of the Coulomb blockade is given in terms of the charging model, also known as the constant-interaction model (for a recent review see Alhassid, 2000). It is assumed in this model that the added, $`(N+1)`$-th, electron occupies the $`(N+1)`$-th single particle state in the dot, with energy $`ϵ_{N+1}`$. The interaction between the added electron and each of the $`N`$ electrons in the dot is described by a constant potential $`V_C`$. These assumptions result in the following expressions for the addition spectrum $`\mu (N)`$ and peak spacing $`\mathrm{\Delta }_2(N)`$:

$`\mu (N+1)=V_CN+ϵ_{N+1},`$ (2)
$`\mathrm{\Delta }_2(N)=V_C+(ϵ_{N+1}ϵ_N).`$ (3)
It follows from Eq. (2), that fluctuations in $`\mathrm{\Delta }_2(N)`$ are determined by the fluctuations in the single-particle level spacing which, for a weakly disordered dot, are given by the Wigner-Dyson distribution. Furthermore, if $`N`$ is odd, then the $`(N+1)`$-th particle occupies a level which is alredy occupied by a particle with opposite spin, so that $`\mathrm{\Delta }_2=V_C`$. Thus, according to the charging model, the distribution $`P(\mathrm{\Delta }_2)`$ should exibit a bimodal structure, composed of a $`\delta `$-function at $`V_C`$ plus a Wigner-Dyson function (shifted by $`V_C`$). Experiments (Sivan et al., 1996, Patel et al., 1998, Simmel et al., 1999) , however, show an approximately Gaussian distribution for $`\mathrm{\Delta }_2`$, with a width considerably larger than the average single-particle level spacing $`\mathrm{\Delta }`$ and with no sign of a bimodal structure. The disagreement between experiment and the charging model led to a number of theoretical studies of $`P(\mathrm{\Delta }_2)`$, based on exact diagonalization of the microscopic Hamiltonian (Prus et al., 1996, Berkovits, 1998) as well as on various approximations such as RPA (Blanter et al., 1997, Berkovits and Altshuler, 1997), the extreme classical limit (Koulakov et al., 1997) and the self-consistent Hartree-Fock scheme (Walker et al., 1999, Cohen et al., 1999, Levit and Orgad, 1999), or the somewhat related configuration interaction method (Benenti et al., 2000). In this paper we propose, and study, a simple model for the addition spectrum and the peak spacing of a weakly disordered (or chaotic) quantum dot. The electron spin is ignored, except for a few comments at the end of the paper.

## II The Model

$`N`$ interacting electrons in a quantum dot are described by the Hamiltonian

$`H={\displaystyle \underset{i=1}{\overset{N}{}}}[{\displaystyle \frac{\mathrm{}^2}{2m}}_i^2+V(\stackrel{}{r}_i)]+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i<j}{}}U(\stackrel{}{r}_i\stackrel{}{r}_j),`$ (4)
where $`m`$ is the electron effective mass, $`V(\stackrel{}{r}_i)`$ includes the confining potential as well as the random potential of impurities, and $`U(\stackrel{}{r}_i\stackrel{}{r}_j)`$ is the interaction potential. Let us discuss first the two limits, of weak and strong interactions. In the first limit $`\mathrm{\Delta }_2`$ is equal to the single-particle level spacing and, thus, obey the Wigner-Dyson statistics. Small corrections, due to interactions, can be treated in perturbation theory (Prus et al., 1996).

In the opposite case of very strong interactions (the ”extreme classical limit”) the potential energy dominates , and the charges settle down in a certain configuration which minimizes the potential energy of the system. For some special values of $`N`$, called “magic numbers”, one can expect particularly stable, i.e. low energy, configurations (Koulakov and Shklovskii, 1998, Morris et al., 1996). Specific values of the magic numbers depend on the shape of the dot, the boundary conditions and the type of interactions. For instance, for a dot of an equilateral triangle shape these numbers are $`N=\mathrm{}(\mathrm{}+1)/2`$, where $`\mathrm{}`$ is an integer. Indeed, for these values of $`N`$ one can view the magic configuration inside the triangle as being cut out of an infinite Wigner crystal, which is known to crystallize into a triangular lattice. In what follows, we consider a square dot, with periodic boundary conditions, often used in numerical calculations. The triangular Wigner crystal, characteristic for an isotropic interaction, is clearly incompatible with the square shape of the dot, and finding out the magic configurations for this geometry becomes a formidable problem. However, for anisotropic interactions, or in the presence of an underlying periodic potential, the charges can crystallize into a square lattice. In this case the magic configurations should occur at $`N=\mathrm{}^2`$. In order to keep the picture as simple as possible, without seriously damaging the forthcoming reasoning, we shall stick to this example.

Consider, thus, a magic configuration of $`N=\mathrm{}^2`$ electrons, in a square dot of size $`L`$. The electrons inside the dot form a square lattice, with a lattice constant $`b=n^{1/2}`$, where $`n=NL^2`$ is the electron concentration. The lattice is pinned by the weak disorder in the dot. The main property of a magic configuration is its stability, i.e., it can accomodate few extra particles without a major rearrangement. An extra particle can cause some polarization effects, by pushing apart the charges in its neighborhood, but, basically, it has to reside in the center of a cell, formed by 4 lattice charges. Furthemore, the extra particles can somewhat lower their energy by making quantum hops between neighboring cells. Thus, the added particles can be viewed as ”defects” moving on the intersticials of a magic configuration. For an infinite perfect crystal this kind of defects have been discussed long ago by Andreev and Lifshitz (Andreev and Lifshitz, 1969) in their theory of quantum crystals (see also Tsiper and Efros, 1997). They proposed that defects (or vacancies) can move through the crystal coherently and can be treated as a gas of excitations. For a disordered mesoscopic quantum dot this picture of ”defects” becomes more arguable. It is clear, in particular, that when too many particles are added to a given magic configuration it is not appropriate to speak in terms of many defects moving on top of that particular configuration: one should rather speak about vacancies (holes) in the next magic configuration. For instance, for $`N=102`$ one can envisage two defects in the magic configuration $`N=100`$. However, for $`N=119`$ it is clearly more appropriate to envisage two vacancies in the next magic configuration, $`N=121`$ (compare to shells in atomic or nuclear physics). Moreover, somewhere in between, say, at $`N=110`$ a major rearrangement is likely to occur, and at this point it is hardly possible to speak in terms of noninteracting defects or vacancies. Inspite of its limitations, we will use this picture of defects, embedded in the corresponding magic configuration, in order to formulate the forthcoming phenomenological model.

The main point of the above discussion is that in both limits, of small and large $`U`$, the motion of a particle added to an interacting system is described in a single particle language: a bare electron, in the first case, and a ”defect” in a magic configuration, in the second. Next, we interpolate between the two limits and arrive to the following, single particle, model for the addition spectra in quantum dots. The model is defined on a lattice with $`2N`$ sites. This lattice originates from the strong interaction limit, as described in the previous paragraph. In this limit one sublattice ($`A`$) is occupied by the $`N`$ electrons of a magic configuration, while the other one ($`B`$) represents sites available for extra particles (Fig. 1). When an extra particle is placed on site $`i`$ of sublattice $`B`$, the energy of the system increases by $`W_i+U_i`$, where $`W_i`$ is the disorder related site energy, and

$`U_i={\displaystyle \underset{ji}{}}U(\stackrel{}{r}_i\stackrel{}{r}_j)`$ (5)
describes the interaction of the added electron with the electrons in the dot. Consider, as an example, a model potential (Walker et al. 1999) consisting of a short range part $`V_0\delta _{i,i+\eta }`$ ($`\eta `$ designates the nearest neighbors of site $`i`$) plus a constant term $`V_C`$, which describes the charging energy due to the long range part of the interaction. For this example $`U_i=NV_C+4V_0`$. More generally, one can replace $`4V_0`$ by a phenomenological parameter $`U`$, whose relation to the pair potential remains unspecified. Thus, in this ”extreme classical limit”, the $`(N+1)`$-th particle, added to the $`N`$-particle magic configuration, will reside on the site with the lowest value of $`W_i`$, on sublattice $`B`$. The corresponding increase in energy is $`\mu (N+1)=V_CN+U+W_i`$. The $`(N+2)`$-th particle will go into the site with the next lowest value of $`W_i`$, etc. Similarly, one can remove particles from a magic configuration and compare energies $`E_0(N)`$, $`E_0(N1)`$, etc. For instance, when one adds the $`N`$-th particle to the $`(N1)`$-th particle system, one fills in a vacancy in sublattice $`A`$. The main point is that the energy $`U_i`$, corresponding to a site on sublattice $`A`$, differs from that on sublattice $`B`$. This is because a site on sublattice $`A`$ $`(B)`$ is surrounded by empty (occupied) sites. We assign values $`U`$ and zero to sites on sublattices $`B`$ and $`A`$, respectively.

The full model for the fluctuations in addition spectra is obtained by ”switching on” the kinetic energy and, thus, interpolating between the limits of small and large $`U`$. It is given by the matrix

$`H_{ij}=W_i\delta _{ij}t_{ij}+U_i\delta _{ij}`$ (6)
where $`t_{ij}`$ accounts for hopping between sites $`i`$ and $`j`$. Energies $`U_i`$ assume two values: $`U`$ for $`iB`$ and $`0`$ for $`iA`$. The model concentrates on the fluctuations, i.e., the charging energy $`V_C`$, due to the long range part of the interaction, is not included (it can be added to $`\mathrm{\Delta }_2`$ at the end). The matrix $`H_{ij}`$ is of size $`2N\times 2N`$. Its $`(N+n)`$-th eigenvalue describes the energy (more precisely, its fluctuating part) of the $`n`$-th particle added to the system of $`N`$ particles, so that $`\mathrm{\Delta }_2`$ (with $`V_C`$ subtracted) is given by the difference between a pair of consecutive eigenvalues. Let us stress that, although the picture in the limit of very strong interactions, discussed above, was essential for formulating the model, this limit is quite remote from the experimentally relevant situation: the latter corresponds to the ”liquid phase”, rather than the solid one. In the liquid phase the ”magic numbers” do not play any essential role and the predictions of the model are expected to be quite robust, i.e., independent of the specific assumptions (e.g. a square lattice, rather than a triangular one) used in the ”extreme classical limit”. (If the number of electrons is varied in too wide a range, one should allow for a weak dependence of the parameters of the model on the number of electrons.)

In the weak interaction limit, when $`U0`$, the model reduces to the Anderson model for a particle propagating in a random potential. In this limit the lattice has no physical meaning and, simply, provides discretization of the continuous space. Eq. (5) defines a simple, single-particle model for addition spectra. The model captures the right physics in the two limits, of small and large $`U`$. The weakness of the model is that it is no more than a simple interpolation between the two limits. Such an interpolation is, clearly, a big leap — especially, since some new quantum phases might exist between the two limits (Chakravarty et al., 1999, Benenti et al., 1999). The model, however, is not intended to describe subtle correlations in the ground state of $`N`$ interacting particles. Its only purpose is to describe the energetics of added particles. In fact, a rather similar picture — with the same characteristic feature of two energy bands in the large-$`U`$ limit, merging into one band for smaller $`U`$ — is often used in qualitative discussions of the Hubbard model (Mott, 1990). Below we study the addition spectra within the model of Eq. (5).

## III Numerical Results

We assume that site energies, $`W_i`$, are uniformly distributed in the interval $`[W/2,W/2]`$ and take $`t_{ij}=t`$ for nearest neighbors (and zero otherwise). $`t`$ is related to the particle mass by $`t=\mathrm{}^2/2ma^2`$ where $`a=(2n)^{1/2}`$ is the lattice constant. This relation insures the correct value, $`\mathrm{\Delta }=2\pi \mathrm{}^2/mL^2`$, for the single particle level spacing. The disorder strength is measured by the parameter $`W/t`$ and the interaction strength by the value of $`U/t`$.

Let us discuss in more detail the limit of strong interaction and weak disorder, i.e., $`W<t<<U`$. The $`U`$-term in Eq. (5) produces two levels, separated by energy $`U`$. Each of the levels is $`N`$-fold degenerate. The disorder and the kinetic energy remove this degeneracy and broaden the levels into bands. The lower band is occupied by the $`N`$ electrons, so that extra electrons must go into the upper band. They reside in the low energy tail of that band, so that the conductance peak spacings, $`\mathrm{\Delta }_2`$, are given by the level spacings in the tail. Since the density of states in the tail is much smaller than that in the middle of the band, one can expect large (in units of $`\mathrm{\Delta }`$) fluctuations in $`\mathrm{\Delta }_2`$. Moreover, since the tail states do not obey the Wigner-Dyson statistics, the same must be true for the statistics of $`\mathrm{\Delta }_2`$. Similarly, removing electrons from the $`N`$-particle system amounts to creating holes in the high energy tail of the lower band. Note that the width of the two bands is determined not by the hopping amplitude $`t`$ in the absence of interactions but, rather, by an effective hopping amplitude $`\stackrel{~}{t}t^2/U`$. This is because hopping between sites on sublattice $`B`$, such as sites 1 and 2 in Fig. 1, requires visiting an intermediate site (e.g. site 3), which belongs to the (occupied) sublattice $`A`$. Such virtual visits are costly and contribute only in second order of the perturbation theory, producing the effective hopping amplitude $`\stackrel{~}{t}`$ for particles (holes) in the upper (lower) band. The estimate $`\stackrel{~}{t}`$ for the band width holds only as long as $`\stackrel{~}{t}>W`$, i.e. $`t>\sqrt{UW}`$. In the opposite case when $`W<t<\sqrt{UW}`$, the $`t`$-term in Eq. (5) can be neglected. In this case the bands have width $`W`$ and all eigenstates are localized.

Thus, the model in Eq. (5) brings out in a simple way the origin of large fluctuations in the peak spacings $`\mathrm{\Delta }_2`$. The point is that, for large $`U`$, the model exhibits two bands and the statistics of $`\mathrm{\Delta }_2`$ is controlled by the states deep in the tails of those bands. For small $`U`$ the two bands merge into one band, and the statistics of $`\mathrm{\Delta }_2`$ is determined by the single-particle levels in the middle of that band. This qualitative picture is supported by numerics. As an example consider a dot with a number of electrons near $`N`$=64. The situation is described, within our model, by a lattice of 128 sites, i.e. by a matrix of size $`128\times 128`$. Seperation between eigenvalues number $`(N+i)`$ and $`(N+i1)`$ gives the spacing between the corresponding conductance peaks (the eigenvalues are ordered in energy). Let us look, for example, at the seperation $`\mathrm{\Delta }_2`$ between eigenvalues number 67 and 66. We take $`t=1`$ so that disorder and interactions are described by the numbers $`W`$ and $`U`$ respectively. The distribution of $`\mathrm{\Delta }_2`$, for an ensemble of 3000 matrices, for $`W=4`$ and various values of $`U`$ is shown in Fig. 2 ($`\mathrm{\Delta }_2`$ is normalized to $`\mathrm{\Delta }`$, so that the $`x`$-axes is labeled by $`\delta `$=$`\mathrm{\Delta }_2`$/$`\mathrm{\Delta }`$). For $`U=0`$ the distribution is close to Wigner-Dyson. This fact demonstrates that $`W=4`$, although not small compared to unity, corresponds to the diffusive regime (to reach diffusion behavior with small $`W`$ one would need larger matrices, i.e., larger number of electrons in the dot). For $`U`$ near 2 significant deviations from Wigner-Dyson statistics start to show up, and between $`U2`$ and $`U3`$ the distribution rapidly changes from Wigner-Dyson towards a more symmetric, broader distribution. For still larger values of $`U`$ (roughly, larger than 10) the distribution gradually approaches a limiting Poissonian shape. The origin of this limiting distribution is clear: as was explained above, for sufficiently large $`U`$ the $`t`$-term in Eq. (5) can be neglected, and one ends up with an ensemble of diagonal matrices $`H_{ij}=(W_i+U_i)\delta _{ij}`$, with statistically independent eigenvalues. Histograms similar to those in Fig. 2 can be obtained also for other pairs of consecutive eigenvalues, not too far from the 64-th eigenvalue, $`ϵ_{64}`$ (otherwise one should change the size of the matrix, adjusting it to the nearest magic number). This eigenvalue is somewhat exceptional, in the sence that it is between $`ϵ_{64}`$ and $`ϵ_{65}`$ that a gap opens up, for a sufficiently large $`U`$ (see Fig. 3, for $`U`$=5). Therefore, for large $`U`$, the distribution of $`\mathrm{\Delta }_2`$ for this case is shifted by $`U`$ (Fig. 4). The other special feature of this particular case is that, for very large values of $`U`$, the distribution assumes a semi-Poisson shape (a semi-Poisson distribution is defined as $`p(x)=xe^x`$).

The average density of states, for the same ensemble of matrices, is plotted in Fig. 3. For $`U=0`$ one has the single band of the Anderson model. When $`U`$ increases, a gap starts to form in the middle of the band and, eventually, two separate bands emerge. There is a clear correlation between formation of the gap and deviation of the peak spacing statistics from the Wigner-Dyson shape. Note that significant deviations from Wigner-Dyson statistics develop long before a genuine gap oppens up in the spectrum, i.e., standard metallic conductivity can coexist with ”non-metallic” thermodynamic properties.

## IV A Random Matrix Model

.

One can further simplify the model, replacing the Anderson matrix, i.e., the first two terms in Eq. (5), by a $`2N\times 2N`$ Gaussian orthogonal ensemble (GOE) (Mehta, 1991). This leads to a matrix model

$`H=G+A`$ (7)
where $`G`$ is taken from the GOE and $`A`$ is a $`2N\times 2N`$ (non-random) diagonal matrix, with half of the eigenvalues equal to $`U`$ and the other half equal to zero. It is interesting that the same matrix model can be arrived at by entirely different, more abstract arguments. Let’s start with a standard Gaussian matrix, $`G`$, which is known to give a good description of the single particle level statistics in a chaotic or disordered quantum dot. The matrix ensemble $`G`$ has a $`O(M)`$ symmetry, where $`M`$ is the number of states (orbitals) spanning the relevant Hilbert space. Consider now $`N`$ interacting particles in the dot and employ some self-consistent description, e.g. the Hartree-Fock scheme. The outcome will be $`M`$ self-consistent orbitals, $`N`$ of which are occupied and the rest, $`MN`$, are empty. One can make orthogonal transformations seperately among the $`N`$ occupied and the $`MN`$ empty orbitals, but not between the two types of orbitals. Thus, the symmetry is broken from $`O(M)`$ to $`O(N)\times O(MN)`$. The matrix $`A`$ in Eq.(6) implements such symmetry breaking and introduces the two types of orbitals. Thus, the matrix model of Eq. (6) seems to emerge from some very general arguments.

Let us study the statistics of $`\mathrm{\Delta }_2=ϵ_{N+i}ϵ_{N+i1}`$, where $`ϵ_n`$ is the $`n`$-th (ordered) eigenvalue of $`H`$ and $`i`$ is a small integer. Qualitatively, the picture is similar to that outlined above, i. e., under the increase of $`U`$, the distribution $`P(\mathrm{\Delta }_2)`$ crosses over from the Wigner-Dyson shape to a broader, more symmetric distribution (Fig. 5). The crossover is related to the formation of a gap in the spectrum of $`H`$ (Fig. 6). Eventually, for sufficiently large $`U`$, the Wigner semicirle of the GOE splits into two well separated semicircles, and $`P(\mathrm{\Delta }_2)`$ is controlled by the tail states, where the density of levels is extremely low (by a factor $`N^{1/3}`$ smaller than in the middle of the band (Mehta, 1991)). It is interesting to note that Brezin and Hikami (Brezin and Hikami, 1998) have studied the matrix model in Eq. (6), but with $`G`$ being the Gaussian unitary ensemble, instead of the GOE. They considered the ”critical situation”, when the two bands just touch one another, and showed that a new universal level statistics emerge in the large-$`N`$ limit.

## V Conclusions

A simple model for addition spectra is proposed and studied. It can be viewed as an extension of the charging model. The latter is based on the assumption that the added particles spread over the entire dot, charging it uniformly. This assumption is, clearly, inadequate in the limit of strong interactions, when the charge density in the dot exibits strong spatial changes. The proposed model tries to account for this feature in the most crude phenomenological way, by introducing the parameter $`U`$. This turns out to be sufficient for obtaining a qualitatively correct behavior for the conductance peak spacings. The model, or rather its version in Eq.(6), is also of interest in the random matrix theory. The point is that the standard (orthogonal or unitary) Gaussian ensembles describe well the addition spectra in a disordered dot only in the absence of interactions. The interaction breaks the $`O`$\- or $`U`$-invariance of a Gaussian ensemble. This is seen most clearly in the ”extreme classical limit”, when the electrons become localized on particular site-orbitals. Matrix $`A`$ in Eq.(6) implements such symmetry breaking. There are, of course, other ways to break this symmetry, some of which lead to a crossover between the Wigner- Dyson and the Poisson statistics for level spacings (for a review see Shapiro, 1996). It would be, perhaps, of some interest to investigate the connection between such generalized matrix ensembles and the addition spectra in disordered quantum dots.

The model in Eq. (5) was designed for spinless electrons. For spinfull electrons the extra particles can hop on occupied sites, so that there is no need for the second sublattice. The ”minimal” model is defined on a $`N`$-site lattice, corresponding to the the $`N`$-electron magic configurations, which is the same as for the spinless case. However, each site of the lattice now contains 2 orbitals. They correspond to single and double occupation, respectively, and differ by energy $`U`$. For zero $`U`$ one has a single band, with doubly degenerate levels. The statistics of $`\mathrm{\Delta }_2`$ exibits in this case a bimodal structure, characteristic for non-interacting electrons with spin. For large $`U`$ the picture is essentially the same as for the spinless case: extra particles (holes ) go into the tail of the upper (lower) Hubbard band (Mott, 1990). Since there is no spin degeneracy in these bands (they correspond to different site-occupation nubmers), there is no trace of a bimodal structure. The absence of such structure, as well as other spin effects, have been studied earlier, both for disordered and clean quantum dots (Berkovits, 1998, Hirose and Wingreen, 1999, Baranger et al., 2000, Cha and Yang, 2000, Jacquod and Stone, 2000).

## Acknowledgements

I benefited greatly from discussions with A. Auerbach, J. Feinberg, M. Janssen, A. Kamenev, Y. Levinson, G. Montambaux, J.-L. Pichard, R. Pnini, M. Pollak and U. Sivan. These discussions were instrumental in shaping the final version of the model. I am greatful to D. Gangard and, in particular, to A. Retzker for the help with numerics and figures. The research was partially supported by the DIP project on ”Quantum Electronics in Low-Dimensional Systems” and by the Fund for promotion of research at the Technion.

## References

* Alhassid, Y., 2000, Rev. Mod. Phys., 72, 845.
* Andreev, A.F., and Lifshitz, A.M., 1969, Zh. Eksp. Teor. Fiz. 56, 2057 \[Sov. Phys. JETP 29, 1107\]
* Baranger, H.U., Ullmo, D., and Glasman, L.I., 2000, Phys. Rev. B61, R2425.
* Benenti, G., Waintal, X., Pichard, J.-L., and Shepelyansky, D., 2000, Europ. Phys. Jour. B17, 515.
* Benenti, G., Waintal, X., and Pichard, J.-L., 1999, Phys. Rev. Lett. 83, 1826.
* Berkovits R., 1998, Phys. Rev. Lett. 81, 2128.
* Berkovits, R., and Altshuler, B.L., 1997, Phys. Rev. B55, 5297.
* Blanter, Ya.M., Mirlin, A.D., and Muzykantskii, B.A., 1997, Phys. Rev. Lett. 78, 2449.
* Brezin, E., and Hikami, S., 1998, Phys. Rev. E58, 7176.
* Cha, M.C., and Yang, S.R.E., 2000, Phys. Rev. B61, 1720.
* Chakravarty, S., Kivelson, S., Nayak, C. and Voelker, K., 1999, Phil.Mag. B79, 859.
* Cohen, A., Richter, K. and Berkovits, R., 1999, Phys. Rev. B60, 2536.
* Hirose, K. and Wingreen, N.S., 1999, Phys. Rev. B59, 4604.
* Jacquod, P., and Stone, A.D., 2000, Phys. Rev. Lett. 84, 3938.
* Koulakov, A.A., Pikus, F.G. and Shklovskii, B.I., 1997, Phys. Rev.

B55, 9223.
* Koulakov, A.A. and Shklovskii, B.I., 1998, Phys. Rev. B57, 2352.
* Levit, S. and Orgad, D., 1999, Phys. Rev. B60, 5549.
* Mehta, M.L., 1991, Random Matrices (Academic Press).
* Morris, J.R., Deaven, D.M., and Ho, K.M., 1996, Phys. Rev. B53, R1740.
* Mott, N.F., 1990, Metal-Insulator Transitions (Taylor$`\&`$ Francis).
* Patel, S.R., Cronenwett, S.M., Stewart, D.R., Huibers, A.G., Markus, C.M., Doruoz, C.I., Harris, J.S., Campman, K., and Gossard, A.C., 1998, Phys. Rev. Lett. 80, 4522.
* Prus, O., Auerbach, A., Aloni, Y., Sivan, U., and Berkovits, R., 1996, Phys. Rev. B54, R14289.
* Shapiro, B., 1996, Int. Jour. Mod. Phys. B10, 3539.
* Simmel, F., Abusch-Magder, D., Wharam, D.A., Kastner, M.A., and Kotthaus, J.R., 1999, Phys. Rev. B59, R10441.
* Sivan, U., Berkovits, R., Aloni, Y., Prus, O., Auerbach, A., and Ben-Yoseph, G., 1996, Phys. Rev. Lett. 77, 1123.
* Tsiper, E.U., and Efros, E.L., 1997, J. Phys. Cond. Matter 9, L561.
* Walker, P.N., Montambaux, G., and Gefen, Y., 1999, Phys. Rev. B60, 2541.


text:
# 1 Introduction

## 1 Introduction

The rare decays $`BK\mathrm{}^+\mathrm{}^{}`$ and $`BK^{}\mathrm{}^+\mathrm{}^{}`$, where $`\mathrm{}`$ is either an electron or muon, are highly suppressed in the Standard Model and are expected to occur via electroweak penguin processes. The internal loops in such processes lead to a $`bs`$ transition via an effective flavor-changing neutral current interaction. Standard Model predictions for these decays are dependent to some extent on the modeling of the strong interactions in the hadronic $`BK^{()}`$ transition, but calculations indicate that $`(BK\mathrm{}^+\mathrm{}^{})6\times 10^7`$, while $`(BK^{}\mathrm{}^+\mathrm{}^{})2\times 10^6`$ . These processes provide a possible window into physics beyond the Standard Model, since new, heavy particles such as those predicted by SUSY models can also enter the loops, resulting in significant changes to both the decay rates and kinematic distributions .

Experimentally, the small expected rates make searches for these modes difficult. In addition, backgrounds from tree-level processes, such as $`BJ/\psi K^{()}`$, with $`J/\psi \mathrm{}^+\mathrm{}^{}`$, complicate the analyses. Searches from CDF  and CLEO  have so far yielded only upper limits, although in the case of $`BK^{}\mu ^+\mu ^{}`$ the limit is only about a factor of two above the Standard Model prediction.

In this paper, we report the results of a preliminary analysis to investigate the backgrounds and the ability of the BABAR detector to reject them. We have analyzed an on-resonance data sample of 3.2$`\text{ fb}^1`$, representing about a third of the current BABAR $`\mathrm{{\rm Y}}(4S)`$ integrated luminosity.

As in many other searches for very rare decays, there are important issues that we face in designing event selection criteria in a manner that does not bias the determination of the signal yield, or, at worst, create a false apparent signal in a small statistics event sample. The main goal of our study is to test the performance of a “blind” analysis in which the event selection is optimized without use of the signal or sideband regions in the data. The events in these regions are therefore selected without bias, and the sidebands can be used for background estimation. The decays $`B^+J/\psi K^+`$ and $`B^0J/\psi K^0`$, which have the same topology as the signal, provide a useful control sample to compare data with Monte Carlo simulation. These events can be used to validate the efficiencies predicted by the full Geant321 Monte Carlo simulation. An additional control sample is provided by events with $`K^{()}e^\pm \mu ^{}`$ candidates. Such events monitor the level of combinatorial backgrounds present in the sample. This analysis provides a basis for extending the measurement to the much larger data samples that BABAR will obtain in the future. The 3.2 $`\text{ fb}^1`$ data sample is comparable in size, however, to that used in an earlier CLEO analysis, and we are able to obtain results of comparable sensitivity for the $`B^+K^+\mathrm{}^+\mathrm{}^{}`$ modes.

We analyze four decay modes: $`B^+K^+e^+e^{}`$, $`B^+K^+\mu ^+\mu ^{}`$, $`B^0K^0e^+e^{}`$, and $`B^0K^0\mu ^+\mu ^{}`$. The $`K^0`$ is reconstructed in the $`K^+\pi ^{}`$ final state. In each case, we include the charge conjugate mode as well.

## 2 The BABAR detector and data sample

The BABAR detector is described in detail elsewhere . All components of the detector are used for this study. Of particular importance for this analysis are the charged tracking system and the detectors used for particle identification. At radii between about 3 cm and 14$`\mathrm{cm}`$, charged tracks are measured with high precision in a five-layer silicon vertex tracker (SVT). This device comprises 52 modules built from double-sided AC-coupled silicon strip detectors. Tracking beyond the SVT is provided by the 40 layer drift chamber, which extends from 23.6$`\mathrm{cm}`$ to 80.9$`\mathrm{cm}`$. Just outside the drift chamber is the DIRC, which is a Cherenkov ring-imaging particle identification system. Cherenkov light is produced by charged tracks as they pass through an array of 144 five-meter-long fused silica quartz bars; the light is transmitted to the ends of the bars by total internal reflection, preserving the information on the angle of the light emission with respect to the track direction. Electrons are detected using an electromagnetic calorimeter comprising 6580 thallium-doped CsI crystals. Muons are detected in the Instrumented Flux Return (IFR), in which resistive plate chambers (RPCs) are interleaved with the iron plates of the flux return.

The data used in the analysis were collected with the BABAR detector at the PEP-II storage ring. We analyze a subsample of this data consisting of 3.2$`\text{ fb}^1`$ taken on the $`\mathrm{{\rm Y}}(4S)`$ resonance and 1.2$`\text{ fb}^1`$ at a slightly lower energy to obtain a pure continuum sample. For the purpose of optimizing the event selection, the following simulated event samples were used:

* 4.2 million generic $`\mathrm{{\rm Y}}(4S)B\overline{B}`$ events generated using the Geant321 package,
* 4.2 million continuum $`c\overline{c}`$ events (Geant321),
* 4.8 million continuum $`u\overline{u}`$, $`d\overline{d}`$, and $`s\overline{s}`$ events (Geant321),
* 50$`\text{ fb}^1`$ of a fast, parametrized Monte Carlo simulation with the correct mixture of $`B\overline{B}`$ and continuum events,
* 5,000 events for each signal mode (Geant321).

## 3 Kinematic properties of signal events

The properties of $`\mathrm{{\rm Y}}(4S)B\overline{B}`$ events containing a $`BK^{()}\mathrm{}^+\mathrm{}^{}`$ decay are modeled with a full Geant321 Monte Carlo simulation using the event generator EvtGen . We have implemented a generator based on the matrix elements computed in Ref. . Although there are significant uncertainties due to strong interaction effects, the matrix elements can be rigorously parametrized in terms of form factors, which are Lorentz-invariant functions of $`q^2=M_\mathrm{}^+\mathrm{}^{}^2`$.

In the regions of $`q^2`$ near the $`J/\psi `$ and $`\psi (2S)`$ resonances, interference occurs between the penguin amplitude and that for charmonium production, since the $`J/\psi `$ and $`\psi (2S)`$ can decay into $`\mathrm{}^+\mathrm{}^{}`$. This interference enhances the rate on the low-mass side of the resonance and suppresses it on the high side. These effects are modeled in Ref. , but not in our Monte Carlo generator. However, since $`BJ/\psi K^{()}`$ and $`B\psi (2S)K^{()}`$ are major backgrounds in our analysis, we eliminate most of these mass regions, as described in the following section. Small residual interference effects are expected for masses close to our exclusion regions, but we have ignored these in the present study. Further enhancements in the rate, along with related interference effects, are expected from charmonium resonances with masses above the $`\psi (2S)`$. These effects are relatively small compared to the expected electroweak penguin contribution, and they are ignored in this study.

The Dalitz plot and $`q^2`$ distributions resulting from these matrix elements are shown in Fig. 1. For Dalitz plot variables we have chosen $`q^2`$ and $`E_{\mathrm{lep}}`$, where the lepton energy is given in the $`B`$ meson rest frame, and the positive lepton is selected for the case $`B^+K^+\mathrm{}^+\mathrm{}^{}`$ and $`B^0K^{}\mathrm{}^+\mathrm{}^{}`$ (the negative lepton is selected for the charge-conjugate states). For $`BK\mathrm{}^+\mathrm{}^{}`$, the $`q^2`$ distribution peaks at low values, corresponding to rapid $`K`$ recoil in the $`B`$ rest frame. For $`BK^{}\mathrm{}^+\mathrm{}^{}`$ the $`q^2`$ distribution peaks at large values, corresponding to slow $`K^{}`$ recoil in the $`B`$ rest frame, but there is a pole in the amplitude at $`q^2=0`$, where the photon becomes real. The decay $`BK\mathrm{}^+\mathrm{}^{}`$ does not have a pole at $`q^2=0`$ because $`BK\gamma `$ is forbidden by angular momentum conservation.

## 4 Overview of the analysis method

The analysis method exploits the strong signature and kinematic constraints for the $`BK^{()}\mathrm{}^+\mathrm{}^{}`$ decay and the broad capabilities of the BABAR detector. Lepton and kaon identification, good momentum measurement, and the constraint associated with the known beam energies are essential for the search. Because even low-level backgrounds can be problematic, we have to be concerned about a possible “cocktail” of small contributions. Some of the most important sources are:

* $`BJ/\psi K^{()}`$ or $`B\psi (2S)K^{()}`$ with $`J/\psi `$ or $`\psi (2S)\mathrm{}^+\mathrm{}^{}`$. If one of the leptons radiates a photon, the mass of the dilepton system can fall below the resonance veto region. These events must be rejected with high efficiency, since they have precisely the same topology as signal events.
* Combinatorial backgrounds from random leptons and kaons in $`B\overline{B}`$ and continuum processes ($`e^+e^{}q\overline{q}`$, where $`q=u`$, $`d`$, $`s`$, or $`c`$). In general, the lepton candidates arise from semileptonic decays, although leptons from photon conversions or from hadrons faking a lepton signature also contribute at a low level.
* Certain backgrounds can peak in the signal region, but these are quite rare and easily rejected. An example is $`B^+\overline{D}{}_{}{}^{0}\pi _{}^{+}`$ with $`\overline{D}{}_{}{}^{0}K^+\pi ^{}`$. If both pions are misidentified as muons, this process can be confused with the signal. This type of three-body background is easily rejected with suitable reassignment of particle identification hypotheses and a veto on the $`D`$ mass region.

Our goal in the present search is to gain an understanding of background levels and distributions and to set limits. For this purpose, we assume that each event in the signal region is potentially due to the signal process and do not perform a background subtraction. However, we do carry through a background estimation procedure using the sidebands in data and compare the results with those obtained by applying the procedure to the full Geant321 simulation.

The signal region is defined as a rectangle in the plane defined by two kinematic variables, the beam-energy substituted mass of the $`B`$ candidate $`m_{\mathrm{ES}}`$, and the energy difference $`\mathrm{\Delta }E`$ , where

$`m_{ES}`$ $`=`$ $`\sqrt{\left({\displaystyle \frac{\sqrt{s}}{2}}\right)^2\left({\displaystyle \underset{\alpha =1}{\overset{n}{}}}𝐩_\alpha ^{}\right)^2},`$
$`\mathrm{\Delta }E`$ $`=`$ $`{\displaystyle \underset{\alpha =1}{\overset{n}{}}}\sqrt{m_\alpha ^2+|𝐩_\alpha ^{}|^2}\sqrt{s}/2.`$
The index $`\alpha `$ is over the particles that make up the candidate $`B`$ meson system. $`m_\alpha `$ are the masses of the particles and $`𝐩_\alpha ^{}`$ are their momenta measured in the $`\mathrm{{\rm Y}}(4S)`$ center of mass frame. $`\sqrt{s}`$/2 is one half of the center of mass energy. It is important to distinguish between the laboratory frame, in which the beams have energies $`E_e^{}9.0`$$`\mathrm{GeV}`$ and $`E_{e^+}3.1`$$`\mathrm{GeV}`$, and the center of mass frame, in which the beam energies are equal and the $`\mathrm{{\rm Y}}(4S)`$ is produced at rest. In a symmetric-energy machine, these frames are identical, but at PEP-II the center of mass frame is moving with respect to the laboratory with a boost factor $`\beta \gamma 0.56`$.

## 5 Event selection

We select events with at least five good quality tracks, with 20 or more hits in the drift chamber and originating from the beam spot. Initially, we also require either two loosely identified, oppositely-charged electrons, muons, or an electron-muon pair (for combinatoric background studies). In addition, we require $`R_2<0.8`$, where $`R_2`$ is the ratio of the second and zeroth Fox-Wolfram moments . This requirement provides a first suppression of continuum events, which have a more collimated (“jet-like”) event topology than $`B\overline{B}`$ events. At this stage, $`R_2`$ is evaluated in the center of mass frame using charged tracks only. Events passing these loose event selection criteria are further selected by requiring the electron and muon momenta in the lab frame to satisfy $`p_e>0.5\mathrm{GeV}/c`$ and $`p_\mu >1.0\mathrm{GeV}/c`$. Electrons are identified on the basis of several quantities, primarily $`E/p_e`$, where $`E`$ is the energy measured in the CsI electromagnetic calorimeter and $`p_e`$ is the momentum determined by the tracking system. This ratio is near one for electrons or positrons. Muons are identified primarily by the number of interaction lengths of iron penetrated by the charged track through the magnet flux return, which is instrumented with layers of resistive-plate chambers. The events are then required to lie within a large, rectangular region in the $`\mathrm{\Delta }E`$ vs. $`m_{\mathrm{ES}}`$ plane: $`m_{ES}>5.0`$ $`\mathrm{GeV}/c^2`$and $`|\mathrm{\Delta }E|<0.8`$ $`\mathrm{GeV}`$.

After these basic event selection criteria, we apply a tighter set of particle identification requirements. In addition, electrons and positrons are required to pass the conversions veto, which suppresses the leptons that are likely to have come from photon conversions in the material. The kaon identification is based on combining the energy loss information ($`\mathrm{d}E/\mathrm{d}x`$) from the silicon vertex detector and the drift chamber for momenta $`p`$ $`<`$ 0.6 GeV/c, with the knowledge of the Cherenkov angle from the DIRC for $`p`$ $`>`$ 0.6 GeV/c. For the $`BK\mathrm{}^+\mathrm{}^{}`$ modes we require that the kaon have $`p_K`$ $`>`$ 0.6 GeV/c in the lab frame.

For the $`B^0K^0\mathrm{}^+\mathrm{}^{}`$ channels, we reconstruct the $`K^0`$ in the $`K^+\pi ^{}`$ final state. The kaon candidate is required to be identified as a kaon, while there are no particle identification requirements on the pion candidate. The mass of the $`K^+\pi ^{}`$ pair is required to be within 75$`\mathrm{MeV}/c^2`$ of the $`K^0`$ mass.

The next step is to veto $`BJ/\psi K^{()}`$ and $`B\psi (2S)K^{()}`$. The requirements used for this purpose are somewhat complicated and are shown in Fig. 2. We remove events with dilepton masses consistent with those of the $`J/\psi `$ or $`\psi (2S)`$. These veto regions are shown as pairs of vertical lines in the $`\mathrm{\Delta }E`$ vs. $`M_\mathrm{}^+\mathrm{}^{}`$ plane. However, bremsstrahlung or track mismeasurement can result in a large departure of the dilepton mass from the $`J/\psi `$ or $`\psi (2S)`$ masses. Such departures are also accompanied by a shift in $`\mathrm{\Delta }E`$. We, therefore, remove these events by applying a correlated selection in this plane. Finally, a pair of horizontal lines is shown. There are veto regions above and below these lines, while the signal region lies between them. This region between the horizontal lines corresponds to the requirement that $`\mathrm{\Delta }E`$ be consistent with zero. This requirement is not part of the charmonium veto *per se* but is still useful to see on this figure, because it further restricts the allowed region in this plane. Veto regions for the dimuon final state are defined in a similar way, but the definitions are somewhat different, since there is less bremsstrahlung.

The $`BJ/\psi K`$ modes can also pass this veto if the kaon is misidentified as a lepton (most often a muon). In a similar way $`BD^0\pi `$, where $`D^0K^{}\pi ^+`$, can pass our selection criteria if both of the leptons are fake. These effects can be suppressed by re-assigning the particle masses and excluding mass combinations around the $`J/\psi `$ and the $`D^0`$.

Continuum background is suppressed by using a four-variable Fisher discriminant. The variables are the $`R_2`$, which at this stage is computed using both charged and neutral candidates, the cosine of the angle between the $`B`$ candidate and the $`z`$ axis in the center of mass frame, the cosine of the angle between the thrust axis of the signal and that of the rest of the event in the center of mass frame, and the invariant mass of the kaon and the oppositely charged lepton. The last variable helps discriminate against the $`D`$ semileptonic decays, as it will tend to be below the $`D`$ mass. Both continuum and $`B\overline{B}`$ combinatorial background is further suppressed by requiring that the combined $`B`$ vertex probability be greater than 0.1%.

Finally, we select the signal box in the $`\mathrm{\Delta }E`$ vs. $`m_{ES}`$ plane: 5.272 $`<`$ $`m_{ES}`$ $`<`$ 5.286 GeV (3$`\sigma `$) and $`0.10<\mathrm{\Delta }E<0.06`$ $`\mathrm{GeV}`$ ($`|\mathrm{\Delta }E|<0.06`$ $`\mathrm{GeV}`$) for the electron (muon) channels.

The net effect of our analysis selection criteria on the Dalitz plot is shown in Fig. 3. The requirements to remove $`BJ/\psi K`$ and $`B\psi (2S)K`$ decays are evident as horizontal gaps, while those on the minimum lepton energy in the lab frame result in a loss of acceptance around the edges of the Dalitz plot.

## 6 Physics results

Figure 4 shows a large $`\mathrm{\Delta }E`$ vs. $`m_{ES}`$ region (the “grand sideband”) and a small box indicating the signal region for each of the four modes. These scatter plots show a trend of decreasing background as $`\mathrm{\Delta }E`$ increases from negative to positive values; this is a characteristic of $`B\overline{B}`$ background, where the energy of three or four randomly selected tracks is not as high as the nominal $`B`$ energy in the center of mass frame. Events from the continuum have a more uniform distribution in $`\mathrm{\Delta }E`$ and are more important for positive $`\mathrm{\Delta }E`$ values.

Figure 5 shows the data distribution in $`m_{ES}`$ after applying all other event selection criteria (including the requirement that the events lie within the $`\mathrm{\Delta }E`$ signal region). The distributions are fit to the ARGUS function , whose shape is derived from the large statistics fast Monte Carlo simulation, and are used to estimate the background contribution to the $`m_{ES}`$ signal region. For this study, however, we do not subtract this estimated background from the yield in the signal region; for the purpose of setting the limits we assume that all events in the signal region might be due to the signal events.

Table 1 lists the signal efficiencies, total yield, the expected background, and the 90% C.L. upper limits on the branching fractions. The signal efficiencies were determined from the signal Monte Carlo events. The efficiencies include the branching fractions for the $`K^0`$ modes. The table also lists the total systematic error, determined as described in Section 7. It is important to note that because we decided before the measurement that our result would be expressed as an upper limit, we do not need to apply the Feldman-Cousins procedure (Ref. ), for a two-sided confidence interval, but instead use the original procedure specified in the Particle Data Book , in which zero observed events corresponds to an upper limit of 2.3 events.

## 7 Systematic studies

As we are setting a conservative limit by assuming that all events observed in the signal region are in fact signal, the only systematic effects that affect the limits are those due to the modeling of our signal efficiency and the number of produced $`B`$ mesons in our sample. Table 2 summarizes the systematic uncertainties that we have considered. The uncertainty on the tracking efficiency is 2.5% per track. Lepton and kaon identification efficiencies are measured in control samples in data (e.g., a Bhabha sample for electrons), where a measure of the uncertainty is obtained by comparing different control samples. We assign an error of 2.0% per electron, 2.5% per muon, and 3.0% per kaon. We have also studied the effects of momentum smearing and shifting on the $`\mathrm{\Delta }E`$ selection criteria, and assign an error of about 3.0% for this effect. Comparing data and Monte Carlo efficiencies for the $`m_{ES}`$ requirement in the $`J/\psi `$ control region, we assign a 3.0% error for the modeling of this requirement. We have compared the data and Monte Carlo efficiencies for the $`B`$ vertex probability and the Fisher discriminant selection criteria. For the former, we assign a systematic error of 3.0% for the $`B^+K^+\mathrm{}^+\mathrm{}^{}`$ modes and 4.0% for the $`B^0K^0\mathrm{}^+\mathrm{}^{}`$ modes. For the latter, we assign an error of 3.0% for all of the four modes. Furthermore, as a check that we understand our signal efficiencies, we have compared the yields in data and Monte Carlo events for the $`B^+J/\psi K^+`$ and $`B^0J/\psi K^0`$ (with $`K^0K^+\pi ^{}`$) channels. We have used the world average  branching fractions for $`B^+J/\psi K^+`$ and $`B^0J/\psi K^0`$ and found agreement within the statistical errors between the data and Monte Carlo samples.

## 8 Summary

We have searched for rare $`B`$ decays $`BK^{()}\mathrm{}^+\mathrm{}^{}`$ in a sample of $`3.7\times 10^6`$ $`B\overline{B}`$ events. We find no observable signal for any of the four modes considered, and set preliminary 90% C.L. upper limits on the branching fractions of

$`(B^+K^+e^+e^{})`$ $`<`$ $`12.5\times 10^6,`$
$`(B^+K^+\mu ^+\mu ^{})`$ $`<`$ $`8.3\times 10^6,`$
$`(B^0K^0e^+e^{})`$ $`<`$ $`24.1\times 10^6,`$
$`(B^0K^0\mu ^+\mu ^{})`$ $`<`$ $`24.5\times 10^6.`$
The limits for the $`B^+K^+\mathrm{}^+\mathrm{}^{}`$ modes are comparable to those set by other experiments, while those for $`B^0K^0\mathrm{}^+\mathrm{}^{}`$ are less sensitive with this data sample. We plan to analyze substantially more data in the near future.

## 9 Acknowledgments

The authors thank JoAnne Hewett and Gudrun Hiller for their advice on theoretical issues related to the analysis.

We are grateful for the contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible. We acknowledge support from the Natural Sciences and Engineering Research Council (Canada), Institute of High Energy Physics (China), Commissariat à l’Energie Atomique and Institut National de Physique Nucléaire et de Physique des Particules (France), Bundesministerium für Bildung und Forschung (Germany), Istituto Nazionale di Fisica Nucleare (Italy), The Research Council of Norway, Ministry of Science and Technology of the Russian Federation, Particle Physics and Astronomy Research Council (United Kingdom), the Department of Energy (US), and the National Science Foundation (US). In addition, individual support has been received from the Swiss National Foundation, the A. P. Sloan Foundation, the Research Corporation, and the Alexander von Humboldt Foundation. The visiting groups wish to thank SLAC for the support and kind hospitality extended to them.