Asymmetric Heat Transport in Ion Crystals

TL;DR

This paper demonstrates heat rectification in ion chains with graded trapping frequencies using a fast algebraic approach to find stationary states, advancing understanding of thermal transport in ion trap systems.

Contribution

It introduces a novel, efficient algebraic method to compute stationary states for heat transport in ion chains, enabling faster analysis of thermal rectification.

Findings

Heat rectification observed in ion chains with graded traps

Algebraic method significantly reduces computation time

Potential applications in thermal management and quantum devices

Abstract

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.

Equations

$$H(\mathbf{x},\mathbf{p})=\sum _{n=1}^N\left[\frac{p_n^2}{2m}+\frac{m{\omega }_n^2}{2}(x_n-x_n^{(0)}{)}^2\right]+V_{int}(\mathbf{x}),$$

$$V_{int}(\mathbf{x})=\frac{1}{2}\sum _n\sum _{l\ne n}{V}_C(\left|x_n-x_l\right|),$$

$$\begin{array}{rl}{\dot{x}}_n& =\frac{1}{m}{p}_n,\\ {\dot{p}}_n& =-m{\omega }_n^2(x_n-x_n^{(0)})-\frac{\partial V_{int}}{\partial x_n}-\frac{{\gamma }_n}{m}{p}_n+{\xi }_n(t),\end{array}$$

$$\begin{array}{rl}{\gamma }_n& =-4\hslash k_{L,R}^2\left(\frac{I_{L,R}}{I_0}\right)\frac{2{\delta }_{L,R}/\Gamma }{{\left[1+(2{\delta }_{L,R}/\Gamma {)}^2\right]}^2},\\ D_n& ={\hslash }^2{k}_{L,R}^2\left(\frac{I_{L,R}}{I_0}\right)\frac{\Gamma }{1+(2{\delta }_{L,R}/\Gamma {)}^2},\\ n& \in L,R,\end{array}$$

$$T_{L,R}=-\frac{\hslash \Gamma }{4k_B}\frac{1+(2{\delta }_{L,R}/\Gamma {)}^2}{(2{\delta }_{L,R}/\Gamma )},$$

$$H_n=\frac{1}{2m}{p}_n^2+\frac{1}{2}m{\omega }_n^2{\left(x_n-x_n^{(0)}\right)}^2+\frac{1}{2}\sum _{l\ne n}{V}_C(\left|x_n-x_l\right|).$$

$${\dot{H}}_n=\frac{p_n}{m}\left[{\xi }_n(t)-{\gamma }_n\frac{p_n}{m}\right]-\frac{1}{2m}\sum _{l\ne n}\frac{\partial V_C(\left|x_n-x_l\right|)}{\partial x_n}(p_n+p_l).$$

$$\text{expectationvalue}{j}_n^B+\text{expectationvalue}{\dot{H}}_n^{int}=0,$$

$$\begin{array}{rl}{J}_L& =\sum _{n\in L}\text{expectationvalue}{j}_n^B,\\ J_R& =\sum _{n\in R}\text{expectationvalue}{j}_n^B,\end{array}$$

$$\text{expectationvalue}\xi (t)\varphi (t)={\int }_0^{t}d{t}^{\prime }\text{expectationvalue}\xi (t)\xi (t^{\prime })\text{expectationvalue}\frac{\delta \varphi (t)}{\delta \xi (t^{\prime })},$$

$$\begin{array}{rl}{x}_n(t)& =x_n(0)+\frac{1}{m}{\int }_0^{t}ds{p}_n(s),\\ p_n(t)& =p_n(0)+{\int }_0^{t}ds\left[-\frac{\partial H}{\partial x_n}(s)-\frac{{\gamma }_n}{m}{p}_n(s)+{\xi }_n(s)\right].\end{array}$$

$$\text{expectationvalue}{j}_n^B=\frac{1}{m}\left[D_n-{\gamma }_n\frac{\text{expectationvalue}{p}_n^2}{m}\right].$$

$$R=\frac{J_{\to }-J_{\leftarrow }}{max(J_{\to },J_{\leftarrow })}.$$

$$V(\mathbf{x})\approx \frac{1}{2}\sum _{n,m}{K}_{nm}(x_n-x_n^{eq})(x_m-x_m^{eq}),$$

$$K_{nm}=\{\begin{array}{ll}m{\omega }_n^2+2\left(\frac{q^2}{4\pi {\epsilon }_0}\right){\sum }_{l\ne n}\frac{1}{{|x_n^{eq}-x_l^{eq}|}^3}& \text{if }n=m\\ -2\left(\frac{q^2}{4\pi {\epsilon }_0}\right)\frac{1}{{|x_n^{eq}-x_m^{eq}|}^3}& \text{if }n\ne m\end{array}.$$

$$\begin{array}{rl}{\dot{y}}_n& =\frac{1}{m}{p}_n,\\ {\dot{p}}_n& =-\sum _l{K}_{nl}{y}_l-\frac{{\gamma }_n}{m}{p}_n+{\xi }_n(t).\end{array}$$

$$\begin{array}{rl}\text{expectationvalue}{p}_n{p}_l-{\gamma }_l\text{expectationvalue}{y}_n{p}_l-\sum _m{K}_{lm}\text{expectationvalue}{y}_n{y}_m& =0,\\ \sum _m\left[K_{nm}\text{expectationvalue}{y}_m{p}_l+K_{lm}\text{expectationvalue}{y}_m{p}_n\right]+\frac{1}{m}\left({\gamma }_l+{\gamma }_n\right)\text{expectationvalue}{p}_n{p}_l& =2{\delta }_{nl}{D}_n.\end{array}$$

$$\begin{split}\bm{\eta}=\Big{[}&\expectationvalue{y_{1}y_{1}},\expectationvalue{y_{1}y_{2}},...,\expectationvalue{y_{N}y_{N}},\\ &\expectationvalue{p_{1}p_{1}},\expectationvalue{p_{1}p_{2}},...,\expectationvalue{p_{N}p_{N}},\\ &\expectationvalue{y_{1}p_{2}},\expectationvalue{y_{1}p_{3}},...,\expectationvalue{y_{N-1}p_{N}}\Big{]}^{T}\end{split}$$

Full text

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 Starr (1936). 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}$ Terraneo et al. (2002). This paper sparked a substantial body of research to this day Pereira (2019) (see Fig. 1 in Roberts and Walker (2011)).

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 Roberts and Walker (2011); Li et al. (2012). There are fascinating proposals to engineer thermal logic circuits Ye and Cao (2017) in which information, stored in thermal memories Wang and Li (2008), would be processed in thermal gates Wang and Li (2007). Such thermal gates, as their electronic counterparts, will require thermal diodes and thermal transistors to operate Li et al. (2006); Joulain et al. (2016). Heat rectifying devices would also be quite useful in nano electronic circuits, letting delicate components dissipate heat while being protected from external heat sources Roberts and Walker (2011).

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 Chang et al. (2006). 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 Terraneo et al. (2002), on a chain segmented into two or more regions with different properties, but using other lattice models such as the Frenkel-Kontorova (FK) model Li et al. (2008); Hu et al. (2006). The fundamental ingredient for having rectification was attributed to nonlinear forces in the chain Zeng and Wang (2008); Katz and Kosloff (2016); Li et al. (2008); Hu et al. (2006); Benenti et al. (2016); Li et al. (2012), 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 Terraneo et al. (2002); Li et al. (2004). Later, alternative mechanisms have been also proposed which do not necessarily rely on anharmonic potentials Pereira (2017); Pons et al. (2017).

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 Hu et al. (2006). 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 Wang et al. (2012); Chen et al. (2015); Romero-Bastida et al. (2017); Yang et al. (2007); Romero-Bastida and Arizmendi-Carvajal (2013); Dettori et al. (2016); Pereira (2010, 2011); Ávila and Pereira (2013). The second one uses particles with long range interactions (LRI), such that all the particles in the lattice interact with each other Chen et al. (2015); Bagchi (2017); Pereira and Ávila (2013). 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 Chen et al. (2015). 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 Freitas et al. (2016); Ruiz et al. (2014); Pruttivarasin et al. (2011); Ramm et al. (2014) and interesting phenomena like phase transitions have been investigated Freitas et al. (2016); Ruiz et al. (2014); Pruttivarasin et al. (2011). The idea of using locally-controlled traps is already mentioned in Freitas et al. (2016) 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 Cirac and Zoller (2000); Krauth et al. (2014); Schmied et al. (2009). 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 $\omega_{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)

[TABLE]

where $\{x_{n},p_{n}\}$, position and momentum of each ion, are the components of the vectors $\bm{x},\bm{p}$, $x_{n}^{(0)}=na$ are the centers of the harmonic traps, and $V_{int}$ is the sum of the Coulomb interaction potential between all pairs of ions,

[TABLE]

with $V_{C}(\left|x_{n}-x_{l}\right|)=\frac{q^{2}}{4\pi\varepsilon_{0}}\frac{1}{\left|x_{n}-x_{l}\right|}$. 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}$ Lepri et al. (2003); Dhar (2008). The equations of motion of the chain, taking into account the baths and the Hamiltonian, are

[TABLE]

where $\gamma_{n}$ and $\xi_{n}(t)$ are only non-zero for the ions in the end regions, in contact with the left and right baths in the sets $L=\left\{1,2,...,N_{L}\right\}$ and $R=\left\{N-(N_{R}-1),...,N-1,N\right\}$, see Fig. 1. The $\gamma_{n}$ are friction coefficients and $\xi_{n}(t)$ are uncorrelated Gaussian noise forces satisfying $\expectationvalue{\xi_{n}(t)}=0$ and $\expectationvalue{\xi_{n}(t)\xi_{m}(t^{\prime})}=2D_{n}\delta_{nm}\delta(t-t^{\prime})$, $D_{n}$ being the diffusion coefficients. These Gaussian forces are formally the time derivatives of independent Wiener processes (Brownian motions) $\xi_{n}(t)=\sqrt{2D_{n}}\frac{dW_{n}}{dt}$ Toral and Colet (2014); Ruiz et al. (2014) and Eq. (3) is a stochastic differential equation (SDE) in the Stratonovich sense Toral and Colet (2014).

The baths are physically implemented by optical molasses consisting of a pair of counterpropagating Doppler-cooling lasers Ruiz et al. (2014). The friction and diffusion coefficients for the ions in contact with the baths are given by Arimondo et al. (1993)

[TABLE]

where $k_{L}$ ($k_{R}$) and $I_{L}$ ($I_{R}$) are the wave vector and intensity of the left (right) laser. $\delta_{L}$ ($\delta_{R}$) is the detuning of the left (right) laser with respect to the angular frequency $\omega_{0}$ of the atomic transition the laser is exciting, and $\Gamma$ 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}\ll 1$. In this bath model, the friction term in Eq. (3) comes from the cooling action of the laser and the white noise force $\xi_{n}(t)$ corresponds to the random recoil of the ions due to spontaneous emission of photons Metcalf and Van der Straten (1999); Metcalf and van der Straten (2003); Cohen-Tannoudji (1990). Using the diffusion-dissipation relation $D=\gamma k_{B}T$ Chee Kong (2010), the temperature of the optical molasses baths are given by

[TABLE]

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 $\delta=-\Gamma/2$ the optical molasses reach their minimum temperature possible, the Doppler limit $T_{D}={\hbar\Gamma}/({2k_{B}})$.

III Calculation of the stationary heat currents

The local energy of each site is defined by

[TABLE]

Differentiating $H_{n}$ with respect to time we find the continuity equation

[TABLE]

Two different contributions can be distinguished: $j^{B}_{n}\equiv\frac{p_{n}}{m}\left[\xi_{n}(t)-\gamma_{n}\frac{p_{n}}{m}\right]$, which is the energy flow from the laser reservoir to the ions at the edges of the chain (only for $n\in L,R$), and $\dot{H}_{n}^{int}\equiv-\frac{1}{2m}\sum_{l\neq n}\frac{\partial V_{C}(\left|x_{n}-x_{l}\right|)}{\partial x_{n}}(p_{n}+p_{l})$, which gives the “internal” energy flow due to the interactions with the rest of the ions. In the steady state $\expectationvalue{\dot{H}_{n}}=0$, and therefore

[TABLE]

where $\langle\cdot\!\cdot\!\cdot\rangle$ stands for the expectation value with respect to the ensemble of noise processes $\bm{\xi}(t)$ ($\bm{\xi}$ represents a vector with components $\xi_{n}$). Equation (8) implies that, in the steady state, the internal rates $\dot{H}_{n}^{int}$ vanish for the inner ions of the chain because $j^{B}_{n}=0$ for $n\notin L,R$. In chains with nearest-neighbor (NN) interactions, $\langle\dot{H}_{n}^{int}\rangle$ 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 Chen et al. (2015), but this approach is not realistic in the current system so we define instead the heat currents for the left and right baths as

[TABLE]

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\equiv\absolutevalue{J_{L}}=\absolutevalue{J_{R}}$. $J$ is defined as $J_{\rightarrow}$ when the hot bath is on the left and $J_{\leftarrow}$ when it is on the right.

To compute the average heat fluxes of the baths $\expectationvalue{j^{B}_{n}}$ in Eq. (9) we need the averages $\expectationvalue{p_{n}(t)\xi_{n}(t)}$. Instead of explicitly averaging $p_{n}(t)\xi_{n}(t)$ over different realizations of the white noise we use Novikov’s theorem Novikov (1965); Ma and Dudarev (2011); Toral and Colet (2014). Novikov’s theorem states that the ensemble average (over the realizations of the noise) of the product of some functional $\phi(t)$, which depends on a Gaussian noise $\xi(t)$ with zero mean value, $\expectationvalue{\xi(t)}=0$, and the noise itself, is given by

[TABLE]

where ${\delta\phi(t)}/{\delta\xi(t^{\prime})}$ is the functional derivative of $\phi(t)$ with respect to the noise, with $t^{\prime}<t$. When the noise is $\delta-$correlated, $\expectationvalue{\xi(t)\xi(t^{\prime})}=2D\delta(t-t^{\prime})$, and Eq. (10) reads $\expectationvalue{\xi(t)\phi(t)}=D\expectationvalue{{\delta\phi(t)}/{\delta\xi(t^{\prime})}}|_{t^{\prime}\to t^{-}}$. 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 $\xi_{n}(t)$,

[TABLE]

Equation (11) implies that the functional derivatives are ${\delta x_{n}(t)}/{\delta\xi_{m}(t^{\prime})}|_{t^{\prime}\to t^{-}}=0$ and ${\delta p_{n}(t)}/{\delta\xi_{m}(t^{\prime})}|_{t^{\prime}\to t^{-}}=\delta_{nm}$ ($\delta_{nm}$ is the usual Kronecker delta symbol). Thus we have $\expectationvalue{x_{n}(t)\xi_{m}(t)}=0$ and $\expectationvalue{p_{n}(t)\xi_{m}(t)}=\delta_{nm}D_{m}$, which gives for the heat flow from the baths

[TABLE]

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

[TABLE]

$R$ values may go from $-1$ to $1$. If there is no rectification $J_{\rightarrow}=J_{\leftarrow}$ and $R=0$. For perfect rectification in the right (left) direction, $J_{\rightarrow}\gg J_{\leftarrow}$ ($J_{\rightarrow}\ll J_{\leftarrow}$), and $R=1$ ($R=-1$).

III.1 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 $\bm{\xi}(t)$. In that way the ensemble averages $\expectationvalue{p_{n}(t\to\infty)^{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 $\expectationvalue{x_{n}x_{m}}$, $\expectationvalue{x_{n}p_{m}}$, $\expectationvalue{p_{n}p_{m}}$ (second order moments) without integrating any SDE Särkkä and Solin (2019). The idea is to impose the condition $\frac{d\expectationvalue{\cdot\cdot\cdot}}{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), $V(\bm{x})=V_{int}(\bm{x})+m\;\sum_{n}\omega_{n}^{2}(x_{n}-x_{n}^{(0)})^{2}/2$, by its harmonic approximation around the equilibrium positions $\bm{x}^{eq}$, defined by \frac{\partial V(\bm{x})}{\partial\bm{x}}\Big{|}_{\bm{x}=\bm{x}^{eq}}=0. The approximate potential (ignoring the zero-point energy) is

[TABLE]

with K_{nm}=\frac{\partial^{2}V(\bm{x})}{\partial x_{n}\partial x_{m}}\Big{|}_{\bm{x}=\bm{x}^{eq}} being the Hessian matrix entries of $V(\bm{x})$ around the equilibrium configuration James (1998)

[TABLE]

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 $\omega_{n}$ and the non-zero off-diagonal elements of the $K$ matrix, respectively. In the following we will use $y_{n}=x_{n}-x_{n}^{eq}$ to simplify the notation. The linearized dynamics around the equilibrium positions are given by

[TABLE]

Now, we set ${d\expectationvalue{\cdot\!\cdot\!\cdot}}/{dt}=0$ for all the moments. Using Eq. (16) and applying Novikov’s theorem we get

[TABLE]

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 $\expectationvalue{y_{n}p_{l}}=-\expectationvalue{y_{l}p_{n}}$, which follows from Eq. (16) and $d\expectationvalue{y_{n}y_{m}}/dt=0$. Since there are $\frac{1}{2}N(N-1)$ independent $\expectationvalue{y_{n}p_{l}}$ moments, we choose the ones with $n<l$. Similarly, the moments $\expectationvalue{y_{n}y_{l}}$ and $\expectationvalue{p_{n}p_{l}}$ contribute with $\frac{1}{2}N(N+1)$ independent variables each and we choose the ones with $n\leq m$. Thus there are in total $\frac{1}{2}N(3N+1)$ independent moments that we arrange in the vector

[TABLE]

There are the same number of independent equations as independent moments: $N^{2}$ equations correspond to the first line in Eq. (17), and $\frac{1}{2}N(N+1)$ equations to the second line because of the symmetry with respect to $n,l$. The system of equations (17) may be compactly written as $\mathbf{A}\bm{\eta}=\mathbf{B}$, where $\mathbf{A}$ and $\mathbf{B}$ are a $\frac{1}{2}N(3N+1)$ square matrix and vector.

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.1. We also check that the results coincide with solving Eq. (3) for many different realizations of the noisy forces $\bm{\xi}(t)$ and averaging. The code for all the numerical simulations have been written in the language Julia Bezanson et al. (2012, 2017). In particular, to solve the Langevin equation, we used Julia’s package DifferentialEquations.jl Rackauckas and Nie (2017).

To model the baths and the chain we use atomic data taken from ion trap experiments Leupold (2015); Lo (2015). We consider 15 24Mg+ 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}\rightarrow 3p^{2}P_{1/2}$, with angular frequency $\omega_{0}=2\pi\times 1069$ THz and excited state line width $\Gamma=2\pi\times 41.3$ MHz Ruiz et al. (2014). 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 $\delta_{L},\,\delta_{R}$ with respect to the atomic transition. We fix two values for the detunings, $\delta_{H}$ and $\delta_{C}$, such that $T_{H}>T_{C}$ (hot and cold baths, also source and drain) and we define $J_{\rightarrow}$ ($J_{\leftarrow}$) 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 $\frac{\Delta\omega}{N-1}$ from one trap to the next. If the frequency of the leftmost trap is $\omega_{1}$, the frequency of the $n$th trap will be $\omega_{n}=\omega_{1}+\Delta\omega\frac{n-1}{N-1}$. The frequency of the rightmost trap is $\omega_{1}+\Delta\omega$. In the segmented chain, the left half of the chain has trapping frequencies $\omega_{1}$ while the other half has $\omega_{1}+\Delta\omega$.

IV.1 Evolution to steady state

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\,\mu$m and a trapping frequency $\omega_{1}=2\pi\times 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 Leupold (2015); Lo (2015). To set the trap distance we fix first the characteristic length $l=\left(\frac{q^{2}}{4\pi\varepsilon_{0}}\frac{1}{m\omega_{1}^{2}}\right)^{1/3}$ 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\,\mu$m ($a=1.29\,l$). The detunings of the hot and cold lasers are $\delta_{H}=-0.02\,\Gamma$ and $\delta_{C}=-0.1\,\Gamma$ which gives temperatures $T_{H}\approx 12$ mK and $T_{C}\approx 3$ mK. We fix the value $\Delta\omega=0.5\,\omega_{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 $\bm{\xi}(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 under-damped regime. In contrast, this increase does not affect the algebraic method.

IV.2 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.1 for 24Mg+ ions with the same parameters for the baths used before. We fix the trapping frequency of the leftmost trap to $\omega_{1}=2\pi\times 1$ MHz, which sets a characteristic length $l=5.25\,\mu$m, and a trap spacing $a=4.76\,l$ ($25\,\mu$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 $\Delta\omega\approx 0.1\omega_{1}$. The fluxes cross so there are some points where the rectification is exactly zero, besides the trivial one at $\Delta\omega=0$, at $\Delta\omega=0.05\,\omega_{1},\;0.3\,\omega_{1},\;1.3\,\omega_{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. [math]-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 Freitas et al. (2016); Enderlein et al. (2012); Bermudez et al. (2013); Schneider et al. (2010).

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 Li et al. (2012); Riera-Campeny et al. (2019), or different nonlinearities in onsite forces.

The calculation of the steady state has been performed with an algebraic approach much faster than the time-consuming 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 Pereira (2017). We underline that our linear model is, arguably, even simpler than some linear “minimalist, toy models” in Pereira (2017) 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.

Acknowledgements.

The authors would like to thank Joseba Alonso and Daniel Alonso for fruitful discussions and comments. M. A. Simón acknowledges support by the Basque Government predoctoral program (Grant no. PRE-2018-2-0177). We thank funding from Grant No. PGC2018-101355-B-I00 and Fis2016-80681-P (MCIU/AEI/FEDER,UE), and Basque Government (Grant No. IT986- 16).