Phases of cold atoms interacting via photon-mediated long-range forces

TL;DR

This paper explores the phase behavior of cold atoms interacting via photon-mediated long-range forces in a two-mode cavity, revealing phases analogous to the Generalized Hamiltonian Mean Field model and discussing ensemble inequivalence.

Contribution

It demonstrates that a two-mode cavity system exhibits phases of the Generalized Hamiltonian Mean Field model, linking atomic spatial order to cavity light emission.

Findings

Identification of paramagnetic, nematic, and ferromagnetic phases in atomic spatial arrangements.

Correlation between cavity emitted light and atomic phase states.

Potential to study ensemble inequivalence in this optical system.

Abstract

Atoms in high-finesse optical resonators interact via the photons they multiply scatter into the cavity modes. The dynamics is characterized by dispersive and dissipative optomechanical long-range forces, which are mediated by the cavity photons, and exhibits a steady state for certain parameter regimes. In standing-wave cavities the atoms can form stable spatial gratings. Moreover, their asymptotic distribution is a Maxwell-Boltzmann whose effective temperature is controlled by the laser parameters. In this work we show that in a two-mode standing-wave cavity the stationary state possesses the same properties and phases of the Generalized Hamiltonian Mean Field model in the canonical ensemble. This model has three equilibrium phases: a paramagnetic, a nematic, and a ferromagnetic one, which here correspond to different spatial orders of the atomic gas and can be detected by means of the light emitted by the cavities. We further discuss perspectives for investigating in this setup the ensemble inequivalence predicted for the Generalized Hamiltonian Mean Field model.

Equations

$$H=\sum _j\frac{p_j^2}{2}+\frac{1}{2N}\sum _{i,j}\left(1-\Delta \cos {\theta }_{ij}-(1-\Delta )\cos 2{\theta }_{ij}\right),$$

$$\frac{\partial }{\partial t}\hat{\rho }=\frac{1}{i\hslash }\left[{\hat{H}}_0,\hat{\rho }\right]-\sum _{n=1,2}{\kappa }_n\left({\hat{a}}_n^{†}{\hat{a}}_n\hat{\rho }+\hat{\rho }{\hat{a}}_n^{†}{\hat{a}}_n-2{\hat{a}}_n\hat{\rho }{\hat{a}}_n^{†}\right),$$

$${\hat{H}}_0=\sum _{i=1}^N\frac{{\hat{p}}_i^2}{2m}$$

$$f_{\mathrm{tot}}(\mathbf{x},\mathbf{p},t)=\int \frac{d^{N}y}{(2\pi \hslash {)}^N}{\mathrm{e}}^{-\frac{i}{\hslash }\mathbf{y}\cdot \mathbf{p}}\mathrm{Tr}\left\{|\mathbf{x}-\frac{1}{2}\mathbf{y}⟩⟨\mathbf{x}+\frac{1}{2}\mathbf{y}|\hat{\rho }(t)\right\},$$

$$\epsilon =\frac{k\Delta p/m}{|{\kappa }_n+i{\Delta }_n|}$$

$$\frac{\partial }{\partial t}f(\mathbf{x},\mathbf{p},t)+\left\{f(\mathbf{x},\mathbf{p},t),H(\mathbf{x},\mathbf{p})\right\}=$$

$$\sum _{i,j=1}^N\sum _{n=1,2}\frac{\partial }{\partial p_i}\left[\sin (k_n{x}_i)\sin (k_n{x}_j)\left(D_n\frac{\partial f(\mathbf{x},\mathbf{p},t)}{\partial p_j}-{\Gamma }_n{p}_{j}f(\mathbf{x},\mathbf{p},t)\right)\right]$$

$$+\sum _{i,j=1}^N\sum _{n=1,2}\frac{\partial }{\partial p_j}\left[{\eta }_n\sin (k_n{x}_i)\sin (k_n{x}_j)\frac{\partial }{\partial x_i}f(\mathbf{x},\mathbf{p},t)\right].$$

$$H(\mathbf{x},\mathbf{p})=\sum _{i=1}^N\frac{p_i^2}{2m}-N\sum _{n=1,2}{\gamma }_n{\Theta }_n^2,$$

$${\Theta }_1=\frac{1}{N}\sum _{i=1}^N\cos (k{x}_i/2);{\Theta }_2=\frac{1}{N}\sum _{i=1}^N\cos (k{x}_i).$$

$${\mathcal{E}}_n=\frac{N{S}_n{\Theta }_n}{{\Delta }_n+i{\kappa }_n}.$$

$${\alpha }_n=\frac{4N{S}_n^2{\Delta }_n^2}{({\Delta }_n^2+{\kappa }_n^2{)}^2},$$

$${\beta }_n=\frac{-4{\Delta }_n}{\hslash ({\Delta }_n^2+{\kappa }_n^2)}.$$

$$D_n$$

$${\Gamma }_n$$

$${\eta }_n$$

$$\frac{{\Delta }_1}{{\Delta }_1^2+{\kappa }_1^2}=\frac{{\Delta }_2}{{\Delta }_2^2+{\kappa }_2^2}.$$

$${\beta }^{-1}=\frac{\hslash ({\Delta }_n^2+{\kappa }_n^2)}{-4{\Delta }_n}.$$

$$f_{\mathrm{st}}=C\exp (-\beta H).$$

$$Z=\frac{1}{h^N}{\int }_0^{\lambda }d{x}_1\dots {\int }_0^{\lambda }d{x}_N{\int }_{-\infty }^{\infty }d{p}_1\dots {\int }_{-\infty }^{\infty }d{p}_N{e}^{-\beta H(\text{x},\text{p})},$$

$$Z=\frac{N\sqrt{{\alpha }_1{\alpha }_2}}{\pi {\left(\pi \hslash {\omega }_r\beta \right)}^{\frac{N}{2}}}{\int }_{-\infty }^{\infty }d{y}_1{\int }_{-\infty }^{\infty }d{y}_2\exp \left[-N\left\{{\alpha }_1{y}_1^2+{\alpha }_2{y}_2^2-\ln \left(\mathcal{I}(y_1,y_2)\right)\right\}\right],$$

$$\mathcal{I}(y_1,y_2)={\int }_0^{2\pi }dq\exp \left[2\left({\alpha }_1{y}_1\cos (q)+{\alpha }_2{y}_2\cos (2q)\right)\right].$$

$$\mathcal{F}(y_1,y_2)=\frac{\ln \left(\pi \hslash {\omega }_{\mathrm{r}}\beta \right)}{2\beta }+\frac{1}{\beta }\underset{y_1,y_2}{\inf }\left\{{\alpha }_1{y}_1^2+{\alpha }_2{y}_2^2-\ln \left(\mathcal{I}(y_1,y_2)\right)\right\}.$$

$$y_n^{\ast }=\frac{{\int }_0^{2\pi }dq\cos (nq)\exp \left[2({\alpha }_1{y}_1^{\ast }\cos (q)+{\alpha }_2{y}_2^{\ast }\cos (2q))\right]}{{\int }_0^{2\pi }dq\exp \left[2({\alpha }_1{y}_1^{\ast }\cos (q)+{\alpha }_2{y}_2^{\ast }\cos (2q))\right]},$$

$$ϵ=\frac{1}{2\beta }-\frac{1}{\beta }\left({\alpha }_1{\Theta }_1^2+{\alpha }_2{\Theta }_2^2\right),$$

$${\alpha }_1<1,{\alpha }_2<1.$$

$${\Theta }_2=\frac{I_1\left(2{\alpha }_2{\Theta }_2\right)}{I_0\left(2{\alpha }_2{\Theta }_2\right)},$$

$${\alpha }_1<\frac{1}{1+{\Theta }_2},$$

$$1<$$

Full text

1 Introduction

Atomic ensembles in optical resonators offer a promising platform for studying the physics of long-range interacting systems [1]. The long-range interaction here originates from multiple scattering of cavity photons, which carry the information about the positions of the scattering atoms and thus mediate an optomechanical interparticle potential [2]. In a single-mode cavity the photons are coherent over the cavity mode volume, which makes the interaction range as large as the system size. Therefore, the energy is non-additive like in gravitational and Coulomb systems in two or more dimensions [1, 3, 4].

In equilibrium statistical mechanics, consequences of non-additivity are for instance the super-linear scaling of thermodynamic quantities with the system size and the inequivalence of the statistical ensembles [4], one manifestation of which are negative specific heats in the microcanonical ensemble [4, 5, 6]. Differing from these systems, however, the dynamics of atomic gases in optical cavities is typically dissipative and non-trivial effects can only be observed if either the atoms or the cavity are pumped by light [2, 7]. The steady state, when it exists, results from the dynamical interplay between drive and losses and its properties thus depend on the drive and on the cavity parameters. It is therefore often not possible to draw a direct connection with equilibrium statistical mechanics of long-range interacting systems.

In this context it is remarkable that for some parameter regimes the dynamics of atoms’ spatial selforganization in an optical resonator can be mapped to long-range interacting systems at equilibrium [2, 8, 9]. Selforganization of the atomic gas in ordered spatial patterns occurs in a single-mode standing-wave resonator when the atoms are driven by lasers whose intensity exceeds a threshold value, which depends also on the cavity decay rate [2, 8, 9, 10, 11]. By suitably tuning the laser frequency, moreover, a stationary state exists which is characterised by a Maxwell-Boltzmann distribution of the atomic momentum. In Ref. [11] it was shown that the stationary dynamics can be mapped to the one of the Hamiltonian-Mean-Field model in a canonical ensemble [4, 12] and in particular that the transition to spatial order can be described in terms of a Landau second-order phase transition. The dynamics leading to equilibrium, moreover, exhibits a slow relaxation that is due to the interplay between the conservative and dissipative cavity-mediated long-range forces [7].

In this paper we extend the model of Ref. [11] and consider a gas of cold atoms that interact with two cavity modes and are transversally driven by lasers. A possible setup is illustrated in Fig. 1. We determine the parameter regimes where the dynamics asymptotically tends to a stationary state and show that its phase diagram as a function of the lasers’ and of the cavity parameters can be mapped to the one of the Generalized Hamiltonian Mean Field model (GHMF) in a canonical ensemble [13, 14, 15]. This model describes the dynamics of $N$ particles with canonically conjugated variables $p_{j}$, $\theta_{j}$ constrained on a circle that interact via competing long-range forces. In the form studied in Ref. [13, 15] its Hamiltonian reads

[TABLE]

where $\theta_{ij}=\theta_{i}-\theta_{j}\in[0,2\pi)$ and $\Delta$ is a dimensionless parameter that can vary continuously in the interval $[0,1]$. The phase diagram as a function of the temperature and $\Delta$ is characterised by (i) a paramagnetic, (ii) a nematic, and (iii) a ferromagnetic phase, with first and second-order transitions. In our case the phases correspond to density modulations of the atoms at different periodicity and can be detected through the light emitted by the cavity. Our motivation draws from ongoing experimental investigations [16]. Theoretical studies of this system focused on the dynamics leading to equilibrium and reported the existence of several metastable states [17]. These properties are at the basis of proposals for using these systems to simulate a quantum Hopfield associative memory scheme [18, 19]. The determination of the condition for a stationary state of the setup in Fig. 1 and the analysis of the corresponding phase diagram is the main result of the present manuscript. The mapping to the GHMF model shows that Cavity Quantum Electrodynamics (CQED) setups offer a versatile platform for studying the statistical mechanics of systems with long-range interactions.

This article is organised as follows. In Sec. 2 we introduce the physical model and sketch the derivation of a Fokker-Planck equation governing the dynamics of the atoms’ external variables in the semiclassical limit. We then determine the parameters’ regime for which the Fokker-Planck equation allows for a stationary state which is a Maxwell-Boltzmann distribution. In Sec. 3 we define the free energy, which we can associate to the stationary state, and introduce an appropriate thermodynamic limit. We then show that the free energy can be mapped to the one of the GHMF model. We determine the phase diagram as a function of the system’s parameters and identify the observables which allow one to measure the predicted phases. In Sec. 4 we discuss possible implementations of the setup which could correspond to the realization of the microcanonical ensemble of the GHMF model, for which ensemble inequivalence has been predicted [15], and identify the parameters regimes for which it could be measured.

2 Semiclassical dynamics of an atomic gas in an optical cavity

The system we consider consists of a gas of $N$ cold atoms of mass $m$, whose motion is confined to occur along one dimension parallel to the unit vector ${\bf e}_{x}$. We denote by $\hat{x}_{j}$ and $\hat{p}_{j}$ the canonically conjugated position and momentum operators ($j=1,\ldots,N$), such that $\left[\hat{x}_{i},\hat{p}_{j}\right]=i\hbar\delta_{ij}$. The atoms experience the optomechanical forces due to the interaction with the lasers and with the quantized fields of two high-finesse cavities in the setup of Fig. 1. Specifically, the axes of the two cavities are in the $x-y$ plane, the wave vector ${\bf k}_{1}$ of cavity 1 forms an angle $\varphi$ with the $x$-axis and thus the force the atoms experience is the projection of the mechanical force along $x$, while the wave vector ${\bf k}_{2}$ of cavity 2 is parallel to ${\bf e}_{x}$. The laser fields are linearly polarised and propagate in the direction orthogonal to the plane, they pump the cavity fields by means of coherent scattering via the atoms. The amplitude of coherent scattering, in turn, is maximal when the atoms form Bragg gratings, whose stability depends on the mechanical forces of the cavity light. As we will show, a nematic phase corresponds to a stable Bragg grating which supports the build-up of the field of only one cavity mode. In the ferromagnetic phase, instead, the atoms form stable Bragg gratings for both modes.

Below we describe the setup in detail and introduce the master equation for the density matrix $\hat{\rho}$ of atoms and cavity fields which governs the system’s dynamics. We then sketch the derivation of a Fokker-Planck equation for the motion of the atoms’ external degrees of freedom, which is valid when the atomic variables can be treated as semiclassical variables and the cavity fields can be eliminated from the equations of motion in a coarse-grained time scale. We finally determine the stationary state of the atoms and identify the regime in which it is a thermal state.

2.1 Master equation

The state of the atoms’ external degrees of freedom and of the cavity modes is described by the density operator $\hat{\rho}$, whose dynamics is governed by a Born-Markov master equation of the form:

[TABLE]

where $\hat{H}_{0}$ is the Hamiltonian of the system, which we introduce below, and the second term of the right-hand side describes photon emission by the cavity modes at rate $\kappa_{n}$ ($n=1,2$). Here, $\hat{a}^{\dagger}_{n}$ and $\hat{a}_{n}$ denote the creation and annihilation operators of a photon of the standing-wave cavity mode $n$ ($n=1,2$), with wave vector ${\bf k}_{n}$, frequency $\omega_{c,n}=c|{{\bf k}_{n}}|$, and linear polarization in the $x-y$ plane. The commutation relations are $\left[\hat{a}_{i},\hat{a}^{\dagger}_{j}\right]=\delta_{ij}$.

Hamiltonian $\hat{H}_{0}$ describes the optomechanical coupling between atoms’ degrees of freedom and cavity modes. It is an effective Hamiltonian derived in the limit where the atoms’ internal degrees of freedom can be adiabatically eliminated, such that the scattering processes are coherent and the relevant parameters of the atomic internal structure is the atoms’ polarizability [20]. Here, cavity mode $n$ couples with the electronic transition $|g\rangle\to|e_{n}\rangle$ at frequency $\omega_{a,n}$ with vacuum Rabi frequency $g_{n}$. We further assume that the coupling of mode 1 (2) with $|g\rangle\to|e_{2}\rangle$ ($|g\rangle\to|e_{1}\rangle$) is off-resonance by orders of magnitude and can be discarded (nonetheless, the wave numbers are assumed to be $|{\bf k}_{1}|\approx|{\bf k}_{2}|=k$). In this limit cavity $n$ is pumped by coherent scattering of the laser, which couples to the transition $|g\rangle\to|e_{n}\rangle$ with Rabi frequency $\Omega_{p,n}$ and frequency $\omega_{p,n}$. The condition for adiabatic elimination of the internal excited state is given by the inequality $|\omega_{a,n}-\omega_{p,n}|,|\omega_{a,n}-\omega_{c,n}|\gg\Omega_{p,n},g_{n}\sqrt{\bar{n}_{c,n}},|\Delta_{n}|$, where $\bar{n}_{c,n}$ is the mean intracavity photon number in cavity $n$ and $\Delta_{n}=\omega_{p,n}-\omega_{c,n}$ is the detuning of the laser from the cavity mode it pumps [3]. In this regime $\hat{H}_{0}$ reads:

[TABLE]

and is here reported in the frame where each atomic transition and cavity mode rotates at the corresponding laser frequency. Beside the kinetic energy of the atoms, it contains the resonators’ energy, which is shifted by the dynamical Stark shift with amplitude $U_{n}$ induced by the coupling between cavity mode and the atoms at the position $x_{i}$ within the cavity spatial mode function $\cos(k_{n}\hat{x}_{i})$. This term is also a periodic potential for the atoms whose depth is a dynamical variable. The last term on the right-hand side, finally, describes coherent scattering by the atoms between laser and cavity mode with coupling strength $S_{n}=g_{n}\Omega_{p,n}/(\omega_{p,n}-\omega_{a,n})$. It is an effective pump of the resonator whose amplitude is maximal when the atoms form Bragg gratings maximizing the expectation value of the operator $\sum_{i=1}^{N}\cos(k_{n}\hat{x}_{i})$.

Note that in Eq. (3) we introduced the notation $k_{1}\equiv|{\bf k}_{1}\cdot{\bf e}_{x}|=k\cos\varphi$ and $k_{2}\equiv|{\bf k}_{2}\cdot{\bf e}_{x}|=k$. In the following we will set $\varphi=\pi/3$, thus $k_{1}=k/2$.

2.2 Fokker-Planck equation for the atoms’ external variables

We now discuss the assumptions at the basis of the derivation of a Fokker-Planck equation (FPE) for the dynamics of the atomic external variables. A semiclassical description of the atoms’ center-of-mass motion is justified when the width $\Delta p$ of the single-atom momentum distribution is much larger than the linear momentum $\hbar k$ carried by a cavity photon $\Delta p\gg\hbar k$ [21]. In this limit it is convenient to use the Wigner function $f_{\rm tot}({\bf x},{\bf p},t)$ for the atomic variables ${\bf x}=(x_{1},\ldots,x_{N})$, ${\bf p}=(p_{1},\ldots,p_{N})$:

[TABLE]

with ${\bf y}=(y_{1},\ldots,y_{N})$. We further assume that the cavity field relaxes very fast to a local steady state depending on the atomic distribution, which is verified when the inequality $k\Delta p/m\ll|\kappa_{n}+i\Delta_{n}|$ is fulfilled, namely, when the dimensionless parameter

[TABLE]

is small. This allows us to identify a coarse-grained time scale $\Delta t$ that is infinitesimal for the external degrees of freedom but over which the cavity degrees of freedom can be eliminated from the equations of the atomic dynamics. In particular, $f_{\rm tot}({\bf x},{\bf p},t)=f({\bf x},{\bf p},t)+f_{\rm na}({\bf x},{\bf p},t)$, where $f({\bf x},{\bf p},t)$ is the term in zero order in the retardation effect, corresponding to the cavity field following adiabatically the atomic motion, and $f_{\rm na}({\bf x},{\bf p},t)$ represents the non-adiabatic corrections scaling with $\varepsilon$. The latter can be expressed in terms of $f({\bf x},{\bf p},t)$ using perturbation theory [22, 23]. The derivation is lengthy but is also a straightforward extension of the derivation for a single-mode cavity, which is extensively reported in Ref. [23]. We thus refer the interested reader to this work and present here the resulting FPE, which reads

[TABLE]

In detail, Hamiltonian (7) results from the adiabatic component of the dynamics and describes coherent long-range, two-body interactions which are mediated by the cavity photons:

[TABLE]

where $\gamma_{n}$ is a scalar and

[TABLE]

The quantities $\Theta_{n}$ are order parameters for spatial selforganization. In fact, they vanish for homogeneous spatial distributions, they are both different from zero when the atomic density forms spatial gratings with periodicity $4\pi/k$, while $\Theta_{1}=0$ and $\Theta_{2}\neq 0$ when the spatial grating has periodicity $2\pi/k$, as illustrated in Fig. 2. The cavity mode field amplitudes $\mathcal{E}_{n}=\langle\hat{a}_{n}\rangle$, in turn, are proportional to $\Theta_{n}$ and in leading order in the expansion in $\varepsilon$ read [11]

[TABLE]

The Bragg gratings can be stable provided that $\gamma_{n}>0$. The sign of $\gamma_{n}$ is here controlled by the detuning $\Delta_{n}$, and hence by the frequency of the pumping laser. For later convenience we write $\gamma_{n}=\alpha_{n}/\beta_{n}$ with

[TABLE]

Friction and diffusion are instead due to retardation effects between atoms and cavity dynamics and describe cross-correlations between the atoms, which can play a relevant role in stabilizing the system in non-thermal metastable states [7]. Their coefficients take the form

[TABLE]

and are here reported in the limit $|\Delta_{n}|\gg NU_{n}$, where we neglected the contribution of the dynamical Stark shift to the dynamics.

2.3 Existence of a stationary state

The FPE (6) allows for a stationary solution satisfying the condition $\partial_{t}f_{\mathrm{st}}({\bf x},{\bf p},t)=0$. We first consider two limiting cases, in which only one cavity mode is pumped. These situations correspond to the dynamics of atoms in a single-mode standing-wave cavity investigated in Refs. [3, 11, 23, 24].

Let us first assume that $S_{1}=0$ but $S_{2}\neq 0$. In this case a stationary state exists provided that $\Delta_{2}<0$ and the stationary distribution reads $f_{\mathrm{st}}=C_{2}\exp(-\beta_{2}H|_{\gamma_{1}=0})$, with normalization constant $C_{2}$ and $\beta_{2}=-\frac{m\Gamma_{2}}{D_{2}}$ given in Eq. (11) 111In this discussion we neglect the terms of Eq. (6) that scale with $\eta_{n}$. This is exact if $\Delta_{n}=-\kappa_{n}$. In general, these terms give rise to corrections to the coherent dynamics that scale with $1/N$ in the thermodynamic limit we use in Sec. 3, see Ref. [24] for details.. Vice versa, when $S_{1}\neq 0$, $S_{2}=0$, and $\Delta_{1}<0$, the stationary state is $f_{\mathrm{st}}=C_{1}\exp(-\beta_{1}H|_{\gamma_{2}=0})$ with $\beta_{1}=-\frac{m\Gamma_{1}}{D_{1}}$ from Eq. (11).

When both resonators are pumped, a stationary solution can be found provided that $\Delta_{1},\Delta_{2}$ are negative and $\beta_{1}=\beta_{2}$, namely:

[TABLE]

This is the situation we consider in the following. In particular, $\beta_{1}=\beta_{2}=\beta$, where

[TABLE]

The detunings and cavity loss rates thus determine an effective temperature $T_{\rm eff}=k_{B}/\beta$ characterizing the stationary state. The stationary state is given by

[TABLE]

3 Mapping to the generalized Hamiltonian Mean Field model

We now consider the stationary state of Eq. (17) and draw a formal analogy to a canonical ensemble at equilibrium. For this purpose we define the thermodynamic limit, according to which $\alpha_{n}\propto NS_{n}^{2}$ is constant. This assumption warrants that the energy is extensive, it is thus equivalent to Kac scaling [4] and physically corresponds to scaling the cavity mode volume linearly with the number of particles [10, 25]. In this thermodynamic limit we obtain an explicit expression of the free energy per particle which allows us to perform a mapping of the steady state in Eq. (17) to the canonical ensemble realization of the GHMF.

We determine the free energy per particle $\mathcal{F}$ using the relation $\mathcal{F}=-\ln(Z)/(N\beta)$, where $Z$ is the “canonical” partition function and reads:

[TABLE]

with $\lambda=4\pi/k$. We integrate over the momenta p and apply the Hubbard-Stratonovich transformation to eliminate $\Theta_{n}^{2}(\textbf{x})$, obtaining

[TABLE]

with $q_{i}=kx_{i}/2$, $\omega_{r}=\hbar k^{2}/(2m)$, and

[TABLE]

In the thermodynamic limit $N\rightarrow\infty$ we perform a saddle-point approximation, which leads to the expression for the free energy

[TABLE]

This expression coincides, apart for irrelevant constants, with the mean free energy of the GHMF model in the canonical ensemble [15]. In particular, the extrema $y_{n}^{*}$ of the free energy fulfill the relation

[TABLE]

their values lie in the interval $y_{n}^{*}\in[-1,1]$ and they can be identified with the variables $\Theta_{n}$: $y_{n}^{*}=\langle\cos nq\rangle$. This shows that $\Theta_{1},\Theta_{2}$ are analogous to the magnetization in the GHMF model. By means of this mapping, moreover, we connect the lasers and cavity parameters with the dimensionless parameter $\Delta$ of the GHMF in Eq. (1): $\alpha_{1}/\beta\to\Delta$ and $\alpha_{2}/\beta\to(1-\Delta)$. Therefore, varying the lasers’ amplitudes would allow one to span over the values of $\Delta$ in Eq. (1), while the effective temperature can be tuned varying the detunings $\Delta_{1}$ and $\Delta_{2}$, and thus the lasers’ frequencies (provided condition (15) is fulfilled).

The phase diagram is obtained following the same analysis as in Ref. [15] and it is illustrated in Fig. 3. The system exhibits second and first order phase transitions (see A for details), which are found for the same corresponding values of the phase diagram as in Ref. [15]. The phases can be measured by monitoring the amplitude and the phase of the fields at the cavity output, since these are proportional to the order parameters $\Theta_{1}$ and $\Theta_{2}$, as visible in Eq. (9).

4 Discussion and outlook

The mapping of the stationary dynamics of Eq. (6) to the canonical GHMF shows that cavity QED can be a versatile platform for studying equilibrium statistical mechanics of long-range interacting systems. In this perspective it is important to identify the parameter regimes for which this setup could simulate the microcanonical GHMF model. This would allow one to experimentally investigate the ensemble inequivalence that has been predicted for the GHMF [15].

Within the validity of the semiclassical description, the microcanonical GHMF could be realised in the regime where the parameter $\varepsilon$, Eq. (5), becomes smaller than the small parameter $\hbar k/\Delta p$ of the semiclassical expansion. This would require one to choose the detunings $|\Delta_{n}|\gg\kappa_{n}$ [25]. In this limit there is a well defined time scale over which the dynamics is coherent and solely dominated by Hamiltonian (7), while the right-hand side of the FPE (6), which scales with $\varepsilon$, can be discarded. Moreover, in order to prepare the system in a microcanonical ensemble, the atomic gas shall be in the asymptotic Boltzmann-Gibbs distribution of the corresponding Hamiltonian dynamics [4, 14]. By identifying the mean energy per particle in the canonical ensemble $\langle E/N\rangle=\partial(\beta\mathcal{F})/\partial\beta$ with the constant energy $\epsilon$ in the microcanonical ensemble, we obtain the relation

[TABLE]

which allows us to directly compare the results obtained in both ensembles. Ensemble inequivalence is predicted in the regions where the three phases meet and is indicated by the shaded area in Fig. 3. In an experimental realization it would become evident by detecting different phases in the canonical and microcanonical realizations for the same values of $\alpha_{1}$ and $\alpha_{2}$ [15]. It would be interesting to identify observables of this system which provide a measure of the specific heat, thus allowing one to determine whether it may become negative in the microcanonical ensemble. This is an open question which we will address in future work.

The authors acknowledge discussions with Shamik Gupta, Julian Léonard, Stefan Schütz, Valentin Torggler, and Tobias Donner. This work was supported by the German Research Foundation (DFG DACH ”Quantum crystals of matter and light”) and by the European Commission (ITN network ”ColOpt”).

Appendix A Determination of the canonical phase diagram

For the calculation of the phase diagram in Fig. 4 we numerically calculated the global minimum of the free energy, Eq. (20), for each pair $\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2})$. The global minimum describes a paramagnetic phase when $\Theta_{1}(\boldsymbol{\alpha})=\Theta_{2}(\boldsymbol{\alpha})=0$. The phase is nematic when $\Theta_{1}(\boldsymbol{\alpha})=0$ and $\Theta_{2}(\boldsymbol{\alpha})\neq 0$ and ferromagnetic when $\Theta_{1},\Theta_{2}\neq 0$. A phase transition occurs where the properties of the minimum change by varying $\boldsymbol{\alpha}$. The order of the transition is determined by calculating numerically the first derivatives of $\mathcal{F}(\boldsymbol{\alpha})$ with respect to $\alpha_{1}$ and $\alpha_{2}$: If they are discontinuous at the phase transition the transition is of first order, while if they are continuous the transition is of second order. This determines the phase diagram in Fig. 4.

We also checked our results by analytically calculating the Hessian matrix of the free energy, Eq. (20), at the extrema where Eq. (21) is fulfilled. We note that $\Theta_{1}=\Theta_{2}=0$ is always a solution of Eq. (21). It is a minimum when the following inequalities hold:

[TABLE]

The thin black dashed line in Fig. 4 delimitates the region $\alpha_{1},\alpha_{2}<1$, where the paramagnetic phase is a local minimum of the free energy. We now consider the nematic phase, set $\Theta_{1}=0$ in Eq. (21) and obtain the equation for $\Theta_{2}$:

[TABLE]

where $I_{n}(x)$ is the modified Bessel function of order $n$. The Hessian matrix is positive definite when the following inequalities are fulfilled:

[TABLE]

One can show that imposing $\Theta_{2}\neq 0$ in Eq. (24) is equivalent to Eq. (26). Inequality (25) determines an upper threshold $\alpha_{1,c}=1/(1+\Theta_{2})$ on $\alpha_{1}$ above which the nematic configuration is no longer a minimum of the free energy. The thin black solid line in Fig. 4 shows $\alpha_{1,c}$ in the $\boldsymbol{\alpha}$-plane. Note that $\alpha_{1,c}<1$ ($\alpha_{1,c}>1$) for $\Theta_{2}>0$ ($\Theta_{2}<0$). In the limit $\alpha_{2}\to\infty$, Eq. (24) delivers $\Theta_{2}\to 1$ and $\Theta_{2}\to-1$, giving $\alpha_{1,c}\to 1/2$ and $\alpha_{1,c}\to\infty$, respectively.

We note that the conditions we determine analytically do not overlap for all values of $\alpha_{1},\alpha_{2}$ with the numerically calculated phase transition lines (see Fig. 4). For instance there is an area where the inequalities (23) hold but the global minimum is in the ferromagnetic phase. Moreover the condition in Eq. (25) predicts a parameter region where the nematic phase with $\Theta_{2}<0$ is a local minimum but the global minimum is a ferromagnetic phase.