Quantum transport of non-interacting Fermi gas in an optical lattice combined with harmonic trapping

We have considered non-interacting Fermi gas in a combined harmonic and periodic potential. We calculate the energy spectrum and simulate the motion of the gas after sudden replacement of the trap centre. For different parameter regimes, the system presents dipole oscillations, damped oscillations around the replaced centre as well as localization. The behaviour can be explained by a change in the energy spectrum from linear to quadratic.

Experimentally, optical lattices are highly controllable and open a wide range of parameters for studying the dynamics. For bosonic atoms, quantum transport, such as Bloch oscillations, has been observed [2]. Such experiments on fermions in optical lattices have been realized very recently [3]- [6].
Unlike electrons in bulk solids, atomic gases feel, in addition to the lattice, an overall confining potential, usually harmonic. The harmonic trapping can be chosen to be weak and the system can be considered to be nearly homogeneous. On the other hand, it can be made significant to induce interesting quantum transport experiments: a shift in the trap centre may cause oscillatory motion of the gas [3]- [5]. Single-particle states and the spectrum in a combined harmonic and periodic potential have been studied in [7,8], and the effect of the potential on quantum transport has been analysed in [5,9] using semiclassical analysis. We present in this paper the numerical treatment, which gives a complementary point of view for understanding the problem and provides a method applicable for a wide range of parameter values. Since we consider non-interacting fermions, the only physical assumptions we need are that the individual atoms obey the Schrödinger equation and that the equilibrium of the many-particle system can be described by the grand canonical ensemble. The physical model is considered in more detail in section 2. Details of the numerical implementation are discussed in section 3. Finally, in section 4, results of the simulations are presented and discussed and are compared with the experiments.

Physical model
Let us consider a one-dimensional optical lattice combined with a three-dimensional quadratic potential. The quadratic potential is weak in the direction of the lattice (axial direction) and tight in the remaining two (radial) directions. The full three-dimensional Hamiltonian operator for a single particle is where ω a and ω r are the axial and radial frequencies, respectively, U is the lattice height and λ is the lattice wavelength. We will consider first the one-dimensional system in the axial direction and later deal with the radial degrees of freedom. Choosing the units of energy and length to be hω a and √ h/mω a , respectively, and using the notation A = U/2 and k = 4π/λ, the one-dimensional Hamiltonian can be written as a sum of the lattice potential H L = A cos(kx) and the oscillator Hamiltonian Let φ n and n be the eigenfunctions and eigenenergies of H, respectively. Using the annihilation and creation operators a n , a † n associated with the eigenfunctions, we may write the many-particle Hamiltonian in the fermionic Fock space aŝ H = n n a † n a n .
Initially, the system is in statistical equilibrium and can be described by the grand canonical ensemble where µ is the chemical potential and N = n a † n a n is the number operator, and β = 1/k B T . The time evolution of the system is given by where H d is the Hamiltonian governing the dynamics of the system. Obviously, H d must be different from the Hamiltonian H used to calculate the initial equilibrium state (0), otherwise the system will be stable. In the experiment, the displacement d of the harmonic trap gives rise to a new Hamiltonian which starts the time evolution of the system. Displacement of the trap centre in the absence of the lattice leads to dipole oscillations that are characteristic of a harmonic potential. We aimed to calculate the time evolution of the expectation value of the position operator. The task can be greatly simplified by the fact that we can represent any many-particle state by a single-particle state ρ( ) when calculating expectation values of single-particle observables. This reduction is discussed in more detail in the appendix. The single-particle state ρ( (0)) corresponding to the grand canonical ensemble (0) is a diagonal operator in the energy eigenbasis. The diagonal elements are given by the Fermi-Dirac distribution When the Hamiltonian operator is itself a single-particle operator, as for a non-interacting Fermi gas, the reduced state obeys the ordinary Schrödinger equation. Hence, the final formula for the motion of the Fermi gas is simply where X is the position operator.
As the final step, we consider dimensional reduction from three dimensions to one. Even though we are interested in the motion of the Fermi gas only in one dimension, the physical system is naturally three-dimensional. However, we can simply trace out the radial degrees of freedom since they are independent of the observable that we are measuring. In practice, this is done by replacing the Fermi-Dirac distribution by the effective one-dimensional energy distribution where the multiples of hω r exhaust the energy spectra of the independent oscillators in the radial direction.

Numerical implementation
The time evolution of the centre of mass of the non-interacting Fermi gas trapped in an optical lattice was calculated using formula (2). The first task was to diagonalize the Hamiltonian. Actually, there were two Hamiltonians, H and H d , to be considered, corresponding to the two positions of the trap. However, we made the simplifying assumption that the displacement of the trap was a multiple of the lattice wavelength (this is justified when the latter is much smaller than the former). Then, by a symmetry argument, the two sets of eigenfunctions were obtained from each other by a simple translation. According to the equation (1), the initial equilibrium state (0) is diagonal in the eigenstate basis of the initial Hamiltonian H, and the diagonal elements are determined by the corresponding eigenvalues. The chemical potential µ was fitted to match the particle number. By a change of basis, the state matrix was then transformed to the eigenstate basis of the displaced Hamiltonian H d . The position operator X was also expressed in this basis. By definition, the Hamiltonian H d itself is diagonal in this basis, so that the time evolution could be easily calculated using equation (2).
The main difficulty was to find a basis suitable for diagonalizing the Hamiltonian numerically. We used the Hermite functions where H n is the nth Hermite polynomial. With this choice, the harmonic oscillator Hamiltonian H O becomes diagonal with matrix elements and even the matrix elements of the lattice part H L can be calculated analytically. To be exact, the matrix elements are recovered as H L nm = A Re(I nm ), where We were able to derive explicit closed expressions for I nm (k) as finite polynomials of k. However, they are rather cumbersome and not very useful in practice. In numerical calculations, the integrals can be most easily obtained from the formula where h nm satisfies the recursion equation with the initial conditions h 00 = 1 and h nm = 0 if either of the indices is negative. We proved equations (3) and (4) using generating function methods. The recursion formula is efficient, but numerically unstable. Therefore high-precision arithmetics had to be used. The Fourier transform F is a unitary transform on the space of square-integrable functions and Hermite functions are eigenfunctions of F, Fψ n = i n ψ n ; therefore, we can use equations ( 3) and (4) to calculate the convolutions Apart from being interesting in their own right, these convolutions yield the matrix elements of the change of basis from (ψ n (x)) to the translated basis (ψ n (x − d )), which is exactly what is needed to transform the state matrix from the eigenstate basis of H to that of H d . Hence, equations (3) and (4) were used twice in the implementation.
For some values of the parameters, it was beneficial to use the scaled Hermite functions ( √ λψ n (λx)) with λ > 1 in the diagonalization for improving the resolution. Even in this case, the matrix elements needed are easily derived from the above formulae. In the simulations, we used less than 3000 basis vectors and calculated typically a few hundred eigenstates. This was probably not enough for all the cases, but at this point our computer resources became a limiting factor.

Results
We used the above methods to simulate the motion of the free Fermi gas of potassium atoms in an optical lattice confined in a magnetic potential and we compared our results with the experimental results described in [3]- [5]. The average temperature of the gas in the experiment described in [5] was T = 90 nK, the lattice wavelength was λ = 863 nm, the axial and radial frequencies were ω a = 2π × 24 s −1 and ω r = 2π × 275 s −1 , respectively, and the average atom number was 25 000. Recall that k = 4π/λ. The depth U = 2A = sE R of the lattice, expressed in terms of the recoil energy E R = h 2 /2mλ 2 , was varied. Initially, the trap was replaced by 15 µm to induce oscillations. Results of the simulation are shown in figure 1. The energies of the lowest-lying single-particle eigenstates are also plotted against the ordinal number in figure 2 to get some qualitative insight of the system. In the absence of the lattice potential (s = 0), the system is reduced to a harmonic oscillator and the average position oscillates freely, showing the characteristic dipole oscillation in a harmonic potential. None of the eigenstates are localized and the corresponding energy spectrum is linear. The energy quantum or the difference between the consecutive eigenvalues gives the frequency of the oscillations.
When the height of the lattice is increased to s = 3, both localized and non-localized eigenstates play a role. The oscillations of the system are damped, which can be understood as follows. According to equation (2), the centre of mass of the cloud is a sum of periodic terms whose frequencies are given by the differences, (ε n − ε m )/h , between the eigenvalues ε n of the dynamical Hamiltonian. When the eigenvalues are not evenly spaced, destructive interference causes dephasing.
Increasing the lattice height also increases the offset of the dipolar oscillations since more and more particles become localized. When the height of the lattice is increased to s = 8, localized eigenstates dominate. The system is almost frozen to its initial position. The energy spectrum becomes almost quadratic. The spectrum corresponds to states that are localized in individual lattice potential sites. Neighbouring sites are shifted in energy according to the superimposed harmonic potential; therefore the quadratic form of the spectrum (the difference in energy of the two localized particles depends essentially on their locations in the lattice). Actually, each energy is doubly degenerate (the symmetric and the antisymmetric states formed from lattice sites on opposite sides of the centre), which can be resolved by a closer look at the spectrum. Therefore, to be accurate, localized particles correspond to a superposition of these doubly degenerate states. Representative examples of localized and non-localized eigenstates are shown in figure 3. It is noteworthy that in the intermediary region, the lowest-lying states are not localized. This is indicated by the linear parts of the spectra near the origin. Hence, the atoms do not tend to localize when the system is cooled down to ultra-low temperatures. Note that higher energy states also may not be localized, i.e. have a linear spectrum, as is shown for one set of parameter values in figure 2. This is associated with the existence of higher bands of the optical lattice potential. Such states may affect significantly the localization and oscillation behaviours of the gas at higher temperatures.
The results are in very good qualitative and quantitative agreement both with the experiment and the semiclassical approximation presented in [5]. We also compared our simulations with two earlier experiments by Modugno et al [3,4]. In these cases the experimentalists observed oscillations with a considerably higher amplitude and a higher frequency compared with what was predicted by the simulations. We believe that this was simply because the physical parameters used in these experiments were much more demanding for our numerical method and would have required more eigenstates to be calculated than what is currently possible. The higher the displacement, number of particles and temperature, the higher the number of eigenstates needed for reliable results. The computer resources at our disposal were rather modest and hence we have hardly pushed the method to its limits.

Conclusions
We have analysed the energy spectrum and motion of a Fermi gas in a combined potential. The harmonic potential alone has a linear spectrum and shows dipole oscillations of the gas when the trap centre is displaced. Introducing the periodic potential transforms the spectrum gradually into a quadratic one and the gas becomes localized. Intermediate cases show damped oscillations around the displaced centre position. Our approach is complementary to the semiclassical analysis presented in [5,9], and is in good qualitative agreement with the experiments. For the parameter values that are not too challenging from the numerical point of view, the results agree also quantitatively both with the experiments and the semiclassical analysis.
It is important to note that the quantum transport in the combined harmonic and periodic potential is qualitatively quite different from the Bloch oscillation or Wannier-Stark ladder behaviour in tilted periodic potentials. The two cases approach each other only when the gas moves in a region where the harmonic trap can be approximated by a linear potential. This requires a much shallower trap than that used in the experiments of [3]- [5]. Understanding the behaviour of the non-interacting Fermi gas sets the basis for investigating the effects of interactions [4,10] and, eventually, superfluidity [11,12].
as the position x. In mathematical terms, the phase space is the cotangent space T * C of the configuration space and it comes equipped with a canonical non-degenerate two-form which is known as the symplectic form, whereas the pair (M, ω) is referred to as a symplectic manifold. The symmetries of the phase space are the diffeomorphisms of M preserving ω. Accordingly, infinitesimal symmetries are the vector fields A on M whose Lie derivatives satisfy L A ω = 0. Then, each infinitesimal symmetry generates a one-parameter group of symmetries at least locally, where the correspondence is given by the flow of the vector field.
Under some mild conditions, we can associate a classical observable or a function H A on the phase space M with an infinitesimal symmetry A. The observable is defined by the equation which is required to be valid for all vector fields B on M. In particular, since the time evolution of a classical mechanical system preserves the symplectic form, it is generated by an observable, known as the Hamiltonian. In the example above, it is clear that both translations and rotations are symmetries of the phase space. An infinitesimal translation is just a three-vector t ∈ R 3 , and the corresponding observable is Similarly, an infinitesimal rotation can be given by a three-vector r as r · R, where the components of R are the infinitesimal rotations about the coordinate axes. The associated classical observable is then given by the relation The general pattern emerging is that, given a Lie algebra g of infinitesimal symmetries, the collection of associated observables is most conveniently expressed as a map ρ from the phase space M to the dual g * , i.e. the space of linear maps g → R. In this way, we associated the momentum p with the translations and associated the impulse moment x × p with the rotations. Therefore the map ρ : M → g * , defined by the equation for all A ∈ g and x ∈ M, is known as the moment map. The Lie group G corresponding to the Lie algebra g acts naturally on both g and g * , and the moment map is equivariant, in the sense that ρ(gx) = gρ(x) for all g in G.
Returning to quantum mechanics, let H be a Hilbert space, on which a compact group G acts unitarily. This is the general framework for a quantum mechanical system with the symmetry group G. Now, the projective space of pure states is canonically a symplectic manifold on which G acts preserving the symplectic form. Hence, from the point of view of symplectic geometry, all that we have said previously applies.
Let g be the Lie algebra of G. The observable associated with an infinitesimal symmetry A ∈ g is calculated to be the expectation value H A (v) = v|Â|v , whereÂ stands for the representation of A in H. There is also the moment map ρ : P(H) → g * from the projective space to the dual of g defined by the equation Hence, the expectation value ofÂ in the vector state v is given by Â = (ρ(v), A). We can extend the moment map to any state by linearity. Then, the moment map image ρ( ) is a kind of reduced state that can be used to calculate expectation values for all observables in g. The point is that ρ( ) may be a much simpler object than the state we started with. Furthermore, the equivariance of the moment map means that ρ(Û Û −1 ) = Uρ( )U −1 for every element U of G. In particular, if the Hamiltonian H itself belongs to g, then the time evolution of the reduced state ρ( ) is given by the usual Schrödinger equation.
In the application we have in mind, G is the group of unitary transformations of the singleparticle Hilbert space V, and g may be identified with the single-particle observables or the Hermitian operators on V. Using the pairing (B, A) = Tr(BA), we may identify g with its dual. Let H be the fermionic Fock space associated with V. Then, G acts naturally on H and by the general theory Â = Tr( Â ) = Tr(ρ( )A) (A.1) for any many-particle state and single-particle observable A. It is easy to check that ρ( ) is a positive operator with trace Tr ρ( ) = Tr(ρ( )I ) = Tr( Î ) = Tr( N ) = N , since the number operator N is indeed the representation of I ∈ g in the fermionic Fock space.
In conclusion, we can represent any many-particle state by a single-particle state ρ( ) when calculating the expectation values of single-particle observables. It is easy to calculate ρ( ) explicitly in a given basis (v n ) of V using equation (A.1). In terms of the annihilation and creation operators associated with this basis, the matrix elements of ρ( ) are ρ( ) nm = Tr( a † m a n ) = a † m a n . Applying the above procedure to the grand canonical ensemble (1) yields where f(x) = 1 1 + e βx is the Fermi-Dirac distribution, which is usually derived by the partition function approach. The above discussion is not just an appealing generalization but is given to show that the reduction to the single-particle state formalism gives not only the energy distribution, but also a state matrix obeying the Schrödinger equation.