Delocalization and heat transport in multidimensional trapped ion systems

TL;DR

This paper investigates how the structural phase transitions in trapped ion systems influence heat transport, revealing that delocalization and spatial configurations significantly affect thermal conductivity in multidimensional Coulomb crystals.

Contribution

It demonstrates, through numerical simulations, how structural phase transitions in trapped ion systems can be used to tune heat transport properties.

Findings

Heat transport is affected by the system's structural phase.

Delocalization correlates with temperature gradient emergence.

Optimal configurations enhance heat flux.

Abstract

We study the connection between heat transport properties of systems coupled to different thermal baths in two separate regions and their underlying nonequilibrium dynamics. We consider classical systems of interacting particles that may exhibit a certain degree of delocalization, and whose effective dimensionality can be modified through the controlled variation of a global trapping potential. We focus on Coulomb crystals of trapped ions, which offer a versatile playground to shed light on the role that spatial constraints play on heat transport. We use a three-dimensional model to simulate the trapped ion system, and show in a numerically rigorous manner to what extent heat transport properties could be feasibly tuned across the structural phase transitions between the linear, planar zigzag and helical configurations. By solving the classical Langevin equations of motion, we analyze the steady state spatial distributions of the particles, the temperature profiles and total heat flux through the various structural phase transitions that the system may experience. The results evidence a clear correlation between the degree of delocalization of the internal ions and the emergence of a non-zero gradient in the temperature profiles. The signatures of the phase transitions in the total heat flux as well as the optimal spatial configuration for heat transport are shown.

Equations

$$H=\sum _{i=1}^N\left[\frac{{\mathbf{p}}_i^2}{2m}+\frac{1}{2}\sum _{j\ne i}^N{U}_{ij}(|{\mathbf{q}}_i-{\mathbf{q}}_j|)+V_i({\mathbf{q}}_i)\right].$$

$${\dot{q}}_{\alpha ,i}=\frac{p_{\alpha ,i}}{m}\text{mbox}fori=(1,\dots ,N),$$

$${\dot{p}}_{\alpha ,i}=g_{\alpha ,i}+\sum _{j\ne i}^N{f}_{\alpha }^{(ij)}-\frac{{\eta }_{\alpha ,i}^L}{m}{p}_{\alpha ,i}+{\epsilon }_{\alpha ,i}^L(t)\text{mbox}fori=(1,\dots ,N_L),$$

$${\dot{p}}_{\alpha ,i}=g_{\alpha ,i}+\sum _{j\ne i}^N{f}_{\alpha }^{(ij)}\text{mbox}fori=(N_L+1,\dots ,N-N_R),$$

$${\dot{p}}_{\alpha ,i}=g_{\alpha ,i}+\sum _{j\ne i}^N{f}_{\alpha }^{(ij)}-\frac{{\eta }_{\alpha ,i}^R}{m}{p}_{\alpha ,i}+{\epsilon }_{\alpha ,i}^R(t)\text{mbox}fori=(N-N_R+1,\dots ,N),$$

$$⟨{\epsilon }_{\alpha ,i}^{L,R}(t)⟩=0,$$

$$⟨{\epsilon }_{\alpha ,i}^{L,R}(t){\epsilon }_{\beta ,j}^{L,R}(t^{\prime })⟩=2D_{\alpha ,i}^{L,R}{\delta }_{\alpha ,\beta }{\delta }_{i,j}\delta (t-t^{\prime }),$$

$${\eta }_{\alpha ,i}^{L,R}=\frac{1}{2k_B{T}^{L,R}}{\int }_{-\infty }^{\infty }⟨{\epsilon }_{\alpha ,i}^{L,R}(t){\epsilon }_{\alpha ,i}^{L,R}(t+\tau )⟩d\tau =\frac{D_{\alpha ,i}^{L,R}}{k_B{T}^{L,R}},$$

$$d{q}_{\alpha ,i}=\frac{p_{\alpha ,i}}{m}dt\text{mbox}fori=(1,\dots ,N),$$

$$d{p}_{\alpha ,i}=\left(g_{\alpha ,i}+\sum _{j\ne i}^N{f}_{\alpha }^{(ij)}-\frac{{\eta }_{\alpha ,i}^L}{m}{p}_{\alpha ,i}\right)dt+\sqrt{2D_{\alpha ,i}^L}d{W}_{\alpha ,i}^L\text{mbox}fori=(1,\dots ,N_L),$$

$$d{p}_{\alpha ,i}=\left(g_{\alpha ,i}+\sum _{j\ne i}^N{f}_{\alpha }^{(ij)}\right)dt\text{mbox}fori=(N_L+1,\dots ,N-N_R),$$

$$d{p}_{\alpha ,i}=\left(g_{\alpha ,i}+\sum _{j\ne i}^N{f}_{\alpha }^{(ij)}-\frac{{\eta }_{\alpha ,i}^R}{m}{p}_{\alpha ,i}\right)dt+\sqrt{2D_{\alpha ,i}^R}d{W}_{\alpha ,i}^R\text{mbox}fori=(N-N_R+1,\dots ,N),$$

$$P(\mathbf{q})=\frac{1}{{\tau }_{ss}}{\int }_t^{t+{\tau }_{ss}}d\tau \left[\frac{1}{N{\sigma }^3(2\pi {)}^{3/2}}\sum _{i=1}^N{e}^{-{\left[\mathbf{q}-{\mathbf{q}}_i(\tau )\right]}^2/2{\sigma }^2}\right],$$

$$\Theta_{i}(\alpha)\,\equiv\,\Theta_{i}(\alpha_{l})\,=\,\frac{1}{\tau_{ss}}\,\bigg{\langle}\,{\displaystyle\int_{t}^{t+\tau_{ss}}}d\tau\,{\displaystyle\int_{\alpha_{l}-\Delta/2}^{\alpha_{l}+\Delta/2}}\,d\alpha\,\delta\left(\,\alpha\,-\,q_{\alpha,i}(\tau)\,\right)\,\bigg{\rangle}\,.$$

$$T(\alpha)\,\equiv\,T(\alpha_{l})=\frac{2}{3\,k_{B}}\,\bigg{\langle}\,\frac{\sum\limits_{i=1}^{N}{\displaystyle\int_{t}^{t+\tau_{ss}}}d\tau\,{\displaystyle\int_{\alpha_{l}-\Delta/2}^{\alpha_{l}+\Delta/2}}\,d\alpha\,\delta\left(\,\alpha\,-\,q_{\alpha,i}(\tau)\,\right)E_{k}({\bf p}_{i}(\tau))}{\sum\limits_{i=1}^{N}{\displaystyle\int_{t}^{t+\tau_{ss}}}d\tau\,{\displaystyle\int_{\alpha_{l}-\Delta/2}^{\alpha_{l}+\Delta/2}}\,d\alpha\,\delta\left(\alpha\,-\,q_{\alpha,i}(\tau)\,\right)}\,\bigg{\rangle}\,,$$

$$h(\mathbf{q})=\sum _{i=1}^N{h}_i\delta (\mathbf{q}-{\mathbf{q}}_i),$$

$$h_i=\frac{{\mathbf{p}}_i^2}{2m}+V_i({\mathbf{q}}_i)+\frac{1}{2}\sum _{j\ne i}^N{U}_{ij}(|{\mathbf{q}}_i-{\mathbf{q}}_j|).$$

$$\frac{\partial h(\mathbf{q})}{\partial t}=\sum _{i=1}^N\left[\frac{1}{2}\sum _{j\ne i}^N\frac{1}{m}\left({\mathbf{p}}_i+{\mathbf{p}}_j\right)\cdot {\mathbf{f}}^{(ij)}+J_i^B\right]\delta (\mathbf{q}-{\mathbf{q}}_i)-\nabla \cdot \sum _{i=1}^N\frac{{\mathbf{p}}_i}{m}{h}_i\delta (\mathbf{q}-{\mathbf{q}}_i),$$

$$\displaystyle J^{B}_{i}=\left\{\begin{array}[]{cl}\sum\limits_{\alpha}\frac{{p}_{\alpha,i}}{m}\left(-\frac{\eta_{\alpha,i}^{L}}{m}\,{p}_{\alpha,i}+\varepsilon_{\alpha,i}^{L}(t)\right)&\mbox{for}\hskip 8.5359pti=(1,\dots N_{L})\,,\\ 0&\mbox{for}\hskip 8.5359pti=(N_{L}+1,\dots N-N_{R})\,,\\ \sum\limits_{\alpha}\frac{{p}_{\alpha,i}}{m}\left(-\frac{\eta_{\alpha,i}^{R}}{m}\,{p}_{\alpha,i}+\varepsilon_{\alpha,i}^{R}(t)\right)&\mbox{for}\hskip 8.5359pti=(N-N_{R}+1,\dots N)\,,\\ \end{array}\right.$$

$$\frac{\partial h(\mathbf{q})}{\partial t}+\nabla \cdot \sum _{i=1}^N\frac{{\mathbf{p}}_i}{m}{h}_i\delta (\mathbf{q}-{\mathbf{q}}_i)-\frac{1}{4}\sum _{i=1}^N\sum _{j\ne i}^N\frac{1}{m}\left({\mathbf{p}}_i+{\mathbf{p}}_j\right)\cdot {\mathbf{f}}^{(ij)}\left[\delta (\mathbf{q}-{\mathbf{q}}_i)-\delta (\mathbf{q}-{\mathbf{q}}_j)\right]=\sum _{i=1}^N{J}_i^B\delta (\mathbf{q}-{\mathbf{q}}_i),$$

$$\delta (\mathbf{q}-{\mathbf{q}}_j)-\delta (\mathbf{q}-{\mathbf{q}}_i)=\nabla \cdot {\int }_0^{1}d\lambda \frac{d\mathbf{G}(\lambda )}{d\lambda }\delta (\mathbf{q}-\mathbf{G}(\lambda )),$$

$$\frac{\partial h(\mathbf{q})}{\partial t}+\nabla \cdot {\mathbf{j}}_h(\mathbf{q})={\sigma }_h(\mathbf{q}),$$

$${\mathbf{j}}_h(\mathbf{q})=\sum _{i=1}^N\frac{{\mathbf{p}}_i}{m}{h}_i\delta (\mathbf{q}-{\mathbf{q}}_i)+\frac{1}{4}\sum _{i=1}^N\sum _{j\ne i}^N\frac{1}{m}\left({\mathbf{p}}_i+{\mathbf{p}}_j\right)\cdot {\mathbf{f}}^{(ij)}{\int }_0^{1}d\lambda \frac{d\mathbf{G}(\lambda )}{d\lambda }\delta (\mathbf{q}-\mathbf{G}(\lambda )),$$

$${\sigma }_h(\mathbf{q})=\sum _{i=1}^N{J}_i^B\delta (\mathbf{q}-{\mathbf{q}}_i).$$

$$\mathbf{J}(t)=\sum _{i=1}^N\frac{{\mathbf{p}}_i}{m}{h}_i+\frac{1}{4}\sum _{i=1}^N\sum _{j\ne i}^N\frac{1}{m}\left({\mathbf{p}}_i+{\mathbf{p}}_j\right)\cdot {\mathbf{f}}^{(ij)}({\mathbf{q}}_i-{\mathbf{q}}_j).$$

$$\big{\langle}\,\frac{{\bf p}_{i}}{m}\,h_{i}\,\big{\rangle}_{ss}\,=\,-\big{\langle}\,{\bf q}_{i}\,\frac{dh_{i}}{dt}\,\big{\rangle}_{ss}\,=\,-\frac{1}{2}\sum_{i=1}^{N}\sum_{j\neq i}^{N}\big{\langle}\frac{1}{m}\left({\bf p}_{i}+{\bf p}_{j}\right)\cdot{\bf f}^{(ij)}\,{\bf q}_{i}\,\big{\rangle}_{ss}\,-\,\sum_{i=1}^{N}\big{\langle}\,J^{B}_{i}\,{\bf q}_{i}\,\big{\rangle}_{ss}\,,$$

$$⟨\mathbf{J}{⟩}_{ss}=-\sum _{i=1}^N⟨J_i^B{\mathbf{q}}_i{⟩}_{ss}=\frac{1}{m^2}\sum _i{}^{\prime }\sum _{\alpha }{\eta }_{\alpha ,i}^{B_i}⟨{\mathbf{q}}_i{p}_{\alpha ,i}^2{⟩}_{ss}-\frac{1}{m}\sum _i{}^{\prime }\sum _{\alpha }⟨{\mathbf{q}}_i{p}_{\alpha ,i}{\epsilon }_{\alpha ,i}^{B_i}(t){⟩}_{ss}.$$

$$\langle\,q_{\eta,i}\,p_{\alpha,i}\,\varepsilon_{\alpha,i}^{B_{i}}(t)\,\rangle_{ss}\,=\,\sum_{\nu}\,\int\,dt^{\prime}\,\langle\,\varepsilon_{\alpha,i}^{B_{i}}(t)\,\varepsilon_{\nu,i}^{B_{i}}(t^{\prime})\,\rangle\,\bigg{\langle}\,\frac{\partial(q_{\eta,i}\,p_{\alpha,i})}{\partial\varepsilon_{\nu,i}^{B_{i}}}\,\,\bigg{\rangle}_{ss}\,=\,D_{\alpha,i}^{B_{i}}\,q_{\eta,i}\,,$$

$$⟨\mathbf{J}{⟩}_{ss}=\frac{1}{m^2}\sum _i{}^{\prime }\sum _{\alpha }\left[{\eta }_{\alpha ,i}^{B_i}⟨{\mathbf{q}}_i{p}_{\alpha ,i}^2{⟩}_{ss}-m{D}_{\alpha ,i}^{B_i}⟨{\mathbf{q}}_i{⟩}_{ss}\right].$$

$$V_i({\mathbf{q}}_i)=\frac{1}{2}m{\omega }_x^2(q_{x,i}^2+n_y^2{q}_{y,i}^2+n_z^2{q}_{z,i}^2),$$

$$U_{ij}(|{\mathbf{q}}_i-{\mathbf{q}}_j|)=\left(\frac{Q^2}{4\pi {ϵ}_0}\right)\frac{1}{|{\mathbf{q}}_i-{\mathbf{q}}_j|},$$

$${\eta }_{\alpha ,i}^{L(R)}=-4\hslash {\left(k_{\alpha ,i}^{L(R)}\right)}^2\left(\frac{I_{\alpha ,i}^{L(R)}}{I_0^{L(R)}}\right)\frac{\left(2{\delta }_{\alpha ,i}^{L(R)}/\Gamma \right)}{{\left[1+4{\left({\delta }_{\alpha ,i}^{L(R)}\right)}^2/{\Gamma }^2\right]}^2}$$

$$D_{\alpha ,i}^{L(R)}={\hslash }^2{\left(k_{\alpha ,i}^{L(R)}\right)}^2\left(\frac{I_{\alpha ,i}^{L(R)}}{I_0^{L(R)}}\right)\frac{\Gamma }{\left[1+4{\left({\delta }_{\alpha ,i}^{L(R)}\right)}^2/{\Gamma }^2\right]}.$$

Full text

I Introduction

The downsizing of electronic devices to the nanometric scale, driven by the rapid progress of micro-electronic technology, has made the problem of thermal conduction increasingly important because of the need to find ways to dissipate a significant amount of energy in a shrinking compact space dubi11 . In this sense, the divergence of thermal conductivity with the size of the materials of reduced dimensionality would allow the construction of nano-materials capable of dissipating heat efficiently. This would solve one of the fundamental problems arising from the miniaturization of electronic and optical devices.

According to Fourier’s law, given a system connected to different heat reservoirs in two separate regions, the amount of heat transferred per unit area and time unit has a linear dependence on the imposed temperature gradient, the thermal conductivity that characterizes the material being the constant of proportionality between both magnitudes. Although this law can be verified in a simple way for three-dimensional systems, it is known that heat conduction exhibits anomalous behavior in systems of reduced dimensionality, such as carbon nanotubes, silicon nanowires or molecular junctions dubi11 ; sevik11 ; lee17 . This issue keeps many important and fundamental questions open in the field of nonequilibrium statistical physics lepri03 ; dhar08 ; gaspard08 .

From a theoretical perspective, the study of heat transport through nanoscale devices typically requires making tradeoffs between the size of the system and the completeness of the model. Most of the work on low-dimensional systems have considered simple, yet nontrivial, models that incorporate elements crucial to heat conduction, such as anharmonicity and disorder li01 ; dhar08b ; lepri03 ; dhar08 . It has been shown that some one-dimensional systems, such as the Frenkel-Kontorova model hu98 or the Lorentz model alonso99 , the temperature gradient is uniform and the thermal conductivity is a constant, independent of the size of the system. This indicates that these systems obey Fourier’s law. In contrast, in the case of one-dimensional integrable systems, such as a harmonic chain and the Toda mono-atomic model, a temperature gradient is not established rieder67 ; dhar16 ; toda79 ; kundu16 . There are also one-dimensional nonintegrable systems, such as the Fermi-Pasta-Ulam model kaburaki93 ; lepri97 ; lepri98 or the Toda diatomic chain hatano99 , for which the thermal conductivity diverges with increasing system size. In two dimensions, anomalous conductivities exhibiting a logarithmic divergence with the size of system have been reported lepri03 ; dhar08 ; yang06 . In the case of polygonal channels with zero Lyapunov exponents, numerical simulations have shown transport properties ranging from normal to anomalous conduction, depending on the system parameters li02 ; alonso02 ; alonso04 . Although these mathematical models have shed some light on the underlying mechanisms for normal heat conduction, understanding heat transport at the microscopic level remains a central topic of current research.

Coulomb crystals of ions confined in electromagnetic traps and manipulated with laser beams ghosh95 provide a versatile platform to study a broad range of intriguing physical phenomena emerging in system driven out of equilibrium, in particular energy transport in both classical and quantum regimes pruttivarasin11 ; lin11 ; bermudez13 ; ramm14 ; ruiz14 ; freitas15 ; cormick16 ; blatt12 ; bermudez16 ; cosco17 . Due to the unique control in the preparation, manipulation, and detection of the electronic and motional degrees of freedom, they have also become a promising candidate for quantum information and computation, quantum networks, quantum simulations, and quantum metrology haffner08 ; schneider12 ; pezze18 . From the theoretical point of view, Coulomb crystals are particularly interesting since the typical distances between ions of several micrometers, with densities several orders of magnitude smaller than the standard crystalline materials, simplifies to a great extent the analysis of many issues related to their very rich and highly non-trivial static and dynamical properties.

It has been shown that the thermodynamic behavior of trapped ion systems strongly depends on the crystal structure, resulting from the interplay between Coulomb repulsion and the trapping potential, and which can be experimentally controlled to an exquisite degree birkl92 ; waki92 ; dubin99 ; mavadia13 ; yan16 . For a highly anisotropic trapping potential the ions exhibit an inhomogeneous alignment in the axial direction ghosh95 . Due to the approximately harmonic confinement the center of the chain presents a higher density, and therefore a higher Coulomb repulsion. The decrease of the radial anisotropy can trigger structural phase transitions that start at the center of the chain and extend towards the edges as the anisotropy decreases birkl92 ; waki92 ; dubin99 ; mavadia13 ; yan16 . In particular, a second order structural phase transition from the linear chain to a planar zigzag spatial configuration has been well characterized schiffer93 ; piacente04 ; fishman08 . Also, the formation of planar concentric ellipses yan16 and three-dimensional helical configurations waki92 ; hasse90 ; nigmatullin16 has been reported.

Theoretical studies of heat transport in crystals of trapped ions connected to heat baths at different temperatures in two separate regions have revealed anomalous heat conduction, with nonlinear temperature profiles and thermal conductivities increasing with system size lin11 ; bermudez13 ; ruiz14 ; freitas15 . It has been shown that the linear chain exhibits almost flat temperature profiles characteristic of harmonic systems. A non-zero temperature gradient can be induced by engineered on-site disorder due to spin-vibron couplings, and dephasing through noisy modulations of the trap frequencies bermudez13 . Also, a small amount of induced disorder using site-specific dipole forces can be used to control the transition from anomalous to normal heat transport in two- and three-dimensional crystals freitas15 . In this context, a proposal for the experimental implementation of local dephasing noise on the vibrational excitations by means of fluctuating optical potentials has been described cormick16 .

Most of the previously mentioned studies of heat transport in trapped ion systems have considered the potential interaction in harmonic approximation about the equilibrium positions of the ions in the different spatial configurations lin11 ; bermudez13 ; freitas15 . Therefore, any possible normal heat conduction behavior had to be necessarily induced by the artificial inclusion of mechanisms such as disorder and dephasing. However, it has been shown that non-zero temperature gradients naturally arise in models that fully take into account the many-body Coulomb interaction, due to non-linearity and axial-transverse mode coupling effects arising in the proximity of the structural phase transitions ruiz14 .

The aim of this work is to study the correlation between the degree of atomic delocalization in the steady nonequilibrium dynamics of classical systems and their heat transport properties. To illustrate the analysis, we focus on a system formed by a Coulomb crystal of trapped ions, across the various structural phase transitions. We shall study a three-dimensional model corresponding to the design of a possible experiment to measure an energy current through the system. This theoretical model will be used to numerically simulate the classical dynamics of the system in contact with laser beams that emulate two heat reservoirs at different temperatures in two separate regions, until reaching the nonequilibrium steady state. Then, heat transport properties, such as temperature profiles and the total heat flux, can be obtained from the dynamical variables in such a state. We will show that the internal ions can exhibit a strong delocalization, and that such behavior is correlated with the onset of Fourier-type temperature profiles. A proper treatment of such delocalization will require a continuous description of the transport properties.

The paper is organized as follows. In Section II we describe the general model considered to study the nonequilibrium dynamics and heat transport in classical systems with atomic delocalization. The nonequilibrium dynamics in the steady state is analyzed in terms of spatial probability densities, which will evidence the degree of atomic delocalization. We shall consider a continuous description to define the steady temperature profile and the total heat flux in terms of dynamical variables. In Section III we particularize the model to simulate three-dimensional Coulomb crystals of trapped ions, and set the different parameters corresponding to a possible experimental setup. We show the numerical results concerning the steady state nonequilibrium dynamics of the ions for various anisotropies of the trapping potential, in particular the spatial probability densities of the entire systems and the spatial distributions of the individual ions. The temperature profiles and the total heat flux through the various structural phase transitions that modify the effective dimensionality of the trapped ion system are also shown. Finally, Section IV summarizes the main conclusions.

II Nonequilibrium dynamics and heat transport in classical systems with atomic delocalization

We consider a classical three-dimensional system composed of $N$ particles of mass $m$, whose motional degrees of freedom are described by the position coordinates ${\bf q}_{i}=(q_{x,i},q_{y,i},q_{z,i})$ and their conjugate momenta ${\bf p}_{i}=(p_{x,i},p_{y,i},p_{z,i})$, with $i=1,\dots,N$. We assume that the particles interact with each other through a central interaction potential $U$, and that the entire system remains confined within a finite volume due to an external trapping potential $V$. Then, the dynamics of the system can be described by the Hamiltonian

[TABLE]

A main goal of this work is to analyze the response of this system to a imposed temperature gradient, as a function of the trapping potential $V$. The variation of such potential can be used to induce structural phase transitions that modify the effective dimensionality of the system, and therefore to control the nonequilibrium dynamics and the corresponding heat transport properties.

II.1 Nonequilibrium dynamics

To induce a heat current across a given direction, we consider that $N_{L}$ particles on the left end along this direction and $N_{R}$ on the right one are connected to thermal reservoirs at different temperatures. We shall analyze a regime in which the localization induced by the interaction with the thermal reservoirs keeps these extreme particles on both ends within the spatial regions where such reservoirs are acting, whereas the remaining internal particles may be delocalized within the intermediate region. Assuming Langevin thermal reservoirs, the equations of motion for the $\alpha=(x,y,z)$ components of the position and momentum coordinates can be expressed as:

[TABLE]

where $g_{\alpha,i}=-\partial V_{i}({\bf q}_{i})/\partial q_{\alpha,i}$ is the external force along the $\alpha$-direction, and $f_{\alpha}^{(ij)}=-f_{\alpha}^{(ji)}=-\partial U_{ij}(|{\bf q}_{i}-{\bf q}_{j}|)/\partial q_{\alpha,i}$ the force that the $j$th particle exerts on the $i$th particle along such direction. The action of the Langevin reservoirs is characterized by the friction coefficients $\eta_{\alpha,i}^{L,R}$ and the stochastic forces $\varepsilon_{\alpha,i}^{L,R}(t)$. This force is assumed to correspond to a Gaussian white noise that satisfies the statistical relationships

[TABLE]

where $\langle\dots\rangle$ denote the average over an ensemble of stochastic trajectories, and $D_{\alpha,i}^{L,R}$ are the diffusion coefficients. These are related to the friction coefficients $\eta_{\alpha,i}^{L,R}$ according to the fluctuation dissipation theorem kubo66

[TABLE]

with $T^{L,R}$ the temperature of the corresponding thermal reservoir.

The equations of motion (2) can be rewritten in terms of both the friction and diffusion coefficients as the following stochastic differential equations:

[TABLE]

where $dW_{\alpha,i}^{L,R}$ denote the Wiener processes associated with the interactions with the laser reservoirs.

The heat transport properties of the trapped system in contact with the two thermal reservoirs are dictated by the steady state solution of the equations of motion (5). To elucidate the underlying nonequilibrium dynamics in such state we shall analyze the probability density of particles in the spatial domain ${\bf q}=(x,y,z)$. This local distribution is obtained from the set of positions that the particles have visited in their dynamics during a sufficiently long time interval $\tau_{ss}$, as:

[TABLE]

with $t$ an arbitrary time value within the steady state, and $\sigma$ a small parameter giving the width of the three-dimensional Gaussian kernel. The purpose of introducing the smoothing in $P({\bf q})$ is to highlight the spatial regions with high probability of finding the particles.

In addition to the spatial probability density $P({\bf q})$ of the entire system, valuable information of the dynamics can be extracted from the analysis of the steady spatial distributions of the individual particles. In order to visualize such individual distributions, we divide the spatial region occupied by the entire system along a given $\alpha$-direction in a series composed of $c_{\alpha}$ cells, centered on the positions $\alpha_{l}\,(l=1,\dots,c_{\alpha})$ and with size $\Delta$ alonso99 . During a sufficiently long time interval $\tau_{ss}$ in the nonequilibrium steady state, we monitor the passage of the particles through the various cells, and determine the time spent within each of them on successive visits. In this way, we obtain the spatial distribution of the $i$th particle along the $\alpha$-direction as:

[TABLE]

In order to obtain a quasi-continuous distribution $\Theta_{i}(\alpha)$, a small enough value of $\Delta$ shall be considered.

As we will show below, strong trapping potentials $V$ lead to point like spatial distributions in which the individual particles can be clearly distinguished, whereas spatially extended distributions emerge in weaker confinements in which the internal particles can become highly delocalized.

II.2 Heat transport properties

In this section we focus on the study of the temperature profiles and the total heat flux through a selected direction, obtained from the position $\{{\bf q}\}_{N}=({\bf q}_{1},\dots,{\bf q}_{N})$ and the momentum $\{{\bf p}\}_{N}=({\bf p}_{1},\dots,{\bf p}_{N})$ coordinates of the particles in the nonequilibrium steady state reached under the action of the thermal reservoirs. We are particularly interested in analyzing the behavior of these heat transport properties through the various structural phase transitions that the system may experience due to controlled variations of the trapping potential $V$.

II.2.1 Temperature profiles

Taking into account that the internal particles may be delocalized for some configurations of the trapping potentials, we shall consider a continuous description to define the steady local temperature $T(\alpha)$ across a given $\alpha$-direction. To proceed, here again we consider a series composed of $c_{\alpha}$ cells along such direction, and monitor the passage of the particles through each cell during a sufficiently long time interval $\tau_{ss}$ in the nonequilibrium steady state alonso99 . We then make use of the equipartition theorem to write the temperature of the $l$th cell in terms of the kinetic energy of the particles as:

[TABLE]

where $E_{k}({\bf p}_{i})$ is the kinetic energy of the $i$th particle, and $k_{B}$ the Boltzmann constant.

II.2.2 Total heat flux

To obtain the heat flux across the trapped system connected to different heat reservoirs in two separate regions, we continue using a continuous description and consider the energy balance equation in local form lepri03 ; dhar08 ; zubarev74 . To proceed, we start by introducing the dynamical variable corresponding to the local energy density

[TABLE]

with

[TABLE]

The time derivative of $h({\bf q})$, taking into account the equations of motion (2), is given by:

[TABLE]

where

[TABLE]

with $\alpha$ running over the components $\{x,y,z\}$, is the energy per time unit exchanged between the $i$th particle and the heat reservoir to which it is connected.

Equation (11) can be rewritten as:

[TABLE]

where ${\bf f}^{(ij)}=-{\bf f}^{(ji)}$ has been used. To transform the last expression into a local energy balance equation, it is necessary to rewrite the third term on the right hand side in the form of a divergence. For this purpose, we take into account that

[TABLE]

with ${\bf G}(\lambda)$ an arbitrary function satisfying the conditions ${\bf G}(1)={\bf q}_{i}$ and ${\bf G}(0)={\bf q}_{j}$ piccirelli68 . Then, we get the balance equation

[TABLE]

with the energy flux density vector

[TABLE]

and the energy source term

[TABLE]

Finally, the integral of ${\bf j}_{h}({\bf q})$ over all space gives the total heat flux

[TABLE]

In the steady state, considering that

[TABLE]

where $\langle\dots\rangle_{ss}$ indicates the steady state average, the expected value of the total heat flux can be expressed as:

[TABLE]

The prime denotes the sum over the particles that are connected to heat reservoirs, and $B_{i}\equiv L$ for $i=1,\dots,N_{L}$ and $B_{i}\equiv R$ for $i=N-N_{R}+1,\dots N$. According to the Novikov’s theorem novikov65 , and taking into account the stochastic relationships (3) and the equations of motion (2), the average of the terms including the stochastic forces are given by:

[TABLE]

with $(\eta,\alpha,\nu)$ running over the components $\{x,y,z\}$. Then, the steady total heat flux becomes

[TABLE]

The requirement that the ensemble average $\langle\,{\bf J}(t)\,\rangle$ of the total heat flux (21) over a large enough number of stochastic trajectories coincides with the result given by (25) provides a good criterion for checking convergence to steady state in the numerical simulations.

III Trapped ion systems

In this section we deal with a three-dimensional system composed of $N$ ions of mass $m$ and charge $Q$, confined within an electromagnetic trap. We focus on the motional degrees of freedom of the ions, described by the position coordinates ${\bf q}_{i}=(q_{x,i},q_{y,i},q_{z,i})$ and their corresponding momenta ${\bf p}_{i}=(p_{x,i},p_{y,i},p_{z,i})$. We consider a system with a small number of ions, and use the pseudopotential theory to replace the trapping potential by a time-independent harmonic potential ghosh95 . We assume that this approximation, which neglects the rapid micromotion, captures the essence of the secular motion. Specifically, we assume that the ions are confined by a harmonic trap with the axial frequency $\omega_{x}$, and the transverse (radial) frequencies $\omega_{y}$ and $\omega_{z}\,$. To characterize the radial anisotropy of the trap we introduce the parameters $n_{\beta}=\omega_{\beta}/\omega_{x}$, with $\beta=(y,z)$. Then, the dynamics of the system is ruled by both the interaction with the harmonic trap

[TABLE]

and the Coulomb repulsion

[TABLE]

with $\epsilon_{0}$ the vacuum permittivity.

In this trapped ion system we study the steady state nonequilibrium dynamics as well as heat transport, through the various structural phase transitions induced by the variation of anisotropy of the trapping potential. We focus on the study of the spatial distributions of the ions, and the temperature profiles and the total heat flux across the axial direction. To induce a heat current across the system, we consider that the $N_{L}$ leftmost ions along the $x$-direction and the $N_{R}$ rightmost ones are connected to laser beams that simulate two thermal reservoirs at different temperatures. In order to resolve the dynamics we assume that such laser beams can be modeled as Langevin heat baths. In addition, considering that the typical separations between the ions are generally of the order of micrometers, we adopt a classical description of the nonequilibrium dynamics. Then, such dynamics can be described by the equations of motion (5).

For small laser intensities, the friction coefficients $\eta_{\alpha,i}^{L,R}$ and the diffusion coefficients $D_{\alpha,i}^{L,R}$ can be obtained from the Doppler cooling expressions phillips92 :

[TABLE]

and

[TABLE]

The ratio $I_{\alpha,i}^{L(R)}/I_{0}^{L(R)}$ denotes the normalized intensity of the laser beam acting on the $i$th ion along the $\alpha$-direction, $k^{L(R)}_{\alpha,i}$ is the corresponding laser wavelength, $\delta_{\alpha,i}^{L(R)}=\omega_{\alpha,i}^{L(R)}-\omega_{0}$ is the detuning of the laser frequency $\omega_{\alpha,i}^{L(R)}$ with respect to the frequency $\omega_{0}$ of a selected atomic transition in the ion, and $\Gamma$ is the natural linewidth of the excited state in such transition. In this work we select the atomic transition $3s^{2}S_{1/2}\rightarrow 3p^{2}P_{1/2}$ of the 24Mg+ ions, with $\omega_{0}/2\pi=1069$ THz and $\Gamma/2\pi=41.296$ MHz nist99 . To induce a heat current through the trapped ion system we shall consider that the extreme ions on both ends are subjected to laser beams with different detunings $\delta_{\alpha,i}^{L}\neq\delta_{\alpha,i}^{R}$. From now on we set $\delta_{\alpha,i}^{L}=-0.02\Gamma$ for the $N_{L}$ leftmost ions and $\delta_{\alpha,i}^{R}=-0.1\Gamma$ for the $N_{R}$ rightmost ones, and the same laser intensity $I_{\alpha,i}^{L(R)}/I_{0}=0.08$ on both ends. These values lead to the friction coefficients $\eta^{L}=6.76\times 10^{-22}\,$kg$/$s and $\eta^{R}=3.13\times 10^{-21}\,$kg$/$s, and the diffusion coefficients $D^{L}=1.16\times 10^{-46}\,$kg2m2$/$s3 and $D^{R}=1.11\times 10^{-46}\,$kg2m2$/$s3. For reference, the corresponding limit Doppler temperatures, obtained in the case of the Doppler cooling of a single isolated ion, would be $T^{L}=D^{L}/k_{B}\eta^{L}=12.41\,$mK and $T^{R}=D^{R}/k_{B}\eta^{R}=2.58\,$mK. Thus the trapped ion system shall be connected to an effective hotter bath on the left end, and to a colder bath on the right end.

III.1 Spatial configurations of the trapped ions

Before dealing with the study of the heat transport properties, we analyze the underlying nonequilibrium dynamics of the trapped ions according to the anisotropy of the confining potential. In this section we present a detailed analysis of the steady state spatial distribution of the trapped ions as a function of the parameters of the radial anisotropy $(n_{y},n_{z})$. We are particularly interested in identifying the values of these parameters leading to the different structural phase transitions.

To perform the numerical analysis, from now on we shall consider a system composed of $N=30$ 24Mg+ ions, with $N_{L}=N_{R}=3$ of them connected to the laser reservoirs on both ends along the axial direction. We set the axial frequency $\omega_{x}/2\pi=50\,$kHz, and study the dynamics for different transverse frequencies $(\omega_{y},\omega_{z})$. In order to ensure that the extreme ions that are connected to the lasers reservoirs remain within the space region where the beams are focussed, the radial frequencies must be sufficiently high. For a fixed value of $\omega_{x}$, this imposes a lower limit to the values of $(n_{y},n_{z})$. In the case of the selected value of $\omega_{x}$, this requires considering values of both anisotropy parameters above approximately 6.7.

In the numerical simulations we set the initial conditions with the ions at rest, and positions randomly distributed in the close vicinity of the equilibrium positions of the linear configuration along the axial direction. Then the laser reservoirs that are focused on the selected ions are switched on instantaneously. We perform the time evolution until the system reaches the nonequilibrium steady state from which the energy transport properties are extracted. Figures 1, 2 and 3 show the steady state spatial probability densities (6) of the ions for different anisotropies of the trapping potential.

Figure 1 illustrates the change of the spatial probability distribution through a linear-zigzag structural phase transition. It shows how for a fixed axial frequency $\omega_{x}$ such transition spreads from the center to the edges as the radial anisotropy $n_{y}$ is lowered, whereas the high radial anisotropy $n_{z}$ confines the ions in the $xy$-plane.

As Fig. 2 shows, in the case of a symmetrical trap with equal radial frequencies, their decrease gives rise to a transition from the linear configuration to a three-dimensional one, in which the ions distribute over a series of rings contained in the transverse $yz$-plane and centered along the axial direction $x$. As occurs with the zigzag configuration, the rings arise at the center of the system and extend towards the ends as the radial frequencies decrease. This three-dimensional configuration of the trapped ions corresponds to the helical arrangement reported in previous studies waki92 ; hasse90 ; nigmatullin16 . From now on we will refer to it as the helical configuration.

Outside the high frequency domain corresponding to the linear configuration, the variation of one of the transverse frequencies through the helical configuration $(n_{y}=n_{z})$ results in a rotation of the zigzag configuration around the axial axis, to be confined again in the plane perpendicular to transverse direction with the highest anisotropy. As an illustration, Fig. 3 shows the transition from the zigzag configuration confined in the $xz$-plane, given by $(n_{y},n_{z})=(9,8)$, to the zigzag configuration in the $xy$-plane, given by $(n_{y},n_{z})=(9,10)$, through the helical configuration corresponding to $(n_{y},n_{z})=(9,9)$.

So far we have used the spatial probability densities of the entire trapped ion system to illustrate the different structural phase transitions that it may experience. Notice that while each of the $N_{L}=N_{R}=3$ outermost dots, as well as all dots in the linear chain, can be assigned to specific ions, we will now show that this not necessarily the case for the internal dots of the zigzag configuration, nor for the various rings of the helical configuration. To elucidate the dynamics performed by the individual ions in the three previously shown spatial configurations, we now analyze the steady spatial distributions (7). Figure 4 illustrates some of such individual distributions along the axial direction, and Fig. 5 along the two transversal directions.

As expected, in the linear chain $\Theta_{i}(\alpha)$ exhibits a single peak centered at the corresponding equilibrium position along the axial axis, for all ions $i=1,\dots N$ and directions $\alpha=(x,y,z)$, see Figs. 4 and 5. In this configuration each ion is strongly confined, and can only perform small oscillations around its equilibrium position.

In the zigzag and helical configurations the distributions $\Theta_{i}(x)$ of the internal ions, see Fig. 4, exhibit a series of peaks, which evidence the delocalization of these ions along the axial direction. The intensities of the various peaks vary from one ion to another, but their positions are the same, as they are determined by the minima of the global potential energy surface. In the zigzag configuration the different peaks correspond to equilibrium positions located on both sides of the axial axis, see for example Fig. 1. Whereas in the helical configuration they give the location of the series of rings shown in Fig. 2. Thus, the lower confinement of the internal ions allows them to move throughout the system and exchange their positions along the axial direction. This axial displacement can occur along practically all the region covered by the helical configuration, through very fast jumps between the different rings. In the zigzag configuration such displacement is more local, as it is restricted to rapid jumps between neighboring equilibrium positions.

The radial distributions $\Theta_{i}(y)$ and $\Theta_{i}(z)$ clearly distinguish between the linear, zigzag and helical configurations, see Fig. 5. In the zigzag configuration the distribution of the radial coordinate corresponding to the highest trapping frequency presents a single peak centered at zero, while the other radial coordinate exhibits two peaks of similar intensity arranged symmetrically around zero. The probability of presence practically zero in between the two peaks indicates that the ions are jumping very rapidly through the axial axis, staying most of the time in the close vicinity of any of the minima that the potential has on either side of this axis. In the helical configuration both radial coordinates show nearly identical bimodal distributions, again arranged symmetrically around zero. But in this case there is a significant probability of presence in between both maxima, which evidences the distribution of the ions within the rings shown in Fig. 2.

Taking into account the delocalization that the internal ions can exhibit, we shall characterize the configuration changes using an alternative order parameter to those considered in previous studies schiffer93 ; freitas15 ; yan16 . Concretely, we have shown that the two radial distributions $\Theta_{i}(y)$ and $\Theta_{i}(z)$ are highly sensitive to changes in the global potential energy surface leading to the structural phase transitions. We now employ the locations of their maxima as a criterion to construct an anisotropy map that shows the values of $(n_{y},n_{z})$ at which the different spatial configurations occur. To be more precise, the maxima of both distributions are located at zero for an ion aligned with the chain axis, only one of the two distributions has the maximum at zero for an ion that is in a zigzag configuration, whereas neither of them presents a maximum at zero for an ion arranged in a helical configuration. The anisotropy map obtained for the ion distributions in the three different spatial configurations is given in Fig. 6.

As expected, all the ions are located along the axial direction in trapping potentials with high radial anisotropy. In the case of a chain composed of $N=30$ 24Mg+ ions, we observe that this occurs for values of both $n_{y}$ and $n_{z}$ above approximately $11.6$, see Fig. 6. Outside such region the trapped ion system exhibits predominantly zigzag configurations in the plane perpendicular to the radial direction with the largest anisotropy, and with a decreasing number of external ions along the axial axis as the radial trapping frequencies become smaller. The helical configuration does not occur unless the two radial trapping frequencies become practically equal. Thus, it emerges as a distinctive feature of traps with symmetrical anisotropy, where $\omega_{y}=\omega_{z}\,$, in contrast to the ubiquitous zigzag configurations for radially asymmetric traps, in which $\omega_{y}\neq\omega_{z}$, provided at least one of the two radial frequencies is sufficiently small.

III.2 Temperature profiles and total heat flux

Now we analyze the steady state temperature profiles and the total heat flux across the axial direction, for various anisotropies of the trapping potential. In particular, we focus on their behavior through the various structural phase transitions shown in the previous section.

Figure 7(a) shows the temperature profile and the total heat flux across a linear-zigzag structural phase transition, Fig. 8(a) in a linear-helical transition, and Fig. 9(a) in a transition between two perpendicular zigzag configurations through a helical configuration.

The presence of separate segments across the various temperature fields is due to the strong spatial confinement of certain ions around their equilibrium positions. In agreement with the results presented in the previous section, such confinement persists for all ions in the linear string, while in the zigzag and helical configurations the delocalization of the internal ions leads to quasi-continuous central regions in the temperature profiles.

The small size of the system leads to significant boundary effects in the temperature profiles, mainly in the regions occupied by the ions that are connected to the thermal baths and their nearest neighbors. In between those regions, the temperature field corresponding to the internal ions becomes very sensitive to the structural phase changes. In the linear chain it shows a nearly perfect plateau. As is known, a flat temperature profile evidences the anomalous heat transport in one-dimensional harmonic crystals, which have infinite conductivity and can not, therefore, support a temperature gradient rieder67 ; dhar16 .

The structural phase transitions from the linear to the zigzag configuration and from the linear to the helical configuration introduce a non-zero temperature gradient with the axial coordinate, whose magnitude increases as the transitions spread from the center to the edges. While in the zigzag configuration the overall gradient is not uniform, in the helical one it remains almost constant, giving rise to a quasi-linear temperature profile consistent with Fourier’s law.

These results elucidate a correlation between the amount of delocalization of the internal ions and the temperature profiles across the axial direction. While an almost complete delocalization in the helical configuration leads to a Fourier like behavior, the strong confinement in the linear chain results in the nearly flat profiles characteristic of an anomalous heat transport in harmonic systems. The more restricted delocalization in the zigzag configuration corresponds to an intermediate behavior, with non fully linear temperature profiles.

As Fig. 10(a) summarizes, the temperature profiles corresponding to the complete linear, zigzag and helical configurations can be clearly differentiated from each other. The analysis of the temperature gradient across the axial region occupied by the internal ions clearly shows the different behaviors of the three configurations, see Fig. 10(b). As expected, this gradient remains very close to zero in the linear chain. In the helical configuration it presents very small variations around a non-zero value, which is consistent with Fourier’s law. While the zigzag configuration clearly departs from the behavior predicted by this law, by exhibiting a significant variation of the temperature gradient.

The temperature profiles obtained for the linear and zigzag configurations are in agreement with those previously reported using a discrete description in a two-dimensional model ruiz14 , as in such configurations the delocalization of the internal ions is absent or remains restricted. However, as shown in this work, a proper analysis of the helical configuration necessarily requires a continuous description that takes into account the displacement of the internal ions across the axial direction.

The spectral analysis of the steady state evolution of the coordinates of the internal ions has shown that the formation of a non-zero temperature gradient, as the system experiences a structural phase transition from the linear chain to configurations of higher dimensionality, can be assigned to non-linearities and the coupling between radial and axial modes, which are absent in the harmonic picture ruiz14 . Such effects become also visible in the total heat flux across the axial direction, which exhibits signatures of the structural phase transitions, see Figs. 7(b), 8(b) and 9(b).

On approaching the transitions from the high anisotropy domain corresponding to the linear chain, the increasing contribution of the transverse modes and the growing level of fluctuations due to thermal motion of the ions may assist transport, leading to a progressive increase of the heat flux ruiz14 . As Fig. (7)(b) and (8)(b) show, once the transition has already emerged and the chain buckles, a further decrease of the radial anisotropy leads to a reduction of the total heat flux as the transition spreads from the center to the edges and the internal ions jump off the axial axis.

The reduction of the heat flux can be attributed to the decreasing interaction between neighboring ions as they arrange in configurations of higher dimensionality, in which the inter-ion distances become larger. While in the linear-helical transition such reduction is nearly uniform, in the linear-zigzag transition it tends to stabilize at low anisotropies, once all the internal ions have jumped off the axial axis. In the case of a system composed of $N=30$ 24Mg+ ions, the total reduction of the heat flux through the transition from the linear chain to the complete zigzag and helical configurations reaches around 34%.

The analysis of the linear-zigzag and linear-helical transitions has shown that heat transport is optimal in the linear configuration, in the proximity of the onset of the structural phase transition. According to Fig. 9(b), in the case of the transition between two perpendicular zigzag configurations, the heat flux exhibits a maximum at the intermediate helical configuration. Thus the coupling between the axial and transverse modes becomes more detrimental to heat transport in the zigzag configuration than in the helical configuration with a similar number of internal ions located outside the axial axis.

IV Conclusions

We have used Coulomb crystals of trapped ions to get deeper insight into the connection between heat transport properties and the underlying nonequilibrium dynamics in systems that are in contact with different thermal baths in two separate regions. We have considered an intrinsically non-linear three dimensional model, which fully takes into account the many-body Coulomb interaction, and analyzed the response of the system through the various structural phase transitions induced by the controlled variation of the anisotropy of the trapping potential.

The results, obtained from the numerical resolution of the classical Langevin equations of motion, have shown a correlation between the degree of delocalization of the internal ions and the temperature gradient across the axial direction. The strong confinement of the ions around their equilibrium positions in the linear chain leads to nearly flat temperature profiles characteristic of the anomalous heat conduction one-dimensional harmonic systems. While the extended delocalization of the ions in the helical configuration is associated with global quasi-linear temperature profiles, as predicted by Fourier’s law. The planar zigzag configuration corresponds to an intermediate situation, in which the more restricted delocalization results in non-zero, but non-uniform temperature gradients.

Considering that Fourier’s law is a macroscopic consequence of ordinary diffusion at the microscopic level, the observed delocalization of the internal ions across the entire helical configuration could correspond to a ordinary Brownian diffusion. This is not expected in the case of the zigzag configuration, for which the temperature profiles do not fully exhibit a Fourier-type behavior.

While previous studies based on harmonic models have shown that the onset of nonzero temperature gradients across the trapped ion system can be artificially induced by engineered dephasing and disorder, in this work we show that the transition from anomalous to normal transport arises naturally in an intrinsically non-linear model. The interplay between the many-body Coulomb interaction and the external substrate potential of the trap leads to a very rich underlying nonequilibrium dynamics, which is ultimately responsible for the strong dependence of the transport properties on the structural phase transitions that modify the dimensionality of the system.

We have shown that the total heat flux across the axial direction is highly sensitive to changes in the effective dimensionality of the trapped ion system. Heat transport is optimal in the linear configuration, in the proximity of the onset for the structural phase transitions to configurations of larger dimensionality. This may be attributed to an increasing contribution of the transverse modes to transport, and the increasing thermal motion of the ions. Upon further decrease of the anisotropy of the trapping potential, the spread of the structural phase transitions across the axial direction results in a progressive decrease the total heat flux, as the larger inter-ion distances in the zigzag and helical configurations reduce the interaction between neighboring ions.

The transition through a helical configuration in between two perpendicular zigzag configurations results in a local maximum of the total heat flow. Thus, the non-linear effects that arise in the dynamics during the transition to the planar zigzag configuration are more detrimental to heat transport than those corresponding to the helical configuration.

Interesting issues to consider in future work include the study of additional transport properties, such as the conductance, and their correlation with the atomic delocalization in systems of different dimensionality and size. In the case of Coulomb crystal of trapped ions, the study of heat transport properties in systems with different ion species or with structural defects are also of great interest.

Acknowledgements.

We thank J. P. Palao and J. Onam González for fruitful discussions, and J. González-Platas for his guidance on the use of the VESTA software. This project was funded by the Spanish MICINN and European Union (FEDER) (FIS2017-82855-P).