Superradiant phase transition in the ultrastrong coupling regime of the two-photon Dicke model

TL;DR

This paper investigates the two-photon Dicke model, revealing a superradiant phase transition in the ultrastrong coupling regime before spectral collapse, and compares its critical behavior with the one-photon case.

Contribution

It demonstrates the existence of a superradiant phase transition in the two-photon Dicke model and extends mean-field analysis with quantum fluctuations.

Findings

Superradiant phase transition occurs before spectral collapse.

Quantum fluctuations modify the critical behavior.

Comparison with one-photon case highlights unique features.

Abstract

The controllability of current quantum technologies allows to implement spin-boson models where two-photon couplings are the dominating terms of light-matter interaction. In this case, when the coupling strength becomes comparable with the characteristic frequencies, a spectral collapse can take place, i.e. the discrete system spectrum can collapse into a continuous band. Here, we analyze the thermodynamic limit of the two-photon Dicke model, which describes the interaction of an ensemble of qubits with a single bosonic mode. We find that there exists a parameter regime where two-photon interactions induce a superradiant phase transition, before the spectral collapse occurs. Furthermore, we extend the mean-field analysis by considering second-order quantum fluctuations terms, in order to analyze the low-energy spectrum and compare the critical behavior with the one-photon case.

Tables

Table 1. Table 1: Critical exponents of the excitation energy E e ​ x ​ c subscript 𝐸 𝑒 𝑥 𝑐 E_{exc} and of the quadratures X ^ a = a + a ^ † subscript ^ 𝑋 𝑎 𝑎 superscript ^ 𝑎 † \hat{X}_{a}=a+\hat{a}^{\dagger} and X ^ d = d + d ^ † subscript ^ 𝑋 𝑑 𝑑 superscript ^ 𝑑 † \hat{X}_{d}=d+\hat{d}^{\dagger} . The one-photon interaction term is taken to be g N ​ ( a ^ + a ^ † ) ​ ( J ^ + + J ^ − ) 𝑔 𝑁 ^ 𝑎 superscript ^ 𝑎 † subscript ^ 𝐽 subscript ^ 𝐽 \frac{g}{N}(\hat{a}+\hat{a}^{\dagger})(\hat{J}_{+}+\hat{J}_{-}) . For the two-photon case, we took β > 0 𝛽 0 \beta>0 .

	Two-photon case	One-photon case
$E_{exc}$	$\frac{1}{2}$	$\frac{1}{2}$
$\Delta {\hat{X}}_d$	$-\frac{1}{4}$	$-\frac{1}{4}$
$\Delta {\hat{X}}_a$	$0$	$-\frac{1}{4}$

Equations

$$\hat{H}=\omega {\hat{a}}^{†}\hat{a}+\frac{{\omega }_q}{2}\sum _{j=1}^N{{\hat{\sigma }}_j}^z+\frac{g}{N}\sum _{j=1}^N{{\hat{\sigma }}_j}^x({\hat{a}}^2+{\hat{a}}^{†2}),$$

$$\hat{H}=\omega {\hat{a}}^{†}\hat{a}+{\omega }_q{\hat{J}}_z+\frac{g}{N}({\hat{J}}_{+}+{\hat{J}}_{-})({\hat{a}}^2+{\hat{a}}^{†2}).$$

$${\hat{J}}_{+}={\hat{b}}^{†}\sqrt[]{N-{\hat{b}}^{†}\hat{b}},{\hat{J}}_{-}=\sqrt[]{N-{\hat{b}}^{†}\hat{b}}\hat{b},{\hat{J}}_z={\hat{b}}^{†}\hat{b}-\frac{N}{2},$$

$$\hat{b}=\beta +\hat{d},\beta =⟨GS|\hat{b}|GS⟩,[\hat{d},{\hat{d}}^{†}]=1.$$

$$\hat{H}=\omega {\hat{a}}^{†}\hat{a}+g_{\beta }({\hat{a}}^2+{\hat{a}}^{†2})+{\omega }_q|\beta {|}^2-\frac{{\omega }_{q}N}{2},$$

$$g_{\beta }=\frac{g}{N}(\beta +{\beta }^{\ast })\sqrt[]{N-|\beta {|}^2}.$$

$$r_a^{(MF)}=\frac{1}{2}\mathrm{arctanh}(\frac{2g_{\beta }}{\omega }).$$

$$E_G=\frac{1}{2}\frac{\cosh (2r_a^{(MF)})}{\omega }({\omega }^2-4g_{\beta }^2)+{\omega }_q|\beta {|}^2-{\omega }_q\frac{N}{2}-\frac{\omega }{2}.$$

$$\beta ={\beta }^{\ast }=\pm \sqrt{\frac{N}{2}}{\left(1-\sqrt{\frac{1-\mu }{4{\mu }^2{\lambda }^2-\mu }}\right)}^{1/2}=\pm {\beta }_0,$$

$$\lambda =\frac{\omega }{2{\omega }_{q}N}\ge 0;\mu =\frac{4g^2}{{\omega }^2}\ge 0.$$

$$\hat{H}={\hat{H}}_0+ϵ{\hat{H}}_c,$$

$${\hat{H}}_{JC}={\omega }_0{\hat{a}}^{†}\hat{a}+\frac{{\Omega }_{spin}}{2}{\hat{\sigma }}_z-ϵg(\hat{a}{\hat{\sigma }}_{+}+{\hat{a}}^{†}{\hat{\sigma }}_{-}).$$

$${\hat{K}}_0=\frac{1}{2}({\hat{a}}^{†}\hat{a}+\frac{1}{2}),$$

$${\hat{K}}_{+}=\frac{1}{2}{\hat{a}}^{†2},$$

$${\hat{K}}_{-}=\frac{1}{2}{\hat{a}}^2,$$

$$\hat{H}=-\frac{{\omega }_{q}N}{2}+{\omega }_q{\hat{d}}^{†}\hat{d}-\frac{g^2}{N\omega }(\hat{d}+{\hat{d}}^{†}{)}^2+O\left(\frac{\omega }{N\sqrt{N}}\right),$$

$$r_s^{(2)}=-\frac{1}{4}\ln \left(\frac{1+\frac{{\alpha }^2}{1-2{\alpha }^2}\left(3+\frac{{\alpha }^2}{1-{\alpha }^2}\right)-\frac{1-2{\alpha }^2}{1-{\alpha }^2-(16g^2/{\omega }^2){\alpha }^2(1-{\alpha }^2{)}^2}}{1+\frac{{\alpha }^2}{1-2{\alpha }^2}}\right),$$

$$A(g\to g_t)\propto |\frac{g-g_t}{g_t}{|}^{{\gamma }_A},$$

$$\hat{H}=\omega {\hat{a}}^{†}\hat{a}+{\omega }_q{\hat{J}}_z+\frac{g_2}{N}{\hat{J}}_x({\hat{a}}^2+{\hat{a}}^{†2})+\frac{g_1}{N}{\hat{J}}_x(\hat{a}+{\hat{a}}^{†}).$$

$$\hat{H}=\omega {\hat{a}}^{†}\hat{a}+g_2^{\beta }({\hat{a}}^2+{\hat{a}}^{†2})+g_1^{\beta }(\hat{a}+{\hat{a}}^{†})+{\omega }_q|\beta {|}^2-\frac{{\omega }_{q}N}{2},$$

$$g_1^{\beta }=\frac{g_1}{N}(\beta +{\beta }^{\ast })\sqrt[]{N-|\beta {|}^2},$$

$$g_2^{\beta }=\frac{g_2}{N}(\beta +{\beta }^{\ast })\sqrt[]{N-|\beta {|}^2}.$$

$$\hat{H}=\omega {\hat{c}}^{†}\hat{c}+g_2^{\beta }({\hat{c}}^2+{\hat{c}}^{†2})+{\omega }_q|\beta {|}^2-\frac{{\omega }_{q}N}{2}-\frac{(g_1^{\beta }{)}^2}{\omega +2g_2^{\beta }},$$

$${\hat{H}}^{\prime }=e^{-\hat{S}}\hat{H}{e}^{\hat{S}}=\sum _{k=0}^{\infty }\frac{1}{k!}[\hat{H},\hat{S}{]}^{(k)},$$

$${\hat{H}}_d^{\prime }=\sum _{k=0}^{\infty }\frac{1}{2k!}[{\hat{H}}_0,\hat{S}{]}^{(2k)}+\sum _{k=0}^{\infty }\frac{1}{2k+1!}[ϵ{\hat{H}}_c,\hat{S}{]}^{(2k+1)},$$

$${\hat{H}}_a^{\prime }=\sum _{k=0}^{\infty }\frac{1}{2k+1!}[{\hat{H}}_0,S{]}^{(2k+1)}+\sum _{k=0}^{\infty }\frac{1}{2k!}[ϵ{\hat{H}}_c,\hat{S}{]}^{(2k)}.$$

$$\hat{F}(j)+\hat{F}(j-2)+\hat{F}(j-4)+...=\sum _{p<\frac{j}{2}}\hat{F}(j-2p),$$

$$\hat{F}(i)=\sum _{n=0}^{\infty }{\rho }_n^i|n{⟩}_K⟨n+i{|}_K+\sum _{n=i}^{\infty }{\chi }_n^i|n{⟩}_K⟨n-i{|}_K.$$

$$[M_{(i)},M_{(j)}{]}_{op}\subseteq M_{(i+j)}.$$

$$\hat{H}=-\frac{1}{2}\omega -\frac{{\omega }_{q}N}{2}+2\omega {\hat{K}}_0+{\omega }_q{\hat{d}}^{†}\hat{d}+\frac{2g}{\sqrt{N}}({\hat{K}}_{+}+{\hat{K}}_{-})(\hat{d}+{\hat{d}}^{†}).$$

$${\hat{H}}_1$$

$${\hat{H}}_1={\hat{H}}_1^{(0)}+{\hat{H}}_1^{(1)},$$

$${\hat{H}}_1^{(0)}={\hat{K}}_0+\frac{{\omega }_1}{N}{\hat{d}}^{†}\hat{d},$$

$${\hat{H}}_1^{(1)}=\frac{{\omega }_2}{\sqrt{N}}(\hat{d}+{\hat{d}}^{†})({\hat{K}}_{+}+{\hat{K}}_{-})=\frac{1}{\sqrt{N}}{\hat{V}}^{(1)},$$

$$\hat{S}=\frac{1}{\sqrt{N}}\hat{P}+\frac{1}{N}\hat{Q}+\frac{1}{N}\hat{R}+O\left(\frac{1}{N\sqrt{N}}\right),$$

Full text

††preprint: APS/123-QED

I Introduction

In quantum optics, the superradiant phase transition Mallory69 ; Hepp75 (SPT) is the abrupt change in the behavior of the ground state properties of a quantum many-body system, while a physical parameter is continuously varied. It is a quantum phase transition, i.e. it can be accessed at zero temperature and it is due to quantum fluctuations. The archetypal model known to display such a transition is the Dicke model Gelfand54 ; Brandes05 (DM), which describes the interaction of an ensemble of two-level quantum systems, or qubits, with a single bosonic mode. In the limit of a large number of qubits, the DM undergoes a SPT Emary03 in the ultrastrong coupling (USC) regime, where the collective light-matter coupling becomes comparable to the qubit and field bare frequencies Ciuti05 . Although the DM is commonly used to describe atomic and solid-state systems, whether it provides a reliable description of the system ground state when approaching the critical coupling is still the subject of debate Rzazewski75 ; Knight78 ; Gawedzki81 ; Keeling07 ; Deliberato14 ; Vukics14 ; Griesser16 . In particular, the presence of the so-called diamagnatic term is expected to prevent the SPT. The debate has been recently extended to the framework of circuit QED Chen07 ; Lambert09 ; Nataf10 ; Viehmann11 ; Jaako16 ; Bamba16 , where the USC regime has been experimentally achieved Niemczyk10 ; Forn-Diaz10 ; Forn-Diaz16 ; Chen16 ; Yoshihara16 .

However, a compelling way to circumvent no-go theorems consists in using driven systems to engineer effective Hamiltonians. Indeed, the SPT has been observed in driven atomic systems which effectively reproduce the DM Baumann10 ; Baden14 ; Klinder15 . In general, driven atomic or solid state systems represent a powerful tool to access the USC regime of quantum optical models, both in few Crespi12 ; Ballester12 ; Pedernales2015 ; Langford16 ; Felicetti16 and many-body physics Dimer07 ; Carusotto13 . In the USC regime, even apparently simple models entail a very complex physics. This is the case for the quantum Rabi model Braak11 ; Rossatto16 , which corresponds to a single-qubit DM. Furthermore, the engineering of effective Hamiltonians allows to implement generalized quantum optical models Genway14 ; Zou14 , including anisotropic couplings or two-photon interactions. In the case of anisotropic couplings, reaching the USC regime leads to parity-symmetry breaking Xie14 ; Cui15 and to a rich phase diagram in the many-body limit Baksic14 .

Similarly, the two-photon Rabi model has highly counter-intuitive spectral and dynamical features. Its spectrum collapses into a continuous band for a specific value of the coupling strength Emary02 ; Dolya09 ; Travenec12 ; Peng16 . In the transition from the strong to the USC regime of the two-photon Rabi model, a continuous symmetry breaks down into a four-folded discrete symmetry, identified by a generalized-parity operator Felicetti15 . However, so far there are no known results on the ground state of two-photon models in the many-body limit.

In this work, we perform first a mean-field analysis of the two-photon Dicke model and we find that the system exhibits a phase transition in the thermodynamic limit. This transition is superradiant in the sense that it is characterized by a macroscopic change in the average photon number. The boundary of the phase transition is set by the critical value of the coupling strength, which depends on the qubit and field energies. Interestingly, the coupling strength for which the spectrum collapses depends only on the field frequency. For larger values of the coupling strength, the Hamiltonian is not bounded from below and the model is not well defined. We define the parameter regime where the SPT could be accessed within the validity region of the model, that is, where the critical coupling is smaller than the collapse coupling strength. Finally, we go beyond mean-field by including second-order quantum fluctuations. This lets us characterize the system phases and analyze the differences with the SPT of the standard DM. We find fundamental differences in the critical scaling of the bosonic field.

The two-photon DM could be implemented using trapped ions Felicetti15 , which have been used to realize spin systems composed of hundreds of qubits Britton12 . Similar schemes can be conceived for other atomic or solid state systems. Particularly promising are superconducting devices, where bosonic modes have been coupled to increasingly large spin ensembles Macha14 ; Kakuyanagi16 . In any implementation, the critical issue would be the number of qubits that can be effectively coupled to a single bosonic mode. In our case, we have considered the thermodynamic limit $N\gg 1$, with $N$ the number of qubits. We show that, in this limit, the superradiant phase transition of the two-photon DM entails squeezing and spin-squeezing properties. These features provide a signature of the phase transition that could be observed Sorensen99 ; Vitagliano14 ; Saideh16 for a smaller number of qubits.

II Mean-Field Analysis

We consider an ensemble of N qubits interacting with a bosonic mode via two-photon interaction, as sketched in Fig.1. The system Hamiltonian is given by:

[TABLE]

with $\hbar=1$, ${\hat{\sigma}_{j}}^{z}$ and ${\hat{\sigma}_{j}}^{x}$ are Pauli operators describing the j-th ion, and $\hat{a}$ and $\hat{a}^{\dagger}$ are bosonic ladder operators. In Ref.Felicetti15 , Felicetti et al. have proposed to implement this model using a chain of ions in a trap illuminated by two lasers, with the motional degree of freedom of the chain playing the role of the bosonic field. The coupling strength $g$, qubit energy spacing $\omega_{q}$ and bosonic frequency $\omega$ are then effective and tunable parameters which depend on the frequencies and amplitudes of the two lasers. In principle, it is then be possible to reach the ultrastrong coupling (USC) regime: $g\sim\omega$. In the following, we will consider this regime of parameters, as well as the thermodynamic limit $N\gg 1$. In the USC regime, the two-photon Dicke model exhibits a spectral collapse Travenec12 ; Felicetti15 : for $g=\frac{\omega}{2}$, the energy levels of the system collapse into a continuum. Beyond this limit, the ground state of (1) is no longer defined, which renders the notion of phase transition meaningless. Thus, the goal of this work is to study the existence of a phase transition for $g<\frac{\omega}{2}$. For this purpose, it is convenient to describe the ensemble of qubits by pseudospin operators: we define $\hat{J}_{z}=\frac{1}{2}\sum_{j}{\hat{\sigma}_{j}}^{z}$, $\hat{J}_{x}=\frac{1}{2}\sum_{j}{\hat{\sigma}_{j}}^{x}$, $\hat{J}_{\pm}=\sum_{j}{\hat{\sigma}_{j}}^{\pm}$ , which gives us:

[TABLE]

In the following, we will speak only in terms of the fluctuations and polarization of this "spin", which are physically interpreted in terms of population and coherence of the states of the qubits. We study the phase diagram of our model by using a mean-field approach, inspired by the analysis performed in Emary03 and Baksic14 . The point here is to determine how the properties of the ground state evolve when $g$ increases. First of all, we use the Holstein-Primakoff (HP) transformation to turn our pseudospin operators into bosonic operators:

[TABLE]

with $[\hat{b},\hat{b}^{\dagger}]=1$. Here and in the following, we restrict ourselves to the eigenspace of $\overrightarrow{J}^{2}$ associated with eigenvalue $\frac{N}{2}(\frac{N}{2}+1)$ (that is, the maximal angular momentum eigenspace $\mathcal{H}_{\frac{N}{2}}$) note subspace .

Next, we shift $\hat{b}$ and $\hat{b}^{\dagger}$ with respect to their mean values in the ground state $|GS\rangle$ of (2):

[TABLE]

As a zeroth order approximation, we neglect altogether the spin fluctuations $\hat{d}$ and $\hat{d}^{\dagger}$, which gives:

[TABLE]

In other terms, we replace the spin operators by their classical mean values. This Hamiltonian is quadratic in $\hat{a}$, and can therefore be diagonalized by Bogoliubov transformation. The ground state is a squeezed vacuum state with squeezing parameter:

[TABLE]

The corresponding ground state energy is given by:

[TABLE]

The final step consists in minimizing this energy in order to determine the value of $\beta$ the system adopts in its ground state. $\beta$ then plays the role of order parameter for our system: a change in its behavior leads to a modification of the qualitative properties of the ground state, which indicates a phase transition. We find that, for $g<g_{t}=\sqrt{\frac{\omega\omega_{q}N}{4}}$, $E_{G}$ is minimal for $\beta=\beta^{\ast}=0$. For $g>g_{t}$, we have two degenerate minima:

[TABLE]

where we have defined:

[TABLE]

Thus, our system exhibits two phases: in the first, $\beta$ and the squeezing parameter are zero for all $g$, meaning that the field $\hat{a}$ is in the vacuum state and the pseudospin $\hat{J}_{z}$ is polarized along the $z$-axis: $\langle\hat{J}_{z}\rangle=\lvert\beta\rvert^{2}-\frac{N}{2}=-\frac{N}{2}$, $\langle\hat{J}_{x}\rangle=\langle\hat{J}_{y}\rangle=0$. In the second phase, the ground state is twice degenerate. This degeneracy comes from the fact that $\beta$ can be positive or negative. The field $\hat{a}$ is in a squeezed vacuum state, the direction of squeezing depending on the sign of $\beta$: (i) For $\beta>0$, the quadrature $\hat{X}_{a}=\hat{a}+\hat{a}^{\dagger}$ is squeezed, the fluctuations of $\hat{P}_{a}=i(\hat{a}^{\dagger}-\hat{a})$ are amplified. (ii) For $\beta<0$, the squeezed quadrature is $\hat{P}_{a}$. In both cases, when $g$ increases, the squeezing parameter increases and the pseudospin polarization evolves towards the $x$-axis: $\lvert\langle\hat{J}_{z}\rangle\rvert$ decreases until it hits zero for $g=g_{collapse}=\frac{\omega}{2}$, and $\lvert\langle\hat{J}_{x}\rangle\rvert$ increases at the same time. This phase diagram is reminiscent of the one-photon Dicke model and its superradiant phase transition Emary03 ; Baksic14 . Here, only the pseudospin acquires macroscopic mean value in the second phase; the mean value of $\hat{a}$ remains zero. However, the average number of photons becomes nonzero at the transition; thus, we argue that our transition may still be dubbed "superradiant". Note also that the zero mean value of the bosonic field comes from the fact that we have considered only quadratic terms for $\hat{a}$ in the Hamiltonian. If a linear term is present, $\hat{a}$ will acquire nonzero mean value at the mean-field level (See Appendix A).

Let us finally notice that, for the phase transition to occur before the spectral collapse, the following condition must be satisfied: $\omega_{q}N<\omega$. Since $\omega$ and $\omega_{q}$ are adjustable effective parameters, it is possible in principle to meet this condition. From this point on, we will consider that the order of magnitude of the parameter $\omega_{q}$ is $\frac{\omega}{N}$ ($\omega_{q}=O(\frac{\omega}{N})$). We can notice that the approximation that we have made to obtain (5b) amounts simply to keeping only terms of order $O(\omega)$ in the Hamiltonian, and to neglecting higher-order terms note zero omegaq .

III Beyond Mean-Field

Let us now take the fluctuations of $\hat{J}_{z}$ into account, in addition to the mean value. We will expand the Hamiltonian keeping terms of order $O(\omega)$, $O(\frac{\omega}{\sqrt{N}})$ and $O(\frac{\omega}{N})$. Then, we will make use of a technique that Hwang & al. applied to the Rabi and Jaynes-Cumming models article Hwang JC ; article Hwang Rabi . This method is inspired by the Schrieffer-Wolff transformation ref Schrieffer-Wolff ; ref adiabatic elimination , and can be summarized as follows: one starts with an Hamiltonian describing a spin operator $\hat{\sigma}_{z}$ and another degree of freedom, that can be written as

[TABLE]

with a small parameter $\epsilon\ll 1$, $\hat{H}_{0}$ an operator that does not couple the $\hat{\sigma}_{z}$ eigenstates, and $\hat{H}_{c}$ that does. In article Hwang JC for instance, this method was applied to a Jaynes-Cummings Hamiltonian:

[TABLE]

The method consists in decoupling the spin eigenspaces up to a certain order of $\epsilon$. This is done by finding a transformation $e^{\hat{S}}$, with $\hat{S}$ an anti-hermitian operator such that $\hat{H}^{\prime}=e^{-\hat{S}}\hat{H}e^{\hat{S}}$ commutes with $\hat{\sigma}_{z}$ up to a certain order in $\epsilon$. It is then possible to project $\hat{H}^{\prime}$ in one of the spin eigenspaces, which allows to effectively suppress the spin degree of freedom and to diagonalize the Hamiltonian more easily.

In our case, it would seem natural to decouple the eigenspaces of the pseudospin $\hat{J}_{z}$. Instead, however, we are going to apply the HP transformation once more, and shift the operators $\hat{b}$ with respect to their mean-field expectation value $\beta$.

We define the following operators:

[TABLE]

which obey spin-like commutation relations $[\hat{K}_{0},\hat{K}_{\pm}]=\pm\hat{K}_{\pm}$, and $[\hat{K}_{+},\hat{K}_{-}]=-2\hat{K}_{0}$. Note that the commutation algebra here is SU(1,1) instead of SU(2), hence these are not spin operators even if the commutation relations are the same. In order to apply the method described above, we take profit of these commutation relations to decouple the eigenspaces of these pseudospin operators up to a certain order of the small parameter $\frac{1}{\sqrt{N}}$. Then, we project out the $\hat{K}$ degree of freedom. This gives us an effective Hamiltonian describing the low-energy fluctuations of $\hat{b}$ above the ground state; the detailed calculations can be found in Appendix B, C and D.

While not very intuitive at first glance, this manipulation can be justified by two arguments. Firstly, it is interesting to use the HP transformation on $\hat{J}_{z}$ again in order to make a link to our previous study and to re-use our results. Hence, only $\hat{a}$ is available to play the role of the pseudospin whose eigenspaces we seek to decouple. Next, as we mentioned previously, we consider the following regime of parameters: $\omega_{q}=O(\frac{\omega}{N})$. Since $N\gg 1$, this means that the typical energy separation between the $\hat{a}$ eigenstates, given by $\omega$, is much bigger than the typical energy separation for the $\hat{J}_{z}$ (or $\hat{b}^{\dagger}\hat{b}$) eigenstates, given by $\omega_{q}$. As a consequence, the eigenspaces associated with the pseudospin (12) are well separated in energy, but the eigenspaces of $\hat{J}_{z}$ are not. A schematic representation of the energy levels is shown in Fig.2. Thus, it is possible to project the Hamiltonian into the lowest-energy eigenspace of $\hat{K}_{0}$ and to study the fluctuations of $\hat{J}_{z}$ while staying inside this subspace, but not the opposite.

In the first phase, after having decoupled the $\hat{K}_{0}$ eigenspaces and projected the system into the lowest-energy one, we obtain the following Hamiltonian:

[TABLE]

where $\hat{d}$ has been defined in Eq.(4). The associated ground state for $\hat{d}$ is a squeezed vacuum state, with squeezing parameter $r_{s}^{(1)}=\frac{1}{4}\ln(1-\frac{4g^{2}}{N\omega\omega_{q}})$ (let us note that this parameter is negative, meaning that the squeezed quadrature is $\hat{P}_{d}$ instead of $\hat{X}_{d}$). This constitutes a piece of information about the fluctuations of $\hat{J}_{z}$, while only mean values were accessible in our first analysis. The $\hat{a}$ field is found to be in a coherent vacuum state, with no modifications with respect to the mean-field analysis.

In the second phase, we also find squeezing properties for the $\hat{b}$ ground state, with squeezing parameter:

[TABLE]

with $\alpha=\frac{\beta}{\sqrt{N}}$. The results for the $\hat{a}$ field are identical to the mean-field case with a slight correction, i.e. we have a squeezed vacuum state with squeezing parameter $r_{a}^{(2)}$ which differs from $r_{a}^{(MF)}$ by a correction of order $\frac{1}{N}$ (see Appendix D). Going back to the definition of the $\hat{b}$ field (3), the behavior for the $\overrightarrow{J}$ operators can be summarized as follows: in the first phase, $\overrightarrow{J}$ is polarized along the $z$-axis: $\langle\hat{J}_{x}\rangle=\langle\hat{J}_{y}\rangle=0$. When $g$ increases, the fluctuations of $\hat{J}_{y}$ are damped, while the fluctuations of $\hat{J}_{x}$ are amplified proportionally. The amplification factor goes to infinity when approaching the transition: even though the approximations we have used break down near the critical point, the divergence of the $\hat{J}_{x}$ fluctuations in our model may have observable consequences in an actual experiment. In the second phase, the $\overrightarrow{J}$ polarization will gradually evolve from the $z$-axis to the $x$-axis as $g$ increases. Near the spectral collapse $g=g_{collapse}=\frac{\omega}{2}$, $\overrightarrow{J}$ will be polarized along the $x$-axis, and we will have squeezing properties for the fluctuations in the $y$ and $z$ directions. As of the $\hat{a}$ field, we have a coherent vacuum state in the first phase; then, at the transition, the field acquires squeezing properties. In the second phase, the squeezing parameter increases with $g$ and diverges at the spectral collapse. The behavior of spin squeezing fluctuations is schematically depicted in Fig.3. We can also characterize the ground state energy and the $\hat{b}$ excitation energy in both phases: the results for the excitation energy are displayed on Fig.4.

Finally, using the effective Hamiltonian in both phases, we can compute the critical exponents of several observables that exhibit critical behavior at the transition, that is,

[TABLE]

$\gamma_{A}$ being the critical exponent of $A$. We can compare those results to the critical exponents in the one-photon case Emary03 ; Vidal06 , with an interacting term of the form $\frac{g}{N}(\hat{a}+\hat{a}^{\dagger})(\hat{J}_{+}+\hat{J}_{-})$ (note that this convention differs slightly from the one usually found in the literature, $\frac{g}{\sqrt{N}}(\hat{a}+\hat{a}^{\dagger})(\hat{J}_{+}+\hat{J}_{-})$). In the one-photon model, the atomic and photonic excitations are hybridized to form polaritons. The energy spectrum exhibits two polaritonic branches, one of which goes to zero at the transition. In the limit $\omega_{q}\ll\omega$, however, the excitations cease to hybridize; the lowest- and highest-energy polaritons are purely atomic and photonic, respectively, which corresponds to our low-energy $\overrightarrow{J}$ excitations and high-energy $\hat{a}$ ($\hat{K}_{0}$) excitations. Notice that this situation is similar to what is considered in Ref.Mottl12 , where an atomic ensemble interacts with a strongly detuned optical cavity. We compare the critical exponents of the one- and two-photon Dicke models in the regime $\omega_{q}\ll\omega$, in a such a way that the polaritonic branches have the same atomic and photonic components in both cases. The results (computed for $\beta>0$) are summarized in Table 1.

Notice that the variance of $\hat{X}_{a}$ remains constant near the transition in the two-photon case, but diverges in the one-photon case, which is a marked difference. However, the diverging terms in the one-photon case are high-order terms; near the transition, $\Delta\hat{X}_{a}\sim 1+O(\frac{1}{\sqrt{N}}\lvert\frac{g-g_{t}}{g_{t}}\rvert^{-\frac{1}{4}})$. Thus, as long as one does not get too close to the transition, the behavior of $\hat{X}_{a}$ remains identical in the one-photon and the two-photon cases. When $\frac{g-g_{t}}{g_{t}}=O(N^{-2})$, though, the results begin to differ.

On the other hand, the critical exponents are the same in both cases for $E_{exc}$ and $\hat{X}_{d}$. Nevertheless, for these quantities, some of the higher-order terms we have neglected in our analysis may diverge faster near the transition that the terms we took into account, leading to the breakdown of our analysis when $g$ comes close enough to the critical value $g_{t}$. More precisely, the scaling parameter for $\Delta\hat{X}_{d}$ is no longer valid when $\frac{g-g_{t}}{g_{t}}$ becomes comparable with $\frac{1}{N}$.

IV Summary and Outlook

As a conclusion, we have shown the presence of a superradiant phase transition of the two-photon Dicke model in the ultrastrong coupling regime. We have characterized the behavior of the qubits and bosonic field in both phases, and we have studied some of their critical properties near the transition. With respect to the one-photon case, fundamental differences arise in the behavior of the bosonic field $\hat{a}$, which does not acquire macroscopic occupation in the second phase, and whose fluctuations do not diverge at the critical point, at least at the order considered. The pseudospin $\overrightarrow{J}$, on the other hand, does exhibit diverging fluctuations at the transition that could lead to observable phenomena which would mark the transition in experimentally accessible situations. Indeed, as an extension of this work, it would be interesting to analyze the spin-squeezing properties of the ground state of the system for a finite number of qubits. Quantum-phase transitions and the ground-state properties could be analyzed also modifying the symmetries of the model Baksic14 , considering anisotropic couplings or including multiple bosonic fields. Finally, it would be interesting to study dynamical features and the emergence of quantum chaos Emary03 ; Alvermann12 ; Bakemeier13 in the two-photon Dicke model.

ACKNOWLEDGMENTS

This work was supported by the French Agence Nationale de la Recherche (SemiQuantRoom project ANR-14-CE26-0029-01) and from University Sorbonne Paris Cite EQDOL contract. S.F. acknowledges funding from the PRESTIGE program, under the Marie Curie Actions-COFUND of the FP7. E.S. and I.L.E. acknowledge funding from Spanish MINECO/FEDER FIS2015-69983-P, Basque Government IT986-16, UPV/EHU UFI 11/55.

Appendix A Linear term expansion of the Hamiltonian

Let us consider an extension of our model, obtained by adding a linear term in $\hat{a}$:

[TABLE]

This gives, in the mean-field approximation:

[TABLE]

We can shift the $\hat{a}$ operator: $\hat{a}=\hat{c}-\alpha$. A proper choice of $\alpha$, namely, $\alpha=\frac{g_{1}^{\beta}}{\omega+2g_{2}^{\beta}}$, allows to suppress the linear term in $\hat{c}$. This gives us:

[TABLE]

which is just Hamiltonian (5b) with an additional constant term. Hence, $\hat{a}$ in the ground state is in a displaced squeezed state, with squeezing factor $r_{c}=\frac{1}{2}\operatorname{arctanh}(\frac{2g_{2}^{\beta}}{\omega})$ and $\langle\hat{a}\rangle_{GS}=-\frac{g_{1}^{\beta}}{\omega+2g_{2}^{\beta}}$. There is also a correction on the values of $\beta$ and $E_{G}$. Importantly, the presence of the term $\frac{(g_{1}^{\beta})^{2}}{\omega+2g_{2}^{\beta}}$ in (18), which is not invariant under the transformation $\beta\rightarrow-\beta$, lift the degeneracy between $\beta>0$ and $\beta<0$.

Appendix B Description of the method used to go beyond mean-field

We are now going to detail the method used in article Hwang JC and article Hwang Rabi . We begin with (10), and we assume that we have a two-dimensional pseudospin operator $\hat{\sigma}_{z}$. $\hat{H}_{0}$ is then a block diagonal operator with respect to the $\hat{\sigma}_{z}$ eigenbasis, and $\hat{H_{c}}$ is block anti-diagonal with respect to this eigenbasis (in the following, we will simply talk about block diagonal and block anti-diagonal operators, and about pseudospin eigenstates and eigenspaces: the reference to $\hat{\sigma}_{z}$ is implicit). We apply the unitary transformation $\hat{U}=e^{\hat{S}}$ to $\hat{H}$, which yields:

[TABLE]

with $[\hat{H},\hat{S}]^{(k)}=[[\hat{H},\hat{S}]^{(k-1)},\hat{S}]$ and $[\hat{H},\hat{S}]^{(0)}=\hat{H}$. We impose $\hat{S}$ to be block anti-diagonal. If we split $\hat{H}^{\prime}$ into its block diagonal $\hat{H}^{\prime}_{d}$ and block anti-diagonal $\hat{H}^{\prime}_{a}$ parts, we get:

[TABLE]

These equations stem from the following relations: consider $D_{1}$ and $D_{2}$ two arbitrary block diagonal operators, and $O_{1}$ and $O_{2}$ two arbitrary block anti-diagonal operators. Then, $[D_{1},D_{2}]$ and $[O_{1},O_{2}]$ are block diagonal and $[O_{1},D_{1}]$ is block anti-diagonal.

At this point, the idea is to expand $\hat{S}$ with respect to $\epsilon$: $\hat{S}=\epsilon\hat{S}_{1}+\epsilon^{2}\hat{S}_{2}+...$, with $\hat{S}_{i}$ block anti diagonal for all $i$. By properly choosing the $\hat{S}_{i}$ for $i$ from $1$ to $p$ with $p$ arbitrary, it is possible to cancel $\hat{H}^{\prime}_{a}$ up to order $\epsilon^{p}$, thus decoupling the pseudospin eigenspaces up to this order.

Our case, however, is slightly different, because we have an infinite number of eigenstates instead of only two. When SU(1,1) is presented in terms of quadratic combinations of creation/annihilation operators the relevant representations are those with Bargmann parameter $q=\frac{1}{4}$ and $q=\frac{3}{4}$. The first one corresponds to even number of $\hat{b}$ excitations, while the second one is associated with an odd number article Duan Bargmann representation and Rabi analytic solution . They are not connected by any of the operators of the algebra. Since we will focus on the ground state, we concentrate on the $q=\frac{1}{4}$ case. Then, we have an infinity of eigenstates $|n\rangle_{K}$, with $\hat{K}_{0}|n\rangle_{K}=(n+\frac{1}{4})|n\rangle_{K}$. This means that we have to modify slightly the method above if we want to apply it to bigger Hilbert spaces.

For this, let us consider the operators diagonal in the $|n\rangle_{K}$ basis, which can be written as $\sum_{n=0}^{\infty}\alpha^{0}_{n}|n\rangle_{K}\langle n|_{K}$ with arbitrary $\alpha^{0}_{n}$. Let us call the ensemble of these operators $M_{(0)}$; $\hat{K}_{0}$, for instance, belongs to this ensemble. Then, let us call $M_{(1)}$ the ensemble of all operators of the form $\sum_{n=0}^{\infty}\alpha^{1}_{n}|n\rangle_{K}\langle n+1|_{K}+\sum_{n=1}^{\infty}\beta^{1}_{n}|n\rangle_{K}\langle n-1|_{K}$, with arbitrary coefficients; $\hat{K}_{+}+\hat{K}_{-}$ is an element of $M_{(1)}$. $M_{(2)}$ contains the operators of the form $\sum_{n=0}^{\infty}\alpha^{2}_{n}|n\rangle_{K}\langle n+2|_{K}+\sum_{n=2}^{\infty}\beta^{2}_{n}|n\rangle_{K}\langle n-2|_{K}+\sum_{n=0}^{\infty}\gamma^{2}_{n}|n\rangle_{K}\langle n|_{K}$; and in a general way, we define the ensemble $M_{(j)}$ that contains operators that can be written as:

[TABLE]

Let us note here that there is a redundancy in this definition: an element of $M_{(j)}$ also belongs to $M_{(j+2)}$, for all $j$. We will use the following property: for all $\hat{A}$ belonging to $M_{(i)}$ and $\hat{B}$ belonging to $M_{(j)}$, $\hat{C}=[\hat{A},\hat{B}]$ belongs to $M_{(i+j)}$, which we denote symbolically by

[TABLE]

The idea in the following will be to isolate the elements of the various $M_{(i)}$ and to cancel those that couple the eigenstates of $\hat{K}_{0}$ to a certain order.

Appendix C Study of the first phase

As indicated in the main text, we perform the Holstein-Primakoff transformation, as well as the transformation (12-14), and we consider the case $\beta=0$: $\hat{b}=\beta+\hat{d}=\hat{d}$. This gives us for the Hamiltonian:

[TABLE]

With a redefinition of parameters $\omega_{1}=\frac{\omega_{q}N}{2\omega}=O(1)$ and $\omega_{2}=\frac{g}{\omega}=O(1)$, we get:

[TABLE]

which is explicitly of the form (10); $\hat{H}_{0}=\hat{K}_{0}+\frac{\omega_{1}}{N}\hat{d}^{\dagger}\hat{d}$ does not couple the $\hat{K}_{0}$ eigenstates, but $\hat{H}_{c}=\omega_{2}(\hat{d}+\hat{d}^{\dagger})(\hat{K}_{+}+\hat{K}_{-})$ does. In our case, the small parameter is $\frac{1}{\sqrt{N}}$.

The generalized method goes as follows: all the operators implied in our calculations can be (not uniquely) decomposed into elements belonging to the different $M_{(i)}$: $\hat{A}=\sum_{i=0}^{n_{max}}\hat{A}^{(i)}$, where $n_{max}$ can take an arbitrary high value. We then propose the following decompositions:

[TABLE]

[TABLE]

with $\hat{P}$ and $\hat{Q}$ that belong to $M_{(1)}$ and $\hat{R}$ that belongs to $M_{(2)}$. Using (22) and (19), we then obtain the following decomposition for $\hat{H}_{1}^{\prime}=e^{-\hat{S}}\hat{H}_{1}e^{\hat{S}}$:

[TABLE]

since $\hat{H}_{1}^{\prime(1)}$ couples the eigenvalues of $\hat{K}_{0}^{\prime}$, we are going to cancel it at order $\frac{1}{N}$. For this purpose, we need $\hat{V}^{(1)}+[\hat{K}_{0},\hat{P}]=0$ and $[\hat{K}_{0},\hat{Q}]=0$. We note that these conditions can be met by the following choice of operators:

[TABLE]

As of the $\hat{H}_{1}^{\prime(2)}$ term, we find, using the above expressions for $\hat{P}$ and $\hat{Q}$:

[TABLE]

In general, we should choose $\hat{R}$ so as to cancel the non-diagonal terms of this operator. Here, though, $-\frac{4\omega_{2}^{2}}{2N}(\hat{d}+\hat{d}^{\dagger})^{2}\hat{K}_{0}$ is already diagonal; it is thus sufficient to set $\hat{R}=0$. We speculate that, for higher order expansions, if $\hat{H}$ can be made diagonal at order $\left(\frac{1}{\sqrt{N}}\right)^{i-1}$, we can make it diagonal at order $\left(\frac{1}{\sqrt{N}}\right)^{i}$ by adding terms of the form $\left(\frac{1}{\sqrt{N}}\right)^{i}\hat{T}_{p}$ to the expansion of $\hat{S}$, with $\hat{T}_{p}$ that belongs to $M_{(p)}$ and $p\leq i$. It could be interesting in a future study to consider higher-order terms in a more systematic way to confirm or to invalidate this conjecture.

In the end, all those manipulations yield:

[TABLE]

This Hamiltonian does commute with $\hat{K}_{0}$; projection in the $\hat{K}_{0}$ ground state gives $\hat{K}_{0}\rightarrow\frac{1}{4}$ (according to the definition (12)). Restoring the constants in the Hamiltonian yields the operator described in the main text:

[TABLE]

This effective Hamiltonian describes the $\hat{b}$ fluctuations near the ground state. Being quadratic in $\hat{d}$, it can be diagonalized by using a Bogoliubov transformation of parameter $r_{s}^{(1)}=\frac{1}{4}\ln(1-\frac{4g^{2}}{N\omega\omega_{q}})<0$, which corresponds to a squeezed vacuum state as described in the main text. The state of $\hat{a}$ is given by the ground state of $\hat{K}_{0}$; here it is just the coherent vacuum state, as in the mean-field scenario. Finally, the ground state energy and the $\hat{d}$ excitation energy are computed: we obtain

[TABLE]

Appendix D Study of the second phase

We proceed in a similar fashion, using HP transformation once more, and setting $\hat{b}=\beta+\hat{d}$; but this time $\beta\neq 0$. As a starting point, we use the value of $\beta$ obtained by our mean-field analysis. For ease of notation, we define the following parameters:

[TABLE]

Expanding the Hamiltonian gives:

[TABLE]

Where we have defined new operators and parameters:

[TABLE]

[TABLE]

Let us note that the definition (35) amounts to a Bogoliubov transformation of parameter $r_{a}^{(2)}$ for the $\hat{a}$ field. We seek to decouple the eigenspaces of $\hat{K}_{0}^{\prime}$; for this, we apply a transformation $e^{-\hat{S}}\hat{H}_{2}e^{\hat{S}}$ with $\hat{S}=\frac{1}{\sqrt{N}}\hat{S}_{1}+\frac{1}{N}\hat{S}_{2}$. Using the procedure described earlier, we propose the following operators:

[TABLE]

[TABLE]

[TABLE]

This gives us an Hamiltonian that commutes with $\hat{K}_{0}^{\prime}$:

[TABLE]

As previously, projection in the ground state of $\hat{K}_{0}^{\prime}$ yields an effective Hamiltonian describing the fluctuations of $\hat{d}$ above the ground state. Note that in contrast to the first phase, $e^{-\hat{S}}\hat{H}_{2}e^{\hat{S}}$ contains a term linear in $\hat{d}+\hat{d}^{\dagger}$. $e^{-\hat{S}}\hat{H}_{2}e^{\hat{S}}$ can be seen as an effective potential for the $\hat{d}$ (or $\hat{b}$) field; adding this linear term shifts its minimum, which amounts to change the value of $\beta$. In our case, the term adds a correction of order $\frac{1}{\sqrt{N}}$ to the value of $\beta$ (keeping in mind that the mean-field value of $\beta$ is of order $\sqrt{N}$). We can thus absorb the linear term by adding this correction to $\beta$; once it is done, $e^{-\hat{S}}\hat{H}_{2}e^{\hat{S}}$ becomes quadratic in $\hat{d}$ and can be diagonalized by a Bogoliubov transformation of parameter $r_{s}^{(2)}$, which gives squeezing properties for $\hat{b}$. Going back to the $\overrightarrow{J}$ operators, we have:

[TABLE]

with $\hat{X}_{d}$ et $\hat{P}_{d}$ the quadratures of the $\hat{d}$ field. Using the squeezing properties of $\hat{d}$ in both phases, we retrieve the properties described in the main text and in Fig.3. As of the state of $\hat{a}$, it corresponds to the ground state of $\hat{K}_{0}^{\prime}$; according to the definition (35), this is just a squeezed vacuum state with squeezing parameter $r_{a}^{(2)}$. This result is identical to what we obtained in the mean-field scenario (6), up to a small correction in the definition of $r_{a}^{(2)}$: $r_{a}^{(2)}-r_{a}^{(MF)}=\frac{1}{2}\operatorname{arctanh}\left(\frac{2g_{\beta}^{(MF)}}{\omega}+O\left(\frac{1}{N}\right)\right)-\frac{1}{2}\operatorname{arctanh}\left(\frac{2g_{\beta}^{(MF)}}{\omega}\right)=O(\frac{1}{N})$ (while $r_{a}^{(MF)}$ is of order O(1)).

Finally, the computation of the $\hat{d}$ excitation energy gives:

[TABLE]

For the sake of simplicity, we do not show here the expression of the ground state energy $E_{G}^{(2)}$ for the second phase. The plot on Fig.4 is obtained by combining $E_{exc}^{(1)}$ and $E_{exc}^{(2)}$.