Symmetries and physically realizable ensembles for open quantum systems

A $D$-dimensional Markovian open quantum system will undergo stochastic evolution which preserves pure states, if one monitors without loss of information the bath to which it is coupled. If a finite ensemble of pure states satisfies a particular set of constraint equations then it is possible to perform the monitoring in such a way that the (discontinuous) trajectory of the conditioned system state is, at all long times, restricted to those pure states. Finding these physically realizable ensembles (PREs) is typically very difficult, even numerically, when the system dimension is larger than 2. In this paper, we develop symmetry-based techniques that potentially greatly reduce the difficulty of finding a subset of all possible PREs. The two dynamical symmetries considered are an invariant subspace and a Wigner symmetry. An analysis of previously known PREs using the developed techniques provides us with new insights and lays the foundation for future studies of higher dimensional systems.


Introduction
Given an open quantum system obeying a Markovian master equation (ME), an experimentalist typically has much choice as to the way in which the environment, with which the system interacts, is monitored. Given this experimental freedom, one can ask the question of whether the system state can be made, in the long time limit, to jump between a finite number of pure states, with static dynamics between jumps. That is, can the evolution of the system state be made to follow a Bohr-Einstein model of jumps between quantum states that are, conditional upon measurement results, stationary (although not necessarily the eigenstates of some Hamiltonian)? If it can, how many such states, at a minimum, are required? Alternatively, what is the minimum Shannon entropy of the ensemble? Questions of this type [1] may offer novel insights into the nature of open quantum system dynamics and the quantum-classical epistemic distinction [2,3,4].
It might have been thought that since there is an infinity of ways to write a state matrix as an ensemble of, in general, non-orthogonal pure states, there would also be an equivalent diversity of ensembles that are realizable via the application of an appropriate measurement scheme. However, this has been shown to be incorrect [5]: for some ensembles there is no way that the experimentalist can know at all (long) times that the system is in one of the states belonging to that ensemble. Ensembles that are realizable are known as physically realizable ensembles (PREs). Such ensembles support an ignorance interpretation in the sense that an individual who is not privy to the experimental results, produced during the ensembles realization, could claim that the system really is in one of the pure states belonging to the PRE.
In Ref. [5], a set of polynomial constraint equations were discovered that govern the existence (or otherwise) of PREs. A solution set of the constraints would describe all possible PREs for a given ME. These constraints will be detailed in the next section, but in our introductory remarks we wish to discuss their generic solubility. Unfortunately, solving a set of nonlinear polynomial equations typically becomes exponentially more difficult as the number of equations and variables increases. In fact, the problem is known to fall into the NP-complete complexity class [6]. This difficulty is acute in the case of the PRE constraint equations in question, as the number of constraints in the set scales as KD 2 , where K is the number of pure states in the ensemble and D is the system dimension. Thus far, a full description of PREs for a given ME -via analytic or numeric means -only exists in the literature for K = D = 2 [1].
In this paper, we introduce methods that can simplify the task of finding example PREs for a ME. This is important as it will make tractable the study of PREs for a broader range of MEs and ensembles, in particular MEs for systems of larger dimension, D, and also PREs containing more members, K. The primary tool that we introduce is that of symmetry; many physical systems of interest possess symmetries and even generic systems possess some symmetry (by our definitions, as will be made clear in the appropriate basis). In Ref. [7], it was stated that PREs of smaller K (than otherwise expected) could perhaps be discovered by exploiting additional structure within the ME: in this paper we develop such ideas in considerable detail, with the conclusion that there is, indeed, much to be gained from an analysis that utilizes system symmetries. Due to the, still considerable, computational difficulty of finding PREs in D > 2, we defer these searches to future work. In the current paper, we are satisfied with an exploration of previously analyzed MEs and PREs from a new perspective that gives increased understanding and, also, experience with our new techniques.
Two types of symmetry are explored, followed by a consideration of their joint application. Firstly, we consider an 'invariant subspace symmetry', by which we mean that the dynamics of the system state are contained to within some region of the space of density matrices, given an initialization within that space. Our definition contrasts with other well-known invariant subspace definitions [8] and allows us maintain the dimensionality, D, of the system but, nevertheless, decrease the size of the space that is searched for the presence of PREs. This reduces the computational task of identifying a PRE. Secondly, we consider 'Wigner symmetries', which are those transformations that preserve the inner product of Hilbert space rays [9]. The consequences of such a symmetry existing are explored for PREs: we find that new PREs can be generated via the application of a Wigner symmetry to an already known PRE and, also, that some fraction of the constraint equations become redundant given a PRE that is Wigner symmetric (a concept that will be defined). Both of these discoveries will make a study of PREs possible in some situations where their finding was previously intractable. Detailed case studies for these symmetries are made for two single-qubit MEs: the resonance-fluorescence ME and the thermal equilibrium (absorption and emission) ME.
The organization of this paper is as follows. Firstly, we provide a more mathematical description of PREs, in Sec. 2, in order to set the grounds for our in-depth study. In Sec. 3, we motivate the need for new techniques to find PREs, based on a discussion of the inherent computational difficulty. Next, in Sec. 4, we consider the invariant subspace symmetry, followed by a discussion of its utility and qubit examples. Then, in Sec. 5, we consider Wigner symmetries and how they can generate both new PREs and be applicable to an individual PRE. Once again we highlight with qubit examples. In Sec. 6, we consider how both these symmetries may co-exist. This is then followed by a conclusion and discussion in Sec. 7.

Mathematical description of PREs
This paper is concerned with autonomous Markovian open quantum systems of finite dimension D. In the absence of measurement (or if the measurement results are ignored) the dynamics of such systems are governed by a Lindblad-form master equation (ME) for the density matrix [10]: whereĤ eff ≡Ĥ − i lĉ † lĉ l /2 andĤ is the Hermitian Hamiltonian. Without loss of generality we take the Lindblad operators {ĉ l } as linearly independent and traceless but otherwise arbitrary, which leads to a bound of L ≤ D 2 − 1 [11]. The separation of Eq. (1) into terms that involveĤ eff and those comprisingĉ l ρĉ † l can be associated with a purity-preserving unraveling of the ME as would arise from perfectly efficient monitoring of the decoherence channels, l. Assuming an initially pure system state, if a detection in channel l is observed at time t then the system jumps from the pre-jump state |ψ(t − ) to the post-jump state |ψ(t) ∝ĉ l |ψ(t − ) . After the jump, the quantum state evolves under the no-jump evolution operatorĤ eff and will not remain stationary unless it happens to be an eigenstate ofĤ eff .
We limit the MEs under consideration to those that produce an impure, but unique, steady-state density matrix, ρ ss , of rank D, defined by Lρ ss = 0. It is always possible to write ρ ss as a mixture of pure states, but because the pure states in this decomposition are not necessarily orthogonal there are an infinite number of ways of doing so. However, only a subset of these ensembles -termed physically realizable ensembles (PREs)can be realized at all sufficiently long times as the states conditioned on the outcomes of some monitoring of the decoherence channels. This might seem obvious due to the requirement that the ensemble members be eigenstates of the no-jump evolution operatorĤ eff . However, that would be a premature judgement, as we need to consider adaptive detection, for which the operatorĤ eff can change, even while Eq. (1) remains fixed. This will be explained fully below.
More than one PRE can be achievable (in separate experiments of course) for a given ME due to the invariance of Eq. (1) under the following joint transformations [12,1] where, with M ≥ L, β is an arbitrary complex M -vector and S is an arbitrary M × L semi-unitary matrix; that is, M m=1 S * ml S ml = δ l,l . By unraveling Eq. (1) with {ĉ m } as the jump andĤ eff =Ĥ − i M m=1ĉ m †ĉ m /2 as the no-jump operators, different β and S thus correspond to different measurement schemes [13,10]. It is important to note that the implied number of detectors (equal to M ) can be greater than the number of Lindblad operators, L, describing the ME.
The above flexibility is most easily understood in a quantum optics context: linear interferometers1, described by S, take the field outputs of the system as inputs while β represents the addition of weak local oscillators (WLOs) to the interferometer outputs prior to photodetection. We are interested in the case where the LOs are weak (that is, |β l | 2 not much larger than Tr[c † l c l ρ ss ]) so that a non-negligible amount of information is gathered concerning the system upon each detection. This means that the quantum trajectory is jump-like rather than diffusive; consequently, the PRE can consist of a finite number of states. We stress that most jump-like unravellings will lead, like diffusive unravellings, to PREs of infinite size. Those that lead to finite PREs, as we are interested in, are exceptional. It was shown in Ref. [5] that an ensemble {℘ k , |φ k } of size K is physically realizable iff there exist real valued transition rates κ jk ≥ 0 (which naturally determine the occupation probabilities ℘ k ) such that For each k ∈ 1...K, Eq. (4) is a matrix of constraints. In general, a ME will allow multiple solutions to Eq. (4) via the experimental freedom described by Eqs.
(2)-(3); see below for further discussion. Most of the difficulty of our research program arises due to the system of non-linear constraints defined by Eq. (4) being difficult to solve, even numerically, when D > 2. The transition rates κ jk play the role of free parameters in Eq. (4), and their possible values are determined by the range of solutions that exist. However, we can make a few observations, without directly solving Eq. (4), based on simple requirements of the PRE. Firstly, it is feasible that some of the κ jk can be zero only if this still allows the entire ensemble to be explored repeatedly in the long time limit. That is, by following the possible non-zero transitions, all ensemble members can be reached from every other ensemble member. A particularly simple arrangement of transitions that satisfies this is of a cyclic nature: the only non-zero transitions are of the form κ mod(k+1,K),k . This is not to say that cyclic PREs will always exist, merely that they do satisfy the principle just discussed. If further non-zero transition rates exist (or the pattern of transition rates is not based on a cycle), then we will refer to the PRE as non-cyclic. Many of the example PREs that we analyze will be of a cyclic nature.
Given a viable PRE [an ensemble {κ jk ≥ 0, |φ k } satisfying Eq. (4)], there must exist an appropriately applied measurement scheme such that the conditioned state of the quantum system will, in the long-time limit, jump between the ensemble members, spending a time in the state |φ k proportional to ℘ k . The measurement scheme in question will, in general, be adaptive, meaning that the experimental setting parameters ( β, S), must be changed according to which state (k) the system is currently in. Adaptive measurements by a controllable local oscillator were first studied in the context of quantum jumps in Ref. [14]. They are more commonly studied for their use in state discrimination [15] and phase estimation [16,17], with the latter having been realized experimentally [18,19,20]. In the body of our paper, we focus on characterizing PREs for MEs possessing a symmetry, and defer to the appendices an analysis of the measurement schemes required to produce the PREs. The reason for this is that finding the measurement scheme for a given PRE is a comparatively simple computational task, whereas we are concerned with ameliorating the difficulty of finding PREs in the first place.

Feasibility of finding PREs
As this paper's purpose is to make easier the solution of Eq. (4), it is important for us to motivate this as being necessary. Further details will be provided in Sec. 4.2, but let us here consider the difficulty of finding PREs in D = 3, which is the smallest dimension for which it is not possible, in general, to construct an example PRE via analytic means. For D = 3, it is argued, in Ref. [1], that the minimum expected PRE size is K = 5. Thus, Eq. (4) presents as 5 coupled matrix equations, each of which represents D 2 = 9 constraints, giving 45 total constraints. The RHS of Eq. (4) consists of terms cubic in the variables describing the transition rates (linear contribution) and the pure state vectors (quadratic contribution).
How difficult is it to solve a system of 45 coupled, mostly cubic, polynomials? The answer is that it is very difficult, at least as judged by current computational resources. In Ref. [6] it is suggested that even when using the highly efficient Gröbner basis Faugère F4 algorithm [21] (a variant on the Buchberger algorithm [22]) quadratic systems of size larger than 15 equations are very difficult. The constraints of Eq. (4) are cubic (except for K quadratic normalisation constraints) and are, resultantly, even more difficult. As another example, the computational algebra software MAGMA [23] has timings [24] provided for several benchmark problems, with the Cyclic 9 problem (see [25] and references therein) over a rational field listed as being solved in 4.6 days. The Cyclic 9 problem has only 9 equations but has a maximum degree of 9. Our experience with Gröbner basis techniques using MAGMA was that an example system of 16 equations was soluble in about 13hrs [26] whereas size 20 systems were out of reach (limitations of 200Gb of Ram were exceeded after several days runtime).
In the current paper, we defer the numerical task of solving Eq. (4) for D > 2 and, instead, after acknowledging its difficulty, focus on developing techniques that will allow such problems to be tackled in future work (where we will also provide more details on the available numerical methods). To test our new methods, we consider in detail mostly D = 2 examples and view these from the new perspective provided by an analysis of the available symmetry.

Definitions
Before defining what is meant by an L-invariant subspace, we introduce some notation. Let B (H) represent the set of bounded linear operators on the Hilbert space H. Density matrices, ρ, being trace-one positive-semidefinite operators, form a convex set D (H) ⊂ B (H), with extreme points in the set corresponding to one-dimensional projectors, |ψ ψ|, that represent pure states. An operator basis for ρ, sayσ σ σ = {σ j } D 2 j=1 , can be chosen that consists of D 2 orthogonal Hermitian matrices, D 2 − 1 of which are traceless. The generalized Pauli matrices, suitably normalized, provide an example of one such basis [27]. In the chosen basis, the state of the system can then be represented by D 2 coordinates, comprising a real-valued vector, r r r = {r j } D 2 j=1 ∈ R D 2 , as follows: Commensurately, the Lindblad dynamics of Eq. (1) can be stated aṡ r r r = Lr r r, with L being a real, D 2 × D 2 matrix representation of the Lindbladian superoperator, L (L should not be confused with the number of decoherence channels, L). We can simplify further by setting the traceful operator of our basis to beσ D 2 = 1, as theṅ r D 2 = 0 and r D 2 = 1. The dynamics of Eq. (6) reduces to the D 2 − 1 dimensional subspace (x x x = {r j } D 2 −1 j=1 , corresponding to the traceless operator basis components), givingẋ where L 0 is the restriction of the Lindbladian matrix representation, L, to the first which are called generalized Bloch vectors [29], occupy the region of R D 2 −1 that is mapped to D (H) via Eq. (5) with r D 2 = 1. This region is contained within a closed (D 2 − 1)-ball of radius D(D − 1)/2. (All geometric 'balls' -and 'discs' -in our work should be understood as 'closed'.) In D = 2, we have the familiar Bloch ball, every point of which is a valid density matrix. However, for D > 2, the inverse mapping of D (H), to the vectors x x x, gives only a convex compact subregion of the ball (that is still of dimension D 2 − 1). Moreover, for D > 2, the boundary of this convex subregion will comprise both extremal (pure) and non-extremal (mixed) states. Extremal states are those that cannot be formed from a combination (mixture) of other states in the region; a collection of such states is called an extremal set.
Recall that we assume L to have a unique steady state, ρ ss , of rank D, where Now let us translate to coordinates u u u = x x x − x x x ss , so that For these new coordinates we obtaiṅ u u u = L 0 u u u. (11) and the new steady-state, u u u ss , is at the origin. As the vectors u u u are merely translations of x x x, we will talk about them in the same geometric terms as the generalized Bloch vectors. The region of R D 2 −1 that they occupy will be referred to as U 0 , which is, in turn, mapped to D (H) via Eq. (10). The boundary of the space U 0 (as viewed in the new coordinate system) will lie within (or on) that of a sphere of radius D(D − 1)/2 translated by −x x x ss from the origin. Many of the example systems in this paper will have D = 2; in such cases, to make it clear which coordinate system is being discussed, we will use x, y, z as the components of x x x (which lies in the Bloch ball) and u, v, w as the components of u u u (which lies in U 0 ).
We define an L-invariant subspace D I ⊂ D (H) to be a convex space containing ρ ss , and having the property that the image of D I under e Lt for any t ≥ 0 is ⊆ D I 3.
Equivalently, this L-invariant subspace can be represented by a convex subregion I 0 ⊂ U 0 , comprising vectors of the form u u u, such that the image of I 0 under e L 0 t is contained within I 0 for t ≥ 0. That is, I 0 characterizes the invariant subspace and we will term the dimension of I 0 as the invariant subspace dimension.
For the invariant subspace definition to be an interesting one in the context of PREs, we require that D I corresponds to an I 0 that is an N -dimensional projective subspace of U 0 , with D − 1 ≤ N < D 2 − 1. The lower bound on dimensionality is necessary so that it is possible for I 0 to represent a PRE (which must contain D linearly independent pure states). The upper bound will allow a simplification in our search for PREs, due to the reduced dimensionality. That I 0 is a projective subspace indicates that it consists of lines through the origin; this ensures that the steady state is included (the origin in U 0 space) and that the extremal set of I 0 is in the boundary of U 0 . However, as discussed earlier, the boundary of U 0 is not necessarily pure, so we additionally require that I 0 has an extremal set in the extremal set of U 0 .
To span the entire density matrix space, which is represented by U 0 , requires the complement to I 0 , an orthogonal projective subspace of U 0 (we also include the origin) that is of dimension at most D 2 − D. We term this space R 0 and note that I 0 ∪ R 0 = U 0 and I 0 ∩R 0 = 0 (the origin, not the empty set, ∅). Let us now move to a new orthonormal basis for U 0 such that each basis vector lies wholly within either I 0 or R 0 . Additionally, we write the coordinates of u u u with those N degrees of freedom corresponding to I 0 first, 3 It is important to realize that our notion of invariance is different from that of Ref. [8]. We consider invariant subspaces of D (H), with dim(H) = D. We require that the invariant space contains rank D density matrices (meaning that there will exist density matrices that are formed by an ensemble of D linearly independent pure state projectors). This is in contrast to Ref. [8], where the invariant subspace, of Hilbert space dimensionD, is D(H) for H =H ⊕ H R (with, here, R being the remainder space). In other words, their invariant space consists of all possible density matrices with support solely in some Hilbert subsystem,H, whereas that considered in this paper is a compact subregion D I ⊂ D (H). Note that in the long time limit, under ME dynamics, the state necessarily relaxes to the unique ρ ss , a single point within the invariant subspace.
followed by those of R 0 . The purpose of this is to clearly expose the invariant subspace in the representation L 0 , which will take on the block form where L I 0 and L R 0 are square with length N and (D 2 − N − 1) respectively and . It is clear, from Eq. (12), that u u u will not be mapped outside I 0 as long as it is initialized inside it. The upper-right block, L I 0 ,R 0 , of the representation L 0 in Eq. (12) is non-zero in general, reflecting that R 0 may not itself be invariant.
To complete this subsection, a procedure for finding an invariant subspace is described. Firstly, we choose a basis for the D 2 − 1 traceless, orthogonal, Hermitian operatorsσ σ σ = {σ j } D 2 −1 j=1 (we have already setσ D 2 = 1). This directly leads to expressions for the matrix representations L and L 0 (see, for example, Ref. [28]). We now examine L 0 and calculate its right-eigenvectors -these eigenvectors will define directions in I 0 ⊂ U 0 . In general, the right-eigenvectors will be complex-valued, and appear as complex conjugate pairs, but we can form a linearly independent real-valued pair of vectors by taking their real and imaginary parts respectively. Each of these real-valued pair of vectors defines a plane, that will form an invariant subspace. Additionally, each real eigenvalue of L 0 will define an invariant subspace via its corresponding real eigenvector. Obviously we can then form larger invariant spaces by combining any of the smaller invariant spaces -this will be necessary, in general, as we require I 0 to be of dimension at least D − 1. The eigenvectors of L 0 will not typically be orthogonal, so that they will not provide an orthogonal basis for U 0 . Consequently, the final step is to find such an orthogonal basis, with each basis vector belonging solely either to I 0 or R 0 . As mentioned, there will typically be eigenvectors of L 0 that are linearly independent from, but not orthogonal to, I 0 -these will not then span R 0 , meaning that R 0 is not itself an invariant subspace. In the final discussed basis, L 0 will take on the form given in Eq. (12).

Utility
Once an appropriate invariant subspace has been identified, it is natural to look for solutions to Eq. (4) that lie entirely in this subspace. That is, where, as a reminder, {|φ k } K k=1 are the PRE members (if the PRE exists). As |φ k φ k | are pure states, they are extremal points of D I . In terms of the space U 0 , Eq. (13) is equivalent to where {u u u k } K k=1 describe the PRE members. Note that the requirement that {u u u} K k=1 correspond to pure states is enforced for For D > 2 further constraints are imposed [27].
Assuming that we have made a change of basis such that L 0 is of the form Eq. (12), the form of Eq. (4) in U 0 space is ∀k, The utility of assuming u u u k ∈ I 0 is now made clear as, since u u u k has only non-zero components in its first N = dim(I 0 ) coordinates, Eq. (15) simplifies to with u u u k | I 0 being the restriction of u u u k to the domain of I 0 . In other words, the constraint equations of Eq. (4) relating to R 0 are trivially satisfied, thus reducing the size of the polynomial system that defines the PRE. Also, the polynomials are more sparse than if the full space were being considered. Both of these considerations reduce the computational cost of finding solutions to the constraints. There are, in general, multiple solutions to Eq. (4) for a given K. Some fraction of these (which may be zero, one or in between) will lie entirely within D I -it is these solutions that we are focusing on here. For N not too large it should be possible to do an exhaustive search to find all of them, if they exist, or prove that there are none. Whether any are found or not says nothing about whether solutions not belonging to D I exist. That is, solutions {|φ k φ k |} having any non-zero operator support outside of D I will not be found with this approach, a point that will be illustrated later in our work.
An important consideration is whether it is to be expected that solutions of Eq. (4) satisfying Eq. (13) can be found -this consideration being distinct from the fact that they will be easier to find if they do exist. That is, is the dimension of D I sufficient to provide enough free variables to satisfy all the constraint equations of Eq. (4)? In Ref. [1], where invariant subspace considerations were not made, it was argued heuristically, via the counting of free parameters and constraints of Eq. (4), that typically ensemble states are required for a PRE to be possible. This was based on D 2 − 1 constraints per ensemble member, K 2 − K transition rates and 2D − 1 unknowns to describe each ensemble state. The restriction of {|φ k φ k |} to some subset of the extremal set of D (H), as appropriate when considering an invariant subspace, can be achieved by placing constraints on the 2D − 1 state variables. Because of the quadratic dependence of the number of constraints upon the size of D I it is possible to reduce their number faster (in the case of a qubit, at an equal rate) than the those of the free parameters when an invariant subspace (in the sense defined above) is considered. In other words, by reducing the scope of the PRE search (to that of an invariant subspace) it becomes apparent that it can be easier to satisfy a parameter and constraint counting heuristic for PRE existence. What may have been an intractable task (due to the computational complexity of finding PREs) is made possible, with the cost being that some PREs (those, if any, that are not contained in D I ) are not discoverable in this way. An example of the simplification provided by invariant subspaces, for the finding of PREs, is provided in the next section.

An important class, and its constraint and parameter count
An important invariant subspace is the conceptually (and computationally) simple one of real-valued density matrices. That is, some predetermined basis exists in which D I = R D×D ∩D (H). This subspace provides an example of the minimum sized invariant subspace that can be formed that is consistent with Eq. (13) without constraining the Hilbert space dimension of span({|φ k }) to be less than D. This is because the extremal subset of D I = R D×D ∩ D (H) still contains D orthogonal states -pure 'redit' states that are each representable as a ray in a D-dimensional real Hilbert space. A minimal invariant subspace (in the sense described above) is desired in order to reduce, to the greatest extent possible, the polynomial system represented by Eq. (4). Of course D I = R D×D ∩ D (H) will only be immediately relevant when L possesses such an invariant subspace, but the concepts -of choosing an as small as possible interesting invariant subspace and then parametrising this space -are more general than this. Given D I = R D×D ∩ D (H), we then search for an example PRE that is real valued in the utilized basis. The imaginary constraints of Eq. (4) are then automatically satisfied.
Let us now give the new inequality (analogous to Eq. (17)) that must be satisfied, for the number of parameters to equal or exceed the number of constraints, when a real-valued invariant subspace is assumed as well as a real-valued PRE. Considering a generic ME 4, there are K 2 − K transition rates as before, but now only D − 1 free parameters and D 2 − D constraints per ensemble state because we are restricting to real state-vectors and matrices. Solving for the integer minimum ensemble size, K, such that the number of free parameters is at least as large as the number of constraints gives Comparing with Eq. (17), we see that the scaling with D 2 is still present, although the coefficient is reduced by a factor of two. Hence, for large D, a search for PREs of roughly half the expected ensemble size indicated by the heuristic of Ref.
[1] may be justified. In Table 1, the minimal PRE size, K, is given, for small D, with and without the use of a redit invariant subspace. Of course, for small D, the difference is less pronounced, but it may still be of great importance given the numerical difficulty of finding PREs. For D = 3, we find, from Eq. (18) (or Table 1), that we can reasonably hope to find PREs with an ensemble size of K = 4, compared with K = 5 from Eq. (17). It also follows that there are only 6 constraint equations (including normalisation constraint of the state) per ensemble member, giving a (square) polynomial system of size 24. This is  Table 1. The heuristically argued minimum number of PRE members, K, required for the number of parameters to equal or exceed the number of constraints, is provided for small dimension, D. The comparison being made is between MEs that have real-valued density matrices as an invariant subspace (redits) and those that do not (qudits). The parameter counting heuristic suggests that a solution of Eq. (4) may be possible for K equalling or exceeding these values. We are considering MEs that are generic, apart from the presence of the invariant subspace symmetry. a vast improvement upon the size 45 system, that presents in the absence of an invariant subspace, and gives hope that future numerical investigations of PREs can extend to D = 3 and beyond. The case of D = 2, describing a 'rebit', actually leads to PREs of the same size as that predicted by Eq. (17). However, even in the case of D = 2, the analysis above leads us to a better understanding of the nature of the PREs; in our upcoming discussion we highlight the symmetry that exists in the D = 2 PREs of Refs. [1,30].
Note the importance of the fact that D I = R D×D ∩ D (H) contains pure states (all pure redit states, as stated). This is what makes it an appropriate invariant subspace in which to search for a PRE. A simple way to engineer such an invariant subspace is to choose a Lindbladian which has a real-valued matrix representation in a basis where it is acting on column-stacked density matrix elements. Interestingly, this will lead to the condition L I 0 ,R 0 = 0, which means that R 0 will also be invariant. However, for the case of D > 2, the latter space does not contain any pure states and is, therefore, not capable of supporting a PRE. This way to engineer D I = R D×D ∩ D (H) is not, however, the only way -we will provide an example in the following subsection for which R 0 is not invariant.

Qubit examples
4.4.1. Resonance fluorescence We explore further the nature of PREs possessing an invariant subspace symmetry -specifically, that of a real subspace D I = R D×D ∩D (H), discussed in Sec. 4.3 -by re-examining the D = 2 (qubit) resonance fluorescence system that was analyzed, in detail, in Refs. [1,7]. The ME is of the form Eq. (1), witĥ H = Ω 2σ x and a single Lindblad jump operatorĉ = √ γ|0 1|. However, to make explicit the D I = R D×D ∩ D (H) symmetry, we perform a change of basis relative to [1,7], such that |0 → i |0 . The Hamiltonian then becomesĤ = Ω 2σ x = iΩ (|0 1| − |1 0|) /2 and it is clear that a density matrix that is initially real-valued will stay real-valued (note that the rebit plane here is the y-z plane of the Bloch ball -σ x has been chosen as the imaginary Pauli operator).
Working in the Pauli basis, {σ x , σ y , σ z } (with the just described, non-standard, representation), we find that the expressions for L 0 and b b b (appearing in Eq. (7)) are from which we can determine the steady state x x x ss = (0, 2γΩ, −γ 2 ) T /(γ 2 + 2Ω 2 ). This defines ρ ss via Eq. (8) (with D = 2). The three right-eigenvectors of L 0 are: e e e 1 = (1, 0, 0) T and e e e ± = (0, γ ± γ 2 − 16Ω 2 , 4Ω) T (the latter being unnormalized). When all eigenvectors are real (γ 2 −16Ω 2 ≥ 0), each eigenvector gives a one-dimensional invariant subspace. When γ 2 − 16Ω 2 < 0, the real and imaginary components of e e e + (equivalently e e e − ) are used to form a two-dimensional subspace, and the one-dimensional subspace corresponding to e e e 1 remains also. The structure of the eigenspaces is such that the space formed by e e e 1 is always orthogonal to the other spaces (be they one-or twodimensional). However, the other one-dimensional spaces (when they exist) are not respectively orthogonal, but together span the space orthogonal to e e e 1 . The connection to the space U 0 (being a displaced 3-ball parameterized by u, v, w) and I 0 is as follows. When γ 2 − 16Ω 2 ≥ 0, there are three one-dimensional I 0 spaces that are defined by the following directions: the u-axis (due to e e e 1 ) and two other rays lying in the u = 0 great disc (due to e e e ± ). When γ 2 − 16Ω 2 < 0, I 0 is either (along) the u-axis or the full u = 0 great disc. If I 0 is taken to be the u-axis, there is a duality, in that R 0 itself is invariant (being the u = 0 great disc). This could be inferred from Eq. (19) as it has the form of Eq. (12), but with L I 0 ,R 0 = 0. This duality, L I 0 ,R 0 = 0, is not present in the case when both γ 2 − 16Ω 2 ≥ 0 and the invariant subspace I 0 is chosen as one of the two rays lying in the u = 0 great disc. The non-orthogonality of e e e ± ensures that a state initialized in R 0 will not, in general, remain in R 0 The reader should remember that U 0 (in which the vectors u u u live) is not the Bloch ball (centered at the origin), but rather a unit ball with origin x x x ss .
Recalling our requirement that I 0 be of dimension at least D − 1 = 1 and that it contain pure states (when mapped to density matrices via Eq. (10)), we see that all of the eigenspaces discussed above are interesting, and we can search for PREs contained fully in each of them, respectively. Note that all the extreme states of U 0 correspond to pure states, a feature unique to D = 2.
First, we consider the one-dimensional I 0 subspaces in the context of PREs. When γ 2 −16Ω 2 ≥ 0 each of the three I 0 meets the surface of U 0 in two locations, corresponding to two pure states on the Bloch sphere. Remembering that the PRE must reside in I 0 (as per Eq. (14)), we see that it is only feasible for PREs with K = 2, and each extremal point in I 0 must be an ensemble member.
The case of D = K = 2 was treated analytically in Ref. [1], with expressions describing the PRE resulting. Translating their results to the space U 0 -which we have defined in our work -the PRE states are  ) with Ω/γ = 0.18. The volume of the ball at the tip of each PRE state arrow represents the probability ℘ k that the qubit occupies the corresponding pure state. The steady state, x x x ss , is plotted as a black arrow but it is not often visible as it lies very close to the overlaid PRE states -for example, in subplot (a) it lies almost coincident with the highly occupied PRE member. In subplot (d), the PRE of subplot (a) is reproduced, but in the u = 0 subspace of U 0 , within which it is contained. In the space U 0 , the steady state lies at the origin, marked by a green +. The brown discs show the PRE state locations, with their area proportional to their occupation probability. (The green + is almost concentric with the larger brown disc.) The circle shows the boundary of U 0 , upon which lie pure states. The steady state is only marginally mixed and consequently lies very close to the circle. The secant marks I 0 as a line through the origin meeting the surface at two pure states.
where k = {1, 2} labels the ensemble member, j labels the real eigenvectors, e e e, of L 0 and c is a scalar that depends on x x x ss , e e e j , and k. The important point for our purposes is that the PRE states in U 0 are in the direction of e e e j -in other words, all K = D = 2 PREs lie in I 0 and each one-dimensional I 0 supports a PRE (via its corresponding real eigenvector). We display each of the K = D = 2 PREs in Fig. 1, with subplots (a),(b) associated with e e e ± and subplot (c) with e e e 1 5. In addition to plotting the PREs on the Bloch ball, subplot (d) shows the PRE of subplot (a) in the u = 0 subspace of U 0 . When γ 2 − 16Ω 2 < 0, the one-dimensional spaces corresponding to eigenvectors e e e ± are no longer available and there exists only a single K = 2 PRE which is attributable to the e e e 1 eigenspace. For K = 3, the D I arising from the three one-dimensional I 0 (that pierce the U 0 ball at only two points) cannot hold PREs, so it is natural to look at the two-dimensional I 0 subspace (being the u = 0 great disc). A complete discussion of K = 3 is given in Ref. [7] for the interested reader -in this paper we just wish to point out the symmetry that exists in some but not all of the K = 3 PREs. For example, Fig. 2(a),(c),(d),(f) 6 all contain PREs with u = 0 (corresponding to x = 0 in the Bloch ball figures). That is, they reside in the D I pertaining to the two-dimensional I 0 discussed above. Importantly, there are PREs that do not possess the invariant subspace symmetry (see Fig. 2 subplots (b) and (e)). This illustrates the fact that by searching for PREs in a reduced space [D I compared with the total space D (H)], one may not find all PREs for a given K. This will also be a feature when we consider a different type of symmetry (Wigner symmetries) in the next section. It is important to note that, following Ref. [1], we have limited to cyclic PREs when illustrating our logical points. If non-cyclic PREs exist (necessarily having K > 2), some of them may well possess the invariant subspace symmetry; a search for them would be significantly simplified utilizing the methods of this section.

Absorption and emission ME
Another qubit ME that possesses much symmetry is that which models incoherent emission and absorption. That is, we setĤ = 0, but add a second decoherence channel, as compared with Sec. 4.4.1, so thatĉ 1 = √ γ − |0 1| and c 2 = √ γ + |1 0|. This ME was briefly considered, in the context of PREs, in Ref. [30]; the reader is also referred to the application discussed in Ref. [32] as evidence of its relevance to quantum information. In the Pauli basis we obtain where γ Σ = γ − + γ + and γ ∆ = γ + − γ − . Without loss of generality we assume γ ∆ ≤ 0 7. The steady state is x x x ss = (0, 0, γ ∆ /γ Σ ) T , which lies on the z-axis of the Bloch ball. The three right-eigenvectors of L 0 (all of which are real) are e e e 1,2,3 = {(1, 0, 0) T , (0, 1, 0) T , (0, 0, 1) T } having eigenvalues −γ Σ {1/2, 1/2, 1} respectively. An important difference from the resonance fluorescence case is that there is a degeneracy of eigenvalues, the consequence of which is that every diameter of the w = 0 disc of the 3-ball U 0 is an invariant one-dimensional subspace, I 0 . Of course, the portion of the w-axis within U 0 is also an invariant one-dimensional subspace. The case of K = 2 is discussed in Ref. [30], and it is indeed the case that each of the invariant subspaces contains a K = 2 PRE. An infinite number of PREs are parametrized by the azimuthal angle in the w = 0 disc, and there is one further PRE consisting of the two points where the w-axis and U 0 intersect. In summary, each I 0 space contains a K = 2 PRE and there are no K = 2 PREs that do not possess the invariant subspace symmetry -this is expected as per the discussion below Eq. (20). K = 3 PREs were not investigated in Ref. [30], but it is of interest for us to search for them, making use of the discussed invariant subspace PRE symmetry technique. For simplicity, we once again limit the investigation to cyclic PREs. Of course the onedimensional I 0 are not sufficient for K = 3 PREs as they only correspond to two pure states. Instead, we form two different two-dimensional I 0 within the 3-ball U 0 : firstly, the w = 0 disc and, secondly, the v = 0 great disc. Within the first subspace there are no cyclic K = 3 PREs. This is an analytic result, as the constraints imposed by Eq. (4) can be shown to be algebraically inconsistent for arbitrary, but non-zero, γ ± 8. Within the v = 0 great disc, our numeric results depend on the parameter regime. For γ + > εγ − , with ε ≈ 1/18, no PREs are found (recall that we restrict to γ − ≥ γ + ). But for γ + < εγ − , two cyclic K = 3 PREs exist. Examples of these are shown in Fig. 3.
The results pertaining to the v = 0 great disc are numeric in that we specify the ratio γ + /γ − before applying Gröbner basis techniques [33] either of a algebraic (using the computational algebra software MAGMA [23]) or numeric nature (MATHEMATICA) to solve the polynomial system. That is, we can be sure of the existence of PREs only 7 Because the ME is invariant under |0 ↔ |1 and γ + ↔ γ − , the results for γ ∆ ≥ 0 can be simply obtained from those for γ ∆ ≤ 0. If a ME, parameterized by {γ + , γ − }, possesses a PRE with K members {x k , y k , z k } K k=1 , then the ME parameterized by {γ − , γ + } will possess a PRE of the form {x k , −y k , −z k } K k=1 . 8 A straightforward way to verify this is to formulate the constraint equations in terms of the Bloch vectors (rather than the state vectors of Eq. (4)) then, due to the symmetry of the system, choose say v 1 = 0. By inspection, this will lead to v 2 = v 3 = 0 also. This is a contradiction as it excludes the possibility of three distinct states in the ensemble. Resultantly, we conclude that there are no such PREs contained in the w = 0 disc.  21) -each is comprised of dots of a single colour, either red or black. As they lie on the v = 0 great circle, only that great disc of U 0 is displayed. The relative area of the dots indicates the relative occupation probability amongst the ensemble members of each ensemble respectively. The direction of cycling within each PRE is shown by the arrows. The steady state location is shown by the green cross, which, by definition, lies at the origin in the u, v, w coordinate system. The parameter choice is γ + = 0.05γ − . at discrete γ + /γ − ratios. However, by conducting a fine-grained search (we stepped in 0.001 increments) we can develop a clear picture 9.
In fact, because of the ME having no preference in the u-v plane (x-y plane in the Bloch ball), any plane containing all of the w-axis is actually a two-dimensional I 0 and possesses PREs corresponding to azimuthal rotations of the v = 0 disc PREs 10. That is, we can obtain an infinite number of PREs (forming a 'family') by rotation of each member of the PRE about the w-axis -this is in complete analogy with the family of K = 2 PREs that are related by rotation. In this section, these families appear somewhat ad hoc; in the next section, the theory of such PRE families will be explicated, with their origin attributed to a symmetry different from the invariant subspace symmetry. Following a general theoretical development, we will return to our qubit examples.
For completeness, a search for K = 3 cyclic PREs for the ME of Eq. (21) can be considered. However, by inspection of the constraint equations, it quickly becomes analytically apparent that no further PREs beyond those contained within I 0 are possible. The demonstration of this is as per 8, where, once again, the symmetry of the system allows a simplification that quickly leads to the conclusion that a PRE must lie on a great disc containing the w-axis.

Wigner invariance of L
Wigner transformations act on Hilbert space rays in a way that preserves the Hilbert space inner product, | ψ 1 |ψ 2 |. Their action is consequently well defined upon pure state projectors, and we denote this as T |ψ ψ|. Wigner showed that such transformations are either unitary (and so linear in their action on Hilbert space) or antiunitary (and so antilinear) [34,35]. In this section, we consider those Wigner transformations which leave the Lindbladian, L, invariant: and term these 'Wigner symmetries'. Using the unitary/antiunitary properties of T and the block structure of L implied by Eq. (7) we can gain further insight into Eq. (22), which reduces to where T 0 , analogously to L 0 , is the restriction of the matrix representation, T, of the superoperator T to traceless Hermitian matrices. Given that the ME has, by assumption, a unique steady-state, we can use Eqs. (23)- (24) and Eq. (9) to infer that the steady is invariant under the action of T , Together, Eq. (23) and Eq. (25) provide a useful formulation of Eq. (22).

Wigner-symmetric families of PREs
By using the Wigner symmetry of Eq. (22) in Eq. (4) we can investigate its significance for PREs. After utilizing the symmetry and then premultiplying by T we obtain Note that the κ jk , because they are real valued, commute with all T , including antiunitary transformations. It is clear that a solution to Eq. (4) can be used to generate new PREs under the action of T . Given how computationally difficult it is to find a PRE, this is very valuable. We also note that the set of all T satisfying Eq. (22) forms a group, G. This is apparent since: the identity I satisfies Eq. (22), every Wigner transformation T has an inverse T −1 that is also a Wigner transformation, and if T 1 and T 2 satisfy Eq. (22) then so does T = T 2 T 1 . If G (or a subgroup of G) is a Lie group then there is an interesting consequence, namely that an infinite number of PREs can be found irrespective of whether or not the constraining polynomial system (defined by Eq. (4)) is square. (Note that a Lie group could not include antiunitary symmetries, which are necessarily discrete.)

Wigner-symmetric PREs
In addition to being used to generate families of PREs, as discussed above, the Wigner symmetry described by a group G can also be applied to make it easier to find individual PREs. To do so we first introduce the following notation. Let P be a permutation of {1...K} with k = P (k). Then we define a PRE as having the Wigner invariance for some T ∈ G iff ∃ permutation P such that If this holds, then the action of T (or T −1 ) on any ensemble member generates an existing member, which could be itself (k = k) or a different member (k = k). It is possible that Eqs.
implies that the k th is satisfied also: If we define an equivalence relation, ∼, amongst ensemble members in the presence of the symmetry T as then the constraints on only one element of each equivalence class, [|φ k ] ∼ , need to be tested as the remainder are implied. A PRE being Wigner-symmetric is consistent with the invariance with which we began in Eq. (22) 11. The equivalence class [|φ k ] ∼ will have more than two members if the different elements in K require different permutations P in order to satisfy Eqs. (27)- (28). As an example, this situation will arise if the period of the elements of K is greater than two (that is, T 2 = 1).

Qubit examples
5.3.1. Resonance fluorescence We can, once again, use the resonance fluorescence qubit ME described by Eq. (19), this time to illustrate the Wigner symmetry. It is invariant (in the sense of Eqs. (23)- (24)) under the antiunitary transformation T 0 = diag (−1, 1, 1), defined here using the Pauli operator basis {σ x , σ y , σ z }, where we have restricted to the traceless Hermitian operators. The steady-state is also invariant under the symmetry's action, as required by Eq. (25). Clearly this Wigner symmetry, which takes σ x → −σ x , is a Z 2 symmetry, meaning T 2 = 1. All three of the K = 2 PREs obey the Wigner symmetry, T |φ k φ k | = |φ k φ k |. For two of them -see Fig. 1(a),(b) -the permutation is trivial and each ensemble member is mapped to itself (k = k ). These PREs necessarily lie wholly on the x = 0 great disc which is the L-invariant subspace attributable to e e e ± . The remaining K = 2 PRE, which lies in the L-invariant subspace attributable to e e e 1 -see Fig. 1(c) -has a non-trivial permutation associated with T ; the ensemble members are mapped to each other.
Regarding the K = 3 cyclic PREs, four of the eight possess the Wigner symmetry (see Fig. 2(a),(c)(d),(f)), but only in the trivial manner just described, lying on the x = 0 great disc. Those PREs that are not Wigner symmetric come in pairs, such that one PRE in the pair is obtained from the other under the action of T , in the manner discussed in Sec. 5.1, as seen in Fig. 2(b),(e). The fact that some of the PREs in this example do not possess the Wigner symmetry highlights that, similarly to the case of invariant subspaces of L, discussed in the previous subsection, by imposing PRE Wigner symmetry we may not find the entire solution set.
It is not a coincidence that the Wigner-symmetric K = 3 cyclic PREs of the resonance fluorescence ME, that are found in Ref. [7], have the symmetry present in only a trivial form, with each ensemble member being mapped to itself. The specific Wigner-symmetry in question is T = diag(−1, 1, 1, 1), so that a K = 3 PRE possessing the symmetry, in a non-trivial manner, would have to have two ensemble members of the form |φ 1 φ 1 | ≡ (x, y, z) and |φ 2 φ 2 | ≡ (−x, y, z) for x = 0, and a third ensemble member as |φ 3 φ 3 | ≡ (0, y 3 , z 3 ). Additionally, according to Eq. (28), we must have κ 1,3 = κ 2,3 and κ 3,1 = κ 3,2 . This is inconsistent with the assumed cyclic nature of the PRE and, as a consequence, rules out their existence.
It is feasible to look for non-cyclic K = 3 PREs possessing the Wigner symmetrythe PRE structure just described collapses the polynomial system greatly. However, a difference to the previously discussed numerical searches is that the polynomial system is no longer square. The additional potential transitions beyond cyclicity provide extra variables, which then exceed, in number, the constraint equations (in this case, by one). Generically, one can then expect a positive dimension solution set (here of dimension one), before imposing that we require real-valued solutions and positive transition rates. To avoid this difficulty, we run repeated numerical tests at discrete values of a chosen variable (100 evenly spaced values of x 1 ), thus reducing the number of parameters to be equal to the number of constraints. This collapses the problem to a set of square polynomial systems. Furthermore, we also carry out this procedure for a series of 200 different MEs, parameterized by the value of Ω/γ (varied between 0 and 10). Despite the fact that the obstacle of the previous paragraph is avoided (by allowing non-cyclicity), no non-trivial Wigner symmetric PREs were found. We can be reasonably confident that none exist, but not certain, due to the numerical nature of the search. The reader should not think that the apparent non-existence of non-trivial PREs with K > 2 in this example indicates that none are possible for a qubit. In the ME example of Eq. (21), Wigner symmetric PREs for arbitrarily large K exist, and will be presented in Sec. 6.1.2.

Absorption and emission ME
We now return to the ME of Eq. (21) that describes incoherent absorption and emission, and consider Wigner symmetries. It is not hard to identify that it has the orthogonal Lie group, O(2) (consisting of rotations and reflections), acting on the first two coordinates, as a Wigner symmetry. That is, For clarity, R(2) is a 2 × 2 matrix representation of the group O(2). Note that although an inversion of the z-coordinate would also satisfy Eq. (23), it does not satisfy Eq. (24) (unless γ + = γ − ); thus, the lower right element of T 0 must be unity. It is therefore expected that families of PREs related by the action of T 0 will be obtained. On the Bloch sphere, this means that for a given PRE, we expect a second PRE to exist that is related to the first by rotation about the z-axis and/or the reflection about a plane containing the entire z-axis. As O(2) is a Lie group, an infinite number of PREs (some of which may be indistinct) will belong to each family. For K = 2, the PREs belong to two families: those contained in the z = z ss plane of the Bloch ball and the PRE comprised of the poles of the Bloch ball. The latter PRE has only one PRE in its family as it is mapped to itself. As expected, the actual K = 2 PREs conform to these predictions. A description of cyclic K = 3 PREs was given in Sec. 4.4.2, and we can also now understand the origin of those families of PREs (that exist for γ + γ − /18): they are being generated under the action of T 0 , just as for the K = 2 PREs.
We now consider Wigner-symmetric PREs; that is, each of the members of a given ensemble must be mapped to one another under the action of T , with the transition rates also related by Eq. (28). The PRE contained on the z-axis is a Wigner symmetric PRE, but only under the trivial permutation in which each member is mapped to itself. As the transition rates between ensemble members are different when γ + = γ − (when γ + = γ − the ME has the additional Wigner symmetry under reflection in the z = 0 plane) it is not possible for this PRE to have a Z 2 symmetry. For K = 2, we require that T 2 0 = 1, which implies that R(2) is either a rotation by π or a reflection. Note that K = 2 ensembles whose states are mapped to each other under reflection, but are not antipodal, cannot form PREs as their ensemble average will not lie on the z-axis. Wigner symmetric K = 2 PREs are, therefore, those that consist of the antipodal points of the intersection of the z = z ss plane and the Bloch ball and, additionally, have each member occupied with equal probability. Thus, the family of K = 2 PREs lying in the z = z ss Bloch plane are Wigner symmetric.
As for the cyclic K = 3 PREs, these are not Wigner symmetric -they possess neither a T 2 0 = 1 (with the third ensemble member mapped to itself) nor a T 3 0 = 1 symmetry (which can be seen as this would require them to be mapped outside of the plane containing the z-axis under the action of T 0 ). The theory to investigate some larger (K > 3) and more complicated (non-cyclic) PREs, for this ME, is now in place, but we postpone this to the following section, as the combination of symmetries makes finding some of them particularly easy.

Combining symmetries
The reader will have noticed that many of the example PREs thus far presented have possessed both the Wigner symmetry and the invariant subspace symmetry. In this section we further explore the simultaneous existence of these symmetries at a ME and PRE level.
The requirement that a ME possess the invariant subspace symmetry of Sec. 4 in addition to the Wigner symmetry of Eq. (22) can be formulated as requiring Eq. (23) to be satisfied when L 0 has the block form given in Eq. (12). This needs to be supplemented by Eq. (25) as it is non-trivial that the steady-state will possess the Wigner symmetry. Given the satisfaction of these constraints, then the results of Sec. 4 and Sec. 5 can clearly be applied to search for PREs, for the ME in question, that possess either one of the two symmetries.
The obvious next question is, under what conditions can we look for PREs that possess both symmetries jointly? This would allow us to leverage both symmetries and greatly simplify the finding of (jointly symmetric) PREs. To investigate, we give T 0 the block structure induced by that of L 0 (see below Eq. (12) for further details) Firstly, from Eq. (27), in the basis for which L 0 has the block form of Eq. (12), it is immediate that T R 0 ,I 0 = 0 is necessary, as otherwise T would map PRE members to outside of the invariant subspace. This then implies that T I 0 ,R 0 = 0 as well, due to the unitary/antiunitary nature of T . Secondly, as the PRE is proposed to lie entirely in I 0 , the effect of T on coordinates outside I 0 is irrelevant. We collect these arguments in the equations T 0 x x x ss = x x x ss .
The orthonormal basis for U 0 also defines the basis in which we write x x x ss . Note that Eq. (36) cannot, in general, be written solely in terms of T I 0 because x x x ss may have support on (translated) coordinates outside of I 0 . An example of this is the resonance fluorescence qubit ME, for which the u-axis was an I 0 space but x x x ss has support outside of the x-axis (which is the translated image of the u-axis). Another way of saying this is that I 0 is determined by L 0 but x x x ss is also dependent upon b b b, as per Eq. (9). It is worth pointing out that, in principle, Eqs. (34)-(36) can be satisfied even when both symmetries are not present at the ME level. That is, when searching for PREs that simultaneously have invariant subspace symmetry and Wigner symmetry it is not necessary that the Wigner symmetry apply to the entire D (H); all that is relevant is its action upon D I (in contrast to both symmetries existing at the ME level). In other words, it may be the case that

Qubit examples
6.1.1. Resonance fluorescence As a first example of the combination of symmetries, let us consider the example of Ref. [1]. As a reminder, we have the Wigner symmetry T 0 = diag(−1, 1, 1) and invariant subspaces, I 0 , being, respectively, the u-axis and the u = 0 great disc. It is immediate that T R 0 ,I 0 = 0, so that Eq. (34) holds and it is clear that the steady state is invariant under the action of T . The first invariant subspace has L I 0 = −γ/2 and T I 0 = −1, both scalars, so that Eq. (35) holds. The second invariant subspace, the great disc, has u = 0 so that T I 0 acts as the identity. We also have T −1 LT = L, so that the Wigner symmetry applies across the entire D (H), not just D I . Thus, the Wigner and invariant subspace symmetries can co-exist both on a ME level and at that of the PREs. Indeed a non-trivial manifestation of the Wigner symmetry is found in Fig. 1(c), which shows a PRE constrained to the Bloch ball image of the first mentioned invariant subspace.
6.1.2. Absorption and emission ME Returning to the ME of Eq. (21), we will now exploit the abundance of symmetry that it possesses to find PREs of arbitrary size K. We avail ourselves of the maximum amount of PRE Wigner symmetry; we take it to be of the form Z K , so that This means that of the K matrix equations in Eq. (4) only one of them need be considered -the rest are guaranteed to be satisfied due to the Wigner symmetry. The explicit form of T is an azimuthal rotation of magnitude θ = 2π/K. Consequently, each ensemble member must have the same w coordinate -the entire ensemble lies in the w = 0 invariant subspace of U 0 . Furthermore, if a single PRE exists, of the nature described, then, due to the Wigner family of PRE symmetry, there must exist an infinity of PREs such that every point of the w = 0 circle of U 0 hosts an ensemble member that belongs to atleast one PRE. This allows us to completely fix all coordinates of the state |φ k φ k | and it can be written as, for example, { 1 − γ 2 ∆ /γ 2 Σ , 0, 0} T (any point on the w = 0 circle could have been chosen). Once we have found a single PRE then, of course, the entire family is obtained by rotation. In summary, the potential PRE that we are As the PRE is entirely contained in the w = 0 disc, we show this two-dimensional space. The PRE members are shown as black dots. All transition rates are of the same magnitude (γ Σ /6) and each member has equal occupation probability. In (b) we show a K = 4 PRE, also in the w = 0 plane. All the blue transitions are of equal magnitude and are non-zero. All red transitions are also of equal magnitude but there exists the possibility of them all being zero. An infinity of PREs can be obtained by arbitrary azimuthal rotation.
investigating consists of K evenly spaced states lying on the w = 0 circle and, without loss of generality, we take one of them as lying on the u-axis.
Whether this PRE exists or not is determined by Eq. (4). We insert the PRE states into these equations and the previously difficult to solve system of polynomial equations is reduced to only two linear equations that constrain the transition rates away from the kth state. Taking the one state that determines them all to be labelled as k = 1, we have where κ j,1 denotes the transition rate from state 1 to state j. The above equations can be satisfied for arbitrary K ≥ 2. By setting K = 2 one regains the PREs described in Sec. 4.4.2. It is also clear (in a mathematical sense) that no cyclic K = 3 PREs can possess the Wigner symmetry. In fact, from Eq. (39) and Eq. (28), no cyclic PREs possess this Wigner symmetry for any K > 2. To find actual existing PREs from Eqs. (38)-(39), non-cyclic PREs need to be considered. This is because, generally speaking, the two constraints require two free parameters to have a solution. In this way we can find PREs having a minimum of two outward transitions per member. Specifically, for K = 3, each of the two transition rates is given by γ Σ /6. A PRE from this family is illustrated in Fig. 4(a). If more than two transitions are allowed (possible when K > 3), then the excess parameters, as compared to constraints, gives hyperplanes of solutions in κ space. Specifically, for K = 4, we must have equal transition rates to the two ensemble members that are obtained by a magnitude π/2 azimuthal rotation (κ 2,1 = κ 4,1 ). Additionally, the following constraint must be satisfied κ 2,1 + κ 3,1 = γ Σ /4 (a line of solutions in {κ 2,1 , κ 3,1 } space -note that κ 2,1 = 0 would return us to the K = 2 PRE as the K = 4 PRE would get trapped in a K = 2 subcycle). Additionally, it is of interest that we do have freedom to set some of the transitions to zero for K = 4: κ 3,1 = 0, together with κ 2,1 = κ 4,1 = γ Σ /4 satisfies the constraints. We show a sample K = 4 PRE in Fig. 4(b). In summary, by using our theory of symmetry to the maximum extent, we have found non-cyclic PREs for arbitrarily large K -a task that would naively be undertaken by trying to solve a, typically intractable, set of polynomial equations.

Discussion and conclusion
In this paper we first discussed the difficulty of finding PREs, then provided some new tools to reduce the complexity of this task. That is, we have considered symmetries of MEs that allow a simplified search for PREs that comprise a subset (that may be empty, but in many cases we have shown otherwise) of those that exist.
The invariant subspace symmetry involved finding a subspace (of the space occupied by density matrices) to which dynamics are constrained, given an appropriate initialisation. It was then logical to propose searching within this confined subspace for PREs, a task that can be considerably more simple than searching the entire space. This was highlighted by considering the class of real-valued density matrices and noting that the expected minimum PRE size is halved (for large D) when this symmetry is appropriate. Due to the exponential scaling, in the ensemble size, of the difficulty of finding PREs, this can make previously intractable searches tractable. Qubit examples of this symmetry were then analyzed in detail, giving new insights into previously known PREs. For example, it was shown that many of the known K = 3 resonance fluorescence ME PREs actually possess the symmetry, as do all of the K = 2 PREs.
Next we considered MEs possessing a Wigner symmetry (where the Hilbert space inner product is preserved). This led to the exciting conclusion that new PREs could be generated from existing ones by using the symmetry. Also, perhaps most importantly, it was shown that if the potential PRE is assumed to possess the Wigner symmetry then a fraction of the constraints governing their existence become redundant. This serves to reduce in size the polynomial system that must be solved to find a PRE. This all comes with the proviso that not all PREs can be found in this manner -there may exist unsymmetric PREs that escape the scope of the symmetry-based searches. In fact, for the qubit examples, we find several of these PREs; that is, there exist PREs in both the symmetric and unsymmetric categories.
Finally, we investigated a joint application of the invariant subspace and Wigner symmetries. For the absorption and emission qubit ME, this allowed us to reduce the constraints defining these symmetric PREs to two linear equations, essentially making trivial the finding of arbitrarily large PREs.
Finding the adaptive measurement scheme necessary to produce a PRE is typically a less demanding computational task than finding the PRE itself, as the system of coupled polynomials that must be solved is smaller. For this reason, we have focused on PRE symmetries rather than measurement scheme symmetries. However, a brief discussion of the latter is contained in Appendix A and Appendix B.
Having demonstrated the practicable and insightful use of the developed symmetry tools, we may now consider new problems to which they can be applied. In Refs. [1,7] a number of important open questions regarding PREs are raised. The first question is: are there MEs for which the minimally sized PRE is larger than D? The motivation of this query is that if it is answered in the affirmative then open quantum systems can be said to be harder to track than classical systems, as there exist quantum systems that cannot be tracked by a D-state finite resettable classical memory (in contrast to classical systems). Clearly our techniques that find only a subset of potential PREs cannot be used to rule out their existence, but we can, of course, place an upper bound on the ensemble size by finding a PRE using symmetry. The second question is as follows: is an ensemble size of K = (D − 1) 2 + 1 always sufficient for a PRE to be found? Once again, our techniques do not allow a direct answer, but would provide an upper bound on the difficulty. The third question is: does the ensemble size of K = (D −1) 2 +1 from Ref. [1] reliably predict whether PREs are feasible for a ME of a given form? To this question, our techniques and results are directly relevant as we showed that K = (D − 1) 2 + 1 is potentially larger than necessary if symmetry is present. Specifically, we showed that ensembles of only K = (D 2 − D + 2)/2 pure states can be expected to be sufficient in the case of real-valued density matrices. But rather than seeing this as an immediate negative answer to the third question, we see it as pointing the way to refinement of the heuristic governing PRE existence in the presence of invariant subspace symmetry. Since we have reduced the PRE size of interest (in the symmetric case) we should be able to test the refined heuristic more easily. We expect to provide further refining of the expected PRE size based on other considerations (such as the number of decoherence channels, L) in future work.
As a final topic for future work, we suggest that the study of a D = 4 composite system (two qubits) in terms of the statistics and dynamics of entanglement of discovered PREs would be intriguing. The construction of multi-qubit systems in a symmetric manner may make possible such an investigation. member. Similarly, the no-jump evolution also maps impure states belonging to D I outside of D I . The no-jump evolution increases the excited state occupation probability of such states, as can be inferred from the average evolution preserving D I (with the latter being true by our definition of an invariant subspace).
In contrast, to provide an example of a measurement scheme that does possess the invariant subspace symmetry, consider the K = 3 PREs, arising from e e e ± , that are contained in the two-dimensional invariant subspace consisting of the rebit great disc. The LO amplitude that achieves these K = 3 PREs is purely real, with the implication that both the jump and no-jump evolution map real-valued density matrices to realvalued density matrices, thus preserving D I . Consequently, the measurement scheme possesses the invariant subspace symmetry as we have defined it.
describing the Hamiltonian -Ω in the case of resonance fluorescence -will most commonly be independent from that describing the decoherence, γ for example), it is simple to show that the no-jump constraintsĤ l eff |φ l ∝ |φ l andĤ l eff |φ l ∝ |φ l are automatically satisfied. The constraint, Eq. (B.3), also relates measurement schemes of different PREs in the case that the PREs belong to the same T Wigner-symmetric family (see Sec. 5.1).