Thermodynamics of Permutation-Invariant Quantum Many-Body Systems: A Group-Theoretical Framework

Quantum systems of indistinguishable particles are commonly described using the formalism of second quantisation, which relies on the assumption that any admissible quantum state must be either symmetric or anti-symmetric under particle permutations. Coherence-induced many-body effects such as superradiance, however, can arise even in systems whose constituents are not fundamentally indistinguishable as long as all relevant dynamical observables are permutation-invariant. Such systems are not confined to symmetric or anti-symmetric states and therefore require a different theoretical approach. Focusing on non-interacting systems, here we combine tools from representation theory and thermodynamically consistent master equations to develop such a framework. We characterise the structure and properties of the steady states emerging in permutation-invariant ensembles of arbitrary multi-level systems that are collectively weakly coupled to a thermal environment. As an application of our general theory, we further explore how these states can in principle be used to enhance the performance of quantum thermal machines. Our group-theoretical framework thereby makes it possible to analyse various limiting cases that would not be accessible otherwise. In addition, it allows us to show that the properties of multi-level ensembles differ qualitatively from those of spin ensembles, which have been investigated earlier using the standard Clebsch-Gordan theory. Our results have a large scope for future generalisations and pave the way for systematic investigations of collective effects arising from permutation-invariance in quantum thermodynamics.

Quantum systems of indistinguishable particles are commonly described using the formalism of second quantisation, which relies on the assumption that any admissible quantum state must be either symmetric or anti-symmetric under particle permutations. Coherence-induced many-body effects such as superradiance, however, can arise even in systems whose constituents are not fundamentally indistinguishable as long as all relevant dynamical observables are permutation-invariant. Such systems are not confined to symmetric or anti-symmetric states and therefore require a different theoretical approach. Focusing on non-interacting systems, here we combine tools from representation theory and thermodynamically consistent master equations to develop such a framework. We characterise the structure and properties of the steady states emerging in permutation-invariant ensembles of arbitrary multi-level systems that are collectively weakly coupled to a thermal environment. As an application of our general theory, we further explore how these states can in principle be used to enhance the performance of quantum thermal machines. Our group-theoretical framework thereby makes it possible to analyse various limiting cases that would not be accessible otherwise. In addition, it allows us to show that the properties of multi-level ensembles differ qualitatively from those of spin ensembles, which have been investigated earlier using the standard Clebsch-Gordan theory. Our results have a large scope for future generalisations and pave the way for systematic investigations of collective effects arising from permutation-invariance in quantum thermodynamics.

I. INTRODUCTION
Permutation symmetry plays a fundamental role in many-body quantum mechanics. According to the principle of indistinguishability, quantum states that differ only by interchanges of identical particles cannot be distinguished by any measurement [1,2]. As a result, any dynamical observable of a system of identical particles, including its Hamiltonian, must be invariant under particle permutations. The eigenstates of such observables can be divided into linearly independent sets according to their symmetry type under permutations -for two particles, the eigenstates of a permutation-invariant observable can be chosen as either symmetric or anti-symmetric; for more than two particles, additional, partially symmetric types arise. These symmetry types are preserved under the time evolution of the system. Moreover, they give rise to super-selection rules which forbid transitions between states of different types and imply that coherences between such states cannot be observed.
The symmetrisation postulate is more restrictive than the principle of indistinguishability. It states that, depending on the sort of particle, only quantum states with one specific symmetry type exist: identical Bosons are described by symmetric states, identical Fermions by anti-symmetric ones [2]. This restriction has profound consequences. It implies, for example, that the particles of an ideal quantum gas must be distributed in energy space according to either Bose-Einstein or Fermi-Dirac statistics, which lead to fundamentally different physical properties. On a mathematical level, the symmetrisa-FIG. 1. Permutation-invariant ensembles. A collection of two-level atoms, which can be externally distinguished by their position in a lattice, is coupled to a thermal radiation field whose spatial intensity is represented by the red shaded areas. The system is initially prepared in a fully symmetric state, where all atoms are excited. (a) If the lattice spacing is larger than the typical wave length of the radiation field, the atoms decay independently under the spontaneous emission of photons. At long times, the ensemble settles to a Gibbs state that is a statistical mixture of many-body quantum states with all possible permutation symmetry types. (b) When its typical wave length is comparable to the lattice spacing, the atoms become indistinguishable to the radiation field. As a result, the permutation symmetry of the initial state must be preserved during the decay process, which leads to the collective emission of radiation, i.e., superradiance [3]. The ensemble thereby approaches a non-thermal steady state, which contains only fully symmetric quantum states.
tion postulate is elegantly accommodated by the formalism of second quantisation, which inherently restricts the Hilbert space of a many-body system to its fully symmet-arXiv:2206.12639v2 [quant-ph] 12 Jul 2023 ric or anti-symmetric subspace.
Here we focus on many-body systems that are not subject to the symmetrisation postulate but whose relevant dynamical observables are still permutation-invariant. Specifically, we consider collections of non-interacting multi-level systems which are coupled to an environment that cannot distinguish between their constituents. The particles may be not fundamentally identical as an observer may still be able to tell them apart through an external degree of freedom, see Fig. 1. For clarity, we will refer to this type of many-body system as permutationinvariant ensembles. Such settings can be realised for instance with cold atoms [4,5], ions [6] or artificial atoms [7] and can host remarkable phenomena with superradiance being a prime example [3].
The quantum states of a permutation-invariant ensemble are generally not restricted to a specific symmetry type and span the entire many-body Hilbert space, which renders the formalism of second quantisation inapplicable. Our first major aim is to show that such systems nevertheless admit an efficient theoretical description, which, among other results, provides explicit means to calculate the expectation values of permutation-invariant observables. At the heart of this framework lies Schur-Weyl duality, a powerful tool from representation theory, which makes it possible to endow the Hilbert space of a permutation-invariant ensemble with a universal structure [37,38]. Specifically, this structure consists of a series of invariant subspaces that are associated with both the permutation group and the symmetry group generated by the dynamical variables of the system. For spin systems, these dynamical variables correspond to angular momentum degrees of freedom, which generate the quantum rotation group SU (2). In this special case, which has been studied before [19][20][21][22][23][24][25] and will serve as a reference for our analysis, the Schur-Weyl decomposition reduces to the conventional Clebsch-Gordan series [2].
As a consequence of the super-selection rules implied by the indistinguishability principle, permutationinvariant ensembles do generally not relax to a Gibbs state in a thermal environment. Instead, they settle to a non-trivial steady state that still depends on the initial state in which the ensemble was originally prepared, see Fig. 1. Our second major aim is to systematically characterise these steady states, determine their thermodynamical properties, analyse their internal structure and explore how they can be utilised to enhance the performance of quantum engine cycles. Our mathematical framework does however have a wider scope, which extends beyond steady states. It thus opens an interesting perspective for future research.
A second potential area of application for our mathematical framework is quantum information theory, where permutation-invariant systems have been studied in the context of entanglement and reference frames [39], as well as for describing hidden degrees of freedom in quantum optics [40][41][42]. Similarly, bunching effects in quantum optics have been shown to have implications for energy transfer between light and a mechanical oscillator [43,44]. In general, there is potential for exploring a variety of quantum settings in which the same permutation-invariant structures arise.
The paper is structured as follows. In Section II, we outline the mathematical background necessary for the analysis. Readers familiar with representation theory may skip this section. In Section III, we describe the physical setting and weak-coupling model of thermalisation in open systems, giving the general partially thermalised form of the steady state. We give computable formulas for basic thermodynamical quantities in Section IV and then put these to use in Section V by studying a model of a heat engine, the Otto cycle. Finally, in Section VI, we explore the non-classical implications of higher-order symmetry groups beyond spins, followed by broader perspectives in Section VII.

II. MATHEMATICAL TOOLS
In this section, we introduce the mathematical tools that will be used in our analysis. We provide a brief guide to the Lie group SU(d), its Lie algebra su(d) and their representations. We further discuss Schur-Weyl duality, Young diagrams, group characters and Weyl's character formula. To illustrate these general concepts, we show how they appear in the common theory of angular momentum. Readers already familiar with these mathematical tools may skip to the next section.

A. The special unitary group
The special unitary group of degree d, denoted by SU(d), is the group of unitary d × d-matrices with determinant 1. Any element u of this group can be written in the form u = e ih , where h is traceless Hermitian d×d-matrix. These matrices form a vector space over the real numbers. Hence, upon choosing a basis {x a } in this space, we can express any u ∈ SU(d) as u = exp[i ∑ a θ a x a ] with a = 1, . . . , d 2 − 1 and θ a ∈ R. The matrices x a are called generators of SU(d). They satisfy a set of characteristic commutation relations which define the Lie algebra su(d) with the coefficients f ab c = −f ba c ∈ R being called structure constants (taking values depending on the chosen basis).
In practice, it is often convenient to chose a special non-Hermitian basis for the Lie algebra su(d), which is known as the Cartan basis [45]. To construct this basis, we first define a set of d−1 diagonal Hermitian generators d i , which satisfy These generators form the Cartan subalgebra of su(d).
The remaining d(d − 1) generators, which we denote by e µ , can be chosen such that they obey the commutation relations where the constants v i µ ∈ R form the root vector v µ = (v 1 µ , . . . , v d−1 µ ) associated with the non-diagonal generator e µ . For simplicity we will often refer to µ as a root.
The non-diagonal generators e µ come in Hermitian conjugate pairs. As can be seen from Eq. (3), they obey the symmetry e † µ = e −µ with v −µ = −v µ . In the Cartan basis, the exponential map that connects the group SU(d) with its Lie algebra su(d) reads u = exp i(∑ i α i d i + ∑ µ γ µ e µ ) , where α i ∈ R and γ * µ = γ −µ ∈ C. Using the chain rule for the commutator, one can now derive the remaining commutation relations [e µ , e ν ] = N µν e µ+ν for ν ≠ −µ and (4) which, together with the relations (2) and (3), fully specify the Lie algebra su(d) in the Cartan basis. Here, µ + ν denotes a root with corresponding root vector v µ+ν ∶= v µ + v ν . The coefficients N µν , M µi ∈ R depend on the normalization of the generators, where the N µν obey the symmetry N µν = −N νµ and vanish if there exists no generator with the root vector v µ+ν . The d(d − 1) root vectors v µ are in one-to-one correspondence with the generators e µ and span a (d − 1)dimensional vector space and therefore cannot all be linearly independent. In fact, one can always find a minimal generating set of rootsv µ such that any root vector can be decomposed as v µ = ∑ d−1 ν=1 n µνvν , where the coefficients n µν are integers, which are either all non-positive or non-negative. The vectorsv µ are then referred to as simple roots [45]. A complete set of simple roots has d−1 elements (and there is in general no unique choice).
In physics, the Hermitian generators of SU(d) are related to the dynamical variables of a d-level system. For example, a Hermitian basis of su(2) is given by These matrices represent, up to a factor ̵ h, the dynamical variables of a spin- 1 2 system. Upon setting x 1 = s x , x 2 = s y and x 3 = s z , the generators satisfy the commutation relations (1) with the structure constants f ab c = ε abc , where ε abc denotes the Levi-Civita symbol. The Cartan basis of su(2) is usually constructed by choosing the diagonal generator as d 1 = s z . The non-diagonal generators are then given by the raising and lowering operators e ± = s x ± is y , which obey the commutation relations (3) with the one-dimensional root vectors v ± = ±1; the coefficients appearing in Eq. (5) are given by M ±1 = ±1 2.
Since v − = −v + , either v + or v − can be chosen as a simple root of su(2).
As a second example, we consider su(3), whose two diagonal generators can be identified with the Gell-Mann matrices i.e., d 1 = Λ 3 and d 2 = Λ 8 . In relativistic particle physics, these generators are associated with the dynamical variables Isospin and Hypercharge. A complete set of simple roots for su(3) is given byv +1 = (2, 0),v +2 = (1, √ 3) and the corresponding non-diagonal generators have the matrix form

B. Irreducible representations
A unitary representation of SU(d) on some Hilbert space H associates every u ∈ SU(d) with a unitary operator U (u) such that U (u)U (u ′ ) = U (uu ′ ). If H has no non-trivial subspace that is invariant under the action of all operators U (u), the representation is called irreducible. Any reducible, i.e. not irreducible, representation can be decomposed into a direct sum of irreducible representations, or irreps. That is, there exists an orthonormal basis of H, in which U (u) takes the form with the operators U λ (u) forming irreps of SU(d) on some orthogonal subspaces H λ of H. The label λ thereby represents an ordered partition of some integer n into d non-negative integers λ 1 ≥ λ 2 ≥ ⋯ ≥ λ d . The decomposition (9) may contain multiple terms with the same label λ. We note that the natural representation U (u) = u of SU(d) on H = C d is always irreducible. The concept of representations naturally extends to the generators: any U (u) can be expressed as U (u) = exp[i ∑ a θ a X a ], where θ a ∈ C and the X a are traceless operators satisfying the same commutation relations as the generators x a , i.e.
with the same structure constants f ab c as in Eq. (1). Hence, the operators X a form a representation of the Lie algebra su(d) on H. These operators can be chosen to be Hermitian, in which case the coefficients θ a must be real. In general, a representation can be constructed for any basis of su(d). In particular, we may pick the Cartan basis, whose representation consists of d − 1 Hermitian operators D i and d(d − 1) operators E µ that come in Hermitian conjugate pairs satisfying E † µ = E −µ . Since these operators obey the same commutation relations as the corresponding generators d i and e µ , cf. Eqs. (2)-(5), there exists a basis of H, in which the D i become diagonal and the the E µ assume the roles of raising and lowering operators.
Any irrep of the group SU(d) must derive from an irrep of the Lie algebra su(d). That is, if the operators U (u) do not share a non-trivial invariant subspace, the same must be true for the corresponding representations X a of the generators. The converse statement also holds: any irrep X a of su(d) gives rise to an irrep of SU(d) via the exponential map U (u) = exp[i ∑ a θ a X a ]. Any reducible representation of SU(d) must therefore be associated with a reducible representation of su(d), which can be decomposed into a direct sum of irreps. Specifically, in the basis of H where Eq. (9) holds, we also have with the X λ a forming an irrep of su(d) on the subspace H λ . Since changing the basis of su(d) is a linear operation on both the generators and their representations, Eq. (11) holds for any basis of su(d).
In the previous section, we have seen that the Hermitian generators of SU(d) can be identified with the fundamental dynamical variables of a physical system. This interpretation extends also to higher-dimensional representations of these generators. In particular, a Hermitian irrep of su(2) on C 2s+1 describes a spin-s system. For example, the 3-dimensional analogues of the matrices (6), form an irrep of su(2) on C 3 with X 1 = S 1 x , X 2 = S 1 y and X 3 = S 1 z and represent the dynamical variables of a spin-1 system.

C. Schur-Weyl duality
We now consider the representations of SU(d) on the product space H (n) = (C d ) ⊗n , which is the Hilbert space of an ensemble of n d-level systems. We will refer to the individual systems as particles in the following. Extending the natural representation of SU(d) to H (n) yields the tensor-product representation U (n) (u) = u ⊗n , which is generally reducible. Hence, the operators U (n) (u) can be brought to the block-diagonal form of Eq. (9). A block with dimension d λ = dim H λ can thereby appear multiple times and we denote its multiplicity by m λ . The basis of H (n) that produces the decomposition (9) is called the Schur basis. Finding this basis is generally complicated. However, the dimensions d λ and multiplicities m λ of the individual irreps can be calculated without explicit knowledge of the Schur basis, as we will see in the next section. In this section, we first introduce the mathematical tool that makes these calculations possible.
Schur-Weyl duality asserts that the tensor-product representation U (n) (u) of SU(d) admits the decomposition where λ stands for an ordered partition of n into d integers, U λ (u) is an irrep of SU(d) with dimension d λ and 1 K λ is the identity operator on an m λ -dimensional Hilbert space K λ [37,38]. As will be explained in Sec. II D, each λ indicates a different symmetry type of the wavefunction. More generally, Schur-Weyl duality states that the product space H (n) decomposes as where K λ and H λ carry irreps of the permutation group over n elements S n and the unitary group U (d) respectively [37,38]. In the following, we will therefore refer to K λ as a permutation subspace and to H λ as a unitary subspace. In practical terms, Schur-Weyl duality implies that any operator O on the product space H (n) that is invariant under arbitrary particle permutations takes the form in the Schur basis. The identity (13) then follows by noting that the operators U (n) (u) are permutation-invariant and that every irrep of U (d) remains irreducible when being restricted to SU(d).
We now consider the Lie algebra su(d). From a given basis x a , we can construct the tensor-product representation X Since operators X (n) a , which generate the tensor-product representation of SU(d), are permutation-invariant, they decompose in the same way as the operators U (n) (u). Hence, in the Schur basis, we have where the X λ a form an irrep of su(d) with dimension d λ . Before moving on, it is again instructive to consider the special case of SU (2). The space H (n) is then the joint Hilbert space of n spin-1 2 particles. The operators X , which form a representation of the canonical basis of su(2), correspond to the collective angular momentum operators of the system. Every ordered partition λ = (λ 1 , λ 2 ) of n can now be uniquely associated with an eigenvalue [J(J + 1)] 1 2 of the total angular momentum operator Hence, the decomposition (16) becomes the usual Clebsch-Gordan series, where every irrep of su(2) corresponds to a different value of J.
The Schur basis is given by the collective angular momentum basis { J, m, p⟩}, where m is an eigenvalue of S (n) z and p is an index for the permutation subspace K J . From the theory of angular momenta, we know that m takes the values −J, −J + 1, . . . , J for fixed J [2]. The irrep J thus has dimension d J = 2J + 1. The multiplicities can be determined recursively from the usual rules for the addition of angular momenta. For n = 1, we trivially have J = 1 2 and m 1 2 = 1. For two spins, J takes the values 0 and 1 corresponding to singlet and triplet states; hence, we have m 0 = m 1 = 1. Adding a third spin turns the value 0 into 1 2 , while the value 1 branches into 3 2 and 1 2 , which gives m 3 2 = 1 and m 1 2 = 2. Analogously, we find m 0 = 2, m 1 = 3 and m 2 = 1 for n = 4. This scheme can, in principle, be applied arbitrarily often, although it becomes cumbersome for large n. In the next section, we describe a more efficient method to calculate the multiplicities of individual irreps for SU (2) and SU(d) in general.

D. Young diagrams
An ordered partition λ = (λ 1 , . . . , λ d ) of n can be graphically represented as a set of d left-justified rows made up of λ 1 , . . . , λ d boxes. This representation is called the Young diagram Y λ of λ. For instance, the partition λ = (3, 1) has the Young diagram The boxes of a Young diagram can be considered as containers for indices that label either single-particle states or particles [37]. In the former case, a so-called semi-standard Young tableau is obtained by filling the numbers 1, . . . , d into the Young diagram Y λ such that the sequence of numbers in every row is non-decreasing and the sequence in every column is strictly increasing. The Young diagram (17) thus admits 3 semi-standard Young tableaux for d = 2: Each Young tableau that is constructed in this way represents a many-particle state with a specific permutation symmetry pattern, which depends only on λ. All of these states are linearly independent and states with different symmetry patterns, i.e., states belonging to different Young diagrams, are orthogonal to each other. Since the symmetry patterns are invariant under the action of the tensor-product operators U (n) (u) [37,46], the manyparticle states derived from a given Young diagram Y λ span an invariant subspace of the U (n) (u). Moreover, it can be shown that any invariant subspace of the U (n) (u) can be constructed in this way. It follows that the dimensions of these subspaces equal the dimensions d λ of the corresponding irreps U λ . Thus, d λ can be determined by counting the admissible Young tableaux for a given Young diagram. This argument leads to the general formula [38] For d = 2, we have d λ = λ 1 − λ 2 + 1 and, upon identifying J = (λ 1 − λ 2 ) 2, we recover the result d J = 2J + 1.
To determine the multiplicities of the irreps U λ n , we have to construct a different set of standard Young tableaux by filling the Young diagrams with the numbers 1, . . . , n, which label the particles of the ensemble, such that the sequence of numbers in every row and column is strictly increasing. For the Young diagram (17), this rule yields 3 possible standard Young tableaux given by The number of possible Young tableaux that can be obtained in this way from a given Young diagram Y λ gives the multiplicity of the corresponding irrep U λ n and can be determined by combinatorial arguments, which yield the closed-form expression [38] For the special case d = 2, we thus have, upon substituting λ 1 = n 2 + J and λ 2 = n 2 − J, It is easy to check that this formula reproduces the multiplicities that we found iteratively for n ≤ 4 in the previous section.

E. Characters
The final tool that will become essential to our analysis is the character associated to an irrep of SU(d) [47]. In general, the character χ U (u) is a function on the group that associates to each element u the trace of the corresponding representation, that is, Here, we are specifically interested in the characters of the Cartan subgroup of SU(d), whose elements u c can be expressed solely in terms of the diagonal generators d i . That is, we want to calculate objects of the form for an irrep U λ (u) of SU(d) corresponding to the ordered partition λ of n. A general expression for these characters is provided by Weyl's character formula [48], which does not require explicit knowledge of the representation and can be applied as follows. We first define the d×d-matices c i , whose ith diagonal element is 1 while all other elements are 0 [49]. The diagonal generators d i can now be expanded in these matrices as d i = ∑ j q ij c j with ∑ j q ij = 0, since the d i are traceless. Once the coefficients q ij have been determined, we can formally express the representations of the diagonal generators as D λ i = ∑ j q ij C λ j such that Eq. (24) becomes This quantity can now be calculated by means of the formula [50] where δ was introduced in Eq. (19) and A denotes a d×dmatrix with elements To show how this recipe works, we consider again the special case of SU (2). With d 1 = s z the coefficients q ij become q 11 = 1 2 and q 12 = −1 2. Hence, we haveα 1 = α 1 2 andα 2 = −α 1 2 and Eq. (26) gives where we have substituted J = (λ 1 − λ 2 ) 2 in the second line. We note that this result could have been obtained directly from Eq. (24) upon recalling that D 1 = S J z represents the z-component of a spin-J system and thus has the matrix form For d > 2, however, such a direct approach is no longer feasible. The character formula (26) then becomes a valuable tool, as we shall see in the following.

III. SET-UP
Permutation-invariant systems encompass all those in which the dynamics are generated by Hamiltonians that are unchanged under exchange of any subsystems. Here, our main objects of interest are ensembles of n identical, non-interacting d-level systems, to which we refer as particles. In the following, we show how the dynamics and stationary states that emerge when such an ensemble is coupled to a thermal bath can be systematically described using the tools of the previous section.
A. Multi-frequency systems

Dynamics
We assume that the Hamiltonian of a single particle has the form where Ω sets the overall energy scale, the a i ∈ R are dimensionless constants and the d i are the diagonal generators of SU(d). The Hamiltonian of the ensemble is thus given by with the D i being the diagonal generators of the tensorproduct representation of SU(d) on the Hilbert space H = (C d ) ⊗n . Note that, from here onward, we drop the index (n) that was used in Sec. II C to denote the tensorproduct representation unless it is required for clarity. The ensemble is coupled to a thermal bath at inverse temperature β = 1 k B T , where k B denotes Boltzmann's constant. We assume that the bath cannot distinguish between the particles of the ensemble so that the systembath coupling can be described in terms of the collective interaction Hamiltonian Here, the E µ are the non-diagonal generators of the tensor-product representation of SU(d), the Hermitian operators B µ = B −µ correspond to observables of the bath and ∆ sets the coupling strength. From the commuta- we can now determine the relevant Bohr frequencies ω µ of the system (i.e., the gaps in its spectrum), which are defined by the relation where v i µ is the element i of the root vector associated with the generator e µ .
To describe the dynamics of the ensemble, we apply the standard weak-coupling, Born-Markov and secular approximations. These approximations require that the time scale of the system-bath interaction, which determined be the inverse coupling strength, is much larger than the relaxation time of the bath and the time scale of the bare system, which is determined by its Bohr frequencies [51]. Under these conditions, one can derive the collective weak-coupling master equation where ρ t denotes the state of the ensemble, the E ν and E † µ play the role of jump operators and the Lamb shift is given by where s µν ω is the anti-Hermitian part of the bath correlation matrix [51]. H LS commutes with the system Hamiltonian H and does not enter into the steady state. The sum in Eq. (33) runs over all pairs of indices ⟨µ, ν⟩ for which ω µ = ω ν = ω. The complex coefficients Γ µν ω are determined by the bathcorrelation functions. They form a positive-semidefinite Hermitian matrix and obey the detailed balance condition Γ µν −ω = e β ̵ hω Γ νµ ω [51]. Here, we further assume that the bath couples to all collective modes of the system independently so that the matrix Γ µν ω has full rank for every ω ≠ 0.
Schur-Weyl duality now makes it possible to simplify the dynamics of the system as follows. We first recall that any operator O that is invariant under arbitrary particle permutations takes the block-diagonal form in the Schur basis, where the index λ corresponds to an ordered partition of n into d integers, 1 K λ is the identity operator on the permutation subspace K λ and the operator O λ acts on the unitary subspace H λ . This result implies that any permutation-invariant state of the ensemble can be written as Here, m λ = dim K λ is the multiplicity of the irrep λ and we have applied the normalisation condition Tr ρ λ t = 1 so that the p λ ≥ 0 add up to 1 and can thus be regarded as the probabilities for the system to occupy the diagonal block λ. Since the operators H, H LS and E µ , which enter the collective master equation (33), all decompose according to Eq. (34), the block structure of the state (35) is preserved at any later time. Thus, if the system is initially in a permutation-invariant state, the blockoccupation probabilities p λ t = p λ are conserved and the sub-states ρ λ t follow the master equation (33) with H, H LS and E µ replaced by H λ , H λ LS and E λ µ , respectively. It is now convenient to introduce the reduced statẽ which is obtained by tracing out the permutation subspaces K λ . This state lives on the reduced Hilbert spaceH = ⊕ λ H λ , which contains all degrees of freedom that are available to an observer who has access only to permutation-invariant observables; the expectation values of such observables are given by Tr As we will see in Sec. IV, the reduced state fully determines the operationally accessible thermodynamical properties of the system. Also note that for any general state ρ t , the reduced stateρ t must always be block-diagonal. That is, there is a superselection rule preventing the existence of superpositions of different λ on this space.

Steady states
In the long-time limit t → ∞, the ensemble settles to a steady state ρ ∞ . For n > 1, this state is not unique. It does, however, admit a universal structure. In App. A we prove a general theorem, which shows that, as long as the dissipative part of the master equation contains a set of simple roots, the full-rank condition on the matrix Γ µν ω and the detailed-balance condition Γ µν −ω = e β ̵ hω Γ µν ω imply that all steady states take the block diagonal form in the Schur basis. Here, the operators σ λ , which have trace 1, and the block-occupation probabilities p λ depend on the initial state [52] and is a Gibbs state with respect to the irrep H λ of H, where the partial partition function Z λ β is fixed by the condition Tr γ λ β = 1. That is, every block thermalises independently, subject to λ being an effective conserved quantity. We stress that this result holds for any initial state, regardless of whether it has the block structure (35). Moreover, our theorem holds more generally for any permutation-invariant system Hamiltonian, which may in principle also contain interactions between particles, and even if the detailed-balance condition does not hold, in which case the Gibbs state γ λ β has to be replaced with some general unique state ρ λ -for details see App. A.
A natural choice for the initial state of the ensemble, on which we will focus in the following, is given by the Gibbs state Here, β 0 = 1 k B T 0 is the inverse temperature of some environment, in which the ensemble has initially thermalised, and Z β0 = z n β0 , where z β0 = Tr e −β0h is the single-particle partition function. The block-occupation probabilities are now given by p λ = p λ β0 = m λ Z λ β0 Z β0 . The steady state thus becomes Finally, tracing out the permutation subspaces in Eqs. (39) and (40), respectively, gives the reduced Gibbs and steady states We note that, for β = β 0 , the Gibbs state γ β0 is a stationary solution of the master equation (33), which implies that the ensemble remains in thermal equilibrium, i.e., ρ β,β =γ β . However, for β ≠ β 0 the properties of the steady state can deviate substantially from those of a thermal state, as we will show in Sec. IV.

B. Spin systems
As a reference for our results, we will consider ensembles of non-interacting spin-s particles, which have been analysed in detail in Ref. [20]. The dynamical variables of a single spin-s system, S s x,y,z , form an irrep of the Lie algebra su(2) on the Hilbert space C 2s+1 . The corresponding tensor-product representation on the ensemble Hilbert space H = (C 2s+1 ) ⊗n is given by Hence, with the single-particle Hamiltonian h = ̵ hΩS s z , the ensemble Hamiltonian is H = ̵ hΩD 1 with D 1 = S s(n) z .
Choosing the system and bath to couple instead via a collective spin, the Lamb shift and the jump operators in the master equation (33) are then given by H LS = ̵ h∆(s ++ . Note that, in contrast to general d-level systems, spin systems are characterised by a single Bohr frequency ω ± = ±Ω. For s > 1 2, the irrep S s x,y,z is no longer the natural representation of su (2). Therefore, the Schur-Weyl duality does not apply in this case. Nevertheless, the operators D 1 and E ± still take the block-diagonal form (34) in the collective angular momentum basis J, m⟩, where J is the total angular momentum quantum number and m is the magnetic quantum number. However, the blocks corresponding to each J are subdivisions of those implied by Schur-Weyl duality. Consequently, the results of Sec. III A equally apply to spin-s particles, where the irrep index λ has to be replaced with J. The crucial difference is that we can no longer use the formulas (19) and (22) to calculate the dimensions and multiplicities of the individual irreps. Instead, the dimension of an irrep J is now given by d J = 2J + 1 and the multiplicities have to be determined from the standard rules for the addition of angular momenta, which lead to the recursion relation where n is the number of particles. The theorem in App. A also holds for spin systems, giving rise to steady states of the same form as Eq. (37) but with λ replaced by J. This also settles a conjecture from Ref. [20] that the off-diagonal blocks vanish in the steady state.

IV. THERMODYNAMICAL QUANTITIES
In this section, we show how the steady-state energy, reduced entropy and reduced non-equilibrium free energy can be calculated for the set-up laid out in the last section. To illustrate our results, we compare an ensemble of spin-1 particles with an ensemble of 3-level particles with SU(3) symmetry. We further derive explicit expressions for our thermodynamical quantities of interest in certain limiting cases as functions of s and d for spin-s particles and general d-level particles, respectively.

A. General expressions
For a general state ρ, and some fixed inverse temperature β, the steady-state energy, reduced entropy and reduced non-equilibrium free energy are defined according to the block structure (36) as where S(ρ) = − Tr[ρ ln [ρ]] denotes the von Neumann entropy. The reduced quantitiesS(ρ) andF (ρ) have been constructed so that they do not contain entropic contributions from the degeneracy spaces K λ . This approach is motivated by the assumption that all degrees of freedom that interact with the experimenter's apparatus or the environment are represented by permutationinvariant observables, which do not give access to any information stored in the subspaces K λ . As a result, the quantitiesS(ρ) andF (ρ) can, like E(ρ), be expressed as averages over their respective counterparts on the individual irreps λ with respect to the probability distribution p λ . Note that, for steady states of the form ⊕ λ p λ σ λ ⊗ ρ λ , we can write the reduced entropy as a dif-ferenceS(ρ) = S(ρ)−S(ρ), whereρ = ⊕ λ p λ σ λ is obtained from ρ by tracing out the unitary subspaces.
For the steady state (40), the thermodynamical quantities (45)-(47) become with p λ β0 = m λ Z λ β0 Z β0 . To evaluate these expressions, we have to calculate the partial partition functions Z λ β , the total partition function Z β and the multiplicities of the individual irreps m λ .
To this end, we first observe that Z λ β can be written as where we have replaced a i by −iα i β 0 in the second expression. Since the D λ i are diagonal generators of an irrep of SU(d), this quantity can be calculated by means of analytic continuation of Weyl's character formula, which we have discussed in Sec. II E. For d = 3, we find [50] where x j = λ j − λ j+1 and λ denotes the ordered partition (λ 1 , λ 2 , λ 3 ). Analogous formulas can in principle be derived for any larger d. The single-particle partition function z β is obtained as a special case for λ = (1, 0, . . . , 0). The total partition function is then given by Z β = (z β ) n .
For an ensembles of spin-s systems the partial partition function over the irrep J is given by as can be easily verified by direct computation. If the individual particles are described in terms of the natural representation of SU(d), the multiplicities m λ can be obtained from the combinatorial formula (21), which yields for d = 3 withλ j being defined in Eq. (19). For spin systems with s > 1 2, the multiplicities can be obtained from the recursion relation (44). With these prerequisites, the quantities (48)-(50) are computationally accessible.
B. SU(2) vs SU (3) As a first application of our theory, we investigate how the transition from SU(2) to SU(3) changes the thermodynamical properties of the steady state. For convenience, we now set Boltzmann's constant to 1 and rescale all energies and temperatures with ̵ hΩ so that β, β 0 and h, H, H λ become dimensionless from here onward.
As the arguably simplest setting in which SU(2) and SU (3) can be compared, we consider ensembles of 3-level systems with excited, intermediate and ground states +⟩ , 0⟩ and −⟩. To isolate the effect of the symmetry group, we pick the single-particle Hamiltonian In an SU(2) description, this Hamiltonian corresponds to a spin-1 system with h = d 1 = S 1 z , where S 1 z was defined in Eq. (12). In terms of SU(3), the Hamiltonian (55) describes a ladder system with h = (d 1 + where we have identified the diagonal generators d 1 , d 2 with the Gell-Mann matrices from Eq. (7). Although the single-particle Hamiltonian is identical in both cases, the steady state (41) is not the same for SU(2) and SU(3), since different system-bath couplings (Eqs. (43),(31)) lead to different dissipation mechanisms. For SU(2), the ensemble relaxes via a single dissipation channel described by the Lindblad operators E ± , which represent the non-diagonal generators By contrast, for SU(3), relaxation to the steady state occurs via 3 dissipation channels, whose Lindblad operators E 1± , E 2± and E 3± represent the generators  (39) at β0. For comparison, the dot-dashed lines labelled "distinguishable" show the fully thermal state energy E β,β (top), and the changes in entropy (middle) and free energy (bottom) during a full equilibration process of distinguishable particles from β0 to β. For all plots, we have chosen the single-particle Hamiltonian (55) and set n = 10 and β0 = 2. Note that, for β = β0, the initial Gibbs state γ β 0 is a stationary state of the master equation (33). Therefore, all curves intersect at β = 2, and all entropy and free energy differences are 0 at β = 2 in the lower two plots.
Hence, in the SU(2)-case the Lindblad operators induce coherent superpositions of jumps between intermediate and excited and ground and intermediate states, while for SU(3) all possible transitions in the single-particle system are addressed separately by the bath. This difference changes the structure of the steady state and thus its thermodynamical properties. In Fig. 2, we plot the internal energy E β,β0 , and changes in reduced entropy and reduced free energy during the relaxation process, which are given by From these plots, we immediately note two features. First, as expected, all three quantities deviate from their thermal equilibrium references in the steady state (40). Second, although we use the same single-particle Hamiltonian for both symmetry groups, these deviations are typically stronger for SU(2) than for SU (3). This result suggest that the 3-channel dissipation mechanism, which applies to systems with SU(3)-symmetry, allows the system to come closer to thermal equilibrium than the single-channel mechanism applying to spin systems.

C. Limiting cases
We now return to the general case of ensembles with SU(d)-symmetry. In the following, we derive explicit expressions for the thermodynamical quantities (48)-(50) as functions of d in special limits to show how our theory can be applied in practice. For comparison, we also analyse spin-s ensembles.

Low initial and high bath temperature
In the limit β 0 ≫ 1, the ensemble will almost exclusively occupy the trivial symmetric subspace, which is spanned by permutation-invariant many-particle states. To understand this effect, we consider the probability p g β0 = ⟨g n γ β0 g n ⟩ of finding the ensemble in its nondegenerate ground state g n ⟩ = ⊗ n k=1 g⟩ (k) , where g⟩ (k) is the ground state of the particle k. Since the state g n ⟩ belongs to the symmetric subspace, the probability of finding the system in this subspace is subject to the bounds 1 ≥ p sym β0 ≥ p g β0 . Taking the low-temperature limit β 0 → ∞ gives p g β0 = 1 and thus p sym β0 = 1. Hence, by continuity, we have p sym β0 ≃ 1 for β 0 ≫ 1. Since the block-structure (39) of the state γ β0 is conserved by the master equation (33), so is the probability p sym β0 . Consequently, the steady state of the ensemble is also nearly confined to the symmetric subspace and the thermodynamical quantities (48)-(50) are dominated by contributions from the symmetric irrep λ sym = (n, 0, . . . , 0). This observation makes it possible to derive asymptotically exact expressions for these quantities in the limit β ≪ 1. That is, we consider a situation where the system is initially thermalised in a low-temperature environment and then brought into contact with a high-temperature bath.
For an SU(d)-ensemble with single-particle Hamiltonian h = diag[ε 1 , . . . , ε d ], the partial partition function Z sym β of the symmetric subspaces is given by as we prove in App. B. Here, is the dimension of the symmetric subspace, which can be obtained from Eq. (19), and is the variance of the single-particle energy at infinite temperature. Note that we have ⟨ε⟩ = 1 d ∑ d i=1 ε i = 0, since the single-particle Hamiltonian is traceless by assumption. Using the expansions (59) and (66), we find that the quantities (48)-(50) become in the limit β 0 → ∞ with We now consider an ensemble of spin-s particles. By a similar argument to above, since the ground state belongs to the subspace of maximal angular momentum J = ns, the thermal state of such an ensemble will almost exclusively occupy this subspace in the limit β 0 ≫ 1. The dimension of this subspace is given by d max = 2ns+1, which is strictly smaller than the dimension of the symmetric subspace for s > 1 2 and n > 1. With the single-particle Hamiltonian h = S s z , the corresponding partial partition function becomes [20] Z ns β = (1 + 2ns) 1 + β 2 ns(1 + ns) as can be seen by expanding Eq. (53). The asymptotic expressions for the internal energy, reduced entropy and reduced free energy can now be obtained by replacing d sym with d max and q d n with Q s n = ns(1 + ns) 6 = n(d − 1)(2 + n(d − 1)) 24 (67) in Eqs. (62)- (64). Here, we have expressed s = (d − 1) 2 in terms of the dimension d of the single-particle Hilbert space to facilitate the comparison between SU(d)-and spin-s ensembles.
The expressions (65) and (67) show that, for the SU(d) ensemble as well as for the spin-s ensemble, the leadingorder corrections in β to the internal energy, reduced entropy and reduced free energy of the steady state all feature a quadratic dependence on n for all d. For a quantitative comparison, we choose a ladder Hamiltonian with eigenvalues ε i = −(d+1) 2+i for the SU(d)-ensemble, which yields ⟪ε 2 ⟫ = (d − 1)(d + 1) 12. We then have for n ≫ 1. Hence, in the limit of many particles, the magnitudes of the first order corrections in Eqs. (62)-(64) are suppressed by a factor (d−1) in SU(d)-ensemble compared to the spin-s-ensemble. We stress that this effect arises solely from the different dimensions of the primarily occupied subspaces, which have multiplicity 1 in both cases.

High initial temperature and many particles
A second interesting limit is realised when the ensemble is initially prepared in a high-temperature state, that is, β 0 → 0. In the following, we calculate the steady-state energy and reduced entropy in this limit, respectively for thermalisation with a low and a high-temperature bath. To this end, we first observe that the probability distribution p λ β0 tends to for β 0 → 0. This quantity is known as a Planchereltype measure [53]. In the limit of many particles, this measure can be determined explicitly through methods of asymptotic representation theory. Specifically, upon changing variables from λ i to ζ i = (λ i − n d) √ n, the limiting measure for n → ∞ becomes the function [54,55] One interesting question is how far the mean energy of an initially hot ensemble can be reduced by thermalising with a cold bath. When the bath temperature is low, i.e., if β → ∞, the partial thermal state γ λ β becomes a projector on the ground state within the unitary subspace H λ . As we show in App. C, each subspace has a unique ground state with energy ∑ d i=1 λ i ε i , where ε i are the single-particle energies ordered increasingly. Here, we assume that ε 1 is non-degenerate. Putting this observation together with the limiting distribution (70), we obtain the mean energy for n → ∞. The integral in this expression runs over all ζ 1 , . . . , ζ d ∈ R subject to the constraints ζ 1 ≥ ζ 2 ≥ ⋯ ≥ ζ d and ∑ d i=1 ζ i = 0. If we choose a ladder Hamiltonian with energy levels ε i = −(d+1) 2+i for the single-particle system, this mean energy takes the form for n → ∞. Notably, this result shows that the steadystate energy of the ensemble is sub-extensive in the particle number n. Hence, thermalisation is strongly inhibited by the symmetry constraints of the system; if the system were to instead fully thermalise to its ground state, its steady-state mean energy would be trivially given by For d = 2 and d = 3, the constant E d in Eq. (72) takes the values E 2 = 2 π and E 3 = 9 3 16π, see App. D for details. For larger dimensions, E d must be determined numerically by solving the integral in Eq. (71). These results are plotted in Fig. 3, which shows a non-linear increase of E d with d for d ≤ 7. For completeness, we may again compare the SU(3)-ensemble with a spin-1 ensemble, for which we find E ∞,0 → − 16 3π √ n, see App. E. Thus, low-temperature thermalisation through the 3-channel mechanism of the SU(3)-case leads to a lower steady-state energy than single-channel mechanism of the spin-1 case by a factor of 27 16 ≈ 1.69. In line with our earlier results, this observation quite naturally suggests that a larger number of dissipation channels allows the system to relax closer to a thermal state. Note however, that the scaling of the steady-state energy with n is insensitive to the dissipation mechanism. We now assume that the bath temperature is also high, i.e., β → 0. The ensemble then remains in a hightemperature equilibrium state, whose reduced entropy is given byS By using the distribution (70) for the many-particle limit, we prove in App. D that (74) to leading order in n for n → ∞. Thus, for d = 3, we findS 0,0 → (3 2) ln[n]. In comparison, for a spin-1 ensemble, the reduced entropy is suppressed by a factor of 3, that is we haveS 0,0 → (1 2) ln[n] to leading order in n. More generally, for a spin-s ensemble, it is clear thatS 0,0 is upper-bounded by the contribution S(γ ns β ) = ln[2ns + 1] ≃ ln[n] from the maximal angular momentum subspace corresponding to J = ns, cf. Eq. (49). It follows that, compared to spin ensembles, the reduced entropyS 0,0 of SU(d)-ensembles is asymptotically enhanced by a factor that is at least quadratic in the dimension of the single-particle Hilbert space.

V. ENGINE CYCLES
So far we have analysed how steady-state thermodynamical quantities are affected by the symmetries of permutation-invariant ensembles of n indistinguishable and non-interacting particles. As a next step, we explore in this section what role these symmetries play in thermodynamical processes. To this end, we consider a thermodynamical engine cycle that uses a collective working medium and investigate how the properties of this medium affect its performance.

A. Protocol and output
For simplicity, we focus on the standard quantum Otto cycle [56,57]. We begin by stating the 0+4 strokes that the system undergoes during this cycle. In the zeroth stroke, the ensemble is initialised in the Gibbs state at the inverse temperature β 0 , which determines the occupation probabilities p λ β0 of the individual irreps, cf. Eq. (39). The ensemble then cyclically undergoes the four strokes of the Otto cycle between a hot and cold bath at inverse temperatures β h and β c : (1) Equilibration with the hot bath, The net extracted work per cycle W , which is our main quantity of interest, is given by the sum of the work contributions from the two instantaneous strokes, Here, ρ β h ,β0 and ρ ′ βc,β0 are the steady states of the system with Hamiltonian H and H ′ , respectively, see Eq. (40).
We now assume that the instantaneous strokes are simple compression, that is H ′ = κH, where the compression factor κ must obey the inequalities β h β c ≤ κ ≤ 1 to ensure that a positive amount of work is generated. Under these conditions, the net work extraction with a collective medium becomes with the internal energies given by Eq. (48). To uncover the role of collective effects, we compare this quantity with the work generated in the same cycle with an ensemble of distinguishable particles, which fully thermalises during the isochoric strokes,

B. Collective work enhancement
Two natural questions arise at this point. First, whether it is possible to extract more work from the collective medium than from the distinguishable one, that is, can the ratio W col W dis become larger than 1? Second, whether there is an advantage in moving from spinensembles to ensembles with higher-order symmetry.
Both of these questions can be answered with the help of Fig. 4, where we plot the ratio W col W dis against β h for SU(2)-and SU(3)-symmetric ensembles. To enable a quantitative comparison, we have chosen the same singleparticle Hamiltonian h = diag[1, 0, −1] for both cases. We find that, first, there is indeed a temperature range for which W col W dis > 1, meaning that more work can be extracted with the collective medium than with the distinguishable one; a similar result was found in Ref. [20] for spin-ensembles. Second, Fig. 4(b) shows a regime where W col W dis > 1 while the extracted work is larger for the SU(3)-than for the spin-ensemble. An advantage of SU(3) over SU(2) is also seen in Fig. 4(c), where κβ c ≃ β 0 . Here, however, W col W dis < 1, so the higherorder symmetry group rather mitigates the disadvantageous collective effects.
The first observation can be understood in the limiting regime where the initial temperature is low and both baths temperatures are high, that is, β 0 ≫ 1 and β h , κβ c ≪ 1. Under these conditions, we can use the asymptotic result (62) for the steady-state internal energy of an SU(d)-ensemble to evaluate the work extraction from the collective medium. Using the formulas (76) and (77), we find This result shows that the work advantage grows monotonically in the number of particles n, but decreases monotonically with the order of the symmetry group d for n > 1. For the special case n = 10 and d = 3, we recover the limiting value W col W dis ≃ 3.25 seen in Fig 4(a). For comparison, the analogous result for spin-s-ensembles derived in Ref. [20] is

C. Optimisation
Higher-dimensional symmetry groups offer larger freedom in modeling the single-particle Hamiltonian as additional diagonal generators become available in the corresponding Lie algebra, cf. Eq. (29). We will now explore how this freedom can be exploited to optimise the work output of a collective quantum Otto cycle. To this end, we let the single-particle energies ε i vary within a bounded range, which we choose as ε max − ε min = d − 1 to match a ladder Hamiltonian with unit spacing.
By using the expansion (62), we can find the optimal spectrum in the limit of low initial temperatures, β 0 → ∞, and high bath temperatures, β h , κβ c ≪ 1. As shown in App. B, for even d, the optimal spectrum has two levels, which are d 2-fold degenerate with the maximal gap of d − 1; that is, we have ε 1 , . . . ε d 2 = −(d − 1) 2 and ε d 2+1 , . . . , ε d = (d − 1) 2. For this spectrum, the optimal work output of a collective cycle and the work output generated by an ensemble of distinguishable particles are Hence, the work-advantage ratio is still given by Eq. (78). For odd d, there is a small correction, since the remaining energy level can be either at 0 or d − 1.
For arbitrary temperatures, we approach the optimisation problem numerically for SU(3)-ensembles. Normalising the difference between ground and excited singleparticle energies to 2, we can parameterise the singleparticle Hamiltonian by the gap δ between the ground and first excited levels. That is, we make the ansatz Here, we have identified the diagonal generators with the Gell-Mann matrices (7), i.e., d 1 = Λ 3 and d 2 = Λ 8 . As shown in Fig. 5, the maximum work output is attained with a doubly degenerate excited level, corresponding to δ = 2. Notably, the collective medium outperforms the distinguishable one for any value of δ.
Our results show that the larger freedom in the choice of the single-particle Hamiltonian, which comes with higher-order symmetry groups, can indeed be used to optimise the performance of many-body Otto cycles, in addition to the boost coming from a collective dissipation mechanism. This can be done, at least in principle, even outside the regime of extreme initial and bath temperatures. Furthermore, our analysis seems to suggest that the optimal single-particle Hamiltonian for collective quantum Otto cycles should feature two levels with approximately balanced degeneracy. Corroborating this presumption with a more systematic investigation is beyond the scope of our present work, but provides an attractive problem for future research.

VI. HIGHER-ORDER SYMMETRIES: A CLOSER LOOK
As we have seen in the previous sections, the general tools developed in this paper can be used to calculate thermodynamical quantities for non-interacting spin-s ensembles and ensembles with arbitrary SU(d)symmetry. Beyond these applications, our framework also makes it possible to analyse the anatomy of the steady states of such ensembles on the level of individual irreps.

A. Free energy and dimensions
To illustrate this approach, which can provide further insights into the role of higher-order symmetries, we first consider the change in the reduced non-equilibrium free energy as defined in Eq. (58). For a quantitative comparison between the different symmetry groups, we choose the coefficients a i in the single-particle Hamiltonian h = ∑ i a i d i to reproduce the ladder spectrum ε i = −(d + 1) 2 + i. With this choice we first plot the non-equilibrium free energy change as function of the inverse bath temperature β in Fig. 6. We see that, for all bath temperatures, ∆F increases with the order d of the applied symmetry group.
To understand this effect better, we now focus on the role of the individual irreps. In Fig. 7, we decompose ∆F β,β0 = ∑ λ p λ β0 ∆F λ β,β0 into contributions from the individual irreps λ; the size of each dot in these plots indicates the relative probability of occupying a particular irrep. For spin ensembles, the maximal angular momentum irrep J = ns, which has the largest dimension d max = 2ns + 1, is the predominantly occupied irrep in the limit β 0 → ∞. For SU(d), the symmetric irrep, which does not generally have the largest dimension, is most occupied for large β 0 . Thus, for spin ensembles, the dominant contribution to the free energy change can only shift to a lower-dimensional irrep as β 0 decreases to a moderate value -see Fig. 7. By contrast, for SU(d)ensembles with d > 2, one observes a shift towards a higher-dimensional irrep.

B. Degeneracies
An explanation for the link between thermodynamical properties and irrep dimension can be given in terms of energy level degeneracies. Here, one finds a notable difference between SU(2) and higher-order groups. SU (2) irreps are special in having non-degenerate energy levels -that is, for each J, the energy basis J, M ⟩ varying over M is non-degenerate. For SU(3) and above, we observe two sources of degeneracy. The first is that different sets of occupation numbers n i of the single-particle energy levels ε i may have the same total energy. This can only happen when the ε i are rationally dependent. For example, with the three-level ladder where (ε i ) i = (−1, 0, 1), the occupation numbers (0, 2, 1) and (1, 0, 2) represent three-particle configurations with the same energy of +1. The ability of the higher-order group dynamics to mix between such configurations results from the greater number of dissipation channels, compared with a single channel for SU (2).
The second type of degeneracy occurs within a given set of occupation numbers, and can be seen by considering Young tableaux. Referring to the construction described in Sec. II D of a basis for the SU(d) irreps, we consider the example d = 3, n = 3 and the irrep for λ = (2, 1, 0). Two of the possible Young tableaux as- FIG. 7. Contributions to the free energy change ∆F λ β,β 0 for each irrep λ, such that ∆F β,β 0 = ∑ λ p λ β 0 ∆F λ β,β 0 , as functions of the irrep dimension d λ . The radius of each dot is scaled in proportion with the corresponding occupation probability p λ β 0 . In the first plot, which corresponds to an SU(2)-ensemble with spin s = 1, the maximal J irrep is circled. In the middle and bottom plots, which correspond to SU(3)-and SU (6) These Young tableaux describe two linearly independent, though non-orthogonal, states in the subspace H λ . By orthogonalising, one obtains a pair of basis states with the same occupation numbers (n i ) i = (1, 1, 1) but which are still distinguishable according to the permutation-invariant system dynamics. It is not hard to see that such situations can only occur with n > 2 and d > 2. Moreover, only symmetries other than fully symmetric and antisymmetric can support this degeneracy -with only a single row or column in the Young diagram, the filling order is fixed. This is why (as in Fig. 7) the symmetric irrep is not the one with highest dimension, other than for SU (2). This second type of degeneracy is intrinsically nonclassical in the sense that states of classical indistinguishable particles are labelled only by their occupation numbers. Thus, for more than two particles and dynamics generating a higher-order symmetry group, the effective state space contains degrees of freedom that do not exist in ensembles of classical indistinguishable particles. Notably, the same mechanism underlies the quantum modification recently found in the well-known Gibbs paradox [11].
Are these additional degeneracies useful? Whereas the ground and fully excited states are non-degenerate (see App. C), the degeneracies tend to increase in the middle of the spectrum. The work output of the Otto cycle in Sec. V is dictated by the heat capacity C = V (E) T 2 , where V (E) is the variance of the energy [21]. This can be seen when κβ c is close to β h by relating the energy difference E β h ,β0 −E κβc,β0 to its derivative with respect to bath temperature, using C = ∂E ∂T . Thus the work output is proportional to the energy variance. For high bath temperatures T c and T h , the populations are roughly uniform over the energy spectrum, and so the advantage gained by the SU(2) system is due to its suppression of degeneracies in the middle of the spectrum, thus raising the energy variance. In contrast, for low bath temperatures, most of the population is in the ground state, so raising the degeneracy of higher energy states will increase the variance. Therefore distinguishable particles will perform better in this case. For an intermediate temperature range, one expects SU(3) to perform best, since it has degeneracies at higher energies, though not as many as for distinguishable particles. These observations thus qualitatively explain the behaviour in Fig. 4.
It remains to be explored whether the additional nonclassical degrees of freedom can be further exploited to gain an additional advantage in some setting. For example, one could engineer a Hamiltonian that breaks these degeneracies in order to extract work from them. Unlike superradiance, which has been argued to exist in classical models [58,59], this effect would have no classical analogue.

VII. PERSPECTIVES
In this paper, we have combined methods from representation theory to build a comprehensive theoretical framework for the thermodynamical description of non-interacting quantum many-body systems with permutation-invariant observables. As our main applica-tion of this formalism, we have investigated the structure and properties of steady states that emerge when such systems are weakly coupled to a thermal environment. We have further shown that permutation invariance induces collective effects that can in principle be used to enhance the performance of quantum thermal machines. At every step of our analysis, we have demonstrated that altering the constituents of a permutation-invariant ensemble from spin systems, which are covered by the conventional Clebsch-Gordan theory, to more general multilevel systems, for which we have developed a systematic description in this paper, can lead to qualitative changes in the thermodynamical properties of these ensembles.
Two possible extensions of our approach appear to be promising perspectives for future research. First, it would be interesting to include oscillating driving fields. For fast driving, this extension can be achieved by replacing steady-state master equations with thermodynamically consistent Floquet-Lindblad equations [35,60], which have been used earlier to study permutation-invariant spin systems [19]. In the slowdriving regime, one can instead employ adiabatic master equations, whose generators are obtained by replacing time-independent system parameters with external control protocols [36,[61][62][63]. In this way, it would in particular be possible to explore the role of permutationinvariance in the context of thermodynamic geometry, a topic that is currently attracting much interest in quantum thermodynamics [64][65][66][67][68][69][70][71].
From a mathematical perspective, progress in these directions would require an extension of our theorem on the structure of steady states arising from autonomous permutation-invariant master equations to periodic limit cycles. Since the generators of Floquet-Lindblad equations become time independent in a rotating basis, this generalisation should be straightforward in the fastdriving regime. For the adiabatic limit, a counterpart of Spohn's theorem [72] for steady states that provides uniqueness conditions for periodic limit cycles is available [73]. This result will however still have to be generalised to permutation-invariant many-body systems. Similarly, it is important to generalise the setting to regimes beyond weak coupling, involving techniques such as stochastic Liouville equations [74] and the thermal leads approach [75]. Permutation symmetry can in principle be included in such settings, so the main theoretical challenge is to find simple models under which the steady state is the same partially thermalised state analysed here.
A second key problem is to systematically investigate the role of interactions. Most of our results, in particular our theorem on steady states, depend only on permutation invariance and should therefore be applicable also to interacting systems. However, since interactions typically lead to a dense energy spectrum, they render the rotating wave approximation that underpins the conventional weak-coupling master equation invalid. It would thus be necessary to either identify a relevant class of in-teracting systems that still feature a spectrum with wellseparated Bohr frequencies, or to employ new types of recently derived quantum master equations [76][77][78], whose thermodynamical consistency is however not yet settled. One class of interacting Hamiltonians which may be most tractable is those with a product state energy basis -such as pairwise d i d j interactions -since the eigenvectors of H λ can be taken as the Schur basis.
From an information-theoretical viewpoint, one extension could involve studying the full work statistics of the Otto cycle rather than just the mean -due to the conservation law relating to the symmetry type, each block has its own associated probability and work output, so one is forced to consider the work as a fluctuating random variable. The thermodynamical quantities considered here would then need to be replaced by single-shot versions [79]. For example: what is the minimal work output that is guaranteed per cycle with some small probability of failure? Another line of questioning is what our results have to say about the thermodynamical processing of ensembles beyond the "independent and identically distributed" (i.i.d.) regime, where one typically considers many independent copies of a state subject to arbitrary operations. This setting has been used to collapse the many quantum "second laws" [80] down to just the standard non-equilibrium free energy [81]. Perhaps collective operations as defined here may give rise to different laws.

DATA AVAILABILITY
All plotted data were generated from the equations discussed in the text. No further data were created by the research presented in this paper.
The main result of this appendix is Theorem 1. This theorem specifies the structure of the steady states of a class of collective quantum master equations, which is more general than the one discussed in the main text.
We consider a permutation-invariant ensemble of n d-level particles embedded in a large, not necessarily thermal, bath. The bare ensemble is described by a permutation-invariant Hamiltonian H, which may in principle also contain interactions between particles. We further assume that the system-bath interaction is described by the Hamiltonian where the E µ are the non-diagonal generators of the tensor-product representation of SU(d) on the ensemble Hilbert space and the Hermitian operators B µ = B −µ correspond to bath observables. The operators E µ can now be decomposed into Bohr-frequency components with respect to the ensemble Hamiltonian H, that is, Provided that the conventional weak-coupling, Born-Markov and rotating wave approximations are applicable, the time evolution of the state ρ of the ensemble can then be described in terms of the master equation ∂ t ρ t = L(ρ t ) with the Lindblad generator L being defined in Eq. (A4).
Before we can move on to the proof of our main theorem, we need a lemma that characterises those operators that commute with all E µ . Note that no information is needed about the dimensions of multiplicity spaces K λ , so that the results below apply not just when Schur-Weyl duality is valid, but also to spin systems upon replacing λ by J. Lemma 1. Let {µ} be a set of simple roots of SU(d). Then, for an operator X it holds that where X λ is an arbitrary operator on the multiplicity space K λ .
Proof. The "if" part is obvious. For the "only if" part, we make use of the fact that the simple roots plus their negatives generate the whole root system (see Ref. [82], Eq. 21.20). In other words, for any root ν, there exists a sequence of simple roots or their negatives, µ (1) , µ (2) , . . . , µ (k) (allowing for repetitions) such that If [X, E µ (i) ] = 0 ∀i, we therefore see that [X, E ν ] = 0, as E ν is constructed from products of operators commuting with X.
having used the Jacobi identity in the second line. Since the d − 1 simple roots are a linearly independent set, (A3) can only hold for all µ if all commutators [X, D i ] vanish. Therefore we have found that X must commute with the representation of SU(d).
Finally, we decompose the unitary representation as U = ⊕ λ 1 K λ ⊗ U λ . Considering any block of X λ,λ ′ with respect to the irrep structure, commutation implies Schur's lemma then tells us that the off-diagonal blocks vanish, while the diagonal ones are of the form X λ,λ = X λ ⊗ 1 H λ .

Theorem 1. Consider the Lindblad generator
where the E µ are the collective representations of the non-diagonal generators of the Lie algebra su(d) in the Cartan basis, H is any permutation-invariant Hamiltonian and the coefficients Γ µν ω form a Hermitian matrix for each Bohr frequency ω. Further assume that there exists a set of simple roots R such that, for each ω, the coefficients Γ µν ω form a matrix that has full rank over those roots µ ∈ R for which the frequency components E µ (ω) do not vanish. Then the steady states of the generator L are all of the form where the probabilities p λ and the operators σ λ are arbitrary, while the ρ λ are unique on each irrep λ.
Proof. We adapt the original proof by Spohn that guarantees a single steady state [72], more closely following the presentation of this proof in Ref. [83]. First, we show that a set of ρ λ exists such that all states of the form (A5) are steady states. Due to the permutation symmetry of the dynamics, any such state must evolve in time as ρ t = ⊕ λ p λ σ λ ⊗ ρ λ t . Since the irreps H λ are finite-dimensional, each contains at least one steady state ρ λ .
Next, we still have to prove that this family covers all the steady states. Initially, we relabel the set {E µ (ω)} µ,ω as {E i } i such that i corresponds to a pair (µ i , ω i ). Next we put the generator into diagonal form by transforming to F i = ∑ j u ji E j , where the u ji form a unitary matrix. Then we have E i = ∑ j u * ij F j , so with A kl = ∑ i,j u * ik δ ωi,ωj Γ µiµj ωi u jl . Choosing u to diagonalise the Hermitian matrix A, we obtain We would like to guarantee that a k > 0. To this end, we note that A is unitarily equivalent to the matrix Γ = ⊕ ω Γ ω ∶= ⊕ ω µ,ν∶ Eµ(ω),Eν (ω)≠0 Γ µν ω ω, µ⟩⟨ω, ν .
Note that we use Dirac notation here for convenience although the vectors ω, µ⟩ do not represent quantum states. The aim of this construction is to include only those terms for which E µ has a nonzero frequency component E µ (ω); otherwise, the full-rank condition would not generally be possible to satisfy. Thus we see that the positivity of A is guaranteed by assuming that all Γ ω have full rank. We then make a second transformation by expanding F i = ∑ p i=1 f ik G k in some operator basis such that the G k are Hermitian and the f ik are complex. The G k are chosen to form a basis of the subspace V = span{F i } i = span{E µ } µ with some dimension p. This substitution results in where B kl = ∑ i a i f ik f * ik . Evidently B = ∑ i a i f i f † i > 0, which means we can find some b > 0 smaller than the smallest eigenvalue of B such that B − b1 p > 0. Using this result, we divide L into two parts: such that both parts are valid Lindblad generators. We first analyse the spectral properties of L 2 . For any Hermitian operator X, we calculate the super-operator expectation value This quantity is generally negative and vanishes if and only if [G k , X] = 0 ∀i = 1, . . . , p. This condition is equivalent to [E µ , X] = 0 ∀µ, and by Lemma 1, also equivalent to X = ⊕ λ X λ ⊗ 1 H λ . Let us denote the subspace of such operators by S = {⊕ λ X λ ⊗ 1 H λ X λ † = X λ }; this space has dimension M = ∑ λ m λ (m λ + 1) 2. We have already found a space of steady states with this dimension, so our goal is to show that L has precisely M linearly independent zero eigenvectors.
Next we divide the vector space of operators into two orthogonal parts S ⊕ S ′ . There is then a block decomposition of L on a suitable basis respecting this structure: representing the super-operator as a matrix, and similarly for L 1,2 , where * denotes an unspecified block and L ′ is the projection onto S ′ . From (A12) we see that ⟨L 2 ⟩ X < 0 for all X ∈ S ′ , so the eigenvalues of L ′ 2 all have strictly negative real part. Similarly, for all X ∈ S ′ , since L ′ 1 must have eigenvalues with non-positive real part in order to be a valid Lindblad generator. Therefore L ′ also has eigenvalues with strictly negative real part.
Finally, to learn about the eigenvalues of L, we use Theorem 1.4.10 of Ref. [84]. This theorem tells us that if an eigenvalue λ of some n × n matrix has geometric multiplicity of at least k + 1, then every (n − k) × (n − k) principal submatrix also has λ as an eigenvalue. Suppose for a contradiction that L has more than M zero eigenvectors. So we apply this theorem with λ = 0, k = M and n as the dimension of the full operator space. This observation implies that L ′ has a zero eigenvector, which contradicts what we have just determined about the eigenvalues of L ′ . Therefore L has at most M zero eigenvectors, and we have already found this many.
With a thermal environment, one can assume a standard detailed-balance relation on the bath correlation function and thereby obtain the following result guaranteeing thermalisation within each block: Corollary 1. With the same assumptions as in Theorem 1, together with Γ µν −ω = e β ̵ hω Γ νµ ω , the steady states are all those of the form where γ λ = e −βH λ tr[e −βH λ ] is the thermal state on each irrep, according to the Hamiltonian H = ⊕ λ 1 K λ ⊗ H λ .
Proof. This result is an immediate consequence of Theorem 1 followed by applying, to each diagonal block, Spohn's result [72] that the condition Γ µν −ω = e β ̵ hω Γ νµ ω requires the steady state ρ λ to be the thermal state γ λ .
The mean energy is having used the fact that the uniform distribution of energies is symmetric with respect to permutations of the energy levels. Similarly, = ⟨n 2 1 ⟩ n d ⟨ε 2 ⟩ 1 + ⟨n 1 n 2 ⟩ n d 2 ⟨ε⟩ 2 1 − d ⟨ε 2 ⟩ 1 .
Given that ∑ i n i = n, we have The probability distribution over n 1 is found by counting the number of ways of distributing the remaining n − n 1 particles over the other d − 1 energy levels: from which we get This is to be compared with the full (infinite temperature) thermal state of n particles, which has energy variance n⟪ε 2 ⟫ 1 .
We can now ask what set of single-particle energies gives the maximal variance; due to Eq. (B10), we only need to consider the case of a single particle. In order to make a well-posed question, we constrain the ε i such that there is at least one at the minimum energy ε 1 = 0 and at least one at the maximum ε d = d − 1. (These values are chosen to match the spread of the ladder Hamiltonian -note that the overall shift can be arbitrary.) The variance must always satisfy For even d, this can be saturated by taking half the levels at 0 and half at d − 1. For odd d, the optimum is with either d−1 2 at 0 and d+1 2 at d − 1, or the opposite. Hence, max {εi}, ε1=0, ε d =d−1 , d even, For the evenly spaced ladder, one instead has ⟪ε 2 ⟫ 1 = d 2 −1 12 -so the level optimisation roughly gains a factor of three.