Asymmetric Heat Transport in Ion Crystals

We numerically demonstrate heat rectification for linear chains of ions in trap lattices with graded trapping frequencies, in contact with thermal baths implemented by optical molasses. To calculate the local temperatures and heat currents we find the stationary state by solving a system of algebraic equations. This approach is much faster than the usual method that integrates the dynamical equations of the system and averages over noise realizations.


I. INTRODUCTION
The ideal thermal rectifier, also "thermal diode", is a device that allows heat to propagate in one direction, from a hot to a cold bath, but not in the opposite one when the temperature bias of the baths is reversed. The name is set by analogy to the half-wave rectifiers or diodes for electric current. More generally thermal rectification simply denotes asymmetric heat flows (not necessarily all or nothing) when the bath temperatures are reversed. Thermal rectification was discovered by C. Starr in 1936 in a junction between copper and cuprous oxide [1]. Many years later, a work of Terraneo et al. demonstrated thermal rectification in a model consisting on a segmented chain of coupled nonlinear oscillators in contact with two thermal baths at temperatures T H and T C , with T H > T C [2]. This paper sparked a substantial body of research to this day [3] (see Fig. 1 in [4]).
Research on thermal rectification has gained a lot of attention in recent years as a key ingredient to build prospective devices to control heat flows similarly to electrical currents [4,5]. There are fascinating proposals to engineer thermal logic circuits [6] in which information, stored in thermal memories [7], would be processed in thermal gates [8]. Such thermal gates, as their electronic counterparts, will require thermal diodes and thermal transistors to operate [9,10]. Heat rectifying devices would also be quite useful in nano electronic circuits, letting delicate components dissipate heat while being protected from external heat sources [4].
Most work on thermal diodes has been theoretical with very few experiments. A relevant attempt to build a thermal rectifier was based on a graded structure made of carbon and boron nitride nanotubes that transports heat between a pair of heating/sensing circuits [11]. One of the ends of the nanotube is loaded with a deposition of another material, which makes the heat flow better from the loaded end to the unloaded end. However, rectifications were small, with rectification factors (relative heat-flow differentials) around 7%.
Much effort has been aimed at improving the rectification factors and the features of the rectifiers. Some works relied, as in [2], on a chain segmented into two or more regions with different properties, but using other lattice models such as the Frenkel-Kontorova (FK) model [12,13]. The fundamental ingredient for having rectification was attributed to nonlinear forces in the chain [5,[12][13][14][15][16], which lead to a temperature dependence of the phonon bands or power spectral densities. The bands may match or mismatch at the interfaces depending on the sign of the temperature bias of the baths, allowing or obstructing heat flow [2,17]. Later, alternative mechanisms have been also proposed which do not necessarily rely on anharmonic potentials [18,19].
It was soon realized that the performance of segmented rectifiers was very sensitive to the size of the device, i.e., rectification decreases with increasing the length of the rectifier [13]. To overcome this limitation two ideas were proposed. The first one consists in using graded rather than segmented chains, i.e., chains where some physical property varies continuously along the site position such as the mass of particles in the lattice [20][21][22][23][24][25][26][27][28]. The second one uses particles with long range interactions (LRI), such that all the particles in the lattice interact with each other [21,29,30]. The rationale behind is that in a graded system new asymmetric, rectifying channels are created, while the long range interactions create also new transport channels, avoiding the usual decay of heat flow with size [21]. Besides a stronger rectification power, LRI graded chains are expected to have better heat conductivity than segmented ones. This is an important point for technological applications, because devices with high rectification factors are not useful if the currents that flow through them are very small.
In this article we propose to bridge the gap between mathematical models and actual systems exploring the implementation of a heat rectifier in a realistic, graded system with long range interactions: a chain of ultracold ions in a segmented Paul trap with graded microtraps for each ion. Long-range interactions are due to the Coulomb forces, and the baths at the ends of the chain may be implemented with optical molasses, see Fig. 1. The trapping frequencies of the microtraps are controlled individually in order to create a graded and asymmetric trap-frequency profile along the chain. This asymmetry will lead to a heat flow that depends on the sign of the temperature difference of the baths. Heat transport in trapped ion chains has been studied in several works [ [31][32][33][34] and interesting phenomena like phase transitions have been investigated [31][32][33]. The idea of using locallycontrolled traps is already mentioned in [31] to implement disorder and study its effects. The device we present here may be challenging to implement, but at reach with the current technology, in particular that of microfabricated traps [35][36][37]. The setting is thought for a small number of controllable ions, and we refrain from scaling up the size of the chain to keep the simulations realistic. The rest of the article is organized as follows. In Section II we describe the physical system of trapped ions with graded trap frequencies. We also set the stochastic dynamics due to the action of lasers at the chain edges. In Section III we implement an efficient method to find the steady state using Novikov's theorem and solving an algebraic system of equations. In Section IV we present simulations of this system exhibiting thermal rectification and discuss the advantages/disadvantages of using a graded frequency profile instead of a segmented one. Finally, in Section V we summarize our conclusions and discuss future research.

II. PHYSICAL SYSTEM
Consider a linear lattice of N individual harmonic traps of (angular) trapping frequencies ω n evenly distributed along the x axis at a distance a from each other. Each trap contains a single ion that interacts with the rest via Coulomb potentials. All the ions are of the same species, with mass m and charge q. The Hamiltonian that describes the dynamics of the system is (we consider only linear, one dimensional motion along the chain axis) where {x n , p n }, position and momentum of each ion, are the components of the vectors x, p, x (0) n = na are the centers of the harmonic traps, and V int is the sum of the Coulomb interaction potential between all pairs of ions, The ends of the chain are in contact with two thermal reservoirs at temperatures T L for the left bath and T R for the right bath respectively. The action of the resevoirs on the dynamics of the chain is modeled via Langevin baths at temperatures T L and T R [38,39]. The equations of motion of the chain, taking into account the baths and the Hamiltonian, arė where γ n and ξ n (t) are only non-zero for the ions in the end regions, in contact with the left and right baths in the sets Fig. 1. The γ n are friction coefficients and ξ n (t) are uncorrelated Gaussian noise forces satisfying ξ n (t) = 0 and ξ n (t)ξ m (t ′ ) = 2D n δ nm δ(t − t ′ ), D n being the diffusion coefficients. These Gaussian forces are formally the time derivatives of independent Wiener processes (Brownian motions) ξ n (t) = √ 2D n dWn dt [32,40] and Eq. (3) is a stochastic differential equation (SDE) in the Stratonovich sense [40].
The baths are physically implemented by optical molasses consisting of a pair of counterpropagating Dopplercooling lasers [32]. The friction and diffusion coefficients for the ions in contact with the baths are given by [41] where k L (k R ) and I L (I R ) are the wave vector and intensity of the left (right) laser. δ L (δ R ) is the detuning of the left (right) laser with respect to the angular frequency ω 0 of the atomic transition the laser is exciting, and Γ is the corresponding natural line-width of the excited state. The expressions in Eq. (4) are valid only if the intensities of the lasers are small compared to the saturation intensity I 0 , I L,R /I 0 ≪ 1. In this bath model, the friction term in Eq. (3) comes from the cooling action of the laser and the white noise force ξ n (t) corresponds to the random recoil of the ions due to spontaneous emission of photons [42][43][44]. Using the diffusion-dissipation relation D = γk B T [45], the temperature of the optical molasses baths are given by with k B being the Boltzmann constant. If the laser intensities are low enough, the temperatures of the baths are controlled by modifying the detunings. When δ = −Γ/2 the optical molasses reach their minimum temperature possible, the Doppler limit T D = Γ/(2k B ).

III. CALCULATION OF THE STATIONARY HEAT CURRENTS
The local energy of each site is defined by (6) Differentiating H n with respect to time we find the continuity equatioṅ (p n +p l ), which gives the "internal" energy flow due to the interactions with the rest of the ions. In the steady state Ḣ n = 0, and therefore where ··· stands for the expectation value with respect to the ensemble of noise processes ξ(t) (ξ represents a vector with components ξ n ). Equation (8) implies that, in the steady state, the internal ratesḢ int n vanish for the inner ions of the chain because j B n = 0 for n / ∈ L, R. In chains with nearest-neighbor (NN) interactions, Ḣ int n simplifies to two compensating and equal-in-magnitude contributions that define the homogeneous heat flux across the chain. For long-range interactions this is not so and defining the flux is not so straightforward. A formal possibility is to impose nearest-neighbor interatomic interactions for some atoms in the chain [21], but this approach is not realistic in the current system so we define instead the heat currents for the left and right baths as respectively. In the steady state we must have J L + J R = 0, since the local energies stabilize and internal energy flows cancel. We use either J L or J R to calculate the total energy flow in the chain, always taking the absolute value, i.e., J ≡ |J L | = |J R |. J is defined as J → when the hot bath is on the left and J ← when it is on the right.
To compute the average heat fluxes of the baths j B n in Eq. (9) we need the averages p n (t)ξ n (t) . Instead of explicitly averaging p n (t)ξ n (t) over different realizations of the white noise we use Novikov's theorem [40,46,47]. Novikov's theorem states that the ensemble average (over the realizations of the noise) of the product of some functional φ(t), which depends on a Gaussian noise ξ(t) with zero mean value, ξ(t) = 0, and the noise itself, is given by where δφ(t)/δξ(t ′ ) is the functional derivative of φ(t) with respect to the noise, with t ′ < t. When the noise is δ−correlated, ξ(t)ξ(t ′ ) = 2Dδ(t − t ′ ), and Eq.
To apply Novikov's theorem to our model we need the functional derivatives of the position x n (t) and momentum p n (t) coordinates with respect to the white noises. We integrate Eq. (3) to have its formal solution as a functional depending on the white Gaussian noises ξ n (t), x n (t) = x n (0) + 1 m Equation (11) implies that the functional derivatives are δx n (t)/δξ m (t ′ )| t ′ →t − = 0 and δp n (t)/δξ m (t ′ )| t ′ →t − = δ nm (δ nm is the usual Kronecker delta symbol). Thus we have x n (t)ξ m (t) = 0 and p n (t)ξ m (t) = δ nm D m , which gives for the heat flow from the baths In all simulations we check that |J L | = |J R | within the numerical tolerance of the computer. To measure the asymmetry of the heat currents we use the rectification factor R defined as R values may go from −1 to 1. If there is no rectification J → = J ← and R = 0. For perfect rectification in the right (left) direction, J → ≫ J ← (J → ≪ J ← ), and R = 1 (R = −1).
A. Algebraic, small-oscillations approach to calculate the steady state To find the temperature profiles and heat currents in the steady state the usual approach is to solve the SDE system in Eq. (3) up to long times and for many realizations of the white noises ξ(t). In that way the ensemble averages p n (t → ∞) 2 , necessary for both the temperature profiles and heat currents, are computed. This standard route implies a heavy computational effort, in particular when we want to study the heat transport for several bath configurations, frequency gradients and chain parameters. It is possible to circumvent this difficulty and find ensemble averages like x n x m , x n p m , p n p m (second order moments) without integrating any SDE [48]. The idea is to impose the condition d ··· dt = 0 for all the second order moments and linearize the dynamical equations of the system around equilibrium. A system of linear algebraic equations for the moments results, that can be easily solved without solving the SDE many times.
To linearize the SDE in Eq. (3) we approximate the potential energy of the Hamiltonian in Eq. (1), n ) 2 /2, by its harmonic approximation around the equilibrium positions x eq , defined by ∂V (x) ∂x x=x eq = 0. The approximate potential (ignoring the zero-point energy) is with K nm = ∂ 2 V (x) ∂xn∂xm x=x eq being the Hessian matrix entries of V (x) around the equilibrium configuration [49] (15) Note that this approximation does not modify the two main features of the system, namely asymmetry and long range interactions, which are manifest in the asymmetric distribution of ω n and the non-zero off-diagonal elements of the K matrix, respectively. In the following we will use y n = x n − x eq n to simplify the notation. The linearized dynamics around the equilibrium positions are given bẏ Now, we set d ··· /dt = 0 for all the moments. Using Eq. (16) and applying Novikov's theorem we get p n p l − γ l y n p l − m K lm y n y m = 0, m [K nm y m p l + K lm y m p n ] + 1 m (γ l + γ n ) p n p l = 2δ nl D n .
The system (17) is linear in the second order moments so it can be solved numerically to find the steady-state values of the moments. Besides Eq. (17) we have that y n p l = − y l p n , which follows from Eq. (16) and d y n y m /dt = 0. Since there are 1 2 N (N − 1) independent y n p l moments, we choose the ones with n < l. Similarly, the moments y n y l and p n p l contribute with

IV. NUMERICAL RESULTS
In this section we display the results of our simulations. To find the temperature profiles and the currents in the steady state we use the algebraic method described in section III A. We also check that the results coincide with solving Eq. (3) for many different realizations of the noisy forces ξ(t) and averaging. The code for all the numerical simulations have been written in the language Julia [50,51]. In particular, to solve the Langevin equation, we used Julia's package DifferentialEquations.jl [52]. To model the baths and the chain we use atomic data taken from ion trap experiments [53,54]. We consider 15 24 Mg + ions. The three leftmost and three rightmost ions are illuminated by Doppler cooling lasers. The Doppler cooling lasers excite the transition 3s 2 S 1/2 → 3p 2 P 1/2 , with angular frequency ω 0 = 2π × 1069 THz and excited state line width Γ = 2π × 41.3 MHz [32]. For this ionic species and atomic transition the Doppler limit is T D = 1 mK. The intensities of the laser beams are small compared to the saturation intensity I 0 so that Eq. (4) holds. We take I n /I 0 = 0.08 for the ions in the laser beams, whereas I n = 0 for the rest.
The temperatures T L , T R of the left and right laser baths are controlled with their detunings δ L , δ R with respect to the atomic transition. We fix two values for the detunings, δ H and δ C , such that T H > T C (hot and cold baths, also source and drain) and we define J → (J ← ) as the heat current in the chain when T L = T H and T R = T C (T L = T C and T R = T H ).
In this section we consider two frequency profiles for the traps: graded or segmented. In the graded chain the frequency increases by ∆ω N −1 from one trap to the next. If the frequency of the leftmost trap is ω 1 , the frequency of the nth trap will be ω n = ω 1 +∆ω n−1 N −1 . The frequency of the rightmost trap is ω 1 + ∆ω. In the segmented chain, the left half of the chain has trapping frequencies ω 1 while the other half has ω 1 + ∆ω. To compare the results by solving Eq. (3) and averaging and those from the algebraic method we simulated a frequency graded chain with a lattice constant (intertrap distance) a = 50 µm and a trapping frequency ω 1 = 2π × 50 kHz for the leftmost ion, Fig. 2. The number of ions interacting with the laser beams (three on each bath) is consistent with the lattice constant and typical waists of Gaussian laser beams [53,54]. To set the trap distance we fix first the characteristic length as the distance for which the Coulomb repulsion of two ions equals the trap potential energy for an ion at a distance l away from the center of its trap. If a < l, the Coulomb repulsion of the ions is stronger than the trap confinement which makes the ions jump from their traps. With the parameters used in this section we have l = 38.7 µm (a = 1.29 l). The detunings of the hot and cold lasers are δ H = −0.02 Γ and δ C = −0.1 Γ which gives temperatures T H ≈ 12 mK and T C ≈ 3 mK. We fix the value ∆ω = 0.5 ω 1 for the frequency gradient.
The results of the two methods are in very good agreement. In the scale of Fig. 2 (a) the calculated local temperatures are undistinguishable. In the calculation based on solving the dynamics we had to integrate (3) for N trials = 1000 realizations of white noise ξ(t). The method based on the system of moments shortened the calculation time with respect to the dynamical trajectories by a factor of 1/700. In fact, the time gain is even more important because the dynamical method requires further processing, performing a time averaging to compute the stationary flux in addition to noise averaging, see Fig. 2 (b).
Additionally, the relaxation to the steady state slows down when the frequencies of the traps increase since the deterministic part of the Langevin equation dominates the dynamics over the stochastic part, entering an underdamped regime. In contrast, this increase does not affect the algebraic method.

B. Rectification in frequency graded chains
In this section we display the results of some simulations that demonstrate rectification. We used the method described in section III A for 24 Mg + ions with the same parameters for the baths used before. We fix the trap- ping frequency of the leftmost trap to ω 1 = 2π × 1 MHz, which sets a characteristic length l = 5.25 µm, and a trap spacing a = 4.76 l (25 µm). Figure 3 shows the results with these parameters in a graded chain. Figure 3 (a) shows that both the right-going flux (dashed line) and the left-going flux (solid line) decrease rapidly as the frequency gradient is increased. The rectification reaches its maximum value for a frequency gradient of ∆ω ≈ 0.1ω 1 . The fluxes cross so there are some points where the rectification is exactly zero, besides the trivial one at ∆ω = 0, at ∆ω = 0.05 ω 1 , 0.3 ω 1 , 1.3 ω 1 . At these points the direction of rectification reverses, presumably as a consequence of the changes in the match/mismatch of the temperature dependent local power spectra. The change of rectification direction occurs for all the choices of parameters, as displayed in Fig. 4. Figure 4 gives the rectification factor for different trap distances and frequency gradients. 0-rectification curves separate regions with different rectification direction. The second region in Fig. 4 (starting from the left) would be the most interesting one to build a thermal diode, since rectification reaches its largest values there. We have also compared the performance of the graded thermal diode against a segmented version of it. Even though the optimal rectification in Fig. 5 (a) for the segmented chain is larger than for the graded chain, the fluxes are generally much larger for the graded chain, see Fig. 5 (b), which makes it more interesting for applications.

V. CONCLUSIONS
In this article we have numerically demonstrated heat rectification in a chain of ions trapped in individual microtraps with graded frequencies, connected at both ends to thermal baths created by optical molasses. An alternative to implement a graded frequency profile in the lab could be combining a collective Paul trap for all the ions with on-site dipolar laser forces [31,[55][56][57].
A goal of this article is to connect two communities, ion trappers and researchers on heat-rectification models. The results found are encouraging and demonstrate the potential of a trapped-ion platform to experimentally investigate heat rectification schemes. Trapped ions are interesting to this end because they are highly controllable, and may easily adopt several features to enhance rectification, such as the ones explored here (long-range interactions and an asymmetrical gradation), or others such as time dependent forces [5,58], or different nonlinearities in onsite forces.
The calculation of the steady state has been performed with an algebraic approach much faster than the timeconsuming integration and averaging over noise and time of the dynamical equations. The algebraic approach linearizes the forces around equilibrium positions which, in this system and for the realistic parameters considered is well justified and tested numerically. The results found provide additional evidence that simple linear models may rectify heat flow [18]. We underline that our linear model is, arguably, even simpler than some linear "minimalist, toy models" in [18] that showed rectification (our on-site forces are already linear from the start and the temperature dependence is only in the coefficients of the Langevin baths), with the important bonus of being also realistic.