Universal finite-time thermodynamics of many-body quantum machines from Kibble-Zurek scaling

We demonstrate the existence of universal features in the finite-time thermodynamics of quantum machines by considering a many-body quantum Otto cycle in which the working medium is driven across quantum critical points during the unitary strokes. Specifically, we consider a quantum engine powered by dissipative energizing and relaxing baths. We show that under very generic conditions, the output work is governed by the Kibble-Zurek mechanism, i.e., it exhibits a universal power-law scaling with the driving speed through the critical points. We also optimize the finite-time thermodynamics as a function of the driving speed. The maximum power and the corresponding efficiency take a universal form, and are reached for an optimal speed that is governed by the critical exponents. We exemplify our results by considering a transverse Ising spin chain as the working medium. For this model, we also study how engine parameters like efficiency and power vary as the engine becomes critical.


I. INTRODUCTION
Advances in quantum science and technology have made possible the laboratory implementation of minimal quantum devices such as heat engines and refrigerators using a variety of platforms that include trapped ions [1][2][3], nitrogen vacancy centers [4], and nuclear magnetic resonance experiments [5]. Quantum engines (QE) transform heat and possibly other resources into some kind of useful work [6]. Their study paves the way to identify quantum effects in their performance. In particular, one may wonder whether there exist scenarios exhibiting a quantum advantage with no classical counterpart [4,7,8].
To a large extent, the study of quantum engines has been restricted to single-particle systems [9]. Such devices already display nontrivial features when their operation involves quantum synchronization [10], nonthermal coherent and squeezed reservoirs [11][12][13][14] and quantum measurements [15][16][17], in the presence of quantum coherence over sustained many cycles [18], or in the small action limit, when different cycles become thermodynamically equivalent [19].
Quantum thermal machines with many-body working mediums (WMs) may allow us to harness many-body effects, such as entanglement and other quantum correlations for operation with enhanced power and efficiency [7]. Quantum statistics can boost the performance of Szilard engines [20,21]. Similarly, the performance of quantum Otto Cycles in both the adiabatic [22] and finite-time operation [7] can exhibit an enhancement due to bosonic quantum statistics, while a detrimental one has been predicted in the fermionic case. Other manyparticle effects that can be harnessed for the engineering of QE include super-radiance [23] and many-body localization [24], while novel configurations become feasible, e.g., by using spin networks [25]. Many-particle QE are also required for scalability and the possibility of suppressing quantum friction during their finite-time operation [26][27][28][29][30] which has been explored in the laboratory with trapped Fermi gases [31,32].
Quantum criticality may offer new avenues to boost the performance of heat engines, as a result of the diverging length and time-scales close to a phase transition [33]. The enhancement of microscopic fluctuations to approach Carnot efficiency in finite time was proposed in [34]. Further, the scaling theory of second-order phase transitions has been used to show that the ratio between the output power and the deviation of the efficiency from the Carnot limit can be optimized at criticality [35]. In adiabatic interaction-driven heat engines, quantum criticality has also shown to optimize the output power [36].
In this work, we introduce a quantum Otto cycle with a working medium that exhibits a quantum phase transition. In particular, we consider the family of freefermionic models that include paradigmatic instances of critical spin systems such as the quantum Ising and XY chains, as well as higher dimensional models. As a result, our setting is of direct relevance to current efforts for building many-particle QEs, with e.g., trapped ions. We explore how signatures of universality in the critical dynamics of the working medium carry over the finitetime thermodynamics of the heat engine. Remarkably, we show that the work output of such quantum engines depend on the universal critical exponents via Kibble-Zurek scaling, thus paving the way for the hitherto unexplored field of universal finite-time thermodynamics describing quantum machines driven through quantum critical points. To the best of our knowledge, such a connection has not been explored before, and is the focus of our paper.
In Sec. II, we introduce the model of a many body Otto cycle using a free-Fermionic WM. We discuss Kibble-Zurek scaling and its connection to output work and power of quantum Otto cycles in Section III A, while Section III B introduces an efficiency bound depending on dynamical critical exponent. We focus on the particular example of a transverse Ising spin chain WM in Sec. IV which is further divided into two subsections depending upon the different phases the WM explores during unitary strokes. We also provide analytical expressions for the energies exchanged in each stroke and compare them with numerics. Finally we conclude in Sec. V.

II. MANY BODY OTTO CYCLE
The use of spins as WM opens a wide range of opportunities recognized early on [37,38]. Recent experiments have implemented single-spin quantum heat engine [3,5] and test fluctuation theorems in single strokes [39,40]. WM composed of interacting spins such as multiferroics have been proposed [41,42] whereas it is shown that WMs with cooperative effects boost engine properties [43]. Quantum critical spin systems in quantum thermodynamics have also been considered under adiabatic performance [44,45], shortcuts to adiabaticity [46] and the limit of sudden driving [47,48]. Such settings preclude the study of signatures of universality associated with the quantum critical dynamics in the finitetime protocols, which is our focus.
We consider an Otto cycle with a many-body WM, described by the Hamiltonian with σ + = (σ x + iσ y )/2, σ − = (σ x − iσ y )/2, and σ x , σ y , σ z being the usual Pauli matrices. HereH k is a 2 × 2 matrix in a basis given by Ψ † k = (c † 1k , c † 2k ) where c jk , c † jk (j = 1, 2) are fermionic operators for the k-th momentum mode. Such a Hamiltonian includes widely studied models, such as the transverse-field Ising and XY chains [49][50][51][52][53], and the two dimensional Kitaev model [54][55][56], through suitable choices of λ, a k and b k . This Hamiltonian exhibits a quantum critical point (QCP) at λ = λ c , when the energy gap ∆ = 2 (λ c + a k ) 2 + |b k | 2 between the ground state and first excited state vanishes, for the critical mode k = k c . The density matrix of such a system can be written in a basis consisting of |0, 0 , |1 1k , 0 , |0, 1 2k , |1 1k , 1 2k where the first index corresponds to presence (1) or absence (0) of c 1k fermion. Similarly, the second index corresponds to c 2k fermions. It is to be noted that the unitary dynamics as per the Hamilto-nianH k mixes |1 1k , 0 , and |0, 1 2k only. As we shall see later, the non-unitary dynamics allows mixing along the other two basis too [57,58]. We denote the full 4 × 4 Hamiltonian matrix by H k .
Before dwelling on the dynamics in Fourier space, let us briefly discuss its real space counterpart. One of the prominent instances within the family of Hamiltonians in Eq. (1) is that of transverse Ising model and the XY model in transverse field (in a ring geometry) which takes the real-space form Here i denotes the site index, c i , c † i are Fermionic annihilation and creation operators, respectively, and M i , N i , R i are scalars [50]. Such a Hamiltonian can be generated ,for example, using a WM consisting of interacting-Fermions in an optical lattice setup [59]. If M i , N i and R i are site independent, one can perform Fourier transform of the Hamiltonian which has the same form as given in Eq. (1). We shall discuss more on this Hamiltonian in Sec. IV.
The quantum Otto cycle alternates between unitary and nonunitary strokes. We now describe below the four general stages of the Otto cycle in details (see Fig. 1): FIG. 1. Schematic diagram of a many-body quantum Otto cycle(a) Schematic diagram of a quantum Otto cycle with a many-body system as the working medium. We get a net output work in the heat engine regime (shown by the glowing light bulb). (b) In the equivalent momentum space, the interacting many-body QE can be represented by independent quantum thermal machines corresponding to the different decoupled Fermionic k modes, each acting as a heat engine (shown by glowing light bulbs), or as a refrigerator (shown by the snow-flake), or even as a heat distributor (not shown here).

Stroke 1 (A → B):
The WM is subjected to a constant Hamiltonian (Eq. (1)) with λ = λ 1 , while being coupled to a dissipative energizing bath B E for a time τ E , thus resulting in non-unitary dynamics.
In general the dissipative dynamics undergone by the density matrix ρ(t) is given by with set to unity, and D[ρ] is the non-unitary part of the dynamics generated due to the interaction of the system with the bath. The exact form of D[ρ] depends upon the nature of the bath and its interaction with the system. Here we consider baths with unique steady states. This can be achieved, for example, by coupling the WM to a thermal bath at a finite temperature.
Alternatively, one can consider Markovian baths coupled locally to the Fermionic modes shown in Eq.
(2), with D[ρ] given by Hereκ i is related to system-bath coupling strength for the site i, L i are local Lindblad operators that describe the interaction of the Fermion at site i with the bath. For L i = c i (see Eq. (2)) and site independentκ i , it can be shown that the Fourier transform of D[ρ] does not mix different modes, so that we arrive at mode-dependent non-interacting local baths in the free-Fermionic representation in momentum space as well. The existence of noninteracting Fermionic modes implies the state ρ(t) of the many-body WM can be written as ρ(t) = ⊗ k ρ k (t), with the time-evolution of ρ k (t) given by the differential equation [57,58] The WM is coupled to B E at A of Fig. 1a, and we assume τ E is large enough so that the system reaches the steady state at B. Here κ E j (j = 1, 2, 3, 4) are positive constants related to the energizing bath, which depend on the coupling strength between the WM and the bath. The energy exchanged in this stroke is denoted as Q in .
2. Stroke 2 (B → C): The system is decoupled from the bath at B and λ is varied linearly in time as t/τ 1 from λ 1 (at B) to λ 2 (at C) in a time interval τ 1 , such that the WM undergoes a unitary dynamics described by We consider λ 1 > λ 2 in this paper. Work is done on or by the system in this stroke.

Stroke 3 (C → D):
The WM is now coupled to a relaxing bath B R at C of Fig. 1a, for a time duration τ R , at a constant λ = λ 2 . The evolution equation will be similar to that given in Eq. (5) with appropriate couplings κ R 1 . . . κ R 4 related to B R . Quantum critical dynamics are more pronounced for systems close to their ground states. Consequently, if the aim is to explore the Kibble Zurek mechanism (explained later), we should consider a relaxing bath which takes the WM to its ground state in this stroke. In principle, we can tune the relaxing bath coupling parameters such that it either takes the system to its ground state or to some steady state corresponding to the bath parameters. We denote the energy exchanged in this stroke with Q out .

Stroke 4 (D → A):
The system is decoupled from B R and λ 2 (at D) is varied back to λ 1 (at A) linearly in a time interval τ 2 as t/τ 2 . Once again, work is done on or by the system in this stroke. One can operate the QE in a steady state cycle, by repeating the above described cycle. It is to be noted that depending upon the values of λ 1 and λ 2 , we may or may not cross the critical point. We consider both of these possibilities in this paper.
At the end of any stroke, the energy of the system is calculated using We choose the parameters κ E 1 , . . . , κ E 4 , κ R 1 , . . . , κ R 4 , λ 1 and λ 2 such that energy Q in is absorbed when coupled to B E in stroke 1, while a smaller amount Q out is released when coupled to B R in stroke 3, such that the setup operates as a heat engine, with a net output work W = − (Q in + Q out ). We assume the following sign convention for the energy flows: Q in , Q out , W are positive (negative) if the energy enters (leaves) the WM. For the Otto cycle to operate as a heat engine, we need Q in > 0, Q out < 0, W < 0. On the other hand, Q in < 0, Q out > 0, W > 0 corresponds to a refrigerator, and Q out < 0, W > 0 denotes a heat distributor [60]. We look at the performance of the heat engine in terms of its efficiency η which is given by and in terms of its power (P) output where τ total = τ E +τ R +τ 1 +τ 2 being the total cycle time.
For a WM which can be decoupled into non-interacting momentum modes as shown in Eq. (1), we have where Q in (k), Q out (k) denote the energy flows corresponding to the k-th mode. We note that even if the complete setup acts as a QE, the individual fermionic modes may act as QE, refrigerator, or heat distributor, depending on the details of the operation and WM, see Fig. 1b.

A. Universal Kibble-Zurek scaling in output work
Two of the strokes of the Otto cycle perform unitary dynamics during which a quantum critical point may be crossed depending upon λ 1 and λ 2 . The universal dynamics in terms of excitations produced due to diverging relaxation time at the critical point is a well studied subject [61][62][63], and can be explained through the adiabaticimpulse approximation [64]. Starting from the ground state corresponding to a given set of parameters of the Hamiltonian, if the Hamiltonian has a time dependence such that the critical point is crossed linearly as t/τ , the amount of density of defects (excitations) with respect to the ground state corresponding to the Hamiltonian at final time, follows a universal power-law with the rate of variation 1/τ . The exponent of the power-law is dependent on the equilibrium critical exponents of the quantum critical point crossed. This power-law relation is known as Kibble-Zurek scaling, after its proponents T. W. B. Kibble and W. H. Zurek, and is given by [27,65] where n ex denotes the density of excitations, d is the dimensionality of the system and ν, z are the correlation length and dynamical critical exponents, respectively. The density of excitations n ex in turn gives rise to the excitation energy E ex , i.e., the energy of the system above the instantaneous ground state, which can also be expected to scale with the rate of quench τ [57,65,66].
In order to extend these results to quantities related to quantum Otto engines described before, we need to operate the QE under the following very generic conditions: • The relaxing bath B R takes the WM close to its ground state. This is one of the important conditions in order to arrive at the scaling derived below.
• The WM is driven at a finite rate across a quantum critical point (or points) during the unitary stroke D → A, i.e., τ 2 is finite.
• The energizing bath B E takes the WM to a unique steady state with high entropy.
The last condition of the WM being in a high-entropy state at B can be realized, for example for a thermal energizing bath B E with temperature T h much larger than the energy scales associated with the WM (eg. L −z and |λ − λ c | νz ), such that the corresponding steady state of the WM is close to a maximum entropy state, with all the energy levels equally populated. The unitary stroke B → C cannot change the entropy of the WM, and consequently the states of the WM at C and B remain approximately equal, for any value of τ 1 . We note that, the states of the WM at C and B can also be approximated to be equal if the WM is quenched rapidly across the quantum critical point during the unitary stroke B → C, i.e., τ 1 → 0, for any form of B E or of the state of the WM at B.
The work done is given by where E A , E B and E C are the energies of the WM at A, B and C, respectively, E G A and E G D are the ground state energies of the WM at A and D, respectively, while E ex,A denotes the excitation energy of the WM at A. The implementation of the engine ensures that E G A , E B , E C and E G D are independent of τ 2 , while the Kibble-Zurek mechanism manifests itself through the presence of E ex,A in the output work: Here is the work output in the limit τ 2 → ∞, which depends only on λ 1 , λ 2 , and the steady-state of the bath B E . Remarkably, as seen above (Eq. (13)), the output work shows the same scaling with τ 2 as the excess energy, upto an additive constant. In particular, for systems and quench protocols in which the excess energy is proportional to the density of defects, such as the examples we consider below, one arrives at a universal scaling form of the output work: The above result (Eqs. (13) and (14)) is the highlight of our paper, and to the best of our knowledge, it shows for the first time, the connection between Kibble-Zurek mechanism, which has been traditionally studied in the context of cosmology [67][68][69] and quantum phase transitions in closed quantum systems [52,53,61,62], and the field of quantum thermodynamics. One can easily extend these results to quantum engines involving nonlinear quenches across quantum critical points as well, following the results reported in Ref. [70]. The importance of this connection stems from the identification of universal signatures in the finite-time thermodynamics of critical QE, as well as the optimization of their performance, to which we now turn our discussion.
To this end, we note that in the limit of τ total ≈ τ 2 τ 1 , τ E , τ R , one can use Eq. (14) to arrive at the scaling relation for the output power P where R is a constant independent of τ 2 , such that W − The optimal quench rate τ −1 delivering the maximum power can be found from the condition while the efficiencyη at maximum power iŝ Further, using Eqs. (8) and (14), one can see that the presence of E ex in W as well as in Q in renders the corresponding efficiency η independent of τ 2 , for large τ 2 as seen in the inset of Fig. 2. We note that Q in would involve Kibble-Zurek scaling through Eqs. (11) and (12), even if the states of the WM at B and C are not equal. However, in such a case, the net output work W might involve Kibble-Zurek scaling arising due to the passage from B to C as well, for example if the WM remains close to its ground state at B and τ 1 is finite.
It is worth pointing out that the universal scaling (14) would be modified in the case of sudden quenches or quenches that start or end at the critical point [71], or in presence of disorder [66].
Later, we shall exemplify the novel results Eqs. (13) and (14) with the transverse Ising model as a working medium which has a well studied quantum critical point.

B. Efficiency bound
One can arrive at a maximum efficiency bound η max of the QE, by defining a maximum possible temperature T max and a minimum possible temperature T min . We design the QE such that the maximum (minimum) possible energy gap ∆ max = {∆(λ 1 , k)} max (∆ min = {∆(λ 2 , k)} min ) between two consecutive energy levels is realized at λ 1 (λ 2 ), where the maximum (minimum) is taken over all the k modes and energy gaps. For sufficiently large λ 1 (i.e., (λ 1 − λ c ) νz k z ∀ k), ∆(λ 1 , k) is independent of k, and is a function of λ 1 alone. In analogy with a thermal bath, we define T max through the following relation [72]: Similarly, one can define an analogous minimum possible temperature through the relation: Here we have assumed κ 1 /κ 2 = κ 4 /κ 3 for both the energizing as well as the relaxing bath. The net efficiency of the spin-chain QE is given by where η(k) is the efficiency corresponding to the k-th mode. Therefore defining η max = {η(k)} max we get For dissipative baths acting as thermal baths with mode dependent temperatures, the second law demands that each η(k) should abide by the Carnot bound of maximum efficiency, with the hot and cold temperatures depending on the mode k. Consequently, one can arrive at η max through T max and T min defined above: The minimum possible non-zero energy gap ∆ min between two consecutive energy levels arise at the QCP (i.e., λ 2 = λ c ), when it assumes the value for a WM with length L [33]. Consequently, for a QE operating between a λ 1 and λ 2 = λ c , we get and As can be seen from Eq. (25), η max increases with increasing system size L, thus showing a possible advantage offered by many-body quantum engines over few-body ones.
Interestingly, as discussed above, η max is maximum for an Otto cycle operating between a λ 1 and the QCP λ 2 = λ c . However, we note that η max in general does not provide a tight bound. The equality in Eq. (21) can be expected to hold only in the limit of a WM with mode independent energy gaps. Furthermore, as shown in the example of Ising spin chain in presence of a transverse field WM below, contrary to the behavior of η max , the actual efficiency of the QE, even though bounded by Eq. (22), may peak slightly away from the QCP.
We note that the effective temperatures defined in Eqs. (18) and (19), and consequently also the efficiency bound (25), depend crucially on the condition that the annihilation and creation operators c jk , c † jk cause transitions between adjacent energy levels for the j = 1, 2 fermions, such that the dissipative baths act as thermal baths with mode-dependent temperatures for each mode k.

IV. A TRANSVERSE ISING SPIN CHAIN WORKING MEDIUM
We now exemplify the generic setup proposed above with the exactly solvable transverse Ising spin chain as the WM, and calculate the efficiency and power close to, as well as away from criticality. The Hamiltonian of transverse Ising model (TIM) in spin space can be written as where σ α i denotes the Pauli matrix in the direction α, acting at the site i, and L is the total number of sites or length of the system. Without any loss of generality, we set J to unity. The Hamiltonian (26), when written in terms of Jordan Wigner fermions c i followed by its Fourier transform c k can be rewritten as . Clearly, λ in Eq. 1 corresponds to the transverse field h, a k = 2 cos k and b k = 2 sin k. The QCP where the gap ∆ between the ground state and first excited state vanishes for this Hamiltonian is given by h = ±1 with the critical mode k c = π and 0, respectively. [49][50][51].
Comparing Eq. (1) and Eq. (27), we find that c 1k = c k , and c 2k = c † −k . As before, the four basis corresponds to |0, 0 , |1 k , 0 , |0, 1 −k and |1 k , 1 −k . The full 4 × 4 Hamiltonian matrix H k is given by The above equation (29) reminds us of a multilevel system coupled with a thermal bath, albeit with a mode dependent temperature [72]. However, we stress that in reality the dissipative baths are not thermal, since they are coupled locally to the WM in the momentum space. Also, let us denote µ s related to energizing bath B E with subscript E and that of relaxing bath B R with subscript R.
The QE with TIM as the WM undergoes an Otto cycle, with λ 1 (λ 2 ) replaced by h 1 (h 2 ). In stroke 2, let h(t) is changed linearly from h 1 to h 2 with time t as h(t) = h 1 + (h 2 − h 1 )t/τ 1 with 0 < t < τ 1 , where τ 1 is related to the speed with which h is varied. In the reverse direction during the stroke 4, h is varied as h 2 +(h 1 −h 2 )(t−τ R −τ 1 )/τ 2 for τ R +τ 1 < t < τ R +τ 1 +τ 2 . One can use the state ρ k and the Hamiltonian H k at the end of each stroke and for each k, to calculate the efficiency and the power using Eq. (8) and Eq. (9), respectively. Depending upon the values of h 1 and h 2 , we can have a heat engine exploring different phases of the phase diagram of the working medium. For example, with h 1 , h 2 > 1, we have the engine exploring the paramagnetic phase only without crossing any critical point. When h 1 1 and −1 < h 2 < 1, the engine crosses one critical point and explores the paramagnetic and ferromagnetic phases. On the other hand, for h 1 1 and h 2 −1, the unitary strokes crosses two critical points separating Paramagnetic-Ferromagnetic-Paramagnetic boundaries.
We shall divide the work into (i) Para-Para QE where it crosses two critical points, (ii) Para-Ferro QE with one critical point crossed, and (iii) Proximity of the critical point, and start the discussion with the engine of the first type. Each of these engines have different purposes for study which we describe below in detail. We will also study the operation of the QE as it approaches the quantum critical point.

A. Para-Para QE
A Para-Para QE can be realized with h 1 1 and h 2 −1. As we explain below, the work done in a Para-Para QE admits a closed form expression, which can directly be connected to Kibble-Zurek scaling. In order to explore the Kibble-Zurek scaling in heat engines, it is important that one of the unitary dynamics start from the ground state of the Hamiltonian. We choose the parameters of the relaxing bath such that it takes the system closest to its ground state. The unitary dynamics from D to A will then show the Kibble-Zurek scaling. For this, we fix µ R = 1 and µ R = 0, since it is the µ term which brings the system to the ground state for negative field values. We choose energizing bath B E parameters as µ E < 1 and µ E = 1. Also, we choose h 1 1 and h 2 −1 so that both the critical points h = ±1 are crossed. The ground state for both the field values is paramagnetic where c− particles are also the quasiparticles. This will help in getting closed analytical expressions for Q in , Q out and work done, and finally their dependence on criticality.

Analytical calculations
Our analytical expressions for the various energy values below are obtained for h 1 1, h 2 −1, |h 1 | |h 2 |. We further consider a high entropy steady state at B or τ 1 small, or both, so that one can write the density matrix at C. We first note that the basis |0, 0 , |1 k , 0 , |0, 1 −k , |1 k , 1 −k are also the eigen basis of the Hamiltonian for large |h|; see Eq. (28). To calculate energies E B , E C , E D and E A , at B, C, D and A respec-tively, using Eq. (7), we need to write the density matrix at each of these points. One can see that at B, when the system has reached its steady state after connecting to the energizing bath B E with µ = µ E and µ = µ E , the density matrix takes the form where P B 1 , P B 2 , P B 3 , P B 4 are the populations in the energy levels E 1 , E 2 , E 3 and E 4 of the Hamiltonian with E 1 < E 2 = E 3 < E 4 for h 1 1. Clearly, the order reverses for h −1. The symbol B in superscript represents point B of the cycle. We shall use the symbol D for quantities related to point D for similar reasons. These probabilities can be obtained using the steady state condition of the master equation which gives where µ E = 1 as discussed before. Also, from the normalisation condition we have From (31) and (32), we get the populations in the energy levels when connected to the B E as Using these expressions, we can write the steady state density matrix of the system at B in terms of µ E . It is to be noted that the density matrix is independent of k in these limits. As mentioned before, we choose an energizing bath which results in a high-entropy state at B, or small τ 1 , or both so that ρ C = ρ B . As shown in Appendix A the energy at B (E B ), and C (E C ) can now be written as Since the decay bath takes the system very close to the ground state E D = −N |h 2 | for h 2 −1. We write E A as E G A + E ex,A where E G A is the ground state energy corresponding to the Hamiltonian at A and is equal to −N h 1 . E ex,A is the excess energy, which will show the Kibble-Zurek scaling. The work done W by the system is −(Q in +Q out ), which can be simplified using the above discussion, and can be written as We verify this scaling in Fig. 2, which establishes the novel relation (13) between the work done in a QE and the universal critical exponents of the quantum critical point crossed. For numerical calculations, the initial density matrix is evolved as per the Eq. (29) when connected to bath, whereas in the unitary stroke it is simply The energies at A, B, C, D are calculated to obtain Q in , Q out , and the work done. This work done, upto an additive constant as discussed in the text, shows the universal scaling and is plotted in Fig. 2. We also show the efficiency of the engine as a function of τ 2 in the inset of Fig 2 which approaches a constant value for large τ 2 , also discussed in the text. We now present the other interesting quantity, namely, power as a function of τ 2 in Fig. 3. As discussed in section III A, it shows a peak at τ 2 = τ opt . For the set of parameter values given in Fig. 3, the analytical value of τ opt = 40 whereas the numerically obtained value is 55. We also calculate the analytical as well as numerical values of efficiency at maximum power,η, which are 0.81 and 0.83, respectively. The figure, which is in log-log plot also captures the 1/τ 2 behavior of power for large τ 2 which can be explained using Eq. (15).

B. Para-Ferro QE
We realize Para-Ferro QE by considering h 1 1 and −1 < h 2 < 1, such that only the paramagnetic -ferromagnetic critical point is crossed during the unitary strokes. We consider an energizing bath of the form shown in Eq. (29), and a relaxing bath B R which takes the system close to its ground state. Similar to the previous case of Para-Para engine, the work done W, up to some constant additive factor, will show Kibble-Zurek scaling, as long as the conditions given in Sec. III A are satisfied. This is presented in Fig. 4.

C. Proximity of the critical point
Next we aim at elucidating how the engine parameters change as the critical point is approached. With h 1 1 and h 2 > 0, we choose the energizing bath B E parameters to be µ E = 1 and µ E < 1 and the relaxing bath B R parameters to be µ R = 1 and µ R < µ E . This set of parameters related to the relaxing bath will take the system to some steady state which is not the ground state at D. One can obtain analytical expressions along the same lines as in the para-para section, also presented in Appendix.
Here, we shall focus on behavior of efficiency η and the power P of the engine for different values of h 2 . The final expressions in the limit of large h 1 and h 2 are: In Fig. 5, we present the behavior of power P as a function of h 2 for fixed h 1 , µ E , µ R τ 1 , and different τ 2 values. Clearly, there is a better agreement between numerical and analytical values of P for larger h 2 . Deviation between the two are more pronounced as the critical point is approached, as excitations generated with the crossing of the critical point are not included in the analytical calculations (36). The power P shows a sharp fall for a QE driven across the phase transition. This behavior can also be attributed to the excitations produced in the WM close to criticality, which in turn results in diminishing Q in , and finally, small output power. Also, |0, 0 , |1 k , 0 , |0, 1 −k , |1 k , 1 −k stops being the eigenbasis of the WM for small h 2 . It is to be noted that the different points in Fig. 5 correspond to different engines.
In the inset of Fig. 5, we present the behavior of efficiency as a function of h 2 for different τ 2 values. As in the previous case, there is a good match between the analytical and numerical calculations when away from the QCP, in the paramagnetic phase. On the other hand, analytical calculations of Sec. IV A 1 fail to explain the numerical results obtained close to the QCP and in the ordered ferromagnetic phase. As is expected from Eq. (36), the efficiency is independent of τ 2 when the operation is confined inside the paramagnetic phase, for large h 1 , h 2 . However, for a QE driven across a phase transition (h 1 > 1, h 2 < 1), as shown in Fig. 5 (inset), the results can be expected to depend non-trivially on τ 2 , owing to the dependence of the non-adiabatic excitations on the rate of driving τ −1 2 across the QCP [52,53,[61][62][63].
The variation of |P| and η as a function of τ 2 is shown in Fig. 6. Clearly, slow quench (large τ 2 ) is beneficial for high efficiency at the cost of low magnitude of power, while fast quench (small τ 2 ) results in large output power. This is a manifestation of the so-called tragedy of finitetime thermodynamics leading to a trade-off between efficiency and power.

V. CONCLUSION
We have studied the effect of quantum criticality in quantum thermodynamics, by considering a many-body quantum machine operating close to a phase transition. As a WM for the Otto cycle studied here, we have considered interacting Fermions coupled to local dissipative baths, which in the Fourier-transformed space, can be treated as non-interacting Fermions coupled to local noninteracting Fermionic dissipative baths. This property makes the setup analytically solvable in many regimes. Earlier studies on dynamics of closed many-body systems driven across quantum critical points have shown the existence of universality in different parameters, including defect density [61][62][63] and its fluctuations [73][74][75], fidelity susceptibility [71,76], and other parameters, owing to diverging length and time scales close to a quantum critical point. On a similar note, we have shown for the first time, the existence of such universality in quantum thermodynamics close to phase transitions, in the form of Kibble-Zurek scaling [61,62] in the work output, and the operation of quantum engines close to criticality. Furthermore, we have defined a maximum efficiency bound η max , which scales with the dynamical critical exponent close to quantum criticality, and increases with increasing system size, thus showing the advantage of developing many-body quantum engines.
We have exemplified our generic results using the model of Ising spin chain in presence of a transverse field. Our analytical and numerical results show that, as expected, the work output indeed shows the Kibble-Zurek scaling form, upto an additive constant, for a quantum engine driven across quantum critical points (h 1 1, h 2 −1 or h 1 1, −1 < h 2 < 1). On the other hand, for a quantum engine confined to the paramagnetic phase, the power attains a maxima close to the QCP (h 1 1, h 2 > 1), rapidly decreasing once the WM approaches the QCP (h 1 1, h 2 → 1 + ), diminishing close to zero when the efficiency is maximum. The loss of power can be attributed to the generation of excitations close to quantum criticality.
While we have mainly focussed on Fermionic baths, our results can be expected to be valid for any kind of bath, as long as the conditions stated in Sec. III A are satisfied. In addition, by considering thermal instead of Fermionic bath, our setting can be readily adapted to the characterization of quantum refrigerators.
The class of quantum machines studied here provides an opportunity to scale up quantum technologies to the macroscopic regime, with a complete understanding of their performance. Experimental implementations can be envisioned in an optical lattice setup [59]. Our results should also be of relevance to the scaling of quantum machines using trapped ion chains as a working medium [1][2][3] in which a quantum Ising chain can be emulated [77,78] and in which universal critical dynamics has been studied [79][80][81], with experiments reported to date probing it in the classical regime [82][83][84]. More generally, and beyond specific implementations, our results advance the study of universal critical phenomena in quantum thermodynamics.
Note: After the completion of our work, reference [85] reported the universal Kibble-Zurek scaling in the cumulants of the work distribution associated with the critical dynamics in a single unitary stroke.

ACKNOWLEDGMENTS
It is a pleasure to thank Fernando J. Gómez-Ruiz for useful discussions and comments on the manuscript. UD acknowledges DST, India for INSPIRE Research grant. UD also acknowledge the hospitalities of Donostia International Physics Center, Spain, and IISER Berhampur, India, during her visits. VM acknowledges Amit Dutta for fruitful discussions, SERB, India for Start-up Research Grant SRG/2019/000411 and IISER Berhampur for Seed grant. (A3) For a system of size L, there are L/2 positive k modes so that Since the Hamiltonian is changed suddenly (small τ 1 ) from h 1 to h 2 , the density matrix is not able to evolve resulting to ρ C = ρ B and thus energy at C (E C ) is Since D is in the ground state, we have E D = −L|h 2 |. We write energy at A to be E A = E G A + E ex,A with E G A = −Lh 1 . This gives and The work done, W = −(Q in + Q out ), is thus or equivalently (A10)

Proximity of the critical point
Clearly, there is no change in E B so that it is given by Eq. (A4). For large h 2 with h 2 < h 1 , there will not be any population change in B to C, and hence the density matrix ρ B will be same as ρ C so that E C is also given by Eq. (A5).
The energy at D would be different since the relaxing bath parameters are so chosen that it need not take the system to the ground state. It can be calculated as follows: Here, ρ D k would be similar to ρ B k as given in Eq. (30) with µ E replaced by µ R . Similar calculations give Now, the Q in and Q out for each k mode is Efficiency of the total system can be calculated using η = k Q in (k) + k Q out (k) k Q in (k) (A20) = 2α k [(h 1 + cos k) − (h 2 + cos k)] 2α k (h 1 + cos k) (A21) Power for the total system is defined as P = net work done by the system total cycle time (A24)