Interaction-induced conducting-nonconducting transition of ultra-cold atoms in 1D optical lattices

The study of time-dependent, many-body transport phenomena is increasingly within reach of ultra-cold atom experiments. We show that the introduction of spatially inhomogeneous interactions, e.g., generated by optically-controlled collisions, induce negative differential conductance in the transport of atoms in 1D optical lattices. Specifically, we simulate the dynamics of interacting fermionic atoms via a micro-canonical transport formalism within both mean-field and a higher-order approximation, as well as with time-dependent DMRG. For weakly repulsive interactions, a quasi steady-state atomic current develops that is similar to the situation occurring for electronic systems subject to an external voltage bias. At the mean-field level, we find that this atomic current is robust against the details of how the interaction is switched on. Further, a conducting-to-nonconducting transition exists when the interaction imbalance exceeds some threshold from both our approximate and time-dependent DMRG simulations. This transition is preceded by the atomic equivalent of negative differential conductivity observed in transport across solid-state structures.

U=0 t=0 U=0 t>0 I U>0 Figure 1. Schematic plot of the experimental set up we simulate: Non-interacting ultra-cold fermions are loaded into the lowest energy state of a 1D homogeneous optical lattice. At t = 0, a focused laser beam then induces on-site repulsive interactions of the Hubbard type on the left half of the lattice. Due to the imbalance of interactions, a mass current is induced. We emphasize that there is no need to "mechanically" set the atoms out of equilibrium by tilting a potential nor is it necessary to introduce additional species of atoms or dissipation. Moreover, this setup may also be implemented to study quench dynamics [17]. The gray dots emphasize that atoms may be in a superposition of different quantum states.

I. INTRODUCTION
Advances in experimental studies of quantum transport of ultra-cold atoms in optical lattices [1][2][3][4][5] draw attention to different aspects of systems out of equilibrium. In conventional condensed matter settings, an electrical current is induced by applying a voltage difference across the sample. The interaction of the charge carriers are usually homogeneous, except near impurities. In contrast, ultra-cold atoms are charge neutral and so far the motion of an atomic cloud has been generated by placing the cloud away from its equilibrium position [1] or by a sudden shift of the minimum of its trapping potential or its distortion [2][3][4][5]. The interactions are generated by either adding another species of atoms [1] or tuning magnetic-field controlled collisions [3,4].
In this paper we instead suggest a method to drive a mass current of atoms via local tuning of interactions. This, in turn, reveals interesting phenomena otherwise difficult to observe in conventional solid state systems. A key to realize interaction-induced transport is controllable inhomogeneous interactions -a novel possibility offered by optically tunable collisions of ultra-cold atoms. The optical Feshbach resonance (OFR) [6][7][8][9] is a promising technique for controlling interactions locally using focused laser beams. For instance, the spatial modulation of density using the OFR of bosons [7] demonstrates the power of controllable inhomogeneous interactions. Other exotic equilibrium structures may also be generated by using this tool [10].
We also point out that these optically controlled interactions can be time dependent. Suddenly turning on the interactions, for instance, will drive the system out of equilibrium. Thus, instead of using external driving mechanisms such as a voltage bias, it is possible, with suitable initial conditions and patterns of interactions, to induce a mass current. We show that for reasonably weak interactions the mass current of fermions is similar to a charge current in electronic systems. Further, increasing the interaction strength leads to the atomic equivalent of negative differential conductance, where the current reaches a maximum value and then subsequently decreases. A mean-field approximation predicts that this decrease leads to a conducting-to-nonconducting transition after the interaction strength exceeds some threshold (dependent on the filling). This transition may be explained by a mismatch of energy spectra between the interacting and non-interacting parts of the system, but this argument does not rule out other possible current-carrying states.
Time-dependent density-matrix renormalization group (td-DMRG) simulations [11], as well as a higher-order approximation, demonstrate that indeed a conducting-nonconducting transition should exist in the strongly interaction-imbalanced regime when the initial state is not far away from a band insulator. The threshold value for the transition, however, differs from different approaches. Importantly, the many-body negative differential conductance remains at all levels of approximation. This is the counterpart to negative differential conductivity observed in solid-state structures [12,13], where changes in the carrier density or subdivisions of the Brillouin zones cause non-monotonic dependence of the current on the external field strength. We note that, recently, negative mobility of cold atoms in temporally modulated optical lattices [14] and negative differential conductance in fermion-boson mixture in optical lattices [15] have been discussed. We propose, alternatively, that one may observe negative differential conductance in a setup where the optical lattice potential is static and by tuning the interactions in real space. We focus on fermions but notice that, for the Bose-Hubbard model, Ref. [16] considers the sudden connection of a superfluid and a Mott insulator with different chemical potentials and found a mass current as well.

II. INTERACTION-INDUCED TRANSPORT
In our approach -which explicitly follows the dynamics of mass transport -we implement the micro-canonical formalism (MCF) [18,19] as in Refs. [20,21] which monitors the evolution of the correlation matrix. When applied to a system of noninteracting atoms on a 1D optical lattice suddenly losing the atoms on the right half of the lattice, MCF shows very different dynamics for bosons and fermions [20]. While the bosonic current decays to zero in the thermodynamic limit, the fermionic current exhibits a quasi steady-state current corresponding to a plateau in the current as a function of time. Importantly, the MCF is designed for finite closed quantum systems so it is particularly suitable for studying the dynamics of ultra-cold atoms. Although our proposed setup can be applied to bosons, to explore their dynamics theoretically one has to consider also the effect of the quasi-condensate which renders the analysis more involved. Therefore, in this paper we focus only on fermions.
We consider N pσ (σ =↑, ↓) two-component fermions, with N p↑ = N p↓ , loaded into a 1D optical lattice of size N , see Fig. 1. This type of finite-length 1D lattice may be generated by inserting a thin optical barrier [22] into a ring of optical lattice [23]. The resulting C-shaped lattice is geometrically identical to the setup here. This confining potential creates a uniform lattice and does not require a background harmonic trapping potential parallel to the lattice. The background interactions are assumed to be negligible. The system is initially described by the tight-binding tunneling Hamiltonian H 0 = −t ij ,σ c † iσ c jσ , where ij , t, c † iσ (c iσ ) denote nearest-neighbor pairs, the hopping coefficient, and the creation (annihilation) operator of site i, respectively. The unit of time is t 0 ≡ /t and is about few ms for reasonable lattice depths [20].
The initial state corresponds to the lowest energy state of H 0 , with no correlations between up and down spin states (see, e.g., Refs. [20,21]). The system is then set out of equilibrium by introducing interactions among atoms on the left half of the lattice by, e.g., optically-induced collisions using a focused laser beam. Here we concentrate on moderate repulsive interactions and non-integer filling. This allows us to focus on the dynamics induced by the interactions and avoid unnecessary confusions about possible equilibrium phase transitions such as the Bardeen-Cooper-Schrieffer instability or the Mott-insulating phase in uniform and static interacting systems [24]. We model an instantaneous switch-on of the interactions with a sharp interface between the interacting and non-interacting regions and relax the former condition later on.
The Hamiltonian generating the dynamics is Heren iσ = c † iσ c iσ , U is the onsite repulsive coupling constant, L denotes the left half of the lattice, andσ is the opposite of σ. This model should be appropriate for moderate lattice depth [25]. Sincet can be tuned by the lattice depth and U can be tuned by OFR, U/t can span a broad range so our study should be relevant to experiments.
The equations of motion for the correlation matrix are where [·, ·] denotes the commutator of the corresponding operators. The explicit expression can be found using standard anti-commutation relations. After some algebra, one gets Here iσ c j−1,σ and n iσ = n iσ . We solve Eq. (2) both at the mean-field level and by adding higher-order correlations.
The mean-field level solution can be obtained by Wick decomposition of c † iσ c iσ c † iσ c jσ as c † iσ c iσ c † iσ c jσ since no spin-flip mechanisms are present. This is the standard Hartree-Fock approximation. The equations are closed after this approximation. Thus the dynamics after the interaction is switched on can be monitored by integrating Eq. (2) with the initial condition c ijσ (t = 0) given by the lowest energy state of the non-interacting Hamiltonian. We update Eq. (2) in a symmetric fashion so that c † iσ c jσ = c † iσ c jσ at any time. Later on we will also go beyond mean-field by developing a higher-order approximation and by performing td-DMRG simulations.
In addition to the current I σ = 2tIm c † N/2,σ c N/2+1,σ , we also evaluate the particle number on the left half lattice, N Lσ = N Lσ and its fluctuations ∆N 2 Lσ = N 2 Lσ − N Lσ 2 , which is another property that can be determined experimentally. Herê N Lσ = i∈Ln iσ . Explicitly, Right before the interactions are switched on, there can be initial number fluctuations due to the wave nature of the initial quantum state. We thus present ∆N 2 L(σ) − ∆N 2 L(σ) (t = 0) which reflects the change of the number fluctuations due to the interaction-induced dynamics. Figure 2 shows the particle number and its fluctuations of the left half lattice with N = 128 and N pσ = 64 for each species.

III. RESULTS AND DISCUSSIONS
Since the system size is finite, there is reflection of the current off the boundary after a revival time but we focus on dynamics   before this revival occurs. As U/t increases, one can see that instead of inducing more fermions to move from the left to the right, above a certain threshold value the transport is no longer observed within a reasonable time scale (e.g., t ≤ (N/2)t 0 ). This corresponds to a conducting-nonconducting transition and it can be found at other ratios of f ≡ N pσ /N . In general the threshold U c /t decreases as f increases. We will show that this is because for higher filling the onsite interaction energy dominates more easily. One can see this transition even better in the average current itself. We note that it has been shown that with suitable initial conditions non-interacting fermions can develop a quasi steady-state  current (QSSC) [20,26]. This is so also for the present interacting case. The fermionic current is shown in Figure 3(a) for N = 128 and N pσ = 64 for each species. The current may be measured by the protocol of Refs. [20,27]. The plateaus shown in I σ (t) indicate the existence of QSSCs. To smooth over small oscillations, we define an averaged current as I σ ≡ (1/30t 0 ) 40t0 10t0 dtI σ (t) and plot I σ as a function of U/t for N pσ = 32, 64, 96 in Fig. 3(b). For small U/t we clearly see that the averaged current increases linearly with U/t. Then the dependence deviates from a linear form before the conductingnonconducting transition. For even larger U/t, I σ decreases until no finite averaged current can be found. The non-monotonic dependence of I σ on U resembles the negative differential conductivity observed in conventional solid-state devices [12,13]. In the strongly interaction-imbalanced regime, we found oscillations in the density profile and further studies may connect those oscillations to the charge domains in solid-state devices exhibiting negative differential conductivity.
The interaction term thus acts as a bias on the left half lattice if U/t is small. Here we investigate the validity of this statement. By extracting the slope of the linear part of I σ for small U/t, one can study how efficient the interaction can induce a current. The inset of Fig. 3(b) shows the effective conductance G ef f ≡ f −1 d I σ /dU as a function of N pσ for a fixed N , where f = N pσ /N is the initial filling factor. As N pσ → N the system approaches a band insulator [28] and there is only a tiny range of U/t where a finite I σ exists. This could contribute to some inaccuracy of the extraction of the slope close to N pσ → N . Nevertheless, as shown on the inset of Fig. 3(b) the effective conductance is close to the ideal quantized conductance G 0 = (2π ) −1 for arbitrary initial filling. This establishes the feasibility of using a weak inhomogeneous interaction as an efficient driving force for inducing a current and its correspondence with a bias, ∆µ = f U .
In the weak interaction regime, we found another interesting phenomenon. By considering different ways of switching on the interaction, one may wonder whether the system will reach the same magnitude of the QSSC. We simulate different switchingon scenarios by introducing a time scale t m such that the interaction grows from zero to its full value during t m and remains a constant after that. Figure 4 shows that despite different time-dependences of the interaction, the system always reaches the same magnitude of the QSSC. We have also checked other more complicated functional forms for switching on the interaction and found the same QSSC. Importantly, even if the system is over-excited by a spike in the time dependence of U , when U comes back to a constant value the height of the QSSC coincides with the case of a sudden switch-on of the interaction as shown on Fig. 4. Thus the magnitude of the QSSC is robust (memory-less) against different ways the interaction is switched on. This feature makes the interaction-induced transport appealing for making devices since a steady output (the QSSC) is insensitive to the transient behavior of the driving force (the interaction).
When the optically-induced interactions are switched on, it is possible that the onsite potential could be shifted due to the increase of the kinetic energy of the atoms interacting with the incoming photons. We study possible effects by adding an extra term V i∈L,σn iσ to H e of Eq. (1) to simulate this effect with a positive V . This shift of the onsite potential acts like a bias so one might expect that more current will flow to the other side. However, the negative differential conductance and the mean-field conducting-nonconducting transition are robust against this positive potential shift. The energetic mechanism described below will make it clear why a positive potential shift will not affect the results qualitatively, thus showing that it is not sensitive to this type of modification of the Hamiltonian.
We now discuss the reason behind the negative differential conductance and the conducting-nonconducting transition, which is energetic in nature. In order to obtain a "macroscopic" steady-state current, αN particles must be transferred from the left to the right halves, where α is some small proportionality constant. However, this entails a change in energy, as the particle density decreases on the site with interactions, and this effect must be compensated by an energy change due to rearrangement of the single particle states. To be more concrete, immediately after interactions are turned on, the system is still in the ground state of H 0 . However, the energy has shifted to where Wick's theorem has been applied since the system is in the ground state of H 0 . The subscript of · · · denotes the time when the expectation is taken. For half-filling, this gives an energy E = E 0 + U N/8. For arbitrary filling f , it is where the expression is approximate since not all site occupation numbers will be exactly f (for each spin component) . The energy at some later time will still be E (since no further changes have been made to the Hamiltonian), but will have the form in the mean-field approximation. Assuming the αN atoms are taken uniformly from the left half, the energy is Taking the difference of Eq. (7) and Eq. (5), and using energy conservation, gives H 0 t − E 0 ≈ αf U N . However, the quantity H 0 t −E 0 can not be arbitrary, since we have a finite bandwidth. Regardless of whether one has an interacting or noninteracting state, the largest energy change due to the H 0 contribution is 4tαN , corresponding to taking the αN particles from the lowest part of the band (on the left half) to the highest part (on the right half). This gives and thus results in a threshold value of U . Beyond this value of U , a macroscopic current can not flow within the mean-field approximation and when the particles are taken uniformly from the interacting lattice. We remark that here we consider the lowest energy band of the optical lattice. If the interaction energy exceeds the band gap separating different bands, the system may re-enter a conducting state. We also remark that energy mismatches can results in other interesting phenomena as discussed in Refs. [29,30].
To demonstrate this energetic mechanism, Fig. 3(c) and (d) show the energy spectra of the left (L) and right (R) half lattices when U is switched on for N = 128 with N pσ = 64. When U/t = 7 (panel (c)), there is still an overlap between the two spectra and exchanging atoms in an energy-conserving fashion is possible. For U/t = 8 (panel (d)), there is no overlap between the two spectra and the system evolves into a nonconducting state. This also explains why the phenomenon of negative differential conductance is robust against an additional positive onsite bias. Any further increase of the onsite energy separates the two energy spectra of the left and the right sides even farther so the nonconducting state remains. Importantly, this energetic mechanism applies to other types of setups for studying transport phenomena. For example, by introducing a step-function bias to the lattice potential and tuning the bias, there is also a conducting-nonconducting transition of the same nature (see Appendix  A).
However, other possibilities may occur in interacting systems for large U/t when a homogeneous conducting state is no longer favorable. From Eq. (6) and the discussion below Eq. (8), states could be generated that have an inhomogenous density which would give enough "energetic relief" to allow a macroscopic current to flow [31]. In the mean-field solution, however, such current-carrying state could not be found. Here we emphasize that the nonconducting state is dynamically generated and thus is very different from the widely-discussed Mott insulating phase at integer fillings, which is an equilibrium ground state. We now consider a higher-order approximation and also td-DMRG simulations.

IV. HIGHER-ORDER CORRELATIONS AND TD-DMRG
In order to determine the accuracy of the mean-field simulations above, we also perform simulations with higher-order correlations and with td-DMRG [11]. In the higher-order approximation, we keep two-particle correlations, rather than truncating at the single-particle level. This procedure is described in Appendix B. The details of the td-DMRG simulations are given in their respective figure captions. We note that td-DMRG has been applied to study the dynamics of different cold-atom systems [32]. Figure 5 shows the results for U/t = 2, N = 40, and N pσ = 20 for the three different methods. They all agree qualitatively, giving rise to a quasi steady state. Keeping track of the higher-order correlations significantly improves the accuracy of the  (labeled next to each curve) from td-DMRG with MPS dimension 500. Simulations with higher MPS dimension indicate that these are accurate to within about 5%. The averaged current is computed from t = 5.0t0 to t = 10.0t0. Similar to Figure 3, the average current initially increases with U and then decreases, thus giving negative differential conductance. The td-DMRG simulations indicate that a conducting-nonconducting transition can exist for large U for fillings larger than half filling.
simulation. Furthermore, the results suggest that for small U/t the transient time to reach a quasi steady state is only a few t 0 while the duration of the QSSC can extend much longer. Thus for reasonably large lattices, the observation of an interactioninduced QSSC and the conducting-nonconducting transition is feasible. Figure 6 shows the average current versus U . We note that inclusion of higher order correlations limits our analysis to shorter times than above. This means that the averaging partially truncates oscillations and the effect of this partial truncation is to give a finite value to I σ even when a longer average would give zero. This will be more prominent for shorter averages, i.e., in the absence of transient effects, one expects that the truncation will give a residual contribution of I σ ∝ 1/U , as U gives the frequency scale of oscillations (for large U ). The inset of Figure 6 gives this tail: for filling greater than 1/2, the average current is consistent with this interpretation. For filling equal to 1/2, the results are ambiguous. Most importantly, all methods clearly display the negative differential conductance, and, thus, this phenomenon should be observable in experiments.
If the initial state is not far from a band-insulator (see the cases with N = 40 and N pσ = 30, 35), the td-DMRG simulations show a conducting-nonconducting transition for large U . Thus this dynamically generated conducting-nonconducting transition should be observable in experiments. For the half-filling initial state, a current carrying state is observed within our simulation range, but the magnitude of the current in this state is decreasing as U is increased as happened with the higher order correlation case.

V. CONCLUSION
In summary, we have shown that one may use time-dependent inhomogeneous interactions of ultra-cold atoms to explore non-equilibrium physics not easily realizable in conventional solid-state setups. For weak-interactions, the response of the system is similar to condensed matter systems where a bias is applied. This response is robust against possible transient behavior of the driving force. The application of inhomogeneous interactions gives rise to negative differential conductance, which is a many-body, atomic analog of this phenomenon in solid-state systems. Furthermore, a dynamically generated conductingnonconducting transition is predicted from different simulations and its observation in experiments could provide another example of nonequilibrium phase transitions. Our work sheds light on the physics of complex systems out of equilibrium and may help in the design of devices in the thriving field of atomtronics [33].
CCC acknowledges the support of the U. S. DOE through the LANL/LDRD Program. MD acknowledges support from the DOE grant DE-FG02-05ER46204 and UC Laboratories.  As with the interaction case, energy conservation prohibits particles on the left with energy higher than the highest energy level of the right from entering the right half of the lattice. This occurs when µ B is larger than the bandwidth of the original uniform lattice, which is 4t. The particle distributions of the two half lattices when the interaction is switched on can be found as follows. The initial correlation matrix c ij (t = 0) of the whole lattice is Np p=1 U † 0 (i, p)U 0 (p, j), where U 0 is the unitary matrix that diagonalizes the hopping Hamiltonian. After finding U L and U R that diagonalize the Hamiltonian of the left and the right half lattices with eigenvalues E L and E R , the particle distributions are given by n L/R (k) = i,j∈L/R (U † L/R ) kj U ik c ij (t = 0), where k = 1, · · · , N/2 labels the energy levels E L/R (k). Fig. 7 (b) and (c) show the particle distributions for the left and right half lattices for µ B /t = 3 and 4 with N p = 64 and 96. One can see that when the bottom of the left half energy band moves above the top of the band of the right, no current is expected since energy is conserved.
Taking the same definition for I σ for the interaction-induced case, we plot I σ as a function of U in Fig. 8. At the meanfield level there is also a conducting-nonconducting transition, but the threshold value depends on the filling N pσ /N . The left half lattice has Hamiltonian H L = ij ∈L,σ c † iσ c jσ + U i∈Ln iσniσ . Heren i = c † iσ c iσ . To find the energy levels of H L at t = 0, we replacen iσ by its expectation value, which is known from c ij,σ (t = 0) and we assume the particle number of the two species are equal. After diagonalizing H ef f L = ij ∈L,σ c † iσ c jσ + U i∈L n iσniσ , we found E L (k) with k = 1, · · · , k and the unitary matrix U L . The right half lattice remains noninteracting and one can find E R (k) and U R . The distribution functions are then given by n L/R (k) = i,j∈L/R (U † L/R ) kj U ik c ij (t = 0). Fig. 8 (b) and (c) show the distributions of the right and left half lattices for N pσ = 64 and 96 and N = 128. The nonconducting states emerge slightly above U/t = 7 (U/t = 5) for N pσ = 64 (N pσ = 96). The distributions for U/t below the threshold value still show observable overlaps between the filled states of the left half lattice and the empty states of the right half lattice. When U/t is above the threshold value, the distributions of the two half lattices are separated so energy conservation forbids transport at the mean-field level. We thus observe a nonconducting state due to energy-level mismatches of the two half lattices in bias-induced transport, where they are exact solutions, and in interaction-induced transport at the mean-field level.