Table Of Content

Dicke simulators with emergent collective quantum computational abilities Pietro Rotondo1, Marco Cosentino Lagomarsino2, and Giovanni Viola3,4 1 Dipartimento di Fisica, Universita` degli Studi di Milano and INFN, via Celoria 16, 20133 Milano, Italy 2Sorbonne Universités and CNRS, UPMC Univ Paris 06, UMR 7238, Computational and Quantitative Biology, 15 Rue de l’Ećole de Médecine, Paris, France 3 Department of Microtechnology and Nanoscience (MC2), Chalmers University of Technology, SE-412 96 Gothenburg, Sweden 4 Institute for Quantum Information, RWTH Aachen University, D-52056 Aachen, Germany (Dated: January 29, 2015) Using an approach inspired from Spin Glasses, we show that the multimode disordered Dicke model is equivalent to a quantum Hopfield network. We propose variational ground states for the 5 1 system at zero temperature, which we conjecture to be exact in the thermodynamic limit. These 0 groundstatescontaintheinformationonthedisorderedqubit-photoncouplings. Theseresultslead 2 totwointriguingphysicalimplications. First,oncethequbit-photoncouplingscanbeengineered,it shouldbepossibletobuildscalablepattern-storingsystemswhosedynamicsisgovernedbyquantum n laws. Second, we argue with an example how such Dicke quantum simulators might be used as a a solver of “hard” combinatorial optimization problems. J 8 PACSnumbers: 2 ] The connection of experimentally realizable quantum n n systems with computation contains promising perspec- - tives from both the fundamental and the technological s viewpoint [1, 2]. For example, quantum computational i d capabilities can be implemented by “quantum gates” [3] . t and by the so-called “adiabatic quantum optimization” a m technique [4–6]. Today’s experimental technology of highlycontrollablequantumsimulators,recentlyusedfor - d testingtheoreticalpredictionsinawiderangeofareasof n physics[7–9],offersnewopportunitiesforexploringcom- o puting power for quantum systems. c [ In the case of light-matter interaction at the quantum 2 level, the reference benchmark is the Dicke model [10]. FIG. 1: In the Dicke model, photons (yellow lines) mediate v Studiesofitsequilibriumpropertieshavepredictedasu- a long range interaction between qubits (green circles). The 9 perradianttransitiontooccurinthestrongcouplingand drawingsketchesschematicallyasixqubitssystemwithinits 0 lowtemperatureregime[11–13]. Thesuperradiantphase fully-connected graph and its internal level structure. In the 2 standard single-mode Dicke model the exchange coupling is ischaracterizedbyamacroscopicnumberofatomsinthe 2 fixed at the same value for every pair of qubits. In systems excited state whose collective behaviour produces an en- 0 where both many modes and disorder are present, the ex- 1. hancement of spontaneous emission (proportional to the changecouplingsarequbit-dependentandtaketheformgiven 0 number of cooperating atoms in the sample). Crucially, by Eq. (3). 5 this phenomenology is in direct link with experimen- 1 tally feasible quantum simulators. Recently, Nagy and : coworkers [14] argued that the Dicke model effectively in the study of BEC in multimode cavities [28, 29]. Re- v i describestheself-organizationphasetransitionofaBose- cent works [30, 31] discussed a multimodal-Cavity QED X Einstein condensate (BEC) in an optical cavity [15, 16]. simulator with disordered interactions. The authors ar- r Additionally, Dimer and colleagues [17] proposed a Cav- guethatthesesystemscouldbeemployedtoexplorespin- a ity QED realization of the Dicke model based on cavity- glass properties at the quantum level [30–32]. However, mediated Raman transitions, closer in spirit to the orig- the possible quantum computation applications of this inal Dicke’s idea. Evidence of superradiance in this sys- new class of quantum simulators remain relatively unex- temisreportedin[18]. Animplementationofgeneralized plored. Dicke models in hybrid quantum systems has also been In this Letter, we consider a multimode disordered put forward [19]. More generally, Dicke-like Hamiltoni- Dicke model with finite number of modes. We calcu- ans describe a variety of physical systems, ranging from late exactly (in the thermodynamic limit) the free en- CircuitQED[20–24]toCavityQEDwithDiracfermions ergy of the system at temperature T = 1/β and we in graphene [25–27]. Additionally, disorder and frustra- findasuperradiantphasetransitioncharacterizedbythe tionoftheatom-photoncouplingshaveanimportantrole same free-energy landscape of the Hopfield model [33] in 2 the so-called “symmetry broken” phase, with the typi- mode Dicke Hamiltonian: cal strong-coupling threshold of the Dicke model. From M N N,M the theoretical standpoint, our results generalize to the (cid:88) (cid:88) (cid:88) (cid:16) (cid:17) H = ωa†a +∆ σz+ g˜ a +a† σx, (2) case of quenched disordered couplings the remarkable k k i ik k k i analysis performed by Lieb et al. [11–13]. The choice k=1 i=1 i,k=1 of frozen couplings is compatible with the characteris- effectively modelling quantum light-matter interaction tic time scales involved in light-matter interactions. The of N two-level systems with detuning ∆ and M elec- calculation of the partition function leads us to suggest tromagnetic modes supposed to be quasi-degenerate at variational ground states for the model, which we con- the common frequency ω and with couplings that we √ jecture to be exact in the thermodynamic limit. parametrize for future convenience as g˜ = Ωg / N, ik ik The physical consequences of this analysis are fasci- whereΩistheRabifrequencyandthedimensionlessg ’s ik nating: once the multimode strong-coupling regime is are both atom and mode-dependent. In Cavity QED re- reached and qubit-photon couplings are engineered, it alizations, ω represents the detuning between the cavity should be possible to build a pattern-storing system frequencyandthepumpingfrequencyandcouldbeboth whose underlying dynamics is fully governed by quan- positive or negative. A possible choice of the couplings tumlaws. Moreover,Dickequantumsimulatorshereana- is g =cos(k x ), being k the wave vector of the photon il l i l lyzedmaybesuitabletoimplementspecificoptimization and x the position of i-th atom [31]. i problems, in the spirit of adiabatic quantum computa- We are interested in the thermodynamic properties of tion [4–6]. We point out a non-polynomial optimization this system in the limit M (cid:28) N, and thus in evaluat- problem [4, 5, 34], number partitioning, which could be ing the partition function Z =Tre−βH. This evaluation implemented in a single mode cavity QED setup with can be performed rigorously in the thermodynamic limit controllable disorder. Computing applications based on (N → ∞) using the techniques introduced in Refs. [11– cavity mediated interactions might owns the advantage 13]. We first consider the fully-commuting limit ∆ = 0. tobeaviablewaytogenerateentagledmany-bodystates In this case the evaluation of the partition function is withremarkablescalabilityproperties,asrecentlyshown straightforward(seeSupplementarymaterial)andweob- in Ref. [35]. tain Z = Z Z , where Z is a free boson partition FB H FB Hopfield’s main idea [33] is that the retrieval of stored function and Z is a classical Ising model with local H information, such as memory patterns, may emerge as a quenched exchange interactions of the form: collectivedynamicalpropertyofmicroscopicconstituents (in“fnoerucreodnosr”)wweahkoesneeidnttehrrcooungnhecattiornaisn(in“gsypnhaapssees(”e).ga.rHeerbe-- Jij =−ΩN2 (cid:88)M gikωgjk . (3) bian learning [36, 37]). This is achieved in his model k=1 through a fictitious neuronal-dynamics whose effect is to The physical interpretation of this result is that photons minimize the Lyapunov cost function: mediate long range interactions among the atoms, re- sulting in an atomic effective Hamiltonian described by N 1 (cid:88) E =− T S S , (1) a fully-connected Ising model (see Fig. 1). The role of 2 ij i j the couplings g can be understood from Eq. (1) in the i,j=1 ik context of the Hopfield network. They are the memory whereN isthenumberofneurons,S =1ifthei-thneu- pattern stored in the system. By computing exactly the i ronisactive,and−1otherwise,andthepstoredpatterns free energy of the model, we will show that this inter- ξ(k) = ±1 (k = 1,··· ,p) determine the interconnections pretation stays unaltered in the more complicated case i T through the relation: T = 1/N(cid:80) ξ(k)ξ(k) −pδ . ∆(cid:54)=0. ij ij k i j ij Wenowproceedtotheevaluationofthequantumpar- The analysis in Ref. [33] shows that the long-time dy- titionfunction. WeusethemethodofWangandcowork- namics always converges to one of the p stored patterns, ers [13, 40] (proved to be exact in the thermodynamic i.e. theseconfigurationsaretheglobalminimaofthecost limit for M/N →0 [12]). We introduce a set of coherent function (1). The interpretation of this result is that a states |α (cid:105) with α =x +iy , one for each electromag- suitable choice of the interconnections allows to store a k k k k netic mode k, and we expand the partition function on given number of memory patterns into the neural net- this overcomplete basis: work. Data retrieval is achieved through an algorithm that minimizes the energy function (1). A phase transition to a “complex” phase marks the intrinsic limitation Z =(cid:90) (cid:89)M d2αk Tr (cid:104){α}|e−βH|{α}(cid:105) , (4) on the number of patterns p that can be stored. If p ex- π A k=1 ceeds the critical threshold p ∼ 0.14N many failures in the process of retrieval occur [38, 39]. where Tr is the atomic trace only. The only techni- A In this manuscript we consider the following multi- cal complication is the calculation of the matrix element 3 in (4). This turns out to be equal, apart from non- extensive contributions, to the exponential of the oper- ator in Eq (2) with the replacements a , a† → α , α∗ k k k k [12, 13]. At this stage the trace over the atomic degrees of freedom can be easily performed. The integral over the imaginary parts of α ’s give an overall unimportant k constant. Finally, defining the M-dimensional vectors x = (x ,x ,··· ,x ) and g = (g ,g ,··· ,g ), and 1 2 M i i1 √i2 iM with the change of variables m = x/ N, the partition functionassumesasuitableformforperformingasaddle- point integration, i.e. Z = (cid:82) dMme−Nf(m). Here f is the free energy FIG.2: Orderparameterm(1)inEq. (9)asafunctionof∆in N 1 (cid:88) the ultrastrong coupling regime at T =0. On the theoretical f(m)=βm·m− logG(m,g ), (5) N i curve (dashed line) are superposed the results (green dots) i=1 of a naive numerical optimization algorithm minimizing the with: G(m,g )=2cosh(cid:104)β(cid:0)∆2+Ω2(g ·m)2(cid:1)12(cid:105). ground state energy ansatz in Eq. (11). The parameters of i i the simulation are Ω=4, N =100, M =10. Green dots are The order parameter m describes the superradiant the result of a single disorder realization. This implies the phase transition. Physically, it gives the mean number correctness of our self-averaging hypothesis, which is a good approximation already for N = 100. The classical Hopfield of photons in every mode [41]. Its value is determined model is recovered in the limit ∆→0. by minimizing the free energy in Eq. (5). Solutions of thisoptimizationproblemare,inprinciple,g -dependent, i butinthethermodynamiclimitboththefreeenergyand thesesolutionsandtheirstabilityundertemperaturede- the saddle-point equation are self-averaging [38]. Thus crease. Forthisanalysis,webothconsideredtheHessian we conclude that the free energy and the saddle point matrix ∂2f/∂m ∂m (see Supplementary materials for k l equations are given by its explicit expression) and numerical optimization (Fig- ure 2) . The key point, as mentioned above is that in f(m)=βm·m−(cid:104)logG(m,g)(cid:105)g , this “symmetry-broken” phase the system takes 2M de- Ω2 (cid:28)(g·m)g (cid:29) generategroundstates(aswellasmanymetastablestates m= tanh(βµ(g)) , (6) energetically well separated from the ground states). In 2 µ(g) g other words, also in this fully quantum limit the free- with: µ(g) = (cid:0)∆2+Ω2(g·m)2(cid:1)12 and (cid:104)···(cid:105) represent- energy landscape still closely resembles that of the Hop- g field model [38]. ingtheaverageoverthedisorderdistribution. Eq.(6)re- The ground state solutions have the explicit form: duces to the mean-field equations for the Hopfield model for ∆ → 0 [38]. Thus, ∆ may be intended as a quan- m =m(1)(0,0,··· ,0, 1,0,··· ,0 ). (8) k tum annealer parameter. To fully specify the model, the (cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125) k−1times M−k+1times probability distribution for the couplings is needed. In the following we will assume Equation (6) for the order parameter m(1) reduces to: 2µ(m(1)) = Ω2tanh(cid:0)βµ(m(1))(cid:1), where µ(m(1)) = (cid:89)M (cid:18)1 1 (cid:19) (cid:112)∆2+Ω2(m(1))2. In the zero temperature limit the or- P(g)= δ(g −1)+ δ(g +1) , (7) 2 k 2 k der parameter can be evaluated exactly: k=1 (cid:114) Ω2 ∆2 but we have verified that the results are qualitatively ro- m(1) =± − . (9) 4 Ω2 bust as long as the disorder is not too peaked around zero in accordance with the classical results of Ref. [38] At zero temperature the most interesting state is the To locate the critical point it suffices to expand in Tay- ground state (GS) of the Hamiltonian (2). Inspired by lor series Eqs. (6). As in the conventional Dicke model, the calculation above we propose the variational ansatz a temperature-independent threshold Ω2 =2∆ emerges. for the GS: c For Ω<Ω , the phase transition is inhibited at all tem- c |GS(cid:105)=|α ,α ,··· ,α (cid:105)|spin(α ,··· ,α )(cid:105) , (10) 1 2 M 1 M peratures. Whenever the magnitude of the coupling exceeds this threshold value, the critical temperature is lo- where |α1,α2,··· ,αM(cid:105) is the product of M coherent cated at T =∆/arctanh(cid:0)2∆/Ω2(cid:1). states and the spin part is factorized. The mean value of c Above the critical temperature T the only solution to the energy in this GS is given by: c (6)isaparamagneticstate,withmk =0forallk. Below (cid:68)(cid:112) (cid:69) E (m)=m·m− ∆2+Ω2(g·m)2 . (11) T , differentsolutionsappear. Wenowsetouttoclassify GS c g 4 This expression exactly equals the free-energy computed controllabledisorderasfollows. Numberpartitioningcan previously in the limit β → ∞, which leads us to con- be formulated as an optimization problem [44]: given jecture that our factorized variational ansatz is exact. a set A = {a ,a ,...,a } of positive numbers, find a 1 2 N The quantum phase transition is located at the critical partition, i.e. a subset A(cid:48) ⊂ A, such that the residue: √ (cid:80) (cid:80) coupling Ω=Ω = 2∆, at which the paramagnetic so- E =| a − a | is minimized. A partition c aj∈A(cid:48) j aj∈/A(cid:48) j lution becomes unstable. In the symmetry-broken phase can be defined by numbers S = ±1: S = 1 if a ∈ A(cid:48), j j j we have 2M degenerate ground states of the form S = −1 otherwise. The cost function can be replaced j by a classical spin hamiltonian, whose ground state is |GS(cid:105) =|0,0,··· ,0,±m(1), 0,··· ,0 (cid:105)|spin(±m(1))(cid:105) , k k equivalent to the minimum partition: (cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125) k−1times M−ktimes (12) N (cid:88) with k =1,···M. The spin wave function is also factor- H = aiajSiSj. (14) ized |spin(±m(1))(cid:105) =(cid:81) |s (cid:105) , with i,j=1 k i i k (cid:32) (cid:112) (cid:33) In a single mode cavity QED network couplings can 1 −∆+ ∆2+β2 |s (cid:105) = − ik |e (cid:105)+|g (cid:105) , (13) be chosen as g = cos(kx ) [31]. By the definition of i k N β i i i i ik a = max a and a˜ = a /a, it is possible to engineer A j j j the g ’s in such a way to implement a given instance of whereN anormalization,β =g m(1) and|e (cid:105),|g (cid:105)are i ik ik i i the problem provided that the cavity is in the “blue” σz’seigenstates. Itisworthnotingthat,asexpected,the i detuned regime to ensure the appropriate sign for the ground state energy is a self-averaging quantity, whereas couplings,seeEq. (3). Withasuitableadiabaticanneal- the ground states are not, being disorder-dependent also ingoftheatomicdetuning∆,thesystemshouldcollapse in the thermodynamic limit. on qubit configurations that are good solutions of the The above calculation shows that in the superradiant corresponding optimization problem. phase the ground state of the system is a quantum su- In conclusion, this Letter provides the first rigorous perposition of the 2M degenerate eigenvectors given by analysis of the multimode disordered Dicke model, valid Eqs. (12,13). Their explicit form suggests that at fixed beyond the weak-coupling regime and exact in the ther- disorderandmodenumbertheinformationaboutthedis- modynamic limit. The equivalence between multimodal ordered couplings belonging to the k-th mode is printed disordered Dicke model and a quantum Hopfield net- on the atomic wave function. Moreover, the photonic work [45], together with the proposal of a cavity QED partsofthewavefunctionsareallorthogonalfork (cid:54)=k 1 2 setupimplementinganonpolynomialoptimizationprob- in the thermodynamic limit. This implies that in princi- lem, demonstrates the possibility of quantum computa- ple a suitable measure on the photons-subsystem causes tional abilities of this new class of quantum simulators. the collapse over one of the 2M ground states and gives Our proposal is conceptually complementary to a stan- thus the possibility to retrieve information (“patterns”) dard quantum computation perspective [46, 47]. Indeed, storedintheatomicwavefunction. Asmentionedabove, theinformationcanbe“written”onthequbitsthrougha a single-mode Dicke model has been recently realized quantumannealingonthedetuning∆,similarlytowhat with cavity-mediated Raman transitions in cavity QED happens for adiabatic quantum computation [4–6]. with ultracold atoms [18]. A Multimode cavity QED Acknowledgements.— We are grateful to B. Bassetti, setupsupportingdisorderedcouplingshasbeenproposed S. Mandrà, G. Catelani, M. Gherardi, S. Caracciolo, in refs. [30, 31], and preliminary evidence of superradi- L. Molinari, F. Solgun and A. Morales for useful dis- ance in this system is found in [42]. A setup operating cussions and feedback on this manuscript. GV was sup- in multimode regime has been recently suggested also in portedbyAlexandervonHumboldtfoundationandKnut CircuitQED[43]. Wearenotawareofasetup(inCavity och Alice Wallenbergs foundation. orCircuitQED)thatimplementsbothmultimodestrong coupling regime and controllable disorder, a necessary condition for the quantum pattern-retrieval system that we conjecture here. We surmise that Multimode Dicke quantum setups [1] R. P. Feynman, International Journal of Theoretical with controllable disorder could be used beyond stor- Physics 21, 467 (1982). age, to simulate specific optimization problems. In- [2] M. A. Nielsen and I. L. Chuang, Quantum computation deed, finding the ground state of classical spin models and quantum information (Cambridge university press, 2010). with disordered interactions is equivalent, in most cases, [3] A. Barenco, C. H. Bennett, R. Cleve, D. P. DiVincenzo, to finding solutions of computationally expensive non- N.Margolus,P.Shor,T.Sleator,J.A.Smolin, andH.We- polynomial (NP) problems [34]. For example, the sim- infurter, Phys. Rev. A 52, 3457 (1995). plest NP-hard problem, number partitioning, could be [4] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lund- implemented in a single-mode cavity QED setup with gren, and D. Preda, Science 292, 472 (2001). 5 [5] G. E. Santoro, R. Martonˇák, E. Tosatti, and R. Car, [27] F.M.D.Pellegrino,L.Chirolli,R.Fazio,V.Giovannetti, Science 295, 2427 (2002). and M. Polini, Phys. Rev. B 89, 165406 (2014). [6] V. Bapst, L. Foini, F. Krzakala, G. Semerjian, and [28] S.Gopalakrishnan,B.Lev, andP.Goldbart,Nat.Phys. F. Zamponi, . 5, 845 (2009). [7] I. Buluta and F. Nori, Science 326, 108 (2009). [29] S.Gopalakrishnan,B.Lev, andP.Goldbart,Phys.Rev. [8] I.M.Georgescu,S.Ashhab, andF.Nori,Rev.Mod.Phys. A 82, 043612 (2010). 86, 153 (2014). [30] S.Gopalakrishnan,B.Lev, andP.Goldbart,Phys.Rev. [9] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, Lett. 107, 277201 (2011). A. Sen, and U. Sen, Advances in Physics 56, 243 (2007). [31] P. Strack and S. Sachdev, Phys. Rev. Lett. 107, 277202 [10] R. Dicke, Phys. Rev. 93, 99 (1954). (2011). [11] K. Hepp and E. Lieb, Ann. Phys. 76, 360 (1973). [32] P. Rotondo, E. Tesio, and S. Caracciolo, Phys. Rev. B [12] K. Hepp and E. Lieb, Phys. Rev. A 8, 2517 (1973). 91, 014415 (2015). [13] Y.K.WangandF.T.Hioe,Phys.Rev.A7,831(1973). [33] J. J. Hopfield, PNAS 79, 2554 (1982). [14] D.Nagy,G.Końya,G.Szirmai, andP.Domokos,Phys. [34] A. Lucas, Frontiers in Physics 2 (2014), Rev. Lett. 104, 130401 (2010). 10.3389/fphy.2014.00005. [15] K.Baumann,C.Guerlin,F.Brenneke, andT.Esslinger, [35] C.Aron,M.Kulkarni, andH.E.Tu¨reci,ArXive-prints Nature (London) 464, 1301 (2010). (2014), arXiv:1412.8477 [quant-ph] . [16] K.Baumann,R.Mottl,F.Brennecke, andT.Esslinger, [36] D. O. Hebb, Archives of Neurology and Psychiatry 44, Phys. Rev. Lett. 107, 140402 (2011). 421 (1940). [17] F. Dimer, B. Estienne, A. S. Parkins, and H. J. [37] D. O. Hebb, Brain Mechanism and Learning (1961). Carmichael, Phys. Rev. A 75, 013804 (2007). [38] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Phys. [18] M. P. Baden, K. J. Arnold, A. L. Grimsmo, S. Parkins, Rev. A 32, 1007 (1985). and M. D. Barrett, Phys. Rev. Lett. 113, 020408 (2014). [39] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Phys. [19] L. Zou, D. Marcos, S. Diehl, S. Putz, J. Schmiedmayer, Rev. Lett. 55, 1530 (1985). and P. Rabl, Phys. Rev. Lett. 113, 023603 (2014). [40] J.LarsonandM.Lewenstein,NewJournalofPhysics11, [20] Y. Zhang, L. Yu, J.-Q. Liang, G. Chen, S. Jia, and 063027 (2009). F.Nori,ScientificReports4,4083(2014),arXiv:1308.3948 [41] C. Emary and T. Brandes, Phys. Rev. Lett. 90, 044101 [quant-ph] . (2003). [21] P. Nataf and C. Ciuti, Nature Communications 1, 72 [42] A. Wickenbrock, M. Hemmerling, G. R. M. Robb, (2010). C. Emary, and F. Renzoni, Phys. Rev. A 87, 043817 [22] O. Viehmann, J. von Delft, and F. Marquardt, Phys. (2013). Rev. Lett. 107, 113602 (2011). [43] D. J. Egger and F. K. Wilhelm, Phys. Rev. Lett. 111, [23] J. A. Mlynek, A. A. Abdumalikov, C. Eichler, and 163601 (2013). A. Wallraff, Nature Communications 5, 5186 (2014). [44] S. Mertens, Phys. Rev. Lett. 81, 4281 (1998). [24] A. Mezzacapo, U. Las Heras, J. S. Pedernales, L. Di- [45] Y. Nonomura and H. Nishimori, eprint arXiv:cond- Carlo, E. Solano, and L. Lamata, Scientific Reports 4, mat/9512142 (1995), cond-mat/9512142 . 7482 (2014). [46] D. Loss and D. P. DiVincenzo, Phys. Rev. A 57, 120 [25] L. Chirolli, M. Polini, V. Giovannetti, and A. MacDon- (1998). ald, Phys. Rev. Lett. 109, 267404 (2012). [47] R. Raussendorf, D. Browne, and H. Briegel, Phys. Rev. [26] D. Hagenmu¨ller and C. Ciuti, Phys. Rev. Lett. 109, A 68, 022312 (2003). 267403 (2012). 6 Dicke simulators with emergent collective quantum computational abilities supporting material Pietro Rotondo, Marco Cosentino Lagomarsino, and Giovanni Viola In this material, we give more details on the derivations of the results presented in the main text. DERIVATION OF EQUATION (3) We begin with the derivation of Eq. (3), of the main text. We consider the partition function Z =Tre−βH with H given in Eq. (2) for ∆ = 0. In this fully commuting limit we can evaluate the partition function straightforwardly. We introduce new set of bosonic operators: N N Ω (cid:88) Ω (cid:88) b† =a† + √ g σx, b =a + √ g σx. (S.1) k k ω N ik i k k ω N ik i i=1 i=1 with [b ,b†]=δ . By means of those, H is written as the sum of two commuting operators: k(cid:48) k k,k(cid:48) (cid:88)M Ω2 (cid:88)N (cid:88)M H =ω b†b − g g σxσx (S.2) k k Nω ik jk i j k=1 i,j=1k=1 As a byproduct we obtain the factorization of the full partition function [S1]:Z =Z Z , where Z is an overall free B H B boson partition function that we can safely ignore in the thermodynamic limit. On the other hand   (cid:88) ZH =Trσexpβ Jijσiσj (S.3) i,j is an Ising contribution with both spin and mode dependent couplings J of the form given in Eq. (3) of the main ij text. In Eq.(S.3) Tr indicates the trace on the spins only. σ DERIVATION OF EQUATION (5) InthissectionwereportthederivationofEq. (5)ofthemaintext,whichessentiallyfollowsthederivationofWang and Hioe [S2] proved to be rigorous by Hepp and Lieb [S3]. The authors of Ref. [S2] have shown explicitly, in the termodynamic limit, that the convenient way to calculate the trace on the Hilbert space of bosons, in the partition function, is to evaluate it on the set of the coherent states |{α}(cid:105). The photonic matrix element in the partition function of Eq. (4) equals in the thermodynamic limit (α =x +iy ): k k k (cid:32) M N N M (cid:33) (cid:88) (cid:88) Ω (cid:88)(cid:88) (cid:104){α}|e−βH|{α}(cid:105)(cid:39)exp −βω (x2 +y2)−β∆ σz−β√ g x σx . (S.4) k k i ik k i N k=1 i=1 i=1k=1 The atomic trace thus factorizes and it can be calculated:  (cid:118)  Z =(cid:90) (cid:89)M dxkπdyk e−βω(cid:80)Mk=1(x2k+yk2)(cid:89)N coshβ(cid:117)(cid:117)(cid:116)∆2+ ΩN2 (cid:32)(cid:88)M gikxk(cid:33)2=(cid:90) (cid:89)M dmπke−Nf(m,β). (S.5) k=1 i=1 k=1 k=1 Inthelasttermweintroducedthevectorialnotationdefinedinthemaintext. Thefinalexpressionforthefreeenergy is: f(m,β)=βm·m− 1 (cid:88)N logcosh(cid:32)β(cid:114)∆2+ Ω2 (g ·m)2(cid:33) , (S.6) N N i i=1 7 By minimizing the free energy above and using the self averaging property of Eq. (S.6), we obtain the exact mean field equations: Ω2 (cid:28)(g·m)g (cid:29) m= tanh(βµ(g)) (S.7) 2 µ(g) g To locate the critical point it suffices to expand in Taylor series Eqs. (S.6),(S.7): f(m)−f(0)=β(cid:18)1− Ω2 tanh(β∆)(cid:19)m·m+O(cid:0)m4(cid:1) , m = Ω2 tanh(β∆)m +O(m3). (S.8) 2∆ k k 2∆ k k As in the conventional Dicke model, a temperature-independent threshold Ω2 =2∆ emerges. For Ω<Ω , the phase c c transition is inhibited at all temperatures. Whenever the magnitude of the coupling exceeds this threshold value, the critical temperature is located at T =∆/arctanh(cid:0)2∆/Ω2(cid:1). Solutions to Eq. (S.7) can be classified according to the c number of non-zero components n of the order parameter m [S4]: m =m(n)(1,1,··· ,1,0,0,··· ,0), (S.9) n (cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125) ntimes M−ntimes where all permutations are also possible. In particular, solutions for n = 1 are the ones with the lowest free energy. There are 2M of such solutions, corresponding to the Z ×S symmetry breaking of our multimode Dicke model. 2 M For completeness we report the explicit expression for the Hessian matrix of the free energy, omitted in the main text for space imitations: ∂2f (cid:28)g g (cid:29) (cid:28)(g·m)2g g (cid:20)tanh(βµ(g)) (cid:21)(cid:29) =2δ −Ω2 k l tanh(βµ(g)) +Ω4 k l +βsech2(βµ(g)) . (S.10) ∂m ∂m kl µ(g) µ(g)2 µ(g) k l g g [S1] P. Rotondo, E. Tesio and S. Caracciolo, Phys. Rev. B 91 014415 . [S2] K. Wang and F. T. Hioe, Phys. Rev. A 7 831. [S3] K. Hepp and E. Lieb, Phys. Rev. A 8, 2517. [S4] D. J. Amit, H. Gutfreund and H. Sompolinsky, Phys. Rev. A 32 1007.