text:
# Kondo Effect of Quantum Dots in the Quantum Hall Regime

## Abstract

We report on the Kondo effect of quantum dots involving the precursor of the Landau level filling factor $`\nu =1`$ state in the quantum Hall regime. We argue that pairs of degenerate single Slater determinant states may give rise to a Kondo effect which can be mapped into an ordinary Kondo effect in a fictitious magnetic field. We report on several properties of this Kondo effect using scaling and numerical renormalization group analysis. We suggest an experiment to investigate this Kondo effect.

The Kondo effect, one of the most extensively studied phenomena in condensed-matter physics Hewson93a , has recently enjoyed a revival in mesoscopic systems. Examples include quantum dots Goldhaber-Gordon98a ; Goldhaber-Gordon98b ; Cronenwett98a ; vanderWiel00a ; JiY00a , quantum point contacts Cronenwett02a , and carbon nanotube coupled to superconductors Buitelaar02a . The main attraction of the Kondo effect in such systems is its tunability, which makes it possible to test various aspects of the Kondo effect that cannot be directly investigated in bulk solids. For example, a scattering phase shift at the Kondo resonance in a quantum dot has been measured using a two-path interferometer JiY00a . In mesoscopic systems the role of magnetic impurities embedded in bulk solids is played by nanoscale “artificial atoms”, like quantum dots or carbon nanotube (as a whole) coupled to electrodes. The Kondo effect arises essentially from the spin-degenerate energy levels associated with a single unpaired electron in the artificial atom, and is accompanied by the Kondo resonance at the Fermi level. In quantum dots a Kondo effect may arise when an external magnetic field induces energy degeneracy Pustilnik00a ; Eto00a ; Izumida01a ; Sasaki00a , while in bulk systems a magnetic field lifts the spin-degeneracy and leads to splitting of the Kondo resonance peak Meir93a ; Wingreen94a .

In this paper we report on the Kondo effect of quantum dots involving the precursor of the Landau level filling factor $`\nu =1`$ state in the quantum Hall regime. An attractive feature of these quantum dots is that single-Slater-determinant states are exact groundstates of the many-body Hamiltonian under certain circumstances YangSRE93a ; MacDonald93a ; YangSRE02a . The Slater-determinant states, denoted by $`|N_{},N_{}`$, have $`N_{}`$ spin-up and $`N_{}`$ spin-down electrons, see Eq. (2). In this work we focus on the parameter region where $`|N,0`$ and $`|N1,1`$ become nearly degenerate ground states YangSRE93a ; MacDonald93a (see Fig. 1). The degeneracy of $`|N,0`$ and $`|N1,1`$ results from the interplay of many-body interaction, parabolic confinement energy, and (real) Zeeman energy. These states are stable ground states in wide regions of the parameter space YangSRE93a ; MacDonald93a  and are easily probed in experiments Oosterkamp99a ; Klein95a . We show that this pair of degenerate many-body ground states give rise to a Kondo effect which can be mapped into an ordinary Kondo problem with a *fictitious* magnetic field. Also our investigation suggests that the effects of the fictitious field may be removed by moving off the phase boundary. Our analysis indicates that a similar effect can arise at each phase boundary between $`|N_{},N_{}`$ and $`|N_{}1,N_{}+1`$ ground states. We report on several properties of this Kondo effect using scaling and numerical renormalization group (NRG) analysis. We also suggest an experiment to investigate this Kondo effect.

We first describe briefly our model for a quantum dot in a strong magnetic field. It consists of two dimensional electrons confined to a finite area that is coupled to two leads. The confining potential is regarded parabolic, $`V(x,y)=\frac{1}{2}m^{}\mathrm{\Omega }^2(x^2+y^2),`$ where $`m^{}`$ is the effective electron mass and $`\mathrm{\Omega }`$ is the frequency scale. A magnetic field $`B`$ is applied perpendicular to the quantum dot, i.e., in the $`z`$-direction. The field is assumed to be so strong that the Landau level spacing $`\mathrm{}\omega _c`$, where $`\omega _c=eB/m^{}c`$, is sufficiently large compared with the confinement energy $`\mathrm{}\mathrm{\Omega }`$ and the electron-electron interaction energy $`e^2/ϵ\mathrm{}`$ ($`\mathrm{}\sqrt{\mathrm{}c/eB}`$ and $`ϵ`$ is the dielectric constant). In this limit, only the lowest Landau level is relevant. In the symmetric gauge, the single-electron orbits of an isolated quantum dot are labeled by the orbital angular momentum $`m`$ ($`m=0,1,2,\mathrm{}`$) and spin index $`\sigma =,`$. An electron in the angular momentum state $`m`$ is located in a ring of the width $`\mathrm{}`$ and a radius $`R_m=\sqrt{2(m+1)}\mathrm{}`$. The single-particle energy of an orbit with $`m`$ and $`\sigma `$ is given by $`ϵ_{m,\sigma }=\gamma (1+m)+\frac{1}{2}g\mu _BB\sigma `$ where $`\gamma \mathrm{}\omega _c(\mathrm{\Omega }/\omega _c)^2`$ and $`g\mu _BB`$ is the Zeeman splitting.

The many-body Hamiltonian of an isolated quantum dot in a strong magnetic field is then given by

$$H_0=\underset{m,\sigma }{}\epsilon _{m,\sigma }d_{m,\sigma }^{}d_{m,\sigma }+\underset{m_1,m_1^{},m_2,m_2^{}}{}\underset{\sigma _1,\sigma _2}{}U_{m_1^{},m_2^{};m_1,m_2}d_{m_2^{},\sigma _2}^{}d_{m_1^{},\sigma _1}^{}d_{m_1,\sigma _1}d_{m_2,\sigma _2},$$
(1)
where $`d_{m,\sigma }^{}`$ ($`d_{m,\sigma }`$) creates (annihilates) an electron in the state $`m,\sigma `$, and $`U_{m_1^{},m_2^{};m_1,m_2}`$ are the matrix elements of the electron-electron Coulomb interaction. As mention above, in quantum dots in a strong magnetic field the single-Slater-determinant state of the form

$$|N_{},N_{}=d_{N_{}1,}^{}\mathrm{}d_0^{}d_{N_{}1,}^{}\mathrm{}d_0^{}|0,$$
(2)
so-called *maximum density droplet states*, can be an exact eigenstate of the many-body Hamiltonian, Eq. (1), in a wide range of parameters $`\stackrel{~}{\gamma }\gamma /(e^2/ϵ\mathrm{})`$ and $`\stackrel{~}{g}=g\mu _BB/(e^2/ϵ\mathrm{})`$ (see Figs. 1 and 2 in Ref. YangSRE93a, ).

In the presence of the coupling to the leads, the many-body states in the quantum dot is hybridized with the conduction bands of the leads. The coupling can be considered within the tunneling model

$$H_T=\underset{k,m,\sigma }{}V\left(d_{m,\sigma }^{}c_{k,\sigma }+c_{k,\sigma }^{}d_{m,\sigma }\right),$$
(3)
where $`c_{k,\sigma }^{}`$ and $`c_{k,\sigma }`$ are operators for conduction electrons, and $`V`$ is the tunneling amplitude (we ignore the $`m`$-dependence of $`V`$ for simplicity, see below). We assume that there are $`N`$ electrons in the quantum dot in equilibrium, ant that the Fermi level $`E_F`$ of the leads lies between the successive electro-chemical potentials $`\mu _{N1}`$ and $`\mu _N`$ of the quantum dot vanHouten92a : $`\mu _{N1}<E_F<\mu _N.`$ Here notice that the electro-chemical potential of the quantum dot includes the contribution from the gate voltage $`V_g`$ applied to the quantum dot, which is given by $`\mu _NE_{N+1}^0E_N^0+eV_g`$ where $`E_N^0`$ is the $`N`$-electron ground-state energy. Below we consider the Kondo effect in the limit $`\mathrm{\Gamma }E_F\mu _{N1}`$ and $`\mathrm{\Gamma }\mu _NE_F`$ ($`\mathrm{\Gamma }2\pi \rho _0|V|^2.`$) We emphasize, however, that the spins involved here is not real spins but pseudo-spins corresponding to the degenerate ground states $`|N,0`$ and $`|N1,1`$.

Given the setup prescribed above, transport of the conduction electrons through the quantum dot via sequential tunneling is highly suppressed. The hybridization of the dot levels are possible only through the virtual tunneling processes. To simplify the discussion, we will further assume that $`\mu _NE_FE_F\mu _{N1}`$. The dominant contributions then come from the virtual processes involving the $`(N1)`$-electron states $`|N1,0`$ and $`|N2,1`$, while the processes involving the virtual states with $`(N+1)`$ electrons in the dot can be ignored. \[In the NRG study to be discussed below, however, we have taken into account all the virtual processes with $`(N1)`$\- and $`(N+1)`$-electron states.\] Effectively, the tunneling Hamiltonian takes the form (see also Fig. 1)

$$\begin{array}{c}H_T=\underset{k}{}V\left(c_k^{}\right|N1,0N,0|\hfill \\ \hfill +c_k^{}|N1,0N1,1|\\ \hfill c_k^{}|N2,1N1,1|)+h.c.\end{array}$$
(4)
Each term in Eq. (4) leads to one of the virtual processes listed in Fig. 2. Within the spirit of the Schrieffer-Wolf transformation taking into account these virtual processes, one can show that for small $`\mathrm{\Gamma }`$ and at low temperatures, the impurity model $`H=H_D+H_c+H_T`$ is equivalent to a Kondo-like model:

$$\begin{array}{c}H=\underset{k,q}{}J_{}\left(S^+c_{q,}c_{k,}^{}+S^{}c_{q,}c_{k,}^{}\right)\hfill \\ \hfill \underset{k,q}{}S^z\left(J_{}c_{q,}c_{k,}^{}J_{}c_{q,}c_{k,}^{}\right),\end{array}$$
(5)
where we have adopted the notation $`S^+=|N,0N1,1|`$, $`S^{}=|N1,1N,0|`$, and $`S^z=(|N,0N,0||N1,1N1,1|)/2`$ to emphasize the roles of the states $`||N,0`$ and $`||N1,1`$ as pseudo-spin components. The coupling constants $`J_{}`$, $`J_{}`$, and $`J_{}`$ in Eq. (5) are given by

$$J_{}=\frac{|V|^2}{E_F+E_{N1,0}E_{N,0}eV_g},$$
(6)
$$\begin{array}{c}J_{}=\frac{|V|^2}{E_F+E_{N1,0}E_{N,0}eV_g}\hfill \\ \hfill \frac{|V|^2}{E_F+E_{N2,1}E_{N1,1}eV_g},\end{array}$$
(7)
and

$$J_{}=\frac{|V|^2}{E_F+E_{N1,0}E_{N1,1}eV_g}.$$
(8)
Before going further, it will be useful to discuss here an important difference between our impurity problem, the model studied in Ref. Pustilnik00a, , and the usual Anderson model. In the Anderson model, the degenerate states are spin states associated with a single *orbital* in the impurity. The impurity spin is isotropic. In particular, the spin components $``$ and $``$ are equal in the contributions to the transport. In our case, the pseudo-spin component $`|`$ ($`|N1,1`$) allows for one additional virtual-process channel compared with $`|`$, namely, the one depicted in Fig. 2 (e), coming from the third term in the parenthesis in Eq. (4). This additional channel is responsible for the difference between the coefficients $`J_{}`$ and $`J_{}`$ in Eq. (5). Physically, this additional channel allows the electrons in the $`|`$ states for more kinetic energy gain than those in the $`|`$ states. The difference $`\mathrm{\Delta }_zJ_{}J_{}`$ leads to an *effective* Zeeman splitting between $`|`$ and $`|`$ states. It is reminiscent of the original Kondo impurity in an external magnetic field Meir93a ; Wingreen94a ; corresp . It should be emphasized, however, that in our case the fictitious field $`\mathrm{\Delta }_z`$ arises *intrinsically*. Moreover, as we see below, the effects of $`\mathrm{\Delta }_z`$ can be removed by detuning the levels of $`|N,0`$ and $`|N1,1`$ (i.e., going off the phase boundary in the parameter space). The two degenerate states involved in the Kondo effect studied in Ref. Pustilnik00a,  become degenerate due to competition between the single-particle orbital-level spacing and the Zeeman energy. In the case of the Kondo effect near the singlet-triplet transition studied in Ref. Eto00a, ; Izumida01a, , the degeneracy purely comes from the many-body exchange interaction, and is four-fold (ignoring very small Zeeman energy). In our case, it is a result of interplay between many-body interaction, parabolic confinement energy, and (real) Zeeman effect. We also remark that the tunneling amplitude $`V`$ in Eq. (3) should, in general, depend on the orbital $`m`$. Its effect introduces another contribution to the fictitious Zeeman splitting, and does not affect the results qualitatively.

To investigate the low-energy physical properties of the Kondo-like model, Eq. (5), we first follow the poor man’s scaling approach Hewson93a  and trace out the conduction electrons in the range $`D\delta D<|\epsilon _k|<D`$. The resultant renormalization group (RG) equations are given by

$$\frac{d}{d\mathrm{}}J_{}=2\rho _0J_zJ_{},\frac{d}{d\mathrm{}}J_z=2\rho _0J_{}^2,\frac{d}{d\mathrm{}}\mathrm{\Delta }_z=0,$$
(9)
where $`J_z(J_{}+J_{})/2`$ and $`\mathrm{}\mathrm{ln}(\rho _0D)`$. The first two RG equations in Eq. (9) are exactly the same as those of the usual Anderson model, which exhibit a Kosterlitz-Thouless-type RG flow diagram Hewson93a . This implies that the effective impurity model in question exhibits a Kondo effect. The additional coupling constant $`\mathrm{\Delta }_z`$ in Eq. (9) plays a role of a fictitious magnetic field on the pseudo-spin. Therefore, for $`\mathrm{\Delta }_z`$ larger than the Kondo temperature $`T_K`$, the Kondo resonance peak at the Fermi level will be splitted into two peaks with separation given by $`\mathrm{\Delta }_z`$. The spectral weights of the splitted peaks will be diminished Meir93a ; Wingreen94a . Further, since $`\mathrm{\Delta }_z`$ arises intrinsically in our model, the effect of this fictitious field may be removed by detuning the levels of $`|N,0`$ and $`|N1,1`$. All these arguments are confirmed by the NRG calculations, see below.

The RG equation, Eqs. (9), and the corresponding arguments above are based on a second-order perturbation in $`J_{}`$, $`J_{}`$, and $`J_{}`$. One may question the validity of the perturbative RG analysis, in particular, the marginal behavior of $`\mathrm{\Delta }_z`$ in Eq. (9). We perform a NRG calculation Bulla01a ; Krishna-murthy80b ; Krishna-murthy80a ; Wilson75a , and justify that the above conclusions are qualitatively correct. Figure 3 shows the local density of states of the quantum dot in a strong magnetic field. For the calculations we have chosen $`\mathrm{\Gamma }=0.025D`$ and $`E_F(E_{N,0}E_{N1,0}+eV_g)=0.085D`$. Near the Fermi level, there are two Kondo peaks separated by $`\mathrm{\Delta }_z0.01D`$. Since this fictitious Zeeman splitting comes from the differences in the kinetic energy gains for $`|`$ and $`|`$ states, it can be removed by detuning the degenerate levels $`E_{N,0}`$ and $`E_{N1,1}`$. In Fig. 4, we have detuned the energy levels of $`E_{N,0}`$ and $`E_{N1,1}`$ by amount $`\mathrm{\Delta }_z`$, i.e., $`E_F(E_{N,0}E_{N1,0}+eV_g)=(0.0850.005)D`$ and $`E_F(E_{N1,1}E_{N1,0}+eV_g)=(0.085+0.005)D`$. One can clearly see that the Kondo resonance peak at the Fermi level has been recovered.

Before concluding we briefly discuss possible experiments. The maximum density droplet states have already been probed in transport measurements through vertical Oosterkamp99a  and lateral Klein95a  quantum dots. The ordinary Kondo effects have also been observed in quantum dots Goldhaber-Gordon98a ; Cronenwett98a ; vanderWiel00a . To observe a Kondo effect involving the maximum density droplet states, vertical quantum dots are particularly convenient. One advantage is that they provide a relatively strong confinement potential. Another important advantage of vertical dots is that the leads can be made of bulk electrodes. In bulk electrodes, the densities of states for both spins near the Fermi level are insensitive to the energy (we assume the Zeeman energy is small compared with the Fermi energy) and the spin-polarization effects can be safely neglected vanderWiel00a . We estimate the typical values of the relevant parameters to be $`N15`$, $`B5\mathrm{T}`$ ($`g\mu _BB0.15\mathrm{meV}`$), $`\mathrm{}\mathrm{\Omega }3\mathrm{meV}`$, $`T_K0.05\text{}0.1\mathrm{K}`$, and $`\mathrm{\Delta }_z0.1\mathrm{meV}`$.

We have reported on the Kondo effect of quantum dots in the quantum Hall regime, involving the precursor of the Landau level filling factor $`\nu =1`$ state. Unlike the ordinary Anderson impurity model, we find that the Kondo resonance splits due to an internally generated fictitious magnetic field. We showed that by detuning the energies of states involved in the Kondo effect, the resonance peak at the Fermi level may be restored. We have argued that these effects may be observed in vertical dots.

We thank Hangmo Yi for very useful discussions. M.-S.C. acknowledges the supports from the Swiss-Korean Outstanding Research Efforts Award (SKORE-A) program and by the KOSEF through the eSSC at POSTECH. S.R.E.Y was supported by the KOSEF Quantum-functional Semiconductor Research Center at Dongguk University and by Grant No. R01-1999-00018 from the interdisciplinary Research program of the KOSEF.


format: text

text:
# Sequential Monte Carlo Samplers

## 1 Introduction

Assume we want to sample from a sequence of probability distributions $`\left\{\pi _n\right\}_{n𝒩}`$ defined on a common measurable space $`E`$ where $`𝒩=\{0,1,\mathrm{},p\}`$ or $`𝒩=`$. As a special case, one can set $`\pi _n=\pi `$ for all $`n𝒩`$. Alternatively the distribution can vary across $`𝒩`$. Similarly to simulated annealing, one could be interested in the sequence of distributions $`\pi _n\left(dx\right)\pi ^{\gamma _n}\left(x\right)dx`$ for an increasing schedule $`\left\{\gamma _n\right\}_{n𝒩}`$ so as to maximize $`\pi `$ or $`\pi _n`$ could be the posterior distribution of a parameter given the data collected till time $`n`$. In this paper, we are interested in sampling this sequence of distributions *sequentially*; that is first sampling from $`\pi _0`$ then $`\pi _1`$ and so on. We will further on refer to $`n`$ as the time index.

The tools favoured by statisticians to achieve this are Markov Chain Monte Carlo (MCMC) methods; see for example (Robert and Casella, 1999). To sample from $`\pi _n`$, MCMC methods consists of building an ergodic Markov kernel $`K_n`$ with invariant distribution $`\pi _n`$ using Metropolis-Hastings (MH) steps, Gibbs steps etc. MCMC have been successfully used in many applications in statistics and physics. When the distribution to sample is multimodal, MCMC samplers can be easily stucked in one mode. A standard approach to improve mixing consists of using interacting parallel MCMC/tempering mechanisms where one runs a MCMC chain on an extended space $`E^N`$ with a specified joint invariant distribution admitting $`\pi _n`$ as a marginal (Geyer and Thompson, 1995). However, MCMC are not well adapted to sequential simulation. At index $`n`$, one needs to wait for the Markov chain with kernel $`K_n`$ to reach its stationary distribution $`\pi _n`$.

We propose here a different approach to sample from $`\left\{\pi _n\right\}_{n𝒩}`$. Our approach is based on Sequential Monte Carlo (SMC) methods (Doucet *et al.*, 2001; Liu, 2001). Henceforth the resulting algorithms will be called SMC samplers. SMC methods have been recently studied and used extensively in the context of sequential Bayesian inference and physics (Doucet *et al.*, 2001; Iba, 2001; Liu, 2001). At a given time $`n`$, the basic idea is to obtain a large collection of $`N`$ ($`N1`$) random samples $`\left\{X_n^{\left(i\right)}\right\}_{i=1,\mathrm{},N}`$ named particles whose marginal distribution is asymptotically ($`N\mathrm{}`$) equal to $`\pi _n`$. These particles are carried forward over time using a combination of Sequential Importance Sampling (SIS) and resampling ideas. This approach is very different from parallel MCMC algorithms where one builds a Markov kernel with a specified joint invariant distribution on $`E^N`$.

Standard SMC algorithms available in the literature do not apply to our problem. Indeed, these algorithms deal with the case where each target distribution of interest $`\pi _n`$ is defined on $`E_n`$ with $`E_{n1}E_n.`$ In (Chopin, 2002), a SMC algorithm is proposed to deal with the case $`E_n=E`$. However, this approach restricts severely the way particles can explore the space. The idea in this paper is different and consists of building an artificial sequence of distributions $`\left\{\stackrel{~}{\pi }_n\right\}_{n𝒩}`$ defined on $`E_n=E^{n+1}`$ with $`\stackrel{~}{\pi }_n`$ admitting a marginal $`\pi _n`$. We are then back to the standard SMC framework. More precisely, $`\stackrel{~}{\pi }_n`$ is defined on $`E^{n+1}`$ by $`\stackrel{~}{\pi }_0\left(dx_0\right)=\pi _0\left(dx_0\right)`$ and

$$\stackrel{~}{\pi }_n\left(d(x_0,\mathrm{},x_n)\right)=\pi _n\left(dx_n\right)\underset{k=1}{\overset{n}{}}L_k(x_k,dx_{k1})$$
(1)
where $`\left\{L_n\right\}_{n𝒩\backslash \left\{0\right\}}`$ is a sequence of auxiliary Markov transition kernels.

Our approach has some connections with Annealed Importance Sampling (AIS) (Neal, 2001) and the algorithm recently proposed in (Cappé *et al.*, 2002) which are detailed in Section 2. However, the framework we present here is more general and allows to derive a whole class of principled integration and genetic-type optimization algorithms based on interacting particle systems. Similarly to MCMC, the efficiency of the algorithms is dependent on the target distributions, the proposal and auxiliary kernels. Nevertheless, generally speaking, one can expect that SMC samplers will outperform MCMC when the distributions to sample are multimodal with well-separated modes. Moreover SMC samplers can be used for sequential Bayesian inference problems with static parameters like those addressed by Chopin (2002).

This paper focuses on the algorithmic aspects of SMC samplers. However, it is worth noting that our algorithms can be interpreted as an interacting particle approximating model of a nonlinear Feynman-Kac flow in distribution space. Under additional assumptions, we provide a nonlinear Markov interpretation of the measure-valued dynamic system associated to the flow $`\left\{\pi _n\right\}_{n𝒩}`$. We show that this interpretation is a natural nonlinear version of the MH algorithm. Many convergence results are available for Feynman-Kac flows and their interacting particle approximations (Del Moral and Miclo, 2000; Del Moral and Miclo, 2001) and, consequently, for SMC samplers. However, the Feynman-Kac flow associated to SMC samplers is such that many known estimates on the asymptotic behaviour of these interacting processes can be greatly improved. Several of these results can be found in (Del Moral and Doucet, 2003).

The rest of the paper is organized as follows. In Section 2, we review a generic SMC algorithm to sample from the sequence of distributions (1). Various settings for this algorithm are presented, some extensions and the connections with previous work are outlined. Section 3 describes the distribution flow associated to our interacting particle approximating model and also presents an original nonlinear Markovian interpretation of this flow. Section 4 applies this class of algorithms to a nonlinear regression problem. Finally, we discuss briefly a few open methodological and theoretical problems in Section 5.

## 2 Sequential Monte Carlo Sampling

### 2.1 A Generic Algorithm

We describe here a generic SMC algorithm to sample from the sequence of distributions $`\left\{\pi _n\right\}_{n𝒩}`$ defined in (1) based on a Sampling Importance Resampling strategy; see (Doucet *et al.*, 2001) for a booklength survey of the SMC literature. Alternative SMC algorithms such as the Auxiliary Particle method of Pitt and Shephard (1999) could also be used.

Further on we will use the notation $`X_{0:k}`$ to denote $`(X_0,\mathrm{},X_k)`$. At time $`n1`$, assume a set of particles $`\left\{X_{0:n1}^{\left(i\right)}\right\}`$ ($`i=1,\mathrm{}N`$) distributed approximately according to $`\stackrel{~}{\pi }_{n1}`$ is available, i.e. the empirical measure

$$\widehat{\stackrel{~}{\pi }}_{n1}\left(dx_{0:n1}\right)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\delta _{X_{0:n1}^{\left(i\right)}}\left(dx_{0:n1}\right)$$

is an approximation of $`\stackrel{~}{\pi }_{n1}.`$ At time $`n`$, we extend the path of each particle according to a Markov111The Markov assumption could be relaxed. kernel $`M_n(x_{n1},dx_n)`$. The resulting path is thus approximately distributed according to $`\stackrel{~}{\pi }_{n1}\left(dx_{0:n1}\right)M_n(x_{n1},dx_n)`$. Importance sampling can then be used to correct for the discrepancy between the sampling distribution and $`\stackrel{~}{\pi }_n\left(dx_{0:n}\right)`$, with the importance weight satisfying

$`G_n(x_{n1},x_n)`$ $`={\displaystyle \frac{\stackrel{~}{\pi }_n\left(dx_{0:n}\right)}{\stackrel{~}{\pi }_{n1}\left(dx_{0:n1}\right)M_n(x_{n1},dx_n)}}`$ (2)
$`={\displaystyle \frac{\pi _n\left(dx_n\right)L_n(x_n,dx_{n1})}{\pi _{n1}\left(dx_{n1}\right)M_n(x_{n1},dx_n)}}`$
and being assumed well-defined. Finally, the particles are resampled according to their importance weights; particles with low weights are discarded whereas particles with high weights are multiplied. The resampled particles are given an equal weight. To sum up, the algorithm proceeds as follows.

Sequential Monte Carlo Sampler

Initialization; $`n=0`$.

Sampling step

* For $`i=1,\mathrm{},N`$, sample $`\stackrel{~}{X}_0^{\left(i\right)}v_0()`$.
* For $`i=1,\mathrm{},N`$, evaluate the normalized weights $`W_0^{\left(i\right)}`$

$$W_0^{\left(i\right)}G_0\left(\stackrel{~}{X}_0^{\left(i\right)}\right)=\frac{\pi _0\left(d\stackrel{~}{X}_0^{\left(i\right)}\right)}{v_0\left(d\stackrel{~}{X}_0^{\left(i\right)}\right)}\text{}\underset{i=1}{\overset{N}{}}W_0^{\left(i\right)}=1.$$

Resampling step

* Multiply/Discard particles $`\left\{\stackrel{~}{X}_0^{\left(i\right)}\right\}`$ with respect to high/low weights $`\left\{W_0^{\left(i\right)}\right\}`$ to obtain $`N`$ particles $`\left\{X_0^{\left(i\right)}\right\}`$.

Iteration $`n;`$ $`n𝒩\backslash \left\{0\right\}`$.

Sampling step

* For $`i=1,\mathrm{},N`$, set $`\stackrel{~}{X}_{0:n1}^{\left(i\right)}=X_{0:n1}^{\left(i\right)}`$ and sample $`\stackrel{~}{X}_n^{\left(i\right)}M_n(\stackrel{~}{X}_{n1}^{\left(i\right)},)`$.
* For $`i=1,\mathrm{},N`$, evaluate the normalized weights $`W_n^{\left(i\right)}`$

$$W_n^{\left(i\right)}G_n(\stackrel{~}{X}_{n1}^{\left(i\right)},\stackrel{~}{X}_n^{\left(i\right)})\text{}\underset{i=1}{\overset{N}{}}W_n^{\left(i\right)}=1.$$

Resampling step

* Multiply/Discard particles $`\left\{\stackrel{~}{X}_{0:n}^{\left(i\right)}\right\}`$ with respect to high/low weights $`\left\{W_n^{\left(i\right)}\right\}`$ to obtain $`N`$ particles $`\left\{X_{0:n}^{\left(i\right)}\right\}`$.

In this algorithm, $`v_0`$ is the initial importance distribution. The resampling step can be done using a standard procedure such as multinomial resampling (Gordon *et al.*, 1993), stratified resampling (Kitagawa, 1996) or minimum entropy resampling (Crisan, 2001). All these resampling schemes are unbiased; that is the number of times $`N_i`$ the particle $`\stackrel{~}{X}_{0:n}^{\left(i\right)}`$ is copied satisfies $`E\left(N_i\right)=NW_n^{\left(i\right)}`$. MCMC steps with invariant distribution $`\stackrel{~}{\pi }_n`$ can also be included after the resampling step (Gilks and Berzuini, 1999).

The complexity of this algorithm is in $`O\left(N\right)`$ and it can be parallelized easily. In practice, the memory requirements are in $`O\left(N\right)`$ too and do not increase over time as one does not need to keep in memory at time $`n`$ the whole paths $`\left\{X_{0:n}^{\left(i\right)}\right\}`$ but only $`\left\{X_n^{\left(i\right)}\right\}`$.

The algorithm can be interpreted as an adaptive importance sampling strategy. Initially, $`v_0`$ is used and the particles with the highest importance weights are multiplied whereas the ones with small weights are discarded. At time $`n`$, new “candidate” particles are sampled according to a proposal distribution kernel $`M_n`$. If $`M_n`$ is a random walk, then the new particles can be interpreted as a local exploration of the distribution. The crucial point is that these candidates are weighted by (2) so as to ensure that after the resampling step their distribution is approximately $`\stackrel{~}{\pi }_n`$. The introduction of the auxiliary kernel $`L_n`$ allows the use of importance sampling without having to compute the marginal distribution $`\pi _{n1}\left(du\right)M_n(u,dx)`$ of the particles $`\left\{\stackrel{~}{X}_n^{\left(i\right)}\right\}`$. Indeed, this marginal importance distribution does not typically admit an analytical expression except when $`M_n`$ is a MCMC kernel of invariant distribution $`\pi _{n1}`$. A similar idea is the basis of the conditional Monte Carlo method described by Hammersley (1956).

### 2.2 Particle Estimates

At time $`n`$, we have the following empirical approximations of $`\pi _n`$ before the resampling step

$$\widehat{\pi }_{n,1}\left(dx\right)=\underset{i=1}{\overset{N}{}}W_n^{\left(i\right)}\delta _{\stackrel{~}{X}_n^{\left(i\right)}}\left(dx\right).$$

and after the resampling step it is equal to

$$\widehat{\pi }_{n,2}\left(dx\right)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\delta _{X_n^{\left(i\right)}}\left(dx\right).$$

For any measure $`\mu `$ and function $`f`$, we will denote $`\mu \left(f\right)=f\left(x\right)\mu \left(dx\right)`$. An estimate of $`\pi _n\left(f\right)`$ is given by

$$f\left(x\right)\widehat{\pi }_{n,1}\left(dx\right)=\underset{i=1}{\overset{N}{}}W_n^{\left(i\right)}f\left(\stackrel{~}{X}_n^{\left(i\right)}\right).$$

or alternatively $`f\left(x\right)\widehat{\pi }_{n,2}\left(dx\right)=\frac{1}{N}_{i=1}^Nf\left(X_n^{\left(i\right)}\right)`$ which has higher variance. If $`\pi _n=\pi `$, then the following estimate can be also used

$$\frac{1}{n}\underset{k=1}{\overset{n}{}}\underset{i=1}{\overset{N}{}}W_k^{\left(i\right)}f\left(\stackrel{~}{X}_k^{\left(i\right)}\right).$$
(3)
Though the particles are statistically dependent, one can show under assumptions given in (Del Moral and Miclo, 2000) that this estimate is consistent as $`N\mathrm{}`$.

The algorithm described above can also be used to compute the ratio of normalizing constants. Indeed, typically the sequence of distributions $`\pi _n\left(dx\right)`$ is only known up to a normalizing constant, i.e. say $`\pi _n\left(dx\right)f_n\left(x\right)dx`$. In this case, the unnormalized importance weights one computes are equal to

$$\stackrel{~}{W}_n^{\left(i\right)}=\frac{f_n\left(\stackrel{~}{X}_n^{\left(i\right)}\right)L_n(\stackrel{~}{X}_n^{\left(i\right)},d\stackrel{~}{X}_{n1}^{\left(i\right)})}{f_{n1}\left(\stackrel{~}{X}_{n1}^{\left(i\right)}\right)M_n(\stackrel{~}{X}_{n1}^{\left(i\right)},d\stackrel{~}{X}_n^{\left(i\right)})}W_n^{\left(i\right)}$$

at time $`n`$. It is possible to obtain an estimate of the ratio of the normalizing constants

$$\frac{Z_n}{Z_{n1}}=\frac{f_n\left(x\right)𝑑x}{f_{n1}\left(x\right)𝑑x}$$

using

$$\widehat{\frac{Z_n}{Z_{n1}}}=\frac{1}{N}\underset{i=1}{\overset{N}{}}\stackrel{~}{W}_n^{\left(i\right)}.$$
(4)
Thus an estimate of $`\mathrm{log}\left(Z_n/Z_0\right)`$ is given by

$$\mathrm{log}\left(\widehat{\frac{Z_n}{Z_0}}\right)=\underset{k=1}{\overset{n}{}}\mathrm{log}\left(\widehat{\frac{Z_k}{Z_{k1}}}\right)=\underset{k=1}{\overset{n}{}}\mathrm{log}\left(\underset{i=1}{\overset{N}{}}\stackrel{~}{W}_k^{\left(i\right)}\right)n\mathrm{log}N.$$

If the resampling scheme used is unbiased, then (4) is also unbiased (Del Moral and Miclo, 2000).

### 2.3 Algorithm Settings

The algorithm presented in the previous subsection is very general. There are many potential choices for $`\{\pi _n,M_n,L_n\}_{n𝒩}`$ leading to various integration and optimization algorithms.

*Homogeneous sequences.* A simple choice consists of setting $`\pi _n=\pi ,`$ $`M_n=M`$ and $`L_n=L`$. In this case, the importance weight (2) is the following generalized MH ratio

$$G(x,x^{})=\frac{\pi \left(dx^{}\right)L(x^{},dx)}{\pi \left(dx\right)M(x,dx^{})};$$
(5)
the standard MH ratio corresponds to $`M=L`$. In this simple case, the particles evolve independently according to a proposal distribution $`M`$, their generalized MH ratio is computed and normalized. The particles are then multiplied or discarded with respect to the value of their normalized MH ratio.

*Sequence of distributions* $`\pi _n`$. It might be of interest to consider non homogeneous sequence of distributions either to move “smoothly” from $`\pi _0=v_0`$ to a target distribution $`\pi `$ through a sequence of intermediate distributions or for the sake of optimization. In the case of integration as suggested by Neal (2001), one can select

$$\pi _n\left(dx\right)\pi ^{\gamma _n}\left(x\right)\pi _0^{1\gamma _n}\left(x\right)dx$$

with $`𝒩=\{0,\mathrm{},p\}`$,$`\gamma _0=0`$ and $`\gamma _p=1`$. For the case of optimization, one can select

$$\pi _n\left(dx\right)\pi ^{\gamma _n}\left(x\right)dx$$

where $`𝒩=`$, $`\left\{\gamma _n\right\}_{n0}`$ is an increasing sequence such that $`\gamma _n\mathrm{}`$. In this case, the resulting algorithm is a genetic algorithm where the sampling step is the “mutation” step and the resampling step is the selection step (Goldberg, 1989). However, there is a significant difference with standard genetic algorithms as we know the asymptotic ($`N\mathrm{}`$) distribution of the particles. This makes the analysis of the resulting algorithm easier than in cases where this distribution is unknown such as in (Del Moral and Miclo, 2003). Convergence properties of the algorithm are currently under study.

Finally, another application of this algorithm consists of estimating the sequence of posterior distributions $`\pi _n\left(dx\right)=\pi _n\left(dx|y_{1,\mathrm{}},y_n\right)`$ where $`y_n`$ is an observation available at time $`n`$. As briefly discussed in the introduction, SMC algorithms have been recently proposed in this framework by Chopin (2002) but the SIS framework used is somehow restricted: it only allows $`M_n`$ to be a MCMC kernel of invariant distribution $`\pi _{n1}`$.

*Sequence of proposal kernels* $`M_n`$ *and* *auxiliary kernels* $`L_n`$. Any couple of kernels can be used as long as the ratio (5) is well defined. However, one can only expect good properties of the algorithm if this ratio admits a reasonable variance and also if $`L_n`$ is mixing. Indeed, loosely speaking, the faster $`L_n`$ mixes, the faster the SMC algorithm forgets Monte Carlo errors (Del Moral and Doucet, 2003).

In SMC algorithms (Doucet *et al.*, 2001), it is known that the importance sampling distribution minimizing the conditional variance of the weights at time $`n`$, i.e. $`\left\{X_{0:n1}^{\left(i\right)}\right\}`$ fixed, is given by

$$M_n(x,dx^{})=\frac{\pi _n\left(dx^{}\right)L_n(x^{},dx)}{\pi _n\left(du\right)L_n(u,dx)}.$$
(6)
In this case, the importance weight $`G_n(x,x^{})`$ is given by

$$G_n\left(x\right)=\frac{\pi _n\left(du\right)L_n(u,dx)}{\pi _{n1}\left(dx\right)}$$
(7)
and is independent of $`x^{}`$. This allows the resampling step to be performed before the sampling step.

In standard applications of SMC algorithms, the kernel $`L_n`$ is usually given by the problem at hand whereas in our setup this kernel is arbitrary and can be optimized for a given proposal distribution $`M_n`$. One can alternatively select the kernel $`L_n`$ so as to be able to compute (7); e.g. a MCMC kernel of invariant distribution $`\pi _n`$, and then sample the particles according to (6).

For a fixed $`M_n`$, an alternative natural choice222thanks to C. Andrieu consists of choosing

$$L_n(x,dx^{})=\frac{\pi _{n1}\left(dx^{}\right)M_n(x^{},dx)}{\pi _{n1}\left(du\right)M_n(u,dx)}.$$
(8)
In this case, the associated importance weight $`G_n(x,x^{})`$ is given by

$$G_n\left(x^{}\right)=\frac{\pi _n\left(dx^{}\right)}{\pi _{n1}\left(du\right)M_n(u,dx^{})}.$$
(9)
If $`M_n`$ is a MCMC kernel of invariant distribution $`\pi _{n1}`$, then the weight (9) can be computed easily. If not, numerical integration using the current set of particles can be used to approximate it but the resulting algorithms would be of complexity $`O\left(N^2\right)`$.

### 2.4 Connections to previous work and Extensions

*Connections to previous work.* AIS is a method proposed recently by Neal (2001). Reversing the time index in (Neal, 2001) to be consistent with our notation, AIS corresponds to the case where one considers a finite sequence of distributions, $`M_n`$ is a MCMC kernel of invariant distribution $`\pi _n`$ and

$$L_n(x,dx^{})=M_n(x^{},dx)\frac{\pi _n\left(dx^{}\right)}{\pi _n\left(dx\right)}.$$
(10)
For a given $`M_n`$, one can check that this choice of $`L_n`$ ensures that (6) is satisfied. In this case, one obtains by combining (5) and (10)

$$G_n(x,x^{})=G_n\left(x\right)=\frac{\pi _n\left(dx\right)}{\pi _{n1}\left(dx\right)}.$$

The resampling step is not included in the AIS algorithm. In our framework, we point out that this is a crucial step to include to make the method efficient as established theoretically in (Del Moral and Miclo, 2000) and practically in our simulations. Otherwise the method is just a special instance of SIS and collapses if $`n`$ is too large. In (Godsill and Clapp, 2001), the authors used the AIS algorithm in combination with resampling in the context of optimal filtering.

A more recent work (Cappé *et al.*, 2002) contemporary of (Del Moral and Doucet, 2003) and developed independently is another special case of our framework. In (Cappé *et al.*, 2002), the authors consider the homogeneous case. Their algorithm corresponds to the case where $`M`$ is an MCMC kernel of invariant distribution $`\pi `$ (namely a Gibbs sampler) and $`L(x,dx^{})=\pi \left(dx^{}\right)`$, it follows that

$$G(x,x^{})=\frac{\pi \left(dx^{}\right)}{M(x,dx^{})}.$$

This particular case has limited applications as $`G(x,x^{})`$ would not be defined in most applications; e.g. $`\pi \left(dx^{}\right)=\pi \left(x^{}\right)dx^{}`$ and $`M`$ is an MH kernel.

*Extensions.* The algorithm described in this section must be interpreted as the basic element of more complex algorithms. It is what the MH algorithm is to MCMC. For complex MCMC problems, one typically uses a combination of MH steps where the $`n_x`$ components of $`x`$ say $`(x_1,\mathrm{},x_{n_x})`$ are updated by subblocks (Robert and Casella, 1999). Similarly, to sample from high dimensional distributions, a practical SMC sampler can update the components of $`x`$ via subblocks. There are also numerous potential extensions:

* It is straightforward to develop a version of the algorithm so as to sample distributions defined on an union of subspaces of different dimensions. However, contrary to reversible jump MCMC algorithms (Green, 1995), no reversibility condition is needed.
* As suggested in (Crisan and Doucet, 2000), one can use a proposal kernel whose parameters are a function of the whole set of current particles. This allows the algorithm to automatically scale the proposal distribution based on the previous importance weights.
* In the general case, the sequence of probability distributions $`\left\{\pi _n\right\}_{n𝒩}`$ of interest is such that $`\pi _n`$ is defined on $`E_n`$ and not on $`E`$. We can generalize the algorithm described in this section to this case. We introduce an auxiliary kernel $`L_n`$ from $`E_n`$ to $`E_{n1}`$ and a proposal kernel $`M_n`$ from $`E_{n1}`$ to $`E_n`$. At time $`n1`$, $`N`$ particles $`\left\{X_{n1}^{\left(i\right)}\right\}`$ approximately distributed according to $`\pi _{n1}`$ are available. At time $`N`$ new particles $`\left\{X_n^{\left(i\right)}\right\}`$ are sampled according to $`X_n^{\left(i\right)}M_n(X_{n1}^{\left(i\right)},)`$ and the following importance weights are computed

$$W_n^{\left(i\right)}\frac{\pi _n\left(dX_n^{\left(i\right)}\right)L_n(X_n^{\left(i\right)},dX_{n1}^{\left(i\right)})}{\pi _{n1}\left(dX_{n1}^{\left(i\right)}\right)M_n(X_{n1}^{\left(i\right)},dX_n^{\left(i\right)})}.$$

Then the particles are resampled.

## 3 Feynman-Kac Representation and Particle Interpretations

In this Section, we show that the algorithm presented in Section 2 corresponds to an interacting particle approximation model of a nonlinear Feynman-Kac flow in distribution space. We provide an alternative nonlinear Markovian representation of this flow and its interacting particle approximation. Here the Feynman-Kac flow corresponds to the special case where the so-called potential function is given by the generalized Metropolis ratio (5). The abstract description and the analysis of general Feynman-Kac flows and their particle approximations have been investigated in several recent research articles. Many asymptotic ($`n\mathrm{}`$ and/or $`N\mathrm{}`$) results are available in this field including empirical process convergence, central limit theorems, large deviation principles as well as increasing propagation of chaos estimates and uniform convergence estimates with respect to the time parameter; all of which can be used for SMC samplers. The interested reader is referred to the survey article (Del Moral and Miclo, 2000) and the more recent studies (Del Moral and Miclo, 2001; Del Moral and Miclo, 2003). As mentioned in the introduction, the particular choice of the potential function (5) simplifies the analysis and many known estimates on the asymptotic behaviour of these interacting processes can be greatly improved. Several of these results can be found in (Del Moral and Doucet, 2003).

### 3.1 Feynman-Kac Representation

Define the following distributions on $`E^2=E\times E`$

$$\left(\pi _n\times L_n\right)\left(d(x,x^{})\right)=\pi _n\left(dx^{}\right)L_n(x^{},dx).$$

Using (1) and (2), it is clear that the sequence of distributions $`\left\{\pi _n\times L_n\right\}_{n𝒩}`$ (with the convention $`\pi _0\times L_0=\pi _0\times \pi _0`$) admits the following so-called Feynman-Kac representation

$$\left(\pi _n\times L_n\right)\left(f\right)=\lambda _n\left(f\right)/\lambda _n\left(1\right),$$

with

$$\lambda _n\left(f\right)=𝔼_{v_0,\left\{M_k\right\}}\left(f(X_{n1},X_n)G_0\left(X_0\right)\underset{k=1}{\overset{n}{}}G_k(X_{k1},X_k)\right),$$

where $`𝔼_{v_0,\left\{M_k\right\}}`$ denotes the expectation with respect to

$$v_0\left(dx_0\right)\underset{k=1}{\overset{n}{}}M_k(x_{k1},dx_k).$$

This representation is at the core of the results given in (Del Moral and Miclo, 2000). We give now two “operator-like” interpretations of the sequence $`\left\{\left(\pi _n\times L_n\right)\right\}_{n𝒩}`$. For a measure $`\mu `$ and a Markov kernel $`K`$, we use the standard notation

$$\mu K\left(A\right)=_A\mu \left(dz\right)K(z,dz^{}).$$

Let $`𝒫\left(E^2\right)`$ be the set of probability measures on $`E^2`$. The mapping $`\mathrm{\Psi }_n:𝒫\left(E^2\right)𝒫\left(E^2\right)`$ is defined as

$$\pi _n\times L_n=\mathrm{\Psi }_n\left(\left(\pi _{n1}\times L_{n1}\right)\stackrel{~}{M}_n\right)$$
(11)
where

$$\mathrm{\Psi }_n\left(\mu \right)\left(d(u,v)\right)=\frac{\mu \left(d(u,v)\right)G_n(u,v)}{\mu \left(G_n\right)}$$
(12)
and $`\stackrel{~}{M}_n`$ is a Markov kernel on $`E^2`$ defined as

$$\stackrel{~}{M}_n((u,v),d(u^{},v^{}))=\delta _v\left(du^{}\right)M_n(u^{},dv^{}).$$

Assuming that $`G_n`$ can be upper bounded over $`E^2`$, one can easily check that an alternative representation is given by

$$\left(\pi _n\times L_n\right)=\left(\pi _{n1}\times L_{n1}\right)\stackrel{~}{M}_nS_{{}_{n}{}^{},(\pi _{n1}\times L_{n1})\stackrel{~}{M}_n}$$
(13)
where

$$S_{n,\mu }((u,v),d(u^{},v^{}))=ϵG_n(u,v)\delta _{(u,v)}\left(d(u^{},v^{})\right)+\left(1ϵG(u,v)\right)\mathrm{\Psi }_n\left(\mu \right)\left(d(u^{},v^{})\right),$$
(14)
$`ϵ`$ being chosen such that $`ϵG_n(u,v)1`$ over $`E^2`$.

The kernel $`\stackrel{~}{M}_nS_{{}_{n}{}^{},(\pi _{n1}\times L_{n1})\stackrel{~}{M}_n}`$ is a so-called nonlinear Markov kernel; i.e. the transition kernel is dependent not only on the current state but also on its distribution. A generic nonlinear Markov chain $`\left\{Z_n\right\}_{n0}`$ satisfies

$$Z_nK_{n,\text{Law}\left(Z_{n1}\right)}(Z_{n1},).$$

It is typically impossible to simulate a realization from such a Markov chain as the distribution of the state is not available. However, a particle approximation of it can be used. Consider a Markov chain $`\left\{Z_n\right\}_{n0}`$ taking values in $`E^2`$ with transition kernel $`\stackrel{~}{M}_nS_{{}_{n}{}^{},(\pi _{n1}\times L_{n1})\stackrel{~}{M}_n}`$. This kernel can be interpreted as follows. Given $`Z_{n1}=(U_{n1},V_{n1})\left(\pi _{n1}\times L_{n1}\right)`$, one first sample a candidate $`Z_n^{}=(U_n^{},V_n^{})=(V_{n1},V_n^{})`$ where $`V_n^{}M_n(U_n^{},)`$. With probability $`ϵG_n(U_n^{},V_n^{})`$, one sets $`Z_n=(U_n^{},V_n^{})`$, otherwise $`Z_n\mathrm{\Psi }_n\left(\left(\pi _{n1}\times L_{n1}\right)\stackrel{~}{M}_n\right)`$. By construction, one has $`Z_n\left(\pi _n\times L_n\right)`$. This algorithm can be interpreted as a nonlinear version of the MH algorithm. The main difference being that, when a candidate is rejected, the chain does not stay where it is a new state is proposed according to $`\mathrm{\Psi }_n\left(\left(\pi _{n1}\times L_{n1}\right)\stackrel{~}{M}_n\right)`$.

### 3.2 Particle Interpretations

The first particle interpretation of the flow follows (11)-(12). It corresponds to the standard algorithm which has been described in Section 2. The second alternative algorithm corresponds to a particle interpretation of the flow corresponding to (13). It proceeds as follows.

Iteration $`n;`$ $`n𝒩\backslash \left\{0\right\}`$.

Sampling step

* For $`i=1,\mathrm{},N`$, set $`X_{0:n1}^{\left(i\right)}=\stackrel{~}{X}_{0:n1}^{\left(i\right)}`$ and sample $`\stackrel{~}{X}_n^{\left(i\right)}M_n(\stackrel{~}{X}_{n1}^{\left(i\right)},)`$.
* For $`i=1,\mathrm{},N`$, evaluate the normalized weights $`W_n^{\left(i\right)}`$

$$W_n^{\left(i\right)}G_n(\stackrel{~}{X}_{n1}^{\left(i\right)},\stackrel{~}{X}_n^{\left(i\right)})\text{}\underset{i=1}{\overset{N}{}}W_n^{\left(i\right)}=1.$$
(15)

Resampling step

* $`J=\mathrm{}.`$
* For $`i=1,\mathrm{},N`$, with probability $`ϵG_n(\stackrel{~}{X}_{n1}^{\left(i\right)},\stackrel{~}{X}_n^{\left(i\right)})`$, set $`X_{0:n}^{\left(i\right)}=\stackrel{~}{X}_{0:n}^{\left(i\right)}`$ otherwise set $`J=J\left\{i\right\}.`$
* Multiply/Discard particles $`\left\{\stackrel{~}{X}_{0:n}^{\left(i\right)}\right\}`$ with respect to high/low weights $`\left\{W_n^{\left(i\right)}\right\}`$ to obtain $`\left\{X_{0:n}^{\left(i\right)}\right\}_{iJ}`$.

## 4 Simulation Results

### 4.1 Model

We consider the following harmonic regression model (Andrieu and Doucet, 1999)

$$Y=D\left(\omega \right)\beta +n,$$

where $`Y=(y_0,\mathrm{},y_{m1})^𝚃`$, $`\beta =(\beta _1,\mathrm{},\beta _{2k})^𝚃,`$ $`n=(n_0,\mathrm{},n_{m1})^𝚃`$, $`\omega =(\omega _1,\mathrm{},\omega _k)^𝚃(0,\pi )^k`$ and $`D\left(\omega \right)`$ is a $`m\times 2k`$ matrix where for $`i=0,\mathrm{},m1,j=1,\mathrm{},k`$ .

$$\left[D\left(\omega \right)\right]_{i+1,2j1}=\mathrm{cos}\left(\omega _ji\right),\text{ }\left[D\left(\omega \right)\right]_{i+1,2j}=\mathrm{sin}\left(\omega _ji\right).$$

We assume that $`n|\sigma ^2𝒩(0,\sigma ^2I_p)`$ and we use the following prior $`p(\sigma ^2,\beta ,\omega )=p\left(\omega \right)p\left(\beta |\sigma ^2\right)p\left(\sigma ^2\right)`$ with

$$\sigma ^2𝒢(\frac{\upsilon _0}{2},\frac{\gamma _0}{2}),\text{ }\beta |\sigma ^2𝒩(0,\sigma ^2𝚺_0),$$

where $`𝚺_0^1=\delta ^2D^𝚃\left(\omega \right)D\left(\omega \right)`$ ($`\delta ^2=25`$), $`\upsilon _0=\gamma _0=1`$; $`p\left(\omega \right)`$ is uniform on $`\mathrm{\Omega }=\left\{\omega (0,\pi )^k;0<\omega _1<\mathrm{}<\omega _k<\pi \right\}`$. The posterior density satisfies on $`\mathrm{\Omega }`$

$$p\left(\omega |Y\right)\left(\gamma _0+Y^𝚃PY\right)^{\frac{p+\upsilon _0}{2}}$$

with

| $`M^1=\left(1+\delta ^2\right)D^𝚃\left(\omega \right)D\left(\omega \right),`$ |
| --- |
| $`m=MD^𝚃\left(\omega \right)Y,`$ |
| $`P=I_pD\left(\omega \right)MD^𝚃\left(\omega \right).`$ |

We simulate a realization of $`m=100`$ observations with $`k=6`$, $`\sigma ^2=5`$,

$`\omega `$ $`=(0.08,0.13,0.21,0.29,0.35,0.42)^𝚃,`$
$`\beta `$ $`=(1.24,0.00,1.23,0.43,0.67,1.00,1.11,0.39,1.31,0.16,1.28,0.13)^𝚃.`$
The posterior density is multimodal with well-separated modes.

### 4.2 Algorithms

To sample from $`\pi \left(\omega \right)=p\left(\omega |Y\right),`$ we use an homogeneous SMC sampler with $`N=1000`$ particles where the $`k`$ components are updated one-at-a-time using a simple Gaussian random walk proposal $`M`$ of standard deviation $`\sigma _{RW}`$. We select $`L`$ to be equal to $`M`$ and use the stratified resampling procedure. We compare our algorithm with a MCMC algorithm. The MCMC algorithm updates the component one-at-a-time using a MH step with the proposal kernel $`M`$. In both case, the initial distribution is the uniform distribution on $`\mathrm{\Omega }.`$

We consider the case where $`\sigma _{RW}=0.1`$. Obviously one could come up with a better proposal kernel. We want to emphasize here that the SMC approach is more robust to a poor scaling of the proposal. A similar remark was made in (Cappé *et al.*, 2002). In Figure 1, we present the marginal posterior distributions of $`\omega _1`$ and $`\omega _2`$ obtained using the SMC sampler with $`100`$ iterations. We then run $`12000`$ iterations of the MCMC algorithm so as the computational complexity to be roughly the same for the two algorithms. The MCMC algorithm is more sensitive to the initialization. On $`50`$ realizations of the SMC and the MCMC algorithm, the SMC always explores the main mode whereas the MCMC algorithm converges towards it only $`36`$ times.

We also use an inhomogeneous version of the SMC sampler so as to optimize $`p\left(\omega |Y\right)`$. In this case the target density at time $`n`$ is $`\pi _n\left(\omega \right)p^{\gamma _n}\left(\omega |Y\right)`$ with $`\gamma _n=n`$ and we use $`50`$ iterations. We compare this algorithm to a simulated annealing version of the MH algorithm with $`60000`$ iterations with $`\gamma _n=n/1200`$. In Table 1, we display the mean and standard deviations of the log-posterior density of the posterior mode estimate; the posterior mode estimate being chosen as the sample generated during the simulation maximizing the posterior density. Contrary to the simulated annealing algorithm, the SMC algorithm converges consistently towards the same mode.

| Algorithm | SMC | MCMC |
| --- | --- | --- |
| Mean of the log-posterior values | -326.12 | -328.87 |
| Standard deviation of the | 0.12 | 1.48 |
| log-posterior values |  |  |

Table 1: Performance of SMC and MCMC algorithm obtained over 50 simulations

## 5 Discussion

In this article, we have presented a class of methods to sample from distributions known up to a normalizing constant. These methods are based on SMC algorithms. This framework is very general and flexible. Several points not discussed here are detailed in (Del Moral and Doucet, 2003).

* In the homogeneous case, assume that we do not initialize the algorithm in the stationary regime, i.e. we do not correct for the discrepancy between $`v_0`$ and $`\pi `$. This has to be paralleled with MCMC algorithms which are not initialized in the stationary regime. Under regularity assumptions, it can be shown that the distribution flow still converges towards the target distribution $`\pi `$. Moreover, it converges at a rate only dependent on the mixing properties of $`L`$. This is in contrast with the MH algorithm whose rate of convergence is dependent on $`\pi `$ and $`M`$.
* The algorithm we have presented can be used to simulate a Markov chain with a fixed terminal point. Indeed, one obtains at time $`n`$ samples from (1). By setting $`L_0(x_1,dx_0)=\delta _x\left(dx_0\right)`$ and reversing the time index, one obtains an approximate realization of a Markov process of initial distribution $`\pi _n`$ at time $`0`$, transition $`\left\{L_{n\text{ }}\right\}`$ and terminal point $`x`$ at time $`n+1`$. This has applications in genetics and physics.

There are also several important open methodological and theoretical problems to study.

* Similarly to MCMC methods, one needs to carefully design the various components of the algorithm to get good performance. In particular, it would be of interest to come up with an automated choice for $`L_n`$ given $`M_n`$. For the homogeneous case, one could look at minimizing the variance of (3). It involves a tradeoff between the mixing properties of $`L`$ and the variance of the importance weights (2). This point is currently under study.
* It would be interesting to weaken the assumptions of the results in (Del Moral and Miclo, 2000; Del Moral and Doucet, 2003) which mostly only hold for compact spaces.

## 6 Acknowledgments

The authors are extremely grateful to Manuel Davy for his comments and some of the numerical simulations. We are also grateful to Christophe Andrieu and Elena Punskaya for their comments.


text:
# 1 Calculated values of sin²{2𝜃_{𝜇𝜏}} versus 𝑚_𝑡. The parameter region of 𝑚_𝑡 at GUT scale is within the allowed uncertainty. The horizontal line indecates the experimental lower bound for atmospheric neutrinos.

## Acknowledgements

This work started from the discussion at the research meeting held in Nov. 2002 supported by the Grant-in Aid for Scientific Research No. 09640375. We would like to thank to A. Sugamoto and T. Kugo whose stimulating discussion encouraged us very much. Also we are stimulated by the fruitful and instructive discussions during the Summer Institute 2002 held at Fuji-Yoshida. M. B. is supported in part by the Grant-in-Aid for Scientific Research Nos. 12640295 from Japan Society for the Promotion of Science, and Grants-in-Aid for Scientific Research on Priority Area A “Neutrinos” (Y. Suzuki) Nos. 12047225, from the Ministry of Education, Science, Sports and Culture, Japan. Also we are thankful to T. Takeuchi for his kind comment on our poor English.


text:
# Signatures of chaotic tunnelling

## I Introduction

When studying tunnelling in non-separable systems with more than one degree of freedom, one immediately encounters difficulties which generically can be traced back to the absence of enough constants of motion. Even in the very peculiar case of integrable systems, where continuous symmetries provide as many constants of motion as degrees of freedom, as soon as separability is lost, the analysis of tunnelling is not a simple generalization of what occurs in one-dimensional (1D) autonomous systems. The later case is detailed in textbooks on quantum physics (see for instance Messiah65a ) and it has even been possible to give a comprehensive analytical treatment in term of complex solutions of the Hamilton equations Balian/Bloch74a . However, it is not until the mid-eighties that a satisfactory quantitative approach has been proposed Wilkinson86b ; Wilkinson/Hannay87a  for tunnelling in non-separable integrable systems involving a larger number of dimensions. Moreover, integrability is a property of higher dimensional systems which is not generic. The coupling between several internal degrees of freedom as well as the coupling to an external source usually destroys some global constants of motion. With such a lack of constraints on the dynamics, the classical motion in phase space may become chaotic: it may explore volumes with higher dimensionality and therefore becomes exponentially sensitive on initial conditions. It is not surprising that these deep qualitative differences between an integrable regime and a chaotic one appear at the quantum level too. Some of the properties of a *quantum* system do change when constants of motion are broken and, indeed, it is the very object of quantum chaos to study the signatures of classical chaos at the quantum level (see for instance Giannoni+91a  to see how rich, vivid and successful this domain is).

We define tunnelling as a quantum process which is classically forbidden to *real* solutions of classical equations of motion. In this paper, we consider Hamiltonian systems only and study how the non-dissipative breakdown of continuous symmetries affects tunnelling. We do not consider how tunnelling is modified by dissipation and decoherence of the quantum wave. Of course, this requires a great care in real experiments where making dissipation negligible is always a hard task. This is one of the main reasons why very few real experiments have been done on these questions, though they would definitely help to understand tunnelling in the presence of chaos (as far as we know, the only experiments explicitly made on chaotic tunnelling in the XXth century are those presented in Dembowski+00a  with electromagnetic microwaves instead of quantum waves).

During the last fifteen years however, theoretical and numerical investigations on autonomous 2D and time-dependent 1D Hamiltonians systems highlighted some mechanisms Wilkinson86b ; Wilkinson/Hannay87a ; Lin/Ballentine90a ; Bohigas+93a ; Bohigas+93b ; Tomsovic/Ullmo94a ; Creagh/Whelan96a ; Brodier+02a  and much insight has been gained on the influence of non-separable dynamics. Experimental evidence of such mechanisms would be of great interest especially in the light of the subtle interplay between interferences and disorder. These phenomena lie in the general context of wave transport in complex media where the role of disorder is played by the (deterministic) chaotic dynamics instead of having a statistical random origin. Of course, other important motivations can be found in the numerous domains where tunnelling plays a crucial role as a fundamental quantum process like ionization Zakrzewski+98a , absorption, nuclear radioactivity, molecular collisions, mesoscopic physics etc. More speculatively, studies on tunnelling in high dimensional Hamiltonian systems should provide us with a natural extension of the instanton techniques which deal with quantum field theories reducible to effective 1D autonomous Lagrangian systems.

In 2001, it has been shown both theoretically Mouchet+01a  and experimentally Hensinger+01a ; Steck+01a ; Hensinger02a ; Hensinger+02a  that atom cooling techniques Arimondo+92a  (and possibly molecular physics as well, where formally similar systems have been extensively studied Lawton/Child81a ; Sibert+82a ; Kellman/Lynch87a ; Rose/Kellman96a ) provide systems which fulfill all the severe requirements to study tunnelling in the presence of classical Hamiltonian chaos: accurate manipulation of internal and external degrees of freedom, precise control of dissipation and decoherence and on the preparation/detection set up. For a brief account intended for a large audience see Heller01a ; Mouchet/Ullmo01a ; GossLevi01a ). The present article has the ambition to mark out the future challenging experiments, that remain to be done for reasons which will hopefully appear clear in the following.

This paper is organised as follows. In section II we give a general and informal overview. In section III, we quickly recall the main theoretical apparatus that is needed in the following. We implicitly refer to Mouchet+01a  for details and demonstrations. In section IV, we comment the results of Hensinger+01a  and Steck+01a . In section V, we show in this context, with the help of numerical experiments, the very precise form taken by the phenomenon known as chaos assisted tunnelling. We explain why it has not been observed yet with real atoms and propose how to actually bring it to the fore. Before the concluding remarks in section VII, we give in section VI some more numerical results which illustrate how subtle the signatures of chaotic tunnelling can be.

## II Chaotic tunnelling

The simplest situation to illustrate tunnelling is probably the case of a particle placed in a 1D time-independent symmetric double-well potential (see Fig. 1 (a)).

Starting in one well with an energy that is below the maximum of the potential, a quantum particle can jump into the other well with a non zero probability though it is a forbidden classical process. In addition to the classical time scale $`\tau `$ given by the oscillating period *inside* one well, we therefore have a longer time scale, the tunnelling period $`T\tau `$ of the oscillations *between* the wells. In the eigenenergy spectrum, tunnelling appears as a quasi-degeneracy of the odd and even-symmetry states whose energies are both of the order of $`\mathrm{}/\tau `$ but differ by an exponential energy splitting

$$\mathrm{\Delta }ϵ=\frac{2\pi \mathrm{}}{T}\text{e}^{A/\mathrm{}}$$
(1)
where $`A`$ is a $`\mathrm{}`$-independent typical action and can be interpreted in terms of a unique complex classical trajectory under the barrier Messiah65a ; Heading62a .

In the following we generalize this elementary situation in two ways. First, unlike the parity in the previous example, we can deal with a symmetry which is not necessarily either a spatial one or a two-fold one. In other words, we can have any discrete symmetry group acting on the whole phase space as well as any $`N`$-fold symmetry which lead to bunches of $`N`$-uplets in the energy spectrum (or bands if $`N1`$). In the following we keep $`N=2`$ since we have a two-fold symmetry $`𝒯`$ actually playing the role of parity (see Fig. 1 (b)) and being somehow decoupled from the other discrete symmetries. The classical structure in phase space is globally invariant under $`𝒯`$ and the quantum eigenstates can be classified according to their symmetric or antisymmetric character under the unitary transformation which represents $`𝒯`$ in the Hilbert space of states. Because $`𝒯`$ acts in phase space, it is usually more complicated than a pure spatial transformation. Thus, the two regions of phase space connected by quantum tunnelling, but classically not connected, are in general not separated by a simple potential barrier, but by a more complicated dynamical barrier. In such a case, tunnelling is emphasized to be called “dynamical tunnelling” as suggested by Davis and Heller in Davis/Heller81a . It often happens that the classically unconnected region are associated to the same region of configuration space, with different momenta. A simple study of the density probability in configuration space is then unsufficient to characterize dynamical tunnelling; an analysis of the density probablility in momentum space is required.

The second kind of generalization leads to much more puzzling questions. When dealing with systems with several degrees of freedom or, equivalently, if an external time dependence exists, classical trajectories generically lose their regular behavior, cannot analytically be computed and are organized in a fractal hierarchy that is described by the KAM perturbative scenario. Recently, important progress has been achieved in the understanding of the continuation of these intricate structures in complex phase space and their role at the quantum level, see Kus+93a ; Leboeuf/Mouchet99a  and especially Shudo/Ikeda98a . Anyway, we are therefore led to the following typical quantum chaos question: if one is able to create two symmetric stable islands separated in classical phase space by a chaotic sea whose volume is under control (see Fig. 2), what is the effect of this sea on the (dynamical) tunnelling between the islands?

The “dual” situation where chaos is created *inside* the wells while the dynamical barrier is kept regular has been introduced and studied theoretically and numerically in Creagh/Whelan96a . For a better understanding of what occurs in the energy spectrum when regular wells are separated by a chaotic sea, it has been proposed in Tomsovic98b  to slightly break the tunnelling symmetry. Nevertheless, in the following of the present paper, it must be kept in mind that a discrete symmetry will be always maintained exactly. At last, a third kind of generalization where the Hamiltonian character is destroyed by introducing dissipation and/or coupling to a thermal bath, is beyond the scope of this work Tomsovic98a .

## III Effective Hamiltonian system

Following Mouchet+01a ; Hensinger+01a ; Steck+01a  (see also Averbukh+95a  in a different context) we deal with an effective 1D time-dependent system whose Hamiltonian is

$$H(p,q;t)=\frac{p^2}{2}\gamma (\theta +\mathrm{cos}t)\mathrm{cos}q$$
(2)
in dimensionless units. $`\gamma `$ and $`\theta `$ are two classical real parameters that can be modified in real experiments. In addition, there is also one parameter, namely $`\mathrm{}_{\mathrm{eff}}`$, which fixes the quantum scale and is defined by the usual relation between canonical operators: $`[q,p]=\mathrm{i}\mathrm{}_{\mathrm{eff}}`$. It turns out that $`\mathrm{}_{\mathrm{eff}}`$ is not constant any longer (see section IV below). It can be also experimentally varied via the rescaling factor that is needed in the canonical commutation relation in order to work in dimensionless units used to write (2).

The time dependence breaks the conservation of energy and therefore may generate chaos. In order to deal with such a Hamiltonian, it is crucial to remark that it has both a spatial and a temporal periodicity. The latter implies that the Floquet theorem can be used, which states that the Hilbert space is spanned by an orthonormal eigenbasis of the evolution operator over one period. The corresponding eigenvalues of this unitary operator are distributed on the unit circle and therefore are labelled by their phase, which is conveniently written as $`\mathrm{exp}(\mathrm{i2}\pi ϵ/\mathrm{}_{\mathrm{eff}}),`$ where $`2\pi `$ stands for the period of the modulation and $`ϵ`$ can be interpreted as a quasi-energy, a generalization of the notion of energy level for a time-periodic system.

The spatial periodicity of the Hamiltonian is also extremely important, as it makes it possible to split the Hilbert space into independant components, each component being characterized by the so-called Bloch vector $`k`$ in the $`]0.5,0.5]`$ range: under translation of $`2\pi `$ along $`q,`$ the quasi-energy eigenstates are just multiplied by the phase factor $`\mathrm{exp}(\mathrm{i2}\pi k).`$ One is thus reduced to solve the Floquet-Schrödinger equation in an elementary spatial cell with boundary conditions depending on $`k.`$ One thus generates – for a fixed value of $`k`$ – a discrete quasi-energy spectrum $`ϵ_i(k).`$ When the full range of $`k`$ values is considered, one obtains the familiar (quasi-)energy bands Ashcroft/Mermin76a .

There is an additional discrete symmetry which can be used. The Hamiltonian (2) is invariant under the time-reversal symmetry $`(q,p,t)(q,p,t).`$ In the classical surfaces of section, this implies a symmetry with respect to the $`q`$ axis. In situations like the one in fig. 1b, this implies the existence of pairs of symmetric classically unconnected tori, i.e. a situation where tunnelling could be observed. In the quantum world, the situation is slightly more complicated, because this symmetry connects the $`k`$ subspace to the $`k`$ subspace. In the particular case $`k=0`$ ($`k=0.5`$ could also be used), this implies that the Floquet eigenstates can be split in two subclasses of states which are either even or odd under the symmetry operation. The splitting between a doublet of even and odd states, $`\mathrm{\Delta }ϵ_n=|ϵ_n^+(0)ϵ_n^{}(0)|`$ will be a measure of tunnelling.

We will extensively use the Husimi representation of a quantum state Husimi40a . Such a representation associates to each quantum state $`|\psi `$ a phase space function $`\psi ^\mathrm{H}(p,q)`$ (where $`p`$ and $`q`$ are real numbers) defined by

$$\psi ^\mathrm{H}(p,q)=|z|\psi |^2$$
(3)
where $`|z`$ is the coherent state corresponding to the complex number $`z=(q+\mathrm{i}p)/\sqrt{2\mathrm{}_{\mathrm{eff}}}`$. Because $`|z`$ is a minimal Gaussian wave packet with average momentum $`p`$ and average position $`q`$, the Husimi function $`\psi ^\mathrm{H}(p,q)`$ contains some information about the degree of localization of $`|\psi `$ in phase space, and makes it possible to associate quantum states with classical phase space structures.

## IV Experiments with cold atoms

Under some severe conditions which constrain the experiments, Hamiltonian (2) can be obtained as an effective dimensionless Hamiltonian for cold neutral independent atoms of mass $`M`$ interacting with two counterpropagating laser beams Mouchet+01a ; Hensinger+01a ; Steck+01a . These two beams have two slightly different frequencies at $`\omega _L+\delta \omega /2`$ and $`\omega _L\delta \omega /2`$. The longitudinal coordinate $`x`$ and the dimensionless $`q`$ are related by $`q=2k_Lx`$ where $`k_L=\omega _L/c`$. The rescaling of the momentum is given by $`p=(2k_L/M\delta \omega )p_x`$. $`\gamma `$ and $`\theta `$ are fixed by the intensity of the lasers and the detuning of the laser frequencies with respect to the atomic resonance. The dimensionless time $`t`$ is taken in $`\delta \omega ^1`$ units and the expression of the effective Planck constant is $`\mathrm{}_{\mathrm{eff}}=8\omega _R/\delta \omega `$ where $`\omega _R=\mathrm{}k_L^2/2M`$. Since in fine we want to measure exponentially small tunnelling splittings $`\mathrm{\Delta }ϵ`$, it requires to maintain these conditions for a time at least larger than $`\mathrm{}_{\mathrm{eff}}/\mathrm{\Delta }ϵ`$. Moreover, a very accurate control of the preparation of the initial state and of the analysis of the final state is compulsory.

As shown above, because of the temporal and spatial periodicity of the Hamiltonian, observing the standard signature of tunnelling, that is an oscillation of a quantum state between two classically unconnected regions of phase space requires that a single doublet of Floquet-Bloch eigenstates is initially populated, with well defined values of the parameters $`(\gamma ,\theta ,\mathrm{}_{\mathrm{eff}})`$ but also with a well defined value of the Bloch angle $`k`$. If more than a single doublet is populated, additional frequencies (related to energy differences between the various populated Floquet states) will appear in the temporal evolution. If any parameter is not fixed, the experimental signal will be the superposition of tunnelling oscillations (with different frequencies) for various sets of parameters. This will at best – if the dispersion of the parameter values is reasonnably small – blur the oscillations at long time and at worst will completely destroy the signature of dynamical tunnelling. It is experimentally rather simple to keep an accurate time-periodicity of the driving signal, i.e. to fix $`\mathrm{}_{\mathrm{eff}}.`$ Similarly, the balance between the constant and the oscillatory term, hence the parameter $`\theta `$ is easily controlled. The $`\gamma `$ parameter is proportional to the laser intensity and may thus slightly vary across the atomic cloud (because of the transverse structure of the laser beams). The most difficult part is to control that a single Bloch angle $`k`$ is excited. Indeed, this requires a phase coherence of the initial wavefunction over a large number of laser wavelengths, which is extremely difficult to achieve experimentally Greiner+02a , as will be shown in the following. In any case, the inhomogeneous broadening of the experimental signal because of the dispersion in $`k`$ will be responsible for a decay of the tunnelling oscillations.

### IV.1 NIST experiments Hensinger+01a

In the nist experiments, the two stable symmetric islands are chosen quite close in phase space in order to deal with not too small splittings. Another crucial point of this experiment is that the classical motion of the islands over one period, unlike those in Mouchet+01a  and Steck+01a , always remains trapped in one spatial elementary cell of length $`2\pi `$. The quantum states localized in these islands are consequently only weakly sensitive on the boundary conditions which are governed by the Bloch angle. In other words, the tunnelling period will be only weakly dependent on the Bloch angle $`k`$, which implies that the unavoidable broadening over $`k`$ will not spoil too much the signature of tunnelling. This is a major improvement over the tunnelling described in Mouchet+01a  and Steck+01a , where a very narrow band of Bloch angle is required to observe clear tunnelling oscillations. Moreover, the atoms involved in the tunnelling process stay longer in the region where the laser intensities are uniform.

Indeed, as proposed in Hensinger+01b , (Hensinger02a, , chap. 4 and 5) the two stable symmetric islands are created from a pitchfork bifurcation of the fixed point at $`(p,q)=(0,0)`$. To visualize it (see Fig. 2 (a) and (b)), we extract a one-parameter sequence by varying $`\gamma `$ while $`\theta `$ is fixed to the experimentally chosen value in Hensinger+01a , i.e. $`\theta =1.7`$. When $`\gamma `$ is increased, the pairs of symmetric tori appear at $`\gamma 0.56`$. At the center of each set ot tori, there is a periodic orbit. Over one period of the driving, the periodic orbit is essentially a rotation over the fixed point at $`(p,q)=(0,0),`$ which explains that the whole structure remains trapped in a single spatial cell. For $`0.56\gamma \gamma _c,`$ the tori remain nested in one connected stable island. At $`\gamma =\gamma _c0.625,`$ a chaotic sea separates the symmetric islands which shrink and move away from the central point before being dissolved through a cascade of bifurcations starting at $`\gamma 0.97`$.

In one series of experiments, $`\gamma 0.96`$ and $`\mathrm{}_{\mathrm{eff}}0.8`$, the atoms are prepared in one island and their average momentum $`p`$ is measured stroboscopically every modulation period ($`=2\pi `$ in dimensionless units). Since, in phase space, the islands rotate about the origin with the same period, if no tunnelling occurred no variation in $`p`$ would be noticeable. In fact, starting the measurement sequence when $`p`$ has its maximum value, oscillations are observed which illustrate the back and forth motion of the atoms between the islands due to dynamical tunnelling. The tunnelling period $`T`$ is about 10 modulation periods in this case ($`200\mu \mathrm{s}`$). This is in perfect agreement with the quasi-energy splitting obtained numerically for the two Floquet eigenstates having the largest Husimi functions inside the islands.

It is worth noting that the nist group uses a Bose-Einstein condensate as a preliminary step for preparing atoms in well defined quantum states, especially for achieving a large coherence length for the wavefunction, i.e. a small spreading of the Bloch angle $`k.`$ In order to prepare phase space localized states, an optical lattice is carefully turned on. When the tunnelling experiment starts, the atomic density and the interaction between atoms is sufficiently small, and the experiment can be analyzed as the interaction of individual independent atoms with the laser beams, i.e. using Hamiltonian (2). However, the cloud of atoms remains cold enough, at a sub-recoil temperature, to prevent a large thermal broadening of momentum distribution that would destroy the signal. Because they start from very low temperature, these preparation techniques based on condensate manipulation seem to allow a wider room to manoeuvre than those working with thermal clouds only. Adiabatic switching of the light potentials is not required and one can actually work with values of the classical parameters $`\gamma `$ and $`\theta `$ which are far from the perturbative regime of an integrable system.

By diagonalizing the evolution operator corresponding to eq. (2) over one period, we are not only able to reproduce the oscillatory behavior of $`p(t)`$ (see Fig. 3 (a)),

but also we can study the spoiling effect of the thermal dispersion $`\mathrm{\Delta }p\sqrt{\mathrm{temperature}}`$ and predict the maximum allowed temperature (see Fig. 3 (b) and (c)). 111After the first version of the present paper was written, we learnt that an independent numerical work Luter/Reichl02a  had obtained the same results as the ones in our Fig. 3 (a). But Luter/Reichl02a  do not consider the thermal effects and work always within the Floquet theory at $`k=0`$ instead of the Floquet-Bloch theory. In the present work we clearly demonstrate that thermal effects cannot be neglected and must be studied carefully when experiments are discussed. If $`\alpha `$ denotes the width of the momentum distribution in recoil momentum units, it can be shown (Mouchet+01a, , §6.a) that it corresponds to a statistical mixture of Bloch states with $`\mathrm{\Delta }k=\alpha /2`$. Fig. 3 (a) corresponds to the ideal situation where all atoms are prepared with $`\alpha 1`$ about the $`k=0`$. When a small but non vanishing $`\alpha `$ is introduced, some states of the quasi-energy bands with non vanishing $`k`$ get involved and blur the tunnelling oscillations. For $`\alpha =0.2`$, the oscillation amplitude is reduced by a factor 2 and for $`\alpha =0.4`$ have nearly disappeared. Therefore, in this experiment, having a sub-recoil atom cloud is required.

In the following we want to focus on tunnelling only and we will implicitly keep $`k=0`$.

### IV.2 Austin experiments Steck+01a

For a better understanding of the dynamics, one must go beyond the two-level model involving the symmetric and the antisymmetric states only. Other states must be taken into account and their influence can be enlightened when a classical parameter ($`\gamma `$ or $`\theta `$) or the quantum one ($`\mathrm{}_{\mathrm{eff}}`$) is continuously varied. Two (quasi-)energies may exactly get degenerate if they belong to distinct symmetry classes. If not, they may follow a so-called avoided crossing whose size reflects the direct coupling between the two states (more precisely the off-diagonal matrix element of the coupling perturbation), but also the indirect coupling with other states. One of the keys of the chaotic tunnelling problems is to identify clearly the qualitative nature and the quantitative influence of indirect coupling. This is the background of Austin experiments.

A third level is involved in a non negligible indirect coupling when its quasi-energy approaches the tunnelling doublet energies. This can be understood from perturbation theory as the leading term of the indirect coupling is proportionnal to the inverse of the energy difference. The two generic scenarii of the crossing of the doublet by a third state are shown in Fig. 4.

Aside from the unavoidable ambiguous definition of the splitting, it appears clearly that case a) corresponds to an increase of the tunnelling splitting during the crossing while, conversely, case b) can lead to arbitrary small splittings since an exact degeneracy occurs. Thus, such a crossing by a third state produces a sharp variation of the tunnelling period that can be measured experimentally. This is actually what is observed in the Austin experiments and can be confirmed by numerical experiments as shown in Fig. 5.

When looking at the Husimi representation of the states, Fig. 6,

one can immediately distinguish between the tunnelling doublet and the third state sufficiently far from the crossing. As expected, the tunnelling doublet has Husimi representations localized in the stable islands though they also spread in the chaotic sea. On a classical Poincaré surface of section, it is easy to make the difference between chaotic and regular motion; how to transpose this distinction at the quantum level is not known with the large values of $`\mathrm{}_{\mathrm{eff}}`$ used in both the Austin and the Nist experiments. Some classical structures much smaller than the de Broglie wavelength are possibly present in some of the states 222From time to time, it is claimed Zurek01a  that a simple matching between classical structures and quantum wave sub-Planckian structures was found but, as it was understood long time ago by Balazs/Voros90a , the latter are generally washed out as soon as one try to measure global averages., but just looking at the Husimi representation of the third state for $`\mathrm{}_{\mathrm{eff}}1`$ does not make it possible to attribute any chaotic or regular character to it. It is only for much smaller $`\mathrm{}_{\mathrm{eff}}`$ that chaotic or regular wavefunctions make sense. This is not surprising: the dichotomy between regular and irregular motion is a classical one and, at present, it can be extrapolated at the quantum level within the semiclassical regime only. Anyway, one must keep in mind that tunnelling does only make sense in the semiclassical regime as well. One of the great merit of Austin experiments is to show for the first time a quantum tunnelling effect where an indirect process is involved, but we fill exagerated to attribute any chaotic origin to it.

## V Numerical evidence of a chaotic tunnelling regime

It is one of the first success of quantum chaos to have shown that the energy levels of an integrable system ignore each others because they are localized on different classical tori, while in chaotic systems level repulsion is the rule. Following the discussion in the previous section, one may therefore expect that the average size of avoided crossings is increased when chaos is present. The fluctuations of a tunnelling process which should be narrow and sparse in an integrable regime should be broader, more numerous and possibly involving many states in a chaotic regime. We now illustrate this statement in the framework of the experimental atomic Hamiltonian (2).

There are two different ways of rendering chaos observable by quantum eyes: the first one consists in increasing the volume of the chaotic sea with the help of a classical parameter, the other one fixes the classical dynamics and decreases $`\mathrm{}_{\mathrm{eff}}`$. We will present both ways.

CAT, RAT and all that… One may try to study the tunnelling *fluctuations* separately from the *average behavior*. This average — in a somewhat vague sense — is increased because chaos diminishes the classical dynamical barrier333The question of defining an integrable system with respect to which the chaotic one should be compared is extremely difficult because of the exponential sensitivity of tunnelling on classical parameters. It is a much more serious problem than the way of defining an averaging procedure. and this is the reason why the phenomenon can be called chaos assisted tunnelling. However, as far as only fluctuations are concerned (there can be an enhancement or a decrease as well), the words “chaos assisted tunnelling” (CAT) Tomsovic/Ullmo94a  may lead to confusion and we simply use “chaotic tunnelling”. At last, Brodier, Schlagheck and Ullmo discovered what they called “resonant assisted tunnelling” (RAT) Brodier+01a  to describe an enhancement of tunnelling due to an indirect process which involves one or several quasi-modes localized in the secondary resonances surrounding the main symmetric islands. Up to now, RAT has been studied quantitatively in a purely quasi-integrable case but there are clues that it can be extended far beyond kam theory.

### V.1 $`\gamma `$ change

Let us first take $`\mathrm{}_{\mathrm{eff}}`$ ten times smaller than in the Nist and Austin experiments. It is much easier to do it numerically than experimentally, as it requires to increase the modulation frequency of the laser beams by one order of magnitude. We then follow the quantum states through the classical bifurcation shown in picture 2 and discussed in section IV.1. Again, in order to calculate the tunnelling splitting $`\mathrm{\Delta }ϵ`$, we select the states that have the largest Husimi functions inside the islands.

Figure 7 shows that, after a smooth decrease of $`\mathrm{\Delta }ϵ`$ up to $`\gamma =\gamma _c0.625`$ there is an abrupt change of regime. First, the mean value of $`\mathrm{\Delta }ϵ`$ increases and second, many fluctuations appear which modify the splitting by several orders of magnitude. It is remarkable that this change of regime can be matched on the Poincaré surfaces of section. The smooth tunnelling regime occurs the resonant tori belong to one quasi-integrable island. For $`\gamma <\gamma _c,`$ we are able to reproduce the main features of tunnelling by using an integrable approximation (see appendix). The decrease of tunnelling is directly understood in terms of the lengthening of the dynamical barrier. $`\gamma _c`$ corresponds exactly to the point where the islands get disconnected. In addition, one can follow the bifurcation on the Husimi representations of the tunnelling doublet.

One can hardly detect by eye any regularity in the chaotic regime but four large spikes in the range $`\gamma _c<\gamma <0.7`$. This can be traced back to the crossing by the same third state whose quasi-energy line is folded four times in the Floquet zone centered on the doublet.

### V.2 $`\mathrm{}_{\mathrm{eff}}`$ change

Figure 8 shows $`\mathrm{\Delta }ϵ`$ as a function of $`1/\mathrm{}_{\mathrm{eff}}`$ for four values of $`\gamma `$. Here again, two types of regimes, a quasi-regular and a chaotic one, can clearly be distinguished. For the values of $`\gamma `$ where some substantial chaos is present in between the islands, by decreasing $`\mathrm{}_{\mathrm{eff}},`$ one gets to the chaotic regime. For $`\gamma <\gamma _c`$, i.e. in the quasi-regular regime, one may enter in a chaotic regime but for far much low value of $`\mathrm{}_{\mathrm{eff}}`$ since the chaotic layers between kam tori are so thin that they are not even resolved in the Poincaré surface of section given in the figures.

On this plots, a purely exponential law given by eq. (1) would produce a straight line with negative slope. It gives rather poor predictions in the chaotic regime and should be corrected even in the quasi-integrable regime in order to reproduce the fluctuations. Whenever such a fluctuation is due to a crossing by a quasi-mode (and not a state delocalized in the surrounding chaotic sea), it might be reproduced by the resonance assisted tunnelling techniques.

## VI Statistical signature of chaotic tunnelling ?

In order to reproduce quantitatively the statistics of the tunnelling splitting fluctuations in the chaotic regime, Leyvraz and Ullmo Leyvraz/Ullmo96a  have introduced a random matrix model. The Hamiltonian can be split in two uncoupled components associated with the two even and odd symmetry subspaces. The corresponding matrices are written as:

$$H^{\mathrm{even}}=\left(\begin{array}{cccc}ϵ_0^+& v_1^+& v_2^+& \mathrm{}\\ v_1^+& & & \\ v_2^+& & H_{}^+& \\ \mathrm{}& & & \end{array}\right)$$
(4)
and

$$H^{\mathrm{odd}}=\left(\begin{array}{cccc}ϵ_0^{}& v_1^{}& v_2^{}& \mathrm{}\\ v_1^{}& & & \\ v_2^{}& & H_{}^{}& \\ \mathrm{}& & & \end{array}\right)$$
(5)
where $`ϵ_0^\pm `$ represent the energies of the doublet, $`H_{}^\pm `$ the Hamiltonian in the chaotic sea (modelled by a random Gaussian matrix) and $`v`$ the indirect coupling.

Neglecting direct tunnelling consists in taking $`ϵ_0^+=ϵ_0^{}`$. The central hypothesis is to consider all the $`v`$’s as independent variables with *the same* Gaussian distribution. This is quite natural to treat all the other states on the same footing, as they are assumed to be chaotic states randomly delocalized in the chaotic sea. With these assumptions, the splitting distribution can be calculated and is given by a (truncated) Cauchy distribution, see Leyvraz/Ullmo96a . For Hamiltonian (2), in each chaotic case where it has been tested, the Leyvraz-Ullmo prediction is in agreement with the numerical results (see Mouchet+01a ). More surprisingly, we have found that the Leyvraz-Ullmo law gives correct predictions even when the classical dynamics is quasi-integrable (see Fig. 9 and the corresponding Poincaré surface of section in Fig. 1) and the de Broglie wavelength much larger than the chaotic layers.

When looking at the splitting, one observes numerous and large fluctuations (over several orders of magnitude) which were supposed to characterize the chaotic regime. Conversely, we have not been able to find a regular regime of fluctuations with a classically chaotic dynamics. Therefore one must conclude that even though two different regimes of tunnelling fluctuations can be identified unambiguously, classical chaos appears to be a sufficient but not a necessary condition for having numerous and large fluctuations governed by the Leyvraz-Ullmo law.

These unexpected results are not explained at the present state of the theoretical approaches of chaotic tunnelling. If it appeared, within future numerical or real experiments, that these results are not due to the peculiar properties of our system, it would definitely mean that theoretical studies should be all the more needed.

## VII Conclusion

In this paper we have studied in great details the transition between a regular and a chaotic regime of tunnelling within a classical configuration that can be achieved experimentally. We have shown why $`\mathrm{}_{\mathrm{eff}}`$ has to be small enough if one wants to reach for the first time the chaotic regime in real systems. In recent experiments with cold atoms, it requires to increase the modulation period up to one order of magnitude (that is at least the MHz).

Theoretically, there is also a lot of work to do if we want to understand and therefore predict the fluctuations quantitatively. The semiclassical regime requires to study carefully the dynamics with complex coordinates. However, the present work has shown clearly that the abrupt transition between a regular and a chaotic regime of tunnelling corresponds to a classical transition that can be identified very precisely. Tunnelling being a relevant concept in a semiclassical regime only, it is therefore not surprising that the future investigations on chaotic tunnelling will have to keep track of the classical dynamics in one way or another.

###### Acknowledgements.

We acknowledge A. Shudo, S. Tomsovic, D. Ullmo and W. Hensinger for stimulating discussions. A. M. is grateful towards O. Boebion for computer assistance in the Laboratoire de mathématiques et de physique théorique of Tours and thanks the Laboratoire Kastler Brossel of Paris for kind hospitality. Laboratoire Kastler Brossel de l’Université Pierre et Marie Curie et de l’Ecole Normale Supérieure is UMR 8552 du CNRS. Laboratoire de mathématiques et de physique théorique de l’Université François Rabelais is UMR 6083 du CNRS.

## Appendix A Integral approximation near the pitchfork bifurcation

When one classical parameter is smoothly varied in a non-integrable Hamiltonian system, its periodic orbits follow infinite fractal-like cascade of bifurcations through which they cannot be followed smoothly. However, all the bifurcations can be classified according to a simple set of scenarios: for a given bifurcation of a given periodic orbit, the phase-space dynamics of the original parameter-dependent Hamiltonian can be uniformly approximated by an integrable parameter-dependent Hamiltonian which retains the relevant features only. Of course, because one cannot get rid of the chaotic dynamics, this approximation does make sense only locally, that is near the periodic orbit, in the neighborhood of the bifurcation. It is the object of the Hamiltonian normal form theory to classify the bifurcations and obtain the simplest form of the approximated integrable Hamiltonians (the so-called normal forms). Generically, i.e. when no constraint or symmetry is present, the one parameter Hamiltonian normal forms have been completely classified by Meyer Meyer70a ; Meyer/Hall92a ; Leboeuf/Mouchet99a . However, the bifurcation shown in Fig. 2 is out of the scope of Meyer’s classification precisely because the time reversal symmetry plays a key role. Even if the Poincaré surface of section near the origin (see Fig. 2 (b) cannot be distinguished from the one corresponding to Meyer’s transitional bifurcation (see for instance Fig. 9 (b) in Leboeuf/Mouchet99a ), it is crucial to note that, in our case, two *distinct* $`2\pi `$-periodic orbits have emerged from the origin. In a transitional bifurcation, the two stable islands would correspond to the *same* $`4\pi `$-periodic orbit. Therefore they would be classically connected to each other and would be irrelevant for tunnelling. Let us sketch briefly how to obtain the Hamiltonian normal form in our case.

1. The first step is to find the value $`\gamma _0`$ of $`\gamma `$ at which the bifurcation occurs ($`\theta `$ is kept fixed). Such a bifurcation occur when the trace of the monodromy matrix $`M_\gamma `$ at the origin after one period is 2 and corresponds in the $`(2\gamma ,4\gamma \theta )`$-plane to the border of the even Arnold tongues (see Fig. 20.1 in Abramowitz/Segun65a ) (the transitional case occurs when $`\mathrm{tr}M_\gamma =2`$ at the border of the odd Arnold tongues). More precisely, the case shown in Fig. 2 corresponds to $`\theta =1.724137`$ (the experimental value in the Nist experiment) and the bifurcation takes place at $`\gamma _00.564673`$. In the following, we use $`\epsilon =\gamma \gamma _0`$.

Let us denote by $`𝒴_\epsilon (z)`$ \[resp. $`𝒵_\epsilon (z)`$\] the solution of the Mathieu equation

$$y^{\prime \prime }(x)+(4\gamma \theta +4\gamma \mathrm{cos}(2x))y(x)=0$$
(6)
such that $`𝒴_\epsilon (0)=0`$ and $`𝒴_\epsilon ^{}(0)=1`$ \[resp. $`𝒵_\epsilon (0)=1`$ and $`𝒵_\epsilon ^{}(0)=0`$\]. The prime stands for the derivative with respect to $`x`$. For $`\epsilon =0`$, it is straightforward to show that

$$M_{\gamma _0}=\left(\begin{array}{cc}1& \frac{1}{4\pi }𝒵_0^{}(\pi )\\ 0& 1\end{array}\right)$$
(7)
where $`𝒵_0^{}(\pi )1.480919`$ for $`\theta =1.724137`$.

2. The second step is to make a (linear) $`2\pi `$-periodic canonical change of coordinate that eliminates the time-dependence in the quadratic part of Hamiltonian (2), near the origin and uniformly in $`\epsilon `$. We are then led to the Hamiltonian

$$\begin{array}{c}\hfill \left(\frac{1}{2}\frac{1}{4\pi }𝒵_0^{}(\pi )+\alpha \epsilon \right)q^2+\frac{1}{2}\beta \epsilon p^2+\frac{1}{2}\delta \epsilon p^2\\ \hfill +\text{higher order }2\pi \text{-periodic terms}\end{array}$$
(8)
where $`\alpha ,\beta ,\delta `$ are $`\epsilon `$-independent coefficients. Only $`\beta =\frac{1}{\pi }_\epsilon 𝒴_\epsilon (\pi )`$ evaluated at $`\epsilon =0`$ will be relevant since $`\alpha `$ and $`\delta `$ can be eliminated by a suitable canonical change of coordinates following a method explained in (Leboeuf/Mouchet99a, , section 4.2). In our case $`\beta 2.008/\pi `$. We are therefore led to the following normal form of the quadratic part of the Hamiltonian:

$$\begin{array}{c}\hfill \left(\frac{1}{2}\frac{1}{4\pi }𝒵_0^{}(\pi )\right)q^2+\frac{1}{2}\beta \epsilon p^2\\ \hfill +\text{higher order }2\pi \text{-periodic terms}\end{array}$$
(9)
3. Following the same reasoning that lead to the transitional normal form (Leboeuf/Mouchet99a, , sections and 4.2), all higher order $`2\pi `$-periodic terms, but the resonant terms of the form $`h_k(\epsilon )p^k`$, can be cancelled by a suitable canonical change of coordinates. Because of the time reversal symmetry, the coefficient $`h_3`$ vanishes identically and therefore the leading order normal form is

$$\left(\frac{1}{2}\frac{1}{4\pi }𝒵_0^{}(\pi )\right)q^2+\frac{1}{2}\beta \epsilon p^2\frac{1}{4}h_4(0)p^4.$$
(10)
The explicit calculation of $`h_4(0)`$ is tedious but it can be estimated numerically by fitting the coordinate $`q=0,p=\pm \sqrt{\beta \epsilon /h_4}`$ of the satellite $`2\pi `$-periodic orbits for $`\epsilon >0`$. We obtain $`h_40.0320`$.

It is far from obvious that the quantization of the normal form (10) will give a good approximation of the quasi-energies of the tunnelling doublet. Some discrepancies may rise from the fact that quantum physics is invariant under canonical transformations only at the leading order in $`\mathrm{}`$. When we swap $`p`$ with $`q`$ and change the sign of the energies, the normal form (10) leads to a standard 1D double-well problem whose quantum spectrum can be found by numerically diagonalizing the Hamiltonian written in a harmonic basis. The tunnelling splitting of the ground doublet is given by the dotted line in Fig. 7 and it agrees reasonnably well with the exact numerical result for $`\gamma `$ near $`\gamma _0`$.