Thermoelectric efficiency in momentum-conserving systems

We show that for a two-dimensional gas of elastically interacting particles the thermoelectric efficiency reaches the Carnot efficiency in the thermodynamic limit. Numerical simulations, by means of the multi-particle collision dynamics method, show that this result is robust under perturbations. That is, the thermoelectric figure of merit remains large when momentum conservation is broken by weak noise.


Introduction
Understanding and controlling the behaviour of out-of-equilibrium systems is one of the major challenges of modern statistical mechanics. From a fundamental point of view, the challenge is to understand the origin of macroscopic transport phenomenological laws, such as diffusion equations, in terms of the properties of microscopic dynamics, typically nonlinear and chaotic [1,2]. The problem is extremely complex for coupled flows, so far barely studied from the viewpoint of statistical mechanics and dynamical systems [3,4,5,6,7,8,9,10,11,12]. In particular, it is of primary importance for thermoelectric transport [13,14,15,16] to gain a deeper understanding of the microscopic mechanisms leading to a large thermoelectric efficiency; see Ref. [17] for a review on the fundamental aspects of heat to work conversion.
Within linear response, and for time-reversal symmetric systems ‡ both the maximum thermoelectric efficiency and the efficiency at the maximum output power [22,23,24,25,26] are monotonous growing functions of the so-called figure of merit ZT = (σS 2 /κ)T , which is a dimensionless combination of the main transport coefficients of a material, that is, the electric conductivity σ, the thermal conductivity κ and the thermopower (Seebeck coefficient) S, and of the absolute temperature T . The maximum efficiency reads η max = η C √ ZT +1−1 √ ZT +1+1 , where η C is the Carnot efficiency, while the efficiency at maximum output power P max is given by η(P max ) = η C 2 ZT ZT +2 . Thermodynamics only imposes ZT ≥ 0 and η max → η C , η(P max ) → η C 2 when ZT → ∞. Since the different transport coefficients are interdependent, it is very difficult to find microscopic mechanisms which could provide insights to design materials with large ZT . While for non-interacting models it is well understood that energy filtering [27,28,29] allows us to reach the Carnot efficiency, very little is known for interacting systems [30]. It has been recently shown [31] that the thermoelectric figure of merit ZT diverges in the thermodynamic limit for systems with a single relevant conserved quantity, an important example being that of momentum-conserving systems, with total momentum being the only relevant constant of motion. While the mechanism is generic, it has been illustrated in Ref. [31] only for a toy model, i.e., a diatomic chain of hard-point elastically colliding particles.
In this paper, we show by means of extensive multi-particle collision dynamics simulations that the momentum-conservation mechanism leads to the Carnot efficiency in the thermodynamic limit also in the more realistic case of two-dimensional elastically colliding particles. Furthermore, we show that this mechanism leads to a significant enhancement of the thermoelectric figure of merit even when the momentum conservation is not exact due to the existence of an external noise. This robustness is particularly relevant in experiments for which inelastic or incoherent processes are unavoidable to some extent. In this case, the figure of merit saturates with the size of the system to a value higher, the weaker is the noise. Finally, we discuss the validity range of linear response.
The paper is organised as follows. In Sec. 2, in order to make the paper selfcontained, we review the theoretical argument of Ref. [31] explaining the divergence of the thermoelectric figure of merit ZT in the thermodynamic limit for systems with a single relevant constant of motion. In Sec. 3 we explain our out-of equilibrium multiparticle collision dynamics simulations. Our numerical results are presented in Sec. 4. We finish with concluding remarks in Sec. 5. ‡ Thermodynamic bounds on efficiency for systems with broken time-reversal symmetry are discussed in [18,19,20,21] 2. Theoretical argument

Linear response irreversible thermodynamics
The equations connecting fluxes and thermodynamic forces within linear irreversible thermodynamics read as follows [32,33]: where J ρ and J u are the particle and energy currents, µ the chemical potential and β = 1/T the inverse temperature (we set the Boltzmann constant k B = 1). The kinetic coefficients L ij ( i, j = {ρ, u}), are related to the familiar transport coefficients as where L denotes the (Onsager) matrix of kinetic coefficients and we have set the electric charge of each particle e = 1. Thermodynamics imposes det L ≥ 0, L ρρ ≥ 0, L uu ≥ 0; L uρ = L ρu follows from the Onsager reciprocity relations. The thermoelectric figure of merit reads Furthermore, the Green-Kubo formula expresses the kinetic coefficients in terms of correlation functions of the corresponding current operators, calculated at thermodynamic equilibrium [34,35]: where · = tr ( · ) exp −βH /tr [exp(−βH)] denotes the equilibrium expectation value at temperature T and Ω is the system's volume. Within the framework of Kubo's linear response approach, the real part of L ij (ω) can be decomposed into a singular contribution at zero frequency and a regular part L reg ij (ω) as The coefficient of the singular part defines the generalized Drude weights D ij §, which can be expressed as where in the volume Ω(l) we have explicitly written the dependence on the system size l along the direction of the thermodynamic flows. Non-zero Drude weights, a signature of ballistic transport [37,38,39,40], namely in the thermodynamic limit the kinetic coefficients L ij scale linearly with the system size l. As a consequence, the thermopower S does not scale with l.

Conservation laws
We now discuss the influence of conserved quantities on the figure of merit ZT . Making use of Suzuki's formula [41] for the currents J ρ and J u , one can generalize Mazur's inequality [42] by stating that, for a system of finite size l (along the direction of the flows), where for readability, in the right hand side of the equation we have omitted the dependence on l. The summation in Eq. (2.8) extends over all the M constants of motion Q n , which are orthogonal, Q n Q m = Q 2 n δ n,m , and relevant for the considered flows. That is, J ρ Q n = 0 and J u Q n = 0.
From Eq. (2.8) one can define the finite-size generalized Drude weights as Therefore, the presence of relevant conservation laws directly implies that the finite-size generalized Drude weights are different from zero. If the thermodynamic limit l → ∞ can be taken after the long-time limit t → ∞, so that the generalized Drude coefficients can be written as and moreover D ij = 0, then we can conclude that the presence of relevant conservation laws yield non-zero generalized Drude weights, which in turn imply that transport is ballistic. We point out that, in contrast to Eq. (2.10), one should take the thermodynamic limit l → ∞ before the long-time limit t → ∞. While it remains an interesting open problem for which classes of models the two limits commute ¶, numerical evidence suggests that it is possible to commute the limits for the models considered in Ref. [31] and in the present paper. Let us first consider the case in which there is a single relevant constant of motion, M = 1. We can see from Suzuki's formula, Eq. (2.8), that the ballistic contribution to det L vanishes, since it is proportional to D ρρ D uu − D 2 ρu , which is zero from (2.8) and (2.10). Hence, det L grows only due to the contributions involving the regular part in Eq. (2.6), i.e., slower than l 2 , which in turn imply that the thermal conductivity κ ∼ det L/L ρρ grows sub-ballistically. Furthermore, since σ ∼ L ρρ is ballistic and S ∼ l 0 , we can conclude that for a proof of the commutation of the two limits for a class of quantum spin chains.
Thus ZT diverges in the thermodynamic limit l → ∞.
The situation is drastically different if M > 1, as it would be the case for integrable systems, where typically the number of orthogonal relevant constants of motion equals the number of degrees of freedom. In that case, due to the Schwartz inequality, The equality arises only in the exceptional case when the vectors x ρ and x u are parallel. Hence, for M > 1 we expect, in general, det L ∝ l 2 , so that heat transport is ballistic and ZT ∼ l 0 .

Momentum-conserving gas of interacting particles
In this section we analyse the consequences of our analytical results in a two-dimensional gas of interacting particles. We consider a gas of point-wise particles in a rectangular two-dimensional box of length l and width w. The gas container is placed in contact with two particle reservoirs at x = 0 and x = l, through openings of the same size as the width w of the box. In the transversal direction the particles are subject to periodic boundary conditions.
The dynamics of the particles in the system are solved by the method of Multiparticle Collision Dynamics (MPC) [43], introduced as a stochastic model to study solvent dynamics. The MPC simplifies the numerical simulation of interacting particles by coarse graining the time and space at which interactions occurs. MPC correctly captures the hydrodynamic equations [44,45]. It has been successfully applied to model steady shear flow situations in colloids [46], polymers [47], vesicles in shear flow [48], colloidal rods [49], and more recently to study the steady-state of a gas of particles in a temperature gradient [50].
Under MPC dynamics the system evolves in discrete time steps, consisting on free propagation during a time τ , followed by collision events. During propagation, the coordinates r i of each particle are updated as where v i is the particle's velocity. For the collisions the system's volume is partitioned in identical cells of linear size a. Then, the velocities of the N particles found in the same cell are rotated with respect to the center of mass velocity by a random angle. In two dimensions, rotations by an angle +α or −α with equal probability p(+α) = p(−α) = 1/2 are performed. The velocity updating after a collision event reads i=1 v i is the center of mass velocity andR θ is the two-dimensional rotation operator of angle θ. Furthermore, to guarantee Galilean invariance, the collision grid is shifted randomly before each collision step. It has been shown that for these dynamics, the equation of state of the gas of particles corresponds to that of an ideal gas [43]. Moreover, the time interval between successive collisions τ and the collision angle α tune the strength of the interactions and consequently affects the transport coefficients of the gas of particles. When α is a multiple of 2π, the particles do not interact, propagating ballistically from one reservoir to the other as they cross the system. For any other value of α, the particles interact, exchanging momentum during the collision events. The value α = π/2 corresponds to the most efficient mixing of the particle momenta. Note that by construction, the collision preserve the total energy and total momentum of the gas of particles.
From the reservoir k (k = L, R for the left and the right reservoir), particles of mass m enter the system at rate γ k obtained by integration of the appropriate canonical distribution to give where ρ k and T k are the particle density temperature. Assuming that the particles in the reservoirs behave as ideal gas, the particle injection rate is related to the value of the chemical potential µ k of the reservoir k as with µ 0 an arbitrary constant whose value does not qualitatively modify the results discussed in this paper; hereafter we set µ 0 in such a way that µ = 0 + . Whenever a particle from the system crosses the boundary which separates the system from reservoir k, it is removed (absorbed in the reservoir), i.e., it has no further effects on the evolution of the system.

Discussion of numerical results
We have numerically studied the nonequilibrium transport of the model defined in Sec. 3, coupled to two ideal particle reservoirs. The nonequilibrium state is imposed by setting the values of T and µ/T in the reservoirs to different values, meaning that from each of the reservoirs, the particles are injected into the system at different rates and with a different distribution of their velocities. Out of equilibrium the kinetic coefficients L ij can be computed, in the linear response regime, by direct measurement of the particle and energy currents in the system. Using (2.1), it is enough to perform two nonequilibrium numerical simulations: one with T L = T R and µ L /T L = µ R /T R , and one with T L = T R and µ L /T L = µ R /T R . In the first simulation the reservoirs' temperatures + This arbitrariness is intrinsic in classical mechanics and can only be removed by means of semiclassical arguments, see Ref. [10].
are set to T L = T − ∆T /2 and T R = T + ∆T /2, so that the temperature gradient is given by ∆T /l, while µ L /T L = µ R /T R . Conversely, in the second simulation we set T L = T R = T and using (3.4), we set the particle injection rates γ L and γ R so that In all simulations the mean particle density and mean temperature in the reservoirs was set to n = N/lw = 22.75 (N is the mean number of particles) and T = 1, respectively. We parametrize the gradients in terms of a single parameter by setting, ∆T = ∆ (µ/T ) ≡ ∆ (in units where k B = e = 1). The rotation angle for the collisions in the MPC scheme was set to α = π/2, unless otherwise specified. The length of the collision cells in the MPC scheme was set to a = 0.1 and the time step to τ = 0.25. For these values and small ∆ the system exhibits reasonably linear temperature and chemical potential profiles in the bulk, with some nonlinear boundary layer near the contacts, arising from the fact that the mean free path of the particles near the boundaries is different than in the bulk * , yielding a contact resistance [50]. We performed numerical simulations with up to Ω = 10 3 (l = 500 and w = 2), so that systems with mean number of particles up to N = 4.55 × 10 4 were considered.

Linear response transport
Using the Suzuki's formula (2.8), the current-current correlation functions C ij (t) can be obtained analytically. The particle current is J ρ = N i=1 v x,i and the energy current The other constants of motion, i.e. momentum in the transverse direction, energy and number of particles, are irrelevant since they are orthogonal to the considered flows. Therefore, in this case M = 1.
Applying Eq. (2.8) and integrating over the equilibrium state at temperature T and fixed number of particles N, we obtain that the finite-size correlators are To verify Eq. (2.8) we have numerically computed the equilibrium current-current time correlation functions for the isolated system, averaged over an equilibrium ensemble of initial conditions with N = 1000 particles of mass m = 1 and temperature T = 1. A square container of size l = 2 and periodic boundary conditions in both directions was considered. The results, shown in Fig. 1, verify our the analytical expressions. Note that the initial values C ρρ (0) and C ρu (0) of the time-averaged correlation functions C ρρ (t) and C ρu (t) are equal to their asymptotic values C ρρ (l) = lim t→∞ C ρρ (t) and C ρu (l) = lim t→∞ C ρu (t). On the other hand, it is easy to compute analytically * The MPC collisions at the boundaries are implemented without taking into account the particles in C uu (0) = 6NT 3 /m, and numerical data show that C uu (t) converges algebraically to its asymptotic value C uu (l) = lim t→∞ C uu (t) = 4NT 3 /m. This asymptotic behaviour may be due to the slow decay of the energy hydrodynamic modes. Equation (4.1) also shows that the dependence of the correlations on the size l comes exclusively through the number of particles N. The thermodynamic limit requires keeping the density of particles fixed, so that the number of particles has to scale linearly with the volume of the system: N ∝ Ω(l) = lw. Therefore, Eqs. (4.1) imply that the finite-size generalized Drude weights of Eq. (2.9) do not scale with l. Using Eqs. (2.8)-(2.10) we obtain, for fixed w, the generalized Drude weights As a consequence of the finiteness of the Drude weights, the transport is ballistic, meaning that all kinetic coefficients L ij scale linearly with the size of the system: L ij ∼ l. This prediction is confirmed by the numerical results shown in panel a of Fig. 2.
More importantly, as discussed in Sec. 2.2, due to the conservation of total momentum, the ballistic contribution to the determinant of the Onsager matrix is zero. Indeed, it can be readily seen from Eq. (4.2) that D ρρ D uu − D 2 ρu = 0. Hence a scaling det(L) slower than l 2 is expected. From the nonequilibrium numerical simulations the scaling of the determinant with l is consistent with det(L) ≈ l 1.15 (dotted-dashed curve in Fig. 2-b). It is worthwhile recalling that different analytical methods such as mode coupling theory and hydrodynamics predict, for momentum conserving systems in two dimensions, a logarithmic divergence of the thermal conductivity with the size of the system [1,2]. Therefore, one should expect that det(L) ∼ l log(l). We show in Fig. 2 (dashed curve), that such scaling is also consistent with our numerical results, though deviations are larger than for the algebraic behaviour at small system sizes. Since we have no reason to expect an algebraic sub-ballistic behaviour of the heat conductivity, we will assume in what follows that its behaviour is logarithmic.

Strong enhancement of ZT
From Eqs. (2.2) and Eqs. (4.2) we obtain that the electric conductivity also scales linearly with the size of the system: with A constant. The dependence on l of the Seebeck coefficient cancels out to give, asymptotically in l, Since the ballistic contribution to det(L) vanishes, i.e. D ρρ D uu − D 2 ρu = 0, we cannot derive an explicit expression for the heat conductivity κ. However, as discussed in the previous section, for momentum conserving two-dimensional systems it is predicted that κ diverges logarithmically with respect to the size of the system: κ ∼ log(l).  Eq. (4.4) (asymptotically in l) only in the limit of small forces (in Fig. 3, S is shown for µ = 0). We have found that S converges to the value S = 2 predicted by (4.4) as 1/∆.
The heat conductivity κ exhibits the logarithmic behaviour up to a size l = l ⋆ dependent on the strength ∆ of the thermodynamic forces. For any value of ∆, the heat conductivity grows as κ ∼ log(l), for l > l ⋆ . The smaller the ∆ the larger the range of validity of the logarithmic κ is. We have obtained numerically that the characteristic length l ⋆ grows linearly with 1/∆. Through Eq. (2.3), this characteristic length l ⋆ does also determine the behaviour of the figure of merit ZT . In Fig. 4 we show ZT as a function of l, for different values of ∆. We observe that for any value of ∆, ZT is in reasonable agreement with an initial grow l/ log(l) for l < l ⋆ and for larger sizes saturates to a maximum value (ZT ) max . Our results show that as a consequence of the existence of a single relevant conserved quantity, the values of ZT are greatly enhanced when the system under consideration is large enough. Moreover, ZT does not grow unboundedly, but reaches a maximum value that grows with ≈ (1/∆) 0.9 (see the inset of Fig. 4). The deviations at short sizes are probably due to the slow convergence of the Seebeck coefficient to its asymptotic value 2.
In the above discussion on the behaviour of the transport coefficients and the figure of merit ZT as a function of ∆, we should keep in mind that such coefficients and consequently also ZT are defined in the linear response regime, i.e. in the limit of small thermodynamic forces, formally for ∆ → 0. On the other hand, we numerically computed the kinetic coefficients, for any given ∆, via the fluxes as discussed at the beginning of Sec. 4. That is to say, there is no saturation of ZT within linear response. On the other hand, the numerically observed saturation (as well as the ballistic behaviour of κ for l > l ⋆ ) signals that the range of linear response shrinks with the system size when computing κ and ZT . At first sight, this failure of linear response for a given ∆ and large l appears counterintuitive, since for fixed ∆ larger l means smaller thermodynamic forces, and it is in the limit of small forces that linear response is expected to be valid. There is actually no such problem when computing the kinetic coefficients L ij . As shown in Fig. 3 for the charge conductivity σ = L ρρ /T , data at different ∆ collapse on a single curve, showing that for all values of ∆ in that figure we are within linear response. The problem arises when considering non-trivial combinations of the kinetic coefficients, as in κ ∝ det(L) and consequently in ZT . Our theory predicts the divergence of ZT in the thermodynamic limit and ZT diverges (thus leading to Carnot efficiency) if and only if the Onsager matrix L becomes illconditioned, namely the condition number [Tr(L)] 2 / det(L) diverges (in our model as l/ log(l)) and therefore the system (2.1) becomes singular. That is, the charge and energy currents become proportional, a condition commonly referred to as strong coupling, i.e. J ρ = cJ u , the proportionality factor c being independent of the applied thermodynamic forces. The Carnot efficiency is obtained in such singular limit and it is in attaining such limit that the validity range of linear response shrinks. Therefore, as expected on general grounds, the Carnot efficiency is obtained only in the limit of zero forces and zero currents, corresponding to reversible transport (zero entropy production) and zero output power. It is worthwhile noticing that for our model in the non-interacting limit the momentum of each particle is conserved, meaning that the system is integrable and the number of conserved observables M ∝ l, thus diverging in the thermodynamic limit. As we have discussed at the end of section 2.2, one expects that in such integrable situation, ZT does not scale with the system size. To corroborate this expectation we show in Fig. 5 the dependence of ZT on l for different values of the collisional parameter α. We recall that at the collisions, α = π/2 corresponds to the most efficient mixing of the particle momenta, while α = 0 corresponds to no interaction. As expected, for the non-interacting gas, namely for an infinite number of conserved quantities, ZT does not scale with l, attaining the value 3/2 characteristic of a two-dimensional ideal gas [7]. The enhancement of ZT is observed for any value of α > 0, as then only the total momentum is preserved and M = 1. Our data also suggest a rather weak dependence of (ZT ) max on α.

Systems with noise
The results discussed above show the enhancement of ZT , and thus of the thermoelectric efficiency, in systems with conserved total momentum. In real systems, however, total momentum is never conserved due to the phonon field, the presence of impurities or in general to inelastic scattering events.
In this section we want to explore to what extent the break down of total moment conservation modifies the results obtained above. To address this question numerically, we consider the existence of a source of stochastic noise. From a physical point of view, this noise source may model the interactions of the gas with the walls of the container, or the inelastic scattering from impurities in the material. We model the stochastic noise as follows: after a collision of the particles in a given cell has taken place, with probability ε the velocities of all the particles in the cell are reflected, namely v i → − v i . Therefore for any ε > 0 the total momentum is not longer conserved. If ε is small the momentum conservation is weakly broken and we want to investigate how our results depend on the strength ε of the perturbation.
In Fig. 6 we show the dependence of the transport coefficients on l for fixed ∆ and different strengths of the noise ε. We observe that for sufficiently strong noise all transport coefficients appear to become independent of l, as expected in a diffusive regime in which total momentum is not preserved.
More interesting is the behaviour of ZT shown in Fig. 7. We see that at stronger noise, ZT becomes constant, as expected in the diffusive regime. From a mathematical point of view, the absence of conserved quantities (M = 0) leads to decaying correlation functions and zero Drude coefficients (inset of Fig. 7). Thus the transport coefficients and ZT become size-independent.
More importantly, we see that when the convergence toward the diffusive regime is smooth, meaning that when the conservation of total momentum is only weakly perturbed (small ε), the enhancement of ZT can be significant. This shows that the effect described here is robust against perturbations.

Conclusions
In summary, we have shown that in two-dimensional interacting systems, with the interactions modeled by the multi-particle collision dynamics method, the thermoelectric figure of merit diverges at the thermodynamic limit. In such limit, the Carnot efficiency is obtained with zero output power. When noise is added to the system, ZT saturates at large l, to values higher the weaker is the noise strength.  Figure-of-merit ZT as a function of the length of the system l, for different noise intensities ε = 0 (squares), 0.01 (triangles) and 0.1 (circles). The thermodynamic gradient was fixed to ∆ = 0.0125. The dashed curve corresponds to the expected linear response dependence at ε = 0, i.e. ZT ∼ l/ log(l). In the inset: energy current-current correlation for different noise intensities ε, for the same parameter values as in Fig. 6. From top to bottom: ε = 0, 10 −4 , 10 −3 , 10 −2 , 10 −1 . The dashed curve stands for ∼ 1/t.
Our findings could be relevant in situations in which the elastic mean free path is longer than the length scale over which interactions are effective in exchanging momenta between the particles. Suitable conditions to observe the interaction-induced enhancement of the thermoelectric figure of merit might be found in high-mobility twodimensional electron gases at low temperatures. In such systems very large elastic mean free paths have been reported (for instance, up to 28 µm in Ref. [51]). At low temperatures the inelastic mean free path is determined by electron-electron interactions rather than by phonons. It should be therefore possible to find a temperature window where electron-electron interactions dominate, i.e. are effective on a scale smaller than the elastic mean free path and are dominant over phonon effects. It would be, however, highly desirable to test our arguments in such regime, by means of numerical simulations of quantum systems.