Coherent Control of Trapped Bosons

Coherent Control of Trapped Bosons

Analabha Roy and L.E. Reichl
Center for Complex Quantum Systems
Department of Physics
The University of Texas at Austin, Austin, Texas 78712
September 16, 2019

We investigate the quantum behavior of a mesoscopic two-boson system produced by number-squeezing ultracold gases of alkali metal atoms. The quantum Poincare maps of the wavefunctions are affected by chaos in those regions of the phase space where the classical dynamics produces features that are comparable to . We also investigate the possibility for quantum control in the dynamics of excitations in these systems. Controlled excitations are mediated by pulsed signals that cause Stimulated Raman Adiabatic passage (STIRAP) from the ground state to a state of higher energy. The dynamics of this transition is affected by chaos caused by the pulses in certain regions of the phase space. A transition to chaos can thus provide a method of controlling STIRAP.

1 Introduction

There have been significant advancements in techniques for cooling and trapping ultracold atom gases in recent years, facilitating the experimental realization of Bose-Einstein condensation in dilute alkali gases (specifically and ) in 1995 [1] [2] [3] [4]. Since then, numerous studies of these condensates have been accomplished. In addition, experiments have been conducted that have obtained boson systems in number squeezed states from ultracold gases of alkali atoms in optical traps [5]. Thus, it is now possible to create a mesoscopic two-boson system of Sodium or Rubidium. Individual alkali atoms from a BEC reservoir can be subjected to quantum tweezers (Gaussian lasers exploiting the Landau-Zener tunnelling between the reservoir levels and the levels in the laser beam) [6]. Other recent methods include number-squeezing the BEC itself by ”culling” atoms from a trapped condensate down to a sub-poissonian regime, making the number uncertainty small enough to be ignored [5].

A two-boson system can be subjected to a micrometer-scale double well potential using coherent laser beams. The methods that are applicable in the required length scales are overlapping cross-sections of Gaussian lasers [7]. Numerous methods can be used to control the extent of the overlap, such as an Acousto-Optical Modulator (AOM) vibrating with two sound waves [9]. Each sound wave causes an impinging laser beam with a Gaussian cross section to Bragg-diffract at an angle. The separation between the two diffracted beams can be adjusted with relative ease, and consequently so can the length scale of the double well that is generated by focussing the two diffracting beams into parallel beams (the double well lies in the direction lateral to the propagation) [9]. Other applicable techniques for double well generation are small volume optical traps formed by multiple lasers propagating as Gaussian sheets [5], and blue-detuned far-off-resonant laser light to add a potential hill in the middle of a cigar-shaped magnetic trap [10].

A two-boson system, in an optical trap, can be subjected to stimulated Raman scattering. As we shall show, coherent population transfer from the ground state into one of the excited states can be achieved using radiation pulses. Time-modulated (i.e pulsed) electric fields or laser radiation pulses, applied sequentially in a counter intuitive manner, can be used to facilitate this process. If the time scale of the pulse modulation is sufficiently large, the Raman process is adiabatic (called Stimulated Raman Adiabatic Passage or STIRAP) [11]  [12].

The analysis of traditional STIRAP systems involves relatively weak radiation pulses and the rotating wave approximation on the three involved levels [12]. More recently, multilevel transitions have been performed and Floquet theory has been used to tell us how the population evolves in time [13] [14] [15] [16].

We investigate the appearence of avoided crossings contributed by resonating levels other than the ones that the radiation pulses connect. These avoided crossings in the Floquet eigenphases appear due to level repulsion caused by a loss of symmetry/degeneracy (actual crossings) [17], and affect the statistical properties of the spectra, bringing them close to that predicted by random matrix theory. These are connected with the dynamics of the underlying classical system, which undergo a transition from KAM tori to chaos in this region of the paramter stace [17]. Thus, this work also demonstrates the quantum effects of chaos, induced by the radiation, on multilevel transitions in a 2-boson system.

In the following sections, we describe the behavior of two bosons in a one-dimensional double-well potential which we shall model as detailed in Section 2. In Section 3 we will describe the classical dynamics, where the pseudopotential interaction in one-dimension is approximated by a Gaussian potential. In Section 4, we will discuss the quantum eigenstates of this system and compare the quantum phase space of the eigenstates with the classical phase space. In Section 5, we will proceed to describe the dynamics of excitations of this system which will be driven by two sinusoidal electric field pulses applied in sequence. The frequencies of the fields are chosen to connect specific undriven energy levels of the system. Section 5.1 poses the Schrödinger equation for the quantum dynamics for this system. Section 5.2 introduces Floquet theory and the numerical methods implemented to evaluate the Floquet matrix. Section 6 discusses a specific set of parameters for the general STIRAP dynamics discussed in Section 5. In subsections 6.1-6.3, we look at the eigenvalue spectra of the Floquet matrices for two amplitudes and use them to explain the actual dynamics of the system obtained from ab initio numerical solutions of Schrödinger’s equation with the same parameters. We repeat these numerical methods for a different set of parameters in Section 7, where the wells are slightly deeper and an additional resonance exists in the eigenvalue spectrum which affect the coherent excitations. For large pulse amplitudes, the dynamics of the STIRAP process is also noticeably affected by the presence of chaos in the stronger pulse amplitudes in both cases. Concluding remarks are made on Section 8.

2 The Basic Model

Our system consists of two alkali metal bosons confined to a double-well optical potential. The two bosons can be obtained from a cold-atom system confined in a Magneto-Optical Trap (MOT) by number squeezing, followed by laser culling in an optical dipole trap (ODT). If the optical laser is far-detuned from the internal atomic resonances, they can be treated as point particles. The effective interaction between the bosons, in three dimensions, is obtained in the long wavelength approximation to be


where is Planck’s constant, is the s-wave scattering length and is the displacement of the th particle [18] [19]. Therefore, the energy eigenvalues of the two particle system are given by the Schrodinger equation


where is the mass of the particles, is the optical double-well potential and depends only on a single coordinate, and is the potential profile of the MOT. is a symmetrized energy eigenstate of the interacting two-particle system with energy .

The system can be confined in two spacial (radial) directions so that the essential dynamics occurs in the - direction. This can be achieved by using anisotropic magnetic traps with high aspect ratio [20], where the radial trap frequencies are considerably larger than the axial frequency [21]. In that case, the radial part of the wavefunction is is not affected by the interaction and the energy eigenstate can be decomposed into an axial part (along the x-axis) and a radial part such that


where and denotes the noninteracting ground state of the system in the plane. If we assume the trap to be a radial harmonic oscillator of stiffness then, multiplying Eq. (2) by and integrating in all radial coordinates reduces the Hamiltonian to a form which describes motion along the axial direction,




and the energy contains contributions from the radial ground state. Thus, the system is effectively one dimensional.

We will consider the case of two identical bosons confined to a quartic double well potential. The total Hamiltonian for the system is


where is the momentum of the th particle (= 1,2), is the position of the th particle along the x-axis, determines the depth of the double well potential, and and determine its spatial extent.

It is useful to write the Hamiltonian in terms of dimensionless parameters . We introduce a unit of length, , which (as we will see) is appropriate to the systems we consider here. Then let , , , , , , and where , and . If we now drop the primes on the dimensionless parameters, we obtain


The unit of length, , is two orders of magnitude below typical length scales in double wells created with Gaussian lasers [7], as well as Gaussian sheet lasers, or lasers diffracted off of acousto-optical modulators vibrating with two sound waves [9]. Nonetheless, this length scale should be attainable in these setups. For instances, well dimensions in a double well generated using an AOM are proportional to the difference between the acoustic frequencies (typically MHz  [9]), producing a barrier length in the micrometer range. It should be relatively straightforward to reduce that frequency difference by two orders of magnitude and adjusted it to produce a double well of a barrier length of , or an of , as is required here. The lifetime of the magnetic traps in these experiments is approximately s. If we use Rb alkali atoms as our bosons then the value of is about , making typical trap lifetimes (about ) translate to units of . The characteristic energy scale is very small, and corresponds to photons with frequency of about . Therefore, monochromatic coherent radiation at around this frequency is necessary for STIRAP excitations, and can be obtained using noble gas masers [8]. Fig. 1 shows a plot of the quartic double well for well depth .

3 Classical Dynamics of the Undriven System

In order to study the classical dynamics of the system governed by the Hamiltonian in Eqn. (7) the contact potential can be replaced by a Gaussian shaped potential of suitably chosen width such that


We have noticed no discernable difference in the quantum case between using the Gaussian for the interaction and using the delta function provided that the width of the Gaussian is sufficiently small. Too small a width generates unresolvable errors in the numerical solution to the classical dynamics and so an optimum width was chosen at .

In the absence of a time dependent external field, the system is conservative, has two degrees of freedom, and energy is constant. Therefore, the system is confined to a three dimensional surface in a four dimensional phase space. Fig. 2 shows Poincare surfaces of section of the classical phase space of this system for and (attractive interaction) and fixed. The momentum and displacement of particle-1 are plotted each time the trajectory of particle-2 crosses the point with such that is fixed and therefore they show the behavior of momentum and position of particle-1 for that energy. Energy conservation places bounds on the values of the coordinates so the trajectories are confined to a finite region of the phase space.

3.1 Relevant Energies of the Quantum System

Figs. 2.a through 2.c are surfaces of section for three different energies , and , respectively. The energies were chosen to match quantum energy levels that will be connected by STIRAP pulses. The numerical integration was done by the order implicit Runge Kutta (Prince Dormand) method [22] using the appropriate subroutine in the GNU Scientific Library [23]. The integration was done non-adaptively, with a fixed temporal stepsize of . In Fig. 2.a (), we can see several regions of interest. There are three kinds of dynamics at this energy vis-a-vis the relative energies of each particle. In the case that the particles are in separate wells, trajectories exist where they don’t see each other and are therefore the same as that of a single particle in a double well. They are seen as large KAM (Kolmogorov-Arnold-Moser) tori around in Fig. 2.a.1. These trajectories lie between two chaotic regions. Both chaotic regions are produced by the case when the energy of one particle is set to a positive value (thus taking it above the wells), and the energy of the other particle is decreased so that they both add up to (the particles being kept sufficiently far apart at so as to make the interaction negligible at that time). The resultant dynamics is chaotic due to the interaction experienced by the particles when they approach each other during motion. The cases when the particle being strobed has high momentum cause the chaos at the separatrix coupling both wells. The island seen immediately around in Fig. 2.a.1 is the result of a bifurcation that occurs at lower energies, and will be discussed in the next section.

The case when both particles are in the same well are seen as the highly elongated tori around in Fig. 2.a.2 due to the interaction between the particles. In Fig. 2.b.1, we note the disappearance of the bifurcated island immediately around , the result of a bifurcation in reverse (as energy is increased). The left-right asymmetry observed at lower energies is reduced at the energy increases and (apart from the elongated tori), disappears in Figs. 2.c. The region of chaos seen deep inside the potential wells of Fig. 2.a.1 merges with the region of chaos at the separatrix for higher the value of energy in Fig. 2.c.1.

3.2 Lower Energies

Figs. 3.a through 3.f are surfaces of section for six energies, all less than the quantum ground state energy, and shown in decreasing order in energy. Only the half of the phase space is shown. For sufficiently low energies, the two particles either oscillate in two wells independently or together in the same well. The former case is seen in Fig. 3.f where the periodic motion of one particle about is visible at energy . As we increase the total energy, it becomes possible for one particle to break it’s confinement for certain values of it’s initial coordinates and influence the dynamics of the other particle through the interaction. This produces prominent chaotic behavior as seen in the figures. Increasing the total energy further from to (Figs. 3.d to 3.c) produces a bifurcation as the total energy is increased further (see Figs. 3.c to  3.a).

4 Quantum Mechanics of the Interacting System

In this section, we discuss the quantum mechanics of two interacting bosons whose Hamiltonian is given by Eq. (7). We first discuss the basis used to diagonalize the Hamiltonian. We then show configuration space and phase space plots of the key eigenstates of the system.

4.1 Diagonalization of the Hamiltonian

In order to diagonalize the Hamiltonian in Eq. (7), we use the eigenstates of two noninteracting bosons in a hard-wall box as a representation to formulate the matrix elements. We choose a box of width (in dimensionless units) so that an adequate balance is achieved between truncability and accuracy. The Hamiltonian for a single particle in a box is and the energy eigenstates are given by


where .

The 2-particle boson states of the box system are obtained by symmetrizing the 2-particle states to obtain a complete orthonormal basis of symmetrized 2-boson states:


These states are then used to create a Hamiltonian matrix from Eq. (7). The eigenvalues and eigenvectors of the Hamiltonian matrix were determined numerically using the appropriate subroutine for diagonalizing real symmetric matrices in the GNU Scientific Library [23].

4.2 Energy Eigenstates

Fig. 1 shows a plot of the double well system with well depth . It also shows energy levels of the interacting system with the contact interaction chosen to be attractive (so the amplitude is negative). This can be achieved by tuning a homogeneous magnetic field to the Feshbach Resonance of the alkali metal atoms [24] [19].

Figs. 4.a through 4.c are plots of the ground state, , first excited state, , and third excited state of the interacting system, respectively. Each figure shows a plot of the probability density , as well as cross sections of the wavefunctions for . The bosonic character of the states is evident from the fact that they are symmetric under exchange.

We can also compute the phase space distribution of the energy eigenstates and compare this with the classical surfaces of section in Fig. 2 . A phase space distribution of quantum states was first constructed by Wigner [25]. A smoothed version of the Wigner distribution was introduced by Husimi [26] and has proved particularly useful for comparison of classical and quantum phase space distributions. In the representation, the Husimi Function for a symmetrized 2-particle wavefunction is defined as


where and are the centroids of the Gaussian wave packets in the phase space.

In order to calculate the standard deviations and , we follow the same basic prescription as is normally followed for one-dimensional Husimi functions [28], where the Gaussians appearing in the Husimi function are interpreted as harmonic oscillator coherent states. Therefore the standard deviation is the same as that of the harmonic oscillator, with the modification that the ”stiffness” be where is the double well period of motion, making ( is or ). This is generalized to two particles by choosing the single-particle energies that, when added up, come closest to the 2-particle energy. We then proceed to calculate for each particle as a single noninteracting particle with the chosen energy level. A straightforward integration of the classical double well problem shows [27]


where is the complete elliptic integral of the first kind, is given by


and .

Now that the standard deviations can be calculated, we have the prescription for numerically evaluating the full Husimi function. We cannot sketch the full four- dimensional function realistically, of course. However, we can sketch a ”quantum Poincare map” of the Husimi function by strobing a particular value of given with the energy classically conserved at the quantum eigenvalue. Thus, we can plot


where is determined from the condition at the unperturbed Hamiltonian for a particular eigenstate of energy .

Figs. 5.a through 5.c are Husimi plots of the double well system for states , and . They can be compared with the corresponding classical Poincare sections in Figs. 2.a through 2.c. The Husimi plot of is provided on Figs. 5.a.1 and 5a.2. We notice that the highest probabilities are located in the separatrix region, where there is significant chaos in the classical map (Figs.  2.a). However, there is a significant probability for the system to tunnel from the separatrix region to the interior near the well minima. This is where the bifurcated trajectories occur in Figs.  2.a.1 and 2.a.2. All the Husimis are symmetrical, apart from the interaction resonance in each case, under phase space inversion (). The tunnelling probability from the separatrix into the wells is considerably reduced in the Husimi plot of (shown in Figs. 5.c.1 and 5.c.2).

5 Quantum Dynamics of the Driven System

In order to control transitions between energy states of the two boson system, we drive the system with two sequential pulses of maser radiation with carrier frequency, ( ) for the first (second) pulse. These frequencies are determined by the transitions of interest. The masers are projected in the dimension of particle confinement, thus the spacial dependence of the electric field


( and ) in that direction can be treated as linear, given that the wavelength of the radiation pulse(s) are in the microwave range, and the trapping length scales are in the double well. Ignoring purely temporal terms that only contribute an overall phase, the interaction Hamiltonian simplifies to , where is the atomic dipole moment. Thus, the Hamiltonian of the driven system can be written as


where is the Hamiltonian of the non-driven system in Eqn. (7) and the amplitudes () of the driving fields have Gaussian shape.


where () is the maximum amplitude of the first (second) pulse, and the dipole moment and s have been absorbed into the s. The duration of each pulse is controlled by the parameter , where is a measure of the width of each pulse (similar to standard deviation of the Gaussian). The time at which the maximum of the th pulse occurs is .

The Schrödinger equation for the 2-boson system in the presence of the driving field is


Given our numerical expressions for the energy eigenstates of the non-driven system, we can expand in terms of these states so that


The Schrödinger equation can then be written in the form


where is the probability amplitude to find the system in state at time and denotes dipole matrix elements taken with respect to the exact energy eigenstates of the undriven system. Values of the dipole matrix elements are given in Table 1.

1 2 3 4
1 0 -0.108 0 0
2 -0.108 0 -0.053 -0.008
3 0 -0.053 0 0
4 0 -0.008 0 0
5 0 0 0.015 0.002
6 0 0 0 0
7 0.017 0 0 0.015
8 0 0.003 0 0
9 0 0 0 0.001
10 0 0 0 0
Table 1: Dipole Matrix elements for and . The first values are shown here.

5.1 Floquet States

For the case when the amplitude of the Gaussian pulses changes very slowly relative to the period of the carrier frequencies of the pulses, it is possible to use Floquet theory to study the dynamics of the driven system. As was shown in [13][14][16], one can divide the time over which the pulses act into a sequence of time intervals. During each time interval, the amplitude of the pulses is essentially constant while the carrier waves undergo many oscillations. Consider the time window centered at time . The Hamiltonian describing the dynamics during this time can be written


If the two frequencies and are commensurate so that where and are integers, then the Hamiltonian is time-periodic with a period


Because is time periodic, Floquet theory can be used to analyze the dynamics of the system during the time window centered at .

Let us assume that the Schrodinger equation, has a solution of the form


where the state is time-periodic with period and the phase is real. If this is substituted into the Schrodinger equation we obtain the following eigenvalue equation


The state is the th Floquet eigenstate and is the th Floquet eigenphase. The quantity is a Hermitian operator and is called the Floquet Hamiltonian. The Floquet eigenstates form a complete orthonormal basis and the Floquet eigenphases are conserved quantities [17].

The state of the boson system at time can be expanded in a Floquet spectral decomposition as follows


Because the Floquet eigenstates are time-periodic, the state of the system at time is given by


The Floquet evolution operator is therefore given by where


and is a unitary operator. When the operator acts on the state of the boson system it evolves it forward in time by one period of the driving field.

It is possible to compute the matrix elements of the Floquet evolution operator using energy eigenstates of the undriven system as the basis functions. Thus,


The th eigenvalue of this matrix is and the th eigenvector is given by a column matrix with entries . Since we will only have numerical expressions for the eigenvalues , we can only determine the value of the eigenphases modulo .

The numerical computation of the Floquet matrix is achieved as follows. Each column of the matrix can be constructed by solving the Schrödinger equation (with Hamiltonian ) for one period . Each column of the initial state starts with a single entry for and otherwise. The integration is done times with ranging from to . The numerical integrations were performed using the appropriate subroutine for the order Runge-Kutta-Fehlberg method [22] from the GNU Scientific Library [23]. Each different initial condition yields one column of the Floquet matrix at time . Performing these integrations yields an Floquet evolution matrix. Numerically diagonalizing this matrix gives us the Floquet eigenphases and eigenstates. The numerical diagonalization of the non-Hermitian Floquet matrices were performed using the appropriate routine in the IBM™ Engineering and Scientific Subroutines Library (ESSL) [29]. This process can be performed for each value of and the resulting eigenphases and eigenstates plotted as functions of .

In order to determine the appropriate numerical truncation for the evaluation of the Floquet matrices, we used iteratively increasing values of and checked the components of the Floquet states at the value(s) of when the STIRAP amplitudes were the largest ie at or (see Eq. (17)) until the higher components were too small to contribute to the dynamics.We chose as the final truncation.

In the subsequent section we show that coherent transitions between symmetrized two-particle boson states can be achieved for this system. Because of the sparsity on nonzero dipole matrix elements, the simplest transition, induced by the laser pulses, is from the ground state to the fourth level , via the intermediate state . We show the behavior of the system for three different amplitudes of the radiation pulses.

6 Case 1: STIRAP Ladder, First Pulse , Second Pulse

Fig. 1 shows the energy levels of the double well system for well depth and interaction strength . The value of was chosen so the radiation pulses would have carrier wave frequencies and commensurate with each other and so that () would be equal to the energy spacing (), with a high degree of precision. The energy levels shown in the figure are the exact energy eigenvalues of the undriven two-boson symmetrized system.

We plan to use radiation pulses to induce a coherent transition of the two boson system from its ground state to the excited state . At the first pulse connects the levels and with zero detuning. The second pulse connects the levels and . The ratio to eight decimal places.

The dipole moments of these transitions have very different values. The dipole moment that couples the states and is two orders of magnitude smaller than dipole moment that couples the states and (See Table 1). From Eqs. (16) and (17), the amplitude of the first pulse is given by and the amplitude of the second pulse is given by . Because the dipole coupling of the first pulse is so much smaller than that of the first pulse, we will make the electric field amplitude, , of the first pulse considerably larger than that of the second pulse, so that


The amplitudes for the two radiation pulses are plotted in Fig. 6.

The duration of each pulse can be controlled by varying the pulse width parameter . We let denote the total time over which both pulses act on the system. We choose the following values for the pulse parameters


In the sections below, we will study the effect of these radiation pulses on the boson system for two values of . In both of these cases, we set a value for and set a suitable truncation value for the Floquet evolution matrix.

The Floquet eigenphases lie within a fundamental zone (they are determined modulo ) that is taken to be where . They can be plotted as a function of . For closely spaced values of , Floquet eigenstates belonging to different eigenphases, at neighboring values of , will be orthogonal. This can be exploited to tag and follow the evolution of each eigenstate and eigenphase as a function of .

In subsequent sections, we label each Floquet eigenphase based on it’s dominant dependence on the undriven Hamiltonian eigenstates, at . For the three levels, , and , that are connected by the STIRAP pulses, the corresponding Floquet eigenstates have the following structure and labels:

  1. The eigenphase whose corresponding Floquet eigenstate is dominated by the undriven ground state at is labelled as and the Floquet eigenstate as .

  2. The eigenphase whose corresponding Floquet eigenstate is dominated by the undriven state at is labelled as and the Floquet eigenstate as .

  3. The eigenphase whose corresponding Floquet eigenstate is dominated by the undriven state at is labelled as and the Floquet eigenstate as .

  4. The eigenphase whose corresponding Floquet eigenstate is dominated by the undriven state at is labelled as and the Floquet eigenstate as .

The symmetric and antisymmetric Floquet states are induced by the first radiation pulse (which couples the states and ) even though the amplitude of the first radiation pulse may be very small.

The results of the Floquet analysis described above can be compared to the exact dynamics of the system obtained by solving the full Schrödinger equation in Eq. (20) for the exact state of the driven system . In solving for , we will always start at time with the system in the ground state of the undriven Hamiltonian. We can then plot the probability of finding the system in the undriven energy level as a function of time for various values of . We cannot show strobed Husimi plots of the Floquet states as they evolve across , nor can we show classical Poincare maps of the system during those times, since the system has five degrees of freedom during the STIRAP process.

6.1 Pulse Amplitude

Fig. 7.a we plot the Floquet eigenphases as a function of in units of the total pulse time . The eigenphases of interest are the ones involved in the STIRAP process (ie ,,) which lie in the interval in Fig. 7.a. Fig. 7.b shows a magnification of that region. It is clear that the three levels contribute in a manner characteristic of a traditional STIRAP ladder process approximated by a three-level system [12]. A three-level avoided crossing occurs at , and a coherent population transfer takes place from the ground state to the third excited state. This is further confirmed by Figs. 8.a and 8.b. Fig. 8.a shows the evolution of the dependence of on the undriven energy eigenstates as a function of time . A population transfer occurs from the ground state (labelled ””) to the fourth energy level (labelled “”).

Fig. 8.b shows the actual time evolution of the state of the system obtained by solving the Schrödinger equation as a function of for . In Fig. 8.b, we plot the value of as a function of time . The real time evolution is very close to the evolution of the Floquet eigenstate as a function of . This indicates that the evolution is governed by a single Floquet eigenstate and that the process is adiabatic. The small oscillations and deviations of Fig. 8.b from Fig. 8.a can be attributed to nonadiabatic effects [30].

6.2 Pulse Amplitudes

We now set at a higher value of . Fig. 9.a shows the evolution of the Floquet eigenphases as a function of . Fig. 9.b shows a magnification of the region containing eigenphases , , and . We notice the prominence of a new Floquet state (with corresponding eigenphase profile ) which, at , is displaced in value from , and . At , is dominated by . This occurs due to the near-resonance between the transition and the transition (see Fig. 1).

We also note, from Table 1, that the dipole moment () is an order of magnitude higher than the dipole moment of the connected states . Since the resonance is very close to the resonance, the evolution of the corresponding eigenphase is affected and influences the evolution of the eigenphase . However no measureable avoided crossing with occurs. Therefore, and appear to cross (or under go an avoided crossing very closely spaced) and do not contribute anything significant to the dynamics.

There are three avoided crossings that can affect transitions in the system for . First, appears to undergo an avoided crossing with at (see Fig. 10.a.1). Then the same pair of states avoid each other again at (see Fig. 10.b). Finally, the states , , and undergo the standard STIRAP transition at . The dependence of ,, and on the unperturbed energy eigenstates is shown in Fig. 11. The influence of these avoided crossings on these states is clearly seen. These avoided crossings are manifestations of classical chaos in the quantum dynamics as elaborated in the introduction.

We can now use known properties of avoided crossings to analyze this process in more detail. When two Floquet eigenphases and approach and undergo an isolated avoided crossing, the probability that the system switches from one Floquet state to the other can be calculated from the Landau-Zener formula for two-level systems  [31] (note that use of this estimate for multi-level systems assumes that other levels are not significantly involved in the avoided crossing). In our dimensionless units, the Landau-Zener probability is


where is the (minimum) spacing between and at the avoided crossing and is the magnitude of the rate of change (slope) of the Floquet eigenphases in the immediate neighborhood of the avoided crossing. Thus,


where is the slope of the eigenphase curve in the neighborhood of the avoided crossing. If the system switches between the two Floquet states at the avoided crossing, (if ), then the energy eigenstates of the undriven system which contribute the evolution do not change significantly. On the other hand, if , then the system follows a single Floquet state through the avoided crossing, but there can be significant change in the energy eigenstates of the undriven system that contribute to the dynamics.

The value of depends on the duration of the pulses because that determines the slopes of the Floquet eigenphase curves as they enter and leave the avoided crossing. To make this explicit, we can write


where is the time (normalized to ) at which the avoided crossing occurs. The quantity has very weak dependence on . We can now write the Landau-Zener transfer probability in the form


The quantity has weak dependence on . Thus, the transfer probability will be very small if .

We can compute the Landau-Zener probability for the avoided crossings that occur for this case. We use the information from Fig.  10 to calculate . For the avoided crossing between and , shown in Fig. 10.a at , and . Therefore, and we must have to have a small probability that the system will transfer from one Floquet state to the other. For the second avoided crossing between and shown in Fig. 10.b at , and . This means that for an adiabatic passage (no change in Floquet eigenstate) through the avoided crossing.

Fig. 12.a, which has a relatively small value of (), shows no effect of these first two avoided crossings, although it does show the effect of the three-way avoided crossing that occurs about halfway into the total time. For this case, the pulses appear to leave the system in a superposition of Floquet states and . The effect of the avoided crossing at is also absent in Fig. 12 b, where . However, the effect of the avoided crossing at can be seen in the Fig. A complex mixing of states occurs just after this avoided crossing as the large central three-way avoided crossing comes into play, and in the end the system is again left in a superposition of Floquet states and . The avoided crossing at finally starts to manifest itself in Fig. 12 (c), where and is clearly visible in Fig. 12 d, where . Indeed, in Fig. 12.d, the system follows a single Floquet state through the entire process. This is confirmed by comparing the evolution of the Floquet eigenstate in Fig. 11.b to the exact time evolution in Fig. 12.d. They are essentially identical.

7 Case 2: STIRAP Ladder, First Pulse , Second Pulse and nearly tuned to

We now want to show an interesting effect that can occur in a multilevel system. We adjust the shape of the double-well potential so that there is a resonance between the and transitions that is almost exact (to within units of energy). Fig. 13 shows the energy levels of the double well system for wells that are a little deeper than in Case 1. Here, and interaction strength . The energy levels shown in the figure are the exact energy eigenvalues of the undriven two-boson symmetrized system.

The classical dynamics of the system, for these deeper potential wells, is qualitatively the same as in Case 1, depicted in Figs. 2 and Figs. 3. However, deepening the wells lowers the quantum energies. Nonetheless, the classical dynamics that was observed at energies , and in Case 1 is very similar to that seen for the corresponding energies , and for Case 2. For Case 2, that we consider in this section, the energy will play a significant role. In Fig. 14.a.1 () and 14.a.2 () we show surfaces of section for the classical non-driven interacting system for an energy equal to the seventh energy level of Case 2. The chaos is more spread out and the effect of a new bifurcation can be seen.

We use radiation pulses to induce a coherent transition of the two-boson system from its ground state in a manner similar to Case 1. At the first pulse connects the levels and with zero detuning. The second pulse connects the levels and . In this case, the ratio . As in Case 1, there is a significant difference ( orders of magnitude) between dipole moments and (see table 2). Therefore, the peak amplitudes of the first and second pulses are adjusted in accordance with Eq. (29) (see Fig. 6).

1 2 3 4
1 0 -0.103 0 0
2 -0.103 0 0.043 0.003
3 0 0.043 0 0
4 0 0.003 0 0
5 0 0 -0.007 -0.0004
6 0 0 0 0
7 -0.005 0 0 0.002
8 0 0.003 0 0
9 0 0 0 -0.0002
10 0 0 0 0
Table 2: Dipole Matrix elements for and . The first values are shown here.

In Fig. 15.a we plot the Floquet eigenphases as a function of in units of the total pulse time . The fundamental zone has been set to , where , the commensurate frequency, is given by and is calculated from Eq. (22). In this case, the ratio . The eigenphases of interest are labelled in the same manner as in case 1 ie as ,,, and . These eigenphases lie in the interval in Fig. 15.a. Fig. 15.b shows a magnification of that region.

As we can see from Fig. 15.b, all four eigenphases , , and , participate in a complicated set of avoided crossings. Firstly, an avoided crossing between and very near causes them to switch their supports. Thus, is predominantly supported by after this crossing is avoided. The next avoided crossing of importance is the one between and at , which causes the support of to change from to that of viz. . Thus, a complete population transfer from to is possible. The effect of these avoided crossings and the possible behavior of the system, as the radiations pulses pass through the system, can be seen in Fig. 16. Figs. 16.a through 16.d show the time strobed plots of through , respectively, analogous to Fig. 11 of Case 1. It is clear from Fig. 16.a that this unexpected transition from energy level to should be possible to achieve, producing a marked influence of classical chaos in the quantum dynamics as elaborated in the introduction.

In order to obtain a rough estimate of the pulse time needed to achieve true adiabatic behavior, we apply the Landau Zener formula Eqn. (34) in the same manner as previously done in Case 1, even though it was not meant to be applicable when multiple avoided crossings are involved. For the avoided crossing of with at , shown in Fig. 15.b, and . Therefore and we must have to have a small probability that the system transfers from to .

Fig. 17.a through 17.c show the actual time evolution of the system, starting from the ground state at for different values of . For small values of , below the threshold calculated with the Landau Zener formula (Fig. 17.a), we see a partial coherent transfer to , as demonstrated above. When is well above threshold, as in Fig. 17.c, complete population transfer to is achieved and we appear to have reached approximately adiabatic behavior.

8 Conclusions

We have observed some unusual behavior in the classical and quantum dynamics of two bosons in a double well. Chaos in the separatrix region of the classical version of the coupled system corresponds with regions of high probability in the quantum Poincare map. However, a noticeable tunnelling has been observed from the separatrix into the individual wells. We have also demonstrated the feasibility of a controlled excitation of the system into a higher energy state using STIRAP. The STIRAP pulses destroy symmetries and produce chaos that we can detect by observing avoided crossings in the Floquet eigenphase spectrum. The chaos produced by additional resonances produce avoided crossings that can cause coherent population transfer to higher states. Thus, radiation pulses can be used to exert coherent control of the coupled boson system through chaos assisted adiabatic passages, just as has been recorded for systems with lower degrees of freedom.

9 Acknowledgments

The authors wish to thank the Robert A. Welch Foundation (Grant No. F-1051) for support of this work. A.R. thanks Kyunsun Na and Benjamin P. Holder for useful discussions about Floquet theory and the numerical implementation of Floquet analysis. Both authors also thank the Texas Advanced Computing Center (T.A.C.C.) at the University of Texas at Austin for the use of their high-performance distributed computing grid.


  • [1] C. Monroe, W. Swann, H. Robinson and C. Wieman, Phys. Rev. Lett. 65, 1571 (1990)
  • [2] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman and E. A. Cornell, Science 269, 198 (1995)
  • [3] W. Ketterle, K. B. Davis, M. A. Joffe, A. Martin and D. E. Pritchard, Phys. Rev. Lett. 70, 2253 (1993)
  • [4] K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn and W. Ketterle, Phys. Rev. Lett. 75, 3969 (1995)
  • [5] C.-S. Chuu, F. Schreck, T.P. Meyrath, J.L Hansses, G.N. Price, and M.G. Raizen, Phys. Rev. Lett. , 95 260403 (2005).
  • [6] Roberto B. Diener, Biao Wu, Mark Raizen, and Qian Niu, Phys. Rev. Lett., 89 070401 (2002).
  • [7] Artem M. Dudarev, Roberto B. Diener, Biao Wu, Mark G. Raizen, and Qian Niu, Phys. Rev. Lett., 91 010402 (2003).
  • [8] T. E. Chupp, R. J. Hoare, R.L. Walsworth and Bo Wu, Phys. Rev. Lett., 72, 15, 2363 (1994)
  • [9] Y. Shin, M. Saba, T. A. Pasquini, W. Ketterle, D. E. Pritchard, and A. E. Leanhardt, Phys. Rev. Lett., 92 050405 (2004).
  • [10] M. R. Andrews, C. G. Townsend, H.-J. Miesner, D. S. Durfee, D. M. Kurn, and W. Ketterle, Science 275, 637 (1997)
  • [11] Oreg J, F.T. Hioe, and J.H. Eberly, Phys. Rev. A., 29 69 (1984).
  • [12] Nikolay V. Vitanov, Thomas Halfmann, Bruce W. Shore, and Klaas Bergmann, Annu. Rev. Phys. Chem., 52 763 (2001).
  • [13] Kyungsun Na and L.E. Reichl, Phys. Rev. A, 70 063405 (2004).
  • [14] Kyungsun Na and L.E. Reichl, Phys. Rev. A, 72 013402 (2005).
  • [15] Benjamin P. Holder and L.E. Reichl, Phys. Rev. A, 72 043408 (2005).
  • [16] Kyungsun Na, Christof Jung and L.E. Reichl, J. Chem. Phys., 125 034301 (2006).
  • [17] L.E. Reichl, The Transition to Chaos: Conservative Classical Systems and Quantum Manifestations, 2nd Edition, Chapter 7, (Springer-Verlag, Berlin, 2004).
  • [18] H.J. Metcalf and P. van der Straten, Laser Cooling and Trapping, (Springer-Verlag, New York, 1999)
  • [19] C.J. Pethick and H.Smith, Bose-Einstein Condensation in Dilute Gases, (Cambridge University Press, Cambridge, 2002).
  • [20] M. Olshanii, Phys. Rev. Lett., 81 938 (1998).
  • [21] D. S. Petrov, G.V. Shlyapnikov, and J. T. M. Walraven, Phys. Rev. Lett., 85 3745 (2000).
  • [22] J.R. Dormand and P. J. Prince, A family of embedded Runge-Kutta formulae. J. Comp. Appl. Math., 6(1):19, (1980).
  • [23] M.Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, M. Booth, and F. Rossi, GNU Scientific Library Reference Manual, 2nd edition, (Network Theory Ltd.,Bristol BS8 3AL, United Kingdom, 2003).
  • [24] H. Feshbach. Ann. Phys., 19:287, (1962).
  • [25] E. Wigner, Phys. Rev., 40 749 (1932).
  • [26] K. Husimi, Proc. Phys. Math. Soc. Jpn, 22 248, (1940).
  • [27] Reichl, op.cit., p. 471
  • [28] M.  Novaes, Wigner and Husimi functions in the double-well potential., J. Opt. B:Quantum Semiclass, 5:S342, (2003).
  • [29] Document Number GC23-3836-00, Parallel Engineering and Scientific Subroutine Library Guide and Reference,
  • [30] M.V. Berry, Proc. R. Soc London, Ser. A, 429 61 (1990).
  • [31] G. Zener, Proc. R. Soc. London Ser., A 137 696 (1932).
Figure 1: Plot of the double-well potential experienced by each boson in case 1. All units are dimensionless. The energy levels, , and of the interacting two-boson system (interaction strength ) are also sketched, with wavy arrows denoting the levels connected by the STIRAP pulses. Here, .
Figure 2: Classical Poincare maps (of () for and ) for two interacting particles in the double well potential. All units are dimensionless. Here and . The interaction is approximated by an attractive Gaussian potential of width . (a) Energy . (b) Energy . and (c) Energy . A unit area of the phase space equals
Figure 3: Classical Poincare maps for lower energies with the bifurcating resonance magnified. All units are dimensionless.A unit area of the phase space equals .(a) Energy . (b) Energy , (c) Energy , (d) Energy , (e) Energy , and (f) Energy . Note the increased prominence of the smaller resonance as the energy decreases from , as well as the bifurcation in the other resonance as the energy increases from (e).
Figure 4: Plots of energy eigenfunctions for the two interacting bosons in a double well potential. All units are dimensionless. (a1) Contour plot of the probability density . (a2) The cross-section of the wavefunction at . (a3) The cross-section of the wavefunction at . (a4) The cross-section of the wavefunction at . (b1) Contour plot of the probability density . (b2) The cross-section of the wavefunction at . (b3) The cross-section of the wavefunction at . (b4) The cross-section of the wavefunction at . (c1) Contour plot of the probability density . (c2) The cross-section of the wavefunction at . (c3) The cross-section of the wavefunction at . (c4) The cross-section of the wavefunction at