Evolutionary games of condensates in coupled birth–death processes

Condensation phenomena arise through a collective behaviour of particles. They are observed in both classical and quantum systems, ranging from the formation of traffic jams in mass transport models to the macroscopic occupation of the energetic ground state in ultra-cold bosonic gases (Bose–Einstein condensation). Recently, it has been shown that a driven and dissipative system of bosons may form multiple condensates. Which states become the condensates has, however, remained elusive thus far. The dynamics of this condensation are described by coupled birth–death processes, which also occur in evolutionary game theory. Here we apply concepts from evolutionary game theory to explain the formation of multiple condensates in such driven-dissipative bosonic systems. We show that the vanishing of relative entropy production determines their selection. The condensation proceeds exponentially fast, but the system never comes to rest. Instead, the occupation numbers of condensates may oscillate, as we demonstrate for a rock–paper–scissors game of condensates.

C ondensation phenomena occur in a broad range of contexts in both classical and quantum systems. Networks such as the World Wide Web or the citation network perpetually grow by the addition of nodes or links and they evolve by rewiring. Over time, a finite fraction of the links of a network may be attached to particular nodes. These nodes become hubs and thereby dominate the dynamics of the whole network; they become condensate nodes [1][2][3] . Condensation also occurs in models for the jamming of traffic [4][5][6][7] and in related mass transport models in which particles hop between sites on a lattice 3,8,9 . A condensate forms when a finite fraction of all particles aggregates into a cluster that dominates the total particle flow. Bose-Einstein condensation, on the other hand, is a quintessentially quantum mechanical phenomenon. When an equilibrated, dilute gas of bosonic particles is cooled to a temperature near absolute zero, a finite fraction of bosons may condense into the energetic ground state [10][11][12] . Long-range phase coherence builds up and quantum physics becomes manifest on the macroscopic scale 13,14 .
In both the classical and the quantum mechanical context, condensation occurs when one or multiple states become macroscopically occupied (they become condensates), whereas the other states become depleted 15,16 . However, the physical origins of condensation in the above examples differ from each other. Why and how condensation arises in a particular system remains a topic of general interest and vivid research.
Here we study condensation in two systems from different fields of research: incoherently driven-dissipative systems of noninteracting bosons and evolutionary games of competing agents. As we show below, the physical principle of vanishing entropy production governs the formation of condensates in both of these systems. The entities that constitute the respective system shall be called particles. They may be quantum or classical particles (bosons or agents). The dynamics of these particles eventually lead to the condensation into particular states (quantum states or strategies). Before describing the above two systems, we now introduce the mathematical framework of our study.
On an abstract level, we consider a system of S (nondegenerate) states E i ; i ¼ 1; . . . ; S, each of which is occupied by N i Z0 indistinguishable particles, see Fig. 1a. The configuration of the system at time t is fully characterized by the occupation numbers N ¼ ðN 1 ; N 2 ; . . . ; N S Þ. This configuration changes continuously in time due to the transition of particles between states. The total number of particles in this coupled birth-death process is conserved ðN ¼ P i N i Þ. We are interested in the probability P(N, t) of finding the system in configuration N at time t. The temporal evolution of the probability distribution P(N, t) is governed by the classical master equation 17,18 : where e i 2 Z S denotes the unit vector in direction i (equal to 1 at index i, otherwise 0). The rate for the transition of particles from state E j to E i depends linearly on the number of particles in the departure and in the arrival state: with rate constant r ij Z0 and constant s ij Z0. Condensation in this framework is understood as the macroscopic occupation of one or multiple states 15,16 : We consider a state E i as a condensate when the long-time average of the number of particles in this state scales linearly with the system size (hN i i t $ OðNÞ for large t). Hence, a condensate harbours a finite fraction of the total number of particles for large systems (N ) 1). We refer to a state as depleted when its average occupation number scales less than linearly with the system size. Therefore, the fraction of particles in a depleted state vanishes in the limit of large systems.
Depending on the values of the rate constants r ij , numerical simulation of equations (1) with rates (2) reveals that all states, multiple states or only one state become condensates when the particle density N/S is large enough to detect condensation 19 . Thus far, various questions about condensation have remained elusive for the coupled birth-death process defined by equation (1): Which of the states become condensates? How does this selection of condensates proceed? Is it possible to construct systems that condense into a specific set of condensates?
In the following, we answer these questions by illuminating the physical principle that governs the formation of multiple condensates on the leading-order timescale. We show that the vanishing of relative entropy production determines the selection of condensates (see equations (3) and (4) below). We elaborate how condensate selection is determined by the rate constants r ij . The condensation proceeds exponentially fast into a dynamic, metastable steady state within which the occupation numbers of condensates may oscillate. By applying our general results to systems with many states, we show that the interplay between critical properties of such networks of states 20 and dynamically stable network motifs 21 determines the selection of condensates. The results of our analysis apply to any system whose dynamics  Figure 1 | Condensation into multiple states due to particle transitions between states and mathematics of condensate selection. (a) With respect to condensation in an incoherently driven-dissipative quantum system, each bowl represents a state E i that is occupied by N i noninteracting bosons (filled circles). If indicated by an arrow, bosons may undergo transitions from state E j to state E i at a rate G i'j ¼ r ij (N i þ 1)N j , with rate constant r ij . In the language of evolutionary game theory, the figure depicts the interaction of N i agents (filled circles) playing strategies E i (bowls) for i ¼ 1,...,S. An agent playing strategy E j adopts strategy E i at a rate G i'j ¼ r ij N i N j . The above rate of bosonic condensate selection is recovered if agents may also spontaneously mutate from E j to E i at a rate r ij . (b) A condensate vector c for an antisymmetric matrix A has two properties: its entries are positive for indices for which Ac is zero, and they are zero for indices for which Ac is negative (''-'' signifies the antisymmetry of matrix A). Temporal evolution of the relative entropy of the condensate vector to the state concentrations under the ALVE (3) relates positive entries of the condensate vector to condensates, and its zero entries to depleted states. Generically, positive entries of c represent the asymptotic temporal average of oscillating condensate concentrations according to the ALVE (3), and negative entries of Ac represent depletion rates.
are described by the coupled birth-death processes (1) with rates (2). Before proceeding to the mathematical and numerical analysis of condensation in these processes, we now give a brief overview of such systems.

Results
Non-interacting bosons in driven-dissipative systems. The classical master equation (1) has recently been derived by Vorberg et al. 19 in the study of bosonic systems that are dissipative and driven by external sources. For a system of noninteracting bosons that is weakly coupled to a reservoir and driven by an external time-periodic force (a so-called Floquet system) [22][23][24] , one can eliminate the reservoir degrees of freedom (Born and Markov approximation) 25,26 and the density matrix of the system becomes diagonal (see the Supplement of the work of Vorberg et al. 19 ). The effective dynamics of the bosons become incoherent and are captured on a macroscopic level in terms of the coupled birth-death processes (1) with rates (2)). These non-equilibrium set-ups may not only lead the bosons into a single, but also into multiple condensates 19 .
For the incoherently driven-dissipative systems described above, the state E i denotes a time-dependent Floquet state [22][23][24] . The total rate G i'j for the transition of a boson from state E j to E i depends linearly on the number of bosons in the departure state (N j ) and the arrival state (N i þ 1). The latter factor stems from the indistinguishability of bosons and reflects their tendency to congregate. Although we refer to equation (1) as a classical master equation and coherence does not build up, the quantum statistics of bosons is still encoded in the functional form of G i'j . The rate constant r ij is determined by the microscopic properties of the system and the reservoir.
Condensation in the above set-up is to be distinguished from Bose-Einstein condensation. Typically, studies on Bose-Einstein condensation focus on the existence of long-range phase coherence in thermal equilibrium [10][11][12][13][14][15] , its kinetic formation 13,14,[27][28][29][30][31][32] and the fragmentation of a coherent condensate into multiple condensates (for example, when the equilibrium ground state is degenerate) 13,16 . In contrast, the classical birth-death processes (1) with rates (2) describe condensation in bosonic systems that are externally driven by a continuing supply of energy, dissipate into the environment and exhibit decoherence.
Equations of type (1) may also arise in atomic physics and quantum optics and are known as Pauli master equations [33][34][35] . They describe how the population of S non-degenerate energy levels changes over time when a system harbours N indistinguishable, non-interacting bosonic atoms. Such changes may occur by interactions with a radiation field that induces transitions between energy levels. A theoretical description of these transitions in terms of a Pauli master equation is appropriate if coherence is negligible. As in the previous example, the system then approaches a state in which some of the energy levels are macroscopically occupied (condensates), whereas the others are depleted. More generally, whenever a rate constant r ij governs the transition of a single boson from a state E j to E i , the rates (2) with s ij ¼ 1 for all i and j apply if N non-interacting bosons are brought into the system 19 .
Strategy selection in evolutionary game theory. The classical master equation (1) also occurs in evolutionary game theory (EGT). Historically, EGT was developed to study the evolutionary processes that are driven by selection and mutation 36,37 and seeks to identify optimal strategies for competitive interactions. For example, EGT has been applied in the study of the prominent 'rock-paper-scissors' (RPS) game, which was proposed as a facilitator of species coexistence and has inspired both experimental and theoretical research [38][39][40][41][42] . Furthermore, the 'prisoner's dilemma' game serves as a paradigmatic model to explore the evolution and maintenance of cooperation 43,44 . The interplay between non-linear and stochastic effects underlies the dynamics of such evolutionary games [45][46][47][48][49][50] .
In EGT, one typically considers a system of N interacting agents (classical particles) who repeatedly play one fixed strategy E i out of the S possible choices E 1 ; E 2 ; . . . ; E S . In each succeeding interaction, the defeated agent adopts the strategy of its opponent. Since N j agents playing strategy E j can potentially be defeated by one of the N i agents playing strategy E i , the rate of change is If an agent who plays E j can also spontaneously mutate into an agent who plays E i (with rate m ij ¼ r ij s ij ), one recovers the classical master equation (1) with rates (2).
Thus, there exists a correspondence between condensation in incoherently driven-dissipative bosonic systems and strategy selection in EGT-the transition of bosons between states can be interpreted in terms of the interaction and mutation of agents employing evolutionary strategies. In effect, the states in an incoherently driven-dissipative set-up play an evolutionary game and the winning states form the condensates.
After having introduced the above examples, we now proceed with the mathematical and numerical analysis of the classical master equation (1). We show that the dynamics of condensation change on two distinct timescales. At the leading-order timescale, the dynamics are described by a set of non-linearly coupled, ordinary differential equations (see equation (3) below), which determine the states that become condensates. We identify these states by applying concepts from EGT. After an exposition of the physical principles that underlie the condensation dynamics, implications of our general results for incoherently drivendissipative systems are discussed.
The antisymmetric Lotka-Volterra equation. The total number of particles needed for condensation phenomena to occur is large (N ) 1). To detect macroscopic occupancies, it is also assumed that the particle density N/S is large. Therefore, one may approximate the classical master equation (1) by a Langevin equation for the state concentrations x i (t) ¼ N i (t)/N (details of the derivation are provided in Supplementary Note 1). Originally proposed for Brownian particles suspended in a liquid, the Langevin equation decomposes the dynamics of a sample trajectory of the random process into two contributions-into a deterministic drift and into noise stemming from the discreteness of particle numbers ('demographic fluctuations'). Both the demographic fluctuations and the contribution to the deterministic drift that corresponds to mutations in the EGT setting are suppressed by a small prefactor 1/N. Therefore, these terms change the dynamics only slowly. The deterministic drift that corresponds to interactions between agents is, however, not suppressed. It thus governs the dynamics to leading order.
Hence, we find that the leading-order dynamics of the condensation process (1)-(2) are described by the differential equations: d dt The matrix A is antisymmetric and encodes the effective transition rates between states (a ij ¼ r ij -r ji ). The constants s ij that occur in the definition of the rates (2) do not change the leadingorder dynamics, but they become relevant on subleading-order timescales.
We refer to equation (3) as the antisymmetric Lotka-Volterra equation (ALVE). It provides a description of pairwise interactions that preserve the total number of particles. Therefore, the ALVE finds a broad range of applications in diverse fields of research, in addition to the aforementioned condensation of bosons far from equilibrium. It was first studied by Volterra 51 in the context of predator-prey oscillations in population biology 47,52,53 . In plasma physics, the ALVE describes the spectra of plasma oscillations (Langmuir waves) 54,55 , and in chemical kinetics it captures the dynamics of bimolecular autocatalytic reactions 18,[56][57][58] . In EGT, the ALVE is known as the replicator equation of zero-sum games such as the RPS game 47,59-61 . Table 1 summarizes all of the above analogies.
Despite the simple structure of the ALVE, it exhibits a rich and complex behaviour. In the following, we show how the mathematical analysis of the ALVE explains condensation into multiple states (condensate selection). To this end, we extend an approach for the analysis of the ALVE that was introduced in the context of EGT 59,60 .
Production of relative entropy and condensate selection. Our analysis starts from a theorem in linear programming theory 62 . Given an antisymmetric matrix A, it is always possible to find a vector c that fulfils the following conditions: its entries are positive for indices in I f1; . . . ; Sg and zero for indices in I ¼ f1; . . . ; Sg À I, whereas the entries of Ac are zero for indices in I and negative for indices in I (Fig. 1b). Although several vectors c with these properties may exist, the index set I is unique and, thus, determined by the antisymmetric matrix A. Finding such a 'condensate vector' c is crucial for the understanding of condensate selection and of the condensation dynamics. The condensate vector has the following physical interpretations.
All condensate vectors yield fixed points of the ALVE (3). Because of the antisymmetry of matrix A, a linear stability analysis of these fixed points does not yield insight into the global dynamics (Supplementary Note 1). However, the global stability properties can be inferred by showing that the relative entropy of a condensate vector to the state concentrations, is a Lyapunov function (note that we do not consider the relative entropy of the state concentrations to the condensate vector, but define the relative entropy vice versa). The relative entropy (4) decreases with time and is bounded from below (see Methods and Supplementary Fig. 1). Therefore, the dynamics relax to a subsystem in which the relative entropy production is zero. The relaxation of relative entropy production is reminiscent of Prigogine's study of open systems in non-equilibrium thermodynamics. Indeed, we find that the system, to cite Prigogine's phrase, 'settles down to the state of least dissipation' 63 . This state of least dissipation is characterized even further by the condensate vector c. Considering the definition of the relative entropy (4) and its boundedness, it follows that every concentration x i with iAI remains larger than a positive constant. On the  The concentration of the state associated to the purple disk in a decays exponentially to a concentration of 1.5 Â 10 À 7 before recovering transiently. Numerical integration cannot rule out permanent recovery at later times. Supplementary Fig. 2 demonstrates such a case.
other hand, states with indices in I become depleted for long times (see Methods). Therefore, we find that the condensate vector determines the selection of condensates. Positive entries of c correspond to states that become condensates, whereas zero entries of c correspond to states that become depleted. Both the set of condensates and the set of depleted states are unique (Fig. 1b) and independent of the initial conditions. Generically, the entries of the condensate vector are also unique upon normalization (its entries sum up to one) and yield the rate |(Ac) i | at which a state E i becomes depleted. The condensate selection occurs exponentially fast (see Fig. 2 and Methods). After relaxation, the dynamics of the system are restricted to the condensates. In other words, the condensates form the attractor of the dynamics. However, the dynamics in this subsystem do not come to rest. The state of least dissipation is a dynamic state with a perpetually changing number of particles in the condensates-periodic, quasiperiodic and non-periodic oscillations are observed ( Fig. 2b and Supplementary Fig. 1). In the generic case, the entries of the condensate vector represent the temporal average of condensate concentrations according to the ALVE (3). After condensate selection, the dynamics of these active condensates take place on a high-dimensional, deformed sphere 61 .
An algebraic algorithm to find the condensates. Numerical integration of the ALVE (3) is neither a feasible nor a reliable method for identifying condensates (Fig. 2, Supplementary Figs 1  and 2). Instead, we determine these states by numerically searching for a condensate vector c. To this end, we reformulate the above conditions on c in terms of two linear inequalities 62 : Ac 0 and c À Ac40 : We solve these inequalities with a linear programming algorithm that is both reliable and efficient. The time to find a condensate vector scales only polynomially with the number of states S (see Supplementary Fig. 3 and Methods for details).
Condensation in large random networks of states. We used our combined analytical and numerical approach to study how the connectivity of a random network of states affects the selection of condensates under the dynamics of the ALVE (3). The connectivity specifies the percentage of states between which particle transitions occur 20,64 . After having generated a network with a given connectivity, the strength and direction of an allowed transition between states E j and E i were determined by randomly sampling the corresponding effective rate constant a ij ¼ r ij -r ji .
Our results for condensation in large random networks of states are summarized in Fig. 3. When the connectivity of a network is zero, all of its states are isolated. Particles are not exchanged between states and none of the states becomes depleted. For an increased connectivity, isolated pairs of states are sampled in a random network. One state in an isolated pair is always depleted and the average number of condensates decreases   20 . We observe that, under the dynamics of the ALVE (3), the average number of condensates becomes minimal for a connectivity that also scales inversely with the number of states (see Fig. 3g). We attribute this minimum to the interplay between the criticality of random networks and condensate selection on connected components of the network 61,65 . Embedded directed cycles are a recurring motif 21 in the remaining network of condensates after condensate selection. Above the critical connectivity, a single giant cluster is formed. On average, half the number of states in this giant cluster become condensates once the network is fully connected (C ¼ 1; Fig. 3a) 19,60,66 . Thus, our analysis emphasizes the importance of critical properties of random networks for condensate selection.
Design of active condensates. Our understanding of the condensate selection can be used to design systems that condense into a particular network of states, a game of condensates. We exemplify this procedure by formulating conditions under which a system relaxes into a RPS game of condensates 40 . Three particular states E 1 , E 2 and E 3 in a system become a RPS cycle of condensates if, and only if, the following two conditions are fulfilled (Fig. 4). First, the 'RPS condition' requires that the rate constants between the three states form a RPS network (for example, r 12 4r 21 , r 23 4r 32 and r 31 4r 13 ). Second, the 'attractivity condition' requires that the inflow of particles into the RPS cycle from any other state E k is greater than the outflow to that state E k (for all k ¼ 4; . . . ; S). The values of the rate constants between the states that become depleted are irrelevant. More complex games of condensates can be designed by formulating similar conditions on the rate constants. These conditions are formulated as inequalities that depend on Pfaffians of the antisymmetric matrix A and its submatrices (see Methods) 61 . The flow of particles between states in these systems causes condensate concentrations to oscillate (Fig. 2b).

Discussion
Our findings thus suggest intriguing dynamics of condensates in systems whose temporal evolutions are captured by the classical master equation (1) with rates (2), for example, in drivendissipative systems of non-interacting bosons. Condensates observed on the leading-order timescale are metastable. For longer times, relaxation into a steady state occurs 19,67 . When detailed balance is broken in the system of condensates, the net probability current between at least two states does not vanish and a non-equilibrium steady state is approached 68,69 . The simplest way of designing such condensates is illustrated by the above RPS game. In this game, detailed balance is broken, for example, when the transition of particles is unidirectional (with totally asymmetric rate constants r 12 4r 21 ¼ 0, r 23 4r 32 ¼ 0, and r 31 4r 13 ¼ 0). For non-interacting bosons in driven-dissipative systems, the continuing supply with energy through the external time-periodic driving force (Floquet system) and the dissipation of energy into the environment may, therefore, prevent the system from reaching equilibrium. How such systems may be realized in an experiment poses an interesting question for future research.
The transition of particles between condensates in the herestudied coupled birth-death processes parallels the interaction and mutation of winning agents in evolutionary game theory, reflecting an 'evolutionary game of condensates'. Our results suggest the possibility of creating novel bosonic systems with an oscillating occupation of condensates. Non-interacting bosons in incoherently driven-dissipative systems are promising candidates. Since the antisymmetric Lotka-Volterra equation also arises in population biology, chemical kinetics and plasma physics, all of our mathematical results apply to these fields as well.

Methods
Asymptotics of the ALVE. The asymptotic behaviour of the ALVE (3) can be characterized as follows: for every antisymmetric matrix A there exists a unique subset of states I f1; . . . ; Sg whose concentrations stay away from zero for all times, that is, x i ðtÞ ! ConstðA; x 0 Þ40 for all t ! 0 and for every i 2 I : ð6Þ The set I is the set of condensates. All of the other states with indices in I ¼ f1; . . . ; Sg À I become depleted as t-N, that is, The set of condensates can be determined algebraically from the antisymmetric matrix A and does not depend on the initial conditions x 0 2 D S À 1 ¼ fx 2 R S j x i 40 for all i; P S i¼1 x i ¼ 1g. To show this result, the time-dependent entropy D(c||x)(t) of a condensate vector c ¼ ðc 1 ; . . . ; c n Þ 2 D S À 1 (c i Z0 for all i and P i c i ¼ 1) relative to the trajectory x(t) is considered (that is, the Kullback-Leibler divergence of x(t) from c), see equation (4). A condensate vector is defined via the properties (see Fig. 1b Such a vector can always be found for an antisymmetric matrix 62 . Notably, the index set I is unique although more than one condensate vector may exist. Considering the time derivative of the relative entropy D(c||x)(t) and employing equations (3) and (8) yields: Since ðAcÞ I o0 and x40, it follows that @ t D(c||x)(t)o0 (please note the overbars in subscripts, which may be lost when read in low resolution). Therefore, the relative entropy D(c||x) is a Lyapunov function if c is chosen in accordance with equations (8) and (9). Moreover, D(c||x) is bounded from above by D(c||x)(0) and from below by zero for all times. This can be seen from the definition of D, and  Figure 4 | Conditions for the emergence of a RPS cycle of condensates.
Three particular states E 1 , E 2 and E 3 (blue, red and yellow disks) of a network condense into a RPS cycle if, and only if, two conditions are fulfilled: First, the 'RPS condition' requires that the rate constants r ij between the three states form a RPS network: r i À 1,i þ 1 4r i þ 1,i À 1 (indices are counted modulo 3, for example, r 42 ¼ r 12 (framed arrows denote rate constants that are larger than rate constants for the respective reverse direction). Differences between these rate constants define the entries c i ¼ r i À 1,i þ 1 -r i þ 1,i À 1 of an admissible condensate vector c. Second, the 'attractivity condition' requires that the weighted sum of rates from any exterior state E k (purple disks) into the RPS cycle, P 3 j¼1 c j r jk (framed arrows), is larger than the weighted sum of outbound rates, P 3 j¼1 c j r kj (black arrows). In other words, the inflow of particles into the RPS cycle from any exterior state needs to be greater than the outflow to that state. From the definition of the relative entropy in equation (4), it follows that every concentration x i with iAI remains larger than a positive constant, that is, x i (t)ZConst(A, x 0 )40 for all times t (if x i (t)-0 for iAI, it follows that D-N, which contradicts the boundedness of D).
Furthermore, equation (11) implies that, for every i 2 I and for all t. Therefore, concentration x i is integrable for every i 2 I (x i AL 1 (0,N)) with the bound: 0o Since the derivative of the concentrations is bounded from above and below, , one concludes that x i is uniformly continuous (|| A || N-N denotes the operator norm of A induced by the maximum norm on R S ). Together with the integrability (13), it follows that the states with indices in I become depleted as t-N, that is, x i (t)-0 for i 2 I.
In conclusion, given an antisymmetric matrix A ¼ R-R T via a rate constant matrix R ¼ {r ij } i,j , one finds a condensate vector c that satisfies inequalities (8)- (9). The index set I, for which the entries of c are positive, represents condensates. The index set I, for which entries of c are zero, represents states that become depleted. Moreover, equation (10) implies that the relative entropy becomes a conserved quantity in the subsystem of condensates ( Supplementary Fig. 1).
Temporal average of condensate concentrations. The ALVE (3) is solved implicitly by, with the time average of the trajectory hxi t defined as: It is shown above that 0oConst(A, x 0 )rx i (t)r1 holds for the states that become condensates (iAI). By rearranging equation (14), one thus obtains: Note that Const is used to denote arbitrary positive, time-independent constants. Therefore, the right-hand side of equation (16) vanishes for t-N. On the other hand, x i is integrable for i 2 I (equation (13)). Thus, the corresponding component of the time average converges to zero, Hence, the distance of the time average hxi t to the kernel of the antisymmetric submatrix A I converges to zero (the submatrix A I corresponds to the system of condensates with indices in I).
Structure of a generic antisymmetric matrix. For systems with an even number of states S, the antisymmetric matrix A ¼ R-R T generically has a trivial kernel, whereas for systems with an odd number of states, the kernel of A is generically one dimensional. A higher dimensional kernel of A only occurs if the matrix entries are tuned 52,61,70 . As a consequence, when all of the entries above the diagonal of A are, for example, randomly drawn from a continuous probability distribution (for example from a Gaussian distribution), all 2 S submatrices of A have a kernel with dimension of less than or equal to one. The projection of x 2 R S to the subspace R ðJÞ R S for an arbitrary index set J f1; . . . ; Sg is defined as x J :¼ P J x:¼ (x j ) jAJ . In other words, entries of x J are zero for indices in the complement J. In the following, the short notation A J :¼ P J AP J is also used (see above). Furthermore, the set of antisymmetric matrices whose submatrices have a kernel with dimension r1 is defined: O :¼ A 2 R SÂS j A is antisymmetric and dim ker A J 1 for all J f1; . . . ; Sg È É : The complement O has measure zero with respect to the flat measure dA on antisymmetric matrices (the translation invariant measure, which is sigma-finite and not trivial).
For antisymmetric matrices AAO, the kernel can be characterized as follows 61,70 . If the number of states S is even, the kernel of A is trivial: ker A ¼ {0}. If the number of states is odd, the kernel is one-dimensional: ker A ¼ {v}. This kernel element can be computed analytically in terms of Pfaffians of submatrices of A: v ¼ Pf A1 ð Þ; À Pf A2 ð Þ; . . . ; Pf AŜ À Á À Á : ð19Þ The submatrix Ak 2 R ðS À 1ÞÂðS À 1Þ denotes the matrix for which the k-th column and row are removed from A. For antisymmetric matrices AAO, the normalized condensate vector c with P i c i ¼ 1 and with properties (8), (9) is unique. The latter follows from A I c ¼ 0 (equations (8) and (18)). Therefore, the condensate vector is the unique kernel vector of the subsystem of condensates whose interactions are characterized by the matrix A I . Furthermore, I contains an odd number of elements. To determine the condensate vector for AAO, one can proceed as follows. For each odd-dimensional submatrix A I with I f1; . . . ; Sg, one computes the kernel element v according to equation (19) and defines the vector w 2 R S by setting w I ¼ v and w I ¼ 0. There exists exactly one set I for which ðAwÞ I o0. The corresponding vector w is the unique condensate vector upon normalization.
Temporal average of condensate concentrations (generic case). It was shown above that the temporal average of condensate concentrations hxi t converges to a non-negative kernel element of the antisymmetric matrix A I . In the generic case, the condensate vector c is the unique kernel element of A I upon normalization. Therefore, positive entries of c represent the asymptotic temporal average of condensate concentrations, Exponentially fast depletion of states (generic case). On inserting equation (20) into the implicit solution (14) of the ALVE, the exponentially fast depletion of states with i 2 I can be seen as follows (note that (Ac) i o0 according to the choice of the condensate vector in equations (8) and (9)): x i ð0Þe tÁ ðAcÞ i þ kAðhxi t À cÞk1 ð Þ ð22Þ x i ð0Þe tÁðAcÞ i þ ConstðA;x0Þ ð23Þ and analogously, Therefore, condensate selection occurs exponentially fast at depletion rate |(Ac) i |. The dynamics of cases for non-generic antisymmetric matrices are discussed in Supplementary Note 2.