Quantum Heat Engines with Complex Working Media, Complete Otto Cycles and Heuristics

Quantum thermal machines make use of non-classical thermodynamic resources, one of which include interactions between elements of the quantum working medium. In this paper, we examine the performance of a quasi-static quantum Otto engine based on two spins of arbitrary magnitudes subject to an external magnetic field and coupled via an isotropic Heisenberg exchange interaction. It has been shown earlier that the said interaction provides an enhancement of cycle efficiency, with an upper bound that is tighter than the Carnot efficiency. However, the necessary conditions governing engine performance and the relevant upper bound for efficiency are unknown for the general case of arbitrary spin magnitudes. By analyzing extreme case scenarios, we formulate heuristics to infer the necessary conditions for an engine with uncoupled as well as coupled spin model. These conditions lead us to a connection between performance of quantum heat engines and the notion of majorization. Furthermore, the study of complete Otto cycles inherent in the average cycle also yields interesting insights into the average performance.


Introduction
Thermodynamics originated as an empirical study of steam engines, which blossomed into a framework of exceptional generality and simplicity. Quantum thermodynamics is an emerging research field that aims to extend classical thermodynamics and statistical physics into the quantum realm, offering new challenges and opportunities in the wake of a host of non-classical features. A dominant interest is to understand energy-conversion processes at length scales and temperatures where quantum effects become imperative. Inspired by our enhanced capabilities towards nanoscale design and control, this endeavour is being pursued by scientists from diverse backgrounds, such as statistical physics, quantum information, quantum optics, many-body physics and so on. In order to lay foundations for technological breakthroughs, a variety of fundamental questions are being addressed, ranging from issues of thermalisation of quantum systems to examining the validity of thermodynamic concepts, such as the definitions of work, heat, efficiency and power at the nanoscale. The accord between quantum mechanics and thermodynamics is yet to fully unfold [1][2][3]. Its fundamental implications have inspired numerous proposals for thermal machines based on quantum working media . Two major issues which are addressed in such proposals are as follows: What are the performance bounds of heat engines working in quantum regime and what are the thermodynamic properties of these quantum systems which control these bounds? The performance analysis of various quantum analogues of classical heat engines serve as a test bed to study different extensions of thermodynamic ideas in the quantum world. With the recent development of quantum information technology [51][52][53][54] and a number of interesting results, the study of quantum heat engines (QHEs) has drawn much interest. In fact, the past few years witnessed conducive studies exploring how quantum statistics, discreteness of energy levels, quantum adiabaticity, quantum coherence, quantum measurement and entanglement affect the operation of heat engines and cycles in various experimental set-ups, including trapped ions, transmon qubits and more .
In this paper, we generalise the above model and treat two arbitrary spins which are coupled through XXX interaction. We derive conditions ensuring the operation of Otto engine with an arbitrary spin (uncoupled model) as well as for the coupled model. We are able to analytically show that coupling can enhance thermal efficiency and derive a new spin-dependent upper bound to the Otto efficiency, which generalizes and obtains, as a special case, the previous bound of Reference [6]. The new bound is dependent on the magnitude of the total spin quantum number s = s 1 + s 2 . Given the arbitrary magnitudes of the spins and complex nature of the energy spectrum, we focus on the worst-case or best-case scenarios (WCS/BCS) to approach the problem. By proceeding in this manner, we observe that WCS implies a certain majorization relation between the canonical probability distributions of the working medium relative to hot and cold reservoirs. Thus, we discover a robust connection between the performance of our thermal machine and the property of majorization.
We also establish the consistency of our model with the second law. As a novel tool, we study complete Otto cycles (COCs) in order to characterise the performance of our engine. For a COC, the working medium starts and ends in the same state. We notice a subtle difference in the sense in which second law may be applied at the COC level compared to the average level. We show that COCs that follow the second law under a certain operation (say as an engine) suggest conditions applicable for the average performance of the machine. In this manner, COCs offer a potentially handy tool of estimating parameters for the operating regime of the machine.
It is apparent that as the quantum working medium becomes complex, an exact analysis becomes intractable. This is especially true when the working medium is neither a few-particles system having a simple energy spectrum nor a medium close to thermodynamic limit where some scaling law may aid in mathematical simplicity [117]. Thus, in order to target this intermediate regime, it is pragmatic to formulate heuristics. A significant motivation for our paper is to explore the use of heuristics in view of the complexity of the given problem. Heuristics have been employed in various disciplines such as cognitive science, behavioural economics and computer science to name a few. A heuristic is a rule of thumb providing insights into the behaviour of a system in the face of complexity or uncertainty [118][119][120]. It must be appreciated that the value of a heuristic lies in providing a shortcut method that requires a simpler analysis, thus trading accuracy and completeness for speed.
The paper is organized as follows. In Section 2, we introduce our model of two coupled spins (s 1 , s 2 ) as the working substance of the Quantum Otto engine. In Section 2.1, various stages of the heat cycle are described, and the positive work condition for the uncoupled model is discussed. The proof for the same is sketched in Appendix A. In Section 3, the spins are coupled, and we find the coupling range in which positive work extraction is ensured (proofs are sketched in Appendices B and C), which is related to the notion of majorization in Section 3.1 and further used to order the system's energy levels for J = 0 in Section 3.2. In Section 4, the conditions for maximum enhancement of coupled system's efficiency over the uncoupled model are discussed. An upper bound to engine's efficiency is also calculated in the considered domain of coupling. A detailed proof for the positive entropy production for the coupled system is sketched in Appendix D. In Section 5, an analysis is carried out by using the notion of complete Otto cycles. Finally, we discuss the results of our analysis in Section 6.

Quantum Otto Cycle
The working substance consists of two spins with arbitrary magnitudes, s 1 and s 2 , coupled by 1-D isotropic Heisenberg exchange interaction in the presence of an externally applied magnetic field of magnitude B along the z-axis. The system Hamiltonian in the first Stage of the cycle can be written as follows: where J > 0 is the strength of the anti-ferromagnetic coupling.
2 } are the spin operators for the first and the second spin, respectively. H int is the interaction Hamiltonian, and H 1 is the free Hamiltonian. We have taken Bohr magneton µ B = 1, and the gyromagnetic ratio for both spins has been taken to be two [121]. Let n = (2s 1 + 1)(2s 2 + 1) be the total number of energy levels with |ψ k as the corresponding energy eigenstates. When the system is in thermodynamic equilibrium with a heat bath at temperature T, the density matrix ρ 1 for the working substance can be written as follows: where P k = e −E k /T /Z are the occupation probabilities of the energy levels, and Z = ∑ k e −E k /T is the partition function for the system. We have rendered the Boltzmann constant k B equal to unity. Let us consider the case where one spin is an integer, and the other is a half integer.
Some examples of such spin combinations are 3 2 , 2 , 1 2 , 2 and 5 2 , 4 . The energy eigenvalues of the Hamiltonian H for a general (s 1 , s 2 ) coupling are shown in Figure 1. It is to be noted that a term 8s 1 s 2 J common in all the eigenvalues has been neglected as the physical properties of the system would be independent of it. The ordering of these energy levels would depend upon the conditions on the parameters which the positive work condition for the system would provide, which will be discussed in the coming sections.

The Heat Cycle
The four stages constituting the Otto cycle are as follows.

Stage 1:
The system is at thermal equilibrium with a heat reservoir at temperature T 1 with energy e k for which its occupation probabilities are p k , and the corresponding density matrix is ρ 1 (here, we are considering two non interacting spins with energy eigenvalues denoted by e k and occupation probabilities by p k ). Stage 2: The system undergoes a quantum adiabatic process after it is isolated from the hot bath, and the magnetic field is changed from B 1 to a smaller value B 2 . Here, the quantum adiabatic theorem is assumed to hold according to which the process should be slow enough so that no transitions are induced as the energy levels change from e k to e k . Stage 3: Here, the system is brought in contact with a cold bath at temperature T 2 (<T 1 ). The energy eigenvalues remain at e k , and the occupation probabilities change from p k to p k with the external magnetic field at B = B 2 , and the density matrix of the system is ρ 2 . Stage 4: The system is detached from the cold bath, and the magnetic field is changed from B 2 to B 1 with occupation probabilities remaining unchanged at p k and energy eigenvalues changing back from e k to e k such that only work is performed on the system during this step. Finally, the system is attached to the hot bath again, and the cycle is completed such that the average heat absorbed is q 1,av = Tr[H 1 ∆ρ], and the net work performed per cycle is w av = Tr[(H 1 − H 2 )∆ρ]. Here, Tr[·] denotes the trace operation, and ∆ρ = ρ 1 − ρ 2 . In this paper, we consider the free Hamiltonian of the form H i ≡ 2B i h 0 (i = 1, 2), where h 0 is an operator. We now have w av = 2(B 1 − B 2 )Tr[h 0 ∆ρ]; therefore, the efficiency in the absence of interaction is as follows.
Let us first discuss the positive work condition when s 1 and s 2 are non-interacting. The energy eigenvalues (e k ) of the free Hamiltonian, written in the order of increasing energy (if one spin is integer and the other is half integer), are listed in Table 1, and as it can be observed many energy levels for the non-interacting system are degenerate. There is only one level with energy proportional to −s as well as s, two levels with energy proportional to −(s − 1) as well as (s − 1) and so on (the proportionality constant always being 2B). Therefore, denoting the degeneracy by "g", we have the following from Table 1: g |s| = 1, g |s−1| = 2, g |s−2| = 3, ..., g |s−r| = 2s 1 + 1 (4) such that the total number of energy levels are as follows. n = (2s 1 + 1)(2s 2 + 1) = 2(g |s| + g |s−1| + g |s−2| + ... + g |s−r| ). (n/2 − 2s 1 ), ..., n/2 −2(s − r)B (n/2 + 1), ..., (n/2 + 2s 1 The Stage 1 occupation probabilities are written as p k = e −e k /T 1 /z 1 , where z 1 = ∑ n k=1 e −e k /T 1 is the partition function of the system, which can be expressed as follows.
The average heat exchanged with the hot reservoir is as follows: where the primed probabilities are tabulated at T = T 2 and B = B 2 . The average heat exchanged with the cold bath is as follows.
Thus, the work performed on average is the following.
The explicit expression of v is given by Equation (A1). Since B 1 > B 2 is assumed, the system works as an engine on average if and only if v > 0. We prove in Appendix A that the condition required to satisfy v > 0 is the following: where θ = T 2 /T 1 . Furthermore, as proved in Appendix A, Equation (9) implies z 2 > z 1 as well as the following. p 1 > p 1 , and p n < p n , From the above conditions, we can make the following inferences. Positive work extraction is favoured when the occupancy of ground (top) level is more (less) at the cold bath than at the hot bath, which suggests that heat is absorbed at the hot bath, decreasing (increasing) the occupancy of the ground (top) level, while heat is released at the cold bath, thus increasing (decreasing) the occupancy of the ground (top) level.
Since the working medium returns to its initial state (restoring the Hamiltonian as well as coming to be in equilibrium with the hot reservoir), the net change in entropy ∆S 0,av is due to the entropy changes only in the heat baths. The decrease in the entropy of the hot bath is −q 1,av /T 1 and increase in entropy of the cold bath is q 2,av /T 2 . Thus, the net entropy change in one cycle is the following.
We have seen that w av > 0 or v > 0 requires Equation (9) to hold. Under these conditions, it follows that ∆S 0,av > 0, and the consistency with the second law is established at the level of average performance as an engine. Similarly, we observe that the efficiency satisfies the following. η 0 < 1 − T 2 /T 1 = η C .

The Coupled Model
Let us now couple the two spins, with J > 0 being the anti-ferromagnetic coupling strength. The corresponding energy eigenvalues are shown in Figure 1b, where the ordering of the eigenvalues can be considered when the coupling parameter J is small. Moreover, as the coupling is switched on, the degeneracy of the previously degenerate levels is now lifted. Let us express an energy eigenvalue of the coupled system as E k = m 1 B − 8m 2 J, where m 1 = −2s, ..., +2s and m 2 can only take positive values including zero, as shown in Table A2 in Appendix C. The values m 1 and m 2 depend on the index k, but we have omitted it here for brevity of notation. Now, the average heat absorbed from the hot bath (Q 1,av ), the heat rejected to the cold bath (Q 2,av ) and the average work performed in one cycle, W av = Q 1,av − Q 2,av , are given as follows: where the following is the case.
The spin dependent factors m 1 and m 2 are obtained from the expressions of the equilibrium occupation probabilities of the energy levels E k (shown in Figure 1), which in general are written as follows.
For explicit expressions of P k , refer to Table A1 in Appendix B. Z 1 is the Stage 1 partition function of the system for which its expression may be rewritten as follows: where Z 1 ≡ 2 ∑ s+1/2 k=1 cosh [2(s − k + 1)B 1 /T 1 ]. Similarly, we can define P k , the canonical probabilities due to cold bath, by replacing B 1 → B 2 and T 1 → T 2 in the above expressions for P k .
For the proof of PWC for the coupled model (Appendix B), we show that for the so-called worst case scenario (WCS) is given by the following: P k < P k , k = 2, 3, ..., n, and P 1 > P 1 , along with Equation (9). It follows that X > 0. Consistent with Equations (15) and (A25), we then calculate the strictest condition on the allowed range of J (Appendix C), which is given by the following.
Therefore, we conclude that X > 0 or PWC is satisfied under Equations (9) and (16), with the latter constituting the sufficient condition for the coupled system to work as an engine.

Majorization
Majorization [122] is a powerful mathematical concept that defines a preorder on the vectors of real numbers. It is particularly useful to compare two probability distributions. We will highlight its occurance in the context of the working regime of our engine by comparing the two equilibrium probability distributions. Now, for the uncoupled model, the relevant probability distributions are the canonical probabilities {p k } and {p k }, which, at finite temperatures, are ordered as p n < p n−1 < · · · < p 1 and p n < p n−1 < · · · < p 1 , respectively. In Lemma A2 of Appendix B, we proved that Equation (9) is a necessary condition that ensures w av > 0 in the regime of the so-called worst case scenario (WCS), given by the following: where the equality holds for B 2 /T 2 = B 1 /T 1 . Therefore, the above relations imply the following.
The above set of conditions (M) is summarised by stating that {p k } majorizes {p k } and denoted as {p k } ≺ {p k }. As a powerful tool, majorization can be used to prove other results. Intuitively, it indicates that the distribution {p k } is more mixed than {p k }. Thus, as an important consequence, is the Shannon entropy of the distribution {p} (proportional to the thermodynamic entropy of the working medium in equilibrium with a reservoir). In fact, this is expected, since the flow of heat for the engine is on the average from hot to cold. Then, along with heat, thermodynamic entropy is also lost to the cold reservoir. However, the condition of majorization is more general than the above mentioned relation between the entropies.
Thus, for the coupled model, we can also write down the set of conditions equivalent to Equation (M), and infer that {P k } ≺ {P k }, which implies S(P k ) ≥ S(P k ). In other words, if the Stage 3 equilibrium distribution majorizes Stage 1 equilibrium distribution, then we have positive work extraction from the coupled system.
It is possible to find a range of parameter values which satisfy Equation (17). In Figure 2, we show the behaviour of (P k − P k ) for (1/2, 1) system. It is observed that (P 2 − P 2 ) changes sign within the range [0, J c ], indicating that every condition of Equation (17) may not hold in this range, especially at high bath temperatures. However, we observe that the majorization conditions continue to hold and {P k } ≺ {P k }, even if P 2 > P 2 (see Figure 3).
Here, J c = 1/3. The value of J for which P 2 − P 2 (red curve) changes sign (from positive to negative) approaches J c for lower temperatures (see also Figure 3).  . Majorization conditions shown by positivity of all quantities P 6 − P 6 (purple), ∑ 6 k=5 P k − P k (green), ∑ 6 k=4 P k − P k (blue), ∑ 6 k=3 P k − P k (brown) and ∑ 6 k=2 P k − P k (red) as function of the coupling strength J for (1/2, 1) system of Figure 2. The parameters are set at B 1 = 4, B 2 = 3, with temperatures (a) T 1 = 4, T 2 = 2 and (b) T 1 = 6, T 2 = 3. The point where the red curve intersects the lower curve is where P 2 = P 2 . It is observed that for higher bath temperatures (for a given ratio T 2 /T 1 ), this point shifts to lower J values.

Energy Level Ordering
The actual arrangement of the energy eigenvalues depends on the positive work conditions derived above. As for the relative position of 2sB energy level, it will not change, because it is the highest energy eigenvalue of the system regardless of the coupling strength J. The ground state or the minimum energy state will be decided as follows.
There are two energy levels −2sB and −2(s − 1)B − 8sJ which can possibly form the ground state of the coupled system, and their energy gap is |2B − 8sJ|. Given that B 1 > B 2 and 0 < J < J c , we can check that the following is the case.
The above implies that 2B − 8sJ > 0, thereby making −2sB the lowest energy of the system and −2(s − 1)B − 8sJ the energy of the first excited state. Now, Equation (18) opens different possibilities for the arrangement of other energy levels. For example, the levels −2(s − 2)B − 8sJ − 8(s − 1)J and −2(s − 1)B have an energy gap of | − 2B + 8sJ + 8(s − 1)J|, and either of them can be at higher energy state than the other, and both the arrangements are acceptable. For the sake of concreteness, we assume the condition that there is no level crossing when B 1 is changed to a lower value B 2 . One method of arranging the energy levels, in accordance with Equation (18), is shown in Figure 1, which is assumed for the discussion that follows. The net entropy production in one cycle ∆S av for the coupled system ∆S av = −Q 1,av /T 1 + Q 2,av /T 2 can be written as follows.
In the above expression, due to Equation (9), the first term is always positive, but since T 1 > T 2 , the sign of the second term depends on Y which may not be positive. We will consider the WCS whereby under Equation (15), all terms in the defining sum Y (Equation (12)) are negative, thus making Y negative definite (note that m 2 > 0 for all k). By defining the following: we have ∆S av = aX − bY 1 . The condition, given by Equation (16), on the coupling strength which ensures W av > 0 implies that a > b. Then, for Y 1 > 0, we have shown in Appendix D that PWC for the coupled system encapsulated in Equations (9) and (16) suffice to prove X > Y 1 and hence ∆S av > 0. This establishes the consistency of our engine with the second law in the considered domain.

Efficiency Enhancement and the Upper Bound
In the above, we have established conditions for work extraction in the quantum Otto cycle for the coupled system and verified consistency with the second law. In this section, we explore how the coupling between the spins may enhance the efficiency of the engine.
The heat absorbed from the hot reservoir is given by Q 1,av = 2B 1 X + 8JY, where X and Y are as defined in Equation (12). From the energy levels diagram, it is clear that the contribution 8JY to the exchanged heat comes solely from levels which depend on parameter J apart from the field B. Now, since Q 2,av = 2B 2 X + 8JY, this 'extra' contribution to heat is not available for conversion into work, and it is wasted if 8JY > 0. However, it may be utilized to enhance the efficiency of the cycle if 8JY < 0, thus effectively decreasing the heat absorbed from the hot reservoir. Remarkably, the WCS considered earlier implies that all terms entering the sum for Y are negative, and so we have Y ≤ 0 with J > 0. Thus, the WCS directly results in a regime where we can expect an enhancement of the efficiency. Thus, for the operational regime discussed in previous sections, we can rewrite the expression for efficiency, η = 1 − Q 2,av /Q 1,av as follows: where (18)), we obtain the following: where the second inequality follows due to the permissible range of J (Equation (16)). Thus, the expression of the following: constitutes an upper bound to the system's efficiency, which is tighter than the Carnot efficiency, and is within the coupling range 0 < J < J c . The above expression bounding the efficiency of the Otto cycle is our main result of the paper. This expression is validated with numerical calculations in the discussion section. Note that η ub given by Equation (22) is dependent solely on the field values and the total spin of the two particles, while it is independent of the bath temperatures. This expression generalizes the upper bound derived earlier in Reference [6] for the ( 1 2 , 1 2 ) system. We close this section with a remark on the three possible spin combinations for our (s 1 , s 2 ) system: • When one spin value is a half-integer and the other is an integer; • When both values are half-integer or both are integers; • When both are of the same magnitude (both as half-integer or integer).
In this paper, we have discussed the first case only. The only difference between the present case and the other two cases is that, for the latter, when the spins are uncoupled, an energy level with zero energy and 2s 1 + 1-fold degeneracy occurs but that does not affect the performance of the system. The reason is that after the coupling is turned on between the spins, this energy state splits into 2s 1 + 1 non-degenerate energy levels, which depend only on the coupling factor J. Since J is kept fixed during the cycle, these levels do not shift in a cycle and hence do not contribute to the average work resulting in the same PWC as already derived for the first case. Similarly, it can be observed that these levels do not change the condition for maximal efficiency enhancement, and same upper bound can be obtained whatever the spin combination may be.

Complete Otto Cycles
The working medium for the classical Otto cycle is usually a macroscopic system amenable to thermodynamic treatment. This medium may be a collection of statistically independent, non-interacting individual quantum systems or elements, such as spin-1/2 particles or harmonic oscillators and so on. In the adiabatic step of the Otto cycle, the thermodynamic entropy of the working medium stays constant. This implies that there is no intrinsic control on the transitions experienced by individual elements of the working medium.
On the other hand, the working medium of a quantum Otto engine consists of individual elements. In a quasi-static cycle, the isochoric steps are stochastic while the adiabatic steps are deterministic. The quantum adiabatic step is executed slowly enough such that no transition is induced between energy levels of the element, which continues to occupy its initial state throughout the process. Thus, at the level of the ensemble, the occupation prob-abilities do not change during this process. Thus, such a process imposes maximal control on the evolution of the isolated element, and it is described by a quantum unitary process.
Still, due to the stochastic nature of the contact with the reservoirs, the element may not return to its initial state after the four steps of the cycle. Usually, we are interested in the average properties of the cycle by which the quantities such as heat and work are defined at the ensemble level. In this section, we focus on the complete Otto cycles (COCs) inherent in the average Otto cycle considered in earlier sections. The reason that Otto cycle is so often studied in the quantum thermodynamics literature is that the contributions towards heat and work can be clearly separated into different steps, which helps in the analysis. This distinction also holds at the level of COCs; the interaction of the working medium with a reservoir involves only the exchange of heat with the reservoir, whereas the quantum adiabatic step involves only work.
Consider the COC shown as an engine in Figure 4. If the working medium starts at energy level e i , then by the end of the four stages it is again found at level e i . Such a cycle can either run forward as an engine or backwards as a refrigerator. Analysing the performance of COCs is much easier since we are dealing with only two levels at a time without invoking occupation probabilities of the levels and any average quantities. Let us represent an energy eigenvalue of the uncoupled system as e ≡ m 1 B, where m 1 varies from m 1 = −2s, ..., +2s. Based on the final ( f ) and initial (i) values of m 1 , let us define the quantity x = m 1, f − m 1,i , ranging as x = ±2, ..., ±4s. Let q 1 , q 2 , w denote the heat exchanged with the hot bath, cold bath and the work performed, respectively.
It is clear that for x > 0 (x < 0), a COC runs as an engine (refrigerator). The net entropy change (∆S 0 ) is contributed only by the reservoirs. Thereby, we obtain the following.
Now, for x > 0, the condition B 2 /T 2 > B 1 /T 1 ensures that ∆S 0 > 0, or we may say that the second law is then satisfied at the level of COC. Note that there is a subtle difference in the statement about the second law at the level of a COC versus the average performance level.
In the former case, x > 0 guarantees the operation of an engine, whereas the additional condition (9), i.e., B 2 /T 2 > B 1 /T 1 makes this operation consistent with the second law.
On the other hand, for the average operation as an engine, we require v > 0, which itself requires the condition (9). The latter then automatically ensures consistency with the second law at the level of average performance.
Moreover, note that we do not impose the second law at the level of a COC, and the net entropy change for a COC may be negative, as for instance with x < 0 or a COC operating as a refrigerator if B 2 /T 2 > B 1 /T 1 . Thus, we do not imply that COCs with ∆S 0 < 0 do not happen. These observations result in the following interesting conclusion about the uncoupled model. A consistency with the second law for the average performance as engine ensures consistency with the second law for a COC as engine and vice versa.
Let us next study the effect of coupling between the spins. Now, there are no degenerate levels. Let us express an energy eigenvalue of the coupled system as E ≡ m 1 B − 8m 2 J, where the m 1 , m 2 values are given in Table A2. The levels with same m 1 were originally degenerate in the uncoupled model. For the coupled model, energy levels belong to the same band if they have the same value of m 1 but have different values of m 2 . Furthermore, note that in every band there is one level that stays at the same energy even after the coupling is switched on. Now, for a COC between any two energy levels of the coupled system, the general forms of heat exchanged with the reservoirs, Q 1 and Q 2 and the work performed, W = Q 1 − Q 2 , can be written as follows: The net entropy change in one cycle is the following.
We discuss the possible COCs as below stated below. 1.
x = 0, y = 0: These cycles occur between any two different energy bands having the same m 2 . Therefore, if such a cycle proceeds as an engine (x > 0), its efficiency is (24), this COC is consistent with the second law for B 2 > B 1 θ.

2.
x = 0, y = 0: These cycles are possible between energy levels of the same band, i.e., having same m 1 . The work performed is zero, and the heat exchanged is Q 1 = 8Jy = Q 2 . Thus, for y > 0, the corresponding efficiency is also zero. 3.
x, y = 0 with the same sign: These cycles are possible between different bands for levels with different m 1 and m 2 . If such cycles proceed as engine, i.e., x > 0 (and y > 0), then the corresponding efficiency is the following.
From Equation (24), this type of COC is consistent with the second law for B 2 > B 1 θ, without imposing any further condition on the coupling strength J ≥ 0. Therefore, if the second law allows COCs with η = η 0 , then it also allows COCs with η < η 0 . 4.
x, y = 0 with opposite signs: These cycles occur between energy levels of different bands with different m 1 and m 2 . If x > 0 for such cycles (and y < 0), the corresponding efficiency is as follows.
From Equation (24), this COC is allowed by the second law if B 2 > B 1 θ and the following is the case.
Now, we look for the values of x and y which place the most stringent condition on the second law (Equation (24)) or, in other words, render ∆S as the least positive. This will be the worst-case scenario (WCS) in this context, as other values of x and y would yield a larger upper bound J a . Thus, the range imposed by the WCS will hold for all COCs, making all of them consistent with the second law.
The first term in Equation (24) takes the minimum value if x = 2. For the second term, let y min < 0 denote the minimum value of y. Then, we obtain −y min = [s + (s − 1) + ... + (s − (2s 1 − 1))] = s 1 (2s 2 + 1). Substituting the above values of x and y in Equation (27), we obtain the following range of J.
Therefore, it follows that for B 2 > B 1 θ and within the range 0 < J < J x , all the COCs perform as an engine and satisfy the second law. Now, from the probabilistic or average analysis, we concluded that the conditions B 2 > B 1 θ and the coupling range 0 < J < J c ensure the average performance as an engine. In order to compare the two ranges for J, note that s 1 (2s 2 + 1) ≥ s, where the equality is obtained for s 1 = 1/2 implying that, in general, J x ≤ J c . This has the following important consequence. The range for the parameter J in which the machine behaves as an engine on average subsumes the range for J in which all COCs, performing as an engine, are also consistent with the second law. Conversely, if we restrict to the range 0 < J < J x , allowing all COCs running as engine to follow the second law, then the average operation as an engine in that range of parameters is also consistent with the second law.
Moreover from Section 5, we learn that out of all the possible COCs with η > η 0 , the maximum possible value of efficiency is obtained from Equation (26) for minimum x, i.e., x = 2 and |y min | = s 1 (2s 2 + 1), given by the following.
This cycle is allowed by the second law for the condition B 2 > B 1 θ and in the 0 < J < J x range of coupling. Interestingly, the coupling range required for W av > 0 goes beyond J = J x since J x ≤ J c . The case of J a = J c is obtained when we substitute x = 2, |y| = s in Equation (27), and out of all the COCs allowed in this range, the maximum efficiency is given as follows: η 0 /(1 − 4sJ/B 1 ). The latter value is same as the upper bound, η ub , inferred by analysing the average performance of the system. As it can be observed, η ub ≤ η max . For the special case of (1/2, s 2 ) working medium, J x and J c values coincide irrespective of the value of s 2 , resulting in η max = η ub .

Discussion
We have analyzed the performance of a quantum Otto engine based on a working medium with a complex energy spectrum. An insight into the possible operational regimes is hard to obtain analytically for such a system. Using a heuristic-based approach and employing techniques such as worst-case/best-case reasoning, we have highlighted a regime in which the machine definitely works as an engine on average. These set of conditions can be related to the concept of majorization for the given model. Thereby, we find that majorization serves as a more robust criterion for positive work extraction from our engine.
We also introduced an analysis based on complete Otto cycles (COCs). Compared to the probabilistic analysis, the COC approach is much simpler and straightforward. The latter utilizes much less information than the 'average' analysis, and the conclusions so obtained may not be as general. However, as a starting point, the criteria for COCs may serve as a useful heuristic to gain insight into the average performance of the Otto machine.
As we have observed, there is an interesting correspondence between the COCs and the average Otto cycle with regard to the validity of the second law. One of our main results is an explicit expression for the upper bound of Otto efficiency for the coupled system. This expression reduces to the one found for the (1/2, 1/2) case, with s = 1 [6], or to the case of coupled, effective two-level systems [18]. The dependence of the average efficiency on coupling factor J and validity of the upper bound is demonstrated in Figure 5. In addition to the above analytic approaches, we may also numerically study the implications of using higher spins on the performance of thermal machines. In order to make a few observations, we note that the higher "s" values shift the maximum of work to the weak coupling regimes as shown in Figure 6a. Thus, higher magnitudes of spin may be a useful resource to achieve more work output for weak coupling strengths. Numerical analysis also shows that increasing the bath temperatures may increase the work output by orders of magnitude (see Figure 6b). We also observe an extended regime of positive work extraction from the system at high temperatures and this effect is more pronounced for lower "s" values. Along these lines, variations of the efficiency and work output with the coupling factor J may be studied, where s 1 and s 2 are varied for a fixed s value. Figure 7 shows different cases for the case of s = 7/2. Note that η ub and J c (which depend on s and not on the values of individual spins) are the same for a given s.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Positive Work Condition (PWC) for the Uncoupled Model
The net work extracted from the system when s 1 and s 2 are uncoupled is given as w av = 2(B 1 − B 2 )v. Since we assume B 1 > B 2 , we need to find conditions for v > 0 to hold. In other words, we have the following: p n/2−2s 1 − p n/2−2s 1 ... + p n/2 − p n/2 + p n/2+1 − p n/2+1 + ... + p n/2+2s 1 +1 − p n/2+2s 1 +1 where r = s − 1/2. Let us denote the term in v with the largest coefficient s as follows.
We will show that L < 0 implies v < 0. In other words, if the term with largest coefficient in w av is negative, the system cannot work as an engine.
Let us look at the explicit expression of L ≡ (p 1 − p n ) − (p 1 − p n ): where z i (B i /T i ) is the partition function for the system given by Equation (5). The above expression is of the following form.
Let us observe the function f (B 1 /T 1 ) = p 1 − p n . First, due to canonical form of probabilities, we know that p 1 > p n , and so f (B 1 /T 1 ) > 0. Then, for a given value B 1 , if we increase the temperature T 1 , thereby decreasing B 1 /T 1 , we know that the difference p 1 − p n decreases and vice versa. This implies that f (B 1 /T 1 ) is a monotonically increasing function of B 1 /T 1 . The same is also true for f (B 2 /T 2 ).
Since f (B/T) > 0 is a monotonic increasing function of B/T > 0, so if L < 0, the following condition must hold.
The above condition further implies z 2 < z 1 , and so we have the following.
We now rewrite the expression of v as follows.

(A7)
Now, if L < 0 and Equation (A4) holds, then in the first term above accompanying the coefficient s, we have the following.
Similarly, in the second term of the expression for v, we have the following: and so on until we have the following in the last term.
It is important to note that Equation (A4) does not imply any definite relation between p k and p k for k = 2, ..., n/2. On the other hand, it is clear from Equation (A7) that p k > p k for all k = 2, .., n/2 would favor the case v > 0. Thus, assuming L < 0, we will now consider the BCS (Best Case Scenario), mathematically written as the following: and show that v < 0. The proof is as follows.
Proof. It has been noted earlier that L < 0 implies Equation (A4) and z 2 < z 1 . From inspection of the form of canonical probabilities, this further results in p k < p k , ∀ k = n/2 + 1, ..., n. Using these relations in the normalization condition of probabilities given as ∑ n k=1 (p k − p k ) = 0, we have s ∑ n/2 k=1 (p k − p k ) < 0 along with the following.
Moreover, under BCS, we have the following.
Adding all the above inequalities, we arrive at the result v < 0, thereby proving Lemma A1.
Using the monotonic property of L, it is obvious that if L > 0, the following must hold.
It can be observed that the above condition implies p 1 > p 1 and p n > p n . Equation (A9) favors v > 0 as it results in the following condition: in all the terms in Equation (A7), as m 1 > 0 for all the upper-half levels, i.e., for k = n/2 + 1, ..., n (see Table A2). Moreover, under Equation (A9), positivity of Equation (A7) is always favored irrespective of the relation between p k and p k for all k = n/2 + 1, ..., n.
As for the rest of the occupation probabilities, Equation (A9) does not imply any relation between them except for p 1 > p 1 . However, it is obvious from Equation (A7) that p k < p k for all k = 2, .., n/2 would not favor v > 0. Thus, assuming L > 0, we will now consider the WCS (Worst Case Scenario), mathematically written as the following: and then show that v > 0. This would prove Lemma A2.
Thus, we can write the following.
These inequalities, along with the normalization of each probability distribution, imply the following.
Now, by using Equation (A11) and the normalization of probability distributions, we can write the following: along with the following conditions.
Moreover, under WCS, we have the following.
Adding all the above inequalities, we obtain v > 0, thereby proving Lemma 2 and concluding that L > 0, or Equation (A9), is a necessary and sufficient condition for positive work extraction from the uncoupled spin system.

Appendix B. PWC for the Coupled Model
When the spins are interacting, the work extracted is given as follows: The explicit expressions of occupation probabilities are given in Table A1. With B 1 > B 2 , we need to find the condition for which we have X > 0, where the following is the case.
P n/2−2s 1 − P n/2−2s 1 ... + P n/2 − P n/2 + P n/2+1 − P n/2+1 + ... + P n/2+2s 1 +1 − P n/2+2s 1 +1 . (A14) As shown in Appendix A, for the uncoupled spins case, the term with the largest coefficient (s) must be positive, i.e., L > 0, for the system to run as an engine and that is possible if the system's parameters satisfy Equation (A9). Now, we are interested in seeking additional conditions which ensure positive work extraction for the coupled case, provided that the uncoupled model works as an engine.
For completeness, we first show that the same conditions as (A9) also serve as PWC for the coupled model. In order to prove it, consider the term L X ≡ (P 1 − P 1 + P n − P n ). We will first show that the opposite condition, given by Equation (A4), yields L X < 0 and so X < 0. Assuming Equation (A4), we have Z 2 < Z 1 as well as the following.
Now, it can be observed from the explicit expressions of P 1 and P 1 that for T 1 > T 2 , Equation (A4) implies P 1 < P 1 . Therefore, we conclude that, under Equation (A4), L X is negative definite. Consider now the expression of X, rewritten as the following.

(A16)
As it can be observed, Equation (A4) or L X < 0 implies that the following conditions: hold in all the terms in Equation (A16) since m 1 > 0. Similar to the uncoupled case, the sign of X does not depend on the relation between P k and P k for all k = n/2 + 1, ..., n.
In this manner, all the relations given by Equation (A18) follow from Equations (A4) and (A17). Moreover, as noted above, Equation (A4) implies P n > P n . Therefore, using all these relations in the normalization condition of probabilities, we obtain the following.
Relations (A18) along with P n > P n imply the following.
Under BCS, we have the following.
Adding all the above inequalities, we obtain the result that X < 0. This means that under Equation (A4), L X and X are negative definite.
On the other hand, Equation (A9) implies Z 2 > Z 1 , which further yields P n > P n . However, unlike the case with the uncoupled model, this does not determine the relative magnitudes of the ground state probabilities P 1 and P 1 (explicit expressions of these probabilities are given in Table A1). Therefore, here, we cannot be sure of the sign of the quantity L X . Now, due to Equation (A9), we note that holds in all the terms in Equation (A16). Moreover, note that X > 0 is always favored under this condition irrespective of the relation between P k and P k for all k = n/2 + 1, ..., n. The relation between P k and P k for k = 2, 3..., n/2 is also not apparent under Equation (A9). As in the uncoupled case, here we will also consider the WCS written as follows. P k < P k ; k = 2, 3..., n/2. (A19) Now, WCS results in the following conditions. P k < P k ; k = n/2 + 1, ..., n − 1.
. Stage 1 occupation probabilities of the energy levels E k of the coupled spin system. s 1 is smaller of the two spins in the terms involving the factor 2s 1 + 1.
Now, as shown in Appendix C, the above inequality yields the strictest condition on the permissible range of J, which is obtained for k = 2, and is given as follows.
It implies that for J to be in the above range, all inequalities (A23) hold good. Now, using Equation (A20) and P n > P n in the normalization condition of probabilities, we have the following.
Adding all the above inequalities, we have X > 0. Therefore, we conclude that, for WCS, the following conditions ensure X > 0: B 2 > B 1 θ and 0 < J < J c , where θ = T 2 /T 1 . Let us sum up the above discussion. There are two relevant cases.
(a) B 2 < B 1 θ, which implies the following: 1 L X ≡ (P 1 − P 1 ) + (P n − P n ) < 0; 2 X < 0, thereby proving that under B 2 < B 1 θ, it is not possible for the coupled system to work as an engine at all.
(b) B 2 > B 1 θ, which implies the following: 1 L X does not bear a definite sign. Although the term (P n − P n ) in L X is positive definite, yet the sign of the other term (P 1 − P 1 ) is not definite; 2 Under WCS, we are able to prove X > 0 for B 2 > B 1 θ, thereby implying that it is a necessary condition for W av > 0. However, WCS also demands P 1 > P 1 or 0 < J < J c . Therefore, the latter constitutes a sufficient condition for W av > 0. 3 When P 1 > P 1 does not hold, L X does not have definite sign. Thus, depending on the control parameters, other terms in X can be positive. In this case, we cannot predict the sign of W av .
Therefore, we conclude the following regarding positive work extraction for the coupled model. a If L X < 0 (which happens for B 2 < B 1 θ), then W av < 0. b If L X > 0 (which happens for B 2 > B 1 θ and 0 < J < J c ), then W av > 0. c If no definite sign can be assigned to L X (which may happen even when B 2 > B 1 θ holds, but with no condition on the range of J), the system may or may not work as an engine.
Appendix C. Condition on J from W av > 0 From the conditions given by Equations (A21) and (A22), we obtained Equation (A23), which results in the following.
The above condition yields different possible ranges for J corresponding to different energies E k . Out of these, the shortest range will clearly be permissible for all energy levels. Thus, we will find the strictest condition on J that ensures W av > 0. Let us express an energy eigenvalue as the following.
As can be observed from Figure 1, there are energy bands in the spectrum of the coupled system such that the energy levels corresponding to the same band have an identical value of m 1 , but different values of m 2 . From the spectrum, we observe that m 1 varies from the minimum value of −2s up to 2s, while m 2 can only take positive values (see Table A2). Table A2. Spin dependent factors m 1 and m 2 when the energy eigenvalues of the coupled system are expressed as the following: E k = m 1 B − 8m 2 J. The energy levels "k" which fall within the same band (i.e., having same m 1 ) have also been specified. . . .
2s 0 n Equation (A25) now takes the following form.
Now within one band (fixed value of m 1 ), it is obvious that the highest m 2 value will give the strictest condition on J. Now, by referring to the spectrum, we infer that for m 1 = −2(s − q) where q = 0, 1, 2, ..., the largest value of m 2 denoted by m 2,L is Substituting these values on R.H.S of Equation (A27), we obtain the upper limit on J as the following. 1 2(2s − q + 1) Now, the strictest condition on the range of J will be obtained for the lowest permissible value of q, i.e., q = 1 (since m 2 = 0 for q = 0). Thus, we obtain the following.
Therefore, we conclude that for the above range, we have W av > 0. Note that Equation (A28) is obtained for m 1 = −2(s − 1) and m 2 = s, which corresponds to the first excited state of the coupled system. Appendix D. Proof for X > Y 1 As discussed in the main text, for proving ∆S av > 0 we need to show the following: for a case where all the terms of Y are negative. We will show that the PWCs, given by Equations (9) and (16), derived for the coupled model are enough to show the above relation and hence ∆S av > 0. As already proved in the previous sections that with B 2 > B 1 θ, the condition J < J c is obtained by combining the following set of conditions and then substituting k = 2.
Equation (A29) also implies maximally negative Y. U ≡ X + Y/s > 0 will now be proved using Equation (A29) where X and Y are given by Equation (12).
Before starting the proof, note that all the levels contribute to X but only the J dependent levels contribute to Y/s. The steps followed for proving U > 0 under relations Equation (A29) are as follows: 1.
We first consider the lower half levels. With m 1 being negative for all k = 1, ..., n/2 (see Table A2), the total contribution from these levels to X takes the following form.
Similarly, the coefficients of these terms in Y/s can be calculated from Table A2 as m 2 /s (note that m 2 > 0 holds for all k). Now, we add these to obtain the coefficients of these terms in U, denoted by m 3 ≡ |m 1 | 2 + m 2 s , which have been listed in Table A3. As can be observed, m 3 has a positive part given by "s" and a negative part, say m 4 . The total contribution from the lower half levels to U is therefore written as follows.  With m 4 < 0, the second part is positive because of Equation (A29), and the first part is considered later on.

2.
We now consider the upper half levels. The total contribution of these levels to X and Y/s is considered separately. The former is given as the following.
With m 1 being positive (Table A2) for all k = n/2 + 1, .., n, the above expression is positive because of Equation (A29). As for these levels' contribution to Y/s, it is given as the following.
Note that not all the levels contribute to Y because many levels do not explicitly depend on J. Here, m 5 and m 6 are the positive and negative parts of m 2 /s, respectively(see Table A4). The second part in the above equation is positive because of (A29), and the first part is considered later on.  (m 5 + m 6 )(P k − P k ) From the first and second points, we now have two parts which are yet to be proved positive. Their sum is given as follows.
n/2 ∑ k=1 s(P k − P k ) + n−2 ∑ k=n/2+1 m 5 (P k − P k ) (A30) Using relations such as P n < P n and P n−1 < P n−1 (from Equation (A29)) and P 1 > P 1 in the normalization condition of the following probabilities: n ∑ k=1 (P k − P k ) = 0 we have the following.
As shown below, m 5 < s. Therefore, with P k < P k , we can safely replace s by m 5 in the above inequality, thereby proving U > 0.
Proof for m 5 < s In order to prove m 5 < s, consider as an example the k = n − 5 level, the explicit expression of the occupation probability is the following: P n−5 = e −2(s−2)B 1 /T 1 +8sJ/T 1 +8J(s−1)/T 1 /Z 1 and m 5 = 2 (see Table A4). On carefully observing the energy spectrum, it can be observed that the energy level corresponding to this occupation probability exists only if the sum of spins "s" in the power of the exponent satisfies 2 < s. Similarly, 1 < s holds in P n−2 and 3 < s holds in P n−9 , ..., P n−6 . Thus, this is true for all the energy eigenvalues. Since P n−5 − P n−5 < 0 (Equation (A29)), we therefore have the following.
Similarly, we have the following.
The last inequality follows from the fact s 1 < s 2 . This proves m 5 < s.
Proof. Proof for X > Y 1 As discussed in the main text, we will be proving U = X + (Y/s) > 0 by using Equation (A29) and the condition P 1 > P 1 . Note that all levels contribute to X but only the J dependent levels (E 2 and E 4 ) contribute to Y. The steps followed for proving U > 0 under relations Equation (A29) are as follows.

1.
We first consider the lower half (k = 1, 2, 3) of the levels. With m 1 being negative (see Table A5), the total contribution from these levels to X takes the form of the following.
Similarly, contribution of lower half levels to Y/s is written as follows.
The total contribution of the lower half levels to U is as follows.

2.
We now consider the upper half levels. The total contribution of these levels to X and Y/s is considered separately. The former is given as follows. 6 ∑ k=4 m 1 2 (P k − P k ) = (s − 1)(P 4 − P 4 + P 5 − P 5 ) + s(P 6 − P 6 ) The above expression is positive because of Equation (A29). As for these levels' contribution to Y/s, it is given as the following.
Y/s = m 2 s (P 4 − P 4 ) = (P 4 − P 4 ) This part is negative and will be considered later on. 3.
From points one and two the following terms in U are yet to be shown positive. 3 ∑ k=1 s(P k − P k ) + (P 4 − P 4 ) Using relations such P 6 < P 6 , P 5 < P 5 (from Equation (A29)) and P 1 > P 1 in the normalization condition of probabilities, given as the following: 6 ∑ k=1 (P k − P k ) = 0 we have the following. 4 ∑ k=1 s(P k − P k ) > 0 =⇒ 3 ∑ k=1 s(P k − P k ) + s(P 4 − P 4 ) > 0 With P 4 < P 4 and 1 < s (s = 3/2 for the present case), we can safely replace s by 1 in the the above expression, thereby proving U > 0.