Heat rectification with a minimal model of two harmonic oscillators

We study heat rectification in a minimalistic model composed of two masses subjected to on-site and coupling linear forces in contact with effective Langevin baths induced by laser interactions. Analytic expressions of the heat currents in the steady state are spelled out. Asymmetric heat transport is found in this linear system if both the bath temperatures and the temperature dependent bath-system couplings are also exchanged.


I. INTRODUCTION
Heat rectification, firstly observed in 1936 by Starr [1], is the physical phenomenon, analogous to electrical current rectification in diodes, in which heat current through a device or medium is not symmetric with respect to the exchange of the baths at the boundaries. In the limiting case the device allows heat to propagate in one direction from the hot to the cold bath while it behaves as a thermal insulator in the opposite direction when the baths are exchanged. In 2002 a paper by Terraneo et al. [2] demonstrated heat rectification numerically for a chain of nonlinear oscillators in contact with two thermal baths at different temperatures. Since then, there has been a growing interest in heat rectification [3][4][5][6][7][8][9][10][11][12][13][14][15][16], and the field remains very active because of the potential applications in fundamental science and technology, and the fact that none of the proposals so far appears to be efficient and robust for practical purposes.
Much effort has been devoted to understand the underlying physical mechanism responsible for rectification [3]. In early times some kind of anharmonicity, i.e. non-linear forces, in the substrate potential or in the particle-particle interactions, was identified as a fundamental requisite for rectification [5,[17][18][19][20][21]. This nonharmonic behavior leads to a temperature dependence of the phonon bands. The match/mismatch of the phonon bands (power spectra) governs the heat transport in the chain, allowing it when the bands match or obstructing it if they mismatch [2,22]. However, a work by Pereira et al. [23] showed that rectification can also be found in effective harmonic systems if two requirements are met: some kind of structural asymmetry, and features that depend on the temperature so they change as the baths are inverted. Indeed, in this article we demonstrate rectification in a minimalistic model of two harmonic oscillators where the coupling to the baths depends on the temperature. This will be justified with a particular physical set up with trapped ions and lasers. * miguelangel.simon@ehu.eus † jg.muga@ehu.es The article is organized as follows. In Section II we describe the physical model and its dynamical equations. In Section III we describe the dynamics of the system in terms of a covariance matrix. We also derive a set of algebraic equations that gives as solution the covariance matrix in the steady state. In Section IV we solve the covariance matrix equations and find analytical expressions for the steady-state temperatures of the masses and heat currents. In Section V we relate the parameters of our model to those in a physical set-up of Doppler cooled trapped ions. In Section VI we make a parameter sweep looking for configurations which yield high rectification. We also study the power spectra of the oscillators, which confirm the match/mismatch patterns in cases where there is rectification. In Section VII we summarize our results and present our conclusions.

II. PHYSICAL MODEL
The physical model consists of two masses m 1 and m 2 coupled to each other by a harmonic interaction with spring constant k and natural length x e . Each of the masses m 1 and m 2 are confined by a harmonic potential with spring constants k L , k R and equilibrium positions x L , x R respectively (see Fig. 1). The Hamiltonian de-scribing this model is 2 are the position and momentum of each mass. Switching from the original coordinates x i to displacements with respect to the equilibrium positions of the system q i = x i − x eq i , where x eq i are the solutions to ∂ xi V (x 1 , x 2 ) = 0, the Hamiltonian can be written as This has the form of the Hamiltonian of a system around a stable equilibrium point where is the mass matrix of the system and K is the Hessian matrix of the potential at the equilibrium point, i.e., We shall see later that the generic form (3) can be adapted to different physical settings, in particular to two ions in individual traps, or to two ions in a common trap.
The masses are in contact with Langevin baths, which will be denoted as L (for left) and R (for right), at temperatures T L and T R for the mass m 1 and m 2 respectively (see Fig. 1). The equations of motion of the system, taking into account the Hamiltonian and the Langevin baths areq where γ L , γ R are the friction coefficients of the baths and ξ L (t), ξ R (t) are Gaussian white-noise-like forces. The Gaussian forces have zero mean ( ξ L (t) = ξ R (t) = 0) and satisfy the correlations ξ L (t)ξ R (t ) = 0, . D L and D R are the diffusion coefficients, which satisfy the fluctuation-dissipation theorem: It is useful to define the phase-space vector − → r (t) = − → q , M −1− → p T (note that − → v = M −1− → p is just the velocity vector) so the equations of motion for this vector arė with and − → ξ (t) = (ξ L (t), ξ R (t)) T , Γ = diag(γ L , γ R ). 0 n×n and 1 n×n are the n-th dimensional squared 0 matrix and identity matrix respectively. With the vector notation the correlation of the white-noise forces can be written as with D = diag(D L , D R ).

III. COVARIANCE MATRIX IN THE STEADY STATE
We define the covariance matrix of the system as This matrix is important because the heat transport properties can be extracted from it. In particular, the kinetic temperatures of the masses, T 1 (t) and T 2 (t), are One approach to find the covariance matrix is to solve Eq. (5). However, this requires solving the equations explicitly or simulate them numerically many times to find the covariance matrix for the ensemble of simulated stochastic trajectories. Instead, we proceed by looking for an ordinary differential equation that gives the evolution of the covariance matrix as described in [24][25][26]. Differentiating C(t) with respect to time and using Eq.
The solution of Eq. (9) allows us to find the local temperatures of the masses as a function of the bath temperatures (Eq. (8)) at all times. In particular, we are interested in the covariance matrix in the steady state, i.e., for t → ∞. According to the Novikov Theorem [27] we can write down the covariance matrix in the steady state without having to integrate the differential equation. We now show how to get the steady-state covariance matrix.
In the steady state, the covariance matrix is constant ( d dt C(t) = 0), therefore it satisfies AC s.s. + C s.s. without having to integrate the equations of motion. Using this theorem and the δ-correlation of the noises, we find the ij-th where lim τ →t − is the limit when τ goes to t from below. Evaluation of the functional derivative δr j (t)/δξ k (τ ) for the τ → t − limit gives Now, the algebraic equation that gives the steady-state covariance matrix becomes with B = 2LDL T . By definition, the covariance matrix is symmetric, but there are also additional restrictions imposed by the equations of motion and the steady-state condition, which reduce the dimensionality of the problem of solving Eq. (13) [28]. Since d q i q j /dt = 0 in the steady state, we have Taking (14) into account, the steady-state covariance matrix takes the form The explicit set of equations for the components of C s.s can be found in Appendix A.

IV. SOLUTIONS
In this section we use the solution to Eq. (13) to write down the temperatures and currents in the steady state. We use Mathematica to obtain analytic expressions for the temperatures, where D(k) = where h (n) ≡ γ R m n 1 +γ L m n 2 . The currents from the baths to the masses [28] are with T i given by Eq. (16). Since, in the steady state, J L = −J R we will use the shorthand notation J ≡ J L . Substituting Eq. (16) into Eq. (18) we get for the heat current where κ = k B k 2 γ L γ R h (1) /D(k) acts as an effective thermal conductance, which depends on the parameters of the system, i.e., the masses and spring constants, and also on the friction coefficients of the baths. From Eq. (19) it could be thought that inverting the temperatures of the baths would only lead to an exchange of heat currents. However, since the thermal conductance κ depends on the friction coefficients, the exchange of the baths implies a change in its value. Moreover, it is possible to have temperature-dependent friction coefficients, as it happens in the physical set-up of laser-cooled trapped ions described in Section V.

V. RELATION OF THE MODEL TO A TRAPPED ION SET-UP
As we mentioned, the parameters k, k L and k R can be related to the elements of the Hessian matrix of a system in a stable equilibrium position. In this section we will identify these parameters with the Hessian matrix of a pair of trapped ions. Here we consider two different set-ups: two ions in a collective trap, and two ions in individual traps. In Section VI we focus on two ions in individual traps to illustrate the analysis of rectification.
In both set-ups we assume strong confinement in the radial direction, making the effective dynamics onedimensional. We will also assume that the confinement in the axial direction is purely electrostatic, which makes the effective spring constant independent of the mass of the ions [29]. Additionally, we will relate the temperatures and friction coefficients of the Langevin baths to those corresponding to Doppler cooling.

A. Collective trap
Consider two ions of unit charge with masses m 1 and m 2 trapped in a collective trap. Assuming strong radial confinement and purely electrostatic axial confinement, both ions feel the same harmonic oscillator potential with trapping constant k trap [29]. The potential describing the system is with C = Q 2 4πε0 . The equilibrium positions for this potential are Assuming small oscillations of the ions around the equilibrium positions, the Hessian matrix of the system is Using Eq. (22) we can relate the parameters of this physical set-up to those of the model described in Section II to find

B. Individual on-site traps
We can make the same assumptions for the axial confinement as in the previous subsection but now each of the ions is in an individual trap with spring constants k trap,L and k trap,R respectively. The potential of the system is where x L and x R are the center positions of the on-site traps. The elements of the Hessian matrix in the equilibrium position are Comparing the parameters in Eq. (25) with those in the model described in Section II we identify In this case, the analytic expressions for the equilibrium positions are more complicated. We get for the distance between the equilibrium positions of the ions In this set-up, the coupling between the ions k can be controlled by changing the distance between the on-site traps.

C. Optical molasses and Langevin baths
Trapped ions may be cooled down by a pair of counterpropagating lasers which are red-detuned with respect to an internal atomic transition of the ions. This technique is known as Doppler cooling or optical molasses [30][31][32][33]. The off-resonant absorption of laser photons by the ions exerts a damping-like force that slows them down. The spontaneous emission of the ions produces heating due to the random recoil generated by the emitted photons. Both, the friction and recoil force are in balance, and eventually the ion thermalizes to a finite temperature. Thus the effect of the lasers on the ion is equivalent to a Langevin bath with temperature T molass and friction coefficient γ molass . The temperature and friction coefficients are controlled with the laser intensity I and frequency detuning δ with respect to the selected internal transition by the expressions [31,33,34], where ω 0 is the frequency of the selected internal atomic transition, Γ is the natural width of the excited state, and I 0 is the saturation intensity.

VI. LOOKING FOR RECTIFICATION
We will say that we observe rectification whenever the heat current J for a configuration of the baths changes when we exchange the baths toJ. The important point here is to define what is meant by exchanging the baths. We consider that a bath is characterized, not only by its temperature T but also by its coupling to the system by means of the friction coefficient γ, so, exchanging the baths is achieved by exchanging both the temperatures and the friction coefficients, as summarized in Table I. When implementing temperatures and friction coefficients by lasers, this exchange operation is performed by changing the values of the intensities and detunings acting on each ion (Eq. (29)). The exchange operation is straightforward when the two ions are either of the same species or isotopes of each other, since the only required action is to exchange the values of the detunings of the lasers without modifying the intensities. However, if we deal with two different species, i.e., with two different atomic transitions, the laser wavelengths and the decay rates depend on the species. Then, exchanging the temperatures by modifying the detunings, keeping the laser intensities constant, does not necessarily imply an exchange of the friction coefficients. Nevertheless it is possible to adjust the laser intensities so that the friction coefficients get exchanged and that is the assumption hereafter. The idea of implementing a bath exchange like this follows the same line of thought as [23], since we are adding a temperature dependent feature to the system -the friction coefficients-that changes as the baths are inverted.
To measure rectification, we will use the rectification coefficient R defined as that is, the ratio between the difference of heat currents and the largest one. As defined, R = 0 for no asymmetry of the heat currents and R = 1 when they are maximally asymmetric.

A. Parametric exploration
We have explored thoroughly the space formed by the parameters of the model to find asymmetric heat transport, namely, m 1 , m 2 , k, k L , k R , γ L , γ R . We have fixed the values of some of the parameters to realistic ones while we have varied the rest. We have set the masses to m 1 = 24.305 a.u. and m 2 = 40.078 a.u., which correspond to Mg and Ca, whose ions are broadly used in trapped-ion physics. The temperatures are also fixed and, as Eq. (19) shows, rectification does not formally depend on the temperature in this model, unless we set the friction coefficients as a function of temperature using Eq. (29) explicitly. Figure 2 depicts the values of the rectification after sweeping the k L k R plane for fixed values of k, γ L , and γ R . A remarkable result from this figure is that parallel lines appear alternating minima and maxima of R. With a numerical fitting, we find that the line corresponding to the highest maximum value of R is determined by In a trapped-ion context the condition (31) may be imposed by adjusting the distance of the traps for fixed k L and k R . It is also remarkable that when Eq. (31) is satisfied, the rectification no longer depends on the spring constants of the model. This last result can be found assuming Eq. (31) when calculating the currents with Eq. (19) and R with Eq. (30), where a and g are the mass and friction coefficients ratios The maximal rectification found does not scale with the magnitude of the masses or the friction coefficients, just with their ratios. Besides a high R, it is important to have non-vanishing heat currents [28]. Using again Eq. (31) in the expression for the currents (19), the maximum current J max = max( J , J ) is Now we analyze how the parameters a and g affect the maximum current J max in (34). To do this, we can divide the ag plane in four quadrants by the axes a = 1 and g = 1 (in those axes R = 0). In Eq. (34) the parameter a appears only in the denominator, thus for a higher a, a smaller current is found. The quadrants with a < 1 will be better for achieving large currents. However, g appears both in the numerator and denominator so there is no obvious advantageous quadrant for this parameter. Equation (32) is symmetric upon the transformations a ↔ 1/a and g ↔ 1/g. Using a logarithmic scale for a and g, the resulting R map will be symmetric with respect to the a = 1 and g = 1 axes. We can limit ourselves to analyze the quadrant a > 1, g > 1, as the results in other quadrants will be equivalent upon transformations a ↔ 1/a and g ↔ 1/g. Fig. 3 shows the rectification given by Eq. (32) in terms of a and g. Along any diagonal line (parallel to the solid cyan or the dashed green lines), the maximum value is at the center, that is, when a = g. However, if we fix a, increasing g always increases R. Although we could increase g arbitrarily to get more rectification this is not a realistic option in a trapped-ion set-up. Since g is defined as the ratio between the friction coefficients, increasing it means making either γ L go to 0 or γ R to infinity. Making γ L go to 0 decouples one of the ions from the bath, so the heat current tends to vanish in any direction. Also, increasing γ R arbitrarily is impossible since the Doppler cooling friction coefficient as a function of the laser detuning (Eq. (29)) is bounded. Although Eq. (29) suggests that boosting the laser intensity can also increase the friction coefficient, this is not an option since Eq. (29) is just an approximation for low laser intensities. When going to higher intensities, the emission/absorption of photons by the ion is saturated and the friction coefficient reaches a finite value proportional to the width Γ of the excited state [33]. As a compromise between feasibility and high R, we set the ratio between the friction coefficients g to be equal to the mass ratio a. As shown in Fig. 3, along the solidcyan and dashed-green diagonal lines the maximum R is achieved for a = g. Fig. 4 shows the rectification in Eq. (32) for the line a = g. When both parameters are large enough, the rectification goes to 1.

B. Spectral match/mismatch approach to rectification
The match/mismatch between the power spectra of the particles controls the heat currents in the system [2,22]. A good match between the power spectra of the two ions in a large range of frequencies yields a higher heat current through the system while the mismatch reduces the heat current. If there is a good match between the spectra of the ions (i.e., their peaks overlap in a broad range of frequencies) for a certain baths configuration, and mismatch when the baths exchange, the system will present heat rectification.
We have studied the phonon spectra of our model for several sets of parameters exhibiting no rectification or strong rectification. The phonon spectra of the ions is calculated through the spectral density matrix. For a real-valued stochastic process − → x (t), its spectral density matrix is defined as [24] with − → X (ω) being the Fourier transform of − → x (t) (we are using the convention of multiplying by a factor of 1 and 1 2π for the transform and its inverse operation). A justification of the use of the spectral density matrix to understand heat transport arises from the Wiener-Khinchin theorem [24], which says that the correlation matrix of a stationary stochastic process in the steady state is the inverse Fourier transform of its spectral density matrix Eq. (36) directly connects the spectral density matrix to the steady-state temperature and, therefore, to the heat currents (in Section III we saw that T s.s.

2
= m 2 C s.s. 4,4 /k B ). For the vector process − → r (t) describing the evolution being the Fourier transform of the white noise − → ξ (t). Note that − → Ξ (ω) does not strictly exist, because it is not square-integrable, however its spectral density is S− → ξ (ω) = 2D [24], which is flat as expected for a white noise. Therefore, the spectral density matrix of the system is As we can see in Eq. (37), the imaginary part of the eigenvalues of the dynamical matrix A correspond to the peaks in the spectrum whereas the real part dictates their width. The spectral density matrix of our model is where P A (λ) is the characteristic polynomial of the dynamical matrix A and S L (ω), S R (ω) are the matrix polynomials in the angular frequency ω whose coefficients are defined in Appendix B. Equation (39) gives the full expressions of the spectral densities for the velocities, S 3,3 (ω) = R 3 (ω)R 3 (−ω) for the left ion, and S 4,4 (ω) = R 4 (ω)R 4 (−ω) for the right ion, since they are the elements related to the calculation of the heat current using Eq. (36), , .
(39) Figure 5 depicts a series of plots of the spectra given by Eq. (39) that correspond to two points in Fig. 4. For c = 1 (Fig. 5(a) and (b)) there is no rectification, since the spectra match in the forward (a) and reversed (b) configurations. However, for c = 10 (( Fig. 5(c) and (d))) the picture is very different: there is a good match between the spectra in the forward configuration whereas in the reversed configuration the spectra are less correlated, giving as a result higher rectification (R ≈ 0.8). Figure  5 only shows the elements (3,3) and (4,4) in the diagonal of S but the remaining elements, including off-diagonal ones, exhibit a similar behavior.

VII. CONCLUSIONS
We have studied heat rectification in a model composed of two coupled harmonic oscillators connected to baths. This simple model allows analytical treatment but still has enough complexity to examine different ingredients that can produce rectification. Our results demonstrate in a simple but realistic system that harmonic systems can rectificate heat current if they have features which depend on the temperature [23]. We implement this no-tion of temperature-dependent features by defining the baths exchange operation as an exchange of both temperatures and coupling parameters of the baths to the system. This kind of temperature-dependent features happens naturally in laser-cooled trapped ion set-ups.
We have also studied the phonon spectra of the system, comparing the match/mismatch of the phonon bands, to reach the conclusion that the band match/mismatch description for heat rectification is also valid for systems which are harmonic, as long as there are temperaturedependent features. We hope this article sheds more light into the topic of heat rectification and that encourages more research regarding its physical implementation on chains of trapped ions.
Here we present the full set of equations for the covariance matrix elements in the steady state,  In Section VI we used the characteristic polynomial P A (λ) of the dynamical matrix A for the calculation of the spectral density matrix. P A (λ) is defined as det(A − λ) = λ 4 + λ 3 γ L m 1 + γ R m 2 + λ 2 (γ L γ R + m 2 (k + k L ) + m 1 (k + k R )) m 1 m 2 + λ (γ R (k + k L ) + γ L (k + k R )) m 1 m 2 + k(k L + k R ) + k L k R m 1 m 2 . (B1) We also used the polynomials S L (λ) and S R (λ), which are defined as S L (λ) = 6 n=0 λ n s L,n and S R (λ) = 6 n=0 λ n s R,n . There are 14 different polynomial coefficients, which are 4 × 4 matrices, which makes very cumbersome to include them in the main text. This is the full list of coefficients,