A solvable model for strongly interacting nonequilibrium excitons

TL;DR

This paper introduces a solvable model to study nonequilibrium phase transitions in excitons formed in semiconductor bilayers under intense laser excitation.

Contribution

The paper presents a novel solvable model for nonequilibrium excitons, capturing phase transitions in open quantum systems.

Findings

The model shows a nonequilibrium phase transition out of the Mott insulator phase of excitons.

The steady-state density matrix can be expressed in closed-form under certain conditions.

Phase transitions in the model differ from thermal systems due to non-thermal steady states.

Abstract

Excitons, bound states of electrons and holes, can be formed by shining a laser on a semiconductor. Recent experiments have suggested that sufficiently intense laser excitation applied to a WS2-WSe2 bilayer creates a “Mott insulator” of excitons in which they fill a lattice and can no longer easily move around. Since the excitons have a finite lifetime before decaying, this phase transition is fundamentally “nonequilibrium.” By developing a solvable toy model that goes beyond the standard approach of the literature, we theoretically study and characterize a nonequilibrium phase transition out of the Mott insulator phase and validate some of the experimental heuristics. The results provide a relevant exemplar of a phase transition in an open quantum system.

We study the driven-dissipative Bose-Hubbard model with an all-to-all hopping term in the system Hamiltonian, while subject to incoherent pumping and decay from the environment. This system is naturally probed in several recent experiments on excitons in WS2/WSe2 moiré systems, as well as quantum simulators. By positing a particular form of coupling to the environment, we derive the Lindblad jump operators and show that, in certain limits, the system admits a closed-form expression for the steady-state density matrix. Away from the exactly solvable regions, the steady state can be obtained numerically for 100s to 1,000s of sites. We study the nonequilibrium phase diagram and phase transitions, which qualitatively matches the equilibrium phase diagram, agreeing with the intuition that increasing the intensity of the light is equivalent to changing the bosonic chemical potential. However, the steady states are far from thermal states, and the nature of the phase transitions is changed.

Figures

We perform scaling collapses for the order parameters at the (A) transition between the nB=1 and nB=2 Mott phases and (B) transition between Mott and superfluid phases. As N→∞, the data in (A) [(B)] collapses to the functions fnB,m(x) [fΨ2,s(x)] defined in the text. We find that λm=1 and βm=0 for the former and λs=2=βs+1 for the latter. We have confirmed that this exponent stays the same regardless of which Mott lobes are involved and regardless of whether we tune I0/γ0 or t/U.

Equations

Full text

The Bose–Hubbard model is the paradigmatic model of interacting bosons on a lattice and describes a transition between a Mott insulating and superfluid phase (1). There are now several platforms for experimentally realizing this model such as trapped bosonic atoms (2???–6) or photon modes (7?–9). In the latter case, interactions are induced through circuit quantum electrodynamics (QED) (10, 11), cavity QED (12??–15), semiconductor platforms (16?–18), or Rydberg atoms (19, 20), and a lattice can be produced by interconnecting these platforms (15, 18, 21). Excitons in transition-metal dichalcogenide (TMD) moiré bilayers (22?–24) or artificial lattices (25) provide another emerging solid-state realization of Bose–Hubbard physics, where the interactions are naturally implemented through strong dipole repulsions. All these experimental settings move beyond the study of equilibrium physics as these systems either controllably couple to their environments or have finite lifetimes for the bosons (10, 26).

We are most interested in understanding several recent experiments potentially realizing bosonic Mott physics with excitons in WSe_2_/WS_2_ moiré bilayers (27???–31). The experiments see a finite energy shift in the photoluminescence spectrum interpreted as an on-site exciton–exciton interaction with multiple excitons per site. There is considerable evidence that the excitons enter an interaction-dominated Mott-insulating regime at sufficiently large intensities of light (27, 30). The inherently nonequilibrium nature of these experiments, given the necessity of continuously pumping the system full of decaying excitons, naturally raises the question of how much these reported nonequilibrium Mott phases resemble their equilibrium counterpart, and whether an analog of the superfluid phase exists out of equilibrium.

The Bose–Hubbard model,

$$[eqn]$$

subject to dissipation and driving, is the prototypical model for studying strongly interacting open quantum systems. Here, $[eqn]$ denotes the bosonic annihilation (creation) operator at site $[eqn]$ , and $[eqn]$ is the boson number operator there. The parameters $[eqn]$ , $[eqn]$ , and $[eqn]$ are the on-site Coulomb repulsion, hopping amplitude between sites $[eqn]$ and $[eqn]$ , and chemical potential, respectively. For moiré excitons, $[eqn]$ is simply the excitation gap of a single exciton in the minimum of the moiré potential.

The drive is often coherently added to the total Hamiltonian, $[eqn]$ , with $[eqn]$ being the Rabi frequency associated to a coherent driving field, e.g. electromagnetic wave mode of a cavity, and the decay is implemented through a Lindblad master equation for the density matrix $[eqn]$ (32?–34),

$$[eqn]$$

where $[eqn]$ , $[eqn]$ are the jump operators and $[eqn]$ are the associated damping rates. By including suitable Lindblad jump operators, we can include incoherent driving as the reverse process of decay. In our work, we assume all driving is incoherent and set $[eqn]$ .

This system has been numerically studied with numerous methods: mean-field theory (MFT) (35?????–41), cluster MFT (42, 43), the positive $[eqn]$ representation (37, 44), the corner-space renormalization group method (45), functional renormalization group using a Keldysh path integral (32, 46), quantum trajectories (38, 43, 47), matrix-product-state-based methods (48?–50), and a truncated hierarchy of correlators (51). These studies all assume a simple and decoupled form of Lindblad jump operators, i.e. $[eqn]$ , which works well for describing nearly isolated noninteracting atoms. In general, however, such an assumption does not hold for a strongly correlated system when $[eqn]$ .

In this work, we emphasize that the explicit form of jump operators is determined by the system Hamiltonian itself, as well as the coupling with the environment, which can significantly modify the Lindblad master equation and the properties of the steady state. In order to fully understand the phases and dynamics in these nonequilibrium systems, it is useful to have prototypical solvable models (52). Here, we study the Bose Hubbard Hamiltonian Eq. 1 with an all-to-all hopping for $[eqn]$ sites, i.e. $[eqn]$ for any $[eqn]$ , which may be considered a sort of infinite-dimensional limit. The all-to-all hopping mimics the effect of a mean-field analysis at the equilibrium level and allows us to move beyond mean-field with analytical control. Remarkably, in such an approach we are able to derive the explicit form of the Lindblad jump operators, analytically obtain closed-form solutions of the steady state in certain limits, and numerically solve the system for up to thousands of sites in others.

Our main results are summarized in the phase diagram, Fig. 1, a color map of the superfluid order parameter $[eqn]$ , as a function of the ratio of exciton production and decay rates $[eqn]$ (appearing as damping rates in the Lindblad equation Eq. 2) and $[eqn]$ . $[eqn]$ and $[eqn]$ are pumping and decay rates at energy $[eqn]$ , the detailed structure of which are specified in Section 4. Similar to the equilibrium model, we see that there are two phases: Mott-insulator phases with an integer values of $[eqn]$ and a superfluid phase identified as a nonzero value of $[eqn]$ with a generic value of $[eqn]$ . Although the steady-state density matrices are not of a (thermal) equilibrium form, the nonequilibrium and equilibrium phase diagrams resemble each other lending credence to the qualitative understanding of the experiment as tuning the boson chemical potential. We additionally study the critical properties of the phase transition. Notably, the Mott–Mott phase transition at $[eqn]$ (either from $[eqn]$ to $[eqn]$ , or $[eqn]$ to $[eqn]$ ), which is simply a level crossing at the equilibrium level, becomes a continuous phase transition and the critical exponents for the Mott–superfluid transition are distinct from their equilibrium values.

Derivation of the Master Equation

We will analyze the Bose–Hubbard Hamiltonian, Eq. 1 with all-to-all hopping between the $[eqn]$ sites, $[eqn]$ , and in the presence of coupling to the environment. We will assume that $[eqn]$ as is relevant for the physical system.

We now proceed to derive the Lindblad master equation, Eq. 2, in the weak-coupling limit for this model. We follow the standard procedure (53, 54): First, we start with a general coupling form between the system and environment $[eqn]$ , where $[eqn]$ ( $[eqn]$ ) are Hermitian operators that act on the system (environment), and the index $[eqn]$ labels different operators. Then, we decompose $[eqn]$ by eigenoperators of the system Hamiltonian $[eqn]$ , which are defined by the relation $[eqn]$ . The Lindblad jump operators are nothing but the eigenoperators $[eqn]$ , and the damping rate $[eqn]$ is given by $[eqn]$ . With these operators in hand, the equation of motion for the system’s density matrix can be written out explicitly according to Eq. 2.

Since the excitons emit light when they decay, we know that they interact with the quantized electromagnetic field. A typical dipole interaction takes the form $[eqn]$ , with $[eqn]$ . In momentum space, after the rotating-wave approximation, it can be expressed as (55)

$$[eqn]$$

where $[eqn]$ is related to the polarization of the light mode. We used $[eqn]$ since the relevant wavelength of light ( $[eqn]$ nm) is much larger than the distance between moiré unit cells ( $[eqn]$ 7.5 nm) (56). Therefore, as opposed to the common assumption where each site couples to different modes in the bath, we assume the environment couples with the system in a collective way, reminiscent of the Dicke model (57). Note that $[eqn]$ where $[eqn]$ is the volume available to the photon modes, and we suppress the light’s polarization for simplicity. With such an assumption, we may write the coupling as

$$[eqn]$$

with

$$[eqn]$$

where $[eqn]$ denotes the environment operator. Though we only argued above for the decay of excitons to have such a form, treating the production and decay on equal footing allows us to make theoretical progress.

We do not specify the explicit forms of $[eqn]$ since they could be complicated in general. Specifically, the long-lived interlayer excitons which form the exciton lattice are generated through at least a two-step process, where intralayer excitons are photoexcited and then converted to interlayer excitons through some relaxation process (23). We therefore consider only the incoherent contributions of $[eqn]$ to the exciton system.

We now decompose $[eqn]$ in $[eqn]$ in terms of eigenoperators, which, in general, requires diagonalizing $[eqn]$ . Note that the coupling to the environment and the Hamiltonian itself is invariant under all permutations of the sites; we also assume the system’s initial state (before turning on exciton-producing effects) is the vacuum $[eqn]$ since $[eqn]$ , which has permutation symmetry as well. Therefore, the only states that are accessed are in the fully symmetric sector of this permutation symmetry and it is sufficient to diagonalize $[eqn]$ in this sector.

Similar to the procedure adopted in multilevel systems (58), we parameterize the sector with the normalized states

$$[eqn]$$

where $[eqn]$ is a vector in $[eqn]$ dimensional space, $[eqn]$ is the symmetric group, $[eqn]$ is the permutation operator that shuffles site $[eqn]$ to $[eqn]$ , and the state $[eqn]$ has the first $[eqn]$ sites being unoccupied, the next $[eqn]$ sites having one boson occupied, etc. We have chosen a maximum number of allowed bosons per site, $[eqn]$ , and $[eqn]$ .

It is convenient to define transition operators $[eqn]$ which changes site $[eqn]$ from having $[eqn]$ bosons to $[eqn]$ bosons and acts as the identity on all other sites. Note that $[eqn]$ for all $[eqn]$ . We can then find

$$[eqn]$$

where $[eqn]$ are the unit vectors and $[eqn]$ are the components of $[eqn]$ .

Letting $[eqn]$ be the total boson number operator with eigenvalues $[eqn]$ , we use these relations to write $[eqn]$ in the $[eqn]$ basis:

$$[eqn]$$

where $[eqn]$ . For a given $[eqn]$ and $[eqn]$ , we can numerically find the eigenbasis within each $[eqn]$ sector where only the ratio $[eqn]$ determines the eigenvectors. We denote the eigenvectors $[eqn]$ where each sector has a different number of allowed $[eqn]$ .

We now eigen-decompose the coupling operator $[eqn]$ within the fully permutation-symmetric sector to find

$$[eqn]$$

where $[eqn]$ and $[eqn]$ is a combined index. Note that the operators in the sum have energy $[eqn]$ .

With the eigenoperators in hand, we can immediately write out the interaction-picture Lindblad equation (ignoring the Lamb shift)

$$[eqn]$$

The first (second) line corresponds to decay (pumping) processes. Equality between $[eqn]$ and $[eqn]$ is determined by the condition $[eqn]$ , where $[eqn]$ sets the relaxation timescale of the system. The specifics of the environment enter only through

$$[eqn]$$

and we have assumed $[eqn]$ (and $[eqn]$ is set by $[eqn]$ which are assumed to be the smallest energy scales in the weak-coupling limit).

We now write out Eq. 8 in the eigenbasis, $[eqn]$ . One difficulty in solving the resulting equation is the presence of coherences (i.e. off-diagonal terms of the density matrix in the basis of energy eigenstates). As we derive in SI Appendix, the Lindblad master equation can only generate coherences within energy-degenerate sectors, and therefore, if there is no degeneracy, the steady-state density matrix is diagonal. However, there is a weaker condition that ensures that the $[eqn]$ operators never generate coherences when acting on a diagonal $[eqn]$ (see SI Appendix): $[eqn]$ ( $[eqn]$ ) only when $[eqn]$ or $[eqn]$ ( $[eqn]$ or $[eqn]$ ), respectively. In our analytic or numerical results, we always check that the above condition is satisfied for this system.

We can therefore write $[eqn]$ and

$$[eqn]$$

As a check on our derivation, we can set $[eqn]$ and $[eqn]$ to represent a thermal bath of free bosons. In this case, the Kubo-Martin-Schwinger condition holds allowing us to write $[eqn]$ and $[eqn]$ where $[eqn]$ is a function only of $[eqn]$ and $[eqn]$ is the Bose–Einstein distribution. The function $[eqn]$ is determined by the specifics of the environment [e.g. $[eqn]$ for the photon field (54)]. We numerically find that the steady-state density matrix, $[eqn]$ has the thermal form $[eqn]$ . In general, we consider nonzero $[eqn]$ as resulting from light-induced, rather than thermal, production of excitons.

Key Observables, Phases, and Exponents

Once the steady state satisfying $[eqn]$ from Eq. 10 is determined, there are two key observables to evaluate

$$[eqn]$$

The two phases present in our model are Mott-insulating phases with integer exciton density $[eqn]$ (and $[eqn]$ ) and a superfluid phase with generic $[eqn]$ and $[eqn]$ . In the equilibrium case, the superfluid phase is determined by the nonzero expectation $[eqn]$ (or equivalently the long-range order of $[eqn]$ ), and $[eqn]$ is our analogous observable.

In addition to the phases themselves, we study the phase transitions and extract the critical exponents. We will tune across the phase transitions primarily by adjusting the ratio $[eqn]$ of the exciton production and decay strengths. With $[eqn]$ , we may Taylor expand $[eqn]$ around $[eqn]$ and the transition occurs when $[eqn]$ reaches a critical value $[eqn]$ .

At a continuous (second-order) phase transition, there is an order parameter $[eqn]$ that takes a particular form close to the phase transition

$$[eqn]$$

where $[eqn]$ is a constant, and the correlation length $[eqn]$ of the system will diverge as $[eqn]$ . The exponents $[eqn]$ and $[eqn]$ are universal numbers set by the phase transition itself, independent of the exact microscopic model. For systems of a finite length $[eqn]$ , this diverging correlation length is cut off and will lead to corrections of the critical behavior that depend on $[eqn]$ . Replacing the correlation length by $[eqn]$ leads to the finite-size scaling form for $[eqn]$ , where $[eqn]$ is a finite-size scaling function (see e.g. ref. 59).

The all-to-all hopping in our system makes it difficult to define a correlation length, so we define the exponent $[eqn]$ as the analog for $[eqn]$ in our system. At the Mott–Mott transition from the $[eqn]$ to $[eqn]$ Mott phases, the order parameter is simply $[eqn]$ , which therefore has the form

$$[eqn]$$

where $[eqn]$ , $[eqn]$ are the critical exponents and $[eqn]$ is the scaling function for $[eqn]$ at the Mott–Mott transition. Similarly, at the Mott–superfluid (Mott–Mott) transition, the order parameter is $[eqn]$ and takes the scaling form

$$[eqn]$$

with similarly defined $[eqn]$ , $[eqn]$ and $[eqn]$ . Note that the exponents $[eqn]$ and $[eqn]$ and scaling functions do not depend on the exact definition of the tuning parameter used to move between the two phases, and they can distinguish between different critical points.

Analytic Results

Due to needing to find the eigenvectors of $[eqn]$ , Eq. 7, we cannot make additional analytic progress in most cases. We will treat the general problem numerically in the next section, but there are two cases where more analytic progress can be made.*

t = 0.

3.1.

In this limit, the eigenvectors of $[eqn]$ are simply the $[eqn]$ states defined above. Instead of simplifying Eq. 10), we note that

$$[eqn]$$

and $[eqn]$ . We can therefore write out the equivalent of Eq. 8

$$[eqn]$$

We now substitute $[eqn]$ . Using Eq. 6 and the fact that we start from a vacuum state, one can check that the condition $[eqn]$ when $[eqn]$ , is preserved under time evolution. Therefore, we may instead write $[eqn]$ .

We can explicitly write out the steady-state equation as

$$[eqn]$$

where the above equation holds for all $[eqn]$ and $[eqn]$ denote components of $[eqn]$ . A sufficient “detailed-balance”-like condition for the steady state is that all terms in the sum are zero,

$$[eqn]$$

which remarkably specifies the solution in this case. (In general, the equivalent set of equations will not hold as they overdetermine the solution.) It can easily be checked that the steady state is given explicitly as

$$[eqn]$$

where we have defined $[eqn]$ . This form of the steady state manifestly reveals that we can set the maximum number of bosons per site we need to consider, $[eqn]$ , by the condition $[eqn]$ , when $[eqn]$ .

We now analyze the case $[eqn]$ in detail, which can be mapped to the case of a three-level system with collective decay and pumping, and thus can be potentially implemented in cavity experiments. We can rewrite the density matrix as $[eqn]$ (and we assume that $[eqn]$ ). We can compute the average number of bosons per site $[eqn]$ as

$$[eqn]$$

We can see, in the thermodynamic limit, that our model has a sharp transition from an $[eqn]$ to $[eqn]$ to $[eqn]$ state as we increase the ratio $[eqn]$ and $[eqn]$ . We identify these regions as Mott insulators (and we can numerically check that $[eqn]$ ). We can analyze these critical points by evaluating

$$[eqn]$$

implying that $[eqn]$ has the scaling form $[eqn]$ with $[eqn]$ and $[eqn]$ . We expect that the values of $[eqn]$ are insensitive to increasing $[eqn]$ .

In equilibrium, this transition is first order since it is simply a level crossing of the two competing ground states, $[eqn]$ to $[eqn]$ . However, in the nonequilibrium setting, as we approach the phase transition there are more and more states in the Hilbert space that contribute to the steady-state density matrix, which allows for a continuous phase transition to emerge between $[eqn]$ to $[eqn]$ .

M = 1.

3.2.

If we consider the case of hardcore bosons, i.e. $[eqn]$ , we can make additional progress as well because there is only one state per (fully symmetrized) sector, $[eqn]$ , and therefore the eigenstates of $[eqn]$ can be specified as $[eqn]$ . The eigenenergy of the states is given by

$$[eqn]$$

and we can evaluate

$$[eqn]$$

so $[eqn]$ and these have energies $[eqn]$ . When $[eqn]$ is nonzero, we see that the energy will take on a range of values from $[eqn]$ in the thermodynamic limit.

Since all $[eqn]$ are unique, we still have no coherences and can write down the density matrix equation

$$[eqn]$$

where $[eqn]$ is specifying $[eqn]$ . Again noting that a sufficient condition for the steady-state solution is

$$[eqn]$$

we find

$$[eqn]$$

If $[eqn]$ , we see that $[eqn]$ for all $[eqn]$ , so we quickly derive Eq. 19 for $[eqn]$ .

We will now show that, in addition to the $[eqn]$ and $[eqn]$ Mott phases seen before, there is an intermediate phase. We first assume the pump-to-decay ratio $[eqn]$ decreases monotonically with $[eqn]$ , which fits the intuition that higher energy excitations are prone to decay more quickly. Defining $[eqn]$ [ $[eqn]$ ], we can recover the Mott phases from before by computing $[eqn]$ $[eqn]$ when $[eqn]$ ( $[eqn]$ ), respectively, where $[eqn]$ are the corresponding values of $[eqn]$ (see SI Appendix for additional calculation details for this section).

Next, we consider the case $[eqn]$ . There then exists some index $[eqn]$ where $[eqn]$ . From Eqs. 25 and 26, we can find as $[eqn]$

$$[eqn]$$

where $[eqn]$ from the monotonic decrease assumption.

We can deduce this density matrix represents an additional phase by measuring $[eqn]$ . Since $[eqn]$ , we see that $[eqn]$ . In the three aforementioned limits, we find

$$[eqn]$$

As above, we can observe a critical scaling emerge as $[eqn]$ . Using Eq. 26, we can find

$$[eqn]$$

which has the scaling form $[eqn]$ , implying that $[eqn]$ and $[eqn]$ . A similar finite-size scaling form, with the same value of $[eqn]$ , can be derived for $[eqn]$ .

We note that the $[eqn]$ hardcore boson case can actually be mapped to an Ising model, with corresponding interactions and dissipations. Various cases have been studied (60?–62), but with a single jump operator which corresponds to cavity loss or cavity pumping. When the Lindblad jump operators are derived from the interacting Hamiltonian, the physics becomes more distinct.

Numerical Results

From the analytic results, we can see there are two phases distinguished by $[eqn]$ . If $[eqn]$ , then we are in a Mott insulating phase where $[eqn]$ is approximately an integer. Otherwise, $[eqn]$ takes an intermediate value, and $[eqn]$ is nonzero, like in the superfluid phase. By solving the system numerically, we can explore what happens when $[eqn]$ and more than one boson per site is allowed.

In order to proceed numerically, we must specify $[eqn]$ and $[eqn]$ . We take $[eqn]$ , the correct form for photon-driven decay, and we assume $[eqn]$ for simplicity.† For the moiré TMD system (27), the realistic excitation gap $[eqn]$ eV and the repulsion $[eqn]$ meV. However, we fix $[eqn]$ meV as $[eqn]$ only affects $[eqn]$ [and consequently $[eqn]$ ] which we need to be large to achieve numerical convergence. We let $[eqn]$ set the energy scale, and vary both $[eqn]$ and $[eqn]$ to map out the phase diagram and study the critical exponents.

Once $[eqn]$ and $[eqn]$ , we can only numerically determine the Lindblad operators and energy eigenstates, and the equivalent conditions to Eqs. 18 and 25 cannot all be satisfied. Instead, we construct the Liouvillian superoperator $[eqn]$ satisfying $[eqn]$ from Eq. 10 and search for the steady-state $[eqn]$ that has eigenvalue zero. We can numerically evaluate the $[eqn]$ case up to a large system size of $[eqn]$ . In this case, when $[eqn]$ , the various $[eqn]$ sectors are nondegenerate implying that the steady-state density matrix is diagonal, as we argued above. We check numerically that this nondegeneracy persists when $[eqn]$ .

In addition to measuring $[eqn]$ and $[eqn]$ , we measure the entropy $[eqn]$ and the Liouvillian gap, $[eqn]$ determined from the smallest real part of any of the nonzero eigenvalues of $[eqn]$ (63). These two quantities allow us to better distinguish between the two phases and characterize their critical behavior numerically. We extrapolate these quantities to the $[eqn]$ limit (SI Appendix).

In Fig. 1, we observe the same two phases as above: The Mott phase is signified by integer $[eqn]$ with $[eqn]$ and $[eqn]$ , and the superfluid phase is signified by generic $[eqn]$ and $[eqn]$ . At the phase boundary, $[eqn]$ attains a maximum, clearly demarcating the two phases. The phase diagram in Fig. 1 qualitatively matches the equilibrium phase diagram with its $[eqn]$ Mott lobe. If $[eqn]$ , then we would see additional Mott lobes as well (SI Appendix).

In addition to mapping out the phase diagram, we can confirm the universal scaling forms discussed above. In Fig. 2A, at the transition between two Mott phases, we observe that $[eqn]$ , as expected from Eq. 13, indicating that $[eqn]$ . Within the equilibrium phase diagram, this transition is first order (a level crossing) indicating that the environment can change the nature of the critical point. In Fig. 2B, we similarly see $[eqn]$ , as above, for a transition between the Mott and superfluid phase with $[eqn]$ and $[eqn]$ , as opposed to the equilibrium values of $[eqn]$ (SI Appendix). We confirm that the same exponents are extracted while tuning $[eqn]$ instead of $[eqn]$ across the phase transition.

We perform scaling collapses for the order parameters at the (A) transition between the nB=1 and nB=2 Mott phases and (B) transition between Mott and superfluid phases. As N→∞, the data in (A) [(B)] collapses to the functions fnB,m(x) [fΨ2,s(x)] defined in the text. We find that λm=1 and βm=0 for the former and λs=2=βs+1 for the latter. We have confirmed that this exponent stays the same regardless of which Mott lobes are involved and regardless of whether we tune I0/γ0 or t/U.

Finally, when $[eqn]$ , we can extract the Liouvillian gap, $[eqn]$ . Away from the critical point, $[eqn]$ approaches a constant indicating an effect similar to superradiance (54), which is balanced by superabsorption (see below). At criticality, however, we find $[eqn]$ , indicating a critical slowing down (i.e. a longer time to reach the steady state) near the phase boundary (SI Appendix, Fig. S7). We expect that this same behavior should arise at the phase transitions when $[eqn]$ , but we are unable to compute $[eqn]$ for large enough system sizes.

Discussion and Conclusion

We have introduced a solvable variant of the Bose Hubbard model coupled to an environment. As opposed to the zero-temperature equilibrium case where the ground state determines the density matrix, for the nonequilibrium system we focus on the steady-state density matrix. Due to the lack of energy conservation for the open system, the steady state is instead determined by a balance between the radiative decay and (incoherent) laser pumping. We find two phases that resemble the equilibrium Mott and superfluid phases, but the density matrix does not take a Boltzmann form (see e.g. Eq. 27). We extract critical exponents $[eqn]$ ( $[eqn]$ ) for the transition between a Mott phase and a Mott (superfluid) phase, respectively. These exponents are different from the observed equilibrium exponents derived in SI Appendix (and the Mott–Mott transition is first order in equilibrium). The change in criticality is not surprising—although previous work found that the open system criticality can only change a dynamic exponent (32, 46), universality is typically only insensitive to local perturbations, and the environment acts as a global perturbation.

We derive the Lindblad operators and master equations from an explicit form of the coupling between the system and environment. Our model is solvable because the environment couples in a site-invariant way that allows us to only consider the fully symmetric sector of the permutation group, and the calculation is controlled by $[eqn]$ for large system size. We note that, although the full eigensystem information of $[eqn]$ is difficult to obtain, only the eigenstates that are connected by the derived Lindblad operators are relevant, and we demonstrate this idea for the case $[eqn]$ in SI Appendix. For generic $[eqn]$ , it may be possible to obtain Lindblad operators perturbatively controlled by $[eqn]$ , but we leave this work for the future.

Our formulation of site-invariant coupling between the system and its environment is reminiscent of the Dicke model (57) and will display superradiance (54). However, since we model the production as the reverse process of decay, our model also has superabsorption. This observation explains why, when $[eqn]$ , the steady state has only a finite number of excitons in the thermodynamic limit: if the production rate $[eqn]$ is less than the decay rate $[eqn]$ , an extensive number of excitons cannot be built up.

In the experimental systems (27???–31), there may be an additional local coupling $[eqn]$ responsible for the production of the bosons, but introducing such a term into our formalism would lead to a sum over site-dependent dissipative terms in the Lindblad equation that requires considering sectors outside of the fully symmetric one. Therefore, the formalism in this paper would be insufficient to solve for the steady state for similarly sized systems. Future work should consider the competition between a local production of bosons and their global (or longer-range) decay, as well as a short-range hopping system Hamiltonian, which might be resolved by numerical methods such as dynamical mean-field theory.

Materials and Methods

We implement numerically finding the steady state in Python. We first specify a total number of sites $[eqn]$ , after which we can loop through all possible $[eqn]$ and sort them into their corresponding $[eqn]$ sector using $[eqn]$ . We then specify $[eqn]$ , $[eqn]$ , and $[eqn]$ and construct the Hamiltonian $[eqn]$ in each $[eqn]$ sector using Eq. 7. We diagonalize $[eqn]$ using built-in numpy routines and then compute $[eqn]$ using Eq. 8 and the eigenvectors from each sector.

By specifying $[eqn]$ and $[eqn]$ functions (as in the main text), we then construct a sparse matrix representing the Liouvillian operator, $[eqn]$ , acting on the vector space of all possible density matrices. That is, $[eqn]$ with $[eqn]$ being the right-hand-side of Eq. 10. Using the sparse matrix routines available in scipy, we can find the eigenvalue–eigenvector pairs with eigenvalues close to zero (or we can directly diagonalize the entire matrix if it is small enough). We have then computed the unique $[eqn]$ satisfying $[eqn]$ and the Liouvillian gap $[eqn]$ , as defined in the main text. For $[eqn]$ , we find considerable speed up of the code by only targeting a few eigenvalue–eigenvector pairs which makes the determination of $[eqn]$ difficult. All other quantities we discuss can be computed using $[eqn]$ :

$$[eqn]$$

and

$$[eqn]$$

Supplementary Material

Appendix 01 (PDF)