Fluctuations and photon statistics in quantum metamaterial near the superradiant transition

TL;DR

This paper investigates photon fluctuations and counting statistics near the superradiant transition in the Dicke model, providing analytical and numerical insights into the behavior of photon number fluctuations and their full counting statistics.

Contribution

It introduces an effective Matsubara action formalism for analyzing photon fluctuations in the Dicke model at the superradiant transition, including full counting statistics.

Findings

Analytical expressions for photon number fluctuations and Fano factor.

Numerical simulations confirm the analytical results.

Identification of the Ginzburg-Levanyuk region near the transition.

Abstract

The analysis of single-mode photon fluctuations and their counting statistics at the superradiant phase transition is presented. The study concerns the equilibrium Dicke model in a regime where the Rabi frequency, related to a coupling of the photon mode with a finite-number qubit environment, plays a role of the transition's control parameter. We use the effective Matsubara action formalism based on the representation of Pauli operators as bilinear forms with complex and Majorana fermions. Then, we address fluctuations of superradiant order parameter and quasiparticles. The average photon number, the fluctuational Ginzburg-Levanyuk region of the phase transition and Fano factor are evaluated. We determine the cumulant generating function which describes a full counting statistics of equilibrium photon number. Exact numerical simulation of the superradiant transition demonstrates quantitative agreement with analytical calculations.

Equations

$$\hat{H}=\omega {\hat{\psi }}^{†}\hat{\psi }+\sum _{j=1}^N\frac{{ϵ}_j}{2}{\hat{\sigma }}_j^z+\sum _{j=1}^N{g}_j(\hat{\psi }{\hat{\sigma }}_j^{+}+{\hat{\psi }}^{†}{\sigma }_j^{-}).$$

$${\hat{\sigma }}_j^{+}=\sqrt{2}{\hat{c}}_j^{†}{\hat{d}}_j,{\hat{\sigma }}_j^{-}=\sqrt{2}{\hat{d}}_j{\hat{c}}_j.$$

$$Z=\int \mathcal{D}[\Psi ,\bar{c},c,d]\exp (-S[\Psi ,\bar{c},c,d])$$

$$\begin{gathered} S[\Psi ,\bar{c},c,d]=S_{\mathrm{ph}}[\Psi ]+S_{\mathrm{q}}[\bar{c},c,d]+ \\ +S_{\mathrm{int}}[\Psi ,\bar{c},c,d]+\ln Z_{\mathrm{ph}}{Z}_{\mathrm{q}}. \end{gathered}$$

$$S_{\mathrm{ph}}[\Psi ]=\underset{0}{\overset{\beta }{\int }}{\bar{\psi }}_{\tau }(-G_{\mathrm{ph};\tau -{\tau }^{\prime }}){\psi }_{{\tau }^{\prime }}d\tau ,$$

$$G_{\mathrm{ph};\tau -{\tau }^{\prime }}^{-1}={\delta }_{\tau -{\tau }^{\prime }}(-{\partial }_{{\tau }^{\prime }}-\omega ).$$

$${\psi }_n=T\underset{0}{\overset{\beta }{\int }}{\psi }_{\tau }{e}^{\mathrm{i}2\pi nT\tau }d\tau ,{\bar{\psi }}_n=T\underset{0}{\overset{\beta }{\int }}{\bar{\psi }}_{\tau }{e}^{-\mathrm{i}2\pi nT\tau }d\tau$$

$$G_{\mathrm{ph};n}^{-1}=\underset{0}{\overset{\beta }{\int }}{G}_{\mathrm{ph};\tau }^{-1}{e}^{\mathrm{i}2\pi nT\tau }d\tau =\mathrm{i}2\pi nT-\omega .$$

$$S_{\mathrm{ph}}[\Psi ]=\beta \sum _n{\bar{\psi }}_n(-G_{\mathrm{ph};n}^{-1}){\psi }_n.$$

$$\begin{gathered} S_{\mathrm{q}}[\bar{c},c,d] \\ =\frac{1}{2}\sum _{j=1}^N\underset{0}{\overset{\beta }{\int }}\left[\begin{array}{ccc}{\bar{c}}_j& c_j& d_j\end{array}\right](-{\mathbf{G}}_{j;\tau -{\tau }^{\prime }}^{-1})\left[\begin{array}{c}{c}_j\\ {\bar{c}}_j\\ d_j\end{array}\right]d\tau d{\tau }^{\prime }. \end{gathered}$$

$$-{\mathbf{G}}_{j;\tau -{\tau }^{\prime }}^{-1}={\delta }_{\tau -{\tau }^{\prime }}\left[\begin{array}{ccccc}{\partial }_{{\tau }^{\prime }}+{ϵ}_j& & 0& & 0\\ & & & & \\ 0& & {\partial }_{{\tau }^{\prime }}-{ϵ}_j& & 0\\ & & & & \\ 0& & 0& & {\partial }_{{\tau }^{\prime }}\end{array}\right].$$

$$\begin{gathered} S_{\mathrm{int}}[\Psi ,\bar{c},c,d] \\ =\frac{1}{2}\sum _{j=1}^N\underset{0}{\overset{\beta }{\int }}\left[\begin{array}{ccc}{\bar{c}}_j& c_j& d_j\end{array}\right]{\delta }_{\tau -{\tau }^{\prime }}{\mathbf{V}}_j[{\Psi }_{\tau }]\left[\begin{array}{c}{c}_j\\ {\bar{c}}_j\\ d_j\end{array}\right]d\tau d{\tau }^{\prime }. \end{gathered}$$

$${\mathbf{V}}_j[{\Psi }_{\tau }]=\sqrt{2}{g}_j\left[\begin{array}{ccccc}0& & 0& & -{\psi }_{\tau }\\ & & & & \\ 0& & 0& & {\bar{\psi }}_{\tau }\\ & & & & \\ -{\bar{\psi }}_{\tau }& & {\psi }_{\tau }& & 0\end{array}\right].$$

$$\ln Z_{\mathrm{ph}}=-\mathrm{Tr}\ln (-G_{\mathrm{ph};\tau -{\tau }^{\prime }}^{-1})$$

$$\ln Z_{\mathrm{q}}=\frac{1}{2}\sum _{j=1}^N\mathrm{Tr}\ln (-{\mathbf{G}}_{j;\tau -{\tau }^{\prime }}^{-1}).$$

$$\begin{gathered} S_{\mathrm{eff}}[\Psi ]=S_{\mathrm{ph}}[\Psi ]+\ln Z_{\mathrm{ph}}{Z}_{\mathrm{q}}- \\ -\frac{1}{2}\mathrm{Tr}\ln (-{\mathbf{G}}_{j;\tau -{\tau }^{\prime }}^{-1}+{\delta }_{\tau -{\tau }^{\prime }}{\mathbf{V}}_j[{\Psi }_{\tau }]). \end{gathered}$$

$$\begin{gathered} \mathrm{Tr}\ln (-{\mathbf{G}}_{j;\tau -{\tau }^{\prime }}^{-1}+{\delta }_{\tau -{\tau }^{\prime }}{\mathbf{V}}_j[{\Psi }_{\tau }])=\ln Z_{\mathrm{q}}+ \\ +\frac{1}{2}\mathrm{Tr}\ln (-{\mathbf{G}}_{j;\tau -{\tau }^{\prime }}^{-1}+{\mathbf{V}}_j[{\Psi }_{\tau }]{\mathbf{G}}_{j;\tau -{\tau }^{\prime }}{\mathbf{V}}_j[{\Psi }_{{\tau }^{\prime }}]). \end{gathered}$$

$${\psi }_0\equiv T\underset{0}{\overset{\beta }{\int }}{\psi }_{\tau }d\tau .$$

$${\mathcal{G}}_{j;\tau -{\tau }^{\prime }}^{-1}[{\Psi }_0]\equiv {\mathbf{G}}_{j;\tau -{\tau }^{\prime }}^{-1}-{\mathbf{V}}_j[{\Psi }_0]{\mathbf{G}}_{j;\tau -{\tau }^{\prime }}{\mathbf{V}}_j[{\Psi }_0].$$

$${\mathbf{V}}_j[{\Psi }_0]=\frac{1}{\beta }\underset{0}{\overset{\beta }{\int }}{\mathbf{V}}_j[{\Psi }_{\tau }]d\tau .$$

$${\mathbf{V}}_j[\delta {\Psi }_{\tau }]\equiv {\mathbf{V}}_j[{\Psi }_{\tau }]-\mathbf{V}[{\Psi }_0],$$

$$S_{\mathrm{eff}}[\Psi ]=S_{\mathrm{ph}}[\Psi ]+S_{\mathcal{G}}[{\Psi }_0]+S_{\mathrm{\Sigma }}[\Psi ]+\ln Z_{\mathrm{ph}}.$$

$$S_{\mathcal{G}}[\Phi ]=-\sum _{j=1}^N\ln \frac{\cosh \frac{\sqrt{{ϵ}_j^2+4g_j^2\Phi }}{2T}}{\cosh \frac{{ϵ}_j}{2T}}.$$

$$\begin{gathered} S_{\Sigma }[\Psi ]= \\ \frac{\beta }{2}\sum _{n\ne 0}\left[\begin{array}{cc}{\bar{\psi }}_n& {\psi }_{-n}\end{array}\right]\left[\begin{array}{ccc}{\Sigma }_n[{\Psi }_0]& & {\tilde{\Sigma }}_n[{\Psi }_0]\\ & & \\ ({\tilde{\Sigma }}_{-n}[{\Psi }_0]{)}^{\ast }& & {\tilde{\Sigma }}_{-n}[{\Psi }_0]\end{array}\right]\left[\begin{array}{c}{\psi }_n\\ \\ {\bar{\psi }}_{-n}\end{array}\right]. \end{gathered}$$

$$S_{\mathcal{G}}[\Phi ]\approx \Phi S_{\mathcal{G}}^{\prime }[0]+\frac{1}{2}{\Phi }^2{S}_{\mathcal{G}}^{\prime \prime }[0].$$

$$S_{\mathrm{\Sigma }}[\Psi ]\approx S_{\mathrm{\Sigma }}[\Phi =0,\delta \Psi ].$$

$$S_{\mathrm{\Sigma }}[\Phi =0,\delta \Psi ]=\beta \sum _{n\ne 0}{\Sigma }_n[0]{\bar{\psi }}_n{\psi }_n.$$

$$T\gg \frac{\omega }{N}$$

$${\kappa }_{\mathrm{c}}=\sqrt{\frac{\omega }{NT}}\ll 1$$

$$S_{\mathrm{eff},0}[\Phi ,{\bar{\psi }}_n,{\psi }_n]=S_{\mathrm{zm}}[\Phi ]+S_{\mathrm{fl}}[{\bar{\psi }}_n,{\psi }_n]+\ln Z_{\mathrm{ph}}.$$

Full text

I Introduction

The dynamics of quantum metamaterials Macha et al. (2014); Kakuyanagi et al. (2016); Braumüller et al. (2017); Zagoskin et al. (2016); Lazarides and Tsironis (2018); Jung et al. (2014); Fistul (2017); Shapiro et al. (2015a, b); Greenberg and Moiseev (2019) is a subject of a great interest. These metamaterials are the hybrid systems where cavity photons interact with multi-qubit environment. The behavior of such systems is captured by the Dicke model Emary and Brandes (2003); Brandes (2005); Kirton et al. (2019). The interactions can be characterized by a collective Rabi frequency proportional to a product of the individual qubit-cavity coupling constant and square root of the qubit number. If the Rabi frequency is larger than a certain value then the superradiant phase transition, characterized by an emergence of a large photon number in a cavity and finite order parameter, occurs for temperatures lower than a critical value. The rigorous study of the superradiant phase transition was proposed in the pioneering work of Fedotov and Popov Popov and Fedotov (1988). These authors proposed semi-fermion parametrization of spin operators and described the phase transition in the framework of Matsubara effective action for the photon field. In that work the chemical potential was assumed to be zero and, consequently, the excitations’ number was not constrained. Another case of finite chemical potential in the Dicke model was addressed in Refs. Eastham and Littlewood (2001, 2006) and it was shown that the Bose condensation of polaritons is emerged Popov and Fedotov (1988); Eastham and Littlewood (2006) The Keldysh diagrammatic approach for finite-$N$ corrections, as well as effects of dissipation and external driving, were studied in Refs. Dalla Torre et al. (2013, 2016). Zero temperature description for a limit of large excitations number was obtained in Ref. Pogosov et al. (2017) by means of Bethe-ansatz technique.

Alternatively to the temperature driven transition discussed in Ref. Popov and Fedotov (1988), the superradiance can be turned on by an increase of the interaction strength. It takes place if the Rabi frequency overcomes a critical value. A realization of a control parameter as the interaction energy is possible for quantum metamaterials such as superconducting qubits arrays Macha et al. (2014); Kakuyanagi et al. (2016); Shulga et al. (2017); Zhang et al. (2017) integrated with a GHz transmission line via tunable couplers Srinivasan et al. (2011); Hoffman et al. (2011); Chen et al. (2014); Zeytinoğlu et al. (2015). Also, this may be done in hybrid systems with a controllable amount of nitrogen-vacancy (NV) centers in a diamond sample which interact with an electromagnetic field Dutt et al. (2007); Sandner et al. (2012); Putz et al. (2014); Angerer et al. (2018).

In the present paper we address the situation where the Rabi frequency in quantum metamaterial is varied from weak to ultra-strong coupling domains while the temperature remains constant. We also keep a constant number of qubits $N$ assuming that $N$ is large but finite. It is implied that the loss rate in the cavity is small. The finiteness of $N$ in our consideration means that the superradiant transition is smoothed by the fluctuations of the order parameter and, beside of that, by the thermal fluctuations of polariton quasiparticles. The aim of this work is (i) to describe fluctuations of the above two types and (ii) to formulate a full counting statistics for the photon numbers in this regime.

Our main results are the explicit expressions for the average photons number, its fluctuations and full counting statistics as functions of the collective Rabi frequency. Proposed formalism provides a solution for low temperature $T$ and large $N$ provided they satisfy the conditions $\hbar\omega\gg k_{\rm B}T\gg\hbar\omega/N$ (in this case all qubits are assumed to be in a resonance with the photon mode of the frequency $\omega$). The generalizations for the high-temperature limit, $k_{\rm B}T\gg\hbar\omega$, and dispersive regime, where a spectral density of qubits energies is strongly broadened, are also discussed.

The paper is organized as follows. In the Sec. II we present Matsubara action for the Dicke model, where the qubits degrees of freedom are expressed through the Majorana and complex fermion variables. This is one of possible representations of Pauli operators acting in a Hilbert space of a two-level system. The Sec. III has methodological character. We derive the photon mode’s effective action, which was obtained in previous works Popov and Fedotov (1988); Eastham and Littlewood (2006), by means of the alternative technique with Majorana fermions. In the Sec. IV we present the general expressions for the average photon number and their fluctuations in a resonant limit. In the Sec. V we discuss fluctuational and statistical properties and present a comparison with results of exact numerical simulations at finite temperature and qubit number of the order of ten. In Sec. VI the results are generalized for high temperatures and inhomogeneous broadening in qubit ensemble. In Sec. VII the cumulant generating function for the photon number is derived. In the Sec. VIII we conclude. In the Appendix X we derive the conditions, where our solution based on the Gaussian approximation for thermal fluctuations is strict.

II Path integral formulation

The Dicke Hamiltonian of $N$ qubits reads (we set $\hbar=1$ and $k_{\rm B}=1$ throughout the text):

[TABLE]

Here $\epsilon_{j}$ are the qubits excitation energies, $g_{j}$ are the individual coupling strengths between $j$-th qubit and the photon field in a single-mode cavity. The fundamental frequency of the photon mode is $\omega$. The coupling term is introduced in the standard rotating wave approximation.

In a path integral formulation the photon mode is described by a conventional complex bosonic fields $\bar{\psi},\psi$. The Pauli operators, $\hat{\sigma}^{\pm}_{j},\hat{\sigma}^{z}_{j}$, acting on the $j$-th qubit degrees of freedom, may be represented in path integrals in different ways. It can be bosonic Holstein-Primakoff representation Holstein and Primakoff (1940) or bilinear forms of fermions. Concerning other fermion representations for the Dicke model, techniques based on semi-fermions with an imaginary chemical potential Popov and Fedotov (1988) or auxiliary boson field Eastham and Littlewood (2001) were employed. These representations allow to eliminate the emergent unphysical states and to reduce a Hilbert space to that of a spin-1/2. The semi-fermion representation for spin operators was generalized for Keldysh technique in Ref. Kiselev and Oppermann (2000). Another one fermion representation, which we choose for our calculations, is given by the product of a complex $\hat{c}_{j}\neq\hat{c}^{\dagger}_{j}$ and Majorana $\hat{d}_{j}=\hat{d}_{j}^{\dagger}$ fermion operators Martin (1959); Tsvelik (2007):

[TABLE]

They correspond to three Grassmann fields $\bar{c},c$ and $d$ in a path integral formalism. The use of Majorana fermion allows to avoid auxiliary constraints in the action. Fields $\bar{c},c$ are related to usual complex fermion mode with the excitation energy of two-level system. Field $d$ stands for Majorana zero energy mode with $\langle\hat{d}^{2}\rangle=1/2$. Majorana representation of spin operators has been applied to spin-boson model Schad et al. (2015, 2016) and to a description of spin-spin interaction via helical Luttinger liquid Yevtushenko and Yudson (2018). Recently, this fermionization has been applied to the Dicke model with counter-rotating terms in the interaction Hamiltonian and a regime of quantum chaos has been studied Alavirad and Lavasani (2018). In our studies, which are focused on fluctuation-dominated regime near superradiant phase transition and behavior at finite $N$, Majorana representation appears as a convenient tool.

Below we demonstrate how one can obtain the effective action for photon field with the use of the fermionization (2). The starting point of such consideration is the path integral formulation of the partition function $Z$ in terms of the boson complex fields $\Psi_{\tau}=[\bar{\psi}_{\tau},\psi_{\tau}]$ and fermion fields $\bar{c},c,d$ Kamenev (2011):

[TABLE]

with the action is

[TABLE]

Here $S_{\rm ph}[\Psi]$, $S_{\rm q}[\bar{c},c,d]$ and $S_{\rm int}[\Psi,\bar{c},c,d]$ are the Matsubara actions of the photon mode, qubit environment and their interaction, respectively. The last term $\ln Z_{\rm ph}Z_{\rm q}$ appears due to a normalization of $Z$ to unity at the decoupled limit $g_{j}\to 0$.

Below we consider the terms in (4) in more details. Both of the qubit and photon subsystems are assumed to be in thermal equilibrium at the temperature $T$. The photon mode action, defined on the imaginary time interval $\tau\in[0,\beta]$, where $\beta=1/T$, is

[TABLE]

where the inverse Green function of free photon mode is

[TABLE]

The Fourier transformations from $\tau$ to Matsubara bosonic frequencies $\omega_{n}=2\pi nT$ are defined for the fields and for the Green functions as

[TABLE]

and

[TABLE]

In this representation the photon mode action (5) is transformed into

[TABLE]

The qubit ensemble action is

[TABLE]

The matrix $\mathbf{G}_{j;\tau-\tau^{\prime}}^{-1}$ describes the $j$-th qubit. It contains the inverse Green functions for the $j$th complex fermion and its conjugate with the energies $\pm\epsilon_{j}$, respectively, and the Majorana fermion of zero energy:

[TABLE]

Note, that a corresponding Fourier transformation of the fields $\bar{c}_{\tau},c_{\tau}$ and $d_{\tau}$ and the elements of $\mathbf{G}_{j;\tau-\tau^{\prime}}$ assumes the fermionic frequencies $\omega_{n}=2\pi nT+\pi T$. Bilinear forms $c_{j}d_{j}$ and $\bar{c}_{j}d_{j}$ appear in $S[\Psi,\bar{c},c,d]$ due to the qubit-cavity coupling encoded by the matrix $\mathbf{V}_{j}[\Psi_{\tau}]$:

[TABLE]

This is the matrix which involves the complex boson fields $\psi_{\tau}$, $\bar{\psi}_{\tau}$ as follows:

[TABLE]

The normalization term in (4) is the product of partition functions of non-interacting photon mode and $N$ qubits. The logarithms of their partition functions $Z_{\rm ph}=\int\mathcal{D}[\Psi]\exp(-S_{\rm ph}[\Psi])$ and $Z_{\rm q}=\int\mathcal{D}[\bar{c},c,d]\exp(-S_{\rm q}[\bar{c},c,d])$ are the following:

[TABLE]

and

[TABLE]

The prefactor of $\nicefrac{{1}}{{2}}$ results from the integration over Grassmann variables in the representation (10). The sign “${\rm Tr}$” means the trace taken over the imaginary time variables, or, equivalently, by the Matsubara frequency index $n$; in a case of qubits, an additional trace is taken over the internal $3\times 3$ structure of a matrix $\mathbf{G}_{j}$.

III Effective action

To derive the effective action for the photon field, $S_{\rm eff}[\Psi]$, from the full one $S[\Psi,\bar{c},c,d]$, we start from integration over the the fermion modes $c_{j},\bar{c}_{j}$ and $d_{j}$. As a result, the path integral in the partition function is reduced to $Z=\int D[\Psi]e^{-S_{\rm eff}[\Psi]}$ where the effective action is obtained in the most general form

[TABLE]

Expanding the logarithm in the last term of (16) we obtain that all odd order terms are equal to zero. This follows from the diagonal and non-diagonal structures of $\mathbf{G}_{j}$ and $\mathbf{V}_{j}$, respectively. The resummation back of the non-zero terms of even orders gives the identity:

[TABLE]

A direct first order expansion of the logarithm in the second line of (17) by $\mathbf{V}[\Psi_{\tau}]\mathbf{G}_{\tau-\tau^{\prime}}\mathbf{V}[\Psi_{\tau^{\prime}}]$ provides Gaussian action for all Matsubara modes $\bar{\psi}_{n}$, $\psi_{n}$. As it will be shown in Sec. IV, this expansion results in divergent number of photons at the critical Rabi frequency near the transition into superradiant phase (see Eq. 50). This follows from an infinite occupation of zeroth Matsubara frequency component of the field

[TABLE]

To make correct description of photonic subsystem we should leave $\psi_{0}$ in zero order term of (17) and expand the logarithm by the fluctuations $\delta\psi_{\tau}\equiv\psi_{\tau}-\psi_{0}$. This results in effective regularization of the divergency. Note, that Fourier transformation $\delta\psi_{\tau}$ gives the non-zero Matsubara components $\psi_{n\neq 0}$. The field $\psi_{0}$ is related to the complex amplitude of a superradiant order parameter while $\psi_{n\neq 0}$ are related to thermal fluctuations of polaritonic quasiparticles.

The regularization of the divergency mentioned above assumes a redefinition of the Green function, $\mathbf{G}_{j}\to\mathcal{G}_{j}[\Psi_{0}]$ with $\Psi_{0}=[\bar{\psi}_{0},\psi_{0}]$, as follows:

[TABLE]

Here we introduce the matrix with zero-mode components

[TABLE]

Below we limit our consideration of the fluctuations taking into account bilinear combinations of the fields $\delta\bar{\psi}_{\tau}$ and $\delta\psi_{\tau^{\prime}}$. These are gauge invariant terms $\delta\bar{\psi}_{\tau}\delta\psi_{\tau^{\prime}}$ which provide normal coupling channel between the photons in the dissipative action $S_{\Sigma}$. Oppositely, terms $\delta\bar{\psi}_{\tau}\bar{\delta}\psi_{\tau^{\prime}}$ and $\delta\psi_{\tau}\delta\psi_{\tau^{\prime}}$ are not gauge invariant and provide anomalous type of coupling. At the given step of the derivation we perform the logarithm expansion in (17) around $\mathcal{G}^{-1}$ in second order by the matrix

[TABLE]

which involves the fluctuating parts in $\delta\Psi_{\tau}=[\delta\bar{\psi}_{\tau},\delta\psi_{\tau}]$. We note that the first order contribution by $\mathbf{V}_{j}[\delta\Psi_{\tau}]$ equals zero in this approach. As a result, we obtain

[TABLE]

The first term $S_{\rm ph}$ in (22) is not changed. The second term $S_{\mathcal{G}}[\Psi_{0}]\equiv-\frac{1}{4}\sum_{j}{\rm Tr}\ln\left(\mathbf{G}_{j}\mathcal{G}_{j}^{-1}[\Psi_{0}]\right)$ involves the zero-frequency mode $\Psi_{0}$ only. Note, that in the Dicke model (1) the interaction is limited by the rotating wave approximation which conserves the excitations number. In this case $S_{\mathcal{G}}$ depends on the zero mode’s magnitude squared, $\Phi\equiv\bar{\psi}_{0}\psi_{0}$, and is independent on its complex phase $\varphi\equiv\arg\psi_{0}$. Thus, $S_{\mathcal{G}}[\Psi_{0}]=S_{\mathcal{G}}[\Phi]$ and its explicit expression is

[TABLE]

This result follows from a representation of the Green functions $\mathbf{G}$ and $\mathcal{G}$ in Matsubara frequencies $\omega_{n}$. After that, $S_{\mathcal{G}}$ is reduced to a calculation of infinite product by $n$. The third term in (22) quadratic by quasiparticle fluctuations reads

[TABLE]

This is the dissipative part of the action; it corresponds to effective photon-photon interaction via qubits degrees of freedom. The self-energy operators $\Sigma_{\tau}[\Psi_{0}]$ and $\tilde{\Sigma}_{\tau}[\Psi_{0}]$ provide normal and anomalous channels of the photon-photon interactions, respectively. They result from a summation over the fermionic Matsubara frequencies. From calculations it follows that normal self-energy depends on $\Phi$ only, $\Sigma[\Psi_{0}]=\Sigma[\Phi]$, while the anomalous one depends also on the phase, i.e., $\tilde{\Sigma}[\Psi_{0}]\equiv\tilde{\Sigma}[\Phi,\varphi]$. Their explicit expressions are presented in the Appendix X, see Eqs. (117) and (118). The above results for $S_{\mathcal{G}}$ and $S_{\rm\Sigma}$ are in full correspondence with that derived in Refs. Popov and Fedotov (1988); Eastham and Littlewood (2006) using alternative spin representations.

The action $S_{\rm eff}$ allows to calculate the thermodynamical average value $\langle\Phi\rangle$ which is superradiant order parameter. As shown below, the action indicates the superradiance as a second order phase transition. The quadratic expansion by $\Phi$ in $S_{\mathcal{G}}[\Phi]$ allows to capture this transition (it corresponds to taking into account the non-Gaussian $|\psi_{0}|^{4}$). As a consequence, if the system is in the normal phase or near the phase transition, one can simplify $S_{\mathcal{G}}$ and $S_{\rm\Sigma}$ assuming that the relevant values of $\Phi$ belongs to a certain region near $\Phi=0$. Namely, analytical calculations presented in this work assumes that we apply second order expansion by $\Phi$ in $S_{\mathcal{G}}$ and neglect by non-Gaussian cross terms $\propto\Phi\bar{\psi}_{n}\psi_{n}$ in $S_{\rm\Sigma}$. For the zero-mode part it means

[TABLE]

For $S_{\rm\Sigma}$ part it means one can neglect the dependencies of self-energies on $\Phi$ and assume

[TABLE]

We obtain that this approximation involves the normal coupling only, i.e.,

[TABLE]

It follows from (118) that \tilde{\Sigma}_{n}[\Phi,\varphi]\propto\Phi\ for small $\Phi$ and, hence, the anomalous terms does not appear in (26). Note, that $S_{\rm\Sigma}$ is purely Gaussian in this case because the terms proportional to $\Phi\bar{\psi}_{n}\psi_{n}$ are neglected.

The validity of the approximations (25) and (26) is analyzed in the Appendix X by means of the effective action for the zero Matsubara mode, see Eq. (115). This action is obtained after the Gaussian integration over all non-zero modes $\psi_{n\neq 0}$ in $S_{\rm eff}$ from (22). The linear by $\Phi$ contributions to the self-energies, $\Sigma_{n}[\Phi]\approx\Sigma_{n}[0]+\Phi\Sigma^{\prime}_{n}[0]$ and $\tilde{\Sigma}_{n}[\Phi,\varphi]\propto\Phi$, are investigated as a perturbations for the action (115). It is shown that such perturbations are small and approximations (25) and (26) are strict if the condition for the temperature and qubit number

[TABLE]

holds. It is assumed here that qubits’ and photon mode’s frequencies are of the same order, $\epsilon_{j}\sim\omega$. The condition (28) also provides the range of parameters where one can go beyond the thermodynamic limit and study finite-$N$ effects. The thermodynamic limit, where fluctuations of order parameter are negligible as $1/N$, corresponds to a situation of simultaneous limits $N\to\infty$ and $g\to 0$ with the constraint $\sqrt{N}g={\rm const}$.

As it is also shown in Appendix X, the ratio

[TABLE]

provides the small parameter of this theory near the superradiant transition which allows one to neglect non-Gaussian terms in a controllable way.

To summarize the above, for large enough qubit number dictated by (28), we obtain an effective theory for low temperatures, $\omega\gg T\gg\omega/N$. The high temperature domain corresponds to $T\gg\omega$. Two approximations (26) and (25) yield the effective action $S_{\rm eff,0}$ which provides a description of the normal phase and fluctuational region near the superradiant transition. It is convenient to represent it as

[TABLE]

where the zero-mode terms are collected in

[TABLE]

and that of quasiparticle fluctuations in

[TABLE]

The parameters for a general case are:

[TABLE]

[TABLE]

and

[TABLE]

In the above formulation, the critical point is $A=0$. For $A>0$ the system is in the normal phase and for $A<0$ a superradiant phase with large amount of photons does emerge. In other words, if $A<0$ then $S_{\rm eff,0}$ has a minimum at the stationary point $\Phi=\Phi^{*}$ with

[TABLE]

In terms of the complex photon field $\psi_{0}$ this corresponds to a saddle line which is a circle in the complex plane.

The control parameter of the phase transition is the collective Rabi frequency defined as

[TABLE]

where we denote $\langle\cdot\rangle_{j}$ as the average over the qubit ensemble. The superradiance condition $A<0$ corresponds to the Rabi frequency exceeding a certain critical value $\Omega>\Omega_{\rm c}$. For the homogeneous limit where all qubits have the same energy, $\epsilon_{j}=\bar{\epsilon}$, the saddle point (36) is given by

[TABLE]

The critical Rabi frequency of the phase transition follows from the condition $\Phi^{*}=0$. From (38) one finds that

[TABLE]

We also introduce the action

[TABLE]

where, in contrast to $S_{\rm eff,0}[\Psi]$ in (30), the zero mode’s part (23) is taken into account exactly and its logarithm is not expanded. This action provides an adequate description of the superradiant phase where is $\langle\Phi\rangle$ large. Calculations combine the exact integration by $\bar{\psi}_{n}$ and $\psi_{n}$ and numerical integration by $\Phi$ in this regime. In what follows we focus mainly on the transition between the normal phase and the fluctuational region employing the formalism of $S_{\rm eff,0}[\Psi]$. A behavior in the superradiant phase is briefly discussed below.

IV Photons number and their fluctuations for resonant limit

In this part of the paper we study fluctuational behavior of the superradiance with the use of $S_{\rm eff,0}$ in a limit of full resonance between qubits and photon mode, i.e.,

[TABLE]

The disorder in $g_{j}$, in its turn, is taken into account. The parameters (33, 34) are reduced to

[TABLE]

[TABLE]

We introduced here the collective Rabi frequency renormalized by $T$,

[TABLE]

and the function $f(x)=\frac{\sinh x-x}{1+\cosh x}\coth^{2}\frac{x}{2}$; the parameter $q$ is a ratio between fourth and second moments for coupling parameters, $q=\langle g^{4}\rangle_{j}/\langle g^{2}\rangle^{2}_{j}$. The absence of the disorder in $g_{j}$ corresponds to $q=1$; in disordered case $q>1$; $q=9/5$ for a flat distribution ranging from $g_{\rm min}$ to $g_{\rm max}$ with $g_{\rm max}\gg g_{\rm min}$.

In further consideration the photon number

[TABLE]

is analyzed. Alternatively, it is given by the following identity

[TABLE]

[TABLE]

where $G_{n}$ is $n$-th component of Matsubara Green function. If the quadratic expansion in (17) is applied then one obtains $S_{\rm eff,0}[\Phi,\bar{\psi}_{n},\psi_{n}]$ with $\gamma=0$. This action is fully Gaussian with respect to all Matsubara modes. The following expression for the Green function is obtained for arbitrary $\epsilon_{j}$ and $\omega$ within this expansion:

[TABLE]

For the resonant limit $\epsilon_{j}=\omega$ it reads:

[TABLE]

It is used in the calculations below. This expression holds for any $n$ in the Gaussian approach ($\gamma=0$). After the summation one obtains the average photon number:

[TABLE]

One can see that $\langle N_{\rm ph}\rangle_{\rm Gauss}$ is divergent at the critical value of the renormalized Rabi frequency $\Omega_{T,c}=\omega$ and is negative for $\Omega_{T}>\omega$. This follows from the condition $G_{n=0}=-\frac{1}{\alpha T}$. It is divergent at the critical point where $\alpha=0$. The regularization is provided by the expansion with respect to $\mathcal{G}$ in (17) which involves high order terms by $\Phi$. As we have shown above, quadratic expansion of $S_{\rm eff}$ by $\Phi^{2}$ gives $S_{\rm eff,0}$. Corresponding zero mode’s Green function is changed to $G_{0}=-\frac{\langle\Phi\rangle}{T}$ which is not divergent anymore at the critical point due to $\gamma\neq 0$.

Within the $S_{\rm eff,0}$ action, for non-zero modes the expressions for $G_{n\neq 0}$ are the same as in (49). We refine the definition for the average $\langle N_{\rm ph}\rangle$ as

[TABLE]

The zero mode’ part is written explicitly here. We emphasize thereby that it is calculated within the non-Gaussian (fourth order) approach by $\psi_{0}$. Let us calculate both of the contributions originating from the superradiant order parameter, $\langle\Phi\rangle$, and from the thermal excitations $\langle\bar{\psi}_{n}\psi_{n}\rangle$. As long as there is no explicit dependence on $\varphi$, one has $\iint d(\mathop{\rm Re}\psi_{0})\,d(\mathop{\rm Im}\psi_{0})=\pi\int\limits_{0}^{\infty}d\Phi$. For $\langle\Phi\rangle$ we find

[TABLE]

with the complementary error function is ${\rm erfc}z=1-{\rm erf}z$. Summation over $n\neq 0$ gives the quasiparticle contribution

[TABLE]

Finally, for the average photon number (51) we obtain

[TABLE]

The fluctuations of the photon number are given by the second cumulant $\langle\!\langle N_{\rm ph}^{2}\rangle\!\rangle\equiv\langle N_{\rm ph}^{2}\rangle-\langle N_{\rm ph}\rangle^{2}$. With use of the above notations it is reduced to

[TABLE]

Calculation of the integrals by $\Phi$ and summation over $n$ provide

[TABLE]

We introduced here the second cumulant in Gaussian approximation $\langle\!\langle N_{\rm ph}^{2}\rangle\!\rangle=T^{2}\sum\limits_{n}G_{n}^{2}$. It reads

[TABLE]

Similar to (54), the divergent zero frequency term in the sum is canceled by $1/\alpha^{2}$ in (56).

In the Section V the properties of the photon number and its second cumulant near the phase transition are analyzed in details.

V Phase transition at low temperatures. Resonant limit

V.1 Average photon number near the phase transition

In the following consideration at low temperatures $T\ll\omega$, we should emphasize that there is also a limitation (28) which means that $T$ can not be arbitrary small. Namely, it belongs to the domain

[TABLE]

In such a limit we set $f(\beta\omega)=1$ and $\Omega_{T}=\Omega$ with the exponential by $\beta\omega$ accuracy. In this Section we continue to study the case of full resonance between cavity and all the qubits.

We obtain an analytical expansion of photon number $\langle N_{\rm ph}\rangle$ (54) around the critical point $\Omega_{\rm c}=\omega$. The expansion in series by the dimensionless detuning $\frac{\Omega-\omega}{\omega}$ for $\langle N_{\rm ph}\rangle$ is

[TABLE]

The main contribution to $\langle N_{\rm ph}\rangle$ follows from the $\psi_{0}$ mode as powers of $\sqrt{\frac{NT}{\omega}}$. The prefactors contain the leading term given by zero mode, and small corrections $\delta n$. The corrections follow from the fluctuations of the modes $\psi_{n\neq 0}$. Their expressions might be obtained from the expansion of (53) as

[TABLE]

and

[TABLE]

The zero order term in (59) gives large but finite photon number at the critical point

[TABLE]

The leading term is much higher than unity under the condition $N\gg\omega/T$. If one goes beyond the validity of $S_{\rm eff,0}$ taking a formal limit of $T\to 0$ in (62) then the unphysical value of $-\nicefrac{{1}}{{4}}$ is obtained. This demonstrates that for low temperatures the non-Gaussian fluctuations are needed to be taken into account. It is known that $\langle N_{\rm ph}\rangle=1/2$ at zero temperature limit above the critical point. This is due to the ground state wave function contains $\nicefrac{{1}}{{2}}$ photon on the average. The ground state is changed at the critical point from a direct product of zero photon state and qubits ground state, $|n{=}0;\sigma_{j}{=}-1\rangle$ ($j=1,\ ...\ N$), to an entangled state with a single photon and excited qubits. The field-theoretical approach provided does not allow to describe this limit because it is restricted to finite temperatures $T\gg\omega/N$. Nevertheless, this formalism allows to demonstrate a positive change in the negative constant value in $\langle N_{{\rm ph}}\rangle_{\rm c}$ for very low $T$. For instance, if the third order correction $\propto\Phi\bar{\psi}_{n}\psi_{n}$ is taken into account in the action $S_{\rm eff,0}$. This correction originates from the dependence of the normal part of the self-energy $\Sigma_{n}$ on $\Phi$. In Appendix X it is shown that the integration over all non-zero modes provides a correction $\delta S[\Phi]=\delta\alpha\Phi$ to $S_{\rm zm}[\Phi]$ due to the third-order term. At the critical point and low temperatures $T\ll\omega$ this coefficient reads as $\delta\alpha_{\rm c}=\frac{3\omega}{4TN}$. This results in the positive shift in (62) as $\langle N_{{\rm ph}}\rangle_{\rm c}^{\prime}=\langle N_{{\rm ph}}\rangle_{\rm c}+b$ where $b=\frac{3}{8}(1-2/\pi)\approx 0.1363$.

V.2 Fluctuations and the Fano factor

At the critical point and low temperature limit (58) the fluctuations of photons number are:

[TABLE]

This is the sum of large leading term due to superradiant order parameter fluctuations, $\langle\!\langle\Phi^{2}\rangle\!\rangle$, and small correction $\sim T/\omega\ll 1$ due to weak fluctuations of quasiparticles. The relative value of fluctuations

[TABLE]

is large as $e^{\omega/T}$ in the decoupling limit $\Omega=0$ and decays monotonously due to the decreasing of the second cumulant. It is less than unity above the phase transition. Using the expressions for $\langle N_{\rm ph}\rangle$ and $\langle\!\langle N_{\rm ph}^{2}\rangle\!\rangle$, Eqs. (54) and (56), one obtains the expansion near the phase transition up to the first order by the dimensionless detuning:

[TABLE]

At the critical point ($\Omega=\omega$) the main contribution is due to the zero mode and, consequently, $r_{\rm c}=\frac{\langle\!\langle\Phi^{2}\rangle\!\rangle}{\langle\Phi\rangle^{2}}$. The universal value of the relative fluctuations $r$ at the transition point is

[TABLE]

It follows from the $\Phi$-integrals (52) at $\alpha=0$ and is exact up to the small correction $\sim N^{-1/2}$.

From the expansion (65) the width of the fluctuational Ginzburg-Levanyuk region, $\Omega_{\rm GL}$, near the critical Rabi frequency can be defined. This is a domain where fluctuations and average value of the number of photons are of the same order. This consideration can be applied straightforwardly to the superradiant phase where $\Omega>\Omega_{\rm c}$. The parameter $\Omega_{\rm GL}$ is obtained from the matching conditions

[TABLE]

which give

[TABLE]

Approaching the critical point from the normal phase, i.e., $\Omega<\Omega_{\rm c}$, fluctuations are always greater that average values and the definition (67) is not valid. Instead of (67) we introduce the width $\Omega^{\prime}$ where the superradiant order parameter fluctuations start to grow and become relevant. In this region the contribution to $\langle N_{\rm ph}\rangle$ due to the non-Gaussian fluctuations of $|\psi_{0}|^{4}$ is comparable with the quasiparticle’s part related to $\psi_{n\neq 0}$. We define $\Omega^{\prime}$ through the value of $\Omega=\Omega_{\rm c}-\Omega^{\prime}$ which provides the matching between the average values obtained in the Gaussian and non-Gaussian approaches:

[TABLE]

From (50) and (54) it follows that $\langle N_{\rm ph}\rangle_{\rm Gauss}\sim T/\Omega^{\prime}$ and $\langle N_{\rm ph}\rangle\sim\sqrt{NT/\omega}$. The width $\Omega^{\prime}$ of fluctuation-dominated region appears the same order as in the superradiant phase, i.e.,

[TABLE]

It is rather narrow and is much less than the temperature due to the condition (58).

In order to illustrate the above results we present in Fig. 1 (a, c) the data for $\langle N_{\rm ph}\rangle$ and $\sqrt{\langle\!\langle N_{\rm ph}^{2}\rangle\!\rangle}$ as functions of $\Omega$. We consider the full resonance limit, $\bar{\epsilon}=\omega$, and low temperature regime $T=0.1\ \omega$. The qubits number is $N=100$ on panel (a) and $N=1000$ on panel (c), hence, the constraint (58) is satisfied. Such a qubit number can be realized in contemporary quantum metamaterials. White, thin light blue and light green sectors correspond to normal (N) phase, fluctuational region and superradiant (SR) phase, respectively. The critical point in this low temperature regime is $\Omega_{\rm c}=\omega$. The width of the Ginzburg-Levanyuk fluctuational region is $2\Omega_{\rm GL}\approx 0.063\ \omega$. The red curves are obtained with the use of the action $S_{\rm eff,0}$ and the corresponding analytical results (54) and (56). It is shown that in the normal phase there are exponential dependencies of the photon number and its fluctuations, as follows from linear sectors in the logarithmic scale. Tuning $\Omega$ to the critical value initiates the superradiant transition where photons number is increased rapidly. Further increase of $\Omega$ drives the system into superradiant state. The quadratic expansion for the logarithm in $S_{\mathcal{G}}$, as it should be, does not work well in this phase. A correct description assumes that the use of the action $S_{\rm eff,1}$ from (40). The green curves are the results obtained by means of $S_{\rm eff,1}$ where the integration over $\Phi$ is performed numerically. Dashed parts of the red curves demonstrate the difference between these two approaches. Green curves show that $\langle N_{\rm ph}\rangle$ and its fluctuations in the superradiant phase grow sub-exponentially. It also follows from this plot that in the normal phase the relative fluctuations value $r_{\rm n}>r_{\rm c}$ and $r_{\rm sr}<r_{\rm c}$ in the superradiant phase.

It is also instructive to analyze the Fano factor defined as a ratio between second and first cumulants as

[TABLE]

This is a representative parameter bringing an information about the statistics. The value of $F$ reflects a type of a coherence between the photons: $F=1$ means that they are uncorrelated, $F<1$ and $F>1$ correspond to their negative and positive correlations, respectively. As shown in Fig. 1 (b, d) the dependence $F(\Omega)$ demonstrates rich behavior. The parameters of calculation here are the same as that in the plot (a): $T=0.1\ \omega$, $N=100$ and $\bar{\epsilon}=\omega$. Red and green curves correspond to calculations based on $S_{\rm eff,0}$ and $S_{\rm eff,1}$, respectively. In the decoupling limit $\Omega=0$ the value of the Fano factor is

[TABLE]

In a low temperature limit $F_{0}\approx 1+e^{-\beta\omega}$ which means that photons are weakly correlated. For a finite $\Omega$ there is the entrance into the negative correlations domain where the dependence is non-monotonous with $F_{\rm min}<F<1$. There is a minimum with $F_{\rm min}\approx 0.8$ for an intermediate strength of $\Omega$. It means negative correlation between photons (antibunching effect) in the normal phase due to the interaction between photons through the qubit environment. It is remarkable, that the dependence $F(\Omega)$ in the fluctuational region demonstrates a dramatic change where $F$ grows rapidly and becomes greater than unity. There is a maximum at the critical point which is given by

[TABLE]

The latter means strongly positive coherence between photons near the superradiant transition. The Fano factor shows the decay entering into the superradiant phase if $\Omega$ is further increased. As one can see there is the reentrance to negative correlations with $F_{\rm sr}<1$.

The finite width of fluctuational region and the peak in the Fano factor dependence are finite-size effects. In thermodynamic limit of $N\to\infty$ the Fano factor peak shrinks to a singularity at the critical point. This tendency is seen from a comparison of Figs. 1 (b) and (d) where $N$ is changed by an order.

V.3 Numerical simulation

In Fig. 2 we compare the results obtained in the field-theoretical formalism and that in exact numerical simulations. The qubit number $N=10$ means that the system is beyond from the thermodynamic limit. Despite that $\kappa_{\rm c}$ is not very small compared to unity, $\kappa_{\rm c}\approx\ 0.58$ on Fig. 2 (b), a well quantitative agreement between the numerical and theoretical results is observed. Surprisingly, analytical solution is in a good agreement with numerical calculations even when $\kappa_{\rm c}=1$, as shown in Fig. 2 (b). We represent results for $\langle N_{\rm ph}\rangle$ and $\langle\!\langle N_{\rm ph}^{2}\rangle\!\rangle$ obtained in three different ways. The red dotted curves are obtained with the use of the action $S_{\rm eff,0}$ and analytical expressions (54) and (56). Green dashed curves are derived with the use of $S_{\rm eff,1}$. There is a difference between them in the superradiant phase Blue solid curves represent the results of numerical calculations based on the definitions

[TABLE]

and

[TABLE]

Here the equilibrium density matrix is

[TABLE]

It is block diagonal due to the conservation of total excitations number in the system. This follows from a commutation of the excitations number operator, $\hat{M}=\hat{\psi}^{\dagger}\hat{\psi}+\sum_{j}\hat{\sigma}^{+}_{j}\hat{\sigma}^{-}_{j}$, and $\hat{H}$. In calculations the maximum of excitations number is $M_{\rm max}=50$. This means that $\hat{\rho}$ has $M_{\rm max}$ blocks each of them has the dimension of $2^{M}$, $M=1,\ ...\ ,M_{\rm max}$. For the above parameters the most relevant part of the Fock space belongs to $M$ which covers a region from one to a value around 30. We observe a good correspondence between theoretical curves (red dashed and green dotted) and numerical simulation (blue solid curves) for the range of Rabi frequencies $0<\Omega\lesssim\omega$ which covers the normal phase and fluctuational region. In the superradiant phase, where $\Omega\gtrsim\omega$, the numerical results are in good agreement with more precise calculations based on $S_{\rm eff,1}$.

VI Some generalizations

VI.1 High temperatures

Below we discuss results obtained at the critical point for the high temperature regime $T\gg\omega$. Note that the phase transition at $\alpha=0$ (see Eq. (33)) is given by the increased collective coupling:

[TABLE]

We use (54) and (56) to obtain the leading order expansions for $\langle N_{\rm ph}\rangle$ and $\langle\!\langle N_{\rm ph}^{2}\rangle\!\rangle$ by the large parameter $T/\omega$.

In the Appendix X we discuss that the Gaussian approximation for quasiparticle fluctuations is valid for any $N$ and the corrections due to cross terms $\propto\Phi\bar{\psi}_{n}\psi_{n}$ are always small. This is distinct from $N\gg\omega/T\gg 1$ in the low-temperature limit addressed above.

We obtain that the photon number at the critical point is

[TABLE]

In contrast to the low temperature limit where it scales as $\propto\sqrt{T}$, in the high temperature regime under consideration it grows as $\propto T$. The fluctuations of photons,

[TABLE]

in contrast to (63), contain not only the contributions from $\Phi$ (first term), but also from the non-zero modes $\psi_{n}$ as well (second term). Thus, the high temperature limit is distinct in that sense that there are two domains of $N$ where fluctuations have different contributions. The first domain for $N$ is related to the thermodynamical limit of very large qubit number. It is given by (78) as

[TABLE]

when only superradiant zero mode is relevant. The second one is the intermediate region,

[TABLE]

when contribution of fluctuations of the order parameter can be neglected compared to that of thermal fluctuations of quasiparticles. The relative value at the transition for this intermediate domain,

[TABLE]

shows a deviation from the universal value $\pi/2-1$ due to the second term. Thus, $N\sim\sqrt{T/\omega}$ defines a condition for the crossover between two types of fluctuational behavior. Namely, $N\gg\sqrt{T/\omega}$ corresponds to thermodynamic limit where fluctuations of superradiant order parameter provide the leading contribution to fluctuations of the photons number. In case of $\sqrt{T/\omega}\gtrsim N$ the contribution due to thermal fluctuations of quasiparticles becomes dominant.

VI.2 Inhomogeneous broadening

In the above results for the resonant limit a spread of coupling energies $g_{j}$ yields the prefactor $q^{-1}$ for the qubit number. The inhomogeneous broadening of qubit energies modifies the expressions in a more significant way described below.

We also assume that qubit frequencies are distributed in a certain interval, temperatures are low enough, $T\ll\epsilon_{j}$, and couplings are homogeneous, $g_{j}\equiv g$. We assume that the system is in the critical point, $\alpha=0$, and photons number (54) is contributed by the zero mode only, i.e., $\langle N_{\rm ph}\rangle=\frac{1}{\sqrt{\pi\gamma}}$, and quasiparticles contributions are neglected. In the definition for $\gamma$ (34) the sum over qubit index is replaced by the integral over energies, $\sum_{j}\to N\int\rho({\epsilon})\,d{\epsilon}$ with the density of states $\rho({\epsilon})$ is normalized to unity ($\epsilon$ is a qubit’s energy). We discuss two cases which correspond to flat distributions with finite and very broad widths.

In the first case we consider the distribution with a median energy at $\overline{\epsilon}$ and width $\Delta$, hence, the density of states is

[TABLE]

The photon number is obtained as

[TABLE]

where dimensionless prefactor $z$ is

[TABLE]

In the homogeneous limit, $\Delta\to 0$, this prefactor is unity. Note, that the expression (83) provides the photon number at the critical point for the off-resonant regime, where ${\overline{\epsilon}}\neq\omega$.

In the second case of the very broad distribution, qubits energies belong to the interval from $\epsilon_{\rm min}$ up to large $\epsilon_{\rm c}\gg\epsilon_{\rm min}$ which is spectrum cut-off. This case is considered as a thermodynamic limit where the average level spacing can be introduced, $\delta\epsilon\equiv\epsilon_{\rm c}/N$. Under the assumption

[TABLE]

we find that

[TABLE]

In a physically relevant situation the lower edge of the qubits spectrum $\epsilon_{\rm min}$ may be of the order of the resonator mode frequency, hence, their fraction is order of unity. The logarithm is also not a very large number. Interestingly, in this case we obtain that the photon number is affected mainly by the ratio between the smallest energy scales – the temperature and level spacing.

VI.3 Off-resonant regime

In this subsection we generalize the result for photon number where $\omega$ and $\bar{\epsilon}=\epsilon_{j}$ are out of the resonance. We assume no disorder in $g_{j}$. The value of $\langle N_{\rm ph}\rangle$ is given by the same expression as in Eq. (54) but $\alpha$, $\gamma$ and the Gaussian part are taken in more general form due to $\bar{\epsilon}\neq\omega$. The functional coefficients are

[TABLE]

[TABLE]

The Gaussian part is given by the sum (46) with $G_{n}$ from (48). In the off-resonant case it reads as

[TABLE]

where

[TABLE]

In the resonant limit of $\bar{\epsilon}=\omega$, addressed in the Sec. IV, the expression (89) reproduces (50).

In Fig. 3 (a) and (b) the photon number as the function of $\Omega$ and qubits energies $\bar{\epsilon}$ is plotted (in units of $\omega$). The effective action $S_{\rm eff,1}$ is employed in this calculation. The data shown in (a) demonstrates the behavior at low temperature $T=0.1\,\omega$; (b) demonstrates the behavior at intermediate temperature $T=\omega$. The qubit number $N=100$ in both of the plots. The dark (bright) regions in the maps correspond to normal (superradiant) phases. Red curves depict dependencies of the critical coupling value $\Omega_{\rm c}(\bar{\epsilon})$ from (39) where $\omega$ is kept constant. Curves in insets demonstrate the average photon number as functions of $\Omega/\omega$ for cuts in the plots marked by green dashed lines. Red points in insets stand for the critical Rabi frequency for a given $T$ and cuts of $\bar{\epsilon}$ in (a) and (b). These plots demonstrate typical scales of photons number in normal and superradiant phases for low and intermediate temperature regimes. Red curves corresponding to $\Omega_{\rm c}(\bar{\epsilon})$ relations reproduce asymptotics for low temperatures in (a), where $\Omega_{\rm c}(\bar{\epsilon})\propto\sqrt{\bar{\epsilon}}$, and for high temperatures with $\Omega_{\rm c}(\bar{\epsilon})\propto{\rm const}$ in (b).

VII Full counting statistics

VII.1 Generating action

The effective action for quantum fluctuations (30) allows to we derive the full counting statistics (FCS) for photon numbers. These are cumulant and moment generating functions (CGF and MGF). These are functions of real counting variable $\xi$. In our consideration the generating action is introduced on the imaginary time.

The CGF and MGF are defined as follows through the partition function $Z(\xi)$

[TABLE]

[TABLE]

$\mathcal{T}$-ordering in the imaginary time representation of the path integrals assumes that the photon number, introduced in (45), is defined as

[TABLE]

in a generating term. Alternatively, the generating term can be also represented as a half sum of (93) with $+{\it 0}$ and $-{\it 0}$, which is symmetric under $\mathcal{T}$- and anti-$\mathcal{T}$-ordering. In the Matsubara representation we obtain the generating action in the form of (92) after such a symmetrization. Due to the commutation of photon operators, we include $-\nicefrac{{1}}{{2}}$ in (92).

The photon number moments $\langle N_{\rm ph}^{n}\rangle\equiv\langle(\hat{\psi}^{\dagger}\hat{\psi})^{n}\rangle$ are given by the derivatives

[TABLE]

while the cumulants are defined as

[TABLE]

Path integration in (92) is reduced to the infinite product of Matsubara Green functions involving the counting variable

[TABLE]

The Green function with the counting variable reads as

[TABLE]

Calculation of the integrals and product in (96) yields for the resonant case ($\epsilon_{j}=\bar{\epsilon}=\omega$):

[TABLE]

where the zero mode’s and quasiparticles’ parts are

[TABLE]

and

[TABLE]

With the use of this result for MGF one can obtain the above expressions for the photon number and its fluctuations (54) and (56).

VII.2 FCS at the phase transition

In the thermodynamic limit of large enough $N$, the leading contribution to cumulants is described by that of the zero mode ${\rm MGF}_{0}(\xi)$. Thus, the CGF for the critical point is

[TABLE]

The first six cumulants, which follows from ${\rm CGF}_{0}(\xi)$, are:

[TABLE]

From a numerical calculation it follows that higher cumulants alter their signs, for instance, as it seen from the negativity of the 5th and 6th ones. The non-zero cumulants for $n>2$ is the consequence of that fact that photons’ probability distribution function is half of a Gaussian because of the positively defined variable of integration $\Phi$ in (96).

The Fourier transformation of the MGF provides the probability density to measure $N_{\rm ph}$ photons on average

[TABLE]

Note, that $\mathcal{P}$ is a non-zero function of the continuous variable $N_{\rm ph}$. This is due to that $N_{\rm ph}$ is not an eigenvalue of the Hamiltonian (1). Hence, non-integer values $N_{\rm ph}$ are assumed to be observed as the thermodynamical averages.

As long as the $\psi_{n}$-fluctuations are frozen out if the system is near the critical point and $N$ is large enough, one finds from (92) and (108) that the probability density is identical to the exponent in $Z$ (92) as

[TABLE]

In particular, at the critical point ${\rm MGF}_{0}(\xi)$ from (98) the distribution is

[TABLE]

This is the half of the Gaussian for $N_{\rm ph}>0$, while for unphysical $N_{\rm ph}<0$ it is zero. At the critical point (we assume below that $\Omega_{\rm c}=\omega$), the distribution’s maximum is located at $N_{\rm ph}=0$. In the superradiant phase, the maximum of $\mathcal{P}(N_{\rm ph})$ is shifted to a non-zero value. In other words, for higher values $\Omega\gg\omega$ one obtains from $\ln[{\rm MGF}_{0}(\xi)]$ that in the leading order $\langle N_{\rm ph}\rangle=\frac{N\omega^{2}}{2\Omega^{2}}$ and $\langle\!\langle N_{\rm ph}^{2}\rangle\!\rangle=\frac{NT\omega^{3}}{2\Omega^{4}}$. The higher cumulants are strongly suppressed by the exponent: for instance, the third one is $\langle\!\langle N_{\rm ph}^{3}\rangle\!\rangle\sim e^{-N\frac{\omega}{T}}$.

VII.3 FCS for weak interaction and normal phase

In this part we discuss MGF at the normal phase and weak coupling limit. It is assumed that the system is far away from the fluctuational region, i.e., $\Omega_{T}\ll\omega$ (see Eq. 70). Taking the limit $\gamma\to 0$ in (98) one obtains the MGF for the normal phase of the Dicke model:

[TABLE]

In the decoupled limit, where the Rabi frequency is the smallest scale $\Omega_{T}\ll\{T,\omega\}$, one arrives at the MGF of the free photon mode of the frequency $\omega$

[TABLE]

Note, that it is $2\pi$-periodic function of the counting variable. The discrete Fourier transformation of (111) at the finite interval $[0;2\pi]$ of the single period yields the standard Hibbs distribution probabilities

[TABLE]

Obviously, the infinite integral definition (108) one would obtains delta-peaks in the probability distribution density located at $N_{\rm ph}=n\geq 0$, being the eigenvalues of the free photon mode Hamiltonian, as

[TABLE]

Note that the cumulant generating function for the free mode is

[TABLE]

The cumulants itself are

[TABLE]

One arrives at the mentioned above Fano factor $F_{0}=(1-e^{-\beta\omega})^{-1}$ in (71) and the relative fluctuations parameter $r_{0}=e^{\beta\omega}$.

VIII Conclusions

In this work we addressed to fluctuations near superradiant transition which is driven by an interaction between a single-mode photons and multi-qubit environment. In such consideration the collective Rabi frequency is varied (it can be close to the critical value of superradiant transition), while the temperature $T$ is kept unchanged. We did not assume the thermodynamic limit of infinite qubits number $N$ and consider it as large enough but finite value. Our analysis was focused on two types of competing fluctuations – the thermal one and that of the superradiant order parameter. This regime is opposite to the transition by the temperature studied in Ref. Popov and Fedotov (1988).

We used Majorana fermion representation of qubits’ Pauli operators in order to formulate a path integral approach. Having started from the Dicke Hamiltonian, we demonstrate how one can derive the effective action for the photon mode, obtained by alternative fermionization techniques in Refs. Popov and Fedotov (1988); Eastham and Littlewood (2006). After that we calculated the average photons number and equilibrium fluctuations in terms of the effective action formalism. As a generalization, the full counting statistics, providing higher order cumulants of the photon numbers, was formulated.

Most of the result of this paper address a low temperature regime and a resonance between qubits and photon mode frequency $\omega$. It was shown that the Gaussian approximation for thermal fluctuations is exact and analytical solution can be found, if $\hbar\omega\gg k_{\rm B}T\gg\hbar\omega/N$. In this limit the critical value of the collective Rabi frequency is $\Omega_{\rm c}=\omega$ and the average photon number at this point is $\langle N_{\rm ph}\rangle=\sqrt{Nk_{\rm B}T/(\pi\hbar\omega)}$. The relative fluctuations parameter $r_{\rm c}\equiv\langle\!\langle N_{\rm ph}^{2}\rangle\!\rangle/\langle N_{\rm ph}\rangle^{2}$, where the second cumulant is $\langle\!\langle N_{\rm ph}^{2}\rangle\!\rangle=\langle N_{\rm ph}^{2}\rangle-\langle N_{\rm ph}\rangle^{2}$, is universal at the critical point $r_{\rm c}=\pi/2-1$. A domain near $\Omega_{\rm c}$ in the superradiant phase, where $r$ is not suppressed, corresponds to the fluctuational Ginzburg-Levanyuk region. The width of such frequency range is proportional to $\sqrt{\omega k_{\rm B}T/(\hbar N)}$; this is much smaller than $k_{\rm B}T$ and shrinks at thermodynamic limit. Another characteristic, Fano factor $F\equiv\langle\!\langle N_{\rm ph}^{2}\rangle\!\rangle/\langle N_{\rm ph}\rangle$, decreases from the unity in decoupled limit $\Omega\ll\Omega_{\rm c}$ to a minimum $F_{\rm min}<1$ at $\Omega\lesssim\Omega_{\rm c}$. The latter indicates a negative correlation between photons. The further increase of $\Omega$ up to the critical value results in a significant growth of the Fano factor to a maximum $F_{\rm c}\approx\langle N_{\rm ph}\rangle\gg 1$. This means significantly positive photon-photon correlations at the superradiant transition. There is a reentrance no negative correlations in the superradiant phase as it follows from the decaying of the Fano-factor above the critical point.

As a generalization, for opposite limit of wide spectral distribution of qubits environment we find $\langle N_{\rm ph}\rangle{\sim}\sqrt{k_{\rm B}T/\delta\epsilon}\ln\frac{\epsilon_{\rm c}}{\omega}$, where $\delta\epsilon$ and $\epsilon_{\rm c}$ are the average level spacing and upper cut-off energy of the spectrum, respectively. For high temperatures, $k_{\rm B}T\gg\hbar\omega$, the neglecting of non-Gaussian fluctuations of quasiparticles is valid for any $N$ – in contrast to the low-temperature regime. The finiteness of the qubit number can change a behavior of fluctuations at the critical point. Namely, for $\sqrt{k_{\rm B}T/(\hbar\omega)}\gtrsim N$ the quasiparticle fluctuations become greater than that of superradiant order parameter. This intermediate region shows a non-universal enhancement of $r_{\rm c}$ which reveals a two-level nature of qubits environment.

We believe that the above results can be of an interest in the context of state-of-the-art hybrid systems and quantum metamaterials operating in GHz frequency domain. The coupling constants $g$ in superconducting systems range from MHz to several GHz showing a realization of an ultra-strong coupling regime. The ratios of $g/\omega\sim 0.071$ Bosman et al. (2017), $g/\omega\sim 0.6$ Andersen and Blais (2017); Braumüller et al. (2017) and $g/\omega\sim 0.72-1.34$ Yoshihara et al. (2016) have been demonstrated. Consequently, the critical qubit number $N_{\rm c}=(\omega/g)^{2}$ needed for turning on the superradiant transition can be around $10^{0}$ to $10^{2}$. Another possibility for realization of the phase transition are hybrid systems with NV centers in diamonds. Our estimations are based on Ref. Putz et al. (2014) where individual coupling constant $g\sim 10$ Hz and the number of NV centers $N\sim 10^{12}$. The collective Rabi frequency $\Omega\sim 20$ MHz is two orders less than the critical value $\Omega_{\rm c}\sim 2$ GHz and, according to our consideration, the system is in the normal phase. For the above value of $g$, the number $N$ should be increased by four orders up to the critical $N_{\rm c}\sim 10^{16}$ in order to reach the superradiant phase.

IX Acknowledgments

The research was funded by the Russian Science Foundation under Grant No. 16-12-00095. Authors thank Andrey A. Elistratov for fruitful discussions.

X Appendix

In this Appendix we analyze a role of the leading non-Gaussian correction $\propto\Phi\bar{\psi}_{n}\psi_{n}$, which originates from $\Phi$-dependent self-energies $\Sigma_{n}[\Phi]$ and $\tilde{\Sigma}_{n}[\Phi,\varphi]$. The analysis is based on the low-energy effective action which depends on variables $\Phi$ and $\varphi$. Such an action follows from Gaussian integration in $S_{\rm eff}[\Phi,\varphi,\bar{\psi}_{n},\psi_{n}]$ (22) over all fluctuations related to fields $\bar{\psi}_{n},\psi_{n}$ with $n\neq 0$:

[TABLE]

The correction $\delta S[\Phi,\varphi]$ here is the result of the integration; in the most general form it reads as

[TABLE]

In what follows we analyze a contribution to the action $\delta S[\Phi,\varphi]$ due to the dependency of the self-energies on $\Phi$. The normal and anomalous parts of the self-energies originates from the matrix structure of $S_{\rm fl}$ introduced in (24). Their explicit expressions are:

[TABLE]

[TABLE]

We focus below on the fluctuational region near the superradiance which is second order phase transition. It means that the leading contribution from quantum dynamics of $\Phi$ belongs to a certain region near $\Phi=0$. It follows from (117) and (118) that near $\Phi=0$ the leading $\Phi$-dependent contributions in self-energies are $\Sigma_{n}[\Phi]\approx\Sigma_{n}[0]+\Phi\Sigma^{\prime}_{n}[0]$ and $\tilde{\Sigma}_{n}[\Phi,\varphi]\propto\Phi$. We analyze these $\Phi$-dependent parts as perturbations. It is required that such perturbations give small contributions to unperturbed part of the effective action (115). Similarly to the previous studies, we limit ourselves by the second order expansion in the unperturbed part:

[TABLE]

Note that $\alpha$ can be arbitrary small at the critical point while $\gamma$ is always finite. It means that the linear by $\Phi$ part form $\delta S[\Phi,\varphi]$ is relevant while higher order terms are not. The expansion of the logarithm in $\delta S[\Phi,\varphi]$ up to the first order by $\Phi$ gives

[TABLE]

where

[TABLE]

In other words, the leading perturbation from $\delta S[\Phi,\varphi]$ is contributed by the linear part in the normal self-energy $\Sigma_{n}[\Phi]$. The anomalous part results in second order corrections which provides small corrections to $\gamma$ and are not relevant. The phase $\varphi$ does not appear in this consideration.

We calculate $\delta\alpha$ from (120) for the case of full resonance $\epsilon_{j}=\bar{\epsilon}=\omega$ and absence of disorder in coupling terms, i.e. $g_{j}=g$. At the critical point under consideration we set $\Omega=\omega\sqrt{\coth\frac{\omega}{2T}}$ and the result is

[TABLE]

The expression (121) has the following asymptotics for low temperatures

[TABLE]

Let us estimate typical value of $\Phi^{\prime}$ where the integral over $\exp[-(\alpha+\delta\alpha)\Phi-\gamma\Phi^{2}]$ in the partition function does converge. In the vicinity of the phase transition ($\alpha=0$) it is given by the Gaussian integrand’s width, i.e., $\Phi^{\prime}\sim\gamma^{-1/2}\sim\sqrt{NT/\omega}$. The requirement that the perturbation $\delta S[\Phi]$ is small means that it is much less than the unity at the convergence region, namely, $|\delta S[\Phi^{\prime}]|=|\delta\alpha_{\rm c}\ \Phi^{\prime}|\ll 1$. The latter results in the following condition for low temperatures

[TABLE]

It defines the range of parameters where our approach based on the Gaussian approximation in $S_{\rm fl}$ (32) for fluctuations is valid.

It is important that this calculation provides the small parameter of this theory $\kappa\equiv|\delta S[\Phi^{\prime}]|$. At the critical point of the superradiant transition it is

[TABLE]

Note, that for much lower temperatures $T\ll\omega/N$ the non-Gaussian contributions by $\psi_{n}$ as well as $\Psi_{0}$-dependencies in $\Sigma_{n}$ and $\tilde{\Sigma}_{n}$ can not be neglected in $S_{\rm fl}$. Technically, it means that the higher order terms in the expansion of (17) by $\mathbf{V}[\Psi_{\tau}]\mathbf{G}_{\tau-\tau^{\prime}}\mathbf{V}[\Psi_{\tau^{\prime}}]$ should be taken into account.

In the high temperature limit, $T\gg\omega$, the critical coupling is enhanced as $\Omega_{\rm c}=\sqrt{T\omega}$, the correction and convergence region are

[TABLE]

We obtain the following value of the non-Gaussian correction $|\delta S[\Phi^{\prime}]|=|\delta\alpha_{\rm c}\ \Phi^{\prime}|\sim\beta^{2}\omega^{2}/\sqrt{N}$. The value of $|\delta S[\Phi^{\prime}]|$ is small compared to unity for any $N$, in contrast to low temperature limit (123).