Thermoelectricity modeling with cold dipole atoms in Aubry phase of optical lattice

We study analytically and numerically the thermoelectric properties of a chain of cold atoms with dipole-dipole interactions placed in an optical periodic potential. At small potential amplitudes the chain slides freely that corresponds to the Kolmogorov-Arnold-Moser phase of integrable curves of a symplectic map. Above a certain critical amplitude the chain is pinned by the lattice being in the cantori Aubry phase. We show that the Aubry phase is characterized by exceptional thermoelectric properties with the figure of merit ZT = 25 being ten times larger than the maximal value reached in material science experiments. We show that this system is well accessible for magneto-dipole cold atom experiments that opens new prospects for investigations of thermoelectricity.


I. INTRODUCTION
The phenomenon of Aubry transition describes the transport properties of a chain of particles linked by linear springs in a periodic potential. At a small potential amplitude the chain slides freely while above a certain potential amplitude it is pinned by the potential. This system is known as the Frenkel-Kontorova model [1]. The transition takes place at a fixed incommensurate density of particles per period ν. In fact the equilibrium positions of particles are described by the Chirikov standard map [2][3][4] which represents a cornerstone model of systems with dynamical chaos and symplectic maps. Indeed, a variety of physical systems can be locally described by this map [5]. In the frame of map description a density of particles corresponds to a winding number of an invariant Kolmogorov-Arnold-Moser (KAM) curve. Such curves cover the main part of the phase space at small potential amplitudes (small kick amplitudes K of the map). At large amplitudes the main part of phase space becomes chaotic and the KAM curves are transformed into cantori invariant sets which correspond to the chain ground states at given density as it was proved by Aubry [6]. In addition to nontrivial mathematical properties, the Frenkel-Kontorova model represents a fundamental interest for incommensurate crystals of interacting particles [7].
An experimental realization of linear spring interactions between particles is not very realistic. Thus in [8] it was proposed to consider a chain of Coulomb charges placed in a periodic potential. It was shown that this system of a one-dimensional (1D) Wigner crystal in a periodic potential can be locally described by the Chirikov standard map and the Frenkel-Kontorova model. Thus the Aubry-like transition from the sliding KAM phase to the Aubry pinned phase takes place at a certain critical potential amplitude K c . The dependence of K c on the density ν is obtained by a local description in terms of the Chirikov standard map [8,9]. For an experimental realization of the Aubry-like transition it was proposed to use cold ions placed in both a periodic potential and a global harmonic trap [8]. The experimental studies of such a system had been started in [10,11]. The first signatures of the Aubry transition have been reported by the Vuletic group with up to 5 ions [12]. The chains with a larger number of ions are now under investigations in [13,14].
A significant interest for a Wigner crystal transport in a periodic potential is related to the recent results showing that the Aubry phase is characterized by remarkable thermoelectric properties with high Seebeck coefficient S and high figure of merit ZT [9,15]. The fundamental grounds of thermoelectricity had been established in far 1957 by Ioffe [16,17]. The thermoelectricity is characterized by the Seebeck coefficient S = −∆V /∆T (or thermopower). It is expressed through a voltage difference ∆V compensated by a temperature difference ∆T . Below we use units with a charge e = 1 and the Boltzmann constant k B = 1 so that S is dimensionless (S = 1 corresponds to S ≈ 88µV/K (microvolt per Kelvin)). The thermoelectric materials are ranked by a figure of merit ZT = S 2 σT /κ [16,17] where σ is the electric conductivity, T the temperature, and κ the material thermal conductivity.
At present, the request for efficient energy usage stimulated extensive investigations of various materials with high characteristics of thermoelectricity as reviewed in [18][19][20][21][22]. The request is to find materials with ZT > 3 that would allow an efficient conversion between electrical and thermal forms of energy. The best thermoelectric materials created till now have ZT ≈ 2.6. At the same time, the numerical modeling reported for a Wigner crystal in the Aubry phase reached values ZT ≈ 8 [9,15].
Thus investigations of Wigner crystal transport in a periodic potential can help to understand the conditions favoring high ZT values. At present, a hundred of cold trapped-ions can be routinely kept for hours in experimental device [23] and thus such a system is promising for experimental investigations of thermoelectricity [9]. However, for a typical distance between charges being ∼ 1µm the Coulomb interactions are rather strong and very high amplitudes of optical lattice potential V A ∼ k B × 3K (Kelvin) are required [9]. This is hard to reach experimentally since typical optical potential amplitudes are V A ∼ k B × 10 −3 K [24]. Thus to find a more suitable experimental realization of Aubry transition we study here a chain of magneto-dipole atoms placed in an optical periodic potential. The strength of interactions between nearby magneto-dipole atoms on a distance of 1µm is significantly smaller compared to Coulomb interactions, and thus a significantly smaller amplitude of the optical potential is required for the observation of the Aubry-like transition. Indeed, the experimental investigations of quantum properties of cold magneto-dipole atoms allowed to observe a number of interesting manybody effects (see e.g. [25][26][27]).

II. METHODS
The chain of atoms with magneto-dipole interactions in 1D periodic potential is described by the Hamiltonian: Here x i , P i are conjugated coordinate and momentum of atom i, and V (x) is an external periodic optical potential of amplitude K. The magneto-dipole-dipole interactions are given by the U I ∝ U dd = µ 2 /( /2π) 3 term with µ ≈ 10µ B (for 164 Dy atoms) and µ B is the Bohr magneton (we assume all magnetic momenta to be polarized) [25]. The Hamiltonian is written in dimensionless units where the lattice period is = 2π and atom mass is m = 1. In the following we also take U dd = 1 so that K in (1) is the dimensionless amplitude of the periodic potential and thus the physical interaction is U A = U dd K. In these atomic-type units, the physical system parameters are measured in units of r a = /2π for length, and of a = µ 2 /r a 3 = U dd for energy. For = 1µm the dimensionless temperature T = 1 (or k B T ) corresponds to the physical The thermoelectric properties of model (1) are determined in the framework of the Langevin approach [9,15] with the equations of motion: Here η describes phenomenologically the dissipative relaxation processes, and the amplitude of Langevin force g is given by the fluctuationdissipation theorem g = √ 2ηT ; the random variables ξ i have normal distribution being δ−correlated, v i is atom velocity, f dc is a static force applied to atoms. As in [9,15], we use η = 0.02, the results being not sensitive to this quantity. The computations of S are done as it is described in [9,15]. At fixed temperature T , we apply a static force f dc which creates an energy (voltage) drop ∆V = f dc 2πL and a gradient of atom density ν(x) along the chain with L potential periods and N atoms.
Then, for f dc = 0 within the Langevin equations, we impose a linear gradient of temperature ∆T along the chain, and in the stabilized steady-state regime, we determine the charge density gradient of ν(x) along the chain (see e.g. Fig. 2 in [15]). The data are obtained in the linear regime of relatively small f dc and ∆T values. Then, the Seebeck coefficient, S = ∆V /∆T , is computed using values of ∆V and ∆T for which the density gradient obtained from the application of a voltage ∆V compensates the one obtained from the application of a gradient of temperature ∆T . We used the computation times up to t = 10 8 to achieve the relaxation of the chain and to reach the required statistical accuracy. We assume that the cold atoms are in contact with an external environment which is able to play the role of thermostat, e.g. residual gas, etc. In order to compute Seebeck coefficient, we need a temperature gradient along the chain. In the Langevin equation, we can impose that the temperature T is a function of the atom position x along the chain, T = T (x) = T 0 + gx, where T 0 is the average temperature and g = dT /dx is a small temperature gradient. In cold atom experiments, such a temperature gradient can be setup by multiple laser beams generating a zero average fluctuating force f (x). The average of the square of the force, f 2 (x), should change linearly along the chain. Consequently, these laser beams induced fluctuating forces will create additional ion velocity fluctuations with (δv i ) 2 ∝ f 2 (x) ∝ (T (x) − T 0 ) = gx producing a temperature gradient along the ion chain.

A. Ground state properties
The equilibrium static atom positions are determined by the conditions ∂H/∂x i = 0, P i = 0 [6,8,9]. In the approximation of nearest neighbor interacting atom, this leads to the symplectic map for recurrent atom positions where the effective momentum This map description assumes only nearest neighbor interactions. Below, we show that it well describes the real situation with interactions between all the atoms. This is rather natural since the nearest neighbor approximation worked already well for ions with Coulomb interactions [8,9], moreover for dipole atoms, the interactions drop even faster with the distance between atoms. As in [9,15], the validity of the map (2) is checked numerically by finding the ground state configuration using numerical methods of energy minimization described in [6,8] and taking into account the long range nature of the interactions between all the atoms. The obtained ground state positions, {x 1 , . . . , x N }, of the N atoms allows to recursively determine the N effective momenta, Once obtained, the effective momenta can be used to compute successively N values {g(x 1 ), . . . , g(x N )} of the kick function g(x), by using the first equation of the map (2) (2) is given in Fig. 1 (middle panels). A phase space portrait (green points and curves) is obtained by taking many different initial conditions (x 0 , p 0 ) and computing many of the corresponding successive phase space points (x i , p i ). In Fig. 1 (middle panels), we show the phase space portrait for KAM phase (K = 0.02 < K c , left panel) and for the Aubry phase (K = 0.1 > K c , right panel). We observe that the ground state positions of the atoms (red circles) are located, for K = 0.02 (left panel), on the top of an invariant KAM curve which is a situation corresponding to a regular motion of the fictitious kicked particle, and, for K = 0.1 (right panel), in the chaotic component of the phase space portrait which is a situation corresponding to a chaotic motion of the fictitious kicked particle. All the chaotic dynamics concepts used in the article are well documented in the literature (see e.g. [3,4]). Fig. 1 (bottom panels) shows the hull function h(x) = (x i + π)[mod 2π] − π with x = (2π(i − 1)/ν + π)[mod 2π] − π. Indeed, for K → 0, i.e. for a vanishing optical potential, we have h(x) = x. For K < K c , the hull function h(x) still follows this linear law with smooth deviations, while, for K > K c , we obtain the devil's staircase corresponding to the fractal cantori structure of the chaotic phase space. The transition from the smooth hull function, which is typical for a sliding chain in the KAM phase, to the devil's staircase, which is typical for a pinned chain in the Aubry phase, is visible in Fig. 1 (bottom panels). Thus, with ν = 89/55 ν g = (1 + √ 5)/2 = 1.618..., i.e. 89 atoms (in their ground state) distributed in 55 optical potential wells, the Aubry transition takes place at a certain K c inside the interval 0.02 < K c < 0.1.
The equations of motion can be linearized in a vicinity of equilibrium positions and, in this way, we obtain the phonon spectrum ω(k) of small oscillations with k = i/N being a scaled mode number. The examples of spectrum are shown in Fig. 2 (left panel). At K = 0.02, in the KAM phase, we have ω ∝ k corresponding to the acoustic modes, while at K = 0.1, inside the Aubry phase, we have appearance of an optical spectral gap related to the atomic chain pinned by the potential. Such a modification of the spectrum properties is similar to the cases of the Frenkel-Kontorova model [6] and the ion chain [8,9]. Fig. 2 (right panel) shows that the minimal spectral frequency ω 0 (K) is practically independent of the optical potential amplitude K inside the KAM phase at K < K c (being close to zero with ω 0 ∝ 1/L) and it increases with K inside the Aubry phase being independent of L for K > K c . Thus, the critical K c values can be approximately determined as an intersection of a horizontal line ω 0 = const with a curve of growing ω 0 (K) at K > K c . From these properties, we obtain numerically that K c ≈ 0.019 for the Fibonacci density ν ≈ ν g = 1.618... We obtained also K c ≈ 0.14 for ν = 2.618, and K c ≈ 0.4 for ν = 3.618.

B. Density dependence of Aubry transition
The dependence of Aubry transition point K c (ν) can be obtained from the local description of the map (2) by the Chirikov standard map. For that, the equation of x i+1 in (2) is linearized in p i+1 near the resonant value p r that leads to the standard map form y i+1 = y i − K eff sin x i , x i+1 = x i − y i+1 with y i ∝ p i and K eff = (2π) 5 K/(12ν 5 ) (see details in Appendix A and [9]). The last KAM curve is destroyed at K eff ≈ 1 [2,3] that gives K c ≈ 12(ν/2π) 5 ≈ 0.0136(ν/ν g ) 5 , ν g = 1.618... (3) The numerically obtained dependence K c (ν) is shown in Fig. 3 for 3 different chain lengths L. On average, it is well described by the theory (3) taking into account that K c is changed by almost 4 orders of magnitude for the considered range 0.5 ≤ ν ≤ 3.5. There are certain deviations in the vicinity of integer resonant densities ν = 1, 2, 3 which should be attributed to the presence of a chaotic separatrix layer at such resonances that reduces significantly the critical K c for KAM curve destruction (similar effect has been discussed for the ion chains [9]).

C. Seebeck coefficient
The dependencies of S on potential amplitude and temperature are shown in Fig. 4 for two Fibonacci-like values of density ν ≈ 1.618 and ν ≈ 2.618. The results clearly show that in the KAM phase K < K c we have only rather moderate values of S ∼ 1 being close to those value of noninteracting particles (similar result was obtained in [9,15]). In the Aubry phase at K > K c , we have an increase of S with the increase of K/K c , and a decrease of S with the increase of T /K c . This is rather natural since at T K c the lattice pinning effect disappears due to the fact that the atom energies become significantly larger than the barrier height, and we approach to the case without potential corresponding to the KAM phase. The maximal obtained values are as high as S ≈ 10 − 15 being still smaller those obtained for ion chains [9,15]. Nevertheless, as shown below, we find very large figure of merit ZT values in this regime.

D. Figure of merit
To obtain the value of ZT we need to compute the conductivity of atoms, σ, and their thermal conductivity, κ. The value of σ is defined through the current J of atoms induced by a static force f dc for a chain with periodic boundary conditions: j = νv at /2π where v at is the average velocity of atoms and σ = j/f dc ; for K = 0, we have σ = σ 0 = ν/(2πη). The heat flow J is induced by the temperature gradient due to the Fourier law with the thermal conductivity κ = J/(∂T /∂x). The heat flow is computed as it is described in [9,15] and in Appendix B. The values of σ and κ drop exponentially with the increase of K inside the Aubry phase at K > K c (see Figs. 6 and 7 in Appendix B). The computations performed for various chain lengths at fixed density confirm that the obtained values of S, σ, κ are independent of the chain length for T > K c confirming that the results are obtained in the thermodynamic limit (see Fig. 8 in Appendix C). For T K c , the pinning is too strong and much larger computation times are needed for numerical simulations to reduce the fluctuations.
We note that, for cold atoms in an optical lattice, a static force can be created by a modification of the lat-tice potential or by lattice acceleration. There is now a significant progress with the temperature control of cold ions and atoms (see e.g. [28][29][30]) and we expect that a generation of temperature gradients for measurements of κ and S can be realized experimentally.
With the computed values of σ, κ, S , we determine the figure of merit ZT . Its dependence on K and T is shown in Fig. 5 for ν ≈ 1.618 and 2.618 (additional data are given in Figs. 9 and 10 in Appendices D and E). For ν ≈ 1.618, we obtain the maximal values ZT ≈ 5 being comparable with those of ion chain reported in [15]. However, for ν ≈ 2.618, we find significantly larger maximal values with ZT ≈ 25. We attribute such large ZT values to the significantly more rapid spatial drop of dipole interactions comparing to the Coulomb case, arguing that this produces a rapid decay of the heat conductivity with the increase of K > K c .

IV. DISCUSSION
The obtained results show that a chain of dipole atoms in a periodic potential is characterized by outstanding values of the figure of merit ZT ≈ 25 being by a factor ten larger than the actual ZT values reached till present in material science [22]. Thus, the experiments with cold dipole atoms in the Aubry phase of an optical lattice open new prospects for experimental investigation of fundamental aspects of thermoelectricity.
We note that for a laser wavelength λ = 564nm, the optical lattice period is = λ/2 = 282nm, and thus for ν ≈ 2.6, we have the Aubry transition at the potential amplitude V A /k B = T A = 0.14 a /k B ≈ 200nK. Such potential amplitude and temperature are well reachable with experimental setups at T ≈ 20nK used in [27]. At the same time at T ≈ 200nK the wave length of Dy atoms becomes λ Dy = / 2m Dy k B T ≈ 90nm < being only a few times smaller than the lattice period . Thus the quantum effects can start to play a role. But their investigations require a separate study.
Since such Aubry phase parameters are well accessible for experiments, it may be also interesting to test the quantum gate operations of atomic qubits, formed by two atomic levels, with dipole interaction between qubits. Such a system is similar to ion quantum computer in the Aubry phase discussed in [31]. As argued in [31,32], the optical gap of Aubry phase should protect gate accuracy.

V. ACKNOWLEDGMENTS
We thank D.Guery-Odelin for useful discussions about dipole gases. This work was supported in part by the Pogramme Investissements d'Avenir ANR-11-IDEX-0002-02, reference ANR-10-LABX-0037-NEXT (project THETRACOM). This work was also supported in part by the Programme Investissements d'Avenir ANR-15-IDEX-0003, ISITE-BFC (project GNETWORKS), and puted as described in [9,15]. Namely, it is computed from forces acting on a given atom i from left and right sides being respectively For an atom moving with a velocity v i , these forces create left and right energy flows In a steady state, the mean atom energy is independent of time and J L + J R = 0. But, the difference of these flows gives the heat flow along the chain: Such computations of the heat flow are done with fixed atom positions at chain ends. In addition, we perform time averaging using accurate numerical integration along atom trajectories that cancels contribution of large oscillations due to quasi-periodic oscillations of atoms. In this way, we determine the thermal conductivity via the relation κ = 2πJL/∆T . The obtained results for κ are independent of small ∆T . It is useful to compare κ with its value κ 0 = σ 0 K c . The dependence of κ/κ 0 on K at different T is shown in Fig. 7.  Fig. 6. Here σ0 = ν/(2πη), κ0 = σ0Kc.
Appendix C: Independence of chain length Here in Fig. 8, we present results at different chain lengths for fixed atom density showing that the Seebeck coefficient is independent of the system size. We found similar independence of the chain length for σ, κ and power factor P S = S 2 σ/σ 0 .  Here Kc = 0.019 for ν = 34/21 and Kc = 0.14 for ν = 55/21.