text:
# Superconducting Vortex Logic Antidots

\[

## Abstract

We examine a building block for logic devices in which the positions of superconducting vortices in coupled elongated antidots provide the elementary logic states of $`0`$ and $`1`$. We show analytically and through simulation the maximum operating frequency of a pair of antidots as a function of antidot spacing and elongation. At finite temperatures, a signal can propagate through a series of identically shaped antidots with correctly chosen spacing, with an exponential distribution of switching times for the signal to move over by one antidot.

\]

Recently there has been considerable interest in superconductors with artificial pinning arrays, such as artificial arrays of holes in thin-film type-II superconductors . Under a magnetic field, the flux penetrates the film in the form of individual quantized vortices which become pinned by the holes. Experiments  and simulations  have shown that for fields where the number of vortices equals an integer multiple of the number of holes, peaks or anomalies in the transport and critical current occur. These anomalies are correlated with the formation of highly ordered vortex lattice crystals. Effective periodic pinning arrays can also be constructed by placing magnetic dots on superconductors, which produces strong pinning and commensuration effects , or by fabricating large “blind” holes (or antidots) which do not pass completely through the film but instead modulate the film thickness. In this case multiple vortices can be trapped at an individual antidot  and the vortex positions from one antidot are correlated with the positions in neighboring antidots . Recently experiments have also demonstrated that nanodots of superconducting materials can capture individual vortices . Besides circular pinning sites, it is also possible to fabricate elongated pinning sites . In this case an individual vortex may not necessarily be located at the center of the pinning site but will reside along a line that cuts through the center of the site in the long direction.

Since the vortices have specific arrangements at certain integer matching or fractional matching fields, it should be possible to use the locations of the vortices or the flux configurations as elementary logic states such as 0 and 1. These states can then be measured or propagated by the application of a magnetic field or current. Puig et al.  made one initial proposal along these lines, in which a superconducting island with $`2\times 2`$ plaquettes was considered as a logic element. In this case the minima in the resistance can be associated with different flux configurations.

The quantum-dot cellular automata (QCA) is another system that uses the locations or configurations of the particles as basic logic states . A basic cell consists of four quantum dots containing two localized electrons. Due to the Coulomb repulsion between the electrons, there are two possible ground state configurations with the electrons located at the diagonals of the cell, slanted toward the right or left. These two states provide the fundamental logic units. With different geometrical arrangements of the basic cell, various logic devices can be constructed. In the QCA system the logic states are propagated by quantum-mechanical means. So far, elementary QCA systems have been demonstrated to operate only at extremely low temperatures. Since these systems must be adiabatically switched between logic states, signal propagation times are relatively slow.

In this paper we examine the building blocks for a superconducting vortex logic system  proposed in analogy with the quantum cellular automata. In our model we consider a superconductor that has been nanofabricated to contain elongated pinning sites, where each site captures a single vortex. Our basic unit consists of coupled elongated antidots of identical shape where, due to the repulsion of the vortices, two ground states can be formed. We show that by flipping one vortex in a pair of antidots, the vortex in the other antidot also flips, and we examine the response time of the second vortex as a function of the antidot geometry. We show that at finite temperatures, a signal can be propagated through a series of identically shaped antidots that have been fabricated in a specific geometry.

The individual logic element in our system is a single vortex inside an elongated antidot, a pinning site created through nanolithography techniques. As illustrated in Fig. 1(a), we take the elongation to be in the $`y`$ direction. In the presence of a neighboring antidot, parallel to the first antidot and offset in the $`x`$ direction, the two equilibrium positions of the vortex are at either end of the antidot. Thus, to create our logic element, we define a vortex position at the top of the antidot to be a logic ’1’, as shown in Fig. 1(a), while a vortex positioned at the bottom of the antidot is in a logic ’0’ state. One way in which the logic state of the antidot can be externally controlled is by means of an STM tip made of magnetic material. If this tip is moved in the $`y`$ direction at a given frequency, it should be possible to drag the vortex in the antidot below between the logic ’0’ and ’1’ states. The position of a vortex in a given antidot can also be detected by using an STM in spectroscopy mode. The vortex location can be identified due to the fact that the spectrum at the vortex core is very different from both superconducting or normal spectra .

In order to transmit information through the system, neighboring antidots must be coupled. We first consider an isolated pair of antidots a distance $`a`$ apart, illustrated in Fig. 1(b). The vortices are prepared in an equilibrium configuration where vortex A is in logic state 1, and vortex B is in logic state 0. At time $`t=t_0`$ the position of vortex A is flipped by external means to be in logic state 0. For operation of the vortex logic, we require that vortex B will subsequently flip without external intervention to logic state 1, so that the final state of the system is reversed from the initial state.

In the case of two isolated antidots, vortex B does not need to overcome an energy barrier in order to flip. We can estimate the operating frequency of the two-vortex system by finding the time required to complete the flip. The interaction between the two vortices in a thin film is given by the Pearl interaction , which can be written

$$f^{vv}(r)=\frac{\mathrm{\Phi }_0^2}{2\mu _0\pi \lambda ^2}\frac{d}{r}$$
(1)
in the limit $`r2\lambda ^2/d`$. Here, $`\mathrm{\Phi }_0`$ is the elementary flux quantum, $`\mu _0`$ is the permeability of free space, $`\lambda `$ is the London penetration depth of the superconductor, and $`d`$ is the thickness of the superconducting film. The vortices obey overdamped dynamics given by $`𝐟_i=\eta 𝐯_i`$, where $`\eta =B_{c2}\mathrm{\Phi }_0/\rho _N`$, $`B_{c2}`$ is the upper critical field, and $`\rho _N`$ is the normal state resistivity. We assume that the antidots are spaced $`a`$ apart in the $`x`$ direction, and are of length $`\alpha `$ in the elongated $`y`$ direction. If the position of vortex A is switched and held fixed such that both vortices are in the same logic state, vortex B requires a transit time $`t_{tr}`$ before it reaches the opposite side of the antidot and the system returns to equilibrium. The transit time can be written

$$t_{tr}=_\delta ^\alpha \frac{1}{v(y)}𝑑y$$
(2)
where $`v(y)=𝐟\widehat{y}/\eta `$, the $`y`$ velocity of vortex B at position $`y`$. The integration must start from a small offset $`\delta `$ because if the two vortices are at the same $`y`$ location, they exert no force on each other in the $`y`$ direction and will not move without thermal assistance. This means that in an experiment conducted at low temperatures where thermal fluctuations are insignificant, vortex A must be moved past the $`y`$ position of vortex B by a distance $`\delta `$. Putting in the approximation for the thin film interaction gives

$$\frac{t_{tr}}{t_0}=\frac{1}{2}(\alpha ^2\delta ^2)+a^2\mathrm{ln}\left(\frac{\alpha }{\delta }\right)$$
(3)
where time is measured in units of $`t_0=\eta /f_0^{}`$, with $`f_0^{}=\mathrm{\Phi }_0^2/(2\pi \mu _0\lambda ^2)`$. Distances are measured in units of the film thickness $`d`$. Thus decreasing the spacing $`a`$ or the antidot length $`\alpha `$ will produce faster switching.

We compare this theoretical result to the transit times of the vortex obtained from a two-dimensional numerical simulation with open boundary conditions. The equation of motion for a vortex $`i`$ is

$$𝐟_i=\eta 𝐯_i=𝐟_i^{vv}+𝐟_i^p+𝐟_i^T$$
(4)
The Langevin force from the temperature is $`𝐟_i^T`$ and has the properties $`<f^T(t)>=0`$ and $`<f_i^T(t)f_j^T(t^{^{}})>=2\eta k_BT\delta _{ij}\delta (tt^{^{}})`$. Initially we consider the case $`T=0`$. The force $`𝐟_i^p`$ is from the pinning well, represented by an ordinary parabolic trap that has been split in half and elongated in the $`y`$ direction. There is no $`y`$ direction confining force in the central elongated portion of the pin. The maximum pinning force is $`f_p`$ in both the $`x`$ and $`y`$ directions. The radius of the pin is $`r_p`$, and the total length of the pin in the $`y`$ direction is $`\alpha +2r_p`$. We first consider two wells a distance $`a`$ apart, each containing a single vortex.

Figure 2(a) shows the transit time $`t_{tr}/t_0`$ obtained from simulation (symbols) as a function of $`a`$ with $`T=0`$, $`r_p=0.24\lambda `$, $`f_p=0.4f_0^{}`$, $`\delta =0.24\lambda `$, and $`\alpha /\lambda =3`$, 4, 5, and 6, along with corresponding plots of Eq. 3 (solid lines). In Fig. 2(b) we show $`t_{tr}/t_0`$ as a function of $`\alpha `$ for $`a/\lambda =3`$, 4, 5, and 6. In each case we find excellent agreement between simulation and theory. In the case of a Nb film of thickness 2000 Å with antidots of anisotropy $`\alpha =3\lambda =`$135 nm, and with spacing between the dots $`a=3\lambda =`$135 nm, the transit time of $`t_{tr}=27.2t_0`$ from Eq. 3 with $`\delta =0.24\lambda =`$10.8 nm corresponds to an actual time of 1.4 ns, indicating that the maximum operating frequency for Nb antidots of this size and geometry is 696 MHz. Smaller or more closely spaced dots will operate at higher frequencies.

To create any type of device, the logic states need to be propagated over distances further than a single well. We therefore consider the case of three antidots, A, B, and C, illustrated in Fig. 1(c). All of the antidots are the same shape and size, but the antidot spacing varies, so that antidots A and B are separated by $`a`$, but antidots B and C are separated by $`a^{}`$. If $`a^{}=a`$, then when the vortex in dot A is switched externally, vortex B will not be able to switch because it experiences a potential barrier due to the presence of vortex C. Instead, the new minimum energy for vortex B will be at the center of well B. This is undesirable for logic operations. In order to allow vortex B to switch, we must increase $`a^{}>a`$ to reduce the energy barrier produced by the repulsion from vortex C, and we must also introduce thermal fluctuations to allow thermal activation over the energy barrier. In a low-temperature material such as Nb, the thermal fluctuations may remain prohibitively small below $`T_c`$ unless very small antidots are fabricated. In a high-temperature superconductor such as BSCCO, it is much easier to produce thermal activation even for relatively large antidots.

For the configuration shown in Fig. 1(c), where vortex A has been switched to a new logic state but vortex B has not yet switched, we can write an expression for the $`y`$ position of vortex B at which the net $`y`$ force on vortex B from vortices A and C is zero. We rescale all distances by $`\alpha `$, so that $`\stackrel{~}{y}=y/\alpha `$, $`\stackrel{~}{a}=a/\alpha `$, and $`\stackrel{~}{a}^{}=a^{}/\alpha `$, and we take $`\stackrel{~}{y}=0`$ as the starting position of vortex B and $`\stackrel{~}{y}=1`$ as the final switched position. We obtain

$$2\stackrel{~}{y}_B^33\stackrel{~}{y}_B^2+(1+\stackrel{~}{a}^2+\stackrel{~}{a}^2)\stackrel{~}{y}_B\stackrel{~}{a}^2=0.$$
(5)
After vortex B reaches $`\stackrel{~}{y}_B`$ satisfying this expression, it has crossed the potential barrier and can move freely to the opposite side of the well. Similarly, after vortex B has switched to the new logic state \[Fig. 1(d)\], we can write an expression for the $`y`$ position at which the force on vortex C from vortices A and B in the switched state is zero:

$$2\stackrel{~}{y}_C^33\stackrel{~}{y}_C^2+(1+\stackrel{~}{a}^2+2\stackrel{~}{a}\stackrel{~}{a}^{}+2\stackrel{~}{a}^2)\stackrel{~}{y}_C\stackrel{~}{a}^2=0.$$
(6)
The position $`\stackrel{~}{y}_C`$ marks the end of the potential energy barrier.

Vortices B and C can cross their respective energy barriers by thermal activation. Well C should be placed as close as possible to well B in order to enhance the coupling of vortex C to vortex B, so $`a^{}`$ should be made as small as possible. However, as $`a^{}`$ approaches $`a`$, the coupling between vortices A and B is weakened, and vortex B will never switch initially. Thus $`a^{}`$ should be chosen just slightly larger than the value at which vortex B can no longer switch. For $`\stackrel{~}{a}=0.187`$, vortex B ceases to switch below $`\stackrel{~}{a}_{crit}^{}=0.29`$; thus, for optimal operation, $`\stackrel{~}{a}^{}`$ should be chosen slightly above $`\stackrel{~}{a}_{crit}^{}`$.

The thermal fluctuations must be large enough for vortices B and C to reach positions $`\stackrel{~}{y}_B`$ and $`\stackrel{~}{y}_C`$ at $`\stackrel{~}{a}_{crit}^{}`$ by thermal activation in order for the signal to propagate. The vortices also thermally fluctuate in the $`x`$ direction, and therefore if the antidots are wide in the $`x`$ direction, the value of $`a_{crit}^{}`$ is increased compared to the value obtained from Eqn. 5 since the minimum possible $`x`$ distance $`a_{\mathrm{eff}}^{}`$ between vortices B and C is less than $`a^{}`$. The following relation holds: $`a^{}2r_p<a_{\mathrm{eff}}^{}<a^{}`$. The same is true of $`a`$. Thus in experiments the pins must be fabricated slightly further than $`a`$ or $`a^{}`$ apart. Additionally, since the switching of the vortices is now thermally activated, strict clocking of the signal is no longer possible as it was in the case of two wells. A mechanism to obtain strict clocking has been demonstrated in Ref. .

The distance between the wells required for the operation of the three well system scales with $`\alpha `$. The smaller $`\alpha `$ is, the smaller the distance between pins can be made and still allow operation of all three wells. We consider the case of BSCCO where thermal activation can play a significant role below $`T_c`$. In Fig. 3 we show the distribution of switching times $`P(t_s)`$ obtained from 200 runs with different random temperature seeds. Here $`f_T=1.2f_0^{}`$, $`a=0.5\lambda `$, $`a^{}=0.68\lambda `$, and $`\alpha =2.48`$. Fig. 3(a) shows that $`P(t_s^{(B)})`$ for vortex B is exponentially distributed with mean value $`\overline{t}_s^{(B)}=1078t_0`$. $`P(t_s^{(C)})`$ for vortex C, with time measured from $`t=0`$ before vortex B has switched, is plotted in Fig. 3(b), and the distribution is clearly broader and more heavily weighted toward later times. $`P(t_s^{(C)})`$ for vortex C is merely the product of two exponential distributions, as can be seen from Fig. 3(c), where we plot the switching time of vortex C with time measured from $`t=t_s^{(B)}`$, the switching time of vortex B. $`P(t_s^{(C)}t_s^{(B)})`$ is also exponentially distributed with a mean value of 693$`t_0`$.

More than three wells can be connected, but the spacing between well $`n`$ and $`n+1`$ must always be greater than the spacing between wells $`n`$ and $`n1`$. Thus there is a practical limitation on the total length of device that can be fabricated in this fashion; when the wells are spaced too far apart, the vortices will thermally decouple. Additionally, as illustrated in Fig. 3, the distribution of switching times for the final well will become increasingly broad and approach a Gaussian as the number of wells is increased. Regardless of the number of wells, these logic elements operate only from the narrowly spaced end to the widely spaced end. A signal introduced at the widely spaced end will not be able to propagate to the narrow end. If all the wells are spaced equally, it is still possible for excitations to propagate through the wells under thermal activation if the temperature is high enough. However, these excitations move diffusively in either direction and are not well controlled. It is possible instead to employ a ratchet mechanism to obtain controlled, clocked motion of the signal .

In summary we have examined the basic building blocks for a vortex cellular automata that is analogous to a quantum-dot cellular automata. We consider elongated pinning sites in superconducting samples where there is one vortex per antidot. The vortices can form two ground states with the vortices located at diagonals in order to minimize their interaction energy. We obtained analytically and in simulation the maximum frequency of operation for a two well system. For the simplest pipeline geometry of three wells we find that, for finite temperatures, a change in logic state introduced at the first well can be propagated over specified distances.

Acknowledgments— We thank W. Kwok and V.V. Metlushko for useful discussions. This work was supported by the U.S. DoE, Office of Science, under Contract No. W-31-109-ENG-38. CJOR and CR were supported by the U.S. DoE under Contract No. W-7405-ENG-36. BJ was supported by NSF-NIRT award DMR02-10519. and the Alfred P. Sloan Foundation.


format: text

text:
# Measurement of the branching ratio for the decay 𝐾^±→𝜋^±𝜋⁰𝜋⁰ with the KLOE detector

## Abstract

We have measured the absolute branching ratio BR($`K^\pm \pi ^\pm \pi ^0\pi ^0`$) with the KLOE detector at the DAΦNE $`e^+e^{}`$ collider. We collected $``$$`3.3\times 10^7`$ tagged charged kaons, from the reaction $`e^+e^{}\varphi K^+K^{}`$. We find BR($`K^\pm \pi ^\pm \pi ^0\pi ^0`$)=$`(1.763\pm 0.013_{\mathrm{stat}}\pm 0.022_{\mathrm{syst}})\times 10^2`$.

, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,2, , , , , , , , , , , , , , , , , , ,1, , , 11footnotetext: Corresponding author: Andrea Ventura, e-mail ventura@le.infn.it 22footnotetext: Corresponding author: Margherita Primavera, e-mail primaver@le.infn.it Recently there has been renewed interest in three pion decays of charged kaons . Because of the small energy available in the reaction, $`K3\pi `$ is an ideal process where to apply the notion of the Goldstone-boson nature of the pseudoscalar mesons, by testing the predictions obtained from the chiral lagrangian realization of the $`\mathrm{\Delta }S=1`$ weak interactions .

The branching ratio for $`K^\pm \pi ^\pm \pi ^0\pi ^0`$ ($`\tau ^{}`$$`^3`$33footnotetext: In the following text, this old notation will be often used to refer to $`K^\pm \pi ^\pm \pi ^0\pi ^0`$ decay.) decays is the least well known of the hadronic $`K^\pm `$ decay modes. The most accurate measurement to date of BR($`\tau ^{}`$), performed with a sample of $``$1300 $`K^+\pi ^+\pi ^0\pi ^0`$ decays in flight , has an accuracy of 3.3% and dates back to more than thirty years ago. We present a new determination of BR($`\tau ^{}`$), performed with the KLOE detector  at the Frascati $`\varphi `$–factory DAΦNE . DAΦNE is an $`e^+e^{}`$ collider operated at a CM energy $`W=m_\varphi `$$``$1020 MeV/$`c^2`$. About 50% of the $`\varphi `$ mesons decay to $`K^+K^{}`$ pairs. Detection of a $`K^\pm `$ meson tags the presence of a $`K^{}`$. Samples of $`N_K`$ (positive and negative) tagged kaons can be identified and used to search for the $`\tau ^{}`$ decay. The branching ratio is then given by BR=$`N(\tau ^{})/N_K`$. This procedure allows to determine directly the absolute branching ratio of interest. The one necessary condition is that the trigger of the tagging kaon should be only slightly dependent on the decay mode of the tagged kaon. We ensure this to be the case by verifying that the tagging kaon did in fact trigger by itself, using the complete set of information, recorded for each event, concerning all signals produced and processed by the trigger system . The measurement described in the following is based on data collected at the $`\varphi `$ peak in 2001 and 2002, corresponding to an integrated luminosity of 441 pb-1 or, equivalently, to the production of $``$$`1.5\times 10^9`$ $`\varphi `$ mesons.

The KLOE detector consists of a large cylindrical drift chamber surrounded by a lead-scintillating fiber sampling calorimeter. A superconducting coil outside the calorimeter provides a 0.52 T magnetic field parallel to the beam axis. The drift chamber  is 4 m in diameter and 3.3 m long. It uses aluminium field wires, a 90% helium–10% isobutane gas mixture and is constructed entirely with carbon fiber-epoxy composites. Transparency to photons is thus maximized and multiple Coulomb scattering minimized. Single point resolutions are $`150\mu \mathrm{m}`$ in the transverse plane and $`2\mathrm{mm}`$ longitudinally. Typical momentum resolution is $`\sigma (p_T)/p_T0.4\%`$. Vertices are reconstructed with a spatial resolution of $`3\mathrm{mm}`$.

The calorimeter , divided into a barrel and two endcaps, covers 98% of the solid angle. The readout granularity is $``$4.4$`\times `$4.4 cm2 with 2440 “cells” arranged in five layers. Cells with signals close in time and space are grouped into a “calorimeter cluster”. For each cluster, the energy deposit is the sum of the single cell energies. Arrival time and position are calculated as energy-weighted averages over the fired cells. The energy resolution is $`\sigma _E/E=5.7\%/\sqrt{E(\mathrm{GeV})}`$ and the time resolution $`\sigma _t=54\mathrm{ps}/\sqrt{E(\mathrm{GeV})}50\mathrm{ps}.`$ The KLOE trigger  requires two local energy deposits in the calorimeter above threshold ($`50\mathrm{MeV}`$ in the barrel, 150 MeV in the endcaps) or an appropriate hit multiplicity in the chamber. The trigger system also produces a cosmic-ray veto using the signals from the outer calorimeter layers.

Charged kaons are identified by observation of their two-body decays: $`\mu ^\pm \nu _\mu (K_{\mu 2}^\pm )`$ and $`\pi ^\pm \pi ^0(K_{\pi 2}^\pm )`$, comprising $``$85% of all $`K^\pm `$ decays. The monochromatic charged particle momentum in the kaon CM produces a very clean signature and is exploited to identify $`K^\pm `$ production and thus tag the presence of another charged kaon, whose decay can be investigated.

As mentioned above, since all information on how the trigger was formed is available, we require that the tagging kaon indeed satisfied the trigger requirements. Charged kaon decays are found by looking for a charged particle originating from the interaction point and ending in a decay vertex reconstructed inside the drift chamber. We require that: 1) the radial distance of the decay vertex from the beam axis be between 40 and 150 cm; 2) the kaon momentum at the decay point be between 70 and 130 MeV/$`c`$ and the point of closest approach, $`x_c,y_c,z_c`$, of its track to the beam axis satisfy $`\sqrt{x_c^2+y_c^2}<10`$ cm and $`|z_c|<20`$ cm (the origin of the reference frame is the collider interaction point); 3) the momentum of the decay product $`p_d`$, assumed to be a pion, satisfy $`180<p_d<270`$ MeV/$`c`$ in the kaon CM. In Fig. 1 the $`p_d`$ distribution is shown for data and Monte Carlo: the two well distinct peaks correspond to $`K_{\pi 2}^\pm `$ and $`K_{\mu 2}^\pm `$ decays.

$`K_{\pi 2}^\pm `$ decay identification requires a decay particle momentum in the $`K`$ meson rest frame compatible with the expected value of 205.14 MeV/$`c`$ and detection of a $`\pi ^0`$. The latter requires observation of two energy clusters in the calorimeter, photons, with energies $`>`$15 MeV and times of flight consistent with the decay vertex, which we refer to as “on-time” clusters in the following. In addition, we require $`85<m_{\gamma \gamma }<185`$ MeV/$`c^2`$. $`K_{\mu 2}^\pm `$ decay identification requires the decay particle momentum in the $`K`$ meson rest frame to be compatible with 235.53 MeV/$`c`$. Furthermore, the missing mass $`m_m`$ at the decay vertex must satisfy $`|m_m^2|<5000`$ MeV2/$`c^4`$, which removes residual background from $`\pi ^\pm \pi ^0`$ decays. The main background in the two samples is mostly due to $`K_{e3}^\pm `$ and $`K_{\mu 3}^\pm `$ decays, especially for $`K_{\pi 2}^\pm `$ decays. From Monte Carlo simulation we find: $`f_{\mathrm{bckgnd}}(\pi ^\pm \pi ^0)=(0.37\pm 0.02)\%,f_{\mathrm{bckgnd}}(\mu ^\pm \nu )=(0.21\pm 0.02)\%`$, where the quoted errors include statistical and systematic contributions.

$`\tau ^{}`$ decays are events in which a decay vertex is observed for the tagged kaon, as for the tagging kaon above. In addition, at least four on-time energy deposits with $`E>15`$ MeV must be detected. The maximum charged pion momentum in the $`K^\pm `$ frame is $`132.95`$ MeV/$`c`$ for $`\tau ^{}`$ decays. We accept only pions with $`p_{\pi ^\pm }<135`$ MeV/$`c`$. The on-time condition ensures that the arrival times are consistent with photons originating at the same point, the decay vertex. Since the cluster time $`t_i`$ is well known but the decay instant is not, we construct the differences $`\mathrm{\Delta }_{jk}`$ between the quantities $`t_ir_i/c`$ for all cluster pairs $`j,k`$, with $`j,k=1,\mathrm{},N`$, $`k>j`$ and $`r_i`$ the distance from the centroid of the cluster $`i`$ to the decay vertex. Each difference is normalized using the appropriate error so that its variance is unity:

$$\mathrm{\Delta }_{jk}=\frac{t_jt_k(r_jr_k)/c}{\sqrt{\sigma _{t,j}^2+\sigma _{t,k}^2}}.$$

We require $`|\mathrm{\Delta }_{jk}|<4`$ for all pairs. Accidental on-time background clusters are effectively rejected by this cut. Events with more than four on-time clusters in the final sample, apart from the decay $`K^\pm \pi ^\pm \pi ^0\pi ^0\gamma `$ (BR$`<10^5`$), are due to residual on-time background and photon showers reconstructed as multiple clusters. We finally require that $`E_i<450`$ MeV, for the four most energetic energy deposits. Data from 188 pb-1 are used to search for $`\tau ^{}`$ decays. The remainder of the data is divided into three subsets for the purposes of efficiency evaluation. We find 41896 $`K^+`$ and 41155 $`K^{}`$. We get 52253 and 30798 $`\tau ^{}`$ decays which pass all requirements above, tagged by $`1.9925\times 10^7`$ $`K_{\mu 2}^\pm `$ and $`1.2753\times 10^7`$ $`K_{\pi 2}^\pm `$ decays, respectively.

As background in the $`K^\pm \pi ^\pm \pi ^0\pi ^0`$ sample, we have considered $`K_{\pi 2}^\pm `$, $`K_{e3}^\pm `$ and radiative decays like $`K^\pm \pi ^\pm \pi ^0\gamma `$, in which a “spurious” cluster in the calorimeter has been paired with the cluster of the radiated $`\gamma `$ (the probability for this to happen has been estimated from a $`K_{\pi 2}^\pm `$ data sample to be $`8\%`$). The nuclear interaction of a charged kaon (mainly $`K^{}`$) with the beam pipe and the chamber walls can produce secondaries that simulate the $`\tau ^{}`$ signal. Finally, we include $`K^\pm l^\pm \pi ^0\pi ^0\nu _l`$ ($`K_{l4}^{}`$) decays. The relative contributions of these backgrounds have been estimated by Monte Carlo simulation, while the total amount of background events has been determined by fitting the observed missing mass ($`m_m`$) spectrum at the decay vertex (assumed to be $`K\pi +x`$) with the sum of the simulated signal and background spectra. The background fraction estimated from the Monte Carlo simulation had to be corrected by a factor $`1.13\pm 0.07`$. After applying this correction we find that the average total background fraction is $`f_{\mathrm{bckgnd}}(\tau ^{})=(0.75\pm 0.11)\%`$, where the quoted error includes the contributions from the finite statistics of the Monte Carlo spectra and from the fit. Contamination from other decays not originating from $`K^\pm `$ is negligible.

The overall $`\tau ^{}`$ efficiency $`ϵ(\tau ^{})`$ is given by the product $`ϵ=ϵ_Kϵ_vϵ_\gamma `$. $`ϵ_K`$ is the efficiency to reconstruct the track of the charged kaon that undergoes $`\tau ^{}`$ decay. $`ϵ_v`$ includes the efficiency for decay vertex finding and losses due to the pion momentum cut. $`ϵ_\gamma `$ is the efficiency for observing at least four photons in the calorimeter satisfying the energy and timing requirements above. Each factor has been measured separately for positive and negative kaon charge and for $`K_{\mu 2}^\pm `$ and $`K_{\pi 2}^\pm `$ tags by using special control samples  extracted from data taken in run periods close in time but different from the one used for signal selection. Control samples selected using calorimeter variables have been used to compute the efficiencies involving the drift chamber and vice-versa.

Where necessary, the efficiencies have been corrected for contamination of the control samples, as estimated by Monte Carlo simulation. For all samples the purities are $`>95\%`$. The efficiency $`ϵ_K`$ has been found, within errors, to be the same for the two tags and kaon charges and its average value is equal to $`0.466\pm 0.001_{stat}\pm 0.002_{syst}`$. It is dominated by the geometrical acceptance (more than one third of the charged kaons decay before reaching the drift chamber). The effect of the above mentioned nuclear interactions of $`K^\pm `$ has been taken into account by properly weighting, to get $`N_K`$, the number of tagging kaons with a factor containing the probability for charged kaons to interact with the materials in front of the chamber. These weights have been extracted from data by using $`K_{\pi 2}^\pm `$ and $`K_{\mu 2}^\pm `$ samples: $`w^+=0.9696\pm 0.0036`$ and $`w^{}=0.9970\pm 0.0034`$ respectively for positive and negative tags.
The efficiency for reconstructing the $`K^\pm `$-$`\pi ^\pm `$ vertex, $`ϵ_{vtx}`$, has been parameterized as a function of the charged pion momentum, as shown in Fig. 2 (left). After including the effect of the cut on the pion momentum, the average value $`ϵ_v=0.539\pm 0.003_{stat}\pm 0.002_{syst}`$ is obtained for the signal. The systematic error includes small differences between the two tags.
The cluster efficiency is evaluated from a control sample and corrected by using a Monte Carlo simulation to consider the effects of accidental clusters. Fig. 2 (right) shows the single cluster efficiency $`ϵ_{clu}`$ vs energy. The only appreciable source of the difference in $`ϵ_\gamma `$ between the $`K_{\mu 2}^\pm `$ and $`K_{\pi 2}^\pm `$ tags results from the on-time cluster requirement, due to the different cluster multiplicities in each event type. As average value on positive and negative kaon charge we find: $`ϵ_\gamma =0.645\pm 0.001_{stat}\pm 0.004_{syst}(0.625\pm 0.002_{stat}\pm 0.005_{syst})`$ in the case of $`K_{\mu 2}^\pm `$ ($`K_{\pi 2}^\pm `$) tag.
Finally, $`ϵ(\tau ^{})=0.1621\pm 0.0010_{stat}\pm 0.0016_{syst}(0.1573\pm 0.0011_{stat}\pm 0.0018_{syst})`$ combining positive and negative kaons for $`K_{\mu 2}^\pm `$ ($`K_{\pi 2}^\pm `$) tags.

The residual dependence of our tagging procedure on the tagged kaon decay mode is expressed by the ratio $`R_{tag}`$ of the average (over all channels) tagging efficiencies (including trigger) and the values of the same efficiencies in the case of $`\tau ^{}`$ signal. $`R_{tag}`$ is determined from Monte Carlo and the averages on the two kaon charges are: $`R_{tag}=1.058\pm 0.005(1.099\pm 0.009)`$ for $`K_{\mu 2}^\pm `$ ($`K_{\pi 2}^\pm `$) tags, where the error includes systematic effects due to small differences between data and Monte Carlo, as estimated on a set of variables to which $`R_{tag}`$ is found to be sensitive.

The $`\tau ^{}`$ branching ratio is computed separately for $`K_{\mu 2}^\pm `$ and $`K_{\pi 2}^\pm `$ tags and for the two signs of the kaon charge, from:

$$\mathrm{BR}(K^\pm \pi ^\pm \pi ^0\pi ^0)=\frac{N(\tau ^{})}{N_K}\times \frac{1f_{\mathrm{bckgnd}}(\tau ^{})}{1f_{\mathrm{bckgnd}}(\mathrm{tag})}\times \frac{R_{tag}}{ϵ(\tau ^{})}\times \frac{1}{[BR(\pi ^0\gamma \gamma )]^2}$$

where BR$`(\pi ^0\gamma \gamma )`$=0.98798$`\pm `$0.00032 . We find:

$$\mathrm{BR}(K^+\pi ^+\pi ^0\pi ^0)=(1.760\pm 0.024_{\mathrm{stat}}\pm 0.017_{\mathrm{syst}})\times 10^2(K_{\mu 2}^{}tag)$$

$$\mathrm{BR}(K^{}\pi ^{}\pi ^0\pi ^0)=(1.793\pm 0.024_{\mathrm{stat}}\pm 0.026_{\mathrm{syst}})\times 10^2(K_{\mu 2}^+tag)$$

$$\mathrm{BR}(K^+\pi ^+\pi ^0\pi ^0)=(1.769\pm 0.027_{\mathrm{stat}}\pm 0.020_{\mathrm{syst}})\times 10^2(K_{\pi 2}^{}tag)$$

$$\mathrm{BR}(K^{}\pi ^{}\pi ^0\pi ^0)=(1.724\pm 0.027_{\mathrm{stat}}\pm 0.027_{\mathrm{syst}})\times 10^2(K_{\pi 2}^+tag)$$

yielding the average value:

$$\mathrm{BR}(K^\pm \pi ^\pm \pi ^0\pi ^0)=(1.763\pm 0.013_{\mathrm{stat}}\pm 0.022_{\mathrm{syst}})\times 10^2$$

where the final systematic error takes into account correlations. All contributions to the uncertainty are listed in Table 1. The statistical error is dominated by the uncertainty in determining the efficiencies, due to the limited statistics of the data samples used in their evaluations. The stability of the result has been verified using data sets from different running periods. It should be pointed out that this measurement is fully inclusive of the radiative decays $`K^\pm \pi ^\pm \pi ^0\pi ^0\gamma `$ (there is no cut on the radiated photon energy $`E_\gamma `$), since we require at least four (rather than exactly four) on-time clusters in the signal selection. An upper limit to the contribution of this decay is estimated to be 0.1$`\%`$ of the selected events, by taking into account possible differences in the efficiency for the two channels.

We thank the DA$`\mathrm{\Phi }`$NE team for their efforts in maintaining low background running conditions and their collaboration during all data-taking. We want to thank our technical staff: G.F.Fortugno for his dedicated work to ensure an efficient operation of the KLOE Computing Center; M.Anelli for his continous support to the gas system and the safety of the detector; A.Balla, M.Gatta, G.Corradi and G.Papalino for the maintenance of the electronics; M.Santoni, G.Paoluzzi and R.Rosellini for the general support to the detector; C.Pinto (Bari), C.Pinto (Lecce), C.Piscitelli and A.Rossi for their help during major maintenance periods. This work was supported in part by DOE grant DE-FG-02-97ER41027; by EURODAPHNE, contract FMRX-CT98-0169; by the German Federal Ministry of Education and Research (BMBF) contract 06-KA-957; by Graduiertenkolleg ‘H.E. Phys. and Part. Astrophys.’ of Deutsche Forschungsgemeinschaft, Contract No. GK 742; by INTAS, contracts 96-624, 99-37; and by TARI, contract HPRI-CT-1999-00088.


text:
# Quantum fluctuations of a vortex in an optical lattice

## Abstract

Using a variational ansatz for the wave function of the Bose-Einstein condensate, we develop a quantum theory of vortices and quadrupole modes in a one-dimensional optical lattice. We study the coupling between the quadrupole modes and Kelvin modes, which turns out to be formally analogous to the theory of parametric processes in quantum optics. This leads to the possibility of squeezing vortices. We solve the quantum multimode problem for the Kelvin modes and quadrupole modes numerically and find properties that cannot be explained with a simple linear-response theory.

Introduction.— There is a large body of research on vortices in Bose-Einstein condensates, both before and after their first experimental observation Matthews et al. (1999); Madison et al. (2000). A recent review of this vibrant field was given by Fetter and Svidzinsky Fetter and Svidzinsky (2001). One reason for this interest is that vortices play a crucial role in explaining rotational and dissipative properties of superfluids. They are also important for an understanding of classical fluids and wave excitations of a vortex were therefore already studied by Lord Kelvin. Such excitations have been observed in superfluid helium Hall (1958) and, most recently, in Bose-Einstein condensates by Bretin et al. Bretin et al. (2003). The latter experiment also elucidated the role of the quadrupole modes of a Bose-Einstein condensate in exciting the Kelvin modes of the vortex line. A theoretical understanding of the observed phenomena was provided by Mizushima et al. Mizushima et al. (2003).

The phenomena of superfluidity is closely connected to the existence of quantized vortices and relies on phase coherence across the condensate. Placing the condensate in an optical lattice makes the establishment of phase coherence more difficult and can even induce a transition into a Mott-insulator state. This quantum phase transition was indeed recently observed experimentally Greiner et al. (2002). In a one-dimensional optical lattice the condensate splits into a stack of weakly-coupled disk-shaped condensates, which leads to some intriguing analogs with high-$`T_c`$ superconductors due to their similar layered structure Blatter et al. (1994).

In this Letter we combine the above ideas and study a vortex in a Bose-Einstein condensate in a one-dimensional optical lattice. This problem has not yet been addressed either experimentally or theoretically, but experimental studies are possible with current technical capabilities. Indeed such an experimental setup appears to be a promising way to achieve the quantum Hall regime in a Bose gas, for which we need to have about one vortex per particle Cor . Also outside the quantum Hall regime an optical lattice allows for unique possibilities to study vortices and their quantum fluctuations. In a normal three-dimensional condensate a vortex line is a very complicated dynamical object and it is nearly impossible to tune its properties in a controlled way. Placing the vortex in an optical lattice provides a way, for example, to reduce the stiffness of the the vortex line by increasing the depth of the lattice.

With this goal in mind, we develop a quantum theory for the Kelvin modes in an optical lattice and obtain the Kelvin mode dispersion relation analytically by means of a variational approach. Furthermore, because of their importance for the creation of the Kelvin waves, we study in detail the coupling between the Kelvin modes and quadrupole modes. Our approach leads to a quantum multimode problem with an interaction between Kelvin modes and quadrupole modes that closely resembles squeezing Hamiltonians familiar from quantum optics. We solve this quantum multimode problem and show how the squeezing is reflected in the quantum mechanical uncertainty of the vortex line. In addition, numerical solution of the multimode problem reveals properties of the nonlinear dynamics that cannot be explained with a linear-response theory.

Theory.— The theory we use is based in our previous work van Oosten et al. (2001); Martikainen and Stoof (2003), and we briefly summarize the main points. We consider an elongated Bose-Einstein condensate in a magnetic trap with a radial trapping frequency $`\omega _r`$ that is much larger than the axial trapping frequency $`\omega _z`$. In addition to this cigar-shaped potential the condensate experiences an one-dimensional optical lattice $`V_0\left(𝐫\right)=V_0\mathrm{sin}^2\left(2\pi z/\lambda \right)`$, where $`V_0`$ is the lattice depth and $`\lambda `$ is the wavelength of the laser light. Furthermore, we assume that the lattice is deep enough so that it dominates over the magnetic trapping potential in the $`z`$ direction and that the number of sites are very large. Then we can in first instance ignore the magnetic trap in the $`z`$ direction.

We expand the condensate wave function in terms of wave functions $`\mathrm{\Phi }_n(x,y)\mathrm{\Phi }(z)`$ that are well localized in the sites labeled by $`n`$. We do not yet assume anything about the wave functions $`\mathrm{\Phi }_n(x,y)`$ for the two-dimensional condensates, but for the wave function in the $`z`$ direction, $`\mathrm{\Phi }\left(z\right)`$, we use the ground-state wave function of the harmonic approximation to the lattice potential near the lattice minimum. This harmonic trap has the frequency $`\omega _L=\left(2\pi /\lambda \right)\sqrt{2V_0/m}`$, where $`m`$ is the atomic mass. These considerations lead to an energy functional of the weakly-coupled two-dimensional condensates, which is characterized by the strength of the on-site particle interactions $`4\pi U=4\pi N\left(a/l_r\right)\sqrt{\omega _L/2\pi \omega _r}`$ and by the strength of the nearest-neighbor Josephson coupling $`J=\left(1/8\pi ^2\right)\left(\omega _L/\omega _r\right)^2\left(\lambda /l_r\right)^2\left[\pi ^2/41\right]e^{\left(\lambda /4l_L\right)^2}`$ Martikainen and Stoof (2003). In these formulae $`N`$ is the number of atoms in each site, $`a`$ is the scattering length, and $`l_L=\sqrt{\mathrm{}/m\omega _L}`$. Also, we use trap units so that the unit of energy is $`\mathrm{}\omega _r`$, the unit of time is $`1/\omega _r`$, and the unit of length is $`l_r=\sqrt{\mathrm{}/m\omega _r}`$. Our theory is valid when the lattice depth is much bigger than the chemical potential. Also, the strength of the Josephson coupling should be small, but still large enough to support a superfluid state. The superfluid Mott-insulator transition in one-dimensional optical lattice is expected to occur only for very deep lattices van Oosten et al. (2003), and such a transition has indeed not yet been observed experimentally.

Variational ansatz.— A variational ansatz for the two-dimensional condensates should satisfy two obvious physical criteria. First, the ansatz must predict the vortex precession dynamics correctly. Second, the ansatz must predict the quadrupole mode frequencies with sufficient accuracy. In an earlier publication we used an ansatz that fulfilled the second requirement very well in the full parameter range from the non-interacting limit to the Thomas-Fermi regime Martikainen and Stoof (2003). Unfortunately, this ansatz does not fulfill the first requirement because it does not properly account for the size of the vortex core, which is not important for the collective modes but plays a crucial role in the dynamics of the vortex. Therefore, we use the ansatz

$$\mathrm{\Phi }_n(x,y)\mathrm{exp}\left(\frac{B_0}{2}\left(x^2+y^2\right)\frac{ϵ_n(t)}{2}\left(x^2y^2\right)ϵ_{xy,n}(t)xy+i\mathrm{tan}^1\left(\frac{yy_n(t)}{xx_n(t)}\right)\right),$$
(1)
where $`(x_n,y_n)`$ is the position of the vortex in the $`n`$th site and the variational parameters $`ϵ_n`$ and $`ϵ_{xy,n}`$ are complex. In addition, the square of the wave function is normalized to $`N`$. We have chosen our variational parameters such that only Kelvin modes and quadrupole modes are included. Our theory can be easily adapted to include the monopole mode as well, but we choose to focus on the quadrupole modes due to their experimental importance for the creation of Kelvin waves Bretin et al. (2003). Since our ansatz does not have a core, the average kinetic energy diverges. Therefore, we introduce a small distance cut-off $`\xi `$ as an additional variational parameter. Together with $`B_0`$ this parameter will be determined by minimizing the equilibrium energy functional. This minimization is a simple numerical task that can be done once in the beginning of the calculation. In that manner the condensate with a vortex is slightly bigger than a condensate without a vortex and the vortex core size becomes smaller with increasing atom numbers. Quantitatively, we find to a very good accuracy that $`\xi =\sqrt{B_0}`$.

The ansatz in Eq. (1) fulfills both our requirements with physically realistic parameter values and is still simple enough to be analytically tractable. Our ansatz fails for very weakly-interacting condensates, but provides an accurate description when the interaction strength $`U`$ is larger than about $`10`$. In the other regime our previous ansatz can be used Martikainen and Stoof (2003). In this Letter we do not consider rotation of the trap, but a change to the rotating frame can again be easily included into our approach.

Quantization of the vortex and the quadrupole modes.— If we expand the Lagrangian for our system to second order in the deviations from equilibrium, we are in a position to study the eigenmodes of the system. The eigenmodes we find are the Kelvin modes and two quadrupole modes with quantum numbers $`m=\pm 2`$. To second order there is no coupling between the vortex positions and the quadrupole modes and we can treat the two independently. We transform to momentum space using the convention $`f_k=1/\sqrt{N_s}_ne^{ikz_n}f_n`$, where $`N_s`$ is the number of sites. We then proceed by defining the vortex positions in terms of bosonic creation and annihilation operators as $`\widehat{x}_k=(\widehat{v}_k^{}+\widehat{v}_k)/2\sqrt{NB_0}`$ and $`\widehat{y}_k=i(\widehat{v}_k^{}\widehat{v}_k)/2\sqrt{NB_0}`$ Fetter (1967). With this definition the Lagrangian for the vortex position takes the desired form

$$\widehat{L}_K=\underset{k}{}\left(i\widehat{v}_k^{}\dot{\widehat{v}}_k\omega _K\left(k\right)\widehat{v}_k^{}\widehat{v}_k\right),$$
(2)
where, to leading order in $`B_0`$, the Kelvin mode dispersion is given by $`\omega _K(k)=B_0/2+\mathrm{\Gamma }[0,B_0^2]\left(2J(k)B_0/2\right)`$. In this expression $`J(k)=J\left(1\mathrm{cos}\left(kd\right)\right)`$, $`d=\lambda /2`$, and $`\mathrm{\Gamma }[a,z]`$ is the incomplete gamma function. Our result for the Kelvin mode dispersion relation has some interesting features as seen in Fig. 1 (b). First, the initially negative precession frequency, that indicates that a straight vortex line is unstable without rotation, can change sign if the momentum $`k`$ of the Kelvin mode is large enough. Moreover, the dispersion is quadratic for small momenta. Both features are due to the attractive interactions between neighboring pancake vortices, which is harmonic for small separations and due to the energy cost for having phase differences between sites. The quadratic Kelvin dispersion in a lattice contrasts with the behavior of the vortex line in a three-dimensional bulk superfluid, where the Kelvin mode dispersion relation behaves like $`k^2\mathrm{ln}\left(1/k\xi \right)`$ and has a logarithmic dependence for small $`k`$. Third, the vortex position is “smeared” by quantum fluctuations. In particular we have for the Kelvin mode vacuum state that $`\widehat{x}_n^2=\widehat{y}_n^2=1/\left(4NB_0\right)`$. Therefore, the quantum properties of the vortex become more important in a lattice. This is due to the reduced on-site particle number, as opposed to the total number of particles, and the spreading out of the condensate wave function as the lattice depth is increased.

The quadrupole modes are technically somewhat more complicated, but can be quantized in a similar way as the Kelvin modes. The quadrupole mode frequencies are

$$\omega _{\pm 2}\left(k\right)=\frac{B_0}{A\left(k\right)}+A\left(k\right)\left(\frac{1}{B_0}+\frac{B_0}{4}UB_0\frac{B_0\mathrm{\Gamma }[0,B_0^2]}{2}\right)+2J(k)\left(A\left(k\right)+\frac{1}{A\left(k\right)}\right)\pm B_0,$$
(3)
where

$$A\left(k\right)=\sqrt{\frac{4B_0^2+8B_0J(k)}{4+8B_0J(k)+B_0^2\left(14U2\mathrm{\Gamma }[0,B_0^2]\right)}}.$$
(4)
The Lagrangian for the quadrupole modes is then

$$\widehat{L}_Q=\underset{m=\pm 2}{}\underset{k}{}\left(i\widehat{q}_{m,k}^{}\dot{\widehat{q}}_{m,k}\omega _m\left(k\right)\widehat{q}_{m,k}^{}\widehat{q}_{m,k}\right).$$
(5)
In Fig. 1 (a) we plot the behavior of the mode frequencies at $`k=0`$ as a function of the interaction strength. The quadrupole mode frequencies are in good agreement with the ones calculated using the Bogoliubov-de Gennes equation Martikainen and Stoof (2003) and with realistic interaction strengths the vortex precession frequency is close to the result using a Thomas-Fermi wave-function for the disk-shaped condensate Fetter and Svidzinsky (2001).

Coupled Kelvin modes and quadrupole modes.— We now proceed to apply our previous results to study the coupling between the Kelvin modes and the quadrupole modes. To obtain coupling between the quadrupole modes and the Kelvin modes, we must expand the Lagrangian up to third order. At third order there are contributions from the kinetic energy, the time-derivative term in the Lagrangian, and the Josephson coupling term. Since we assume small $`J`$ we can safely ignore the last one. The remaining terms describe the coupling between the Kelvin modes and both quadrupole modes. The contribution from the time-derivative term in the Lagrangian includes the time-derivative of the Kelvin mode operators. In the rotating-wave approximation we can replace these time-derivatives with the Kelvin mode dispersion. In addition, this eliminates the nonresonant terms from the Hamiltonian, which in particular removes the coupling between the Kelvin modes and the $`m=2`$ quadrupole mode. In this way we obtain the following interaction Hamiltonian for the Kelvin modes and the quadrupole modes

$$\widehat{H}_I=\underset{k,k^{}}{}E_c(k,k^{})\left(\widehat{q}_{2,k^{}}\widehat{v}_{k+k^{}}^{}\widehat{v}_k^{}+\widehat{q}_{2,k^{}}^{}\widehat{v}_{kk^{}}\widehat{v}_k\right),$$
(6)
where $`E_c(k,k^{})=\sqrt{A(k^{})/8N_s}\left(B_0\mathrm{\Gamma }[0,B_0^2]+\omega _K\left(k\right)\right)`$. This expression is accurate as long as the coupling between the Kelvin modes and the quadrupole modes is sharply peaked around some resonant value $`k_0`$. Note that the angular momentum is conserved, since the quadrupole mode has an angular momentum of $`2`$ whereas each Kelvin mode has an angular momentum of $`1`$.

Multimode problem.— Our interaction Hamiltonian shows how the coupling singles out a wavenumber for the Kelvin modes and thus selects the wavelength for the ensuing wiggles in the vortex line Mizushima et al. (2003). We demonstrate this in Fig. 1 (b). The width of the resonance is still unknown, but we can give a simple estimate. Assume that initially we have a straight vortex line and a coherent quadrupole field $`\widehat{q}_{2,k}(0)=q_{2,0}(0)\delta _{k,0}`$. At short times we can assume that the quadrupole field stays constant and then solve the operator equation for the Kelvin modes $`\widehat{v}_{\pm k}`$. It turns out that the solutions are dynamically unstable only when $`|\mathrm{\Delta }\omega \left(k\right)|<|E_c(k,0)q_{2,0}(0)|`$, where $`\mathrm{\Delta }\omega (k)=\omega _2\left(0\right)2\omega _K\left(k\right)`$. Our numerical solutions of the multimode problem confirm this simple estimate with a reasonable accuracy. The number of modes within the interval depends on the number of sites in the lattice and on the quadrupole amplitude, but with realistic parameter values the multimode behavior is more likely than a single-mode behavior.

In linear response we can use Fermi’s golden rule to calculate the transition rate from the initially coherent state $`|\psi _i=|q_{2,0},\left\{v_k\right\}=|q_{2,0},\left\{0\right\}`$. After a short calculation we get the transition rate as

$$\mathrm{\Gamma }_{QK}=\frac{\pi E_c(k_0,0)^2|q_{2,0}\left(0\right)|^2}{2J\mathrm{\Gamma }[0,B_0^2]|\mathrm{sin}\left(k_0d\right)|}.$$
(7)
In Fig. 2 we show an example of the time evolution of the system, when the quadrupole field is treated as a classical field, but Kelvin modes are treated fully quantum mechanically. In the same figure, we also compare the exponential $`q_{2,0}\left(0\right)e^{\mathrm{\Gamma }_{QK}t}`$ with the numerical solution of the nonlinear multimode problem and show the number distribution of Kelvin modes at long times. It should be noted that at these time scales the system appears irreversible due to the large number of modes present. If the number of sites in the lattice is smaller the oscillatory behavior becomes stronger. It should also be realized that, depending on parameters, the result from Fermi’s golden rule can be unreliable since it ignores the finite lifetime of the final states. It is interesting to observe that the number distribution of the Kelvin modes is asymmetric. This is due to the asymmetric nature of the coupling constant $`E_c(k,0)`$ and the detuning $`\mathrm{\Delta }\omega \left(k\right)`$.

The interaction Hamiltonian in Eq. (6) closely resembles the squeezing Hamiltonian encountered for parametric processes in quantum-optics Scully and Zubairy (1997). Therefore, it is not surprising that we can observe a dramatic squeezing under certain conditions. For example, if we only have a single resonant mode then the initial minimum uncertainty state, i.e., the vacuum, will remain a minimum uncertainty state, but can become strongly squeezed. Such a squeezing is reflected in the vortex position distribution. This is shown in Fig. 2 where we also monitor the time-evolution of the deformation $`ϵ(t)=\mathrm{max}\left[\frac{\widehat{x}_n^2(t)\widehat{y}_n^2(t)}{\widehat{x}_n^2(t)+\widehat{y}_n^2(t)}\right]`$ of the vortex position distribution. Because the uncertainty ellipse that characterizes the vortex position distribution precesses in time, the deformation is maximized with respect to the rotation angle of the coordinate system.

Discussion.— We have developed a quantum theory of a vortex in a one-dimensional optical lattice and applied the theory to study the coupled quadrupole-Kelvin mode dynamics. Our solution of the Kelvin mode dispersion relation essentially amounts to a solution of the vortex-line dynamics in an optical lattice. Among other things our study revealed the possibility of squeezed vortex states and asymmetric number distributions for the Kelvin modes. As yet, there have been no experiments on vortices in optical lattices reported. Detailed studies of the Kelvin mode dynamics would require imaging of the three-dimensional vortex line which is difficult, but not impossible Anderson et al. (2001); Bretin et al. (2003). An alternative, but less direct way to probe the system, would be to track the evolution of the quadrupole oscillations. If the resonance is narrow enough our theory predicts revivals of these oscillations. In addition, the multimode nature of the Kelvin mode distribution can cause a decay rate of the quadrupole mode that is quite different from Fermi’s golden rule. Our theory can be applied to a variety of phenomenon. For example, the interplay between the superfluid Mott-insulator transition, the quantum melting of a vortex lattice, and the quantum Hall regime in an optical lattice is a very interesting topic for further research.

This work is supported by the Stichting voor Fundamenteel Onderzoek der Materie (FOM) and by the Nederlandse Organisatie voor Wetenschaplijk Onderzoek (NWO).


text:
# Online Pattern Recognition for the ALICE High Level Trigger

## I Introduction

The ALICE Experiment  at the upcoming Large Hadron Collider at CERN will investigate PbPb collisions at a center of mass energy of about 5.5 TeV per nucleon pair and pp collisions at 14 TeV. Its main tracking detector, the Time Projection Chamber (TPC), is readout by 557568 analog-to-digital channels (ADCs), producing a data size of $``$75 MByte per event for central PbPb collisions and around $``$0.5 MByte for pp collisions at the highest assumed multiplicities .

The event rate is limited by the bandwidth of the permanent storage system. Without any further reduction or compression the ALICE TPC detector can only take central PbPb events up to 20 Hz and min. bias pp events at a few 100 Hz. Significantly higher rates are possible by either selecting interesting (sub)events, or compressing data efficiently by modeling techniques. Both requires pattern recognition to be performed online. In order to process the detector information of 10-25 GByte/sec, a massive parallel computing system is needed, the High Level Trigger (HLT) system.

### I-A Functionality

The HLT system is intended to reduce the data rate produced by the detectors as far as possible to have reasonable taping costs. The key component of the system is the ability to process the raw data performing track pattern recognition in real time. Based on the extracted information, clusters and tracks, data reduction can be done in different ways:

* Trigger: Generation and application of a software trigger capable of selecting interesting events from the input data stream.
* Select: Reduction in the size of the event data by selecting sub-events or region of interest.
* Compression: Reduction in the size of the event data by compression techniques.

As such the HLT system will enable the ALICE TPC detector to run at a rate up to 200 Hz for heavy ion collisions, and up to 1 kHz for pp collisions. In order to increment the statistical significance of rare processes, dedicated triggers can select candidate events or sub-events. By analyzing tracking information from the different detectors and (pre-)triggers online, selective or partial readout of the relevant detectors can be performed thus reducing the event rate. The tasks of such a trigger are selections based upon the online reconstructed track parameters of the particles, e.g. to select events which contain e+e- candidates coming from quarkonium decay or to select events containing high energy jets made out collimated beams of high pt particles . In the case of low multiplicity events such as for pp collisions, the online reconstruction can be used to remove pile-up events from the trigger event.

### I-B Architecture

The HLT system receives data from the front-end electronics. A farm of clustered SMP-nodes ($``$500 to 1000 nodes), based on off-the-shelf PCs and connected with a high bandwidth, low latency network provide the necessary computing power. The hierarchy of the farm has to be adapted to both the parallelism in the data flow and to the complexity of the pattern recognition.

Figure 1 shows a sketch of the architecture of the system. The TPC detector consists of 36 sectors, each sector being divided into 6 sub-sectors. The data from each sub-sector are transferred via an optical fiber from the detector front-end into 216 custom designed *readout receiver cards* (RORCs). Each receiver node is interfaced to a RORC using its internal PCI bus. In addition to the different communication interfaces, the RORCs provide a FPGA co-processor for data intensive tasks of the pattern recognition and enough external memory to store several dozen event fractions. A hierarchical network interconnects all the receiver nodes.

Each sector is processed in parallel, results are then merged in a higher level. The first layer of nodes receive the data from the detector and performs the pre-processing task, i.e. cluster and track seeding on the sub-sector level. The next two levels of nodes exploit the local neighborhood: track segment sending on sector level. Finally all local results are collected from the sectors or from different detectors and combined on a global level: track segment merging and final track fitting.

The farm is designed to be completely fault tolerant avoiding all single points of failure, except for the unique detector links. A generic communication framework has been developed based on the publisher-subscriber principle, which allows to construct any hierarchy of communication processing elements .

## II Online Pattern Recognition

The main task of the HLT system is to reconstruct the complete event information online. Concerning the TPC and the other tracking devices, the particles should ideally follow helical trajectories due to the solenoidal magnetic field of the L3 magnet, in which these detectors are embedded. Thus we model a track by an helix with 5(+1) parameters describing it mathematically. A track is made out of clusters. So the pattern recognition task is extract clusters out of the raw data and to assign them to tracks thereby determining the helix track parameters.

For HLT tracking, we distinguish two different approaches: the *sequential feature extraction* and the *iterative feature extraction*.

The sequential method –corresponding to the conventional way of event reconstruction– first searches the cluster centroids with a *Cluster Finder* and then uses a *Track Follower* on these space points to extract the track parameters. This approach is applicable for lower occupancy like pp and low multiplicity PbPb collisions. However, at larger multiplicities expected for PbPb at LHC, clusters start to overlap and deconvolution becomes necessary in order to achieve the desired tracking efficiencies.

For that reason, the iterative method first determines track candidates using a suitable defined *Track Candidate Finder* and then assigns clusters to tracks using a *Cluster Evaluator* possibly deconvoluting overlapping clusters shared by different tracks. In both cases, a helix fit on the assigned clusters finally determines the track parameters.

In order to reduce data shipping and communicaton overhead within the HLT, as much as possible of the *local* pattern recognition will be done on the RORC. We therefore intend to run the *Cluster Finder* or the *Track Candidate Finder* directly on the FPGA co-processor of the receiver nodes while reading out the data over the fiber. In both cases the results, cluster centroids or track candidate parameters, will be sent from the RORC to the host over the PCI bus.

## III Sequential tracking approach

The classical approach of pattern recognition in the TPC is divided into two sequential steps: Cluster finding and track finding. In the first step the Cluster Finder reconstructs the cluster centroids, which are interpreted as the three dimensional space points produced by the traversing particles. The list of space points is then passed to the Track Follower, which combines the clusters to form track segments. A similar reconstruction chain has successfully been used in the STAR L3 trigger , and thus has been adapted to the ALICE HLT framework.

#### III-1 The Cluster Finder

The input to the cluster finder is a list of above threshold timebin sequences for each pad. The algorithm builds the clusters by matching sequences on neighboring pads. In order to speed up the execution time every calculation is performed on-the-fly; sequence centroid calculation, sequence matching and deconvolution. Hence the loop over sequences is done only once. Only two lists of sequences are stored at every time: The current pad and the previous pad(s). For every new sequence the centroid position in the time direction is calculated by the ADC weighted mean. The mean is then added to a current pad list, and compared to the sequences in the previous. If a match is found, the mean position in both pad and time is calculated and the cluster list is updated. Every time a match is not found, the sequence is regarded as a new cluster.

In the case of overlapping clusters, a crude deconvolution scheme can be performed 111The deconvolution can be switched on/off by a flag of the program. In time direction overlapping sequences are identified by a local minimum in a sequence, and is separated by cutting at the position of the minimum in time direction. The same approach is being used for the pad direction, where the cluster is cut if there is a local minimum of the pad charge values.

The algorithm is inherently local, as each padrow can processed independently. This is one of the main reasons to use a circuit for the parallel computation of the space points on the FPGA of the RORC .

#### III-2 The Track Follower

The tracking algorithm is based on conformal mapping. A space point (x,y) is transformed in the following way:

$$x^{}=\frac{xx_t}{r^2}$$

$$y^{}=\frac{yy_t}{r^2}$$

$$r^2=(xx_t)^2+(yy_t)^2,$$
(1)
where the reference point $`(x_t,y_t)`$ is a point on the trajectory of the track. If the track is assumed to originate from the interaction point, the reference point is replaced by the vertex coordinates. The transformation has the property of transforming the circular trajectories of the tracks into straight lines. Since then fitting straight lines is easier and much faster than fitting circles (if we neglect the changes in the weights of the points induced by conformal mapping), the effect of the transformation is to speed up the track fitting procedure.

The track finding algorithm consists of a follow-your-nose where the tracks are built by including space points close to the fit . The tracks are initiated by building track segments, and the search is starting at the outermost padrows. The track segments are formed by linking space points which are close in space. When a certain number of space points has been linked together, the points are fitted to straight lines in conformal space. The tracks are then extended by searching for clusters which are close to the fit.

#### III-3 Track merging

Tracking can be done either locally on every sub-sector, on the sector level or on the complete TPC. In the first two scenarios, the tracks have to be merged across the detector boundaries. A simple and fast track merging procedure has been implemented for the TPC. The algorithm basically tries to match tracks which cross the detector boundaries and whose difference in the helix parameters are below a certain threshold. After the tracks have been merged, a final track fit is performed in real space.

#### III-4 Tracking performance

The tracking performance has been studied and compared with the offline TPC reconstruction chain. In the evaluation the following quantities has been defined:

* Generated good track – A track which crosses at least 40% of all padrows. In addition, it is required that half of the innermost 10% of the clusters are correctly assigned.
* Found good track – A track for which the number of assigned clusters is at least 40% of the total number of padrows. In addition, the track should not have more than 10% wrongly assigned clusters.
* Found fake track – A track which has sufficient amount of clusters assigned, but more than 10% wrongly assigned clusters.

The tracking efficiency is the ratio of the number of found good tracks to the number of generated good tracks. The identical definitions have been used both for offline and HLT for comparison.

Figure 2 shows the comparison of the integral efficiency of the HLT and offline reconstruction chains for different charged particle multiplicities for a magnetic field of B=$`0.4`$T. We see that up to dN/dy of 2000 the HLT efficiency is $`90`$%, but for higher multiplicities the HLT code becomes too inefficient to be used for physics evaluation. In this regime other approaches have to be applied.

#### III-5 Timing performance

The TPC analysis in HLT is divided into a hierarchy of processing steps from cluster finding, track finding, track merging to track fitting.

Figure 3 shows the foreseen processing hierarchy for the sequential approach. Cluster finding is done in parallel on each Front-End Processor (FEP), whereas track finding and track fitting is done sequentially on the sector level processors. The final TPC tracks are the obtained on the event processors, where the tracks are being merged across the sector boundaries and a final track fit is performed (cmp. to figure 1).

Figure 4 shows the required computing time measured on a standard reference PC 222800 MHz Twin Pentium III, ServerWorks Chipset, 256 kB L3 cache corresponding to the different processing steps for different particle multiplicities. The error bars denote the standard deviation of processing time for the given event ensemble. For particle multiplicity of dN/dy=4000, about 24 seconds are required to process a complete event, or 4800 CPUs are required to date for the TPC alone at an event rate of 200 Hz 333Estimate ignores any communication and synchronization overhead in order to operate HLT. Table I compares the CPU time needed to reconstruct a TPC event (dN/dy=4000) for HLT and offline. For offline, loading the data into memory is also included in the measurement, while the HLT result only included the processing time as memory accesses are done completely transparent by the publisher-subscriber model 444For offline 17 % of the time is needed for data loading..

### III-A Iterative tracking approach

For large particle multiplicities clusters in the TPC start to overlap, and deconvolution becomes necessary in order to achieve the desired tracking efficiencies. The cluster shape is highly dependent on the track parameters, and in particular on the track crossing angles with the padrow and drift time. In order to properly deconvolute the overlapping clusters, knowledge of the track parameters that produced the clusters are necessary. For that purpose the Hough transform is suited, as it can be applied directly on the raw ADC data thus providing an estimate of the track parameters. Once the track parameters are known, the clusters can be fit to the known shape, and the cluster centroid can be correctly reconstructed. The cluster deconvolution is geometrically local, and thus trivially parallel, and can be performed in parallel on the rawdata.

#### III-A1 Hough Transform

The Hough transform is a standard tool in image analysis that allows recognition of global patterns in an image space by recognition of local patterns (ideally a point) in a transformed parameter space. The basic idea is to find curves that can be parametrized in a suitable parameter space. In its original form one determines a curve in parameter space for a signal corresponding to all possible tracks with a given parametric form to which it could possibly belong. All such curves belonging to the different signals are drawn in parameter space. That space is then discretized and entries are stored in a histogram. If the peaks in the histogram exceeds a given threshold, the corresponding parameters are found.

As mentioned above, in ALICE the local track model is a helix. In order to simplify the transformation, the detector is divided into subvolumes in pseudo-rapidity. If one restricts the analysis to tracks originating from the vertex, the circular track in the $`\eta `$-volume is characterized by two parameters: the emission angle with the beam axis, $`\psi `$ and the curvature $`\kappa `$. The transformation is performed from (R,$`\varphi `$)-space to ($`\psi `$,$`\kappa `$)-space using the following equations:

$$R=\sqrt{x^2+y^2}$$

$$\varphi =\mathrm{arctan}(\frac{y}{x})$$

$$\kappa =\frac{2}{R}\mathrm{sin}(\varphi \psi )$$
(2)
Each ADC value above a certain threshold transforms into a sinusoidal line extending over the whole $`\psi `$-range of the parameter space. All the corresponding bins in the histogram are incremented with the corresponding ADC-value. The superposition of these point transformations produces a maximum at the circle parameters of the track. The track recognition is now done by searching for local maxima in the parameter space.

Figure 5 shows the tracking efficiency for the Hough transform applied on a full multiplicity event and a magnetic field of 0.2T. An overall efficiency above 90% was achieved. The tracking efficiency was taken as the number of verified track candidates divided with the number of generated tracks within the TPC acceptance. The list of verified track candidates was obtained by taking the list of found local maxima and laying a road in the rawdata corresponding to the track parameters of the peak. If enough clusters were found along the road, the track candidate was considered a track, if not the track candidate was disregarded.

However, one of the problems encountered with the Hough transform algorithm is the number of fake tracks coming from spurious peaks in the parameter space. Before the tracks are verified by looking into the rawdata, the number of fake tracks is currently above 100%. This problem has to be solved in order for the tracks found by the Hough transform to be used as an efficient input for the cluster fitting and deconvoluting procedure.

#### III-A2 Timing performance

Figure 6 shows a timing measurement of the Hough based algorithm for different particle multiplicities. The Hough transform is done in parallel locally on each receiving node, whereas the histogram adding, maxima finding and merging tracks across $`\eta `$-slices are down sequentially on the sector level. The histograms from the different subsectors are added in order to increase the signal-to-noise ratio of the peaks. For particle multiplicities of dN/dy=8000, the four steps require about 1000 seconds per events corresponding to 200,000 CPUs for 200 Hz event processing rate. It should be noted that the algorithm were already optimized but some additional optimizations are still believed to be possible. However, present studies indicate that one should not expect to gain more than a factor of 2 without using hardware specifics of a given processor architecture.

The advantage of the Hough transform is that it has a very high degree of locality and parallelism, allowing the efficient use of FPGA co-processors. Given the hierarchy of the TPC data analysis, it is obvious that both the Hough transformation and the cluster deconvolution can be performed in the receiver nodes. The Hough transformation is particular I/O-bound as it create large histograms that have to be searched for maxima, which scales poorly with modern processor architectures and is ideally suited for FPGA co-processors. Currently different ways of implementing the above outline Hough transform in hardware are being investigated.

## IV Data modeling and Data compression

One of the mains goals of HLT is to compress data efficiently with a minimal loss of physics information.

In general two modes of data compression can be considered:

* Binary lossless data compression, allowing bit-by-bit reconstruction of the original data set.
* Binary lossy data compression, not allowing bit-by-bit reconstruction of the original data, while remaining however all relevant physical information.

Run-length encoding (RLE), Huffman and LZW are considered lossless compression, while thresholding and hit finding operations are considered lossy techniques that could lead to a loss of small clusters or tail of clusters. It should be noted that data compression techniques in this context should be considered lossless from a physics point of view.

Many of the state of the art compression techniques were studied on simulated TPC data and is presented in detail in . They all result in compression factors of close to 2. However, the most effective data compression can be done by cluster and track modeling, as will be presented in the following.

### IV-A Cluster and track modeling

From a data compression point of the view, the aim of the track finding is not to extract physics information, but to build a data model which will be used to collect clusters and to code cluster information efficiently. Therefore, the pattern recognition algorithms are optimized differently, or even different methods can be used compared to the normal tracking.

The tracking analysis comprises of two main steps: Cluster reconstruction and track finding. Depending on the occupancy, the space points can be determined by a simple cluster finding or require more complex cluster deconvolution functionality in areas of high occupancy (see sec. III and III-A) . In the latter case a minimum track model may be required in order to properly decode the digitized charge clouds into their correct space points.

However, in any case the analysis process is twofold, clustering and tracking. Optionally the first step can be performed online while leaving the tracking to offline, and thereby only recording the space points. Given the high resolution of space points on one hand, and the size of the chamber on the other, would result in rather large encoding sizes for these clusters. However, taking a preliminary zeroth order tracking into account, the space points can be encoded with respect to their distance to such tracklets, leaving only small numbers which can be encoded very efficiently. The quality of the tracklet itself, with the helix parameters that would also be recorded, is only secondary as the tracking is repeated offline with the original cluster positions.

### IV-B Data compression scheme

The input to the compression algorithm is a lists of tracks and their corresponding clusters. For every assigned cluster, the cluster centroid deviation from the track model is calculated in both pad and time direction. This length is quantized with respect the given detector resolution 555The quantization steps has been set to 0.5 mm for pad direction and 0.8 mm for time direction, which is within the range of the intrinsic detector resolution, and represented by a fixed number of bits. In addition the total charge of the cluster is stored. Since the cluster shape itself can be parametrized as a function of track parameters and detector specific parameters, the cluster widths in pad and time is not stored for every cluster. During the decompression step, the cluster centroids are restored, and the cluster shape is calculated based on the track parameters. In tables II and III, the track and cluster parameters are listed together with their respective size being used in the compression.

The compression scheme was applied to a simulated PbPb event with a multiplicity of dN/dy=1000. The input tracks used were tracks reconstructed with the sequential tracking approach (see Fig. 2). The remaining clusters, or the clusters which were not assigned to any tracks during the track finding step, were disregarded and not stored for further analysis 666The remaining clusters mainly originates from very low pt tracks such as $`\delta `$-electrons, which could not be reconstructed by the track finder. Their uncompressed raw data amounts to a relative size of about 20%.. A relative size of 11% for the compressed data with respect to the original set was obtained. In order to evaluate the impact on the physics observables, the data was decompressed and restored cluster processed by the offline reconstruction chain.

In figure 7 the offline tracking efficiency before and after applying the compression is compared. A total loss of $``$1% in efficiency and no significant loss in pt resolution was observed.

However, keeping the potential gain of statistics by the increased event rate written to tape in mind, one has to weigh the tradeoff between the impacts on the physics observables and the costs for the data storage.

For high occupancy events, clusters start to overlap and has to be properly deconvoluted in order to effectively compress the data. In this scenario, the Hough transform or another effective iterative tracking procedure would serve as an input for the cluster fitting/deconvolution algorithm. With a high online tracking performance, track and cluster modeling, together with noise removal, can reduce the data size by a factor of 10.

## V Conclusion

Focusing on the TPC, the sequential approach –consisting of cluster finding followed by track finding– is applicable for pp and low multiplicity PbPb data up to dN/dy of 2000 to 3000 with more than 90% efficiency. The timing results indicate that the desired frequency of 1KHz for pp and 200 Hz for PbPb can be achieved. For higher multiplicities of dN/dy $``$ 4000 the iterative approach using the Circle Hough transform for primary track candidate finding shows promising efficiencies of around 90% but with high computational costs.

By compressing the data using data modeling, results show that one can compress data of up to 10% relative to the original data with a very small impact on the tracking efficiency and the Pt resolution.