Kibble–Zurek Scaling in the Dicke Model at Mesoscopic Scales

Abstract

The Dicke model is a paradigmatic setting for collective light–matter physics and the superradiant phase transition. Yet extracting the critical exponents is challenging at experimentally accessible mesoscopic sizes, due to the slow divergence of the correlation time under all-to-all coupling and a photon-loss-driven crossover to a distinct dissipative universality class. Here, we perform a large-$N$ analysis that identifies distinct coherent and dissipative fixed points for the closed and open Dicke models. We then develop a unified mesoscopic scaling framework that incorporates the leading irrelevant correction and, going beyond static and spectral probes, brings ramping dynamics under the same scaling description. It recovers the corresponding exponents, verifies Kibble–Zurek scaling, and clarifies how finite size, dissipation, and speed compete in the ramping dynamics. Our work thus establishes a unified framework for resolving static and dynamical critical scaling in closed and open quantum systems, with broader applicability to mesoscopic systems with long-range interactions.

Introduction.— Collective light–matter systems realize long-range interactions among quantum emitters mediated by optical modes 46; 35; 53. The Dicke model captures this physics as a paradigmatic model of superradiance in both coherent and dissipative cavity-QED settings 13; 30. The model exhibits a $\mathbb{Z}_{2}$ symmetry-breaking transition 21; 57; 17, near which the low-energy and long-time behavior is governed by universal critical exponents that dictate how observables scale with tuning parameters 6; 49. Extracting these exponents experimentally is a central way to identify the underlying universality class, as has been done in many condensed-matter systems and, more recently, in quantum simulators with a finite number of particles 26; 16; 28; 29; 36; 62.

Despite the simplicity of the Dicke Hamiltonian, a quantitative extraction of its critical exponents has remained challenging. Early cavity-QED realizations reached large atom numbers, but spatial inhomogeneity of the cavity mode drove the light–matter coupling far from the ideal collective limit 3; 4; 31. Recent tweezer arrays offer nearly uniform couplings, but only at mesoscopic atom numbers 50; 22; 58; 10. In this regime, all-to-all coupling makes the critical correlation time grow slowly with $N$ 38; 20; 40; 12, so finite-size corrections readily obscure its asymptotic scaling and bias estimates of the dynamical exponent 36; 22; 8. Unavoidable photon loss 14; 3 is a relevant perturbation that drives the system toward a distinct dissipative universality class 42; 33; 7; at mesoscopic sizes and finite loss, the crossover between the coherent and dissipative limits can be broad, so the apparent exponents need not coincide with those of either asymptotic fixed point.

Figure 1: (a) Phase diagram of the Dicke model in the plane of light-matter coupling $g$ and dissipation $\kappa$. Two distinct universality classes are identified at large $N$: closed (green) and open (red). The ramp protocol is defined by the ramping speed $s$. (b) Crossover regions set by size $N$, dissipation $\kappa$, and speed $s$ in the ramping dynamics. (c) Scaling protocol at mesoscopic scales that encodes the irrelevant exponent $\omega$. The exponents $\nu$, $z$, and $\mu$ are then extracted independently, and the KZ relation $1/\mu=1/\nu+z$ is tested.

Critical properties can be probed beyond the static and steady-state cases. Under a finite-speed ramp through the transition, Kibble–Zurek (KZ) scaling links the critical exponents to the universal dynamical response 27; 64; 63; 44; 15; 48, a relation verified across diverse systems 52; 43; 26; 62; 29. For the Dicke model, the absence of spatial structure eliminates the usual domain-wall picture. Nevertheless, dynamical scaling laws have been formulated in the asymptotic large-$N$ regime for various closed all-to-all models 1; 11; 60. At mesoscopic sizes, however, strong finite-size drift obscures the asymptotic scaling, and a controlled extraction of the exponents from experimentally accessible systems remains lacking. Moreover, the signatures of the coherent and dissipative fixed points become intertwined with the finite ramp speed, making the dynamical critical response difficult to characterize 45; 24 and calling for a protocol that can disentangle these effects.

In this Letter, we develop a unified mesoscopic scaling theory for Dicke superradiant criticality that treats static and driven scaling, as well as coherent and dissipative criticality, on the equal footing. Using a large-$N$ field theory, we identify distinct universality classes for the coherent and dissipative fixed points [Fig. 1(a)]. Because finite size, photon loss, and finite ramp speed introduce competing dynamical time scales, broad crossover regimes emerge in which neither fixed-point response is cleanly visible [Fig. 1(b)]. To overcome this difficulty, we construct a finite-size scaling protocol that retains the leading irrelevant correction 59; 5 and extracts the thermodynamic exponents from mesoscopic static, spectral, and ramped data [Fig. 1(c)]. This analysis recovers the large-$N$ exponents and verifies the KZ relation in both cases, establishing universal driven scaling in a fully connected model without spatial domain formation.

Dicke model: closed and open.— We consider the Dicke model, in which $N$ two-level atoms with energy splitting $\omega_{z}$ couple to a single cavity mode $\hat{a}$ of frequency $\omega_{0}$. The Hamiltonian reads 13

$$\hat{H}=\omega_{0}\hat{a}^{\dagger}\hat{a}+\omega_{z}\hat{J}_{z}+\frac{g}{\sqrt{N}}(\hat{a}+\hat{a}^{\dagger})\hat{J}_{x},$$ (1)

where $\hat{J}_{\alpha}=\frac{1}{2}\sum_{i=1}^{N}\sigma_{i}^{\alpha}$ for $\alpha=x,z$, and $g$ is the light-matter coupling strength. The cavity may additionally undergo photon loss at a rate $\kappa$. For $\kappa=0$, the model is closed and evolves coherently; for $\kappa\neq 0$, it is open and evolves dissipatively according to the Lindblad equation 14

$$\dot{\hat{\rho}}=-\mathrm{i}[\hat{H},\hat{\rho}]+\kappa\left(2\hat{a}\hat{\rho}\hat{a}^{\dagger}-\hat{a}^{\dagger}\hat{a}\hat{\rho}-\hat{\rho}\hat{a}^{\dagger}\hat{a}\right).$$ (2)

In the thermodynamic limit, the model undergoes a superradiant phase transition. The $\mathbb{Z}_{2}$ symmetry-breaking phase is characterized by a nonzero order parameter, $\mathrm{Tr}(\hat{\rho}\,\hat{x})\neq 0$, where $\hat{x}=(\hat{a}+\hat{a}^{\dagger})/\sqrt{2\omega_{0}}$, and the critical point is located at $g_{c}=\sqrt{(\kappa^{2}+\omega_{0}^{2})\,\omega_{z}/\omega_{0}}$.

Large-$N$ field theory.— For both the closed and open Dicke models, the critical exponents can be analyzed within a unified Keldysh field theory 7; 51; 25, with the leading $1/N$ corrections kept explicitly. The theory is written in terms of the real order-parameter field $x(t)$ associated with the Hermitian operator $\hat{x}$ defined above. On the Keldysh contour, this field is doubled into $x_{+}(t)$ and $x_{-}(t)$ on the forward and backward time branches. After rotating to the Keldysh basis, with $x_{\rm cl}=(x_{+}+x_{-})/\sqrt{2}$ and $x_{\rm q}=(x_{+}-x_{-})/\sqrt{2}$, the effective action takes the form

$$S=-\int dt \Big[x_{\rm q}\left(2\kappa\partial_{t}+K\partial_{t}^{2}+Mr_{\kappa}\right)x_{\rm cl}$$ (3)

$$+\frac{u}{N}\left(x_{\rm q}x_{\rm cl}^{3}+x_{\rm q}^{3}x_{\rm cl}\right)-{\rm i}D_{\kappa}x_{\rm q}^{2}\Big].$$

The tuning parameter is the distance from criticality

$$r_{\kappa}=\frac{(g/g_{c})^{2}-1}{\sqrt{\omega_{0}/\omega_{z}+\kappa^{2}/\omega_{0}\omega_{z}}}.$$ (4)

The noise strength $D_{\kappa}$ distinguishes the closed and open cases: it vanishes for $\kappa=0$ and is nonzero for $\kappa\neq 0$. The remaining coefficients $K$, $M$, and $u$ are nonzero in both cases; their definitions, together with the derivation of the action, are given in the End Matter.

This action provides the starting point for power counting. We assign scaling dimensions with respect to the particle number $N$, rather than a linear size $L$, and adopt the convention $[N]\equiv-1$ and $[t]\equiv-z$. A key feature of Eq. (3) is that the near-critical behavior is governed by a single tuning parameter $r_{\kappa}$, rather than by $g$ and $\kappa$ separately, even in the presence of dissipation; see the End Matter for details. The near-critical scaling is therefore controlled by the combination $r_{\kappa}N^{1/\nu}$, which defines the scaling dimension $[r_{\kappa}]\equiv 1/\nu$.

Closed case: $\kappa=0$. In the closed model, the two time branches are not coupled by noise, and the classical and quantum fields are treated on equal footing: $[x_{\rm cl}]=[x_{\rm q}]$. At criticality, $r_{\kappa}=0$, and the coherent kinetic term competes with the leading $1/N$ term. Requiring both terms to be dimensionless gives $\int dt\,x_{\rm q}\partial_{t}^{2}x_{\rm cl}\rightarrow z+2[x_{\rm cl}]=0$, and $\int dt\,\frac{1}{N}\left(x_{\rm q}x_{\rm cl}^{3}+x_{\rm q}^{3}x_{\rm cl}\right)\rightarrow-z+1+4[x_{\rm cl}]=0$. Hence $[x_{\rm cl}]=[x_{\rm q}]=-1/6$ and $z=1/3$. Away from criticality, the mass term $r_{\kappa}$ scales in the same way as $\partial_{t}^{2}$, so $1/\nu=2z$, giving $\nu=3/2$. This is the all-to-all Ising universality class 54.

Open case: $\kappa\neq 0$. In the open model, the noise term couples the two time branches and changes the scaling structure. At long times, the dissipative term $\kappa\partial_{t}$ is more relevant than the coherent term $\partial_{t}^{2}$. The term $x_{\rm q}\partial_{t}x_{\rm cl}$ then implies $[x_{\rm cl}]=-[x_{\rm q}]$. At criticality, the relevant competition is between the dissipative kinetic term, the noise term, and the leading $1/N$ term $x_{\rm q}x_{\rm cl}^{3}/N$ (since $x_{\rm cl}$ is more relevant than $x_{\rm q}$). Power counting these terms gives $[x_{\rm cl}]=-[x_{\rm q}]=-1/4$ and $z=1/2$. Finally, the mass term $r_{\kappa}$ scales in the same way as $\partial_{t}$, so $1/\nu=z$, giving $\nu=2$. This is the Model-A universality class 7.

The above analysis concerns the static tuning parameter $r_{\kappa}$. For ramping dynamics, we take $r_{\kappa}=st$, where $s$ is the ramping speed and define $[s]=1/\mu$; see Fig. 1(a). Since $[s]=[r_{\kappa}]-[t]$, we obtain $1/\mu=1/\nu+z$, which is the KZ scaling relation. In the following sections, we recover these large-$N$ critical exponents using the mesoscopic scaling protocol.

Scaling framework at mesoscopic scales.— We now formulate a scaling framework for extracting critical exponents from mesoscopic all-to-all interacting systems. As a reference point, consider a locally interacting system near criticality, with $r$ measuring the distance from the critical point. In the thermodynamic limit, the singular behavior is governed by a diverging correlation length $\xi\sim|r|^{-\nu}$ and correlation time $\tau\sim\xi^{z}\sim|r|^{-z\nu}$. In a finite system of linear size $L$, this divergence is cut off when $\xi\sim L$, giving the standard scaling variable $rL^{1/\nu}$. If the system is instead driven through the transition at a finite speed, $r=st$, the growth of $\xi$ is cut off by the KZ scale $\xi_{\rm KZ}\sim s^{-\mu}$, where $1/\mu=1/\nu+z$. Thus $\nu$, $z$, and $\mu$ characterize the static, dynamical, and driven scaling properties of the transition, and together provide key fingerprints of its universality class.

For the all-to-all systems considered here, there is no natural notion of a linear length scale. The particle number $N$ therefore replaces $L$ as the finite-size scaling variable. We keep the notation $\nu$, but it should now be understood as the finite-$N$ critical-window exponent: the width of the critical region scales as $N^{-1/\nu}$. A generic observable $Q$ in equilibrium, or in the steady state of an open system, is then expected to obey $Q(N,r)=N^{y_{Q}/\nu}\mathcal{F}_{Q}\!\left(rN^{1/\nu}\right)$, where $\mathcal{F}_{Q}$ is the scaling function associated with $Q$ 6. A dimensionless observable has $y_{Q}=0$ and can be used to extract $\nu$, while a relaxation time has $y_{Q}/\nu=z$ and can be used to extract $z$.

This scaling form is the natural starting point, but it is not sufficient in the mesoscopic regime. In all-to-all models, the approach to the thermodynamic scaling limit can be slow, so leading irrelevant corrections may strongly affect the accessible sizes and, if ignored, produce apparent exponents that differ substantially from the thermodynamic values. The central idea of our protocol is to first determine the leading irrelevant correction and then use the same correction consistently for other observables [Fig. 1(c)].

Stage I: Fix the irrelevant exponent $\omega$. We begin with a dimensionless observable $U$, for example a Binder ratio 34. Let $u$ denote the amplitude of the leading irrelevant scaling field, and let $\omega$ be the corresponding correction exponent. The scaling form becomes

$$U(N,r) =\mathcal{F}_{U}\!\left(rN^{1/\nu},uN^{-\omega}\right)$$ (5)

$$=\sum_{j=0}f_{U,j}\!\left(rN^{1/\nu}\right)N^{-j\omega},$$

where the second line is an expansion in the irrelevant variable, with the coefficients absorbing powers of $u$. At the critical point , this reduces to $U(N,0)=\sum_{j=0}f_{U,j}(0)\,N^{-j\omega}$. Thus the finite-$N$ drift of $U(N,0)$ can be fitted to a truncated series in powers of $N^{-\omega}$, allowing $\omega$ to be determined.

Stage II: Extract $z$, $\nu$, and $\mu$. The remaining exponents are then obtained from independent measurements, all using the same irrelevant correction exponent $\omega$.

$\bullet$  For $z$: The exponent $z$ is encoded in the relaxation time $\tau$, and hence in the inverse gap $1/\Delta$. The relevant gap is the excitation gap in a Hamiltonian setting or the Liouvillian gap in a dissipative setting 41. At criticality, it scales as

$$\Delta(N,0)=N^{-z}\sum_{j=0}f_{\Delta,j}(0)N^{-j\omega},$$ (6)

so the size dependence of the critical gap determines $z$.

$\bullet$  For $\nu$: Differentiating Eq. (5) with respect to $r$ and evaluating at criticality gives

$$\left.\partial_{r}U(N,r)\right|_{r=0}=N^{1/\nu}\sum_{j=0}f_{U,j}^{(1)}N^{-j\omega},$$ (7)

where $f_{U,j}^{(1)}$ denotes the first derivative of $f_{U,j}$ at the critical point. Thus the size dependence of the derivative of the dimensionless quantity determines $\nu$ .

$\bullet$  For $\mu$: The exponent $\mu$ is extracted from driven dynamics. We consider a linear ramp through the critical point with speed $s$. For a closed system, the initial state is the ground state far from the critical point; for an open system, the corresponding preparation is the steady state. During the ramp, we measure the instantaneous dimensionless quantity $U(N,r,s)$, whose scaling form is 19; 9; 32; 39; 23

$$U(N,r,s)=\mathcal{F}_{U}\!\left(rN^{1/\nu},sN^{1/\mu},uN^{-\omega}\right).$$ (8)

We focus on the slow-ramp regime, where the response is controlled by the critical point rather than by fast-quench dynamics. Expanding Eq. (8) in the scaling variable $sN^{1/\mu}$, let $n$ denote the order of the first nonvanishing speed correction. Then

$$\left.\partial^{n}_{s}U(N,r,s)\right|_{r,s=0}\!\!=N^{n/\mu}\sum_{j=0}\tilde{f}_{U,j}^{(n)}N^{-j\omega},$$ (9)

where $\tilde{f}_{U,j}^{(n)}$ denotes the corresponding $n$-th derivative, taken with respect to the speed scaling field, of the $j$-th coefficient in the irrelevant correction expansion. Thus the size dependence of this leading small-speed response determines $\mu$.

This protocol shows that, once $\omega$ is fixed, the exponents $z$, $\nu$, and $\mu$ can be extracted independently. Their agreement with $1/\mu=1/\nu+z$ then provides a nontrivial check that the mesoscopic data analysis recovers the universal KZ scaling, as we verify below for both the closed and open Dicke models.

Numerical results at mesoscopic scales.— In the thermodynamic limit, the superradiant transition is characterized by spontaneous $\mathbb{Z}_{2}$ symmetry breaking and a nonzero order parameter $\mathrm{Tr}(\hat{\rho}\hat{x})\neq 0$. At finite $N$, however, the $\mathbb{Z}_{2}$ symmetry is not truly broken, so the transition must be diagnosed through fluctuations. We use the dimensionless Binder ratio 29; 56

$$U(N)=1-\frac{\mathrm{Tr}(\hat{\rho}\hat{x}^{4})}{3[\mathrm{Tr}(\hat{\rho}\hat{x}^{2})]^{2}}.$$ (10)

Here $\hat{\rho}$ can denote the ground-state density matrix $|\mathrm{GS}\rangle\langle\mathrm{GS}|$, the steady state of the Lindblad equation, or a time-dependent state generated by Eq. (2).

Figure 2: Closed Dicke model. (a) Critical Binder ratio used to fix the leading irrelevant exponent $\omega$. The red dashed line marks extrapolated $U(N\to\infty)$. (b)–(d) Extraction of $z$, $\nu$, and $\mu$ from the critical gap, Binder-ratio derivative, and small-speed response, respectively. Gray dashed curves include corrections to scaling up to $j=2$; green dashed lines show the leading power laws with the fitted exponents. For all figures, we set $\omega_{0}=\omega_{z}=1$.

Figure 3: Open Dicke model. (a) Joint critical Binder-ratio fit for different $\kappa\in[0.1,3]$, used to fix $\omega$. Red dashed curves show extrapolated $U(N)$ vs $\kappa$ at fixed $N=10^{2},10^{3},10^{4}$, and $N\to\infty$, ordered from sparse to dense. (b)–(d) Extraction of $z$, $\nu$, and $\mu$ from the Liouvillian gap, Binder-ratio derivative, and linear small-speed response, respectively. Insets in (b)–(d) show the corresponding $\kappa$ dependence.

Closed case: We first consider the ground state. Figure 2(a) shows the finite-$N$ dependence of the Binder ratio at the critical point $r=0$. Fitting the drift of $U(N)$ with Eq. (5), truncated to $j=2$ in the irrelevant correction, gives $\omega=0.372$ using sizes up to $N=100$. The extrapolated infinite-$N$ value, $U(N\rightarrow\infty)\approx 0.102$, is shown by the red dashed line in Fig. 2(a). Its visible separation from the $N=100$ data highlights the strong finite-size corrections in this mesoscopic regime. With this $\omega$ fixed, the critical gap gives $z=0.33$ [Fig. 2(b)], while the critical derivative of the Binder ratio gives $\nu=1.52$ [Fig. 2(c)]. These values recover the thermodynamic exponents to two significant digits.

We next study ramping dynamics by driving $r$ from deep in the normal phase through the critical point with speed $s$. For the closed system, time-reversal symmetry forbids a linear-in-$s$ correction to $U(s)$. The leading speed dependence is therefore quadratic. We compute $\partial^{2}_{s}U$ and fit its size dependence. This gives $\mu=0.99$, close to the large-$N$ prediction $\mu=1$.

Open case: For the open system, we first analyze the steady state. Although the loss rate $\kappa$ can vary, all $\kappa\neq 0$ cases are expected to belong to the same universality class. We therefore impose that both $\omega$ and the critical value of $U(N\rightarrow\infty)$ are independent of $\kappa$. A joint fit over $\kappa\in[0.1,3]$ gives $\omega=0.363$, as shown in Fig. 3(a).

The dynamical exponent $z$ is obtained from the Liouvillian gap. For $\kappa=1$, exact diagonalization gives $z=0.48$. For generic $\kappa$, however, the small sizes accessible to exact diagonalization, $N\leq 24$, show crossings and rearrangements among low-lying Liouvillian modes, which obscure the asymptotic scaling. To avoid these crossover effects, we perform a numerical large-$N$ expansion, described in the End Matter. For $N\in[10^{5},10^{7}]$, the extracted exponent approaches $z=1/2$ for different $\kappa$, as shown in the inset of Fig. 3(b). The exponent $\nu$ is extracted from the critical derivative of the Binder ratio. At $\kappa=1$, we find $\nu=2.02$ [Fig. 3(c)], while fits at different $\kappa$ fluctuate around the expected value, giving $\nu=1.97(9)$ in the inset.

For ramping dynamics, because dissipation breaks the time-reversal structure, the Binder ratio can acquire a linear correction in the ramp speed. We therefore compute $\partial_{s}U$ and fit its size dependence. At $\kappa=1$, we find $\mu=0.96$; across different $\kappa$, we obtain $\mu=1.02(5)$.

Time scales and crossover.— We now discuss how the finite size $N$, the ramp speed $s$, and the dissipation $\kappa$ compete in the ramping dynamics [Fig. 1(b)]. In the closed model, the relaxation time diverges as $\tau\sim|r|^{-1/2}$ [cf. $z\nu=1/2$]. For a linear ramp $r=st$, the KZ time $t_{\rm KZ}$ is defined by $\tau(r_{\rm KZ})\sim|r_{\rm KZ}|/s$ 64, which gives $t_{\rm KZ}\sim s^{-1/3}$. On the other hand, the critical finite-size gap scales as $\Delta_{N}\sim N^{-1/3}$, yielding the finite-size time scale $t_{N}\sim N^{1/3}$. Comparing $t_{\rm KZ}$ with $t_{N}$ gives the scaling variable $sN$ [cf. $\mu=1$]. Thus $sN\ll 1$ is the finite-size-dominated perturbative regime, while $sN\gtrsim 1$ marks the onset of the conventional KZ regime. Figure 4(a) mainly probes the former: the response is nearly flat at small $sN$, consistent with the absence of a linear-in-$s$ correction.

We next consider how $\kappa$ drives the system away from the closed fixed point. In the large-$N$ long-time limit, any finite $\kappa$ is a relevant perturbation: the coherent closed fixed point is unstable, and the dynamics flows to the dissipative fixed point, with $z$ changing from $1/3$ to $1/2$. Near the closed fixed point, the competition between the coherent kinetic term and photon loss yields $\partial_{t}^{2}\sim\kappa\,\partial_{t}$, which defines a dissipative timescale $t_{\kappa}\sim\kappa^{-1}$. In a finite system, the closed critical dynamics is cut off at $t_{N}\sim N^{1/3}$. Hence, the crossover scaling variable is $\kappa N^{1/3}$: small systems remain closed-like when $\kappa N^{1/3}\ll 1$ and enter the dissipative regime otherwise. This is consistent with scaling theory, where the dissipation rate enters as an additional scaling variable through $Q\sim\mathcal{F}_{Q}(r_{0}N^{2/3},\kappa N^{1/3})$ 61; 47. However, the scaling protocol should not be limited to this weak-dissipation crossover. Since the limit $\kappa\to 0$ is singular, accessing the dissipative fixed point itself requires a separate scaling analysis. This is precisely what our protocol achieves: it extracts the distinct dissipative fixed-point scaling form $\mathcal{F}_{Q}(r_{\kappa}N^{1/2})$ [cf. Fig. 3(c)].

Figure 4: (a) Ramped Binder ratio in the closed Dicke model for different system sizes; red dashed lines indicate the $N\to\infty$ extrapolation. (b) Comparison of the ramped Binder ratio between the closed and open Dicke models for different $\kappa$.

For ramping at nonzero $\kappa$, the KZ time $t_{\rm KZ}\sim s^{-1/3}$ enters. The same comparison gives $t_{\rm KZ}\sim t_{\kappa}$, or $s\sim\kappa^{3}$. Thus fast ramps with $s\gg\kappa^{3}$ do not give the system enough time to relax and can appear closed-like, while slow ramps cross over to open dissipative scaling. This behavior is visible in Fig. 4(b). For large $\kappa$, such as $\kappa=1,2$, the response is already open-like and shows a linear small-speed decay. For small $\kappa$, such as $\kappa=0.1$, the slowest ramps exhibit the open linear response, but increasing $s$ eventually bends the curve toward the flatter, closed-like speed dependence seen in Fig. 4(a). These intertwined crossovers explain why the raw ramp data can look complicated, and why a unified finite-size scaling treatment that includes the irrelevant corrections is needed to extract the correct exponents.

Outlook.— This framework is not tied to cavity-QED platforms. It applies more broadly to long-range and all-to-all quantum simulators, from trapped ions with tunable power-law interactions 36 to Sachdev–Ye–Kitaev-type models related to non-Fermi-liquid criticality 55. In many such systems, a controlled large-$N$ theory may not be available, so extracting the leading irrelevant structure directly from data becomes essential. Our results therefore provide a practical route for probing universality when the thermodynamic limit is clear in principle but experimentally out of reach.

Acknowledgments.— We are grateful to Hui Zhai for stimulating ideas. We also thank Chengshu Li for reading the manuscript and providing valuable feedback. H.L. is supported by Quantum Science and Technology-National Science and Technology Major Project (Grant No. 2025ZD0300400), China National Postdoctoral Program for Innovative Talents (Grant No. BX2026034) , and National Natural Science Foundation of China (Grant No. 12547168). H.W. is supported by China Postdoctoral Science Foundation under Grant No. 2024M751609 and Postdoctoral Fellowship Program of CPSF under Grant No. GZC20231364.

End Matter

Appendix A: Keldysh action.— Here we derive Eq. (3) of the main text. We represent the collective spin degrees of freedom by a Holstein–Primakoff boson $\hat{b}$,

$$\hat{J}_{z} =\hat{b}^{\dagger}\hat{b}-N/2,$$ (A1)

$$\hat{J}_{x} =\frac{1}{2}\left[\hat{b}^{\dagger}\sqrt{N-\hat{b}^{\dagger}\hat{b}}+\sqrt{N-\hat{b}^{\dagger}\hat{b}}\,\hat{b}\right].$$

Substituting this representation into the Dicke Hamiltonian and expanding to leading order in $1/N$, we obtain $\hat{H}=\hat{H}_{2}+\hat{H}_{1/N}$, with

$$\hat{H}_{2}=\omega_{0}\hat{a}^{\dagger}\hat{a}+\omega_{z}\hat{b}^{\dagger}\hat{b}+\frac{g}{2}(\hat{a}+\hat{a}^{\dagger})(\hat{b}+\hat{b}^{\dagger}),$$ (A2)

$$\hat{H}_{1/N}=-\frac{g}{4N}(\hat{a}+\hat{a}^{\dagger})(\hat{b}^{\dagger}\hat{b}\hat{b}+\hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}).$$

We now formulate the theory on the Keldysh contour. Including the Lindblad jump term from Eq. (2), the action reads

$$S=\!\int\!dt \Big[\mathrm{i}(\bar{a}_{+}\partial_{t}a_{+}-\bar{a}_{-}\partial_{t}a_{-})+\mathrm{i}(\bar{b}_{+}\partial_{t}b_{+}-\bar{b}_{-}\partial_{t}b_{-})$$ (A3)

$$-(H_{+}-H_{-})-\mathrm{i}\kappa(2a_{+}\bar{a}_{-}-\bar{a}_{+}a_{+}-\bar{a}_{-}a_{-})\Big].$$

To proceed, we introduce real coordinate and momentum variables for the cavity mode, $x=(a+\bar{a})/\sqrt{2\omega_{0}}$ and $p=\mathrm{i}\sqrt{\omega_{0}/2}\,(\bar{a}-a)$. We similarly introduce the atomic coordinate $y=(b+\bar{b})/\sqrt{2\omega_{z}}$ and its conjugate momentum $q=\mathrm{i}\sqrt{\omega_{z}/2}\,(\bar{b}-b)$.

We first treat the quadratic part $S_{2}$. Performing the Gaussian integration over the momenta $p$ and $q$ gives

$$S_{2}= -\int dt\,\left[x_{\rm q}(\partial_{t}^{2}+2\kappa\partial_{t}+\omega_{0}^{2}+\kappa^{2})x_{\rm cl}-\mathrm{i}D_{\kappa}x_{\rm q}^{2}\right]$$ (A4)

$$-(S_{y}^{+}-S_{y}^{-}).$$

In the first line, the Keldysh rotation for $x$ follows the convention used in the main text. We have omitted the term $\frac{\mathrm{i}\kappa}{\omega_{0}}(\partial_{t}x_{\rm q})^{2}$, which is irrelevant at low frequencies. The second line contains the remaining action for the massive atomic coordinate $y$, with

$$S_{y}=\int dt\,\left[\frac{1}{2}y(\partial_{t}^{2}+\omega_{z}^{2})y+g\sqrt{\omega_{0}\omega_{z}}\,xy\right].$$ (A5)

We now turn to the $1/N$ contribution, which supplies the leading interaction vertex. In terms of $x$, $y$, and $q$, the interaction reads

$$H_{1/N}=-\frac{g\sqrt{\omega_{0}}\,\omega_{z}^{3/2}}{4N}xy^{3}-\frac{g}{4N}\sqrt{\frac{\omega_{0}}{\omega_{z}}}xyq^{2}.$$ (A6)

The $q^{2}$ term gives only irrelevant contributions to the critical theory, so the $xy^{3}/N$ term provides the leading relevant vertex.

We are now ready to integrate out the massive field $y$. Without the $1/N$ correction, one can eliminate $y$ by a Gaussian integration. With the $1/N$ vertex present, it is useful to first shift $y\rightarrow y-\mathcal{G}^{-1}x$, so that the shifted $y$ field has zero mean. The quadratic action (A5) becomes

$$S_{y}\rightarrow\int dt\,\left[\frac{1}{2}y(\partial_{t}^{2}+\omega_{z}^{2})y\right]-\int dt\,\left[\frac{1}{2}g\sqrt{\omega_{0}\omega_{z}}\,x\mathcal{G}^{-1}x\right],$$ (A7)

where, in the low-frequency limit,

$$\mathcal{G}^{-1}=g\sqrt{\omega_{0}\omega_{z}}(\partial_{t}^{2}+\omega_{z}^{2})^{-1}\approx g\sqrt{\frac{\omega_{0}}{\omega_{z}^{3}}}(1-\omega_{z}^{-2}\partial_{t}^{2}).$$ (A8)

This shift removes the bilinear $xy$ coupling. The second term in (A7) renormalizes the coherent kinetic and mass terms of the cavity field.

We then treat the $xy^{3}/N$ vertex perturbatively using the shifted zero-mean Gaussian field $y$,

$$S_{1/N}=-\int dt\,\frac{g\sqrt{\omega_{0}}\,\omega_{z}^{3/2}}{4N}\left\langle x(y-\mathcal{G}^{-1}x)^{3}\right\rangle_{y},$$ (A9)

where $\langle\cdot\cdot\cdot\rangle_{y}$ denotes averaging over the zero-mean Gaussian action. The terms $x\langle y^{3}\rangle/N$ and $x^{3}\langle y\rangle/N$ vanish, while $x^{2}\langle y^{2}\rangle/N$ gives an $\mathcal{O}(1/N)$ mass renormalization and shifts the critical point. The only relevant nonlinear term is therefore the quartic contribution $x^{4}/N$. Expressed on the forward–backward contour as $x_{+}^{4}-x_{-}^{4}$ and rotated to the Keldysh basis, it gives

$$S_{1/N}\approx-\int dt\,\frac{u}{N}\left(x_{\rm q}x_{\rm cl}^{3}+x_{\rm q}^{3}x_{\rm cl}\right),\quad u=\frac{g^{4}\omega_{0}^{2}}{2\omega_{z}^{3}}.$$ (A10)

Collecting the quadratic renormalizations, the effective quadratic action takes the form

$$S_{2}=-\int dt\,\left[x_{\rm q}\left(2\kappa\partial_{t}+K\partial_{t}^{2}+\widetilde{M}\right)x_{\rm cl}-\mathrm{i}D_{\kappa}x_{\rm q}^{2}\right],$$ (A11)

with renormalized coefficients

$$K=1+\frac{g^{2}\omega_{0}}{\omega_{z}^{3}},\quad\widetilde{M}=-(\omega_{0}^{2}+\kappa^{2})\left[(g/g_{c})^{2}-1\right].$$ (A12)

Near the critical point, $\widetilde{M}$ is linear in $g-g_{c}(\kappa)$ and controls the transition. The coefficients $K$ and $u$ are nonsingular at criticality and can be evaluated at $g=g_{c}(\kappa)$.

Appendix B: Universal theory for $\kappa\neq 0$.— We now show that, for any finite loss rate $\kappa$, the low-energy dissipative theory can be written in terms of a single scaled distance from criticality. The key observation is that the explicit $\kappa$ dependence of the quadratic and nonlinear coefficients can be absorbed entirely into a rescaling of the fields and time.

In the dissipative regime, the damping term $\kappa\partial_{t}$ dominates over $\partial_{t}^{2}$ at low frequencies. We therefore drop $\partial_{t}^{2}$ and perform the rescaling

$$\begin{cases}t=\frac{8\kappa N^{1/2}\omega_{z}^{2}}{g_{c}^{3}(\kappa)\omega_{0}}\,\tilde{t},\\ x_{\rm cl}(t)=N^{1/4}\left(\frac{2\omega_{z}}{g_{c}(\kappa)\omega_{0}}\right)^{1/2}\phi_{\rm cl}(\tilde{t}),\\ x_{\rm q}(t)=\frac{1}{2\kappa N^{1/4}}\left(\frac{g_{c}(\kappa)\omega_{0}}{2\omega_{z}}\right)^{1/2}\phi_{\rm q}(\tilde{t}).\end{cases}$$ (B1)

With this choice, the effective action becomes

$$S=-\!\int\!d\tilde{t}\left[\phi_{\rm q}\!\left(\partial_{\tilde{t}}-4r_{\kappa}N^{1/2}\right)\!\phi_{\rm cl}+4\phi_{\rm q}\phi_{\rm cl}^{3}-\mathrm{i}\phi_{\rm q}^{2}\right],$$ (B2)

where the scaled distance from criticality is

$$r_{\kappa}=\frac{(g/g_{c})^{2}-1}{\sqrt{\omega_{0}/\omega_{z}+\kappa^{2}/\omega_{0}\omega_{z}}}.$$ (B3)

Thus the dissipative low-energy theory depends on the microscopic parameters only through the scaling variable $r_{\kappa}N^{1/2}$. This gives the open-system exponent $\nu=2$. Hence, the corresponding dissipative fixed-point scaling takes the form $Q\sim\mathcal{F}_{Q}(r_{\kappa}N^{1/2})$. This should be distinguished from the weak-dissipation crossover scaling near the closed fixed point, $Q\sim\mathcal{F}_{Q}(r_{0}N^{2/3},\kappa N^{1/3})$, where $r_{0}=(g/g_{c})^{2}-1$.

This rescaling also clarifies the form of the mass term in Eq. (3). The bare mass of Appendix A factorizes as $\widetilde{M}=Mr_{\kappa}$, isolating the scaling variable $r_{\kappa}$, with

$$M=-g_{c}^{3}(\kappa)\frac{\omega_{0}}{\omega_{z}^{2}}.$$ (B4)

Appendix C: Hartree–Fock–Bogoliubov analysis of the open-system soft mode.— At the small sizes accessible to exact Liouvillian diagonalization, the lowest nonzero eigenmode can switch between different branches as $N$ and $\kappa$ are varied. This obscures the asymptotic scaling of the critical soft mode.

We treat the $1/N$ interaction in Eq. (A2) at the HFB level. With $n_{b}=\langle\hat{b}^{\dagger}\hat{b}\rangle$ and $m_{b}=\langle\hat{b}\hat{b}\rangle$,

$$\hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}+\hat{b}^{\dagger}\hat{b}\hat{b}\simeq(2n_{b}+\mathrm{Re}\,m_{b})(\hat{b}+\hat{b}^{\dagger}).$$ (C1)

The effective Hamiltonian therefore becomes

$$\hat{H}_{\rm HFB}=\omega_{0}\hat{a}^{\dagger}\hat{a}+\omega_{z}\hat{b}^{\dagger}\hat{b}+\frac{g_{\rm eff}}{2}(\hat{a}+\hat{a}^{\dagger})(\hat{b}+\hat{b}^{\dagger}),$$ (C2)

where $g_{\rm eff}=g[1-(2n_{b}+\mathrm{Re}\,m_{b})/2N].$ Together with the cavity loss in Eq. (2), this gives

$$\dot{\hat{v}}=A\hat{v}+\hat{\xi},\qquad\hat{v}=(\hat{a},\hat{a}^{\dagger},\hat{b},\hat{b}^{\dagger})^{T}.$$ (C3)

The drift matrix is

$$A=\begin{pmatrix}-\kappa-\mathrm{i}\omega_{0}&0&-\mathrm{i}g_{\rm eff}/2&-\mathrm{i}g_{\rm eff}/2\\ 0&-\kappa+\mathrm{i}\omega_{0}&\mathrm{i}g_{\rm eff}/2&\mathrm{i}g_{\rm eff}/2\\ -\mathrm{i}g_{\rm eff}/2&-\mathrm{i}g_{\rm eff}/2&-\mathrm{i}\omega_{z}&0\\ \mathrm{i}g_{\rm eff}/2&\mathrm{i}g_{\rm eff}/2&0&\mathrm{i}\omega_{z}\end{pmatrix}.$$ (C4)

The input noise satisfies,

$$\langle\hat{\xi}(t)\hat{\xi}^{\dagger}(t^{\prime})\rangle=D\delta(t-t^{\prime}),\quad D=\mathrm{diag}(2\kappa,0,0,0).$$ (C5)

The covariance matrix $C_{ij}=\langle\hat{v}_{i}\hat{v}_{j}^{\dagger}\rangle$ satisfies

$$AC+CA^{\dagger}+D=0.$$ (C6)

The moments $n_{b}=\langle\hat{b}^{\dagger}\hat{b}\rangle$ and $m_{b}=\langle\hat{b}\hat{b}\rangle$ extracted from $C$ are inserted back into $g_{\rm eff}$, and this loop is iterated to self-consistency.

Figure 5: Soft-mode spectral gap of the open Dicke model at the critical point, obtained within the Hartree–Fock–Bogoliubov approximation for $\kappa=1$. The calculation is performed over the range $N=10^{5}$–$10^{7}$, where the asymptotic $N^{-1/2}$ decay becomes clearly visible.

After convergence, the homogeneous equation has modes $\delta\hat{v}(t)\propto e^{\lambda t}$, where $\lambda$ are eigenvalues of $A$. The soft-mode decay rate is therefore

$$\Delta_{\rm soft}=-\max_{\mathrm{Re}\,\lambda<0}\mathrm{Re}\,\lambda.$$ (C7)

This quantity equals the Liouvillian gap for a quadratic open bosonic model. Here it is used only as an estimator of the asymptotic critical branch. As shown in Fig. 5, $\Delta_{\rm soft}\sim N^{-1/2}$ for $N=10^{5}$–$10^{7}$, confirming the dissipative exponent $z=1/2$ used in the main text.