Exactness of mean-field equations for open Dicke models with an application to pattern retrieval dynamics

Open quantum Dicke models are paradigmatic systems for the investigation of light-matter interaction in out-of-equilibrium quantum settings. Albeit being structurally simple, these models can show intriguing physics. However, obtaining exact results on their dynamical behavior is challenging, since it requires the solution of a many-body quantum system, with several interacting continuous and discrete degrees of freedom. Here, we make a step forward in this direction by proving the validity of the mean-field semi-classical equations for open multimode Dicke models, which, to the best of our knowledge, so far has not been rigorously established. We exploit this result to show that open quantum multimode Dicke models can behave as associative memories, displaying a nonequilibrium phase transition towards a pattern-recognition phase.

Open quantum Dicke models are paradigmatic systems for the investigation of light-matter interaction in out-of-equilibrium quantum settings.Albeit being structurally simple, these models can show intriguing physics.However, obtaining exact results on their dynamical behavior is challenging, since it requires the solution of a many-body quantum system, with several interacting continuous and discrete degrees of freedom.Here, we make a step forward in this direction by proving the validity of the mean-field semi-classical equations for open multimode Dicke models, which, to the best of our knowledge, so far has not been rigorously established.We exploit this result to show that open quantum multimode Dicke models can behave as associative memories, displaying a nonequilibrium phase transition towards a pattern-recognition phase.
Since its inception [1], the Dicke model has become a paradigm for the study of light-matter interaction and its equilibrium as well as isolated-system dynamical properties have been widely investigated both theoretically and experimentally [2][3][4][5][6][7][8][9][10][11][12][13][14].Nowadays, the interest is in understanding how the presence of an environment, leading to dissipative effects, impacts on the behavior of Dicke models.In this out-of-equilibrium setting, much less is known.Several arguments indicate the persistence of the Dicke superradiant phase transition [15][16][17][18], and this hypothesis is further supported by numerical [19] and experimental [20] evidence.
Particularly intriguing is the possibility that these nonequilibrium spin-boson systems can feature dynamics akin to associative memories [21,22], i.e. they can display pattern-recognition behavior [23][24][25][26][27][28], and implementations of this physics are being explored in realistic experimental setups [29].Couplings between spins and bosons encode different patterns which, in the simplest case, are strings of ±1, see Fig. 1(a).The overlap ξ µ of the spin configuration with pattern µ, which plays the role of an order parameter, is defined by means of a generalized magnetization [c.f.Fig. 1(a)].Assuming the initial configuration to be close to one pattern, two different regimes may emerge.In the first, the state converges -due to dissipation-to a stationary one where all information about the initial time is lost.As sketched in Fig. 1(b), this coincides with a regime where the overlaps ξ µ are all zero.In the other, instead, it converges to a stationary state displaying a finite overlap with the initially stored pattern.In this case, the system "recognizes" the initial condition as a pattern and stores this information in its nonequilibrium steady state.In Dicke models, the observed stationary regime is expected to depend on the spin-boson coupling strength, see Fig. 1(b).
Understanding whether this pattern-recognition be-   Each pattern is associated with a mode.The overlap of the quantum state with the patterns is defined as a generalized magnetization aligned with the coefficients G µ,k .b) As a function of the spin-boson coupling strength, the quantum system passes from a disordered phase, in which it cannot store any pattern, to an "ordered" one, in which it can recognize and protect a pattern.
havior corresponds to a genuine nonequilibrium phase requires the study of quantum systems with large number of bosons and spins.Simulations in fully quantum regimes beyond perturbative approaches [28][29][30] are thus infeasible.Analytically, one may study these sytems relying on so-called mean-field equations, obtained by assuming that expectation values of products of operators factorize [15,19,31].However, a proof of the validity of this assumption in nonequilibrium open Dicke models is still missing, and a widespread belief is that a "full quantum treatment" may lead to different results.
In this paper, we provide a proof of the exactness of the mean-field assumption for open multimode Dicke models.This result is relevant as it solves an open question on the validity of the semi-classical treatment for these systems.Further, it allows us to establish the existence of a nonequilibrium pattern-recognition phase transition in Dicke models.Our proof -which arXiv:2009.13932v2[cond-mat.stat-mech]14 Jul 2021 takes inspiration from Ref. [32]-is of broad applicability: it can be adapted to account for the presence of individual spin dissipative processes [19], to account for time-dependent coefficients in the generator [33][34][35], or even to other models with all-to-all couplings [36][37][38][39][40][41], also with multi-body interactions [42][43][44], Open multimode Dicke models.-OurDicke model consists of an ensemble of N sp spins coupled to M different bosonic modes, described by annihilation and creation operators a µ , a † µ obeying canonical commutation relations [45].Spins are two-level systems, with excited state |• and ground state |• .Transitions between states in the k-th spin are implemented by the Pauli operator σ x .The (Markovian) nonequilibrium dynamics of the spinboson model is implemented by the Lindblad generator Ẋ = L[X] [46][47][48], providing the time-evolution of a generic operator X. Defining n µ = a † µ a µ , we consider the second term appearing on the right-hand side describes boson losses, at rate κ µ for the different modes, while H is the system Hamiltonian.This operator consists of a free contribution for both spins and bosons and of an interaction term The coefficients G µ,k specify the spin-boson interaction.
We consider these to be independent identically distributed random variables assuming the values +1 or −1 with equal probability, as sketched in Fig.
(1)(a).The scaling 1/ N sp -typical for these models-is important to establish a well-defined thermodynamic limit [15] (see also [49] for an application to open systems).For each µ, the string G µ,k forms a pattern which is encoded in the Hamiltonian.A key result of this paper consists in showing that the system can recognize and protect an initially stored pattern, for strong enough spin-boson coupling |g 0 |, see Fig. 1.
Before showing this, we make some considerations which bring the model into a convenient form, see Fig. 2. First, without loss of generality, the first pattern, G 1,k , which is made of ±1, can be brought into a pattern with k < l a t e x i t s h a 1 _ b a s e 6 4 = " s u 9 z 8 2 b z h I all +1, by means of the gauge transformation σ z → −σ z applied to those spins h for which, originally, G 1,h = −1.
Then, we reorder the remaining M − rows of G µ,k .We look at G 2,k : this has ±1 at random positions.We now relabel the spins.We take those with G 2,h = +1 to the left and those with G 2,h = −1 to the right.This reshaping is not affecting the first pattern.In addition, there We then move to G 3,k and we relabel spins as follows.In the subset of spins for which G 2,k = 1, we have values of G 3,k which can be both positive and negative.We thus reorder this subsequence in such a way that all +1 are moved on the left and −1 on the right.The same can be done for the subset of the sequence G 3,k corresponding to values G 2,k = −1.This procedure, sketched in Fig. 2, is then iterated up to the last pattern.This mapping generates 2 M −1 subsets of spins, described by "large-spin" operators and interacting with the bosonic modes.For N sp 1, these subsets are expected to have the same number of spins.This is due to the fact that, given the statistical properties of the G µ,k , in a large enough set of randomly chosen spins there is, at leading order in extensivity of the set, an equal number of +1 and of −1, in their G µ,k .We can thus consider subsets to contain N = N sp /2 M −1 spins.In this representation, the interaction Hamiltonian reads where is the sum of the σ a -spin operators which belongs to the k-th subset, denoted as Γ k [see Fig. 2(a)].In addition, we have defined g = g 0 / √ This representation provides a more compact formulation of the model.This mapping can be extended to consider models whose spin-only part of the dynamical generator is not invariant under the gauge transformation or also to consider generic distributions for G µ,k [50].
Mean-field dynamics.-Asa consequence of the previous mapping, it is sufficient for understanding the behaviour of our nonequilibrium Dicke model to focus on the dynamics of the "large-spin" operators.In this representation, the generator is L N , the same as the one in Eq. ( 1) with Hamiltonian rewritten as where the functional ω represents the initial state, while ω t the time-evolved one.As a consequence, we have We are interested in the "macroscopic" operators [51-54] the first ones are the usual average "magnetization" operators of the spin ensembles, while the rescaled bosonic operators appear typically in superradiant transitions.Indeed, a non-vanishing expectation of these operators implies a macroscopic (∝ N ) bosonic occupation.We want to derive the dynamics of these quantum operators in the thermodynamic limit N, N sp → ∞.We thus compute the action of the generator L N on the operators in Eq. ( 5) and get [50] where abc is the fully anti-symmetric tensor.To make progress, one typically assumes that the dynamics does not generate correlations among the different constituents in the thermodynamic limit, so that expectation values factorize.This leads to the mean-field equations In order to show that they are exact in the thermodynamic limit, we need to prove that lim meaning that the expectation of the operators of Eqs. ( 5) behaves, for large N , as the time-dependent scalar functions m a,k (t), α µ (t) obeying Eqs.(7).To obtain this result, a proper strategy must be identified.In particular, an appropriate "cost function" controlling the above limits is needed.
This quantity is a sum of positive contributions consisting of the expectation of the square of the distance of the operators from their mean-field counterpart.Namely, E N (t) measures the fraction of spins or bosons not behaving as dictated by Eqs.(7).In addition, via Cauchy-Schwarz inequality, one can show that and thus lim N →∞ E N (t) controls the limits in Eq. ( 8), as desired.For physical initial states [52][53][54], with short-range correlations, one has lim N →∞ E N (0) = 0.As we now show, for these states, E N (t) vanishes for large N , implying the exactness of the mean-field assumption for these nonequilibrium multimode Dicke models.
Theorem.With the above definitions, if the initial state of the system is such that lim N →∞ E N (0) = 0 then, for all finite t, we have that lim N →∞ E N (t) = 0. Proof: The full proof is reported in Ref. [50].Here we provide the main steps.The idea is to use Gronwall's Lemma [55,56], which states that if a positive, bounded, and N -independent constant C, such that ĖN (t) ≤ C E N (t), exists then With the assumption lim N →∞ E N (0) = 0, letting N → ∞ in the above relation would prove the theorem.What is missing is to show that such constant C indeed exists.This can be achieved by directly inspecting the time derivative of all terms forming E N (t).They are given by sums of contributions having, for instance, the form ω t (E b,k BA µ ), where B can either be an operator or a scalar from Eqs. (7).In addition, it can be shown that and this gives a way to estimate a suitable constant C.
We thus obtain and we can exploit Gronwall's Lemma to finish the proof of the theorem as already discussed.
Pattern-recognition phase transition.-With the above result, we establish that the semi-classical meanfield equations (7) correctly capture the behavior of our system, in the thermodynamic limit.As such, we can now use these equations to unveil the presence of a nonequilibrium pattern-recognition phase transition.
In the original formulation of the problem, see Eq. ( 2) and Fig. 1(a), we can define the overlap of the quantum state of the spins with the pattern µ as This equation shows that, if the expectaction value of the operator σ z is, for each spin, aligned with the corresponding value of G µ,k , then the overlap |ξ µ | is different from zero (pattern retrieval).Otherwise, ξ µ tends to vanish for N sp → ∞ (pattern not retrieved).In the large-spin representation, the overlaps can be expressed in terms of the coefficients f M µ,k and of the macroscopic operators m N z,k , [c.f.Eq (3) and Fig. 2].In particular, Invoking our theorem, we can thus study the dynamics and the stationary properties of these overlaps through the scalars m z,k , obeying the mean-field equations (7).
To prove the existence of the phase transition, we first show the presence of different stationary solutions to Eq. ( 7), featuring a finite overlap with one of the patterns.Without loss of generality, we consider all rates of the dynamical generator to be positive and, further, that the constant of motion m 2 T,k = a m 2 a,k = 1, ∀k.Then, we take the ansatz solution m z,k = f M ν,k |z|, aligned with pattern ν, and look for conditions ensuring its existence as a stationary solution for Eqs.(7).Note that such ansatz has indeed a finite overlap with pattern ν, since ξ ν = |z| while ξ µ = 0 ∀µ = ν, and that also m z,k = −f ν,k |z| would be valid, with ξ ν = −|z|.
By substituting the ansatz for m z,k in Eqs.(7), taking m y,k = 0 and appropriately fixing the values of m x,k (see Ref. [50] for details) we find that the relation must be satisfied, in order for the assumed stationary solution to exist.This is not always the case; indeed, |z| must be a positive real number, |z| ∈ [0; 1], and this only happens if the argument of the square root is positive.This observation yields a critical value, such that for g 0 ≥ g crit the ansatz solution exists, with |z| given by Eq. ( 12).On the other hand, if g 0 < g crit , we can only have |z| = 0, and we are outside the patternrecognition phase.The critical g depends on the pattern through the parameters Ω ν , κ ν , see also Fig. 3(a-b).
Further, note that a finite stationary overlap corresponds to a macroscopic occupation of the associated bosonic mode.Our theorem indeed implies Discussion.-We have derived two key results for multimode Dicke models.First, we have shown that the mean-field assumption, typically exploited to consider the large-scale behavior of these systems, actually provides an exact description in the thermodynamic limit.Second, we have used this new insight to reveal the presence of a nonequilibrium phase transition from a disordered phase to a pattern-recognition phase in open multimode Dicke models.The stability of stationary solutions, such as the one of Eq. ( 12), for open Dicke models has been shown, for instance, in Refs.[15,19].For the multimode settings investigated here, the agreement of our numerical results with analytical ones [c.f.Fig. 3] suggests that the proposed stationary states, having finite overlap with the patterns, possess stable basins of attraction in the pattern-recognition phase.Interestingly, the critical spin-boson coupling strength depends on the specific pattern through the corresponding bosonic mode parameters.This may allow for intermediate regimes of pattern recognition, where only certain patterns can be stored and retrieved.Following Ref. [37], we remark that the validity of the semi-classical Eqs. ( 5) provides a necessary ingredient to obtain mathematically rigorous results on quantum fluctuations.It would be interesting to exploit it to re-obtain bosonic descriptions [57,58] employed for the investigation of quantum fluctuations in closed Dicke models and to extend these to open systems, via quantum central limit theorems [37,53,59].Contrary to Holstein-Primakoff approximations, these procedures do not assume a conserved total spin operator and are thus more general [15].

Exactness of Mean-Field Equations for Open Dicke Models with an Application to Pattern Retrieval Dynamics
Federico Carollo, 1 and Igor Lesanovsky As briefly mentioned in the main text, the mapping that we have introduced is applicable to more general models than the one we have focussed on in this work.We discuss here two possible extensions.
The first step of the mapping, as presented in the main text, is the gauge transformation σ z → −σ z .This transformation may also, in general, modify the spin Hamiltonian (or a possible dissipative contributions on the spins) or may lead to complications.However, such a step is not necessary.Indeed, instead of applying the gauge transformation to the spins in order to make the first pattern uniform (G 1,k = 1, ∀k), one can directly start acting on the first pattern reordering it by moving all spins with G 1,k = +1 to the left and all spins with G 1,k = −1 to the right.After this is done, one can then proceed to reorder analogously the other patterns.In this way, instead of 2 M −1 subsets of spins, the mapping generates 2 M of them.An illustration of this version of the mapping is given in Fig. S1.+1 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " g n K h h 9  In this version of the mapping, we do not make the initial gauge transformation but rather start reordering the first pattern G 1,k .We relabel spins in such a way that all G 1,k = 1 appear first, and are followed by all the G 1,k = −1.This merely amounts to a relabeling of the spins.We then proceed to reorder also the other patterns in an analogous way to what done in the main text.We notice that, in this way, the mapping produces 2 M subsets of spins Γ k , which are equally coupled to each mode, (Γ1 and Γ8 are highlighted in the figure for clarity).b) These subsets of spins couple to bosons as specified by the matrix f M µ,k .
Another possible extension considers different probability distributions for the coupling coefficients G µ,k .In general, each G µ,k may assume the value x i , where i = 1, 2, . . ., d, with probabilities p i .In this case, one can reorder the first pattern by relabelling all those spins k, associated with a G 1,k = x 1 to move them to the first positions, followed by those with G 1,k = x 2 , and so on till the spins having G 1,k = x d have been accomodated.As for the mapping in the main text, this is just a simple relabelling of the spins.Then, moving to the second pattern G 2,k , we can proceed as follows.We focus, one by one, on the subsets of spins k for which G 1,k = x i .Within these subsets we further reorder the spins, relabelling them in such a way that the ones having G 2,h = x 1 are moved to the first positions in the subset, followed by those with G 2,h = x 2 , till those having G 2,h = x d .For the third, as well as for the remaining patterns the spins are analogously reordered.
Mapping to large spins for general distributions.We present here an example of the mapping for couplings G µ,k which can take three possible values x1, x2, x3 and for M = 2 bosonic modes.The idea is to relabel the spins in such a way that the two patterns G 1,k and G 2,k are reordered as in the figure .In particular, the first pattern will be reordered in such a way that the coupling x1 is associated with the first subset of spins, x2 with the second and x3 with the third.Each of these subsets can be partitioned in smaller subsets.For instance, spins which couple through x1 with the first mode, can be reordered according to their coupling with the second mode.The reordering of this last pattern reveals 3 2 ensembles of spins which are all coupled with both bosonic modes in the same way.Since the probabilities of having the different values xi for the couplings are not uniform, the extensivity of the different ensembles is modulated by the probabilities pi as explained in the main text.
We notice that in this case, the procedure generates d M subsets of spins which are, within each subset, all equally coupled to the different bosonic modes.As such, these ensembles can be treated as collective large spins.Interestingly, the extensivity of these subsets is, in general, not uniform and in fact depends on the probabilities p i of extracting the different couplings.In particular, we have that if the spins in the subset Γ k interact with the µ bosonic mode with the coupling x iµ , then, for large numbers of total spins, N sp 1, the number of spins forming the subset Γ k is given by N k ≈ p i1 p i2 . . .p i M N sp .When all subsets are expected to be different, the proof that we present is still applicable; it is sufficient to redo the same steps using N sp instead of the N introduced in the text, also in the definition of the macroscopic operators (5).The extensivity of the different ensembles will then pose constraints on the modulus of the expectation values of the different macroscopic operators.For instance, in a given ensemble Γ k made of N k ≈ p i1 p i2 . . .p i M N sp spins, expectation values of macroscopic spin operators cannot be larger than p i1 p i2 . . .p i M , in the thermodynamic limit.

PROOF OF MAIN THEOREM
Theorem 1.Given the sequence of generators L N introduced in the main text, the sequences of operators in Eq. (5), and the scalar time-dependent functions appearing in Eqs.(7), the error E N (t), defined as is satisfied.Eq. (S2) implies that the mean-field equations for the scalar quantities α µ and m a,k , for all µ, k and a = x, y, z, correctly describe the dynamics of the expectation values of the limiting operators of the sequences in Eqs.(5).
Proof.In order to prove the theorem, we start considering the time-derivative of E N (t) (S4) Using the results in Lemmata 3 and 4 to bound the modulus of all terms appearing in the sums, we find where c 0 , d 0 are the time-independent and N -independent bounded positive quantities defined in Lemmata 3 and 4.
The above inequality implies, because of Gronwall's Lemma, If the initial state of the system is such that condition in Eq. ( S3) is satisfied, we have lim Now, because of the bounds the statement in Eq. (S2) implies that lim Physically, these relations mean that the mean-field dynamical equations (7) correctly capture the time-evolution of macroscopic spin and boson operators, in the limit N → ∞.
LEMMATA Lemma 1.Given the sequences of operators defined in Eq. ( 5), the scalar time-dependent functions of (7), and the quantity E N (t), defined in Eq. ( 9) and given explicitely in Eq. (S1), the following bounds hold: Proof.Before considering all different cases, we derive a bound which is valid for generic operators.We will then show, case by case via an appropriate choice of the operators, each of the above relations.
We start considering the expectation value ω t A † CB .The state ω t is a positive linear and normalized functional.Thus, we can use Cauchy-Schwarz inequality to obtain In addition, the inequality Altogether, we have To proceed we make the following observation.Consider two positive numbers x, y: the square of their difference is a positive number By turning the above relation around, we get where the second inequality is obtained by removing the factor 1/2. Identifying x = ω t (A † A) and y = ω t (B † B), this means that Using the above finding in Eq. (S13), we have the relation which we use to show all relations in the statement of the Lemma.
To prove the Eq.(S8), we consider the bound in Eq. (S14) with , and then add on the right hand side of the resulting Eq. (S14) all the missing terms to reconstruct E N (t).This can be done since each term forming E N (t) is positive.Notice that if b = a and k = h, then equation (S8) is trivially satisfied.
To prove Eq. (S9), we consider Eq. (S14) with C = m N b,k , noticing that m N b,k = 1.Then, we take † and B = m N a,h − m a,h (t), and add on the right hand side of the resulting Eq. (S14) the remaining terms.
Lemma 2. Given the sequences of operators m N a,k and α µ,N defined in Eq. (5), and the sequence of dynamical generators L N introduced in the main text, we have that Proof.To obtain the relations above, it is sufficient to compute the action of the generator on the considered operators, using their definition and exploiting the commutation relations of the bosonic operators and the algebraic rules The tensor abc is the fully antisymmetric tensor, while the Kronecker delta appears because operators of the subset k commute with operators of the subset h.
Lemma 3. Given the sequences of operators m N a,k defined in Eq. ( 5), the sequence of dynamical generators L N introduced in the main text and the scalar time-dependent functions in (7), we have that where Proof.First of all, we recall that the scalar functions m a,k , α µ are solution to the mean-field equations (7), which are obtained via the factorization assumption of expectation values of the Heisenberg equations emerging from the results presented in Lemma 2. With this in mind, we start defining and consider explicitly the time-derivative: The first term on the right-hand side of the above equation is obtained by taking the derivative on the state functional ω t and using Eq. ( 4); the second term is, instead, emerging from the time-derivative applied to the term in the square brackets of Eq. (S15).
Focussing on the action of the Lindblad map on the operator, we have that

(S16)
Notice that the terms m a,k (t), ṁa,k (t) can be safely pulled inside or outside of the state expectation, since they are scalar quantities.Collecting the first two terms of the above equation, as well as the second two terms, we obtain Since the second term of the above equation is the complex conjugate of the first one, we just focus on the latter.We define it as Exploiting Lemma 2 and the differential equations (7), we have that (we leave summation indeces implicit) (S17) We need to reshape the last contributions to the above equation in a way that we can rewrite them in terms of the difference between macroscopic operators and their corresponding scalar values.To this end, we consider that The task is now to find proper bounds for each of these terms.We consider that for the expectation value in the first and the last term we can use Eq.(S8) in Lemma 1, while for the expectation value in the second term we use Eq.(S9) and for the third Eq.(S10).Considering also the extensions of the summations, we obtain the following bound In the above relation we have introduced the quantity β t defined as for t ≥ 0, needed to bound the scalar term in the fourth term of Eq. (S19).To achieve a meaningful bound we need to show that β t is finite.To this end, we take advantage of the formal solution of the mean-field equations (7) to get α µ (t) = e −(iΩµ+κµ/2)t α µ (0) − ig t 0 ds e −(iΩµ+κµ/2)(t−s) We then take the modulus of the above relation where we have further considered that g = g 0 / √ 2 M −1 .This relation can be satisfied only when the right hand side is positive.For concreteness, we take s = 1.This means that the proposed solution is possible only if

z i < l a t e x i t s h a 1 _
b a s e 6 4 = " A 4 t d + H o p f B V P O 3 f B p C T R 7 Q K 7 M f Y = " > A A A C L 3 i c b Z D L S g M x F I Y z 3 q 2 3 q k s 3 w S p U 0 D I j i m 4 E Q V B X o m C r 0 N Q h k 6 Y 1 T J I Z c h H r M G / k x l d x I 6 K I W 9 / C 9 L L Q 6 g + B n + + c w 8 n 5 o 5 Q z R j D E 0 e F 3 Y x z x L F s c w q R Z m 2 B w 4 e b r B x v 5 E j 1 6 F p Y L P k V v y f 4 1 w Q D U w I D n Y f F Z 9 R M i B V U G s K x 1 v X A T 0 0 j w 8 o w w m l e Q F b T F J M Y t 2 n d W Y k F 1 Y 2 s d 2 8 O 1 x 1 p w l a i 3 J M G 9 u j P i Q w L r T s i c p 0 C m 1 s 9 X O v C / 2 p 1 a 1 r 7 j Y z J 1 B o q S X 9 R y 3 J o E t g N D z a Z o s T w j j O Y K O b + C s k t V p g Y F 3 H B h R A M n / z X 1 L Y r w U 5 l 9 2 K 7 d O g P 4 p g C K 2 A V l E E A 9 s A h O A X n o A o I e A T P 4 A 2 8 e 0 / e i / f h f f Z b R 7 z B z D L 4 J e /r G 8 3 a q X s = < / l a t e x i t > z < l a t e x i t s h a 1 _ b a s e 6 4 = " g 3 S x e e + M a U B n e W S W j 0 3 F H s p 3F F k = " > A A A B 8 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L D a C p 5 I U R Y 8 F L x 4 r 2 F Z o Q t l s N + 3 S 3 U 3 Y 3 Q g 1 9 G 9 4 8 a C I V / + M N / + N 2 z Y H b X 0 w 8 H h v h p l 5 U c q Z N p 7 3 7 Z T W 1 j c 2 t 8 r b l Z 3 d v f 2 D 6 u F R R y e Z I rR N E p 6 o h w h r y p m k b c M M p w + p o l h E n H a j 8 c 3 M 7 z 5 S p V k i 7 8 0 k p a H A Q 8

FIG. 1 .
FIG.1.Pattern recognition in Dicke models.a) Patterns -strings of ±1-are encoded in the couplings G µ,k between Nsp spins and M bosonic modes.Each pattern is associated with a mode.The overlap of the quantum state with the patterns is defined as a generalized magnetization aligned with the coefficients G µ,k .b) As a function of the spin-boson coupling strength, the quantum system passes from a disordered phase, in which it cannot store any pattern, to an "ordered" one, in which it can recognize and protect a pattern.

FIG. 2 .
FIG. 2. Mapping to large spins.a) Example of the mapping for M = 3 patterns and Nsp = spins.The original coupling between the µ-th mode and the k-th spin is encoded in G µ,k .To perform the mapping we first apply a gauge transformation making G 1,k = 1, ∀k.Then, we reorder G 2,k to put all +1 first.Finally, also the last pattern is reordered by moving the +1 towards the right and the −1 towards the left in each sub-block identified by the new G 2,k .In this way, 2 M −1 subsets of spins Γ k , equally coupled with each mode, are identified (Γ1 and Γ4 are highlighted in the figure for clarity).b) These subsets of spins are described by "large-spin" operators and couple to bosons as specified by the matrix f M µ,k .

2 ,FIG. 3 .
FIG.3.Pattern-recognition phase transition.Comparison between theoretical prediction (solid lines) and numerical simulations of the mean-field equations (circles).We consider M = 4. a) Each curve corresponds to the stationary overlap |ξµ| computed from the initial condition ξµ = 1 as a function of g0, for Ω = 0.5, Ωµ = Ω(µ + 2).Rates are in units of κ.Different colors correspond to values of µ growing as indicated by the arrow.Both theoretical and numerical results display a nonequilibrium phase transition, as shown by the behavior of the overlap.b) Same parameters and same order for the curves as in a).The occupation of the µ-th bosonic mode becomes macroscopically occupied when the corresponding pattern is stored in the stationary state.

H 8 P
U p / J 9 0 y p Z d t S 5 v y 6 U m W s W R B y f g F J w D G 9 R B E 9 y A F m g D A m L w A J 7 A s 6 G M R + P F e F 2 2 5 o z V z D H 4 A e P t E 6 W W k C c = < / l a t e x i t > G 2,k < l a t e x i t s h a 1 _ b a s e 6 4 = " n M G 8 J y X Y J W V T G v i / z t r C l u I a M j 0 = " > A A A B 8 n i c d V D L S s N A F J 3 4 r P V V d e l m s B F c S J i k t b W 7 g g t d V r A P S E O Z T C f t 0 M m D m Y l Q Q j / D j Q t F 3 P o 1 7 v w b J 2 0 F F T 1 w 4 X D O v d x 7 j 5 9 w J h V C H 8 b K 6 t r 6 x m Z h q 7 i 9 s 7 u 3 X z o 4 7 M g 4 F Y S 2 S c x j 0 f O x p J x F t K 2 Y 4 r S X C I p D n 9 O u P 7 n K / e 4 9 F Z L F 0 Z 2 a J t Q L 8 S h i A S N Y a c k 1 z e t B 5 p x P Z q Y 5 K J W R V a v a D a c O k Y V Q r V p r 5 M S u o I o N b a 3 k K I M l W o P S e 3 8 Y k z S k k S I c S + n a K F F e h o V i h N N Z s Z 9 K m m A y w S P q a h r h k E o v m 5 8 8 g 6 d a G c I g F r o i B e f q 9 4k M h 1 J O Q 1 9 3 h l i N 5 W 8 v F / / y 3 F Q F l 1 7 G o i R V N C K L R U H K o Y p h / j 8 c M k G J 4 l N N M B F M 3 w r J G A t M l E 6 p q E P 4 + h T + T z q O Z V e t i 1 u n 3 E T L O A r g G J y A M 2 C D O m i C G 9 A C b U B A D B 7 A E 3 g 2 l P F o v B i v i 9 Y V Y z l z B H 7 A e P s E h Z m Q E Q = = < / l a t e x i t > G 1,k < l a t e x i t s h a 1 _ b a s e 6 4 = " + e t N 4 W H M v 6 d J d j 9 H W m S + H t l Q w q M = " > A A A B 8 n i c d V D L S g M x F M 3 4 r P V V d e k m 2 B F c y J A Z 2 9 r u C i 5 0 W c E + Y D q U T J q 2 o Z n M k G S E M v Q z 3 L h Q x K 1 f 4 8 6 / M X 0 I K n r g w u G c e 7 n 3 n j D h T G m E P q y V 1 b X 1 j c 3 c V n 5 7 Z 3 d v v 3 B w 2 F J x K g l t k p j H s h N i R T k T t K m Z 5r S T S I q j k N N 2 O L 6 a + e 1 7 K h W L x Z 2 e J D S I 8 F C w A S N Y G 8 m 3 7 e t e 5 p 6 P p 7 b d K x S R g 9 x a B Z U h c r w L V K 1 W D C l 5 5 V I N Q d d B c x T B E o 1 e 4 b 3 b j 0 k a U a E J x 0 r 5 L k p 0 k G G p G e F 0 m u + m i i a Y j P G Q + o Y K H F E V Z P O T p / D U K H 0 4 i K U p o e F c / T 6 R 4 U i p S R S a z g j r k f r t z c S / P D / V g 2 q Q M Z G k m g q y W D R I O d Q x n P 0 P + 0 x S o v n E E E w k M 7 d C M s I S E 2 1 S y p s Q v j 6 F / 5 O W 5 7 g l p 3 z r F e t o G U c O H I M T c A Z c c A n q 4 A Y 0 Q B M Q E I M H 8 A S e L W 0 9 W i / W 6 6 J 1 x V r O H I E f s N 4 + A Z G 4 k B o = < / l a t e x i t > f M µ,k < l a t e x i t s h a 1 _ b a s e 6 4 = " r Tu M N X v E A v h m B L F o z O V s 0 k J x c / U = " > A A A B + H i c b V D L S s N A F J 3 4 r P X R q E s 3 g 4 3 g Q k p S F F 0 W 3 L g R K t g H t L F M p p N 2 6 M w k z E O o o V / i x o U i b v 0 U d / 6 N 0 z Y L b T 1 w 4 X D O v d x 7 T 5 Q y q r T v f z s r q 2 v r G 5 u F r e L 2 z u 5 e y d 0 / a K r E S E w a O G G J b E d I E U Y F a W i q G W m n k i A e M d K K R t d T v / V I p K K J u N f j l I Q c D Q S NK U b a S j 2 3 5 H l x L + t y c z a a P N x 6 X s 8 t + x V / B r h M g p y U Q Y 5 6 z / 3 q 9 h N s O B E a M 6 R U J / B T H W Z I a o o Z m R S 7 R p E U 4 R E a k I 6 l A n G i w m x 2 + A S e W K U P 4 0 T a E h r O 1 N 8 T G e J K j X l k O z n S Q 7 X o T c X / v I 7 R 8 V W Y U Z E a T Q S e L 4 o N g z q B 0 x R g n 0 q C N R t b g r C k 9 l a I h 0 g i r G 1 W R R t C s P j y M m l W K 8 F 5 5 e K u W q 7 5 e R w F c A S O w S k I w C W o g R t Q B w 2 A g Q H P 4 B W 8 O U / O i / P u f M x b V 5 x 8 5 h D 8 g f P 5 A x N N k g I = < / l a t e x i t > k < l a t e x i t s h a 1 _ b a s e 6 4 = " s u 9 z 8 2 b z h I X H o 6 D c G U w 7 B k J D l p g = " > A A A B 7 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L D a C p 5 I U R Y 8 F L x 4 r m L b Q h r L Z b t q l m 9 2 w u x F K 6 G / w 4 k E R r / 4 g b / 4 b t 2 0 O 2 v p g 4 P H e D D P z o p Q z b T z v 2 y l t b G 5 t 7 5 R 3 K 3 v 7 B 4 d H 1 e O T t p a Z I j Q g k k v V j b C m n A k a G G Y 4 7 a a K 4 i T i t B N N 7 u Z + 5 4 k q z a R 4 N N O U h g k e C R Y z g o 2 V A t e d u O 6 g W FIG. S1.Mapping to large spins without gauge transformation.a) Example of the mapping without initial gauge transformation for M = 3 patterns and Nsp = 8 spins.The coupling between the µ-th mode and the k-th spin is encoded in the coefficients G µ,k .In this version of the mapping, we do not make the initial gauge transformation but rather start reordering the first pattern G 1,k .We relabel spins in such a way that all G 1,k = 1 appear first, and are followed by all the G 1,k = −1.This merely amounts to a relabeling of the spins.We then proceed to reorder also the other patterns in an analogous way to what done in the main text.We notice that, in this way, the mapping produces 2 M subsets of spins Γ k , which are equally coupled to each mode, (Γ1 and Γ8 are highlighted in the figure for clarity).b) These subsets of spins couple to bosons as specified by the matrix f M µ,k .

t) 2 =
i H N , m N a,k − m a,k (t) 2 = L N m N a,k m N a,k − m a,k (t) + m N a,k − m a,k (t) L N m N a,k ,which follows from the fact that the generator annihilates the term proportional to the identity (L N [m a,k (t)] = 0) and that L N acts directly on spins only via a Hamiltonian term.We thus writeD t = ω t L N m N a,k m N a,k − m a,k (t) − ω t ṁa,k (t) m N a,k − m a,k (t) + + ω t m N a,k − m a,k (t) L N m N a,k− ω t m N a,k − m a,k (t) ṁa,k (t) .