Boltzmann equation with double-well potentials

Boltzmann equation with double-well potentials

Silvia Chiacchiera CFisUC, Department of Physics, University of Coimbra, P-3004-516 Coimbra, Portugal    Tommaso Macrì Departamento de Física Teórica e Experimental, Universidade Federal do Rio Grande do Norte, 59072-970 Natal-RN,Brazil International Institute of Physics, 59078-400 Natal-RN, Brazil    Andrea Trombettoni CNR-IOM DEMOCRITOS Simulation Center, Via Bonomea 265, I-34136 Trieste, Italy SISSA and INFN, Sezione di Trieste, Via Bonomea 265, I-34136 Trieste, Italy

We study the dynamics of an interacting classical gas trapped in a double-well potential at finite temperature. Two model potentials are considered: a cubic box with a square barrier in the middle, and a harmonic trap with a gaussian barrier along one direction. The study is performed using the Boltzmann equation, solved numerically via the test-particle method. We introduce and discuss a simple analytical model that allows to provide estimates of the relaxation time, which are compared with numerical results. Finally, we use our findings to make numerical and analytical predictions for the case of a fermionic mixture in the normal-fluid phase in a realistic double-well potential relevant for experiments with cold atoms.

51.10.+y, 02.70.Ns, 67.85.Lm

I Introduction

Double-well energy potentials, with two, degenerate or not, minima separated by a maximum, are ubiquitous in several branches of science. These potentials are used to model a variety of processes involving an energy barrier Hanggi90 (), from the computation of rate coefficients in a chemical reaction Calef83 () and the modeling of solid-state junctions Barone82 () to the calculations involving non perturbative instanton configurations in quantum field theory Rajaraman87 (). The importance per se of the study of the dynamics in double-wells potentials stems also from the fact that it is often preliminary to (and useful for) the investigation of more complex dynamical effects in multi-well potentials.

The prototypical problem of the dynamics in a double-well potential is to determine in how much time the particles move from one well to the other one and if (and in how much time) they equilibrate reaching a vanishing , where is the difference of population between the two wells. A source of inter-well motion is of course given by the quantum tunneling Barone82 (); Razavy03 (). When quantum effects are dominant, a single particle tunnels from one well to the other, and for many particles the macroscopic quantum coherence exhibited at low temperature by many systems – including He and He, ultracold atoms and superconductors Tilley (); Annett () – allows for a net current between the wells. At variance, at high temperatures thermal effects give rise to noise-assisted hopping events Hanggi90 (): for ultracold atoms, at temperatures higher than the temperature at which the effects of quantum statistics become relevant, this corresponds to an incoherent flow of atoms and the ensuing thermalization of the two-well systems with macri13 ().

Experiments with cold atoms in double-well potentials give the concrete possibility to explore physical situations in which both phenomena – quantum tunneling and thermal noise-assisted hoppings – are present. The high degree of control of atomic gases at very low temperatures lewensteinbook () allows to design and perform experiments where particles, either bosons, fermions or mixtures of both, are subjected to properly engineered and highly tunable trapping potentials. In this context, superfluid double- and multi-well dynamics have been extensively studied, both theoretically and experimentally, primarily for bosonic atoms and more recently also for fermionic gases. For bosons in the superfluid regime at , which are well described by the Gross-Pitaevskii equation, this has lead to the identification of the atomic analog of the Josephson effect of superconductivity and of the macroscopic quantum self-trapping effect Sme97 (); Alb05 (); Ank05 (), a direct consequence of the nonlinearity of the dynamical equations of motion. The study of superfluid fermions in double-well potentials is more recent and experiments for these systems are in progress. The peculiarity of such fermionic systems is, among others, the tunability of the strength of the inter-particle interactions via the so-called Feshbach resonances, which results in the well known BCS-BEC crossover Zwerger (). In a recent experiment Valtolina15 () Valtolina et al. studied ultracold fermionic Li atoms in two different hyperfine states loaded in double-well potentials, reporting on the observation of the Josephson effect between fermionic superfluids along the crossover.

Regarding thermal effects, another class of experiments focused during the years on the study of i) polarized and two-component fermionic gases across and above the transition from the superfluid to the normal state occurring at the critical temperature ; and ii) bosons at finite temperature, also near and above the Bose-Einstein condensation temperature ( coinciding with for Bose-Einstein condensates). The Boltzmann equation Landau10 () provides a major tool to describe the collective dynamics of cold atoms in the normal-fluid phase. From such studies it emerged that the description with the Boltzmann equation (including the quantum statistics modification of the collision term Haug08 ()) works rather well not only for fermions above the Fermi temperature and for bosons above , but also for two-components fermions at temperature smaller than and larger than Riedl2008 (); Chiacchiera2011 (); Pantel2012 (); Pantel2015 (). We also mention that the dynamics of bosons at finite temperature below has been studied resorting to a Boltzmann-like description of the thermal part of the gas ZNG () – a similar approach based on the Boltzmann description of the non-superfluid part of a two-component fermionic mixture has been discussed in Urban07 (). Most of these studies have been performed in single-well potentials without tunneling between wells, also given the difficulty of computing in a quantitatively reliable way the relaxation time which is determined at finite temperature by rare events of hopping across the wells. Studies of the dynamics of a Bose gas below in a periodic multi-well potential were presented in Nikuni (); Nikuni07 (); McK08 (), while the study of the fermionic transport in optical lattices at finite temperature (above ) was reported in Rosch (), where the Boltzmann equation was investigated in local relaxation time approximation. The center of mass oscillations of a normal Fermi gas in a one-dimensional periodic potential were studied both theoretically and experimentally Pez04 (); Orso04 ().

Although from one side a huge effort has been devoted to the study of quantum tunneling for superfluid atomic gases in double-well potentials at low temperatures, and from the other side the study of reaction-rate theory Hanggi90 () and of the Boltzmann equation are two workhorses of non-equilibrium physics, the study of the Boltzmann equation itself in a double-well potential is to the best of our knowledge a relatively not addressed topic. Motivated by systems of cold atoms in tunable geometries at finite (and possibly) variable temperature, in this paper we therefore study the double-well dynamics of an interacting gas at finite temperature within the framework of the Boltzmann equation with the classical collision term.

Our aim is to understand and describe, both qualitatively and quantitatively, the effect of two-body collisions on the the double-well dynamics. This problem is interesting in view of current and future experiments with cold atoms at finite temperature, since our approach can be applied to study the normal-fluid dynamics in double-well potentials. Here we examine how a gas in a symmetric double-well, prepared with an initial population imbalance, relaxes towards the balanced equilibrium situation and we propose a simple analytical model to describe our numerical findings. Performing a comparison between numerical and analytical results, we analyze the relaxation time as a function of barrier heights and interaction strengths. We observe that, even though we are not going to consider the effect of the quantum statistics on the collision term, we expect that both the numerical computations and the analytical model can be straightforwardly extended to include such correction and that the relaxation time would have a similar dependence on barrier heights and interaction strengths, with qualitative changes of the dependence of on the temperature only for, say, Riedl2008 (); Chiacchiera2011 () (and therefore close to in the unitary limit).

The article is organized as follows. In Sec. II we introduce the formalism briefly describing the numerical method used in the work. We also define the two model external potentials we consider. The first model is a square-well potential with a barrier of vanishing width in the middle that provides a simplified toy model for our numerical and analytical study. The second model, which is directly inspired by experimental work with cold atoms, is the superposition of a harmonic isotropic potential and a gaussian barrier. In Sec. III we introduce a simple analytical model for the study of the dynamics. In Sec. IV and Sec. V we present the numerical results for, respectively, the square double-well potential and the realistic one and compare them with the predictions of the analytical model. Finally, we draw in Sec. VI our conclusions and discuss possible improvements of the present paper for future work, while some useful results are collected in the Appendices.

Ii The Boltzmann equation with a double-well potential

In this Section we briefly recall the Boltzmann equation and introduce the two double-well potentials we consider in the following. Section II.1 is devoted to introduce the Boltzmann equation formalism and remind some of its basic properties used hereafter, while in Section II.2 we briefly sketch the procedure for the numerical solution of the Boltzmann equation based on the test-particle method. Section II.3 is devoted to introduce the two double-well potentials studied in this paper: a “toy” square double-well (SDW) with a filtering wall and a more realistic harmonic-gaussian double-well (HGDW) potential. In Section II.4 we define the various time scales appearing in the double-well problem.

ii.1 Boltzmann equation

We consider a one-component gas of classical interacting particles of mass , in an external potential . The evolution in time of the phase-space distribution function is governed by the Boltzmann equation Landau10 ()


The left-hand side, where and , represents transport. The right-hand side is the collision integral


where the differential cross section and the notation is used as a shortcut for the distribution function evaluated at the same , but with momenta , respectively. The normalization condition is , and the average of any one-body variable is , where the phase-space volume element is .

At equilibrium, the distribution function reads


where is the inverse temperature and the chemical potential. The total collision rate at equilibrium can be computed exactly using


where the factor is needed to avoid double counting (since we are dealing with a one-component gas). Expressions for the collision rate in a single square-well and in a harmonic potential are collected in Appendix A.

In this work, we consider particles interacting with a constant cross section, i.e., with no energy or momentum dependence:


ii.2 Test-particle method

A fully numerical approach to solve the Boltzmann equation is provided by the so-called test-particle method, in which the coordinates of all the particles in phase space are evolved individually. This method, similar to molecular dynamics but with a stochastic component, was developed in the ’80s in the context of nuclear physics Bertsch88 () and recently has been applied also to cold atoms Jackson02 (); Toschi2003 (); Lepers2010 (); Wade2011 (); Goulko2011 (); Pantel2015 (). Formally, the distribution function is discretized as a sum of delta-functions peaked at the position and momentum of each test particle:


where is the number of real particles and is the number of test particles, those entering the simulation. In some cases it is convenient to have : if the number of real particles is very low, one associates to each real particle many test-particles (), to describe the evolution in phase-space with a finer resolution. At variance, in the opposite case of too many real particles to be simulated individually, one test particle represents many real ones () Pantel2015 (). We observe that for a generic , the interactions between test-particles are ruled by a cross-section that is related to the real one as follows


In this work, we always take .

The average value of any one-body observable is obtained as follows


In absence of inter-particle interactions, the phase-space coordinates of each test particle are evolved according to the Hamilton equations


The actual numerical scheme to integrate them depends on the potential. In the case of a simple box, particles are propagated via the Euler method (exact in this case), and the collisions with the walls are implemented reversing the appropriate momentum component. For a generic potential, the particle position and momentum are evolved from the time step to via the velocity Verlet algorithm Allen ()


with the acceleration and an intermediate time step.

For interacting particles, collisions have to be implemented too. The cross section defines an interaction range : in each time step, all the pair of particles (within a certain distance) are checked, and they are collided if: 1) they reach their closest approach Lepers2010 () in the time step; and 2) their distance at closest approach is within the interaction range. A collision takes place randomly assigning new momenta to the participants, with the constraints of energy and momentum conservation. In our simulations the trajectories of the colliding particles are corrected to take into account the fact that the collision takes place at the closest approach point Lepers2010 (). For example, one can check that in the case of a harmonic well this correction is necessary to respect the balance of kinetic and potential energy.

ii.3 Two models for the double-well

We consider two model potentials for the double-well. The SDW potential has a rectangular energy barrier (of negligible width) located at the center of the system, which is in turn chosen to be a square well. This SDW has the advantage to be more simply numerically simulated using the test-particle method and it allows for an analytical treatment of the approximate model we are going to introduce in Section III for the determination of the relaxation times in double-well potentials. The other potential is relevant for cold-atom experiments and it is provided by the sum of a harmonic potential plus a barrier energy, chosen of gaussian form both for simplicity and also because such a barrier could be created by a superimposed blue-detuned potential. We refer to the second double-well potential as the harmonic-gaussian double-well (HGDW) potential.

The first model (SDW) is a toy model corresponding to a cubic box of side partitioned into two regions by a filtering wall: i.e., particles are allowed to pass through it or are reflected depending upon their momentum component orthogonal to the wall (say, ). When a particle during its propagation tries to cross the plane , we check for its momentum along : if , where and the height of the barrier, then the particle it is allowed to pass, as if the barrier were not there; if instead , it is reflected by the barrier. This model can be seen as a finite barrier of negligible thickness , and the corresponding potential reads


This model has the advantage to simplify both numerical and analytical approaches and it allows to study the effect of the barrier without introducing a specific shape for it.

The second potential (HGDW) is more realistic and it represents an isotropic harmonic potential (i.e., a spherical trap) plus a Gaussian barrier in one direction:


This is the structure of a realistic double-well potential in a cold-atom experiment.

Figure 1: Plot of the SDW (top) and HGDW (bottom) potentials. The figures represent , where is the potential energy. Particles are also pictorially shown, with their vertical coordinate representing their energy.

In Fig. 1 we represent schematically the two model potentials at , as a function of and . In this pictorial representation, the vertical coordinate of the particles corresponds to their energy.

In both cases, we denote as the barrier height: it is for the SDW, and , the difference of potential energy between the barrier top and the well minima in the case of the HGDW.

The quantity we use to follow the macroscopic dynamics of the double-well problem is the fraction of particles in the left well at time


The population imbalance is defined as


In the literature, as in the case of the superfluid dynamics in double-well potentials, is also often used as well the relative imbalance .

ii.4 Time scales

There are in general two fundamental time scales that rule the dynamics of a trapped interacting gas in a confining trap: one is related to the inter-particle interactions, the other to the external potential.

The average time between two consecutive collisions experienced by the same particle is the collisional time


On the other hand, the average time a particle needs to travel across the whole trap is


where the root-mean-square velocity at equilibrium is , and the system linear size. The reference length is chosen to be for the SDW model (and as well for the single square-well treated in Sec. III.1) and for the HGDW model (in the latter it is the distance between the two points of the harmonic trap having energy , neglecting the barrier).

Depending upon the frequency of collisions, the gas can be in different dynamical regimes, the two limiting cases being the collision-less regime (very rare collisions) and the hydrodynamical one (very frequent collisions): if the gas is collision-less, whereas it is hydrodynamical if .

For the double-well problem, we are interested in the time needed to smooth out the initial population imbalance, and our aim is to relate it to the basic properties of the system. This defines for double-well potentials a third time scale, related to the relaxation of imbalance, and that we will indicate with . Classically, the particles can cross the barrier only if their momentum perpendicular to it is high enough. Since the momentum is continuously redistributed between particles in collisional events, we expect that the relaxation time depends on . Also, an Arrhenius-type of behaviour is expected in the limit of large barriers ().

These three time scales (, , ) will of course depend on the potential considered. Throughout the next Sections, for different reasons, we deal with a square single-well, a square double-well, a harmonic oscillator and a harmonic-gaussian double-well. To keep as light as we can the notation, we use the following convention: within any subsection, unless otherwise stated and denoted (respectively with , , and ), we intend that the potential-dependent quantities are computed within the potential under consideration.

Iii A simple analytical model

In this Section we develop a simple analytical model to describe the effect of collisions on the relaxation dynamics.

This model takes inspiration from the tight-binding ansatz used for the Gross-Pitaevskii equation in double- or multi-well potentials Sme97 (); Jak97 (); Tro01 (). The Gross-Pitaevskii equation, describing the dynamics of a superfluid, is an equation for a complex wavefunction : one then introduces two degrees of freedom for (phase, , and number, ) and this leads for a double-well potential to the introduction of four degrees of freedom, two per well (say , and , in the left, , and right, , wells). These four degrees of freedom are not independent since the total number of particles is constant. Via these four non-independent degrees of freedom one then obtains a simplified, yet very good, description of the superfluid tunneling dynamics in the weakly coupled regime.

Of course, for a classical gas in a double-well potential, there is no tunneling and the distribution is not a complex number: anyway, one can still think to introduce suitable degrees of freedom per well and after coupling the two wells via the Boltzmann equation one obtains a description of the double-well dynamics. This leads to a simplified model, which allows for an approximate estimate of the relaxation time . It is the choice of the degrees of freedom in the separate wells (effectively one in the model below, or eventually more for a more accurate description) that characterizes the model. Such choice is suggested by the form of the potential and by the properties of the physical system at hand: as an example, in DeWeert88 () a simplified model was introduced to study the nonequilibrium distribution functions for electrons in the electrodes of a metal-insulator-metal junction.

In our case, we proceed to a rather drastic modelization of the Boltzmann equation dynamics taking into account that it is the energy barrier that suggests/sets an energy scale, dividing particles in two classes: the particles having energy larger than (which then can move from one well to the other) and the particles having energy smaller than (which do not). The mechanism for which particles can move from the latter class to the former is provided by interactions: two particles scatter below barrier and as consequence of the scattering one of the two, conserving energy, has an energy larger than . Finer details, such as higher-order scattering process, are neglected. As a result the model is not expected to give a quantitatively accurate prediction of the relaxation time. Nevertheless, since the estimate of the relaxation time is in general a not simple problem, and what in particular is difficult is the determination of prefactors, it provides simple analytical formulas which are in general qualitatively reasonable. In the regime of intermediate barriers () the agreement is found to be also quantitative.

We start with the single-well problem: this allows to classify the different types of collisions that will play a role in the double-well case and quantify their relative importance.

iii.1 Square single-well

Consider particles in a cubic box of volume having size . This is the same potential that will be considered in the next Sec. III.2, except the fact that there a filtering wall is added in .

In view of the double-well potential problem which we are interested in, we choose a reference momentum and a reference energy (remember that here this choice is arbitrary since there is no barrier yet). We define and [] the number of particles of the gas having () at time . Clearly, at any time, and, at equilibrium, collisions maintain them to their equilibrium values (given in Appendix B). However, what is their evolution in time if the system starts from a situation with ?

We make the following assumption for the distribution function:






is a constant, and the chemical potentials are constant (see Appendix B) and such that , so that , as it should. Notice that the choice Eq. (17) implies a distribution that is uniform in space, thermal for and and has a discontinuity in the momentum distribution at . Since the number of particles in the well is fixed there is only an independent parameter in Eq. (17) (given that the ’s are fixed by ).

Let us now consider, among all the possible collisions, those that will alter and . In a collision, each of the two incoming and outgoing particles can have momentum along above or below the reference values: 16 cases are therefore possible. Among these, there are 6 type of processes altering , namely:


For example, with the notation we mean that in the collision the two incoming particles have , whereas one of the outgoing ones (order does not matter) is above and one below reference. Let us indicate as , the rates of each kind of process: the evolution in time of satisfies


and , as it should. Each term is multiplied by an appropriate factor taking into account the changes in the process implies. By a direct computation of the rates it is found that


where . At equilibrium, and the rates of the processes just defined are equal two by two: , , and . The factors , are appropriate equilibrium phase space integrals detailed in Appendix C. They are proportional to the equilibrium collision rate: in fact, they are given by


where the reference energy and are adimensional functions (see Appendix C).

Inserting these results into Eq. (21) we get


or, equivalently, dividing by :


Therefore, according to this model, the relaxation dynamics of the single-well problem depends only upon the parameter ; all the others (, , , , and the volume ) enter only through the combination [see Appendix A, Eq. (48)] with the result of setting a time scale for the problem, i.e. a prefactor entering the time in which the particles move from above (below) to below (above) the reference energy.

The stationarity condition is


The fraction of particles that at equilibrium is above reference (have ) can be easily evaluated and is


In panel a) of Fig. 2

Figure 2: Rates at equilibrium. Panel a): equilibrium rates of the processes , , , , , and , defined in (20), and their sum vs . The rates are equal two by two and correspond to . The simulations (points) were performed for a square single-well potential, with parameters , , . Results with (full circles) and without (open circles) a barrier are shown. The theory prediction is shown for comparison (lines). Panel b): the adimensional functions , and (with linear scale in vertical axis) vs .

we show (lines) the various rates at equilibrium vs . The rates obtained in the simulation at equilibrium are also shown (points) as a check of the numerical algorithm. The empty circles denote the rates in a single-well, and the full circles denote the rates in a double-well. The presence of the filtering wall prevents collisions between particles on different sides, unless both have high enough energy; we also verified that and are separately very similar for the single square-well and the SDW since the effect of the tiny barrier on these bulk quantities is negligible.

From panel a) of Fig. 2 we also see that, out of all the collisions taking place, the processes affecting (full line) are, on the whole, always a percentage smaller than and, for have already decreased to . In our simulations, we will always consider . Concerning the relative importance of the different processes, and are the most probable ones for .

For completeness, we show in panel b) the adimensional functions , ; notice that these quantities, strictly related to those shown in panel a) [see Eq. (22) and Eq. (23)], are shown here using a linear scale in the vertical axis.

Finally, we would like to emphasize that, while of course the total collision rate depends on the trap shape (see Appendix A), the ratios do not (see Appendix C). In fact, the spatial dependence factors out and cancels between numerator and denominator, leaving only integrals over the momenta. Therefore, also the data obtained in a realistic trap would fall on top of the curves of Fig. 2a.

We now discuss the linearization of the dynamical equations Eq. (24) and the drawbacks of the analytical model.

iii.1.1 Linearization and relaxation time

The evolution in time of is ruled by the non-linear equation Eq. (24). If linearized around the equilibrium, it admits an exponentially decaying solution with the single-well relaxation time that (omitting details) reads


In the last equality we used the relations Eq. (15) and Eq. (27).

iii.1.2 Shortcomings of the analytical model

The simple assumption Eq. (17) is expressed in terms of just one dynamical degree of freedom, : this choice does not allow to respect energy conservation exactly. The energy of our gas is purely kinetic and can be computed at any time from the distribution function as


At equilibrium, it yields correctly . However, the assumption Eq. (17) implies the following energy variation




The function is positive for : it starts from and increases for larger values of . So, by construction, energy conservation is violated in the analytical model as soon as . We can expect the model to be better as long as remains small enough during the whole evolution. However, we explored different initial conditions having the same and we observed that even in cases in which become relatively large, the estimate of is good due to a compensation of effects. It also emerged that the estimate of appears to be better when is positive.

The limitations due to the non-conservation of the energy could be overcome by choosing a structure for with at least a further independent degree of freedom: for example, one could introduce a variable giving the amount of energy above reference energy .

iii.2 Square double-well

We now insert an energy barrier inside the box of Sec. III.1, obtaining the SDW potential. To be specific, assume the box is and the barrier located at . The barrier is perfectly transparent for particles with momentum and perfectly reflecting for particles with momentum . The energy associated to the barrier is .

We can now straightforwardly extend the model seen for the single-well to this case: we write


as in Eq. (17) with and denoting the well. Since the role played by is the same as in the single-well model, the ’s depends on in the same way and they do not depend on : . The numbers of particles in the two wells are and . Now we have four variables/components , , , . The equations of motion are


We have denoted with a the quantities referring to particles in a volume : of course the density is the same and therefore the collision rate per particle is the same too ().

We observe that the numbers , , , obey the condition ( is the total number of particles in the double-well) resulting in three independent parameters. However, during the first part of the dynamics, the particles above the barrier rapidly flow from one well to the other practically giving , and therefore the independent parameters in the subsequent dynamics are just two.

With respect to the single-well case Eq. (25), there is a qualitatively new term coupling the and sides of the barrier: it is a diffusion term giving the rate of particles passing from one side of the wall to the other


The coefficient (Arrhenius) is


A global factor in arises from the ratio between the area of the filtering wall and the volume of each partition: it would be replaced, in a more general geometry, by .

iii.2.1 Linearization and relaxation time

Linearizing Eqs. (33) around equilibrium, we find they admit the following solution


where the eigenvalues are found to be


We have defined the single-well rate , with the single-well relaxation time given by Eq. (28). We also denoted the eigenvalues so that is larger than . In the comparison with the numerical results (see next sections), since is associated to the short-time dynamics and to the long-time one, we plot as the analytical prediction for the relaxation time. In the limit of large cross sections, tends to the diffusive limit .

Since is small, one can get more insight approximating the eigenvalues to lowest order in :


Notice that is a sum of two terms, one depending on the collision (to which we may refer as ) and the other not (). The latter term is dominating for large , and for all it is .

We can also write as


where is the length scale associated to interactions.

Summarizing, in the cubic box with a filtering barrier (the DSW model) we get:


where the coefficients are


and is given in Eq. (27).

iii.3 Adapting the analytical model to the harmonic-gaussian double-well

As an approximation for the realistic double-well, we suppose we can still use Eq. (40), but with other coefficients:


with as given in Eq. (27).

In fact, the rates are the same in any potential , because they are global quantities in which spatial dependence cancels out (see Appendix Sec. C). What changes when passing from the toy to the realistic double-well are the rate of collisions, entering , and the characteristic frequency, entering . In the Eqs. (43) and (44) we have used the harmonic approximation of the realistic HGDW well close to its minima, where and is the rate for particles in a harmonic trap of frequencies (see Appendix A). As a consequence of the modification of , also the expression of the diffusive limit of is altered into .

Iv Numerical results for the square double-well

We show in this Section the numerical results of test-particle simulations for the SDW model potential. We consider a system with given , density and temperature (, , ) and vary the barrier height and the cross section . In all the simulations, the initial population is in the left well and in the right one. A comment is in order about units: in this Section we use units in which and .

As an example, we show in Fig. 3

Figure 3: Dynamics in the SDW potential. Time evolution of the well population for two barrier heights. Parameters: , , , . The barrier height is for the lower curve and for the upper one. The simulation curves are obtained averaging over runs with different microscopic initial conditions and smoothing over small time intervals. The fits are done with the function .

the evolution in time of population imbalance for two values of . The curves are obtained averaging over runs with different microscopic initial conditions and smoothing over small time intervals. These two procedures are necessary because the barrier crossing is a rare event and leads to large fluctuations in the well population. Notice that this is not a numerical artifact: also in an experimental realization it would be necessary to average over different realizations to be able to observe the population evolution in time clearly for large barriers. As a consequence of the smoothing, the curves in Fig. 3 start from values different from : for example, in the case of a very low barrier, , the adjustment from to was very fast.

The obtained behaviour can be nicely fitted with a single exponential decay


that allows to extract the relaxation time . For other choices of the cross section or other barrier heights, the evolution in time of is sometimes more complex, showing an oscillatory behavior at the early times or an exponential decay with two clearly separated time scales. Therefore, in the following, to extract the relaxation time from , we will also consider the fitting functions and .

In Fig. 4

Figure 4: SDW potential relaxation times. The inverse relaxation time , in units of , as a function of for different values of the cross section: with . The numerical results are denoted by points. The system is a box of size , with particles, at and different barriers . Each point is obtained from an average over runs ( for the largest value of the cross section ) and smoothing over small time intervals. The fits are done in the whole time interval, using one of the three functions (full circles), (stars) and (full squares) defined in the text. Different colors indicate different values of the cross section. For comparison, we show (solid lines) the results obtained in an approximated analytical solution presented in Sec. III, namely obtained using Eqs. (40), (27), (41) and (42). The diffusive limit is also shown (dashed line). Inset: rescaled in units of collisional time . Colors correspond to the same cross section as the main figure.

we show the results for the relaxation time for an initial imbalance of [], and different values of the cross section as a function of . The relaxation time is shown in units of , that for the toy model it . Any point is the result of a fit done on a curve obtained averaging ( in the case of the largest cross section) runs with different microscopic initializations. The fits are done on the whole available time interval, and the chosen fitting function is the one leading to the smaller . The point colour (on-line) indicates the value of the cross section, the point shape the fitting function: (full circles), (stars), (full squares). In the inset we show again the inverse relaxation time, but in units of the collisional time: in this scale all the points fall in a narrow region (notice that the range of values in the vertical axis of the inset is much smaller than in the main plot), showing that it is the collisional time that gives the major contribution to .

In Fig. 5

Figure 5: SDW model potential. Inverse relaxation time , in units of , as a function of for a given barrier height, . Parameters and the meaning of symbols and lines are the same of Fig. 4. The numerical results are denoted by points and the model prediction by lines.

we focus on a single barrier height and show the dependence of (again in units of ) upon the interaction strength.

iv.0.1 Comparison with the analytical model predictions

In Fig. 4 and Fig. 5 we compare the numerical results with the predictions of the analytical model for the SDW (lines). The model in general predicts a faster relaxation (smaller ) than what found in the simulation, and as the barrier height increases the agreement is deteriorated. In Fig. 4 one sees, however, that for low and intermediate barriers the agreement is rather satisfactory.

A feature that the model captures nicely is the dependence of the relaxation time upon the interaction strength: in Fig. 4 one can see that, as the cross section increases, both types of curves accumulate towards the diffusive limit. This finding is better seen in Fig. 5, where we fix the relative barrier height at an intermediate value () and study the evolution of with the interaction strength.

V Numerical results for the harmonic-gaussian double-well

In this Section we present our results for the realistic HGDW, Eq. (12), that is a combination of a spherical harmonic trap and a gaussian barrier along one direction. We start with an initial configuration with of the particles in the left well and in the right one; their initial momentum distribution is the thermal equilibrium one. To prepare this configuration numerically, we start with the balanced population and then move of the particles from the right to the left well.

We then let the system evolve: results from a combination of a center of mass oscillation due to the harmonic trapping, and the relaxation of population imbalance. We choose to fit such behaviour with the fit function , where the third term represents the damped center of mass oscillation. We have checked this interpretation exciting explicitly the c.o.m. motion, i.e., considering a balanced cloud and displacing it on the whole by a certain amount. The frequency and damping of this c.o.m. mode are in reasonable agreement with those extracted from the oscillation that arises when the system is prepared with an initial population imbalance and not displaced.

We fix the double-well shape (i.e., , , ) and study the relaxation dynamics for different temperatures and interactions. In all the cases we have particles and the initial imbalance is . The curves are obtained averaging over different microscopic realizations and smoothing over small time intervals. Since we are in presence of a harmonic trap, in this Section we use the harmonic oscillator units (or trap units), in which all the dimensional quantities are built from as usual. For example, , and so on. To pass to physical units, values have to be chosen for the mass and the trap frequency .

In Fig. 6

Figure 6: Dynamics in the HGDW. Evolution of with time, for three values of the temperature. The gas contains particles, with an initial imbalance of . The well parameters, in trap units, are , , (therefore ), and the temperatures (panel a), (panel b) and (panel c). To express them in dimensional units, reported in the top -axis, we need to specify the trap frequency and the particle mass: for Li and Hz they correspond to , , , . The values of and fitted from the data for figures (a,b,c) in trap units are respectively and .

we show three typical behaviours of : they correspond to a given cross section and different temperatures, increasing from left to right. The frequency of the c.o.m. oscillation is very close to at high temperatures, where the presence of the barrier does almost not affect the cloud oscillation; at low temperatures, instead, it is reduced to smaller values.

Repeating similar calculations for a set of temperatures and interaction strengths, we obtain the results shown in Fig. 7. In some cases, in which the c.o.m. oscillation is not visible anymore, we use the fitting function to extract the relaxation time. As before, the symbol shape indicates the fitting function used, with the same notation of Fig. 4.

Figure 7: HGDW relaxation time. The inverse relaxation time as a function of , obtained varying the temperature (the trap is fixed as in Fig. 6 with a barrier height ). The numerical results are denoted by points. Different colors correspond to different values of the cross section: . Each point is obtained from an average over runs. The fits are done in the whole time interval, using either (full circles), or (stars). The dimensional quantities reported in the figure are obtained as in Fig. 6. For comparison, we show as solid lines the results obtained in an approximated analytical solution presented in Sec. III, namely obtained using Eqs. (40), (27), (43) and (44). Inset: rescaled in units of collisional time . Colors correspond to the same cross section as the main figure.

Comparing with the analogous plot for the SDW, Fig. 4, we see an analogous accumulation of curves as the cross section increases. The dependence upon seems qualitatively different: the point is that here, along a curve for a given cross section, the global collision rate changes, whereas this was not the case in Fig. 4. If we rescale all the curves by the corresponding equilibrium collision rate (inset), we see that they all fall in the same region and the trend so obtained is therefore similar in SDW and the HGDW.

In panel a) of Fig. 7, together with the numerical results, we show also (full lines) the analytical predictions obtained extending our model to the realistic well case (see Sec. III.3).

Next, in Fig. 8,

Figure 8: HGDW potential. Inverse relaxation time as a function of for the same trap of Fig. 6 and a given temperature (, then ). The fitting function and its parameters are the same as in Fig. 6.

we show the analogous of Fig. 5 for the realistic HGDW well, for : the qualitative behaviour is the same and it is nicely reproduced by the analytical model.

In the figures of this Section (Fig. 6, Fig. 7, Fig. 8), together with the harmonic oscillator units, we show also axes with physical units: they are obtained assuming a reference mass (that of ) and a trap frequency Hz.

v.0.1 Estimates for a two-component Fermi gas

Finally, we can use our results to make an approximate prediction for a balanced two-component mixture of . In fact, in this case we would have two species, equally populated () with only inter-species interactions. In the Boltzmann framework, we would have two distribution functions, but they coincide if the mode and the potential do not depend upon the “spin”: just one distribution normalized to is needed. Our calculations for classical particles correspond therefore to a balanced mixture of fermions. Anyway, notice that we are not including Pauli blocking in the collision term and we are approximating the cross section with a constant. With this in mind, in Fig. 9

Figure 9: HGDW potential. Inverse relaxation time as a function of , for a fixed barrier . The results for are the same of Fig. 7, here shown in a different representation. The Fermi temperature has been computed as : that is, the Fermi temperature for a balanced mixture having in a harmonic well. For details on the parameters and the dimensional units, see the caption of Fig. 6.

we plot the same results of Fig. 7, as a function of , with the Fermi temperature in a harmonic trap, and . Notice that even we are not considering the effect of the quantum statistics on the collision term, both numerical computations and the analytical model can be straightforwardly extended to include it: we expect that, at temperatures (and in any case larger than ), the relaxation time has a similar qualitative dependence on barrier heights and interaction strengths.

Vi Conclusions

In this work we have presented a study of the double-well dynamics of a classical gas that obeys the Boltzmann equation, with the purpose of assessing the role played by collisions in the relaxation of population imbalance. We think that a detailed study of the Boltzmann equation in a variety of double-well potentials is a rewarding subject of interest, not only for its paradigmatic and pedagogical importance, but also to concretely model currently ongoing experiments in a range of temperature and to set a basis for further quantitative theoretical studies of tunneling of ultracold atoms at finite temperature below , in particular near .

Two model potentials have been considered: a toy square double-well (SDW) potential with a filtering wall and a more realistic double-well obtained by combining a gaussian to a harmonic potential (HGDW). In both cases, we have performed numerical simulations (test-particle method) of the relaxation dynamics from an initial imbalanced population of the symmetric wells to the balanced equilibrium, in a range of cross sections values and barrier heights . For convenience, in the toy SDW potential we have fixed the temperature and varied the barrier height, whereas in the realistic HGDW one we have fixed the trap shape, therefore the barrier height, and varied the temperature.

Beside the numerical results, we have also presented a simple analytical model for the dynamics, that allows to compute the relaxation time analytically from the system parameters. The agreement of the model predictions with the results of the simulations is qualitative, and it is quantitatively better for the realistic HGDW potential. In general the analytical findings are in reasonable agreement with the numerical results for low up to intermediate barriers ().

Finally, we have used our results to estimate the relaxation times for realistic values of on-going experiments with a mixture of fermionic cold atoms ().

As a possible continuation of this study, it would be interesting to compare our results with those obtained in the framework of Klein-Kramers equation in presence of a barrier Hanggi90 (): to this end, the friction parameter appearing in the Klein-Kramers equation should be appropriately (we mean, quantitatively) computed.

Another extension of this work, relevant for ultracold atom experiments at low temperature, would be to include, beyond the classical crossing mechanism studied here, the hopping via quantum tunneling: a possible way could be to include an effective ad hoc term in the Boltzmann equation or to couple the Boltzmann equation to the equations for the superfluid dynamics.


S.C. is supported by the “Fundação para a Ciência e a Tecnologia” (FCT, Portugal) and the “European Social Fund” (ESF) via the post-doctoral grant SFRH/BPD/64405/2009. A.T. acknowledges support from the Italian PRIN “Fenomeni quantistici collettivi: dai sistemi fortemente correlati ai simulatori quantistici” (PRIN 2010_2010LLKJBX).

The authors are grateful to D. Davesne, A. Gambassi, A. Laio and I. Vidaña for valuable discussions. We also gratefully thank for many useful discussions people of LENS and QSTAR groups in Florence, in particular A. Burchianti, K. Xhani, G. Roati, A. Smerzi and M. Zaccanti. A.T. acknowledges the University of Coimbra for kind hospitality and the Isaac Newton Institute for Mathematical Sciences, Programme “Mathematical Aspects of Quantum Integrable Models in and out of Equilibrium”, where the final part of this work was completed. S.C. acknowledges CNR-IOM for kind hospitality in Trieste. The authors also thank the Laboratory for Advanced Computing at the University of Coimbra for providing CPU time in the Milipeia and Navigator clusters.

Appendix A Collision rates at equilibrium

In a generic potential and for constant cross-section, the collision rate at equilibrium reads