On the number of Bose-selected modes in driven-dissipative ideal Bose gases

# On the number of Bose-selected modes in driven-dissipative ideal Bose gases

Alexander Schnell Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Straße 38, 01187 Dresden, Germany    Roland Ketzmerick Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Straße 38, 01187 Dresden, Germany Technische Universität Dresden, Institut für Theoretische Physik and Center for Dynamics, 01062 Dresden, Germany    André Eckardt Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Straße 38, 01187 Dresden, Germany
July 14, 2019
###### Abstract

In an ideal Bose gas that is driven into a steady state far from thermal equilibrium, a generalized form of Bose condensation can occur. Namely, the single-particle states unambiguously separate into two groups: the group of Bose-selected states, whose occupations increase linearly with the total particle number, and the group of all other states whose occupations saturate [Phys. Rev. Lett. 111, 240405 (2013)]. However, so far very little is known about how the number of Bose-selected states depends on the properties of the system and its coupling to the environment. The answer to this question is crucial since systems hosting a single, a few, or an extensive number of Bose-selected states will show rather different behavior. While in the former two scenarios each selected mode acquires a macroscopic occupation, corresponding to (fragmented) Bose condensation, the latter case rather bears resemblance to a high-temperature state of matter. In this paper, we systematically investigate the number of Bose-selected states, considering different classes of the rate matrices that characterize the driven-dissipative ideal Bose gases in the limit of weak system–bath coupling. These include rate matrices with continuum limit, rate matrices of chaotic driven systems, random rate matrices, and rate matrices resulting from thermal baths that couple to a few observables only.

## I Introduction

Other than in thermodynamic equilibrium, the density matrix of a gas that is driven into a nonequilibrium steady state does not follow directly from a simple principle like entropy maximization. It depends on the microscopic details of the system, the bath, and the system–bath coupling. Lately, a high interest into out-of-equilibrium steady states which remember the initial condition of an isolated system has developed. For example, such steady states were observed in many-body localized systems Basko et al. (2006); Nandkishore and Huse (2015); Schreiber et al. (2015); Smith et al. (2016); Choi et al. (2016), including discrete time crystals Khemani et al. (2016); von Keyserlingk et al. (2016); Else et al. (2016); Zhang et al. (2017); Choi et al. (2017) that occur in many-body localized time-periodically driven (Floquet) systems. Also nonequilibrium steady states (NESS) of driven-dissipative many-body systems have generated much attention. This includes open Floquet systems Tsuji et al. (2009); Vorberg et al. (2013, 2015); Foa Torres et al. (2014); Seetharam et al. (2015); Dehghani et al. (2015); Goldstein et al. (2015); Iadecola et al. (2015); Shirai et al. (2016); Qin and Hofstetter (2017); Chong et al. (2017), state engineering via dissipation Diehl et al. (2008); Letscher et al. (2017); Labouvie et al. (2016); Schnell et al. (2017) and photonic systems Deng et al. (2010); Carusotto and Ciuti (2013), where coherent phenomena in some limits may be related to lasing but in other limits rather to Bose condensation Szymańska et al. (2006); Kirton and Keeling (2013); Chiocchetta and Carusotto (2014); Klaers et al. (2010); Byrnes et al. (2014); Leymann et al. (2017). In this context, optical microcavitities offer great freedom for designing system and dissipative environment, which allows for tailoring the coherent emission of multiple modes Türeci et al. (2006); Cao (2003) or to control a switching between emission of two (or more) different modes Sondermann et al. (2003); Marconi et al. (2016); Ge et al. (2016); Leymann et al. (2017).

In this paper, we are focusing on driven-dissipative ideal gases of noninteracting bosons that exchange energy with the environment but no particles. They can be driven out of equilibrium, e.g. by periodic driving in combination with the coupling to a heat bath or by coupling the system to two heat baths of different temperature. In such setups the ideal gas will relax to a nonequilibrium steady state (NESS) which is characterized by a finite heat current through the system. It is an interesting question, whether (and if yes when and in which form) a system can show Bose condensation (or other forms of ordering) under such nonequilibrium conditions. Since this NESS does not follow from thermodynamic principles it is not obvious whether such a state will feature Bose condensation or not. It was observed in Ref. Vorberg et al. (2013) that in the quantum degenerate limit of large densities the single-particle states split into two groups; the Bose-selected states, whose occupations increase linearly with the total particle number, much like for the ground state in thermal equilibrium, while the occupations of all other states saturate.

However, so far very little is known about the factors that determine the number of Bose-selected states. The answer to this question is crucial since systems hosting a single, a few, or an extensive number of Bose-selected states will show rather different behavior. While in the former two scenarios each selected mode acquires a macroscopic occupation, corresponding to (fragmented) Bose condensation, the latter case rather bears resemblance to a high-temperature state of matter. Moreover, inducing transitions between several condensate modes can be a very efficient mechanism to exchange energy with the environment, which is not present in systems hosting a single condensate only Vorberg et al. (2013).

In this paper, we investigate how the number of Bose-selected states depends on the properties of the system and its coupling to the environment. To this end, we study several different scenarios, which are reflected in different forms of the rate matrices that characterize the driven-dissipative ideal Bose gases in the limit of weak system-bath coupling. These include rate matrices with continuum limit, rate matrices of chaotic driven systems, random rate matrices, and rates matrices resulting from thermal baths that couple to a few observables only.

## Ii Driven-dissipative ideal Bose gas and Bose selection

In the limit of weak system–bath coupling, the system will approach a nonequilibrium steady state that is diagonal in the eigenstates of the Hamiltonian for an autonomous (i.e. non-driven) system or in the Floquet states for a time-periodically driven system Breuer and Petruccione (2002); Kohler et al. (1997); Vorberg et al. (2015). The mean occupations of these states obey the equation of motion Vorberg et al. (2013)

 ∂t⟨^ni⟩=∑jRij⟨(^ni+1)^nj⟩−Rji⟨(^nj+1)^ni⟩=0. (1)

The rate for a boson to jump from single-particle level  to is given by the single-particle rate multiplied by the bosonic enhancement factor which manifests that bosons favor to “jump” into states that already have large occupation.

It was pointed out in Refs. Vorberg et al. (2013, 2015) that a generalization of Bose condensation is also observed in the NESS, called Bose selection. Here, a whole group of an odd number of single-particle states, the Bose-selected states, can acquire large occupation. As we briefly recapitulate in this section, these selected states are only determined by the rate asymmetry matrix

 Aij=Rij−Rji. (2)

We assume that the gas may exchange heat with an environment of one or more thermal phonon baths. A single bath is described as a collection of harmonic oscillators which are in thermal equilibrium. The corresponding system–bath coupling operator reads with dimensionless system coupling operator , coefficients and coupling strength . Throughout the paper we assume that the baths are Markovian. Thus, the single-particle rate for the autonomous system,

 Rij=2πγ2ℏ|⟨i|^v|j⟩|2g(εi−εj), (3)

is of golden rule type. Here, is the energy of single-particle eigenstate . It enters the bath-correlation function,

 g(ε)=J(ε)eε/T−1, (4)

where is the temperature (measured in units of energy, ) of the bath and is its spectral function . In the following we will consider ohmic baths with a continuum of modes and spectral function .

A nonequilibrium situation is found when e.g. the system is coupled to multiple baths at different temperatures, where the total rate is given by the sum of the rates corresponding to the individual bath ,

 Rij=∑bR(b)ij. (5)

Another possible scenario are time-periodically driven systems coupled to a heat bath. In this case the rates read

 Rij=2πγ2ℏ∞∑m=−∞|vij(m)|2g(εi−εj−mℏω) (6)

with , driving frequency , Floquet states , and corresponding quasienergies .

The equation of motion (1) for the mean occupations depends also on the two-particle density–density correlations, the equations of which depend in turn on three-particle correlations and so on. In this way it establishes a hierarchy, which in the following will be truncated already at the single-particle level by employing the mean-field decomposition , . As a result, the steady state occupations, , follow from the nonlinear equations of motion

 0=∑jAij⟨^nj⟩⟨^nj⟩+Rij⟨^nj⟩−Rji⟨^ni⟩. (7)

Here we have used the rate asymmetries , which for the rates of a single bath read

 (8)

so that they are independent of temperature. Note that it has been shown by comparison to quasi-exact Monte-Carlo simulations Vorberg et al. (2015) that the mean-field equation (7) yields excellent predictions for the mean occupations  for a broad range of models111A possible reason for this good agreement has been pointed out recently Langen et al. (2015); Lange et al. (2017); a driven system with a set of observables which are approximately conserved quantities will relax towards a steady state that is well described by a generalized Gibbs ensemble Due to the weak-coupling limit that we assume, the occupations  of the system’s single-particle states  are almost conserved, so in our context this set is given by . For a state , the mean-field decomposition is exact, , ..

The solid lines in Fig. 1 show the steady-state solutions of Eq. (7) as a function of the total particle number for three different scenarios: Fig. 1(a) shows occupations for a tight-binding chain of sites, coupled to one heat bath only; therefore, the steady state is thermal. Figure 1(b) shows occupations for the same chain, but additionally in contact also with a second, population inverted heat bath described by a negative temperature . Note that negative temperatures have been realized, e.g., in atomic quantum gases by preparing a state at the upper edge of a Bloch band Braun et al. (2013). Figure 1(c) shows occupations for a time-periodically driven quantum kicked rotor with Floquet states coupled to a single bath. The system is in a regime, where the corresponding classical system, the Chirikov standard map Chirikov (1979), is known to be chaotic.

For small total particle number , the bosons behave classically and the occupations are given by the single-particle probabilities to occupy the state  scaled linearly with particle number , . However, at large total particle numbers the bosonic quantum statistics makes itself felt. As a result, we observe Bose selection Vorberg et al. (2013): the occupations of some of the states saturate, while all additional particles gather in a set of selected states whose occupations grow linearly with . In equilibrium, Fig. 1(a), this corresponds to Bose condensation in the ground state. Away from equilibrium several states can be Bose selected.

From Figs. 1(a)-1(c) we already observe that the number of selected states can range from only few up to an extensive number, while the former case corresponds to fragmented Bose condensation, since each of the selected state acquires a macroscopic occupation in the limit , the latter case does not correspond to Bose condensation, since none of the selected states will acquire a macroscopic occupation. To be more precise, in the thermodynamic limit, , , there can only be (fragmented) condensation if the number of Bose-selected states is intensive, i.e. asymptotically independent of . However, if there is an extensive number of states, , the system will behave effectively classically also in the ultra degenerate limit. So far, however, very little is known about how depends on the properties of the system and in particular on the form of rates. The main goal of this paper is to obtain a better understanding of how the number of selected states is determined by the properties of the rate matrix.

Before we begin with our analysis, let us briefly review the equations that determine the set of selected states and what so far has been known about their number. It has been shown that generally the number of these selected states is odd. The starting point for determining the set of selected states is an asymptotic expansion (dashed lines in Fig. 1) of the mean occupations in the limit of large occupation, . For the selected states, , Eq. (7) yields in this limit Vorberg et al. (2013, 2015)

 0=∑j∈SAijηj,∀i∈S, (9)

where are the leading order occupations . Note that for the nonselected states, these contributions vanish, so it must hold . The vector , is thus a nontrivial vector in the kernel of the skew-symmetric matrix . The leading order of the occupations of the nonselected states, , is then given by

 ⟨^ni⟩=−∑j∈SRijηj∑j∈SAijηj+O(N−1). (10)

The physical condition of having positive occupations, for both the selected and the nonselected states, has been shown to uniquely determine the set of the selected states Vorberg et al. (2015). Using equations (9) and (10) this condition can be cast into the simple form

 μ=Aη with {μi=0,ηi>0, % for i∈Sμi<0,ηi=0, for i∉S. (11)

From this condition a few simple statements about the number of selected states have been drawn already. First, we know that is generically odd, since without fine tuning the skew-symmetric possesses a non-trivial kernel for an odd number of selected states only. Second, in case the system possesses a ground-state like state defined by

 Rki>Rik,∀i≠k, (12)

one can immediately see that the problem (11) is solved by , i.e. a single selected state is found. Finally, for uncorrelated random rates it has been observed numerically that the number of selected states follows a binomial distribution. Thus, an extensive number of states (on average half of the states) are selected.

However, a general estimate for is not straight forward. In the following, we will discuss different scenarios in which such estimates can be found; First we will study rates which possess a well-defined continuum limit for . Second we will consider rates that do not have such a continuum limit, discussing the two important cases where rates are truly random, and rates that stem from a chaotic kicked system. Finally, we will discuss rates that are given by a sum of direct products, as they are relevant for autonomous systems that couple to the environment via a few observables only.

## Iii Rates with continuum limit

In this section we discuss systems described by rate matrices that have a smooth continuum limit and which are, thus, strongly correlated. We assume that the quantum numbers can be labeled by a variable

 ki=αiM, with constant α∈R, (13)

that becomes continuous in the limit . Moreover, we focus on one-dimensional systems whose rate matrices shall become smooth in that limit, i.e. that there exists a function , such that

 Rij=R(ki,kj)Δ2 (14)

with . The generalization to higher-dimensional systems described by several continuous quantum numbers is straight forward.

An example for this situation are the rates for a single bath coupled to site of the tight-binding chain described by the Hamiltonian [see Fig. 1(b)]

 HS=−JM−2∑i=1(^a†i+1^ai+^a†i^ai+1), (15)

where is the tunneling constant and is the bosonic annihilation operator at site . In this case we find Schnell et al. (2017)

 R(k,q)=2γ2πℏg(ε(k)−ε(q))sin(kℓ)2sin(qℓ)2 (16)

with dispersion and -space sampling with , and .

We make the ansatz that there is a discrete set of Bose-selected states , such that in the asymptotic limit for large densities the mean occupation density reads

 ⟨n(k)⟩=⟨nn(k)⟩+¯nsδ(k−ks), (17)

where we have introduced the normalization condition

 n=NΔ=∑i⟨ni⟩Δ⇒n=∫α0⟨n(k)⟩dk. (18)

Eq. (9) then translates to

 a(ks)=0,∀ks∈S, (19)

where we have defined the function

 a(k)=∑s∈SA(k,ks)¯ns, (20)

with rate asymmetry function . Note that for smooth rates, the function will be smooth as well. The starting point for our reasoning is the analog of Eq. (10), which predicts that in the continuum limit the asymptotic occupations of all modes read

 ⟨nn(k)⟩=−∑s∈SR(k,ks)¯nsa(k). (21)

Since here the numerator is strictly non-negative, the denominator has to be strictly negative

 a(k)<0,∀k∉S. (22)

In the following, we assume that is two-fold differentiable. By discussing the rate function in the vicinity of the selected states, we will then be able to restrict the possible selected states to a few candidates only. From Eqs. (19) and (22) it follows that is negative almost everywhere, however at local maxima it assumes the value zero, whenever , see sketch in Fig. 2. This implies both

 0=a′(ks)=∑p∈SA′(ks,kp) ¯np, (23) and 0>a′′(ks)=∑p∈SA′′(ks,kp) ¯np, (24)

where we have defined and . Note that these local criteria are necessary, but not sufficient for Bose selection in state . Note also that state space may have a boundary or not. For example for the tight-binding chain in Fig. 1(b), we do not impose periodic boundary conditions in real space, then takes values in the interval . Thus the dispersion is not a periodic function of and our state space possesses boundaries at and . At such boundaries the criteria (23) and (24) do not have to apply, because here the maxima of are no longer characterized by derivatives.

As we will see, the criteria (23) and (24) strongly constrain the set of selected states. To illustrate this fact, we will discuss different scenarios in the following subsection.

### iii.1 Different selection scenarios

If only one state is selected and does not lie at the boundary of the state space, it follows that

 0=∂kA(k,q)|(k0,k0)and0>∂2kA(k,q)|(k0,k0). (25)

Since is skew-symmetric,

 A(k,q)=−A(q,k),∀k,q, (26)

we also find

 0=∂qA(k,q)|(k0,k0)and0<∂2qA(k,q)|(k0,k0). (27)

Since also the gradient of vanishes at , this point must be an extreme point of the asymmetry function . Moreover, we find that the mixed derivative at this point must vanish

 ∂k∂qA(k,q)|(k0,k0)=∂q∂kA(k,q)|(k0,k0)=−∂q∂kA(q,k)|(k0,k0)=−∂k∂qA(k,q)|(k0,k0)=0. (28)

Here we first used Schwarz’ theorem and in the third step renamed the variables. Thus either corresponds to a saddle point on the diagonal of the rate asymetry matrix or it lies at the boundary of state space (if there is one), where both local criteria do not have to apply.

An example is given by the rate-asymmetry function  shown in Fig. 3(a). Here we mark the position of the selected state by black arrows on the side and the relevant matrix element by a red arrow. It lies at a saddle point (having the correct curvature) on the diagonal. Clearly, also the rest of function must remain below the blue plane.

For the vector of the occupations is a homogeneous solution of Eq. (19) and at the same time of Eq. (23). This provides a strong restriction, since both equations generally will not have a common set of solutions. Therefore, selected states will generically occur in two special scenarios.

The first scenario is the following. Since the rate asymmetry is continuous we will naturally find zero lines also away from the diagonal . If the selected states lie at these zero lines, i.e. , the occupations are only determined by Eq. (23).

Note that since the coefficient matrix is not necessarily skew-symmetric, in this case also an even number of selected states may occur in the continuum limit. In the discrete case, an even number of selected states requires fine-tuning in the rate matrix Vorberg et al. (2015) such that for example some of the vanish. However, for continuous rate asymmetry functions it is natural to have zero lines, i.e.  for , such that in the continuous case no fine-tuning is needed to observe an even number . However, if an even number of selected states occurs in the continuous model, a corresponding discrete system will still feature an odd number of selected states. We then typically find pairs of neighboring selected states around at least one222Around which of the states pairs form is, however, not universal and depends on the discrete grid. If one plots the selected states as a function of discretization parameter for example, one can observe how such pairs jump from one to the other quite irregularly (not shown). of the selected states of the continuous model [as in the example in Fig. 3(b)]. To avoid confusion, we refer to the number of selected states in the continuous system as continuum number of selected states .

In Fig. 3(b) we observe Bose selection at zero lines where the continuum selection of states occurs. The asymptotic states lie at a zero line of . Since in the discrete system the number of selected modes is always odd, we find in the discrete system with a pair of neighboring states in the vicinity of state .

Another possible scenario is that the selected states are found such that at the rate asymmetry function possesses saddle points. Then such that the occupations are solely determined by Eq. (19). In this case the corresponding continuum will be odd.

However typically neither of the cases – selection at only zero lines or only saddle points – occurs in a “pure” form. This can be seen in the example in Fig. 3(c), where we observe the selection at both zero lines and saddle points with a continnum selection of states. In this case, some of the relevant points lie on the zero lines of , and others have saddle points in the vicinity of the point . Here, states are selected in the discrete system, with pairs at two continuum wave numbers .

But it is only in continuous rate matrices with few oscillations (like in the examples chosen in Fig. 3) that we find selected states defined by one of the different mechanisms stated above. If we consider systems with more variation in the rate matrix, like in Fig. 4, the points are not clearly relatable to either zero lines nor saddle points anymore. Nevertheless, these points are often still found in the vicinity of zero lines and saddle points.

### iii.2 Random-wave model

The rate functions shown in Fig. 3 and 4 are motivated by a random wave model for chaotic eigenfunctions Berry (1977). They are defined by

 R(k,q)=L∑l=1 Re{clexp[i(κk,lk+κq,lq)]}+C. (29)

It is a superposition of independent plane waves with , , fixed absolute value of the wavenumber and uniformly distributed angles . Amplitudes are drawn from a normal distribution and phases are distributed uniformly. The global constant is chosen such that all rates are non-negative. Its specific value is irrelevant for the set of the Bose-selected states, which are solely determined by the rate-asymmetry function in which will drop out.

This model produces rate functions that are smooth and show oscillations on the length scale in state space, see e.g. Fig. 3. It allows us to investigate the typical behavior of smooth rates that show variations on such a scale. We discretize state space via Eq. (13) where, to unravel the continuum physics, we choose the discretization length to be small with respect to the oscillation length of the model, .

As we increase the wavenumber of the model the structure of will become more and more complex, leading to more and more zero-lines and saddle points in the rate asymmetry function. Thus, we expect that with the characteristic wavenumber , which determines the typical oscillation length of the smooth rate function , also the average number of selected modes will increase. This is confirmed by the numerical results shown in Fig. 5(a). The mean number of selected states increases as soon as exceeds the threshold , where the wavelength of the plane waves is about double the system size. Afterwards, we see a linear increase with as marked by the black and white dashed line. The linear scaling is explained by the above reasoning, since e.g. the number of extrema (and also the number of zeros) of on the interval is of the order of . However the origin of the prefactor of about remains open.

Note that the mean number only depends weakly on the number of components that we use in the random wave model as shown in Fig. 5(b) for discrete states and (c) and three different values of .

The behavior seen in Fig. 5(a) clearly suggests that for smooth rates, the number of selected states is typically on the order of the number of oscillations in the rate asymmetry function . Therefore, for fixed smooth , even if the number of discrete states is large, , as observed in Fig. 5, the number will remain intensive as it is solely the property of the smooth function .

We expect a breakdown of the theory for continuous rates as soon as the oscillation length of the rate function becomes comparable to the discretization length. For the random rates this breakdown occurs for as visible in Fig. 5(a). Interestingly, the mean number of selected states for the random wave model does not saturate at as one would expect for truly random rates (see Sec. IV.1). We observe a first saturation at values that are slightly above : for discrete states we find saturation at about , or for we find .333Note that after this first saturated regime, for the red line , the number drops again at around . This is related to the effect of aliasing occuring when the wavelength of the random wave model is smaller than the discretization length .

## Iv Rates without continuum limit

### iv.1 Uncorrelated random rates

The single particle rates that one observes typically for a fully chaotic quantum rotor have been shown to roughly follow an exponential distribution Wustmann (2010). If we suppose that there are no additional correlations between the rates , we may draw typical rates from a random realization of the exponential distribution, where . Note that the choice of the parameter is irrelevant, since it only determines the time scale on which the system relaxes but not the steady state.

Figure 6(a) shows the distribution of the number of Bose-selected states for random rates connecting single particle states. The odd number of Bose-selected states is given by a binomial distribution,

 p(MS)=⎧⎪⎨⎪⎩0for MS even% 12M−1(MMS)for MS odd (30)

centered around [cf. Fig. 6(c)]. Such a behavior was observed in the literature for random rates, mostly in the context of population dynamics, where similar equations, the Lotka-Volterra equations, appear Fisher and Reeves (1995); Chawanya and Tokita (2002); Allesina and Levine (2011); Vorberg et al. (2013); Knebel et al. (2015). Usually this result is motivated by arguing that for the random rates all states are equal, so that every state has the same probability to be Bose selected or not. Taking into account the additional constraint that the total must be odd, one can in this way motivate the distribution (30) by just counting the number of possible choices of states among the total states.

However, the argument that every state has the same probability to be Bose selected or not must also follow rigorously from the steady state equations, and thus also from Eq. (11). To fill this gap, we will compute the distribution of selected states. To this end, we will first contruct transformations of the rate-asymmetry matrix under which the number of Bose-selected states remains invariant, i.e. the solutions of the problem (11) have the same .

#### iv.1.1 Reordering transformations

Since the indexing of the states is arbitrary, will remain invariant when reordering the states. The matrix describes the transpositions that exchange state and . Since the corresponding transformation takes the form

 A→Tk,lATk,l. (31)

Under this transformation, the matrix remains skew-symmetric, and also the number of Bose-selected states remains invariant.

So for discussing the properties of a general with a given number of Bose-selected states, we can thus assume solutions of the form

 η=(η1,...,ηMS,0,...,0)T. (32)

#### iv.1.2 Rescaling transformations

We are only interested in the number of selected states and not in their occupation. Therefore, let us consider a rescaling of the coefficients and . This is accomplished by a transformation induced by the matrix

 Dλ=diag(λ1,...,λM),with positive λk>0. (33)

As we want to apply this matrix to both the vectors and , it is important that only positive rescaling is allowed. Otherwise we would transform a solution of problem (11) into vectors that do not solve a problem of this type. Since the inverse of this matrix is again diagonal, we find that the rescaling transformation

preserves skew-symmetry of the matrix . Note that we multiply with the inverse from the left and the right. The rescaled quantities with fulfil again Eq. (11), since

#### iv.1.3 Standard matrices

As we know that for any given skew-symmetric the problem (11) has a unique solution Vorberg et al. (2013, 2015), we always find a sequence of transformations such that

 ~η=(1,...,1,MS entries0,...,0)T,~μ=(0,...,0,MS entries−1,...,−1)T (36)

holds. As the transformations (31) and (34) are invertible, it turns out that it suffices to discuss the properties of the set of standard matrices that is constructed such that it gives rise to a solution of the form (36). Then the space of rate-asymmetry matrices with selection number is spanned by a sequence of the two physically intuitive transformations discussed above.

We find that these standard matrices take the form

 AMS=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝x(1)MS+1⋯x(1)MASMSx(2)MS+1−x(1)MS+1⋯x(2)M−x(1)Mx(3)MS+1−x(2)MS+1⋯x(3)M−x(2)M⋮⋱⋮1−x(MS)MS+1⋯1−x(MS)M−x(1)MS+1−x(2)MS+1+x(1)MS+1⋯−1+x(MS)MS+1−x(1)MS+2−x(2)MS+2+x(1)MS+2⋯−1+x(MS)MS+2⋮⋮⋱⋮Aarb−x(1)M−x(2)M+x(1)M⋯−1+x(MS)M⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (37)

with arbitrary real numbers , arbitrary ()-dimensional skew-symmetric matrix and -dimensional skew-symmetric matrix restricted by the existence of the homogeneous solution

 ASMS⎛⎜ ⎜ ⎜ ⎜⎝1⋮1⎞⎟ ⎟ ⎟ ⎟⎠=0. (38)

#### iv.1.4 Consequence for random rates

For a random rate matrix, all entries of the rate asymmetry are random and statistically independent. We aim to find the probability to randomly choose a matrix that is from the class generated by the -matrix (under the transformations (31) and (34)). To this end, we first count the number of degrees of freedom in determining a general matrix . In a second step we then discuss the influence of the transformations.

Let us begin with the block . To construct such a matrix, we start from an arbitrary -dimensional skew-symmetric matrix which has degrees of freedom. We have to distinguish two different cases: For odd this matrix has always a homogeneous solution, for even we have to fine tune one parameter for a homogeneous solution to exist, which reduces the number of degrees of freedom by one. Furthermore we have to subtract degrees because the homogeneous solution is pinned to . Thus the number of degrees of freedom in the subspace of the selected modes is

 doS =12(MS−1)(MS−2), (39) deS =doS−1 (40)

for an odd or even number of selected modes, respectively.

Then, there are rows, each of which has free variables to choose. This contributes

 df=(M−MS)(MS−1) (41)

degrees of freedom.

Also there is still an arbitrary -dimensional skew-symmetric matrix free to choose, which adds another

 darb=12(M−MS)(M−MS−1) (42)

degrees of freedom.

This sums up to

 do =doS+df+darb=12(M−1)(M−2) (43) de =do−1 (44)

for odd and for even respectively. Interestingly for every the number of free parameters of the generating matrix depends on the parity of only. Matrices with an even have one degree of freedom less.

If we now choose a rate-asymmetry matrix randomly, the probability to hit a matrix with a specific number of selected states is proportional to its size in parameter space. The diagonal transformations (34) do not favour any number of selected states, as they always contribute degrees of freedom. For an even number of selected states there is one free parameter less than for an odd number , such that their generating matrices form a set of probability measure zero in parameter space. Therefore all generators of odd selection numbers have equal size in probability space and the even numbers are suppressed.

It is left to discuss the influence of the reordering transformations (31). They allow to distribute the selected states over the states. The number of possible configurations for this is given by the binomial factor

 (MMS). (45)

After normalization, we infer the distribution (30).

### iv.2 Chaotic quantum kicked rotor

We want to compare these results for uncorrelated random rates to chaotic systems.

A paradigm for quantum chaos is the quantum kicked rotor, a one-dimensional rotor governed by the Hamiltonian

 ^H(^φ,^p,t)=^p22+Kcos(^φ)∑n∈Zδ(t−n) (46)

with time-periodic kicks of strength , period and . For , , we can restrict the system to a torus with periodic coordinate and periodic momentum . Note that since the available phase space volume on the torus is , there exist Floquet states on the torus.

These Floquet states are eigenstates of the one-cycle evolution operator

 ^U(1,0)=e−iℏeffKcos(^φ) e−iℏeff^p22 (47)

fulfilling with corresponding quasienergies .

This kicked rotor is coupled to a bath with temperature . We consider the coupling operator

 ^v=sin(^φ)+cos(^φ), (48)

which respects the periodicity of and breaks the parity (such that also even and odd Floquet states are coupled to each other).

In Fig. 6(b) we show the distribution for the number of selected states that we obtain when randomly choosing the kicking strength within the interval (where the classical counterpart of the quantum kicked rotor is essentially fully chaotic) for a rotor with Floquet states. From Fig. 6(b) and (c) it is clear that the random rate model fails to predict the number of selected states for a typical realization of the chaotic quantum kicked rotor. The distribution is centered around a much larger value than expected for random rates [Fig. 6(a)]. It also seems that the distribution may not be fitted with a binomial distribution, which for example for (black crosses) is much broader then the one that is observed. This trend also manifests itself in Fig. 6(c) where the triangles show the mean number of Bose-selected states for the quantum kicked rotor. This number lies well above , for systems of size about 80% of the states are Bose selected. Consequently, we come to the intriguing conclusion that there must be additional correlations among the rates that are responsible for the fact that significantly more states are selected for the chaotic quantum kicked rotor than in the random-rate model. In the Appendix, we describe a model with correlated random rates that shows a similar distribution of selected states as the quantum kicked rotor model. However, the origin of the large number of selected states for our quantum kicked rotor model remains an interesting open question.

## V Rates with product structure

Consider an arbitrary time-independent system with Hamiltonian which is coupled to a positive temperature bath () and a population-inverted bath () described by a negative temperature through coupling operators which obey the form of a projector on a single quantum state,

 ^v(B1)=|f⟩⟨f|,^v(B2)=|g⟩⟨g|. (49)

Note that the states and can also be coherent superpositions of the single-particle eigenstates . In this section we show that the number of selected states remains always smaller or equal than three, , for all system sizes . Note that in case of this specific system-bath coupling even a system with chaotic single–particle dynamics features only a maximum number of three condensates. An example of a system with such a system bath coupling is the one depicted in Fig. 1(b). Here the states and correspond to the local Wannier orbitals at lattice sites and , respectively.

For coupling operators of the form (49), we find the rate asymmetry matrix

 Aij =2πℏJ(εj−εi)(fifj−gigj) (50)

from Eq. (8), where , and . Here we assume that , giving rise to rates . Moreover, let us first consider ohmic baths with spectral density . We find a rate asymmetry matrix having the product struture

 Aij ∝fifjεj−fifjεi−gigjεj+gigjεi. (51)

Now let be the solution of Eq. (11). It then follows that in the subspace of selected states one has

 0=Aη=⎛⎜ ⎜ ⎜ ⎜⎝fi1fi2...fiMS⎞⎟ ⎟ ⎟ ⎟⎠c1−⎛⎜ ⎜ ⎜ ⎜⎝fi1εi1fi2εi2...fiMSεiMS⎞⎟ ⎟ ⎟ ⎟⎠c2−⎛⎜ ⎜ ⎜ ⎜⎝gi1gi2...giMS⎞⎟ ⎟ ⎟ ⎟⎠c3+⎛⎜ ⎜ ⎜ ⎜⎝gi1εi1gi2εi2...giMSεiMS⎞⎟ ⎟ ⎟ ⎟⎠c4 (52)

with

 c1=∑i∈Sεifiηi,c2=∑i∈Sfiηi,c3=∑i∈Sεigiηi,c4=∑i∈Sgiηi. (53)

Since , and (otherwise we can always shift all by some constant), these coefficients are positive,

 ci>0. (54)

Now if or then the four vectors in Equation (52) will be linearly dependent, thus the can be positive as they should. However, if was greater than three, then generally the four vectors will be linearly independent, so that must hold, in contradiction to the assumption .

Therefore we have shown that for two ohmic baths with product coupling to the system, the number of Bose-selected states is restricted to a maximum of three.

There is a straight-forward generalization of this simple algebraic argument to the case where ohmic baths are coupled to an autonomous system via coupling operators of the product form, Eq. (49). In this case, we find an analog to Eq. (52), but with a linear combination of vectors. By similar reasoning the number of selected states is then restricted to .

We expect that there is a generalization of the above arguement also to non-ohmic systems described by arbitrary spectral densities. We have checked functions of the form (remember that must be odd)

 J(ε)∝|ε|dsgn(ε) (55)

with some power (not necessarily integer) and similarly observe .

Fig. 7(a) confirms the result for ohmic spectral densities . It shows the maximum number of selected states for 200 systems that are randomly drawn from the gaussian orthogonal ensemble, GOE, the ensemble of orthorgonal matrices, where the probability to find a matrix is given by . In random matrix theory these Hamiltonians serve as a model for a fully chaotic system with time-reversal symmetry. We choose random Hamiltonians to make sure that the number of Bose-selected states is not additionally restricted by the system dynamics. We couple these systems to ohmic baths, where we randomly choose an index to which the bath is coupled with operator . For half of the baths we choose positive temperature, for the other half negative temperature (the selected states are only determined by the rate asymmetry matrix and thus independent of the absolute values of these temperatures).

In Fig. 7(a) we plot the maximum number of selected states found for an ensemble of realizations of this model. For sufficiently small it equals the predicted upper bound for . However, we observe that for larger values of , the observed maximum number of selected states saturates very quickly at values that are of the order of where is the system size. This bounds the mean number of selected states as shown in Fig. 7(b).

We would like to stress that the behavior of the autonomous GOE is very different from that of the time-periodically driven rotor, although both systems exhibit chaotic single particle dynamics. First, in the limit of large system size, , the number of selected states can be intensive, as long as the number of baths is not scaled with system size, so that we may find fragmented condensation in the thermodynamic limit. However, even if we scale with system size, as shown in Fig. 7(c) for , the mean number of selected states seems to scale strongly sublinear with at maximum , rather than the drastically different extensive scaling that is observed for time-periodically driven systems, cf. Sec. IV.

## Vi Conclusion and outlook

In this paper, we have addressed the effect of Bose selection in nonequilibrium steady states of driven-dissipative ideal Bose gases that exchange energy with an environment they are weakly coupled to. Namely, when the total particle number is increased in such a system, eventually only the occupation of a number of selected modes will increase linearly with , while the occupation of all other modes saturates. So far, only very little was known about the factors that determine how many modes will be selected.

Within this paper, we have investigated the question of how many Bose-selected states will be found by addressing different relevant scenarios leading to different forms of the rate matrix describing the driven-dissipative ideal Bose gas. We have shown that upper bounds for the number of Bose-selected states that are independent of the dimensionality of the underlying single-particle state space appear both for systems whose rates can be understood as the discretization of a continuous function as well as for systems that couple to baths via a few single-particle quantum states only. In these cases Bose selection corresponds to (fragmented) Bose condensation, since each Bose-selected state acquires a macroscopic occupation in the thermodynamic limit at fixed density .

Moreover, we have discussed two scenarios where the number of selected states is found to grow linearly with , so that none of the selected states will acquire a macroscopic occupation in the thermodynamics limit. On the one hand, we have shown that for randomly drawn uncorrelated rates the number of selected states follows a binomial distribution so that on average half of the single-particle states are selected. On the other hand, we have numerically treated a periodically driven system, the quantum kicked rotor, in a regime where the corresponding classical model is fully chaotic. Averaging over various kicking strengths we find that more than half (about 75 percent) of the states become selected. This is an intriguing result. It implies that the chaotic driven system must give rise to correlations among the rates, which are responsible for a significant enhancement of the selected states. While we were able to construct a model of correlated random rates showing similar behavior, the question about the origin and the nature of the large number of Bose-selected states for the quantum kicked rotor remains open.

Further open questions to be addressed in future research include the role of interactions and larger system bath coupling on the effect of Bose selection, an understanding of the factors determining the number of Bose-selected states in photonic systems where Bose selection is induced by the interplay of pumping, particle loss, and thermalization Vorberg et al. (2018), and the investigation of possible experimental platforms.

###### Acknowledgements.
We acknowledge discussions with Daniel Vorberg. This work was supported by the German Research Foundation DFG via the Research Unit FOR 2414.

*

## Appendix A Modified random wave model for the rates in the quantum kicked rotor

The deviations from the random rate model indicate that the rates of the chaotic quantum kicked rotor contain correlations that lead to the Bose selection of more states than in the uncorrelated case. In this appendix we construct a random rate matrix having correlations that lead to similar behavior.

Correlations of the rates associated with the quantum-kicked-rotor model in contact with a heat bath that we considered were discussed in Ref. Wustmann (2010) (where the single-particle problem was studied). The rate and the backwards rate were proposed to obey

 Rji=(1+ξij)Rij, for i>j, (56)

with both and stemming from individual exponential distributions with scale parameters and , where decreases with system size as . However, the rate model (56) leads to steady states with even less then half of the states being Bose selected (data not shown).

Note that in Sec. III.2 we have encountered a model, which also features that typically more than half of the states are Bose selected. It is the random wave model, Eq. (29), with parameter [cf. also Fig. 5]. A suggestive point of view for why the random wave model might be suitable to approximate rates of a chaotic map, is that random waves have been used successfully to model typical chaotic eigenstates Berry (1977); Bäcker et al. (2010); Berry (2002). Note, however, that this vague reasoning is not based on a microscopic picture for the derivation of the rates.

We can see in Fig. 8(a) that the distribution for the random wave model with is much broader than the one we observe for the quantum kicked rotor in Fig. 6(b). Also the rates that result from a random wave model, Eq. (29), do not follow an exponential distribution. Their distribution is rather peaked at some finite value. To correct for this, we introduce a modified random wave model for the rates

 R(k,q)=L∑l=1 |cl|Re{1+ei(κk,lk+κq,lq+αl)}e−λ|k−q|, (57)

where we constrain the waves to positive values and localize them on a length by introducing an exponential factor. In this model we choose , similar to the random wave model with fixed absolute value of the wavenumber , uniformly distributed angles and from a normal distribution.

We observe that the modified random-wave model has rates that are distributed exponentially and the distribution of the shows relatively good agreement with the data from the quantum kicked rotor for and localization parameter . This can be seen for example by comparing the distribution for discrete states in Fig. 8(b) to the distribution of the quantum kicked rotor in Fig. 6(b), although the mean values coincide, the distribution of our model is, however, a bit broader than the one obtained for the quantum kicked rotor. Also as a function of system size, Fig. 8(c), the model seems to reproduce the mean value of selected states in Fig. 6(c) quite nicely.

However, despite the fact that the modified random wave model (57) gives rise to a similar distribution of the number of selected states as the one obtained for the quantum kicked rotor, we have no evidence that the modified random-wave model mimics the physics of the quantum kicked rotor coupled to a heat bath [see Sec. IV.2].

## References

• Basko et al. (2006) D. M. Basko, I. L. Aleiner,  and B. L. Altshuler, Ann. Phys. 321, 1126 (2006).
• Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Ann. Rev. Cond. Mat. Phys. 6, 15 (2015).
• Schreiber et al. (2015) M. Schreiber, S. S. Hodgman, P. Bordia, H. P. LÃ¼schen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider,  and I. Bloch, Science 349, 842 (2015).
• Smith et al. (2016) J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse