Quantum heat engine with long-range advantages

The employment of long-range interactions in quantum devices provides a promising route towards enhancing their performance in quantum technology applications. Here, the presence of long-range interactions is shown to enhance the performances of a quantum heat engine featuring a many-body working substance. We focus on the paradigmatic example of a Kitaev chain undergoing a quantum Otto cycle and show that a substantial thermodynamic advantage may be achieved as the range of the interactions among its constituents increases. Interestingly, such an advantage is most significant for the realistic situation of a finite time cycle: the presence of long-range interactions reduces the non-adiabatic energy losses, by suppressing the detrimental effects of dynamically generated excitations. This effect allows mitigating the trade-off between power and efficiency, paving the way for a wide range of experimental and technological applications.


I. INTRODUCTION
Thermodynamics was born during the industrial revolution in response to the very practical necessity to understand and master the functioning of thermal machines such as heat engines and refrigerators. This has led to the development of thermal machines of increasingly high performance in terms of thermodynamic efficiency and power consumption. Such devices permeate our daily life and are as well crucial for technology advances. Few out of countless examples span from fuel-based vehicles, household air conditioners and refrigerators to the dilution refrigerators for quantum applications [1]. Fault-tolerant quantum computing [2] represents the latest technological challenge, which in recent years has attracted significant research efforts [3][4][5][6][7][8][9]. In particular the route towards fault tolerance of available quantum processors is dictated by the quantum threshold theorem [10][11][12], which states that it is possible to correct errors, even by using noisy gates, provided that the noise level remains below a certain threshold.
Cooling quantum hardwares at sufficiently low temperatures allows, in principle, to achieve this threshold, but the introduction of large classical apparata, i.e. thermal baths, may generate additional sources of decoherence. Thus, the design of microscopic and coherent thermodynamic machines constitutes an urgent technological challenge [13], which has led researchers to study quantum thermal engines, i.e., heat engines and refrigerators whose working substance operates directly in the quantum domain [14]. Much theoretical [15][16][17][18][19][20][21][22] and experimental [23][24][25][26][27][28][29][30] efforts have been devoted to their study, recently showing their potential applicability to existing quantum processors [31,32]. However these devices suffer from the well-known trade off between power and efficiency [22]. This is due to two main contributions: the first one comes from the fundamental limitations imposed by the second law of thermodynamics on irreversible processes, which implies that the thermodynamic efficiency of a heat engine has to be smaller than that of a Carnot engine [33]. Furthermore, any realistic cycle is carried out in a finite time. This limitation is an additional source of losses, due to the dynamic generation of excitations in the system which dissipate energy, worsening the performance of the device. Consequently, an increase in power typically involves a greater dissipation, with detrimental effects on efficiency. Various methods that suppress the non-adiabatic transitions, based on the so-called shortcuts to adiabaticity, have been proposed to overcome this problem [34][35][36][37]. However these techniques typically require switching on additional driving fields, a procedure whose additional energetic cost reduces the actual work output of the device [38,39].
Here, we propose a new paradigm to reduce these detrimental effects, which consists in using, as the working substance of the engine, a long-range interacting quantum system, i.e., a system in which the coupling energy between two of its microscopic constituents J i,j decays as a power law of their distance r = |i−j|: J i,j ∝ r −α , with α > 0 [40,41].
Long-range interacting systems are emerging as promising platforms for quantum technological applications, due to their stability against external perturbations, which allows keeping the impact of dynamically generated excitations under control, therefore mitigating their detrimental effects [42]. An example of the rigidity of long-range interacting platforms against external drivings, is the possibility for such systems to host clean discrete Floquet time crystal phases [43][44][45][46], i.e., robust out-of-equilibrium symmetry broken phases in systems subjected to an external timeperiodic driving [47]. Moreover, long-range interactions are known to alter the universal behaviour of critical systems at and out of equilibrium, generating a wide range of unprecedented phenomena. These include novel form of dynamical phase transitions [48,49], defect formation [48,[50][51][52], anomalous thermalization [53], information spreading [54][55][56] and metastable phases [57,58], whose features have no-counterpart in systems with short-range interactions.
In this work we show how the peculiar properties of long-range interacting quantum systems may be used to boost the performances of a many-body quantum Otto cycle [21,22], which we take as a prototypical example of a quantum thermal device. In particular, we identify multiple advantages stemming from the presence of long-range interactions 1. In the limit of an infinitely slow cycle, enhanced optimal performances of the device are observed, confirming previously reported results [59]. In particular, the enhancement is observed in the engine most useful operation modes: the heat-engine mode and the refrigerator mode.
2. In the realistic case of a finite time cycle featuring the crossing of a quantum critical point, the choice of a longrange-interacting working substance leads to a significant reduction in nonadiabatic energy losses and defect proliferation if compared with its short-range counterpart.
From the experimental point of view, long-range interacting systems may be implemented in typical quantum simulation platforms, such as atomic molecular and optical (AMO) systems [60][61][62][63][64]. Interestingly, trapped ions setups allow to tune the power law exponent α, dictating the decay of the interaction energy with distance, from α 0 to α 3 [60]. Moreover, such platforms have been proven useful in realising quantum computers [65][66][67][68][69][70]. Finally, single atom quantum heat engines have been already implemented in trapped ions setups [24]. Accordingly, the implementation of thermodynamic devices featuring long-range interacting working substances may be feasible in current trapped ion platforms, with possible direct applications to the cooling of the existing ion-based quantum processors [70].
The paper is organized as follows. In Sec. II we introduce the many-body Otto cycle, taking as a prototypical example of a long-range working substance, the long-range version of the celebrated Kitaev chain [71]. In Sec. III we focus on the limit of the infinite-time, adiabatic, Otto cycle, observing a long-range advantage in the optimal regime. Finally, in Sec. IV we estimate the energy losses of more realistic finite-time cycles with respect to the adiabatic case, showing that they are reduced in the critical regime thanks to the presence of long-range interactions.

II. MANY-BODY QUANTUM OTTO CYCLE
A. Many body working substance As a prototypical example of a many-body working substance with, possibly, long-range interactions, we consider a generic model of spinless fermions hopping across the N sites of a linear chain in the presence of pairing interaction, and with a chemical potential h. Assuming periodic boundary conditions, the system Hamiltonian reads where c † j and c j are creation and annihilation operators for fermions at site j, while t r and ∆ r are the hopping and pairing amplitudes, respectively. We choose their dependence on the intersite distance r according to the power laws with the hopping exponent α 1 > 0, the pairing exponent α 2 > 0, and N α = N/2 r=1 r −α the Kac scaling factor [72], which guarantees extensivity of the energy in the case α i < 1, with i = 1, 2. Hereafter, we set J = ∆ = 1 as the energy scale and work in units of = k B = 1. This model, often referred to as long-range Kitaev chain [71,73,74], has been mainly studied in two particular cases: with short-range hopping α 1 = ∞ and generic long-range pairing α 2 [73,75,76], and with same range for pairing and hopping α 1 = α 2 = α [41,57]. As observed in Refs. [77][78][79], the latter case can be related to the quantum Ising model. In particular, in the short-range case with α → ∞, the relation becomes exact through the Jordan-Wigner mapping [80]. In this Section, we review the main properties of the model with general power law couplings as in Eq. (1). For reasons that will become clear later, we restrict ourselves to the regime α 1 , α 2 > 1.
Due to the translational invariant nature of the couplings, it is useful to write the Hamiltonian in terms of the momentum-space operatorsc where k = 2πn/N , with n = −N/2 + 1, . . . , N/2 (In the following we will drop the˜on the c k unless it is ambiguous).
Then we obtain where t k and ∆ k are the Fourier transforms of the hopping and pairing amplitudes, respectively, which in the thermodynamic limit may be written as where Li α (z) denotes the polylogarithm and ζ(α) is the Riemann zeta function. We notice that the Hamiltonian in the Fourier space can be decomposed into the sum of single mode Hamiltonians, where σ (a) k , a = x, y, z are the sigma Pauli operators. Let us notice how the k-th term of the Hamiltonian acts on a different sector of the total Hilbert space, namely the two dimensional subspace spanned by the states |0 k , 0 −k , Then the Hamiltonian is diagonalized via a Bogoliubov transformation, in terms of the fermionic quasiparticle operators γ k = u k c k + v * −k c † −k , with Bogoliubov coefficients where with the spectrum For α 1 , α 2 > 1, the system possesses two different phases separated by two quantum critical points h c = 1, −1 + 2 1−α1 , in correspondence of which the dispersion relation becomes gapless near to the critical mode k c = 0, π, respectively [41,81]. Figure 1 shows the single-particle spectrum ω k for different values of 0 ≤ k ≤ π, plotted as a function of the chemical potential h. In particular Fig. 1a) shows the spectrum for a short-range Kitaev chain, which, as previously stated, is exactly equivalent to the nearest neighbor quantum Ising chain, showing a ferromagnetic and an antiferromagnetic quantum critical point located at h c = 1, −1, respectively. A similar structure is present in the two relevant long-range cases displayed in Fig. 1b) where α 1 = ∞ and α 2 = 1.5, and in Fig. 1c) where α 1 = α 2 = 1.5. Although the mapping with the quantum Ising model is no longer exact for finite α 1 ,α 2 , the two critical points are still referred to as ferromagnetic and antiferromagnetic. However, we notice that the location of the antiferromagnetic critical point is α 1 -dependent since h c = −1 + 2 α1−1 [81]. Moreover, the spectrum dispersion relation in the proximity of the critical modes is affected by the presence of long-range interactions, thus leading to α-dependent critical exponents [41,81]. Further details on the dispersion relations at the quantum critical points and for different values of α 1 and α 2 are provided in Appendix B.

B. Description of the cycle
The quantum Otto cycle [21,22,82], consists of the following four strokes (see Fig. 2): • Initially the system is assumed to be in thermal equilibrium with the hot reservoir at temperature T 1 = 1/β 1 and h = h 1 . Then, in the first stroke (unitary decrease of h, 1 → 2) the system is decoupled from the bath and it undergoes a unitary evolution where the chemical potential changes from h 1 to h 2 .
• In the second stroke (thermalization at fixed h, 2 → 3) the chemical potential is kept fixed at h = h 2 , while the system is coupled with the thermal bath 2 so as to reach equilibrium at the temperature T 2 = 1/β 2 .
• In the third stroke (unitary increase of h, 3 → 4) the system undergoes another unitary driving which brings the chemical potential back to the initial value h 2 → h 1 .
• In the fourth and last stroke (thermalization at fixed h, 4 → 1) the system at fixed h = h 1 is again coupled with bath 1 so as to reach the equilibrium at temperature T 1 , thus closing the thermodynamic cycle.
At the beginning of each stroke of the cycle (points 1, 2, 3, 4 in Fig.2), the state of the system and the corresponding average energy are given by where H i = H(h i ), Z i = Tre −βiHi i = 1, 2, U andŨ are the unitary evolution operators associated to the first and the third stroke respectively. In the following, we shall assume a linear time dependence of the chemical potential during the unitary strokes. Then the driving protocol corresponding to the first step 1 → 2 can be written as where δ = (h 1 − h 2 )/τ is the sweep rate. The driving protocol during the third step of the cycle 3 → 4 is given bỹ The corresponding unitary evolutions then read where T exp denotes the time-ordered exponential. During the second and the fourth strokes the external driving is switched off and the system interacts only with the baths, reaching thermal equilibrium. While long-range interacting systems are known to evade thermalization allowing for quasistationary states [57,83], it can be shown that thermal equilibrium is safely reachable when α 1 , α 2 > 1. The behavior of the long-range Kitaev chain coupled to thermal reservoirs has been recently addressed in Ref. [84]. In particular, one can couple the system to N thermal baths, one for each lattice site, all characterized by the same temperature T . The internal baths dynamics are described by a continuous model of free fermions with Hamiltonian [85] which is coupled to the chain through the linear interaction Hamiltonian Assuming that the Born-Markov and the secular approximations [86] are satisfied, then a Linblad master equation can be derived [87] for the system density operator ρ, reading where H Ls is the Lamb-shift correction to the system Hamiltonian, and the dissipator D has the form [85] f (ω k ) = (1 + e βω k ) −1 being the Fermi-Dirac distribution, and {A, B} = AB + BA representing the anticommutator between the A and B operators, and J (ω) = π dq|g(q)| 2 δ(ω − (q)) the bath spectral density. The solution of the Lindblad master equation (18) shows that the populations of the various normal modes evolve independently from one another towards the steady state, each with a decay rate equal to r k = 2J (ω k ) [84] When the latter is nonzero for all k, the steady state is unique and is characterized by the thermal expectation values [87] γ This justifies the use of the canonical thermal equilibrium state ρ = e −βH /Z, with Z = Tr[e −βH ], in Eqs. (12a) and (12c). Then, all the thermodynamic quantities can be easily computed using Eq. (20) or, equivalently, directly computing the partition function Z. In particular the internal energy reads with Since during the second and the fourth strokes the external driving is switched off and the system interacts only with the baths, then energy is exchanged with them only in the form of heat According to the first law of thermodynamics we also have The three average energy exchanges Q 1 ,Q 2 and W completely characterize the cycle operation.

C. Operation modes
Depending on the signs of Q 1 , Q 2 and W , our engine may operate in any of the following four modes where  [88]. The Kitaev chain can be decomposed into a collection of N non-interacting qubits, each with level spacing ω k , (see Eqs. (II A) and (II A)), which act as independent quantum thermal machines. In fact, the energy exchanges of the many-body device can be written as where Q 1,k , Q 2,k and W k denote the energy flows corresponding to the k-th mode. Let us introduce the transition probabilities 1 − P k , i.e., the probability of a nonadiabatic transition among the levels of the k-th qubit during the unitary stroke of the cycle. Then the energy exchanges take the form where we have introduced the shortcut notation f i,k = tanh(β i ω i,k /2), for i = 1, 2. Let us notice that, as the system undergoes an Otto cycle, different qubits may operate in different regimes, among the [E], [R], [A], or [H] mode. As a consequence, the resulting operation mode of the many-body engine depends non-trivially on the interplay between the different qubit engines.

III. ADIABATIC CYCLE
Let us now analyze the case of an infinitely slow cycle i.e., the limit δ → 0 (τ → ∞). This regime is usually referred to as adiabatic, since the unitary evolution is sufficiently slow for the adiabatic theorem to hold, preventing transitions between the instantaneous eigenstates of the Hamiltonian, and leading to P k ∼ 1 − O(τ −2 ). The adiabatic approximation is known to break down as the energy gap closes. Strictly speaking, however, this happens only in the thermodynamic limit N → ∞. Accordingly, for any finite N , one can choose the driving time scale such that the adiabatic approximation is justified, allowing us to set P k = 1 in Eqs. (27) obtaining (28c) Figure 3 shows the regions of the parameters, β 2 and h 2 , for fixed values of β 1 and h 1 , corresponding to the different operation modes (25). These are compared with corresponding regions in the Einstein approximation, in which the chain is replaced by N identical single qubits with frequencȳ i.e., the average level spacing. The region boundaries obtained in this approximation correspond to the white lines in Fig. 3 and provide a rough estimation for the engine operation mode. Since, in this limit, only one level-spacing is present, one can obtain these boundaries by applying the results of Ref. [82] for the operation modes of a single qubit, namely whereω 1,2 corresponds to Eq. (29) for h = h 1,2 respectively and we have assumed, without loss of generality, β 1 ≤ β 2 .
Let us notice that, as conditions (30) rule out the possibility of a single qubit acting as a heater ([H]), this regime cannot be well described within the mean-spacing approximation in the adiabatic limit.
In the region, h 2 /h 1 > 0, however, in which the heater phase is not present, the operation mode phase diagram is well reproduced by the Einstein approximation, see Eq. (30), regardless of the values of α 1 and α 2 . By comparing the three diagrams in Fig. 3 we notice that the case α 1 = ∞, first studied in Ref. [59], is pretty similar to the nearestneighbours limit, and does not lead to substantial advantage. On the other hand, in the presence of both long-range hopping and long-range pairing (α 1 = α 2 = 1.5) the heater region is greatly reduced, and the Einstein approximation works better even in the h 2 /h 1 < 0 region.
[H]  In conclusion, the presence of long-range interactions may generate a substantial advantage provided proper values of the decay exponents α 1 and α 2 are chosen. Indeed, in the equally long-range case (Fig. 3c) the [R] and [E] regimes, which are the most relevant to technological applications, are enhanced with respect to the nearest neighbor case (Fig. 3a)), and become prevalent in the whole parameter region |h 2 /h 1 | < 1 . In our studies, the most profitable configuration of the Hamiltonian in Eq. (1) corresponds to the realistic case, where the Kitaev Hamiltonian approximately represents the long-range Ising model, i.e., α 1 = α 2 , which will be the focus of the following discussion.

A. Heat Engine operation mode
The purpose of a heat engine is to exploit natural flow of heat from a hot reservoir to a cold one to generate work. Thus, the performance of a device operating in the [E] mode may be optimized by maximizing the work output.  Another estimator of the engine performance is the heat engine efficiency, defined as the ratio between the energy gain in the form of work and the heat extracted from the hot reservoir The second law of thermodynamics imposes this efficiency to be always smaller than the Carnot efficiency η C [E] = 1 − T 2 /T 1 [33]. The functioning of a heat engine is naturally boosted when the difference between the temperatures of the baths is large, in fact in this situation the energy flow and consequently also the work extraction, are favored. Indeed, in this regime, the Carnot efficiency gets close to unity. This basic physical intuition leads us to consider the region of parameters where T 2 T 1 as the most interesting for the [E] operation. More precisely, we take T 2 ω(h) T 1 , withω(h) playing the role of a typical energy scale of the system. Accordingly, the working substance is close to the ground state when it is in equilibrium with the cold bath, and it is close to the maximally mixed state when in equilibrium with the hot bath. As a consequence, the work extracted from the working substance can be written as (see Appendix C for details) which is fully determined by the average level spacingω i , i = 1, 2 for h = h 1,2 . It follows that the optimal work output is reached for for the values of h 1 , h 2 that respectively maximise and minimise the functionω(h) in Eq. 29, namely: where the optimization has to be performed over the values of h compatible with the approximation (32), i.e., such that T 2 ω(h) T 1 . Within the same approximation, the heat engine efficiency reads Remarkably, this choice of h 1 and h 2 allows to maximize both the work output and the cycle efficiency. In Fig. 4a) and b) the exact work output W , and the exact engine efficiency η [E] respectively (dotted lines) are compared with (29), as a function of the chemical potential h. Different colors correspond to different values of α1 = α2 = α. The blue dashed line represents the asymptotic value in the large |h| limitω ≈ 2|h|. The system size is N = 500. The system size is N = 500.
W 0 and η 0 [E] (solid lines) for different values of h 1 , h 2 and with T 1 = 100, T 2 = 0.01 and α 1 = α 2 = 1.5. We notice that the approximation W W 0 breaks down for large values of h 1 , whenω(h 1 ) becomes of the same order of T 1 , while η 0 [E] remains a good estimate of η [E] even in this regime. Finally, let us notice that, regardless of the validity of the approximation, the maxima of η [E] and W are actually close.
In Fig. 5 we plot W and η [E] against h 2 /h 1 , for different values of α, showing that they grow as the range of the interaction increases, signaling a clear advantage of the long-range regime. This advantage can be traced back to the properties of the spectrum of the system, encoded in the average level spacingω. In fact, the minimum ofω(h 2 ), which corresponds to the maximum of both W 0 and η 0 [E] , is affected by the presence of long-range interactions as shown in Fig. 6.

B. Refrigerator operation mode
In the typical situation in which a quantum refrigerator operates we may expect the two temperatures to be pretty similar T 2 T 1 , since we can imagine that also the baths are embedded in the same quantum hardware of the working substance. Moreover, we want the system to be deeply in the quantum regime, accordingly, the temperatures involved in the cycle should be small with respect to the system energy scale, T 2 T 1 ω. Under these assumptions, the heat extracted from the cold reservoir reads Since the above expression is positive definite, we can conclude that within the considered approximation the Otto cycle always operates as a refrigerator. On the other hand, far from the quantum critical points so that Q 2 has an exponentially decaying behavior as T 2 → 0. As h 2 becomes close to h c (α) the spectrum is no longer gapped, so the above considerations do not apply. In this regime, the main contribution to Q 2 comes from the soft modes, resulting in a power law decay in T 2 where z is the dynamical exponent which, in general, depends on α. In particular for the h 2 = 1 ferromagnetic critical point we have while for the h 2 = −1 + 2 α−1 antiferromagnetic critical point, Further details on the derivation of Eqs. (35), (36), (37) are provided in Appendix D. Since close to the ferromagnetic critical point 1/z grows indefinitely as α → 1, the presence of long-range interactions does not result in any advantage at low temperatures. However, close to the antiferromagnetic case, z does not depend on α: in this case the factor K(α) can actually provide an advantage. That this is indeed the case, is confirmed by the data shown in Fig. 7, where Q 2 /N is plotted as a function of h 2 for different values of α and (T 1 = 0.1, T 2 = 0.099). The figure shows a clear advantage as the range of the interaction is increased. Thus, the different low-temperature scalings of Eq. (37) lead to the peaks in Q 2 at h 2 /h 1 = h c (α)/h 1 , corresponding to an enhanced cooling capability at quantum criticality, as shown in Fig. 7. This effect is augmented when long-range interactions are present, leading to larger and larger peaks as α → 1 and showing long-range advantage also in the most significant regime, i.e., the refrigerator operation. It is worth noting that, while the heat engine configuration is optimized by a long-range interacting machine operating close to the ferromagnetic critical point, the refrigerator operates optimally in the vicinity of the antiferromagnetic critical point.

IV. FINITE TIME CYCLE
While in previous sections we have focused on the thermal machines operating in the adiabatic regime, assuming an infinitely slow unitary driving of the working substance,i.e., setting P k = 1 in Eq. (27c), we are now going to consider the effects of a finite time driving and allow the working substance to have dynamically generated defects during the Otto cycle, causing nonadiabatic energy losses. We provide evidence of long-range advantages also in the more realistic case of a finite-time heat engine and refrigerator, by computing the universal scaling of the nonadiabatic energy losses with the driving speed. In contrast with the nearest-neighbor case, such energy losses are universally suppressed when long-range interactions are present.
A. Exact solution of the non adiabatic dynamics and adiabatic perturbation theory In general, the unitary evolution generated by H(h(t)) is such that it only mixes the states |0 k , 0 −k and |1 k , 1 −k , for each value of k. As a consequence, the dynamics of the Kitaev chain can be exactly described by the N independent evolution equations, each restricted to the two-dimensional subspace associated with the corresponding k-mode [89]. These can be cast into a matrix evolution for the Bogolyubov coefficients u k ,v k (9), where H k is given by Eq. (8), with h = h(t). By means of the transformation t = ∆ k (t k − h + δt)/δ, this is mapped onto a Landau-Zener-Stückelberg-Majorana (LZSM) problem [89][90][91][92][93]: where Ω k = δ/∆ 2 k . The exact general solution of Eq. (41) can be written in terms of Weber (or parabolic cylinder) D-functions D ν (z), (see Ref. [89]), leading to with s = (4iΩ k ) −1 , z = √ Ω k t e iπ/4 , and a, b arbitrary complex parameters to be fixed by the initial conditions u k (t i ), v k (t i ). Accordingly, the solution of Eq. (40) reads: where u k (t) = u k (t (t)), v k (t) = v k (t (t)). We can introduce the instantaneous eigenstates of the two-level Hamiltonians H k (t) at time t, given by withū k (h) = cos(θ k (h)/2),v k (h) = sin(θ k (h)/2), where θ k (h) = arctan(∆ k /(h − t k )) is the Bogoliubov angle for a chemical potential h = h(t). The probability for the system to be found in the instantaneous eigenstates during the evolution reads By inserting the expression for u k , v k in Eqs. (42), (43), in the above expression, one obtains an analytical expression for P k (t) [94]. This exact solution, however, is rather cumbersome. Considering the limit of a slow driving protocol δ → 0, with final time τ = |h f −h i |/δ → ∞ allows for a simpler description that captures and better grasp the relevant physics involved in the dynamics. In this regime, the first non-trivial correction to the P k takes the celebrated LZSM form see Ref. [95] for its derivation using adiabatic perturbation theory. Although for finite ∆ k the O(δ 2 ) contributions is leading, as the transition point is crossed, the physics is dominated by the soft modes with small ∆ k . As a consequence, in any relevant thermodynamical quantity, the O(δ 2 ) contribution in the r.h.s. of (48) is negligible with respect to the non-analytic exponential one.
Knowing the probabilities P k is sufficient to compute the energy exchanges with the two reservoirs and with the external drivings, during our finite-time Otto cycle, in Eqs (27). Let us consider a single mode and the corresponding two-level system formed by |0 k , 0 −k , |1 k , 1 −k . It is known [82] that, when P k < 1, a region in parameter space corresponding to the heater [H] appears and it becomes the only possible regime when P k ≤ 1/2, since in this case the energy exchanges become negative definite. This happens for the many-body Kitaev chain as well if the driving is so fast that P k < 1/2 for all the values of k. For any finite-time driving, the presence of finite transition probabilities 1 − P k hinders the engine performance enhancing the irreversible character of the cycle. This can be explicitly proven by computing the so-called entropy production (as shown in Appendix H).
Here, we follow the ideas introduced in Ref. [96], about the dependence of the non-adiabatic energy losses from the velocity δ of a finite time cycle. Indeed, those can be directly related to the universal Kibble-Zureck scaling, i.e. the scaling of the density of excitations generated dynamically as the chemical potential slowly ramps across one of the critical points [97]. Let us consider the heat engine operation in the optimal regime identified in section III A ( T 2 ω T 1 ). In this situation, the first unitary stroke cannot increase the entropy as the system already starts from the maximal entropy state ρ 1 ∝ I. On the other hand, we expect defects to be generated in the fermionic chain during the third stroke of the cycle as the system is almost in its ground state ρ 3 |gs gs| before the unitary evolution takes place. For a slow driving of the system through a quantum critical point, the density of excitations is given by [48] where in the last step we used the result in Eq. 48 and took the N → ∞ limit. Since δ is small, the leading contribution to the integral comes from the soft modes with ∆ k ∼ 0. By considering the asymptotic scaling of the dispersion relation in correspondence of the ferromagnetic critical point ∆ k ∼ |k| min(α2−1,1) (see Appendix B) we find the scaling law [48] The same scaling holds for the nonadiabatic work losses, i.e. the difference between the adiabatic work W ∞ extracted in an infinitely slow cycle, and the work W extracted in the more realistic finite-time case. Indeed, this difference can be expressed as [98] which, in the optimal regime for the heat engine (T 2 ω T 1 ) becomes As ω 1,k remains finite for k → 0, the above expression has the same scaling of Eq. (49), as δ → 0, This result suggests that the actual leading physical mechanism behind the energy losses in a finite time cycle is the defect proliferation induced by the driving. Let us notice that, from the definition of the exponent θ in Eq. (50), as the system is sufficiently long-range (α < 2), then θ LR = 1/(2α − 2) > θ SR = 1/2. This simple observation tells us that, at least in the limit of a slow cycle δ → 0, dynamical excitations are suppressed when long-range interactions are present. This additional source of long-range advantage mitigates one of the main limitations of quantum thermal devices, namely the trade-off between power and efficiency. Moreover, let us stress the fact that the result in Eq. (53) only depends on universal quantities, suggesting that the reported long-range advantage is present in generic long-range interacting systems independently of the microscopic details of the model. Figure 8 shows the nonadiabatic work loss ratio 1 − W/W ∞ as a function of δ, for different values of α. We notice that, excellent agreement is found between the exact numerical data (scatter plots) and the approximated result (52) we used to extract the universal scaling in Eq. (53). Moreover, as predicted by our analytical results, the work losses are widely reduced when α < 2. Finally, a similar reasoning applies to the refrigerator [R] in its most realistic temperature setting T 2 T 1 ω. As discussed in Sec. III B, in this case, the relevant quantity to be optimized is the heat extracted from the cold reservoir Q 2 . Let us then consider the difference between the adiabatic cooling capability Q 2,∞ and the heat extracted in a finite time cycle. In the range of temperatures T 2 T 1 ω this reads In order to determine the scaling of this quantity for a slow cycle (δ → 0) we have to distinguish the case in which h 2 is critical or not. While in the latter case we find the same result of Eq. (53), for h 2 = 1 (i.e. the ferromagnetic critical point) the dynamical scaling is affected by the presence of soft modes in ω 2,k as well, namely ω 2,k ∼ k min(α1,α2,2)−1 (see Appendix B). We find then the two different scaling behaviors Let us notice that, away from criticality, the performances of the thermal machine are enhanced for all α 2 < 2, while for h 2 = 1 we have to require the additional constraint α 1 > 2α 2 − 2 (which is trivially satisfied in the limits α 1 = ∞, α 1 = α 2 < 2). We conclude that we can have long-range advantage for the cooling capability of a finite time cycle as well.

V. CONCLUSION
In this work, we have investigated the performance of a quantum thermal machine consisting of a chain of fermions with power-law decaying interactions, which undergoes a quantum Otto cycle. We exactly computed the energy exchanged during the cycle, which in turn allowed us to provide a fully-fledged characterization of the device. We determined the regions of the parameter space corresponding to the most useful operation modes for quantum technological applications, i.e., the heat-engine [E] and the refrigerator [R] modes. Focusing on these two operation modes we then investigated the role of long-range interactions while optimizing the device performances, detecting several sources of long-range advantage with respect to the corresponding nearest-neighbor case.
Most remarkably, those results in high thermodynamic efficiency even when operating at finite power. Indeed the long-range character of the interactions is known to hinder the proliferation of defects as the critical point is crossed [48]. We show that this mechanism is able to mitigate the energy losses, resulting in an advantage with respect to the short-range counterpart of the device. In particular, we were able to link such losses to the universal scaling of the defects within the Kibble-Zurek picture [97], suggesting this mechanism to be indeed very general. More precisely, as shown in Sec. IV, the nonadiabatic energy losses in the work output (53), for the [E] operation, and in the cooling capability (55), for the [R] operation, are proportional to the density of dynamically generated defects during the evolution. This, in turn, is suppressed by the presence of sufficiently long-range interactions (α < 2). Indeed, the presence of long-range couplings results in a faster power law decay of the density of defects as the speed of the driving goes to zero, see Eq. (50). The universality of the resulting long-range advantage makes systems with long-range interactions promising platforms for the implementation of finite-time many-body quantum thermal cycles with improved performances.
Aside from nonadiabatic losses, long-range couplings in the fermionic chain were found to boost the performance of quantum thermodynamic machines even in the adiabatic regime. In particular, the optimal regime in the [E] operation mode can be obtained by maximizing both the work output and the engine efficiency simultaneously. For a large temperature gradient between the thermal reservoirs, the optimal performance is achieved as the exponent of the power law decaying interaction is decreased. This effect is clearly demonstrated by the plots in Fig. 5, where the work output and the engine efficiency are shown to have a larger optimal value as the exponent α of the power law decaying interaction is decreased. Indeed, long-range interactions generate cusps in the low-energy spectrum of the fermionic quasi-particle that increase the work output of the engine, as follows already from the Einstein approximation, in which only the average level spacing is considered.
At variance, in the refrigerator operation mode [R], a realistic implementation requires the two baths to be close in temperature (T 2 T 1 ) and both deep in the quantum regime (T 1 , T 2 ω). This in turn results in a peak of the cooling capability corresponding to the quantum critical points of the model for any interaction range. The performance of the refrigerator close to the quantum critical points improves even further at small values of α as clearly indicated by Fig. 7, which yields a clear proof of the long-range advantage.
The thermodynamic advantage demonstrated in the present analysis shall become even more prominent in the so-called strong long-range regime (α 1,2 < 1), in which the system loses additivity. However, the thermalization of the system on an accessible time scale is not guaranteed in this regime, as the fermionic dispersion relation becomes gapped [57,58], thus allowing for the presence of long-lived prethermal and quasi-stationary phases [83]. Further investigation is thus needed in order to understand this regime. However, we still may expect to find some source of thermodynamic advantage similar to the one obtained in Ref. [99] where the lack of thermalization is due to the presence of many-body localization.
Finally, our findings could be experimentally checked in nowadays available trapped-ions platforms, which are currently used to realise long-range interacting quantum systems. The presence of a long-range advantage could thus lead to new quantum technological developments, with direct application to the cooling of quantum computers.

Appendix A: Nonadiabatic energy exchanges
In this Appendix we provide the details of the derivation of the explicit expressions for the engine energy exchanges in the generic nonadiabatic case, reported in Eqs. (27) of the main text. Let us start from the definition of the energy exchanges with the two thermal reservoirs involved in the cycle Since E 1 and E 3 are thermal expectations they can be readily expressed as In order to compute E 2 and E 4 , we notice that, as shown in Section IV, the unitary evolution generated by H(h(t)) is such that it only mixes the states |0 k , 0 −k and |1 k , 1 −k , for each value of k. As a consequence, the dynamics of the Kitaev chain can be exactly described by the unitary dynamics of N independent two level systems, i.e., the unitary evolution operator can be decomposed as where the unitary operator associated to each two level system is generated by the time dependent Hamiltonian where |φ ± k are the instantaneous eigenstates of the Hamiltonian defined in Eq. (46). The energies E 2 and E 4 at the end of the two unitary strokes of the cycle can then be written as where E k,i is the energy associated to the kth mode, which can be computed by using the results for the single qubit Otto cycle of Ref. [82], which we briefly summarize here for the sake of completeness. In particular, for the energies E k,i we find where This immediately implies that one of its elements is sufficient to determine all of them, and that the matrix is symmetric: if P k ++ = P k , then P +− = P −+ = 1 − P k , and P −− = P k . Moreover,as shown in Ref. [82], thanks to the fact that H k (t) is invariant under the antiunitary complex conjugation operator we have thatP k = P k . Hence, E k,2 and E k,4 read Then, summing over the k-modes and using the symmetry for k → −k of the Hamiltonian we obtain Finally, inserting the expression for E 1 and E 3 in Eq. (A3) and those for E 2 and E 4 in Eq. (A11) into the definition of the heat exchanges we find which reduce to the expressions in Eqs. (27a) and (27b) once the coefficients f k,i = tan(β i ω i,k /2) are introduced.
Appendix B: Taylor expansion of the spectrum around the critical modes In this Appendix we compute the Taylor expansion of the quasiparticle spectrum (II A) at lowest order in |k − k c |), where k c = 0 at the critical point h = 1, while k c = π at h = −1 + 2 1−α1 . Lets start from the h = 1 critical point, for long-range couplings with 1 < α 1 < 3 and 1 < α 2 < 2, at lowest order in k 0 we find [48] where we have introduced A(α) = sin(απ/2)Γ(1 − α)/ζ(α) and B(α) = cos(απ/2)Γ(1 − α)/ζ(α). While in the short-range case α 1 = α 2 = ∞ we simply have Inserting these expansions in Eq. II A we obtain where α = min(α 1 , α 2 ) and On the other hand in the nearest neighbor case α 1 = α 2 = ∞, we obtain The Taylor expansions around k = π are obtained from the results at k = 0, using the following relation Then applying this property of the polylogarithm to the definition of t k and ∆ k in Eqs. (5), (6), we find The Taylor expansion of t k and ∆ k around k = π follows by applying the expansion around k = 0 to t k and ∆ k with k = 2(k − π) and k = k − π, respectively, then inserting the results in Eqs. (B9), (B10). Then for k π, 1 < α 1 < 3 and 1 < α 2 < 2, we find where D(α 1 ) = (2 3−α1 − 1)ζ(α 1 − 2)/2ζ(α 1 ) and E(α) = (1 − 2 2−α2 )ζ(α 2 − 1)/ζ(α 2 )). While in the short-range case we have Inserting these results in Eq. (II A) we obtain, for α 1 = α 2 = α, and 1 < α < 2 with α = min(α 1 , α 2 ) and Finally, in the short-range case α 1 = α 2 = ∞ we find Appendix C: Work output in the infinite temperature gradient limit In this Appendix we provide the derivation of Eq. (32) for the work output in the limit of infinite hot temperature T 1 ω(h) and infinitely small cold temperature T 2 ω(h), with respect to the typical energy scale of the system which we identify with the average spectrum (29). Let us start from the expression of the energy exchanges of the adiabatic cycle in Eqs. (28a), (28b) , and (28c), then we have to expand the coefficients f k,1 and f k,2 for ω 1,k β 1 = ω 1,k /T 1 → 0 and ω 2,k β 2 = ω 2,k /T 2 → ∞, respectively. In particular, at leading order we find Then, inserting this expressions into Eqs. (28a) and (28b) we obtain Finally, using the first principle of thermodynamics we find the work output with The effect of higher order corrections, leading to the advantage-disadvantage transition outside from the optimal region, is considered in Appendix G.
In this Appendix we provide the detailed derivation of Eq. (37) in the main text for the low-temperature limit of the cooling capability when the device works as a refrigerator. We start from expression (28b) for the heat extracted from the cold reservoir in an adiabatic cycle. Then the leading contribution as T 2 T 1 ω is given by For large system size N 1 we can perform a continuum limit passing from a sum over the k-modes to an integral Since T 2 → 0 (β 2 → ∞), the leading contribution to the integral comes from the modes at which ω 2,k is minimum. Then if h 2 = 1, −1 + 2 1−α1 , the spectrum is gapped and min k [ω 2,k ] > 0. Accordingly, at leading order, we have the exponential behavior On the other hand at the quantum critical points, min k ω 2,k = 0, leading to a power law decay. In particular at the ferromagnetic critical point h 2 = 1, using the dispersion relation in Eq. (B5), we obtain where α = min(α 1 , α 2 ) and the C(α) is given by Eq. (B6). Finally performing the change of variables y = β 2 C(α)k α−1 we find where we have introduced the gamma function Γ(x) = ∞ 0 dyy x−1 e −y . Analogously, at the antiferromagnetic critical point, using the dispersion relation in Eq. (B15), we find Finally, introducing the dynamical critical exponent (38), (39), we obtain the result of Eq. (37), the α dependent prefactor given by Appendix E: Characterization of the long-range advantage in the full parameter space In Sections III A and III B we have seen how long-range interactions lead to a substantial advantage in the optimal regimes of the heat engine [E] and the refrigerator [R] operations, respectively. However, a question may arise on whether such an advantage is preserved in an extended region of the parameters space, or if it is present only for those fine-tuned values. In this Appendix we address this question showing that the long-range advantage actually extends to a connected region of the parameter space before an advantage/disadvantage transition takes place. In particular, we can introduce the advantage parameters which for the [E] operation are defined as where W lr(sr) and η lr(sr) are the engine work output and efficiency in the long-range(lr) and short-range(sr) case, respectively. Then, a long-range advantage is indicated by the conditions ∆W ≥ 0, ∆η ≥ 0, signaling larger work   output and efficiency in the long-range case. Figure 9 shows the ∆W (E1) and ∆η (E2) parameters as functions of the relevant parameters h 2 /h 1 and β 1 /β 2 , in the [E] region (see fig. 3), and for different values of α 1 ,α 2 . In particular, we notice that in the α 1 = α 2 case a connected advantage domain is present (red area in the plots of Fig. 9) near to the optimal region with β 1 /β 2 1. Then, as the temperature ratio is increased, an advantage to disadvantage transition takes place, with the transition point identified by the condition ∆W = 0. On the other hand for short-range pairing α 1 = ∞, only minor advantages are present while the disadvantage region is much more extended with respect to the fully long-range case. This observation justifies our choice to focus mainly on the α 1 = α 2 case in this paper.
Analogously for the [R] operation, we can introduce an advantage parameter as the difference between the cooling capabilities Q 2 , in the long-range and short-range cases This quantity is plotted as a function of β 1 /β 2 and h 2 /h 1 and for different values of α 1 , α 2 , in Fig. 10. Also in this case we notice that an extended advantage region is present for long-range pairing and hopping α 1 = α 2 = 1.5 (see Fig. 10a)), while the advantage almost disappears in when the pairing is short-range α 1 = ∞ (Fig. 10b)).
Appendix F: Hopping and Pairing contributions to the long-range advantage In this, Appendix we extend the analysis of Sections III A to generic values of power law decay exponents α 1 and α 2 . In particular, Fig. 12 shows the work output W and the heat engine efficiency divided by the Carnot efficiency , as a function of h 2 /h 1 for the same parameter values of Fig. 5 in this two opposite cases. We notice that, independently from the values of α 1 and α 2 , both the work output and the engine efficiency have the same qualitative behavior, in particular they both present a maximum for h 2 /h 1 0. Moreover, this optimal value is enhanced as the power law decay exponent of the long-range coupling (α 1 , α 2 or both) is lowered, thus leading to a long-range advantage. Accordingly, the analysis of the three extremal cases with long-range hopping and short-range pairing (Fig.12a-c), long-range pairing and short-range hopping (Fig.12b-d), and equally long-range pairing and hopping amplitudes (Fig. 5), suggests that a similar qualitative behavior may be found for any intermediate values of α 1 and α 2 leading to an advantage whenever at least one of the two couplings is sufficiently long-range. On the other hand, the optimal values of W and η [E] , corresponding to their maxima at h 2 /h 1 ≈ 0, are enhanced as the range of both couplings is increased, i.e., when α 1 , α 2 → 1. This is clearly shown in Fig. 13, where W/N and η [E] /η C [E] are plotted as a function of α 1 and α 2 , for h 2 /h 1 = 0.
Analogously, in the refrigerator operation mode, the heat extracted from the cold bath Q 2 shows the same qualitative behavior independently of the values of α 1 and α 2 . Figure 14 shows Q 2 in the same parameter region as in Section III B of the main text, as a function of h 2 /h 1 in the three extremal cases we considered also for the heat-engine. In particular, Q 2 is maximal when h 2 corresponds to one of the quantum critical points h 2 = 1 and h 2 = −1 + 2 1−α1 , which are represented by the vertical dashed lines in the plots. However, we notice that the presence of long-range pairing amplitude α 2 < ∞ is necessary in order to have a clear long-range advantage for the cooling capability.  Moreover, as shown in Fig. 15, the cooling capability Q 2 at the quantum critical points h 2 = −1 + 2 1−α1 is enhanced when both α 1 and α 2 are small. The analysis carried out in this Appendix justifies our choice to restrict the study presented in the main text to the most interesting case corresponding to equally long-range hopping and pairing amplitudes α 1 = α 2 = α.

Appendix G: Advantage-disadvantage transition
In this Appendix, we provide a detailed analysis of the advantage-disadvantage transition in the advantage parameter (E1), as the temperatures ratio T 2 /T 1 departs from the optimal regime. Figure 11 shows ∆W as a function of β 1 /β 2 , for different values of the system size N . We notice that both the maximum advantage, max [∆W ], and the maximum disadvantage − min[∆W ], increase as N grows. However, a finite advantage region with ∆W > 0 is always present at sufficiently small β 1 /β 2 , even in the thermodynamic limit. More precisely the finite size scaling of ∆W is of the form ∆W ≈ δwN (G1) Consequently the transition point between long-range advantage and disadvantage is identified by the condition δw = 0, which happens at β 1 /β 2 = (β 1 /β 2 ) * , independently from the system size, see the vertical black dashed line in Fig. 11. A good estimate of the transition point, (β 1 /β 2 ) * = (T 2 /T 1 ) * , can be obtained by considering the leading finite temperature corrections to Eq. (32), namely W W 0 + W(T 1 , T 2 ).
We notice that W 1 ∝ β 1 = 1/T 1 , independently of the values of h 1 and h 2 . On the other hand for W 2 , as T 2 → 0 (β 2 → ∞),if h 2 = h c , we have the exponential decay W 2 ≈ e −β2 min k ω 2,k .
On the other hand, if the value of h 2 corresponds to one of the quantum critical points h 2 = h c we find a power law decay where z is the dynamical critical exponent (38), (39). It follows that the finite T 2 corrections are more important when h 2 = −1 + 2 1−α1 , since in this case the dynamical critical exponent is z = 1. This corresponds to the advantagedisadvantage transition point with the smallest cold temperature T 2 = T * 2 (see the white dashed line in Fig. 9). In fact, the same finite temperature behavior is found also for the advantage parameter (E1), since in the T 2 ω T 1 limit, we have where and ∆W 2 has an exponential decay for h 2 = h c , while at quantum criticality it has the power law decay In particular, in order to determine the smallest transition temperature T * 2 , we take h 2 = −1 + 2 1−α1 . Accordingly, we can write the advantage parameter as Finally, imposing the transition condition ∆W = 0 we obtain the transition temperature is present also in the infinitely slow cycle and it is somehow unavoidable since it is due to the fact that the two thermalization strokes of the cycle are intrinsically irreversible. Thus Σ ∞ = 0 only exactly at the Carnot point, where however all the energies exchanges are null Q 1 = Q 2 = W = 0. On the other hand, the second contribution, reading is present only when the unitary strokes are performed at a finite velocity. We notice that each term in the sum of Eq. (H2) is proportional to the nonadiabatic transition probability 1−P k , then showing explicitly that they provide an additional source of irreversibility, resulting in worse efficiency performances in finite time cycles. In fact, comparing (H2) and (52), we notice that in the T 2 ω T 1 limit showing the relation between irreversible entropy production and nonadiabatic energy losses.