text:
# Vertically shaken column of spheres. Onset of fluidization

## 1 Introduction

Granular material confined in a vertically vibrating container reveals complex effects such as surface structure formation (e.g. Coulomb:1773 ; MeloUmbanhowarSwinney:1994 ), spontaneous heap formation (e.g. DinkelackerHueblerLuescher:1987 ), convection (e.g. EhrichsJaegerKarczmarKnightKupermanNagel:1994 ) and others. A precondition common to all these effects is that (at least) the particles at the free surface of the granular material lose contact to their neighbours for at least a small part of the oscillation period. If a granular material is agitated in a way so that the particles at the surface separate from their neighbours we will call the material fluidized.

The conditions under which surface fluidization occurs are not clear yet. For sinusoidal vertical excitation of the container most of the literature reports that surface fluidization starts as soon as the acceleration amplitude $`A\omega ^2`$ of the oscillation $`z=A\mathrm{cos}\omega t`$ exceeds gravity $`g`$ or another critical constant. The Froude number is defined as $`\mathrm{\Gamma }=A\omega ^2/g`$ and, hence, it has been reported in numerous publications that surface fluidization or, respectively, effects which require fluidization, occur for $`\mathrm{\Gamma }>1`$.

Few references BarkerMehta:1993a  report, however, that in numerical simulations surface fluidization was observed for $`\mathrm{\Gamma }1`$. So far these results have not been confirmed experimentally.

Whereas a rigid solid body, e.g. a single sphere on a oscillating surface, would certainly start to jump if the acceleration amplitude $`A\omega ^2`$ exceeds gravity $`g`$, there are arguments which may lead to a different conclusion for a shaken amount of granular material, i.e., a many body system:

* The deformation-force-law of contacting bodies may be nonlinear due to geometrical effects, even if the material deformation is small enough to assume linear material properties, i.e., Hooke’s law. For contacting ideally elastic spheres it has been derived by Hertz Hertz:1882  that the interaction force $`F`$ scales with the deformation $`\xi `$ as $`F\xi ^{3/2}`$. This law applies for all particle contacts (under mild conditions) provided the deformation $`\xi `$ is small enough Brilliantov . As a result, nonlinear phenomena like the propagation of solitary waves has been shown to appear in bead chains with Hertz contact law CosteFalconFauve:1994PRE .
* Under the influence of gravity the particles are deformed differently even at rest, depending on their vertical position in the material. Therefore, the speed of sound in the material is not uniform but a function of the vertical coordinate. This property may lead to complicated pulse motion through the material and has been reported also for chains of contacting spheres SinkovitsSen:1995 .
* In a polydisperse granular material the geometrical properties, i.e., the contact network, varies with the applied force. If the material is loaded the particles are deformed and more and more contacts emerge. This effect leads to a more complicated deformation – force law, $`F\xi ^\gamma `$ with $`1.5\gamma 4`$ HerrmannStaufferRoux:1987a .

In a recent theoretical paper BelowG  the granular material was modeled as a vertically shaken column of viscoelastic spheres. The main result of BelowG  is that the column may be fluidized even if the condition $`\mathrm{\Gamma }>1`$ does not hold. Instead a different condition for fluidization was derived: Assuming viscoelastic material properties the force between two contacting spheres of radius $`R`$ at vertical positions $`z_k`$ and $`z_{k+1}`$ reads BSHP

$$F_{k,k+1}=\sqrt{R}\left(\mu \xi _{k,k+1}^{3/2}+\alpha \dot{\xi }_{k,k+1}\sqrt{\xi _{k,k+1}}\right),$$
(1)
where $`\xi _{k,k+1}2R\left|z_kz_{k+1}\right|`$ is the compression and $`\mu `$ and $`\alpha `$ are elastic and dissipative material constants (details see BSHP ). The height of the column of $`N`$ spheres is $`L=2NR`$. It was found that the top sphere separates from its neighbour if

$$\frac{A\omega ^2}{g}>\mathrm{\Gamma }\left(\frac{3}{5}\right)\left|M^{2/5}J_{2/5}\left(2M\right)\right|,$$
(2)
where $`J_{2/5}`$ is the (complex) Bessel function and $`M`$ abbreviates

$$M\frac{(18\pi )^{1/3}}{5}\left(\frac{\mu L^5\rho ^2}{g}\right)^{1/6}\frac{\omega }{\sqrt{3\mu 2i\omega \alpha }},$$
(3)
with $`i=\sqrt{1}`$ and $`g`$ and $`\rho `$ being gravity and material density.

For small frequency $`\omega `$ the Taylor expansion of the rhs. of the condition (2) yields

$$\frac{A\omega ^2}{g}>1B_2\omega ^2+B_4\omega ^4$$
(4)
with

$`B_2`$ $``$ $`{\displaystyle \frac{18^{\frac{2}{3}}}{45}}\left({\displaystyle \frac{\pi ^2\rho ^2L^5}{g\mu ^2}}\right)^{\frac{1}{3}}`$
$`B_4`$ $``$ $`{\displaystyle \frac{1}{18^{\frac{2}{3}}}}\left({\displaystyle \frac{\pi ^2\rho ^2L^5}{g\mu ^2}}\right)^{\frac{2}{3}}({\displaystyle \frac{3}{100}}+{\displaystyle \frac{4324^{\frac{2}{3}}}{405}}{\displaystyle \frac{\alpha ^2}{\pi ^2}}\left({\displaystyle \frac{g\pi ^4}{\mu ^4\rho ^2L^5}}\right)^{\frac{1}{3}})`$
Both coefficients $`B_2`$ and $`B_4`$ are positive definite values, i.e., for all materials there is a range of frequency where the rhs. of the inequality (2) is smaller than one. This means that there exists always a frequency interval where the column fluidizes although the acceleration amplitude of shaking is smaller than gravity, $`A\omega ^2/g<1`$. The full derivation of the sketched theory and the discussion of the frequency range where the condition (2) holds can be found in Ref. BelowG . For our experimental system consisting of steel spheres, of density $`\rho `$=7700 Kg/m3, Young modulus $`Y`$=210 GPa, Poisson ratio $`\nu =0.29`$, in a column of $`L=`$0.6 m, the Taylor expansion Eq. (4) holds for frequencies $`\omega `$ much smaller than 1000 s-1, which covers entirely the experimentally accessible range of frequencies.

In the present paper we want to report experimental results on the onset of fluidization of a vertically shaken column of spheres, i.e., on the same system which was studied theoretically in Ref. BelowG . The dynamical behaviour of a vibrated one dimensional system of spheres in the well fluidized regime was studied e.g. in BernuDelyonMazighi:1994 ; x1 . This regime is explicitely not considered here. In x  a column of particles which interact as linear springs with linear damping, $`\ddot{\xi }+2\mu \dot{\xi }+\mathrm{\Omega }^2\xi =0`$, was studied with Molecular Dynamics. For their sets of parameters the authors found that fluidization occurs above $`\mathrm{\Gamma }1.1`$, however, the authors claim that the onset of fluidization depends on the choice of the parameters $`\mathrm{\Omega }`$ and $`\mu `$. Within their model it might be possible to get also fluidization for $`\mathrm{\Gamma }<1`$, depending on the material parameters. If one considers spheres, however, in the limit of vanishing damping the Hertz law implies that the duration of contact behaves as $`t_cg^{1/5}`$, whereas for the model considered in x  one finds a constant $`t_c=\text{const.}=\pi /\mathrm{\Omega }`$. Hence, the linear spring model certainly does not describe the contact of spheres, not even in the case of pure elastic contact. For the same linear spring system in xx  the detachment effect was reported, the same effect for viscoelastic spheres was discussed in xxx .

## 2 Experimental setup

We expect that the effect to be measured is quite small, i.e., for the critical parameters of shaking when the material starts to fluidize we expect $`1A\omega ^2/g`$ to be a small value. Therefore, we have to adjust the amplitude of shaking $`A`$ and in particular the frequency $`\omega `$ with good accuracy. Moreover, we have to assure that no energy originating from different frequencies than the driving frequency enters the column of spheres, e.g., due to high frequency noise.

The column of spheres was confined in a framework of three aluminum bars, which assures pure one dimensional motion of the spheres at low friction, see Fig. 1.

The framework with the spheres was mounted on a precisely vertically balanced linear bearing which was driven by a stepping motor via a crankshaft with adjustable eccentricity, Fig. 2.

The motor was computer controlled, i.e., in precisely 20,000 impulses the motor axis revolves once. This high angular resolution provides quasi-steady motion. The motion was smoothened additionally using a flywheel (with diameter 160 mm, thickness 32 mm, mass 4.5 kg and moment of inertia $`J=0.0144\text{kg\hspace{0.17em}m}^2`$) which was fixed on the axis and a hard rubber clutch connecting the motor axis and the axis of the shaking device (see Fig. 3). Using an acceleration sensor we checked that the finite step size per computer signal does not influence the sinusoidal motion of the column of spheres. The entire mechanical device was mounted on a massive oscillation damping table. Figure 3 sketches the device and Fig. 4 shows a photograph.

In the experiments the column consists of up to 20 ball bearing balls of radius $`R=12.5`$ mm. There are several possibilities to determine the onset of fluidization, i.e., for fixed amplitude $`A`$ the critical frequency $`\omega `$ when the top sphere separates from its neighbour:

* Direct observation: a comoving camera was mounted to the probe carrier and a LASER pointer was fixed at the opposite side of the column. As soon as fluidization starts at the critical frequency the emerging submillimeter gap between the top particle and its neighbour can be recognized as a flash on the screen.
* Acoustic method: When fluidization occurs there emerges a periodic acoustic signal originating from the collision of the top particles. This signal can be detected either using a microphone or the acceleration sensor which was mounted on the top of the device, see Fig. 3.
* Electric circuit: We attached a very thin wire (diameter $`<0.1`$ mm) to the top particle and its neighbour and observed the contact of the spheres using an oscilloscope.

We checked all three methods and found that the results are very similar. We have chosen the electric method since it turned out to be the most simple and most reliable one. The direct observation requires complicated mechanical tuning and for the acoustic method one has to filter the periodic signal from other sources of sound as, e.g., the sound of the motor.

## 3 Raw Data

For fixed amplitude $`A`$ we want to determine the according smallest frequency $`\omega `$ for which the top sphere separates from its neighbour in each period. Figure 6 shows a sketch of the oscilloscope signal for subcritical (top), critical (middle) and over-critical frequency, all for the same amplitude of shaking. For subcritical frequency the top particle follows the motion of the vibrating table resulting in a sinusoidal signal of the acceleration sensor and a horizontal line for the conductivity (limited by a resistor). At over-critical frequency the conductivity drops to zero for a certain time interval in each period. Close to the critical frequency the electrical contact between the spheres is not in each period interrupted, instead we get a noisy signal.

To collect data points we tuned the adjustable crankshaft to a desired amplitude $`A`$. Then the motor was accelerated up to a rotation speed which is slightly less than the speed where we expect fluidization (initial ramp). At this value all spheres are permanently in contact (top in Fig. 6). From this point on we increased the motor speed in very small steps until we observe separation of the top sphere in each period (bottom in Fig. 6). This procedure was repeated 10 times for all amplitudes independently for 2 and 20 spheres. The critical frequency for each amplitude was determined as the average over the 10 independent measurements. Figure 7 shows the raw data points with error bars.

The somewhat unusual representation of the data as acceleration vs. amplitude is due to the fact that our experimental input is the amplitude, the output, as stated above, is the according critical frequency and hence the acceleration.

In Figure 8 the critical acceleration of the container ($`A\omega ^2`$) vs. frequency $`\omega `$ for 2 and 20 beads is shown. One can see that even for two beads the critical acceleration, which is supposed to be a constant varies significantly with the frequency, i.e., it decreases with increasing frequency. The curve for 20 beads shows a similar behaviour, however, the critical acceleration for 20 beads is smaller than the critical acceleration for 2 beads. This is a first indication to an amplification effect similar as theoretically predicted. The experimental device reveals its own resonances, namely 7.6 Hz for two beads and 11.8 Hz for twenty beads. These resonances are caused by the mechanical properties of our experimental setup but not by the properties of the columns of beads themselves. Data points which belong to these frequencies have been eliminated for further analysis (see Fig. 7).

The data points are affected by large statistical errors, i.e., the data points scatter around a smooth curve which one expects to find if only systematic effects were relevant. Therefore, we have to apply a sophisticated data analysis to extract the amplification effect in order to compare it with the theoretical prediction. This procedure is described in the following sections.

## 4 Discussion of the Measuring Procedure

The measurement is affected by different systematic and random errors:

* Due to the generation of the shaking motion via a connecting rod of finite length ($`l=18`$ cm) there is a systematic deviation of the motion from the ideal sinusoidal shaking. Instead of a sinusoidal oscillation the generated oscillation is

$$z=A\mathrm{cos}\varphi +\sqrt{l^2A^2\mathrm{sin}^2\varphi }l.$$
(5)
For small ratios $`A/l`$ we can write

$$zA\mathrm{cos}\varphi \frac{A^2}{2L}\mathrm{sin}^2\varphi $$
(6)
This leads to a constant shift in the extremal acceleration of the container

$$\ddot{z}_{Extr}\pm A\omega ^2\frac{A^2}{L}\omega ^2=A\omega ^2\left(\pm 1\frac{A}{L}\right)$$
(7)
We see that the absolute value of the acceleration at the upper turning point is too large by $`A/l`$ which can be as much as $`3\%`$.
* Uncertainty in determining the driving frequency and amplitude are very small as described above.
* Due to the noise caused by the motor and the friction of the linear bearings the systematic (ideally sinusoidal) vertical motion of the column is superposed by a spectrum of high frequencies. These waves are source of an undesired extra energy feed which is not controllable and causes errors. Measurements of this noise using the acceleration sensor have shown that it does not depend much on the height of the column, it is determined mainly by the motor and the bearings. This property suggests to measure our effect not directly but to compare the results for a column of 20 steel spheres with the results for a “column” of two spheres (see below).

## 5 Data Analysis

To eliminate systematic errors we compare the result for a column of 20 spheres with the equivalent measurement for a “column” of only two spheres. The column of two spheres serves as a reference. Actually, the second bead is fixed to the bottom and oscillates rigidly while the top bead is free. Effectively there is only one free moving sphere, for one single sphere taking off from a flat wall could not serve as a reference for the case when the top bead of a column of $`N`$ spheres separates from the rest. This approach is based on the assumption that the systematic errors originating from the motor noise and from the bearings affect both systems almost in the same way. With the above mentioned method (see section 2) we determined the critical frequency of driving at which the topmost bead starts to jump for a given shaking amplitude. From this we calculated the critical acceleration. Due to the influence of systematic errors also the curve $`a_{\mathrm{crit}}(\omega )`$ for two beads, which is supposed to be constant, shows a significant dependence on the driving frequency, resp. amplitude. As already mentioned before, in Fig. 8 one can see that the critical acceleration decreases with increasing frequency. The curve for 20 spheres shows a similar behaviour but the critical amplitude is always smaller than the critical amplitude for 2 spheres at the same driving frequency. To identify the amplification effect, i.e., to eliminate the error of the apparatus we divide the critical acceleration for 20 spheres by the critical amplitude for 2 spheres at the same frequency which yields the absolute value for the critical Froude number for 20 spheres. This method is applicable if all significant systematic errors affect the measurements by a factor which is independent on the number of beads. This precondition holds certainly for the influence of the limited length of the connecting rod and holds to a reasonable degree for the influence of the higher harmonics in the driving motion due to the discrete steps in the motor motion. Since we believe these errors to be the most significant ones our method should give reliable results.

Following this idea we fit both sets of data, for 2 and 20 spheres, to the function

$$f(\omega )=\left(a_2\omega ^2+a_1\omega +a_0\right)^1$$
(8)
and take the ratio of these curves. There are other possibilities to construct fit-functions which reflect the desired properties of monotonous decay and significant curvature, it turns out that the final result does not depend sensitively on the choice of this function.

The ratio curve, which represents the predicted amplification effect, can be described in good approximation as a decreasing parabola with maximum at $`\omega =0`$. Its curvature is small, in agreement with the expansion given by Eq. (4). Dissipative effects, represented by the material constant $`\alpha `$, are not observable at these frequencies, as the restitution is high for all velocities (0.90 or more): we estimated the second term of $`B_4`$ being of the order of 10-6 or less. This curve, multiplied by the reference acceleration 9.81 m/s2, is shown in Fig. 10. The ratio between the fitting curves in Fig. 9 is not exactly 1 at $`\omega =0`$, but by about $`3\%`$ smaller. This discrepancy is due to the above mentioned systematic errors, namely due to the higher harmonics in the motion of the stepping motor (see Fig. 11) which can have different magnitude at different load. Therefore the division of the curve for 20 beads by the reference curve alone can not completely remove this systematic error. Other systematic errors such as deviations from an ideal sinusoidal motion due to the finite length of the connecting rod have been removed by the data analysis procedure.

The influence of the higher harmonics can be observed in Fig. 11 showing the output of the acceleration sensor on top of the column. The sinusoidal motion of the container with driving frequency $`10`$ Hz is superimposed by an oscillation of frequency of about $`150`$ Hz. When calculating the maximum acceleration of this motion one gets a value which is up to 20% higher than $`A\omega ^2`$. Thus, by taking $`A\omega ^2`$ as the maximum acceleration of the container we underestimate the actual acceleration, which, since the magnitude of the high frequency oscillation varies with different load, contributes to the difference between the maximum of the ratio curve in Fig. 10 and the theoretical value of 9.81 m/s2.

## 6 Discussion

We investigated experimentally the onset of fluidization of a vertically vibrated column of spheres. It was found that one observes fluidization even if the amplitude of acceleration of the vibrating motion is smaller than the acceleration due to gravity, $`A\omega ^2/g<1`$. This result is in agreement with the theoretical prediction BelowG . The quantitative theoretical result for the critical parameters for the onset of fluidization given by Eqs. (2) and (4) does not contradict the experiment but could also not be conclusively confirmed by the experiment due to insufficiencies of the experimental setup. The most important shortcoming is the existence of high frequency vibrations mainly due to the motor noise and the bearings. While qualitatively agreement is shown, due to these limitations the experimental measurements are not completely conclusive to verify the theoretical results: a quantitative check of the theoretical predictions requires a more sophisticated experimental setup, directed to reduce noise and enhance the effect. The noise can be reduced by improving mechanical isolation: enclosing the column of spheres into a massive block or wall, using very rigid bearings and separating the motor at a good distance from the column. Higher columns or more balls will help to reduce the noise, but one has to work harder on a correct vertical alignment; for this reason the beads should not be very small. The use of a softer material alone would enhance the effect, but then one would like still good restitution, smooth surface and good conductivity (for the electrical method of determination), so the choice of material is not straightforward. Also higher frequency will enhance fluidization at lower Froude numbers, but then controlling amplitude (and noise) will become much more critical.

###### Acknowledgements.

We thank Hans Scholz for discussion and Christine Rosinska for technical aid.


format: text

text:
# Kondo effect out of equilibrium in a mesoscopic device

## Abstract

We study the non-equilibrium regime of the Kondo effect in a quantum dot laterally coupled to a narrow wire. We observe a split Kondo resonance when a finite bias voltage is imposed across the wire. The splitting is attributed to the creation of a double-step Fermi distribution function in the wire. Kondo correlations are strongly suppressed when the voltage across the wire exceeds the Kondo temperature. A perpendicular magnetic field enables us to selectively control the coupling between the dot and the two Fermi seas in the wire. Already at fields of order 0.1 T only the Kondo resonance associated with the strongly coupled reservoir survives.

Mesoscopic devices form a powerful tool for the study of fundamental many-body phenomena. A striking example is the Kondo effect in a quantum-dot device, which consists of a small electronic island connected by tunable tunnel barriers to extended leads . Predicted in 1988 , the Kondo effect in a quantum dot was observed for the first time ten years later , leading to intense research . Quantum dots offer new control to study various aspects of the Kondo effect. One option is the possibility to explore non-equilibrium regimes, e.g. by applying a finite voltage between the leads connected to the dot.

At equilibrium, Kondo correlations give rise to a sharp resonance in the density of states (DOS). The resonance is aligned with the Fermi energy of the leads. For a finite bias voltage the resonance is predicted  to split as shown qualitatively in Fig. 1a. The resulting two peaks are aligned with the Fermi energies of the leads, and become progressively smaller as the bias voltage is increased. This suppression is due to decoherence introduced by the inelastic scattering of electrons from the high- to the low-energy lead .

Measuring differential conductance, $`dI/dV`$, as a function of bias voltage, $`V`$, results in a nonlinear characteristic with a single peak at zero bias . This zero-bias anomaly, observed in a variety of experiments, is a characteristic signature of the Kondo effect. However, a $`dI/dV`$-vs-$`V`$ measurement does not provide enough information to extract the $`V`$-dependent DOS. In particular, whether the Kondo resonance is indeed split and how the resulting double-peak DOS evolves with bias is an issue that has not yet been experimentally investigated. Whether Kondo correlations survive at voltages larger than the Kondo temperature at equilibrium, is a question that has recently raised considerable debate .

Some authors  have argued the possibility of a two-channel Kondo effect that would develop at finite bias. A direct measurement of the DOS would contribute important information to understand Kondo physics out of equilibrium. Recent theoretical papers  have proposed to use three-terminal quantum dots in which one of the leads serves as a weakly coupled electrode to probe the local DOS. Here we address this issue from a different angle.
We couple a quantum dot to the middle of a quasi-ballistic quantum wire that is brought out of equilibrium by an applied bias voltage. The finite voltage creates inside the wire two ensembles of electrons (i.e. right and left movers) whose quasi-electrochemical potentials, $`\mu _{S1}`$ and $`\mu _{S2}`$, are defined by the electron reservoirs connected to the wire. Accordingly, the electron distribution function in the wire develops two steps located at $`\mu _{S1}`$ and $`\mu _{S2}`$  (see Fig. 1b, left side). Our main goal is to show that the Kondo resonance splits for $`\mu _{S1}\mu _{S2}`$ (see Fig 1b, right side).

The sample designed for this experiment was fabricated on a GaAs/AlGaAs heterostructure containing a two-dimensional electron gas (2DEG) 90 nm below the surface. (The experiment was reproduced on a sample with similar geometry, fabricated on a different heterostructure.) As shown in Fig. 1c, the device has a small quantum dot (indicated by a dotted circle), formed inside a 500-nm-wide region defined by shallow etching. Confinement in the longitudinal direction is accomplished by applying negative gate voltages, $`V_{gs}`$ and $`V_{gd}`$. These gate voltages define the tunnel coupling between the quantum dot and its 2DEG leads. A third gate voltage, $`V_g`$, is used to control the electrostatic potential on the dot. On the top side, the dot is connected to the middle of a short quantum wire formed in the 2DEG by applying a voltage $`V_{ch}=0.4`$ V to a large gate coming from above. The wire then has a resistance, $`R_{ch}`$, of about 1.3 k$`\mathrm{\Omega }`$, corresponding to $``$10 spin-degenerate propagating modes. There are in total three regions of 2DEG connected by ohmic contacts to the external circuitry. The bottom 2DEG (drain) is connected to ground via a current-voltage converter that measures the current through the quantum dot. The other two regions, source 1 and source 2, are connected to dc voltages $``$$`V`$ and $``$($`V+\mathrm{\Delta }V`$).

An additional ac voltage $`\stackrel{~}{v}`$ with a rms amplitude of 2.5 $`\mu `$V is applied to both source leads allowing for lock-in measurements. The resistances of the leads, $`R_1`$ and $`R_2`$, are in the 0.5$``$0.8 k$`\mathrm{\Omega }`$ range, including the contribution from the respective 2DEGs. We call $`\mu _{S1}`$ and $`\mu _{S2}`$ the Fermi energies of the source 2DEGs, and $`\mu _D`$ the Fermi energy of the drain 2DEG.
Measurements were done in a dilution refrigerator at 15 mK. The quantum dot parameters ($`V_{gs}`$, $`V_{gd}`$, and $`V_g`$) were tuned to achieve a robust Kondo regime. In the inset to Fig. 2a we show the linear conductance, $`G`$, of the dot as a function of $`V_g`$ ($`V`$ and $`\mathrm{\Delta }V`$ are set to zero). The Kondo regime (see arrow) takes place between the two Coulomb blockade peaks. Here the valley conductance is enhanced. By sweeping $`V`$ for $`\mathrm{\Delta }V=0`$ we measure $`dI/dV`$ as a function of bias voltage $`V`$ across the dot. As in several earlier experiments, we find a narrow peak at zero bias (Fig. 2a) which reflects the existence of a Kondo resonance at the Fermi energy. The peak has a full width at half maximum, $`w_025`$ $`\mu `$V, yielding a Kondo temperature $`T_Kew_0/k_B=0.3`$ K.
In Figs. 2b, 2d and 2f we show on grey scale $`dI/dV`$ vs ($`V`$,$`V_g`$). Each plot is obtained from many $`dI/dV`$ vs $`V`$ traces taken for a set of closely spaced $`V_g`$ values ranging between -65 and -20 mV. Figure 2b, taken at $`\mathrm{\Delta }V=0`$, shows a Coulomb diamond (indicated by dotted lines) with a clearly distinguishable peak at zero bias. A non-zero $`\mathrm{\Delta }V`$ results in a current $`I_{12}`$ between source 1 and 2, and a finite voltage $`(\mu _{S2}\mu _{S1})/e`$ across the wire. Hence the quantum dot faces, at its topside, a non-equilibrium electron distribution with two quasi-Fermi energies, $`\mu _{S1}`$ and $`\mu _{S2}`$, associated with the carriers from the left and right reservoir (see Figs. 1b,c). Figure 2c shows a $`dI/dV`$ vs $`V`$ trace taken in the middle of the Kondo valley for $`\mathrm{\Delta }V=150`$ $`\mu `$V. Strikingly, the Kondo peak has split due to the applied voltage $`\mathrm{\Delta }V`$. The splitting is seen all over the Coulomb diamond as shown in Fig. 2d. A similar result is found for $`\mathrm{\Delta }V=150`$ $`\mu `$V (Figs. 2e,f). In this case, however, the split peaks are not fully resolved . As we show below, the split peaks occur at those voltages, $`V`$, for which $`\mu _{S1}`$ or $`\mu _{S2}`$ line up with $`\mu _D`$. This represents the first observation of a bias-induced splitting of the Kondo resonance.
In Fig. 3a we plot several $`dI/dV`$ vs $`V`$ traces taken again in the middle of the Kondo valley, but now for different $`\mathrm{\Delta }V`$ between -300 and 300 $`\mu `$V. To compensate for the series resistance $`R_1`$, each trace has been shifted horizonally by as much as $`R_1I_{12}=R_1\mathrm{\Delta }V/(R_1+R_2+R_{ch})`$, i.e. $`VV+0.20\mathrm{\Delta }V`$. Then $`V`$ becomes the voltage $`\mu _{S1}/e(\mu _{S1}\mu _D)/e`$ (because of the relatively small current through the drain lead we can neglect the voltage drop on its series resistance and assume $`\mu _D=0`$).

The splitting of the Kondo resonance is proportional to $`\mathrm{\Delta }V`$ as emphasized by the dotted lines in Fig. 3a. One of the split peaks stays at zero bias, corresponding to the lineup condition between $`\mu _{S1}`$ and $`\mu _D`$. The other peak moves with $`\mathrm{\Delta }V`$ following the condition $`V=(\mu _{S1}\mu _{S2})/e=R_{ch}\mathrm{\Delta }V/(R_1+R_2+R_{ch})`$. This corresponds to the alignment of $`\mu _{S2}`$ with $`\mu _D`$.
Increasing $`|\mathrm{\Delta }V|`$ (i.e. moving to the top or bottom of Fig. 3a) results in a suppression of the split Kondo resonance. The suppression is due to inelastic scattering processes in the wire which transfer electrons from the high- to the low-energy Fermi reservoir. On what voltage scale is the Kondo resonance washed out? To address this question we plot normalized peak height as a function of $`\mathrm{\Delta }\stackrel{~}{\mu }(\mu _{S2}\mu _{S1})/w_0`$ (see solid circles in Fig. 3c). Since $`ew_0k_BT_K`$, we are scaling the voltage across the wire in units of the Kondo temperature. For $`|\mathrm{\Delta }\stackrel{~}{\mu }|1`$ (i.e. $`|\mu _{S2}\mu _{S1}|k_BT_K`$) the peak height becomes half of the value for $`\mathrm{\Delta }V=0`$, which means that the Kondo resonance is suppressed on a voltage scale $`k_BT_K/e`$. For $`|\mathrm{\Delta }\stackrel{~}{\mu }|>2`$ we have fitted the data to the analytical function $`A\text{ln}^2(\alpha |\mathrm{\Delta }\stackrel{~}{\mu }|)(1+2\text{ln}^1(\alpha |\mathrm{\Delta }\stackrel{~}{\mu }|))`$, derived in Ref.  in the limit $`|\mathrm{\Delta }\stackrel{~}{\mu }|1`$. Although the latter condition is not fully satisfied we find good agreement for $`A1`$ and $`\alpha 5`$.
We now discuss the effect of a magnetic field, $`B`$, perpendicular to the 2DEG. Already for fields of order 0.1 T, the Lorentz force causes a significant shift in the electron waves traveling through the wire. (We note that the magnetic length for $`B=0.1`$ T is about 80 nm, i.e. a few times smaller than the effective width of the wire.) For $`B<0`$, left moving electrons are pushed away from the dot, whereas right movers travel closer to it (Fig. 4a). The situation is reversed for $`B>0`$ (Fig. 4b).
In Fig. 4c and 4d we show $`dI/dV`$ vs $`V`$ traces taken for $`B=0.2`$ and 0.2 T, respectively. As in Fig. 3a, $`V_g`$ is adjusted to the middle of the Kondo valley, $`\mathrm{\Delta }V`$ is varied between -300 and 300 $`\mu `$V. The voltage is compensated by $`0.20\mathrm{\Delta }V`$ for the series resistance $`R_1`$. As opposed to the results for $`B=0`$, no splitting is observed at finite $`\mathrm{\Delta }V`$. We explain such a different behavior as follows. At $`B=0.2`$ T only right movers couple effectively to the dot, while the coupling for left movers is strongly suppressed (see Fig. 4a). As a result, most of the current through the dot comes from the strongly coupled reservoir (see inset to Fig. 3c and relative caption) and only the Kondo peak corresponding to the lineup between $`\mu _{S1}`$ and $`\mu _D`$ survives. Indeed the Kondo peak stays at $`V=0`$, irrespective of $`\mathrm{\Delta }V`$. For $`B=0.2`$ T, only electrons from the right reservoir couple effectively to the dot giving rise to Kondo effect. The resulting Kondo peak occurs when $`\mu _{S2}`$ is aligned with $`\mu _D`$, and its position shifts proportionally to $`\mathrm{\Delta }V`$.
In the inset to Fig. 4d we show that the coupling of the weakly coupled lead can be restored by making the voltage $`V_{ch}`$ more negative. This reduces the effective width of the wire pushing together left and right movers. The splitting of the Kondo resonance by a finite $`\mathrm{\Delta }V`$ is recovered for $`B=0.2`$ T, similarly to the case $`B=0`$.
Finally, we compare quantitatively results for $`B=\pm 0.2`$ T and $`B=0`$. Figure 3b shows the $`V`$-positions of the Kondo peaks as obtained from Figs. 3a, 4c, and 4d. Split peaks from Fig. 3a (solid circles) match well with the single peaks from Figs. 4c and 4d (open squared and diamonds). This corroborates the interpretation given earlier that split peaks occur when $`\mu _D`$ lines up with $`\mu _{S1}`$ and $`\mu _{S2}`$.

In Fig. 3c we compare the $`\mathrm{\Delta }V`$-dependence of the peak height at zero and finite $`B`$. The suppression of the Kondo resonance becomes considerably less severe at finite $`B`$. This increased robustness stems from the displacement between left and right movers which results in a lower rate for inelastic scattering events. Negative magnetic fields appear to be more effective in this respect .
In conclusion, we have shown that a voltage difference between two Fermi seas coupled to a quantum dot causes the Kondo resonance to split, in agreement with theory. Split peaks are strongly suppressed for voltage differences much larger than $`k_BT_K/e`$. A perpendicular magnetic field has proved a useful tool to control the relative coupling of the leads. This provides a way to obtain a three-terminal quantum dot in which one of the leads is weakly coupled and can be used as a probe for the local DOS .
We thank Yu. V. Nazarov, J. Kroha, Y. Meir, L. I. Glazman, J. E. Mooij, T. Hayashi, B. van der Enden, and R. Schouten for discussions and help. This work was supported by the Specially Promoted Research, Grant-in-Aid for Scientific Research, from the Ministry of Education, Culture, Sports, Science and Technology in Japan, the Dutch Organisation for Fundamental Research on Matter (FOM), the DARPA-QUIST program (DAAD19-01-1-0659), and the EU via a TMR network (ERBFMRX CT98-0180).


text:
# Optical properties of 𝛼'-𝑁𝑎_𝑥𝑉₂𝑂₅

## I Introduction

$`\alpha ^{}`$\- NaxV2O5 belongs to the family of AV2O5 oxides (A=Li, Na, Ca, Mg) which have demonstrated, due to their peculiar crystal structures 1 ; 1a , a variety of low-dimensional phenomena . The nominally pure $`\alpha ^{}`$-NaV2O5 is a mixed valence compound ($`V^{4+}:V^{5+}=1:1`$), with the structure consisting of vanadium-oxygen (VO5) pyramids that are connected via common edges and corners to form layers in the (ab) plane. Its structure can be described as an array of parallel ladders running along the b-axis, Fig. 1. Each rung is made of a V-O-V bond, and contains one valence electron donated by the sodium atoms which are situated between layers. The sodium deficiency does not alter $`\alpha `$’ crystal structure 1b  up to x=0.70, but changes the relative abundance of the $`V^{4+}`$ and $`V^{5+}`$ ions, e.g. introduces ”extra” holes (more $`V^{5+}`$ ions) in the (ab) planes.

$`\alpha ^{}`$-NaV2O5 is one-dimensional antifferomagnetic (AF) insulator at room temperature 1b , despite the fact that the vanadium atoms have uniform valence of $`+4.5`$, which indicates the quarter-filled band structure, and suggests a metallic state. However, the optical measurements 2  revealed the electronic gap of the order of 1 eV, which is formed between linear combinations of the $`3d_{xy}`$ states of the two vanadium ions forming the rungs 3 . The insulating ground state is argued to originate from the strong electron-electron interaction 3 ; 3a  e.g. due to the presence of a large on-site Hubbard repulsion parameter U.

The effects of the hole doping on the optical properties of $`\alpha ^{}`$-NaV2O5 have not been studied in details. This interesting topic certainly deserves attention because of the spectacular effects observed in doped AF insulators, e.g. high-temperature superconductivity. Reduction of the Na concentration with respect to the nominally pure $`\alpha ^{}`$-NaV2O5 has a strong influence on the room temperature Raman spectra 4 ; 4a . The largest effect is the energy shift (about 35 cm-1) of the bond bending ($`VO_3V`$) 448 cm-1 phonon caused by the Na-deficiency. Such a strong renormalization of the phonon frequency is conjectured to arise from the delocalization of the electrons from the ladder legs to the bridge oxygens 4a .

Here we further examine the electron-phonon coupling effect in $`\alpha ^{}`$-NaxV2O5 by measuring Raman spectra using different laser line energies. We also analyze the influence of the Na-deficiency on the optical transitions in the mid-IR, and visible energy range. First, we show that the phonon energy renormalization is not as large as reported in Refs. 4 ; 4a , since the structure around 450-500 cm-1 consists in fact of two modes. Their intensity ratio changes due to the resonant conditions caused by the activation of the hole-doping associated electronic band, peaked at $`2.8`$ eV. Second, we discuss the electronic structure of $`\alpha ^{}`$-NaxV2O5 on the basis of the infrared (IR) reflectivity and ellipsometric measurements.

## II Experimental details

The present work was performed on single crystal plates of $`\alpha ^{}`$-NaxV2O5 (0.85$``$x$``$1.00) with dimensions typically about $`2\times 4\times 0.5`$ mm3 in the a, b, and c axes, respectively. The details of the sample preparation were published elsewhere 1a . The ellipsometric measurements on the (001) surfaces of the crystals, both with the plane of incidence of the light along the a and the b axis. An angle of incidence $`\theta `$ of $`66^o`$, was used in all experiments. Procedure for extraction of the dielectric constants is described in Ref. (5 ). The infrared measurements were carried out with a BOMEM DA-8 FIR spectrometer. A DTGS pyroelectric detector was used to cover the wave number region from 100 to 700 cm-1; a liquid nitrogen cooled HgCdTe detector was used from 500 to 5000 cm-1. Spectra were collected with 2 cm-1 resolution, with 1000 interferometer scans added for each spectrum. Raman spectra were measured in the backscattering configuration using micro-Raman system with DILOR triple monochromator including liquid nitrogen cooled CCD-detector. The Ar and He-Ne lasers were used as excitation sources.

## III Raman spectra

Fig.2 shows the (aa) polarized Raman spectra of $`\alpha ^{}`$-NaxV2O5, measured at room temperature, for several laser line energies. The spectra are scaled along the vertical-axis for clarity. The spectra of nominally pure $`\alpha ^{}`$-NaV2O5, and the Na-deficient spectra measured with 514.5 nm laser line show the characteristic features which were already observed, and discussed in previous reports 10 ; 10a ; 10b . Here, we focus on the 350 - 800 cm-1 energy range, where three $`A_g`$ phonon modes at 418 cm-1, 448 cm-1 and 532 cm-1, and a broad structure centered at 650 cm-1 are found. Two phonons at 448 cm-1 and 532 cm-1 are strongly interacting with underlying continuum, which is manifested through their asymmetric line shapes. In the Na-deficient samples this effect seems to be strongly enhanced 4a . The 448 cm-1 mode corresponds to the $`VO_3V`$ bond bending vibrations, whereas the 532 cm-1 mode corresponds dominantly to the bond stretching $`VO_2`$ vibrations with a small contribution of the $`VO_3V`$ bond bending vibrations 101 . The origin of the 650 cm-1 structure is still under debate. Its possible magnetic origin, due to the AF coupling of the electrons along the ladder direction (b-axis), can be ruled out since the structure is observed only in the (aa) polarized configuration, and because its energy does not fit the estimates 102  of the exchange coupling constant. The electronic band structure calculations 11 ; 11a  show no evidence for the existence of low lying electronic states even though both Raman and IR spectra confirmed their activity. The polaronic scenario, as discussed in Ref. \[ 4a \], is also questionable since the maximum of the electron density of states does not fall above the energy of the highest frequency phonon mode. Nevertheless, Raman spectra indicate the existence of the strong electron phonon coupling, that should be related to the charge transfer of the rung electron. The sodium deficiency introduces the holes in the rungs, so it is natural to conjecture that such charge dissonance is responsible for the bond bending phonon renormalization 4a .

However, the Na-deficiency has very little or no effect on the Raman spectra measured with the 638.2 nm line, see Fig.2. The $`VO_3O`$ bond bending phonon at 448 cm-1 does not exhibit large energy renormalization, and the 650 cm-1 electronic continuum is clearly visible even in the sample with the $`15\%`$ hole concentration. This finding suggests that it is necessary to include the resonant effects into consideration. For example, in Fig.2b we present the (aa) polarized room temperature Raman spectra of $`\alpha ^{}`$-Na0.95V2O5 single crystal measured with several laser line energies. This graph shows a typical evolution of the Raman spectra in $`\alpha ^{}`$-NaxV2O5 under strong resonant conditions. Thus, we fit the Raman spectra using Fano model 12 , which includes the effects of the interaction between discrete (phonon) and continuum (electron) states. For a particular energy range, 350 - 800 cm-1, we took five oscillators with 17 parameters. It is important to note, as it will be evident from our fits, that not all the parameters are crucial for the observed line shape change. The parameters of the first Lorentz oscillator at 418 cm-1 are kept constant, and used for the normalization of the spectra, since this mode does not show any change of intensity, energy, or half width under resonant conditions. For the phonons at 448 cm-1 and 532 cm-1 we adopt the Fano line shapes of the form 12 : $`I(q+ϵ)^2/(1+ϵ^2)`$, with q being the ratio of the Raman tensors for discrete vibronic and continuum electronic scattering. The parameter $`ϵ`$ is $`(\omega \omega _p)/\mathrm{\Gamma }`$, where $`\omega _p`$ and $`\mathrm{\Gamma }`$ represent the real and the imaginary part of the phonon self-energy, when its interaction with a continuum is taken into account. The continuum is represented with a Gaussian. The fifth oscillator appears at about 495 cm-1 and its difference from the 448 cm-1 mode is nicely seen in the Raman spectra of $`\alpha ^{}`$-Na0.95V2O5, measured with 488 nm laser line, see Fig.2b. It is difficult to argue whether 495 cm-1 mode interacts with the continuum or not, so we used Lorentzian for this mode. The results of the fit are shown as solid lines in Fig.2. The fits are in all cases excellent. The values of the parameters for $`\alpha ^{}`$-Na0.95V2O5 are presented in the Table I.

The major change in Raman spectra under resonant conditions is described by only four parameters; intensity (I), and q, of the 448 cm-1 mode, intensity of the 495 cm-1 mode, and the $`\omega `$ of the Gaussian. The values of the $`\mathrm{\Gamma }`$, and $`\omega _p=\omega _M\mathrm{\Gamma }/q`$ ($`\omega _M`$ is the position of the mode maximum) obtained from these fits are, within experimental error, the same for all spectra. Since the $`\mathrm{\Gamma }`$ is squared matrix element of the coupling between continuum and discrete state, this parameter should have no relationship with scattering wavelength 11 . This is to be expected for the pure resonant effect, unless there is a change in the strength of the electron-phonon interaction, caused by the sodium deficiency. However, these parameters are also found to be the same for all samples with different Na concentrations, which indicates that the electron-phonon interaction is not dependent on the hole concentration. The parameter q of the 448 cm-1 mode has a strong dependence on scattering frequency, and its variation increases with increasing sodium deficiency. Similar increase is also found for the intensity of the 495 cm-1 mode. The square of the parameter q represents the ratio of the scattering probability for discrete state to that of the continuum, hence q can exhibit a dependence on the scattering frequency if the two processes have different frequency dependencies. This is confirmed by our fits. The new mode at about 495 cm-1 is probably bond bending vibration of the rung without electron, since its energy is close to the corresponding bond bending vibration in V2O5 (482 cm-1, Ref. 12n ). Under resonance the change of these parameters (more precisely 1/q, and/or intensity) can be related to the change of the optical conductivity or dielectric constant 13 . In the nominally pure $`\alpha ^{}`$-NaV2O5, there is no electronic band in the visible energy region 2 ; 14 , and the Raman spectra show no dramatic transformation induced by the measurements with different laser line energies. However, strong resonant scattering in the Na-deficient samples should originate from the new electronic transitions at energies between 2 and 3 eV. This band, induced by hole doping is found at about 2.8 eV using ellipsometric measurements, as discussed in the next section.

## IV Ellipsometric and infrared measurements

The optical conductivity of $`\alpha ^{}`$-NaxV2O5, obtained from the ellipsometric measurements is given in Fig.3. For the nominally pure $`\alpha ^{}`$-NaV2O5 the bands at the energies 0.9, 1.2, and 3.2 eV for $`\sigma _a`$, and at 1.2, 1.6, and 3.7 eV for $`\sigma _b`$ are found in accordance with previous reports 14 ; 14a ; 14b ; 14c . The sodium deficiency causes the activation of the new optical transitions at about 2.8 eV in $`\sigma _a`$ and at 3.2 eV in $`\sigma _b`$, see Fig.3. As we have already mentioned, the 2.8 eV peak is responsible for the resonant behavior of the (aa) polarized Raman spectra. The intensity of the 495 cm-1 mode, see Inset Fig.3, and the parameter $`q^1`$ of the 448 cm-1 phonon, mimic the change in $`\sigma _a`$, which reinforces our conclusions about resonant effects in the Raman spectra.

The essential features of the $`\alpha ^{}`$-NaV2O5 electronic structure are formed from oxygen 2$`p`$ and vanadium $`3d`$ states. The fully occupied oxygen 2$`p`$ states are located about 3 eV below the lowest of the V$`3d`$ states split by ligand field 15 . The lowest state of the $`V^{4+}`$ 3$`d`$ manifold, 3$`d_{xy}`$ state, is occupied by one electron. The V-O-V rung is formed through the electron hopping via oxygen $`O_3`$ states. Due to this coupling, usually denoted as $`t_a`$, two 3$`d_{xy}`$ orbital levels form symmetric and antisymmetric combinations of levels (bonding and antibonding levels) that are split by 2$`t_a0.7`$ eV 11a . The insulating state appears due to the large on-site Hubbard repulsion, U, which pushes the Fermi level in the middle of the gap between bonding and antibonding states 3 ; 3a . Thus, the main peak in the optical spectrum along the a direction at 0.9 eV seems to originate from the transitions between bonding (B) and antibonding (AB) states of the same rung 5 . The higher energy bands, above 3 eV, are similar to those found in $`V_2O_5`$ 16 , and correspond to the $`O2pV3d`$ type of transitions.

The sodium deficiency changes the relative abundance of the $`V^{4+}`$ and $`V^{5+}`$ ions. In fact, introduction of holes leeds to the formation of rungs without electrons, similar to the $`V^{5+}O_3V^{5+}`$ bonds in V2O5. Therefore, the optical transitions above 2 eV should also be similar in V2O5 and $`\alpha ^{}`$-NaxV2O5. However, previous studies 14 ; 14a ; 14b ; 14c  did not clearly show what kind of the final state is involved in the 3.2 eV ($`\sigma _a`$) transition, and this turns out to be the most important information in order to understand the electronic structure of the Na-deficient samples. If the final state is the bonding V$`3d_{xy}`$ state, as described in Ref. 14c , the $`\alpha ^{}`$-NaV2O5 should be in the metallic state, which is in strong disagreement with the observed insulating behavior. In the insulating state, the optical transitions between $`O2p`$ and the bonding V$`3d_{xy}`$ states should be forbidden. The next possibility is that the final state of 3.2 eV peak corresponds to the nonbonding V $`3d_{xy}`$ state, which forms when there are two electrons in the rung. According to the modified Heitler-London model, proposed in order to account for the observed features in the optical conductivity spectra of $`\alpha ^{}`$-NaV2O5 5 , the nonbonding band lies in the middle of the bonding - antibonding gap. The sodium deficiency is expected to drive the Fermi level inside of the bonding bend, which now activates the transitions between the $`O2p`$ and bonding V$`3d_{xy}`$ states. In this case, the bonding level is partially filled, so this scenario also rely on the formation of the metallic state in $`\alpha ^{}`$-NaxV2O5 which is not observed.

The IR reflectivity measurements show that the material remains insulating up to the maximal achieved hole concentration of about $`15\%`$. The optical conductivity, obtained from the Kramers-Kronig (KK) analysis of the IR reflectivity, is presented in Fig.4. In both $`\sigma _a`$ and $`\sigma _b`$ spectra of $`\alpha ^{}`$-NaxV2O5 we found no Drude behavior at low frequencies. Instead, the dominant feature activated by the Na-deficiency is observed in $`\sigma _b`$ at about 0.44 eV, with a half width at half maximum (HWHM) of about 0.25 eV. Note that its energy corresponds to the half of the energy of 0.9 eV peak.

Clearly, all scenarios discussed so far are not yet adequate to fully describe the optical transitions and the electronic structure of the $`\alpha ^{}`$-NaxV2O5. We argue below that the electronic structure based on the four electronic bands can account for the major features observed in the optical conductivity of $`\alpha ^{}`$-NaxV2O5. Simplified schematic representation of the optical transitions, and the energy levels in $`\alpha ^{}`$-NaV2O5, is shown in the Inset of Fig.3. In the nominally pure $`\alpha ^{}`$-NaV2O5, the arrows in the left side of the Inset of Fig.3 represent the optical transitions which we have already discussed. The energy levels of the sodium deficient samples are shown in the right side of the Inset of Fig.3. In order to account for the appearance of the new modes at 2.8 eV and 0.44 eV in $`\alpha ^{}`$-NaxV2O5 spectra, we have to assume the existence of the impurity band (IB), which forms inside the electronic gap of the pure $`\alpha ^{}`$-NaV2O5. This band is further on considered to be empty, so the optical transitions, denoted as processes ”3” in the Fig.3, may appear in the spectra. We find strong support for our assumption in the electronic structure, and the optical spectra 16  of V2O5. In this compound very narrow electronic band, which is separated from the rest of the conduction band by an additional gap, exists at the energy of about 2.5 eV. Thus, it is possible that impurity band in sodium-deficient samples has its origin from this level. If so, the optical process ”3” corresponds to the transition between oxygen 2p band and IB. Its energy, $`E_3=E_21/2E_12.8`$ eV, indicates that the IB forms almost exactly in the middle of the electronic gap of the pure $`\alpha ^{}`$-NaV2O5. This is further supported by the existence of the mid-IR absorption peak at about $`E_{IB}=1/2E_10.44`$ eV, see Fig.4, which corresponds to the transitions between V $`d_{xy}`$ and IB states. The nature of the IB determines the electronic properties of $`\alpha ^{}`$-NaxV2O5. By increasing sodium deficiency the 0.44 eV peak grows in intensity and its width increases, which indicates that the IB gains dispersion. In order to preserve the preferential insulating ground state (at least up to 15 $`\%`$ deficiency) this process is accompanied with the increase of the gap between vanadium $`3d`$ bands. This is evident from the energy increase of the peak ”1” in $`\sigma _a`$. Because of that, and since this band is indeed very narrow (0.25 eV), IB does not overlap with the occupied dxy band, and the solution is an insulator since the Fermi level is placed in the middle of the gap between bonding $`d_{xy}`$ states and IB. The 0.44 eV peak is found to be much stronger in $`\sigma _b`$ then in $`\sigma _a`$ spectra. Thus, we tentatively assigned it as the inter-rung electron transition along the b axis.

## V Conclusion

In conclusion, we present the analysis of the optical properties of $`\alpha ^{}`$-NaxV2O5 using ellipsometric, IR and Raman scattering techniques. Despite strong hole-doping the material remains insulating up to the doping concentration of about 15 $`\%`$. The insulating state appears as a consequence of very narrow impurity band formed in the middle of the electronic gap of the pure $`\alpha ^{}`$-NaV2O5. The optical absorption peak develops at $`0.44`$, eV and corresponds to the transitions between vanadium V $`d_{xy}`$ and impurity band states. The change of the Raman active multiple phonon structure upon doping, in the energy range between 450 and 550 cm-1, is the consequence of the resonance, and is caused by the activation of the hole-doping associated optical transitions peaked at about $`2.8`$ eV .

Acknowledgments

M.J.K. and Z.V.P. acknowledge support from the Research Council of the K.U. Leuven and DWTC. The work at the K.U. Leuven is supported by Belgian IUAP and Flemish FWO and GOA Programs.


text:
# Explicit finite-difference and direct-simulation-MonteCarlo method for the dynamics of mixed Bose-condensate and cold-atom clouds

## I Introduction

The physics of ultracold atomic vapours under magnetic or optical confinement has been a continuing and ever expanding focus of interest since the realization of Bose-Einstein condensation . Following characterizations of the basic thermodynamic and dynamical properties of condensates , a number of experiments have been addressed to investigations of their phase coherence and superfluidity , to the study of non-linear effects and special spectroscopic features , and to the observation of vortices . Parallel efforts are being made in the study of gases of fermionic atoms  and of boson-fermion mixtures , with the ultimate aim of realizing novel superfluids.

Theoretical studies of the dynamics of these systems, involving the analytical solution of approximate models have very often been successful in explaining or predicting such novel phenomena. However, the interplay of different species or the thermal fluctuations of a condensate are not easily handled by analytical methods. The validity of the approximations may be limited to the extreme collisionless or hydrodynamic regimes, the confinement of the sample being treated within a local-density approach. Thus, while the equations governing these dilute systems remain simple, their numerical solution can be helpful for investigating more complex dynamical problems where an intermediate regime is met or a cold-atom cloud accompanies the condensate.

The atomic interactions in a highly dilute Bose gas at very low temperatures, as is relevant for current experiments, are described by a contact pseudopotential accounting for $`s`$-wave scattering in binary atom-atom elastic collisions. A condensate is then treated within a mean-field approach by solving the stationary or the time-dependent Gross-Pitaevskii equation (GPE). Several types of numerical approaches have been developed: eigenvalue solvers, variational solvers, or explicit solvers  for the ground state, and implicit  or explicit  time-marching schemes for the dynamics.

Methods for studying the dynamics of an isolated cloud of ultracold (bosonic or fermionic) atoms are also well developed. One needs to solve the Vlasov-Landau equation of motion (VLE) for the Wigner distribution function. Various numerical techniques have been used for this purpose, which are based either on an ergodic assumption  or on the inclusion of statistical noise , or else use a direct simulation MonteCarlo (DSMC) approach borrowed from the methods of molecular dynamics .

The general theoretical background is provided by the book of Kadanoff and Baym . These authors developed a Green’s function approach to transport phenomena, which extends the Boltzmann equation to strongly interacting quantum fluids and allows for progressively improved self-consistent approximations. This formalism was extended to a homogeneous Bose-condensed gas at finite temperature by Kane and Kadanoff , within a Beliaev approximation including interactions up to second order. More recently, these methods have been adapted to the theory of transport phenomena in a confined Bose gas within the Hartree-Fock-Bogoliubov approximation , dealing with a Bose-Einstein condensate accompanied by its thermal cloud. Jackson and Adams  have proposed to combine the GPE with a quantum version of the DSMC to numerically evaluate the dynamics of such a fluid. Numerical studies based on a generalized GPE combined with a semiclassical kinetic equation, and including collisions between the condensate and the thermal cloud, are also becoming avaible for dynamical properties of a trapped Bose-condensed gas at finite temperature .

In the present work we proceed along the path traced by Jackson and Adams . We propose a different approach to the solution of the GPE and a different method for preparing the initial equilibrium state, which would be immediately applicable to a multi-component cold-atom cloud. The method is applied to two classes of dynamical problems: the first concerns the ballistic expansion of a fermion cloud and the role played by the presence of a condensate, while the second concerns a Bose-condensed gas in a periodic optical-lattice potential at finite temperature.

After introducing in Sec. II the model for both the equilibrium state and the dynamical evolution of the fluid, we describe in Sec. III the numerical methods that we have used to consistently solve the GPE for the condensate and the VLE for the Wigner distribution function of the cloud. The procedure followed in the actual computations is also outlined in Sec. III, with due emphasis on the preparation of the initial equilibrium state. The physical applications that we have carried out are presented and discussed in Sec. IV. A discussion of computational aspects in Sec. V and some final remarks in Sec. VI conclude the paper.

## II The model: a review of the main concepts

### A Transport in a normal quantum fluid

In the original formulation by Kadanoff and Baym (KB), the problem of transport in a quantum fluids in the normal state is tackled by deriving an analogue of the Boltzmann equation for the Wigner distribution function $`f(𝐩,𝐫,t)`$ from the microscopic equations of motion for the non-equilibrium density matrix $`\rho (𝐫_1,t_1;𝐫_1^{}t_1^{})<\psi _𝒰^{}(𝐫_1,t_1)\psi _𝒰(𝐫_1^{},t_1^{})>`$, which is defined through the particle creation and destruction operators $`\psi _𝒰^{()}`$ in the presence of an external, slowly varying disturbance $`𝒰(𝐫,t)`$. Namely,

$$f(𝐩,𝐫,t)=𝑑𝐱\mathrm{exp}(i𝐩𝐱)<\psi _𝒰^{}(𝐫+𝐱/2,0)\psi _𝒰(𝐫𝐱/2,t)>$$
(1)
where $`𝐫=(𝐫_1+𝐫_1^{})/2`$ and $`𝐱=𝐫_1𝐫_1^{}`$ are the center-of-mass and the relative coordinate of the two particles. The moments of the Wigner function yield observables such as the particle density $`n(𝐫,t)=(2\pi )^3𝑑𝐩f(𝐩,𝐫,t)`$ and the current density $`𝐣(𝐫,t)=(2\pi )^3𝑑𝐩(𝐩/m)f(𝐩,𝐫,t)`$.

Contact with the Boltzmann transport equation is made by performing gradient expansions. As in the conventional Boltzmann-equation approach, the validity of the KB formulation is limited to slowly varying perturbations. On the other hand, the advantage of the KB formulation is that higher-order correlations enter the equation of motion for the density matrix in an explicit manner, and therefore systematically improved approximations which are consistent with the conservation laws are accessible. Examples of such treatments of the correlation term are the Hartree approximation and the Born-collision approximation. In the former case the collisionless Boltzmann equation is recovered, while in the latter the collisional Boltzmann equation is extended to non-dilute systems by including the effect of the external potential on the motion of the particles between collisions.

### B Extension to coupled condensate-noncondensate dynamics

The extension of the KB treatment to gases including a Bose-condensed component has been made by Kane and Kadanoff  and further developed by Griffin and coworkers  and also by Wachter et al.  through a different derivation. The presence of two components and the appearence of off-diagonal elements in the density matrix (the so-called anomalous densities) below the Bose-Einstein condensation temperature requires the introduction of three Wigner distribution functions: $`f_c(𝐩,𝐫,t)`$ for the condensate component described by $`<\psi _𝒰^{()}>`$ and involving $`|<\psi _𝒰^{}>|^2`$, $`f_b(𝐩,𝐫,t)`$ for the noncondensate described by the fluctuations operators $`\stackrel{~}{\psi }_𝒰^{()}\psi _𝒰^{()}<\psi _𝒰^{()}>`$ and involving $`<\stackrel{~}{\psi }_𝒰^{}\stackrel{~}{\psi }_𝒰>`$, and $`f_m(𝐩,𝐫,t)`$ for the anomalous part involving $`<\stackrel{~}{\psi }_𝒰\stackrel{~}{\psi }_𝒰>`$ and its Hermitean conjugate.

We thus have to deal with the density of condensate $`n_c(𝐫,t)=(2\pi )^3𝑑𝐩f_c(𝐩,𝐫,t)`$, the density of noncondensate $`n_b(𝐫,t)=(2\pi )^3𝑑𝐩f_b(𝐩,𝐫,t)`$, and the anomalous density $`m(𝐫,t)=(2\pi )^3𝑑𝐩f_m(𝐩,𝐫,t)`$. Analogous expressions holds for the current densities.

As to the consistency of the approximations with the conservation laws, the same general remarks as for normal systems apply. However, the appearence of the condensate introduces an additional principle of gauge invariance, leading to the requirement that the excitation spectrum to be gapless . It is well known  that approximations capable of simultaneously accommodating the conservation laws and the gaplessness condition are hardly available, so that a choice has to be made depending on the specific conditions of density and temperature of the system.

In the regime that we address in the present work the anomalous densities can be neglected, resulting in the gapless Hartree-Fock-Bogoliubov-Popov approximation . Thus, the equation of motion for the condensate wavefunction $`\mathrm{\Phi }(𝐫,t)<\psi _𝒰^{}(𝐫,t)>`$ is

$$i\mathrm{}\frac{\mathrm{\Phi }(𝐫,t)}{t}=\left[\frac{\mathrm{}^2}{2m}\mathbf{}^2+V_{c,b}^{eff}(𝐫,t)\right]\mathrm{\Phi }(𝐫,t)$$
(2)
and is coupled to the collisionless Vlasov equation for the noncondensate Wigner function $`f_b(𝐩,𝐫,t)`$,

$$\frac{f_b(𝐩,𝐫,t)}{t}+\frac{𝐩}{m}\mathbf{}_𝐫f_b(𝐩,𝐫,t)\mathbf{}_𝐫V_b^{eff}(𝐫,t)\mathbf{}_𝐩f_b(𝐩,𝐫,t)=0.$$
(3)
The mean-field potentials in Eqs. (2) and (3) are

$$V_{c,b}^{eff}(𝐫,t)=V_b^{ext}(𝐫)+𝒰(𝐫,t)+U_g[n_c(𝐫,t)+2n_b(𝐫,t)]$$
(4)
and

$$V_b^{eff}(𝐫,t)=V_b^{ext}(𝐫)+𝒰(𝐫,t)+2U_g[n_c(𝐫,t)+n_b(𝐫,t)],$$
(5)
including the time-dependent driving potential $`𝒰(𝐫,t)`$ and an axially symmetric confining potential given for a harmonic trap by $`V_b^{ext}(𝐫)=m_b\omega _b^2(r_{}^2+ϵ_b^2z^2)/2`$. In Eqs. (4) and (5), $`n_c(𝐫,t)=|\mathrm{\Phi }(𝐫,t)|^2`$ and $`U_g=4\pi \mathrm{}^2a_{bb}/m_b`$ is the boson-boson interaction parameter in terms of the $`s`$-wave scattering length $`a_{bb}`$ and of the boson mass $`m_b`$.

Once the algorithm to solve Eqs. (2) and (3) is implemented, it is easily extended to a mixture of condensed bosons and a fermionic cloud in the collisionless regime. In this case, the effective mean-field potentials become

$$V_{c,f}^{eff}(𝐫,t)=V_b^{ext}(𝐫)+𝒰(𝐫,t)+U_gn_c(𝐫,t)+U_fn_f(𝐫,t)$$
(6)
and

$$V_f^{eff}(𝐫,t)=V_f^{ext}(𝐫)+𝒰_f(𝐫,t)+U_fn_c(𝐫,t).$$
(7)
Here, $`V_f^{ext}(𝐫)=m_f\omega _f^2(r_{}^2+ϵ_f^2z^2)/2`$ and $`𝒰_f(𝐫,t)`$ are the external trapping and driving potentials acting on the fermions and $`U_f=2\pi \mathrm{}^2a_{bf}/m_r`$ is the boson-fermion coupling constant with $`a_{bf}`$ the boson-fermion $`s`$-wave scattering length and $`m_r=m_bm_f/(m_b+m_f)`$, $`m_f`$ being the fermion atomic mass. Notice that fermion-fermion collisions in the $`s`$-wave channel are effectively suppressed by the Pauli principle in a dilute gas of spin-polarized fermions, as is relevant to current experiments on boson-fermion mixtures.

### C Validity of the model

To summarize, Eqs. (2) and (3) describe the coupled dynamics of a Bose-Einstein condensate and a bosonic or fermionic cold-atom cloud. In the former case the potentials $`V_{c,b}^{eff}`$ and $`V_b^{eff}`$ are used, and in the latter these are replaced by $`V_{c,f}^{eff}`$ and $`V_f^{eff}`$. The range of validity of this approach is in principle limited to: (i) slowly varying space- and time-dependent external potentials, which allows a low-order gradient expansion of the equations of motion for the one-body density matrix; (ii) a collisionless regime, which allows expansion of the self-energies to first order in the strength of the atomic interactions; and (iii) not too low temperature, so that the anomalous averages can be neglected.

We implement below a numerical method to solve Eqs. (2) and (3) and test it against dynamical behaviors in a one-component fermionic cloud and in clouds of either fermions or thermal bosons accompanied by a condensate. Of course, the role of the statistics enters at this level only from the initial distribution of the particles in phase space. We thus turn to discuss the equations that we use to determine the initial conditions for the subsequent time evolution.

### D Semiclassical description of the equilibrium state

Before proceeding to present the algorithm for our dynamical simulations we briefly recall the basic steps that we take in preparing the initial state of the gas in thermodynamic equilibrium and in evaluating the corresponding densities of the condensate and of the fermionic or bosonic cloud. We refer the reader to Ref.  for the details of the theory and for a discussion of the excellent agreement that it yelds with thermodynamic data on Bose-Einstein condensed gases, under conditions of temperature and dilution that will be verified in the calculations reported in Section IV below.

The equilibrium condensate density is calculated within the Thomas-Fermi approximation, which amounts to neglecting the kinetic energy term in the GPE. Its validity is ensured whenever the average mean-field energy $`U_gn_c`$ is much larger than the typical confining energy. In the case of harmonic confinement the condition $`N_ca_{bb}/a_{ho}1`$ is required, where $`N_c`$ is the number of atoms in the condensate, $`a_{ho}=(\mathrm{}/m_b\overline{\omega }_b)^{1/2}`$ is the harmonic oscillator length and $`\overline{\omega }_b=\omega _bϵ^{1/3}`$ is the geometric average of the trap frequencies. The equilibrium density profile of the condensate is given by

$$n_c(𝐫)=\frac{1}{U_g}\left[\mu _bV_b^{ext}(𝐫)k_{b,f}n_{b,f}(𝐫)\right]\theta (\mu _bV_b^{ext}(𝐫)k_{b,f}n_{b,f}(𝐫))$$
(8)
where $`k_b=2U_g`$, $`k_f=U_f`$ and $`\mu _b`$ is the chemical potential for the bosons.

For the equilibrium cloud density of bosons or fermions we adopt the semiclassical Hartree-Fock scheme . This choice is justified as long as the gas is in a very dilute regime, a condition which is usually met in current experiments. In this approximation we have

$$n_{b,f}(𝐫)=\frac{\mathrm{d}𝐩}{(2\pi )^3}\left\{\mathrm{exp}\left[\frac{1}{k_BT}\left(\frac{p^2}{2m_{b,f}}+V_{b,f}^{eff}(𝐫)\mu _{b,f}\right)\right]1\right\}^1,$$
(9)
with $`V_b^{eff}(𝐫)`$ and $`V_f^{eff}(𝐫)`$ determined by the confining potentials supplemented by static mean-field interaction term as in Eqs. (5) or (7).

The chemical potentials $`\mu _b`$ and $`\mu _f`$ are determined from the total numbers of bosons and fermions. In the case of a bosonic thermal cloud, $`\mu _b`$ is fixed by the relation

$$N_b=d𝐫\left[n_c(𝐫)+n_b(𝐫)\right].$$
(10)
For a fermionic cloud the chemical potential is determined from the total number of fermions,

$$N_f=d𝐫N_f(𝐫).$$
(11)
These equations complete the self-consistent closure of the model in the initial equilibrium state.

## III The numerical method

In this section we present the numerical procedure that we have used to solve the system of Vlasov-Landau and Gross-Pitaevskii equations. Since most experimental setups are invariant under rotation in the azimuthal plane, we use cylindrical coordinates $`\{r,z\}`$. The wavefunction $`\mathrm{\Phi }`$ is discretized on a two-dimensional grid of $`N_r\times N_z`$ points, which are uniformly distributed in a box of sizes $`r_{max}\times 2z_{max}`$: that is, $`\mathrm{\Phi }_{jk}\mathrm{\Phi }(r_j,z_k)`$ with

$$\{\begin{array}{cc}r_j=(j1)\mathrm{\Delta }r\hfill & \hfill (j=1,\mathrm{},N_r)\\ z_k=z_{max}+(k1)\mathrm{\Delta }z\hfill & \hfill (k=1,\mathrm{},N_z)\end{array}$$
(12)
$`\mathrm{\Delta }r`$ and $`\mathrm{\Delta }z`$ being the steps in the two space variables. The particle distributions are discretized by means of a set of $`P`$ computational particles,

$$f(𝐩,𝐫,t)f_P(𝐩,𝐫,t)\underset{i=1}{\overset{P}{}}\delta (𝐫𝐫_i(t))\delta (𝐩𝐩_i(t))$$
(13)
where the $`6P`$ phase-space coordinates $`𝐫_i(t)`$ and $`𝐩_i(t)`$ represent the actual numerical unknowns entering the VLE. They obey Newton’s equations,

$$\{\begin{array}{cc}\frac{d𝐫_i}{dt}=𝐩_i/m\hfill & \\ \frac{d𝐩_i}{dt}=𝐅_i\hfill & \end{array}$$
(14)
where $`𝐅_i`$ collects all forces acting on the $`i`$-th computational particle. As a result, at each time step we have to solve for $`N_G=N_rN_z`$ grid unknowns coupled to $`6P`$ discrete particle coordinates. In order to keep the systematic signal reasonably above the noise level, each grid cell contains of the order of $`10100`$ computational particles. As a result, $`P10N_G100N_G`$.

### A Propagation of the Gross-Pitaevskii equation

The GPE is advanced in time by an explicit finite-difference method based on a non-staggered variant of the Visscher method . The full details are given in original papers  and here we shall just review the main ingredients of the algorithm. The basic idea is to advance the real and imaginary parts of the wavefunction, say $`A`$ and $`B`$, in alternating steps as follows:

$$\{\begin{array}{cc}A_{jk}^{n+1}A_{jk}^{n1}=2[K_{jk}(B)+V_{jk}^n(B)]\mathrm{\Delta }t\hfill & \\ B_{jk}^{n+1}B_{jk}^{n1}=2[K_{jk}(A)+V_{jk}^n(A)]\mathrm{\Delta }t\hfill & \end{array}$$
(15)
where $`n=1,\mathrm{},N_T`$ labels the discrete time sequence $`t_n=n\mathrm{\Delta }t`$. The scheme is initiated as follows. The initial conditions specify the values of $`A`$ and $`B`$ at $`n=0`$ and subsequently a first-order Euler step provides their values at $`n=1`$. With these values available, all later steps $`n=2,3,\mathrm{},N_T`$ are taken by using the above set of equations.

In Eq. (15) $`K_{jk}`$ is the kinetic energy operator,

$$K_{jk}(\mathrm{\Phi })=\frac{\mathrm{}^2}{2m}[\frac{\mathrm{\Phi }_{j1,k}}{\mathrm{\Delta }r^2}+\frac{\mathrm{\Phi }_{j+1,k}}{\mathrm{\Delta }r^2}+\frac{\mathrm{\Phi }_{j,k1}}{\mathrm{\Delta }z^2}+\frac{\mathrm{\Phi }_{j,k+1}}{\mathrm{\Delta }z^2}2(\frac{1}{\mathrm{\Delta }r^2}+\frac{1}{\mathrm{\Delta }z^2})\mathrm{\Phi }_{j,k}]$$
(16)
and

$$V_{jk}(\mathrm{\Phi })=V^{eff}(r_j,z_k;n_{jk})\mathrm{\Phi }_{jk}$$
(17)
is the potential energy operator, including both external and self-consistent interaction terms. The self-consistent potential requires the specification of the particle density $`n_{jk}`$ at each node of the spatial grid. This is obtained by convoluting the discrete particle distribution,

$$n_{jk}(t)=\underset{iC_{jk}}{}W_{jk,i}f_i(t)$$
(18)
where $`C_{jk}`$ is the grid cell centered in $`r_{j+1/2}=(r_j+r_{j+1})/2`$, and $`z_{k+1/2}=(z_k+z_{k+1})/2`$, while $`f_i=1`$ if the particle belongs to $`C_{jk}`$ and $`f_i=0`$ otherwise. The factor $`W_{jk,i}`$ weights the contribution of particle $`i`$ to the density at the grid-point $`(r_j,z_k)`$.

We adopt a bilinear ”cloud-in-cell” (CIC) interpolator (see fig. 1a), which yields the scattering (particle-to-grid) rule

$$W_{jk,i}=e_i(r_j)e_i(z_k)(j=j_i,j_i+1,k=k_i,k_i+1).$$
(19)
Here $`(j_i,k_i)`$ identifies the lower-left corner of the grid cell to which the $`i`$-th particle belongs. In Eq. (19) $`e_i(r)`$ and $`e_i(z)`$ are one-dimensional piecewise linear splines centered on the discrete particle position $`(r_i,z_i)`$,

$$e_i(r)=\{\begin{array}{c}\hfill 1+\frac{(rr_i)}{\mathrm{\Delta }r}(r_i\mathrm{\Delta }r<r<r_i)\\ \hfill 1\frac{(rr_i)}{\mathrm{\Delta }r}(r_i<r<r_i+\mathrm{\Delta }r)\end{array}$$
(20)
and similarly for $`e_i(z)`$.

### B Propagating the Vlasov-Landau equation

The VLE is advanced in time by a standard Particle-in-Cell (PIC) method . In a modified Verlet time-marching scheme, we obtain the following set of discrete algebraic equations:

$$\{\begin{array}{ccc}r_i(t+\mathrm{\Delta }t)=r_i(t)+v_{ri}(t)\mathrm{\Delta }t+a_{ri}(t)\mathrm{\Delta }t^2/2\hfill & & \\ (r_i\theta _i)(t+\mathrm{\Delta }t)=(r_i\theta _i)(t)+v_{\theta i}(t)\mathrm{\Delta }t+a_{\theta i}(t)\mathrm{\Delta }t^2/2\hfill & & \\ z_i(t+\mathrm{\Delta }t)=z_i(t)+v_{zi}(t)\mathrm{\Delta }t+a_{zi}(t)\mathrm{\Delta }t^2/2\hfill & & \\ v_{ri}(t+\mathrm{\Delta }t)=v_{ri}(t)+a_{ri}(t)\mathrm{\Delta }t\hfill & & \\ v_{\theta i}(t+\mathrm{\Delta }t)=v_{\theta i}(t)+a_{\theta i}(t)\mathrm{\Delta }t\hfill & & \\ v_{zi}(t+\mathrm{\Delta }t)=v_{zi}(t)+a_{zi}(t)\mathrm{\Delta }t\hfill & & \end{array}$$
(21)
Here $`a_r,a_\theta ,a_z`$ and $`v_r,v_\theta ,v_z`$ are the three components of the acceleration and of the velocity along the $`(r,\theta ,z)`$ coordinates. The advantage of the modified Verlet time marching is that fourth-order accuracy is preserved while synchronously keeping the coordinate and momentum degrees of freedom on the same sequence of discrete times.

The algorithm is standard except for the specification of the self-consistent coupling, namely the force due to the density gradients. In the first place, we form density gradients from auxiliary values of the density field $`n_{j+1/2,k+1/2}`$ at cell centers,

$$\{\begin{array}{cc}Gr_{ij}=[n_{j+1/2,k1/2}n_{j1/2,k1/2}+n_{j+1/2,k+1/2}n_{j1/2,k+1/2}]/(2\mathrm{\Delta }r)\hfill & \\ Gz_{ij}=[n_{j1/2,k+1/2}n_{j1/2,k1/2}+n_{j+1/2,k+1/2}n_{j1/2,k1/2}]/(2\mathrm{\Delta }z)\hfill & \end{array}$$
(22)
The azimuthal acceleration is zero in cylindrical symmmetry. Next, the grid-forces are evaluated on the discrete particle locations, this being the inverse of the scattering operation discussed in the previous section (see fig. 1b). The grid-to-particle convolution is

$$𝐅_i=\underset{s=0}{\overset{1}{}}W_{i,(j+s,k+s)}𝐅_{j+s,k+s}.$$
(23)
In order to avoid spurious self-forces we use again CIC-interpolation, which amounts to using the same weighting function as for the GPE in the grid-to-particle scattering rule,

$$W_{i,jk}=e_j(r_i)e_k(z_i).$$
(24)
With the force/acceleration field transferred to the particle locations, everything is set to march the VLE in time.

### C Boundary conditions

The conditions imposed on the wavefunction are (i) periodicity along the $`z`$ coordinate, (ii) symmetry (zero radial gradient) at $`r=0`$, and (iii) vanishing at the outer radial boundary $`r=R`$. For the discrete particle distribution we have again periodicity along $`z`$ and specular reflection at the outer radial boundary. Specular reflection means that a particle flying from, say, $`r_\alpha <R`$ to $`r_\beta >R`$ is replaced by a particle at $`r=2Rr_\beta `$ with inverted radial speed. Since the particle trajectories are tracked in a three-dimensional cylindrical coordinate frame, the $`r=0`$ axis requires no special treatment.

#### 1 Time-step considerations

The GPE and the VLE are advanced on the same discrete time sequence. This maximizes simplicity, but implies that the time-step is controlled by the fastest process at work, which usually is the self-consistent potential acting upon the condensate wavefunction. Better efficiency can be achieved by sub-cycling the time-stepper, namely by advancing the slowest equation (say the VLE) only every $`\mathrm{\Delta }t_{VL}/\mathrm{\Delta }t_{GP}`$ steps, $`\mathrm{\Delta }t_{VL}`$ and $`\mathrm{\Delta }t_{GP}`$ being the largest time-steps allowed by the stability conditions on the two equations. The maximum time-step for the GPE solver is estimated from

$$\mathrm{\Delta }t_{GP}(C_1\frac{\mathrm{}}{m\delta ^2}+C_2\mathrm{}\frac{V_M}{\mathrm{}})<1$$
(25)
where $`\delta `$ is a typical mesh size, $`V_M`$ is the maximum value of the potential and $`C_1`$ and $`C_2`$ are two coefficients $`O(1)`$ which depend on geometry and dimensionality. We note here the concurrent effects of quantum diffusion (kinetic energy) and scattering/absorption (potential energy).

The maximum time-step for the VLE solver is estimated from

$$\mathrm{\Delta }t_{VL}\frac{v_{max}}{\delta }<1$$
(26)
where $`v_{max}`$ is the maximum speed in the velocity grid, of the order of the Fermi velocity for fermions and of the thermal velocity for bosons. Under ordinary conditions the kinetic energy contributions dominate over those from the potential energy, so that the condition for the GPE to be the time-limiting section of the code takes the form of a numerical “uncertainty principle”, $`mv_{max}\delta >\mathrm{}`$. For the cases discussed in this work, this inequality is generally fulfilled within a factor ten, so that sub-cycling is not compulsory. The inclusion of collisional interactions would make it mandatory.

#### 2 Radial singularity

A source of potential trouble is the singularity in $`r=0`$, which is known to affect all calculations in polar coordinates. To date, singular factors $`1/r`$ are regularized by a simple numerical cutoff $`1/r1/(r+r_c)`$ with $`r_c0.001\mathrm{\Delta }r`$, with a check that the physical results are virtually insensitive to the specific value of $`r_c`$.

Another undesirable side-effect of the cylindrical geometry is the relative depletion of near-axis cells, which tend to host fewer particles just because of the $`r\mathrm{\Delta }r`$ volumetric effect. On the other hand, this volumetric depletion is often more than compensated by the physical behavior of the radial density, which is generally largest at $`r=0`$. At any rate, the volumetric effect can be readily disposed of by moving to a non-uniform grid along the radial coordinate.

#### 3 Statistical noise

Considerations of statistical accuracy require of the order of a few tens of particles per grid-point (or equivalently per grid-cell) in order to keep the noise-to-signal ratio below an acceptable threshold. A practical consequence of this statistical accuracy requirement is that the VLE part of the computational scheme should be designed in such a way as to evolve these tens of computational particles in approximately the same amount of CPU time that it takes the GPE solver to advance a single grid-point.

Another interesting consequence is that - at variance from ordinary situations in (classical) rarefied-gas dynamics - the number of computational particles in the simulation of Bose-Einstein condensates far exceeds the number of physical atoms, typically by a factor $`10^3`$ in our case. As a result, each single computer simulation performs de facto a built-in ensemble average over a set of about a thousand realizations.

### D Procedure: preparing the initial state

The initial condition for the populations of bosons or fermions in the cold-atom cloud is prepared as follows. Particles are sampled from the probability distribution functions (pdf)

$$f_{b,f}=\frac{e^\eta }{1e^\eta }$$
(27)
where $`\eta \beta (\mu _{b,f}V_{eff}+p^2/2m)`$ with $`\beta =1/k_BT`$ the inverse temperature.

The initializiation procedure starts by assigning to each spatial cell centered about position $`𝐫`$ a corresponding amount of computational particles,

$$\mathrm{\Delta }N(𝐫)=n_{b,f}(𝐫)\mathrm{\Delta }V(𝐫)$$
(28)
where $`\mathrm{\Delta }V(𝐫)=2\pi r\mathrm{\Delta }r\mathrm{\Delta }z`$ is the cell volume. These particles must be sampled in momentum space according to the pdf (27). Owing to the non-separability of the pdf, straightforward sampling based on exact inversion is ruled out and more general - but less efficient - accept/reject methods must be resorted to.

The particle momenta are sampled from the distribution function using the standard Box-Mueller algorithm in three dimensions for cylindrical coordinates,

$$\{\begin{array}{c}p_{}=p_{max}r_1\hfill \\ \xi =2\pi r_2\hfill \\ p_z=p_{max}(2r_31)\hfill \end{array}$$
(29)
where $`r_1`$, $`r_2`$ and $`r_3`$ are sampled from a uniform distribution in the range $`[0,1]`$. The maximum momentum $`p_{max}`$ is taken to be of order $`2p_F`$ for cold fermions and of order $`3m_bv_{th}`$ for bosons, with $`p_F=\mathrm{}k_f`$ being the Fermi momentum and $`v_{th}`$ the average thermal velocity. The particle momentum coordinates in the azimuthal plane are evaluated as

$$\{\begin{array}{ccc}p_r\hfill & =\hfill & p_{}sin\xi \hfill \\ p_\theta \hfill & =\hfill & p_{}cos\xi .\hfill \end{array}$$
(30)
Finally, in order to avoid poor acceptance rates, the standard accept/reject test is performed by comparing the pdf with the maximum value $`f_{max}(𝐫)`$ that it can take in the cell at position $`𝐫`$, that is

$`Iff(p_r,p_\theta ,p_z;𝐫)>r_4f_{max}(𝐫):accept`$ (31)
$`Else:reject`$ (32)
where $`r_4`$ is a random number uniformly distributed in $`[0,1]`$.

## IV Physical applications

### A Expansion of cold fermions

As a first test of the dynamical algorithm, we consider the expansion of a cloud of cold fermions after release of the harmonic trap. The neglect of the collisional integrals is a reliable approximation in this context, since as already remarked the $`s`$-wave scattering between spin-polarized fermions is suppressed by the Pauli principle.

This problem has been analytically solved in Ref.  under the assumption of ballistic expansion. The time-evolution of the mean square radii is found to be

$$<r_{}^2(t)>=\frac{1}{3N_f}E_{rel}\frac{4}{m_f\omega _f^2}(1+\omega _f^2t^2)$$
(33)
and

$$<z^2(t)>=\frac{1}{3N_f}E_{rel}\frac{2}{m_f(ϵ_f\omega _f)^2}(1+(ϵ_f\omega _f)^2t^2).$$
(34)
In Eqs. (33) and (34) $`E_{rel}`$ is the so-called release energy, namely the energy of the system with $`N_f`$ fermions after switching off the trap, which amounts to one half of the total average energy. For a non-interacting fermion gas at temperature $`T>0.2T_F`$ this is best approximated by the classical relation $`E_{rel}3N_fk_BT/2`$, while at lower temperature it is given by $`E_{rel}(3/4)N_fE_F[1+(2\pi /3)(T/T_F)^2]`$, with $`E_F=k_BT_F=(6N_f)^{1/3}\mathrm{}\omega _f`$ being the Fermi energy.

We prepare the initial thermodynamic state (9) for $`N_f=1000`$ $`{}_{}{}^{40}K`$ atoms in an isotropic trap with $`\omega _f=2\pi \times 15.92`$ rad/s at $`T=0.55T_F7.6`$ nK. The chemical potential of the gas is $`\mu _f=0.063E_F=1.14\mathrm{}\omega _f`$. We use a $`N_r\times N_z`$ mesh of $`201`$x$`401`$ points with $`1.610^6`$ representative particles in a box of sizes $`40`$ and $`80`$ in units of $`a_{ho}=\sqrt{\mathrm{}/(m_f\omega _f)}`$ along the radial and axial directions, respectively. The time-step in the dynamical simulation is $`\omega _f\mathrm{\Delta }t=10^4`$.

We then evaluate from the simulation runs the radial width

$$\sigma _r_{}(t)=\sqrt{𝑑x𝑑y𝑑z[(x(t)<x(t)>)^2+(y(t)<y(t)>)^2]n_f(𝐫;t)}$$
(35)
and the axial width

$$\sigma _z(t)=\sqrt{𝑑x𝑑y𝑑z(z(t)<z(t)>)^2n_f(𝐫;t)}$$
(36)
of the cloud as functions of time, after averaging over the density distribution $`n_f(𝐫;t)`$. Here, the center-of-mass coordinates are defined as $`<x(t)>=𝑑x𝑑y𝑑zxn_f(𝐫;t)`$ and similarly for $`<y(t)>`$ and $`<z(t)>`$. Of course, during free expansion the center-of-mass coordinates must remain unchanged: this property is used as a test of the numerical method.

Fig. 2 shows that the calculated $`\sigma _r_{}`$ (circles) and $`\sigma _z`$ (squares) agree with the results from the analytical expressions (33) and (34) (solid lines), where the classical expression has been used for $`E_{rel}`$. Snapshots of the density profiles at selected times are shown as contour plots in fig. 3. The definition of the profile degrades in time because the number of particles per cell drops during the expansion.

After this test of the reliability of the simulational method, we proceed to use it in some original applications.

### B Expansion of a mixture of a condensate and a cold-fermion cloud

As a first novel application we look at the case in which a core of Bose-condensed atoms is present inside the dilute Fermi gas. We prepare a state with $`N_f=1000`$ $`{}_{}{}^{40}K`$ atoms and $`N_c=10^5`$ $`{}_{}{}^{39}K`$ atoms in identical harmonic traps and at the same temperature as for the Fermi gas studied in Sec. IV.A. The scattering length which describes the interactions between the atoms in the condensate is $`a_{bb}=80a_{Bohr}`$, while the interspecies scattering length is $`a_{bf}=40a_{Bohr}`$.

The inclusion of the GPE algorithm at fixed mesh size normally requires shorter time-steps to mantain stability. We make the choice of a thinner $`501`$x$`1001`$ mesh than in Sec. IV.A in order to keep the time-step at $`\omega \mathrm{\Delta }t=10^4`$, all other simulation parameters remaining the same.

The initial state is characterized by $`\mu _b=0.52E_F=9.51\mathrm{}\omega _f`$ and $`\mu _f=0.10E_F=1.83\mathrm{}\omega _f`$. We display in fig. 4 the average widths $`\sigma _r_{}`$ (circles) and $`\sigma _z`$ (squares) for both the fermionic cloud (open symbols) and the condensate (filled symbols). The solid lines are the analytical solution for the ideal Fermi cloud as in fig. 2. It is seen that with the scattering lengths of the <sup>39</sup>K-<sup>40</sup>K mixture the mean-field force of the inner condensate core on the outer fermionic cloud is not strong enough to sizeably affect the expansion of the latter.

Snapshots of the condensate density profiles and of the fermionic cloud at times $`t=0`$, $`8.5`$, $`17`$ and $`25.5`$ ms are displayed as contour plots in figs. 5 and 6. Comparison with those in fig. 3 shows that the reduced number of computational particles per cell tends to increase the statistical noise. This degradation worsens as the simulational time elapses, as is evidenced by the last snapshot in fig. 6.

### C Oscillations of Bose gases inside an optical lattice

Here and in the following subsection we apply our numerical method to study the dynamics of a Bose-Einstein condensate and a thermal cloud of <sup>87</sup>Rb atoms at finite temperature inside a one-dimensional optical lattice. The initial state is prepared by adding to the harmonic trap, described by $`V_{trap}(𝐫)=(1/2)m\omega ^2(r_{}^2+ϵ^2z^2)`$, a periodic potential given by $`V_{latt}(z)=\alpha E_R\mathrm{sin}^2(k_Lz)`$, where $`E_R=\mathrm{}^2k_L^2/(2m)`$ is the recoil energy and $`k_L=2\pi /\lambda `$ is the wave number of the laser beam which creates an optical lattice with period $`\pi /k_L`$ in the axial direction.

Such a system, which has been realized at LENS  and examined numerically at $`T=0`$ by two of us , shows a rich variety of dynamical behaviors. Thus, the study of the sloshing-mode oscillations of an almost pure condensate with $`N=310^5`$ atoms in a lattice with $`\alpha =1.6`$ shows that superfluidity is superseded by dissipation as the initial displacement of the condensate away from the harmonic-trap center is increased. This behavior is quantitatively understood as a gradual destruction of superfluidity via emission of sound waves in the periodically modulated inhomogeneous medium. Below the dissipative threshold, on the other hand, the oscillatory motion of the condensate through the optical lattice can be mapped into the dynamics of superconducting carriers through a weak-link Josephson junction . This implies the possibility of observable resonances and of multimode behavior.

Here we extend the above numerical studies to contrast the oscillations of a condensate with the motions of a thermal cloud. We prepare initial states for the two cases $`\alpha =1`$ and $`\alpha =5`$, both at $`T=0`$ for the BEC  and for the thermal cloud at temperature $`T`$ above the critical temperature $`T_c`$. We give an initial displacement $`\mathrm{\Delta }z=42.6\mu `$m to the trap center and follow the subsequent dynamics with a time-step $`\omega \mathrm{\Delta }t`$ of order $`10^5`$.

The snapshots of the atomic density show that for $`\alpha =1`$ (see Figs. 7 and 8) the condensate behaves as a superfluid executing harmonic oscillations at a frequency equal to the trap frequency, while the thermal cloud at $`T>T_c`$ diffuses away in a quarter of a period. For $`\alpha =5`$ (see Figs. 9 and 10) the condensate breaks instead into fragments as it attempts to perform the first oscillation, and after a period its center of mass becomes localized at the bottom of the harmonic well. In the same setup the thermal cloud becomes localized at the center of the trap in one tenth of a period and spreads out.

Figure 11 gives a clear picture of these behaviors by reporting the axial center position and width of the condensate and of the thermal cloud as functions of time in the two cases.

### D Expansion of a Bose-condensed gas in an optical lattice

In our last study we look at the expansion of a Bose-Einstein condensate and its thermal cloud, which initially reside in a harmonic well and a superposed optical-lattice potential. The external potentials are characterized by parameters typical of an experiment at LENS , namely $`\omega =2\pi 90`$ rad/s, $`ϵ=8.9/90`$, $`2\pi /k_L=795`$ nm and $`\alpha =5`$. The condensate contains $`6935`$ <sup>87</sup>Rb atoms and the thermal cloud is composed of $`3065`$ atoms: the temperature of the gas is $`T=86nK=0.24E_R/k_B`$ and its chemical potential is $`\mu =5.86\mathrm{}\omega =0.14E_R`$. We use a mesh of 111$`\times `$2801 points with 308000 representative particles.

We evolve the gas with a time-step $`\omega \mathrm{\Delta }t=710^5`$ after switching off both the harmonic trap and the periodic potential. Snapshots of the atomic densities of the condensate and of its thermal cloud, taken at the moment in which the potentials are switched off and after 3.5, 7 and 10.5 ms of free expansion, are shown in figs. 12 and 13. The condensate is seen in fig. 12 to develop side bands which separate out of the central cloud, while the thermal cloud in Fig. 13 simply spreads out. These features of our numerical results reproduce those observed in the experiments .

The appearence of side bands in the condensate during expansion is due to Bragg scattering against the periodic potential. In fact, in a long-time simulation run of a one-dimensional model of the expansion we have found that the condensate side-bands move at velocity $`v2\mathrm{}k_L/m`$, corresponding to the momentum associated with the first reciprocal vector of the optical lattice.

## V Computational remarks

We have assessed the computational performance of the numerical method by repeating the test of Sec. IV.A after changing either the number of computational particles or the mesh size. We list in Table 1 the computational times elapsed while running the HPF-PGI-compiled code on a fully dedicated 1GHz Pentium III SCSI.

These data provide the following values for the specific CPU time costs of the GPE and VLE of the code per time-step:

$$\{\begin{array}{cc}t_{GP}1.5\hfill & \mu s/gridpoint\hfill \\ t_{VL}0.6\hfill & \mu s/particle\hfill \end{array}$$
(37)
These figures invite a number of comments. First, they show that the VLE section can evolve just a few computational particles while a single grid-point of the GPE solver is advanced. Since statistical accuracy requires of the order of $`1050`$ particles per cell, we conclude that the VLE solver is a potential computational bottleneck.

Let us nonetheless assume that the VLE and GPE sections can evolve on a one-to-one CPU time basis. We can then focus on the grid part only and estimate the feasibility of large-scale applications to finite-temperature condensates in optical lattices. Covering a simulation span of $`100`$ $`ms`$ in steps of $`0.1`$ $`\mu s`$ requires $`10^7`$ time-steps. At a cost of $`1`$ $`\mu s`$ per time-step and grid-point, a grid with, say, $`10^6`$ grid-points takes of the order of $`10^6`$ seconds, namely almost two weeks of CPU time to complete.

Ways to achieve substantial speed-up are clearly needed. Among others, two promising (and not mutually exclusive) strategies are non-uniform meshes and parallel computing. Both strategies appear conceptually straightforward and will be the object of future work.

## VI Concluding remarks

The increasing complexity and variety of phenomena observed in current studies of the dynamical behavior of normal and superfluid quantum gases at finite temperature motivate the development of suitable numerical tools to assist theoretical understanding.

To this aim, we have combined a particle-in-cell method with an explicit time-marching algorithm to evaluate the time evolution in models of a Bose-Einstein condensate and a cold-atom cloud.

We have tested the method against known analytical results in the simple physical situations offered by the expansion of a collisionless fermionic cloud without and with an inner Bose-condensed core. We have also applied it to simulate novel experimental observations on the dynamical behavior of a condensate with its thermal cloud in a harmonic plus optical lattice potential, where we have found substantial accord with current experiments.

We have also analyzed those computational aspects of the algorithm which are most relevant to applications in large-scale problems. This analysis emphasizes the need for non-uniform meshes and parallel computing. On the physics front, an extension of the method to include the quantum collisional integrals is under way.

###### Acknowledgements.

This work was partially supported by INFM under PRA2001. One of us (MLC) thanks Dr. F. S. Cataliotti and Dr. S. Burger for making their experimental results available prior to publication.