Thermalization and temperature distribution in a driven ion chain

We study thermalization and non-equilibrium dynamics in a dissipative quantum many-body system -- a chain of ions with two points of the chain driven by thermal bath under different temperature. Instead of a simple linear temperature gradient as one expects from the classical heat diffusion process, the temperature distribution in the ion chain shows surprisingly rich patterns, which depend on the ion coupling rate to the bath, the location of the driven ions, and the dissipation rates of the other ions in the chain. Through simulation of the temperature evolution, we show that these unusual temperature distribution patterns in the ion chain can be quantitatively tested in experiments within a realistic time scale.


I. INTRODUCTION
Many-body non-equilibrium dynamics has raised significant interest in recent years, in particular with connection to the atomic experiments where far-fromequilibrium phenomena can be conveniently investigated due to the long relaxation time of these systems [1,2]. For instance, the dynamics after a quantum quench has been studied in several model Hamiltonians [3][4][5][6]. The physics becomes even richer if we add engineered dissipation to the underlying system, which is possible to realize in atomic experiments. Several interesting effects have been analyzed recently from interplay between interaction and dissipation in cold atomic gas [7,8].
Motivated by this line of research, in this paper we study thermalization and temperature distribution in a linear ion chain, with two points of the chain driven respectively by a heating and a cooling thermal bath. In a sense, the configuration here is an analog of a classical example in thermodynamics -the heat propagation in a bar with its two ends fixed at different temperature. In contrast to a simple linear temperature gradient as one sees in a classical bar, we find that the temperature distribution in a driven ion chain shows very rich behaviors in its steady states: first, the temperature distribution depends critically on the ratio between the driving speed and the interaction rate in the ion chain. In the weak driving region, the temperature distribution is non-monotonic across the ion chain and shows a strange mirror effect. In the strong driving region, instead of a linear temperature distribution, all the ions between the two driven ones are stabilized to a constant temperature which is in the middle of the two bath temperatures. Second, we find that the bath which drives the ion with a medium coupling rate plays a more effective role in determining the temperature of the other ions in the system. When the bath drives the ion strongly, it has little influence on temperature of the other ions in the chain, which is pretty counter-intuitive. To see the unusual phenomena predicted in this paper, we discuss the requirements for an experimental observation and calculate the time dynamics to achieve the steady-state temperature dis-tribution. The qualitative features of the steady-state temperature distribution patterns are pretty insensitive to the size of the ion crystal and they show up already in a pretty small ion chain with less than ten ions. The realization time to the steady state is realistic for observation compared with the current experimental time scale.
The paper is arranged as follows: In Sec. II we provide the theoretical model and the calculation method. The main results of the temperature distribution patterns under different circumstances are presented and discussed in Sec. III. In Sec. IV, we discuss the time dynamics to achieve the steady-state temperature distribution. Sec. V summarizes our major findings.

II. MODEL
We consider a chain of ions along the axial z-direction coupled with the Coulomb interaction and driven individually by a thermal bath. The thermal bath can be provided, for instance, by cooling or heating laser beams shined on each ion. Both the driving rate and the effective bath temperature can be tuned by controlling the intensity and the detuning of these laser beams. The oscillation of the ith ion around its equilibrium position is described by the coordinate and the momentum operators x i , p i , with the dynamics determined by the Heisenberg-Langevin equation: where A ij denotes the coupling matrix between the ions, γ i is the driving rate of the bath, and ζ i (t) denotes the corresponding random force from the thermal bath. We consider in this paper the ions' motion along the transverse x direction, so the coupling matrix A ij is given by where ω x is the transverse trapping frequency and z i denotes the ion equilibrium position along the axial direction. We take the ion spacing d 0 as the length unit, e 2 /d 0 as the energy unit, and e 2 /(md 3 0 ) as the fre-quency unit so that every quantities in the matrix A ij become dimensionless. Note that with control of an anharmonic trapping potential along the axial direction, one can make the ion spacing uniform in a scalable trap [9]. For the conventional harmonic axial trap, the ion spacing is not uniform. In that case, the equilibrium positions z i are determined numerically with the given trapping potential, and we take the unit d 0 as the largest spacing in the ion chain. The temperature distribution pattern that will be shown in Sec. III is insensitive to the details of the axial trapping potential, and for most of the calculations in the following, we assume the ion spacing is uniform for simplicity (with exceptions in Sec. III for consideration of the time dynamics in a harmonic trap). For independent Markovian bath, the random force ζ i (t) can be expressed as i is the average phonon number which characterizes the temperature of the bath, and ω i √ A ii is the local oscillation frequency of the ith ion by fixing all the other ions in their equilibrium positions. For typical linear ion traps with strong transverse confinement, we have ω i ≈ ω x . The correlation of ζ i (t) is then given by . The Langevin equation (1) can be solved exactly through diagonalization, with the solution formally expressed as where matrix, which can be diagonalized as U −1 ΩU αβ = λ α δ αβ . From this formal solution, we obtain the variance of the operators x i and p i as where µ = 1, 2, ..., N correspond to the x-operators and µ = N + 1, N + 2, ..., 2N correspond to the p-operators.
The temperature of each ion is similarly characterized by the average phonon number of its thermal oscillation, so the temperature of the ith ion as a function of time t is represented by As long as γ i is nonzero for some ions, each eigenvalue λ α has a positive real part due to coupling of all the ions, and the system approaches to a steady state as t −→ ∞. The steady state temperature for each ion is denoted by

III. STEADY-STATE TEMPERATURE DISTRIBUTION
In analogy to the example of heat propagation in a conducting bar, we consider a long ion chain with the two edge ions driven by different thermal bath with temperature T B 1 and T B 2 (with T B 1 < T B 2 ), respectively. For simplicity, we assume the corresponding driving rate γ 1 = γ 2 = γ and all the ions in the middle have no dissipation. After all the ions attain the steady state, the temperature distribution T s i across the chain is shown in Fig. 1 for N = 100 ions under different driving rates γ. Apparently, the temperature distribution does not follow a linear gradient. There are three regions for the distribution of T s i , depending on the ratio between the driving rate γ and the ion interaction rate. Note that in our unit the interaction energy between the neighboring ions is of the order of unity. If γ ≫ 1, the dissipation is much faster than the energy propagation in the chain, and without surprise the two edge ions have temperature basically fixed by their corresponding bath temperature T B 1 and T B 2 . However, it is surprising that all the other ions in the middle approach almost the same temperature given by T s i ≃ (T B 1 + T B 2 )/2 for i = 2, 3, · · · , 99 in this case. The temperature does not fall down gradually from the hot end to the cold end as in the classical heat propagation problem, but takes a sharp jump right from the driven ion to the next one and then keeps constant over the whole chain. In the opposite limit of weak driving with γ ≪ 1, the phonon propagation is faster than the bath driving, and all the ions approach the same temperature. It sounds that the temperature distribution in this limit resembles a classical thermal equilibrium, however, this picture is not true. The example in Fig. 1 represents an exception instead of a rule where we put the cooling and heating ions exactly at the symmetric positions of an ion chain. As we will see in the following, when we shift the position of one of the ions to break the reflection symmetry, the temperature distribution in the weak driving limit has strange features. It is even not monotonic across the ion chain, in sharp contrast with the distribution from the diffusion process. Between these two limiting regions there lies a transition region with 0.01 < γ < 1, where the temperature of the edge ions gradually approach the corresponding bath temperature. Note that the temperature for the ions next to the driven ones follow a non-monotonic curve as one changes γ, although such variation is pretty small.
To break the position symmetry of the cooling and the heating ions, we fix the cooling ion still at the edge, but move the heating ion inside the chain. Fig. 2(a) shows an example of the temperature distribution where the heating ion is right at the center of the chain. When γ ≫ 1, the heating ion separates two plateaus in the tempera-ture distribution: a lower temperature region between the heating and the cooling ions and a higher temperature region on the other side of the heating ion. As γ decreases, the two plateaus smoothly descend to the same temperature while the heating ion remains as a peak in the temperature distribution. When γ ≪ 1, the temperature distribution has a reflection symmetry with respect to the chain center. A dip in temperature appears on the free edge of the chain (the left side of Fig. 2(a)), mirroring the cooling ion on the other side of the chain (the mirror effect). The temperature distribution is apparently non-monotonic in this case. The mirror effect exists for other positions of the heating ion. For instance, Fig.  2(b) shows the temperature distribution where the distance between the heating and the cooling ions is about one third of the chain length. Both the cooling and the heating ions have their mirror images in the temperature distribution in the weak driving limit.
In the above calculation, we assumed γ 1 = γ 2 for the cooling and the heating ions. When these two driven ions are put at the edge, the middle ions are at a constant temperature T m which is exactly the average of the two bath temperatures. If the driving rate γ 1 = γ 2 , we may expect that the bath which drives the ion more strongly plays a larger role in determining the temperate T m of the middle segment ions. This expectation, however, turns out to be not true. To show that, we map out the middle segment temperature T m in Fig. 3 as a function of γ 1 and γ 2 . The figure shows that γ i ≈ 0.1 is the optimal driving rate for which the corresponding bath has the largest influence on the temperature T m . At this optimal value, the driving rate is comparable with the ion interaction rate in term of the order of magnitude. The middle segment temperature is close to the temperature of the bath that drives the ion at the opti-   Fig. 1, as a function of two independent bath driving rates γ1 and γ2. The dark (light) regions correspond to low (high) temperatures, and the numbers there indicate different temperature contours. The parameters are the same as in Fig. 1 expect that γ1 associated with the cooling bath(at T B 1 = 2) and γ2 associated with the heating bath (at T B 2 = 10) are varied independently. The temperature profiles at different γ bg (cross-section view of (a)). We find that the temperature profile for the main part of the chain, excluding the two edge ions, approaches a linear distribution at an intermediately small γ bg ≃ 10 −2 γ. mal rate. When the driving gets too strong or too weak, the bath plays little role in determining the temperature of the other ions in the chain. This result has important implication for the sympathetic cooling [10][11][12]: instead of fast cooling of the ancilla ions, cooling at a moderate optimal rate is more efficient in reducing the temperature of the computational ions.
So far we have neglected dissipation of the middle segment ions. Now, apart from the cooling and the heating bath (with temperature T B 1 and T B 2 , respectively) attached to the two edge ions, we assume all the middle ions are coupled to a background bath with temperature at a coupling rate γ bg . The final temperature distribution of the ions is shown in Fig. 4 under different background coupling rates γ bg . In this calculation, we fix γ 1 = γ 2 = γ = 0.1. For tiny γ bg , the temperature distribution shows no noticeable difference compared with the case of γ bg = 0. However, as γ bg increases to a moderate value with γ bg /γ ∼ 0.01, the temperature distribution of the middle segment ions has a clear linear spatial gradient, resembling the linear temperature distribution of a classical bar in the heat diffusion problem. As γ bg further increases and gets close to γ, the temperature of the middle ions are pined to the background bath temperature T B bg as one expects and the temperature gradient disappears again in this limit.

IV. TEMPERATURE EVOLUTION
We note that the temperature distributions shown in this paper are insensitive to the size of the ion system. We have checked the temperature distribution of the ions with the size of the chain varying from ten to a few hundreds of ions and noticed no qualitative change of the distribution pattern. To experimentally test these unusual temperature distribution patterns, one can start with a small system of a few ions that are within the reach of the current experimental technology. The final tempera- The initial temperature is assumed to be Ti(t = 0) = 5 for all the ions. The transverse frequency ωx = 10 in the unit of e 2 /(md 3 0 ), where d0 is the largest spacing in the ion chain (between the 1 st and the 2 nd , or the 19 th and the 20 th ions) for the harmonic cases. In terms of real numbers, with a typical choice of d0 = 10µm for 20 ytterbium ions in a harmonic trap, the trapping frequencies for the transverse and the axial traps are given respectively by ωx ≈ 2π × 1.4MHz and ωz ≈ 2π × 76.5kHz, and the driving rate γ ≈ 14kHz. ture (the mean phonon number of the ion motion) can be detected, for instance, by measuring the scattering sidebands of a laser beam through a ion. The asymmetry in the blue and the red sidebands and their ratio give direct inference of the thermal phonon number of the ion motion [13][14][15].
To probe the properties of the steady states, we need to know the time scale to approach these steady states. Fig. 5 shows the relaxation dynamics for 20 ions in either a harmonic or an anharmonic uniform trap [9]. For a uniform trap (in Fig. 5a), the calculation shows that there are two time scales in the thermalization. First, with a time scale t 1 ∼ 1/γ, the temperature of the two edge ions quickly approaches the corresponding bath temperature. During this step, the temperature of the middle segment ions only changes slightly, and the change gets smaller as one moves away from the driven ions. After that, a longer time scale t 2 sets in, representing the interaction driven thermalization process. All the ions gradually approach the steady-state temperature. Note that the temperature of the edge ions (as a well as the other ions that are close to the two driven ones) does not follow a monotonic evolution curve. Instead, it first comes pretty close to the corresponding bath temperature and then is dragged back toward its steady state value. The value of the second time scale t 2 increases with the system size (roughly linearly) and t 2 ∼ 40/γ for 20 ions. For a ytterbium ion ( 171 Yb + ) chain with spacing 10 µm, t 2 is about 3 ms, which is a pretty reasonable time scale for experiments. For a harmonic trap (in Fig. 5b), the basic feature is similar except that the second time scale t 2 gets site dependent. The edge ions and their neighbors approach the corresponding steady-state temperature with a time scale that is comparable with the case of a uniform chain. However, as one moves away from the edge, the thermalization time gets much longer with an exponential increase. For the middle ion, the thermalization is not finished yet with t ∼ 10 10 /γ. With such an extremely long time scale, of course one can not neglect the small dissipation of the other ions to the background bath. If we take into account a small background dissipation, e.g., with a rate γ bg /γ ∼ 10 −3 as shown in Fig. 5c, the long tale in the temperature evolution is completely gone for the harmonic trap. Now, the temperature of all the ions approach their steady values within a time scale that is comparable with the uniform case, although the steadystate temperature of the middle segment ions is dragged down a bit toward the background bath temperature.

V. CONCLUSION
In summary, we have shown that the steady state temperature distribution of a driven ion chain shows surprisingly rich patterns. Many features of these patterns are unexpected, and they are in sharp contrast with the simple linear temperature gradient as one sees in the classical heat diffusion problem. Our calculation is based on exact numerical methods with no unreliable approximations, so we believe all the unusual temperature distribution patterns revealed by this calculation will show up in experiments, although we still lack intuitive explanation of many of these features. We also investigate the relaxation dynamics and the time scale to reach the steady state, and show that these patterns should be observable under a realistic time scale in a small system that is within the current experimental reach.