Nonequilibrium phase transition in a driven-dissipative quantum antiferromagnet

A deeper theoretical understanding of driven-dissipative interacting systems and their nonequilibrium phase transitions is essential both to advance our fundamental physics understanding and to harness technological opportunities arising from optically controlled quantum many-body states. This paper provides a numerical study of dynamical phases and the transitions between them in the nonequilibrium steady state of the prototypical two-dimensional Heisenberg antiferromagnet with drive and dissipation. We demonstrate a nonthermal transition that is characterized by a qualitative change in the magnon distribution from subthermal at low drive to a generalized Bose-Einstein form including a nonvanishing condensate fraction at high drive. A finite-size analysis reveals static and dynamical critical scaling at the transition, with a discontinuous slope of the magnon number versus driving field strength and critical slowing down at the transition point. Implications for experiments on quantum materials and polariton condensates are discussed.

Figure 1

Nonequilibrium phase diagram of driven-dissipative steady states. Steady states as a function of drive strength $g$ and quantum fluctuations, parametrized here by the inverse spin length $1/S$ but controlled in physical systems by many factors, including geometrical frustration. The red section along the vertical axis marks the antiferromagnetically ordered ground state at $T=0$. The black curve separating the ordered (orange) and disordered (green) subthermal phases is obtained by determining the value of $1/S$ at which the staggered magnetization [as defined in Eq. (A10)] vanishes for a given $g$. The gray vertical line at $g=1$ separates the subthermal regime from the superthermal regime, which turns into a thermal distribution plus a $\delta$-function in the interacting system in the thermodynamic limit. The critical end point at $g=1, 1/S=0$ is a specific feature of the Heisenberg antiferromagnet in two dimensions. The gray dashed curve indicates the expected behavior in three dimensions or in the anisotropic $xy$ or $xxz$ (Ising, gapped) regimes in two dimensions.

Figure 2

Nonequilibrium phase transition. (a) Interaction-induced changes in the steady-state magnon occupation $n\left(\omega \right)$. Plotted is $\omega n\left(\omega \right)$ as a function of magnon frequency $\omega$ in order to highlight the difference between subthermal [$\omega n\left(\omega \right)\to 0$ for $\omega \to 0$], thermal [$\omega n\left(\omega \right)\to \text{const}$], and superthermal [$\omega n\left(\omega \right)\to \infty$] regimes. The blue (red) data points show the interacting results for representative subthermal (superthermal), $g=0.5$ ($g=1.5$), cases in comparison with the noninteracting results shown by blue (red) curves. The dark gray solid line indicates a thermal state at $g=1$ and $T=0.6$; the light gray dashed line is a best fit to the high-frequency part of the interacting distribution function at $g=1.5$ and corresponds to a thermal state with an effective temperature $T>\tilde{T}$. Inset: The same results plotted as $n\left(\omega \right)$ versus $\omega$, focusing on the low-frequency part to highlight that the interacting superthermal system shows a low-frequency divergence that is stronger than both the noninteracting system and the best thermal fit. (b) Points: magnon number vs total energy curve defined from Eq. (4) with $g$ as an implicit parameter for both noninteracting and interacting steady states. Solid black line: magnon number vs total energy relation obtained from Bose distribution with chemical potential $\mu =0$ with temperature as an implicit parameter. States below this critical Bose-Einstein condensation line have a lower number of magnons per energy than a thermal state. States above the critical line have a number of magnons that exceeds the maximal number in states with $\omega \left(\mathit{k}\right)>0$ that is compatible with the given system energy in a thermal state, implying the existence of a $\delta$-function contribution at zero energy (condensate fraction) in the thermodynamic limit.

Figure 3

Finite-size scaling analysis revealing the $\delta$-function contribution at $g>1$ in the thermodynamic limit. (a) and (b) Ratio of the magnon density at the lowest frequency and the total number of magnons in the system ${\mathcal{N}}_0/\mathcal{N}$ and (c) and (d) ratio of the magnon density at the second-lowest frequency and the lowest frequency ${\mathcal{N}}_1/{\mathcal{N}}_0$ for $g=1.25$ (left panels) and $g=1.5$ (right panels). Different colors correspond to different values of ${\gamma }_{\text{out}}=0.002,0.02,0.2$ as indicated. Black points correspond to the noninteracting stationary state, gray points show thermal behavior ($g=1$), and red stars correspond to the stationary state to which the interacting, closed system evolves when initialized with the respective noninteracting stationary state at given $g$.

Figure 4

Static and dynamical critical behavior in the interacting driven-dissipative steady state. (a) Rate of change of magnon number $\mathcal{N}$ as a function of $g$. Inset: Scaling behavior with the linear system size collapses the data points onto a single curve. (b) Rate of decay of total magnon number $\mathcal{N}$ towards the stationary state, plotted as a function of $g$ for different system sizes as indicated. Inset: Scaling behavior with the linear system size is consistent with collapse onto a single curve, suggesting critical slowing down as $g\to 1$. (For critical behavior in the strength of the condensate fraction, see Appendix B.)

Figure 5

Condensate fraction ${\mathcal{D}}_0$ as a function of the dimensionless tuning parameter $g$ for different linear system sizes as indicated.

Figure 6

Finite-size scaling analysis analogous to Fig. 4 for the subthermal regime ($g=0.875$). (a) Ratio of the magnon density at the lowest frequency and the total number of magnons in the system ${\mathcal{N}}_0/\mathcal{N}$ and (b) ratio of the magnon density at the second-lowest frequency and the lowest frequency ${\mathcal{N}}_1/{\mathcal{N}}_0$. Different colors correspond to different values of ${\gamma }_{\text{out}}=0.002,0.02,0.2$ as indicated. Black points correspond to the noninteracting stationary state, gray points show thermal behavior ($g=1$), and red stars correspond to the stationary state to which the interacting, closed system evolves when initialized with the respective noninteracting stationary state at given $g$.

Figure 7

Magnetic Brillouin zone (MBZ) for $\ell =8$. The full MBZ (yellow) can be reduced to $({\ell }^2+2\ell )/8$ lattice sites (green) due to the symmetry of the lattice. The multiplicity weights of the reduced lattice vectors that are sufficient to simulate the dynamics in the system are marked as indicated.

Figure 8

Mapping of the momentum grid onto an energy grid for $\ell =8$. The interval $\{0,{\mathrm{\Omega }}_{\max }\}$ into $\ell$ equidistant energy bins (blue and magenta) and for each momentum vector ${\mathit{k}}_{\mathrm{PZ}}$ of the associated bin is determined. The red numbers give the total weight of all vectors within the energy bin, so, for example, energy bin ${\omega }_5$ has momentum vectors $\{{\mathit{k}}_5,{\mathit{k}}_6,{\mathit{k}}_7,{\mathit{k}}_8,{\mathit{k}}_9,{\mathit{k}}_{10}\}$, which have a total weight of 10.