Signatures of associative memory behavior in a multi-mode spin-boson model

Spin-boson models can describe a variety of physical systems, such as atoms in a cavity or vibrating ion chains. In equilibrium these systems often feature a radical change in their behavior when switching from weak to strong spin-boson interaction. This usually manifests in a transition from a"dark"to a"superradiant"phase. However, understanding the out-of-equilibrium physics of these models is extremely challenging, and even more so for strong spin-boson coupling. Here we show that non-equilibrium strongly interacting spin-boson systems can mimic some fundamental properties of an associative memory - a system which permits the recognition of patterns, such as letters of an alphabet. Patterns are encoded in the couplings between spins and bosons, and we discuss the dynamics of the spins from the perspective of pattern retrieval in associative memory models. We identify two phases, a"paramagnetic"and a"ferromagnetic"one, and a crossover behavior between these regimes. The"ferromagnetic"phase is reminiscent of pattern retrieval. We highlight similarities and differences with the thermal dynamics of a Hopfield associative memory and show that indeed elements of"machine learning behavior"emerge in strongly coupled spin-boson systems.

Introduction-Cavity, circuit [1] and waveguide quantum electrodynamics (QED) [2], as well as trapped ions [3,4] provide controllable platforms for quantum simulation. These systems are often described as collections of two-level components (spins) coupled via one or more bosonic degrees of freedom (such as photons or phonons). At the most fundamental level, their properties are captured by the Dicke model [5,6] which, under equilibrium conditions and for large number of spins N → ∞, displays a quantum phase transition from a dark to a superradiant phase: for weak spin-boson coupling the average occupation number of bosons is sub-extensive (it grows slower than N with the number of spins), whereas it is extensive (∝ N ) above a critical threshold [7,8]. This phase transition has been widely investigated [9][10][11], observed in Bose-Einstein condensates in optical cavities [12,13] and recently studied in more general, non-equilibrium settings [14,15].
While originally predicted for a single bosonic degree of freedom, this phase transition was eventually generalized to a multi-mode scenario [16], where each spin interacts with several different bosonic modes. Interesting phenomena can emerge in this setting, in particular when the spin-boson couplings are non-uniform, and depend on both the specific spin and the specific boson [17,18]. In this case the superradiant phase features glassy behavior [19][20][21] as well as an energy landscape that shares similarities with simple neural networks [22][23][24][25].
In this work, we consider a driven-dissipative multimode spin-boson model, sketched in Fig. 1(a), and show that its stationary properties resemble those of the socalled Hopfield neural network (HNN) [26,27], a classical spin model behaving like a basic associative memory. More precisely, we identify a crossover between a disor-FIG. 1. Disordered, dissipative spin-boson system. (a) N spin-1 2 particles are strongly coupled to two bosonic modes with couplings giµ (i = 1, ..., N , µ = 1, 2). The spins are weakly driven at a strength Ω ≪ giµ. Each bosonic mode is subjected to two dissipative processes, gain and loss, occuring at rates κ and γ respectively. Patterns are encoded in the couplings between the spins and the bosons, corresponding to ⃗ Gµ = [sign(g1µ), ..., sign(g N µ )] T (µ = 1, 2). (b) Phase diagram of the overlap mµ = (Gµ ⋅ ⃗ σstat) N between the classical, stationary spin configuration and the pattern as a function of the parameter η = (γ − κ) ω (for details see text). A smooth transition between two phases occurs, these being characterized by the presence of several (ferromagnetic or "retrieval" phase) and one (paramagnetic one) basins of attraction. The insets show sketches of the effective free energy landscape. dered "paramagnetic" phase and an ordered "ferromagnetic" one [see Fig. 1 ; the latter appears to be closely arXiv:2003.01004v1 [quant-ph] 2 Mar 2020 related to the retrieval phase of the HNN, characterised in turn by the ability to recall previously stored information. We furthermore explore the impact of fluctuating spin-boson coupling constants on this crossover. Our results show that analogies between multi-mode spin-boson systems and phases of the HNN are not limited to equilibrium settings [21].
The Hopfield model: a brief recap-In physical terms, the HNN is a fully-connected classical spin model consisting of N Ising variables ("neurons") σ i = ±1, i = 1 . . . N , collated to form configurations ⃗ σ = (σ 1 , . . . , σ N ) T . Information is stored in the form of a number M of special configurations, or "patterns", ⃗ ξ µ , µ = 1 . . . M . In what follows, we employ Roman (Greek) letters for site (pattern) indices. The patterns enter the definition of the energy function via the connectivity matrix J ij = 1 N ∑ M µ=1 ξ iµ ξ jµ . Under the hypotheses that (i) the patterns are "few" M N < 0.14 and (ii) different patterns are approximately orthogonal, or more precisely, that they are chosen in such a way that lim N →∞ ⃗ ξ µ ⋅ ⃗ ξ ν N = δ µν , one can prove that the energy is minimised by the 2M configurations ⃗ σ = ± ⃗ ξ µ (µ = 1 . . . M ) [28]. In other words, the patterns (and their opposites − ⃗ ξ µ ) are the ground states of the system and an ideal annealing procedure to zero temperature would recover them perfectly. In order to quantify such retrieval, it is convenient to introduce the M overlaps ζ µ = ⃗ ξ µ ⋅ ⃗ σ N which measure how much the configuration aligns with any given stored pattern. For sufficiently large N , a frequent choice is to describe the pattern components ξ iµ as independent random variables which take the values ±1 with equal probabilities.
Similarly to a fully-connected Ising model, the HNN undergoes a continuous equilibrium phase transition at inverse temperature β = 1, from a paramagnetic phase (β < 1) in which ζ µ → 0 ∀µ to a retrieval (or ferromagnetic) phase β > 1 where a single component ζμ acquires a non-vanishing value, while ζ µ → 0 ∀µ ≠μ. In the retrieval phase (and in the thermodynamic limit), a thermal dynamics (e.g. heat-bath) implemented on the HNN will eventually lead to the configuration approximately reproducing the features of one of the stored patterns. The specific choice can be dictated by the initial condition: an initial configuration with small, but non-negligible, overlap ζμ will restrict the dynamics to the corresponding basin of attraction and the configuration at long times will fluctuate approximately around ± ⃗ ξμ. This ability to faithfully, up to some noise, reconstruct a pattern from partial initial information is what qualifies the HNN as an associative memory.
Non-equilibrium spin-boson model-We now introduce the open quantum spin-boson system of interest. As sketched in Fig. 1(a), we consider N spin-1 2 particles interacting with M independent bosonic modes, according to the Dicke-like [5][6][7][8] x i . Here,σ x,y,z i are the i-th spin Pauli operators,â µ andâ † µ the annihilation and creation operators of the µ-th bosonic mode, ω µ the corresponding frequency. The parameter Ω drives transitions between spin states and the g iµ s are the spin-boson couplings. For concreteness, we assume these to be independent, identically-distributed, real random variables. Additionally, the system exchanges bosons with a Markovian bath. The system state ρ evolves according to a Lindblad equation [29,30] µ describe independent processes of loss and gain of bosons at rates γ µ > κ µ ⩾ 0.
As for the case of a single boson [31], we can obtain an effective master equation restricted to the spin degrees of freedom only. The technical details are reported in the Supplementary Material [32]; here we summarize the main conceptual steps: for vanishing Ω, eachσ z i is a conserved quantity; it is therefore natural to focus our attention on the z-component eigenbasis ⃗ which is defined in terms of "classical configurations" ⃗ σ. Due to their conservation, each configuration labels, at Ω = 0, a subspace of states disconnected from the others, where the only non-trivial evolution takes place in the bosonic part. In particular, for any fixed ⃗ σ⟩ the bosonic stationary state is a displaced Gaussian state ρ ⃗ σ . This implies that the dynamics has a degenerate stationary space spanned by the 2 N elementary combinations ⃗ σ⟩ ⟨⃗ σ ⊗ ρ ⃗ σ . This "classical" subspace gets dynamically coupled to the remainder of the space by the introduction of the term ∝ Ω; however, as long as Ω is sufficiently small, standard perturbative techniques [33,34] can be employed to project the resulting dynamics back onto the classical subspace. By additionally tracing over the bosonic modes, a reduced spin dynamics is found [35,36]. Up to order Ω 2 , this is encoded in a classical rate equatioṅ with p ⃗ σ the probability to find the spins in the configuration ⃗ σ and W ⃗ σ→⃗ σ ′ the transition rate for switching from configuration ⃗ σ to ⃗ σ ′ . Due to the structure of the perturbative term, the only allowed elementary processes are single spin flips, i.e. W ⃗ σ→⃗ σ ′ ≠ 0 exclusively when ⃗ σ and ⃗ σ ′ differ by a single spin. Restricting for simplicity to γ µ ≡ γ, κ µ ≡ κ, ω µ ≡ ω, ∀µ, the derivation of Eq. (2) can be found in Ref. [32]. We also set ω = 1.
Energy function, overlaps and coupling distribution-Remarkably, transition rates only depend on the spin configuration through the quantity ∆E i = σ i ∑ µ,j≠i g iµ g jµ σ j , which can be regarded as the cost of flipping the i-th spin given an energy function of the form Equation (3) bears a clear resemblance to the Hopfield energy (1) if one replaces the patterns ξ iµ with the spin-boson couplings g iµ . With this analogy in mind, we interpret the spin-boson couplings as "noisy patterns". Specifically, we pick the coupling costants g iµ from a bimodal distribution peaked around ±1, given by the symmetric superposition of two Gaussians N g0,s (g) = (2πs 2 ) 1 2 exp −(g − g 0 ) 2 (2s 2 ). From this perspective, it is natural to introduce, as prospective order parameters, and in analogy to the HNN, the overlaps For the sake of simplicity, we focus on the behaviour of m µ upon varying only the parameter η = γ − κ that quantifies the net bosonic loss rate. As we show further below, this parameter controls, to some extent, effective thermal fluctuations. We further fix κ γ = 0.9, N = 50, M = 2. Finally, as transition rates are ∝ Ω 2 , we rescale the time accordingly and effectively set Ω = 1.
Non-equilibrium dynamics-We simulate the dynamics (2) via kinetic Monte Carlo methods [37][38][39]. For convenience, we call different realizations of the stochastic process "trajectories", reserving "realizations" for different random choices of the couplings.
As a starting point of the analysis we consider that the HNN dynamics within the retrieval [paramagnetic] phase is characterized by the presence of multiple [a single] basins of attractions ζ µ ≈ ±1 [ζ µ ≈ 0]. Two trajectories (at fixed coupling realization) of the stochastic process (2) are shown in Fig. 2(a) for two distinct η regimes. In both cases, the initial configuration is 80% aligned with the first pattern, implying m 1 (t = 0) = 0.6. For ease of visualisation we apply a gauge transformation σ i → G i,1 σ i , G i,µ → G i,1 G i,µ which aligns all the components of the first pattern. The components are then re-ordered to bring all the positive components of the second pattern on one side and all the negative ones on the other. In the small-η regime (left-hand side) the dynamics is trapped for long times in configurations close to either pattern; these long-lasting periods are separated by fast switching from one pattern to the other, a behavior commonly seen in finite-size versions of systems undergoing spontaneous symmetry breaking. Phenomenologically [see Fig. 1(b)], the free energy landscape over the space of configurations breaks into different, approximate basins of attraction (BA); the dynamics tends to remain confined in one such potential well until a rare, sufficiently large fluctuations is able to overcome the barrier, after which the dynamics is trapped again for long times in the new BA.
A similar behavior would also be observed, at finite size, in the thermal dynamics of the HNN within the retrieval phase. For large η (right-hand side), instead, memory of the initial state is quickly lost; the dynamics does not appreciably approach either pattern, implying m 1 ≈ m 2 ≈ 0 and corresponding to a case with a single, trivial BA, analogous to the paramagnetic phase of the HNN.
The emergence of multiple BAs for small η, in which the configuration tends to fluctuate close to one of the recorded patterns, shows that, in spite of the dissipation, the system is still dynamically capable of retrieving part of the stored information. In order to quantify this very capacity, we consider averages over a number N traj of trajectories, which we denote by O(t) for an observable O at time t. Figure 2(b) shows the typical evolution of the average overlaps m µ (t) (grey and black curves) at fixed coupling realization, η = 1 and N traj = 200, where the absolute value is taken to avoid averaging to zero due to the ⃗ σ ↔ −⃗ σ symmetry of the model. These quantities clearly assume finite values at long times. These values, however, are not indicative of the overlaps that can be actually achieved in any given trajectory. This is because, once the dynamics starts exploring different BAs, at most one overlap will be appreciably different from zero at any given time. We clarify this with an example: say that, at a large time t, pN traj trajectories (with p < 1) align with the first pattern giving m 1 (t) = m, whereas the others align with the second one (m 1 (t) ≈ 0). Then, m 1 (t) ≈ pm < m. To obviate this reduction, we also study the averaged maximum M(t) = max( m 1 (t) , m 2 (t) ) which constitutes a one-component order parameter symmetric under sign change ⃗ σ ↔ −⃗ σ and pattern permutation (1 ← 2); this symmetrization makes M(t) a more reliable estimate of the overlap that can be achieved within individual trajectories on either pattern.
Stationary properties-To estimate the extent of the retrieval regime, we now focus on the stationary (longtime) properties and study M(t → ∞) as a function of η. We additionally account for the fact that this quantity depends on the specific realization of the couplings and is thus a random variable. We thereby denote by ⟨M⟩ g its average over the previously-defined doubly-peaked distribution of the g iµ s and plot it against η in Fig. 2(c).
As discussed in Ref. [31], when considering the large and the small η limits, effectively thermal regimes can occur. At large η, an infinite temperature scenario emerges. The finite value taken by the disorder-averaged ⟨M⟩ g in this regime is due to the finite size of the system. As η is decreased, a crossover emerges towards larger and larger values and for η ≈ O(1) it seems to reach a maximum of approximately 0.5, meaning 75% of the configuration spins are aligned with a pattern. In this regime, the error bars get larger as well, implying that the system becomes more sensitive to the specific values taken by the spin-boson couplings. Combining these results with the dynamical ones from the previous section, showing the emergence of several approximate BAs, suggests that the two regimes, η large and small, reproduce to an extent the physics of the paramagnetic and retrieval phases of a HNN, with η playing the role of an effective temperature. Future investigations at variable system size will shed further light on the features of the crossover.
Influence of pattern noise-We recall that, in order to connect the physics of the HNN to our model, we considered the couplings ⃗ g µ as noisy patterns ⃗ G µ . One can thus distinguish between two sources of fluctuations: (i) the randomness of the patterns, whose components can be either ±1 with equal probability , and (ii) the noise due to the finite width s of the distribution peaks around ±1.
It is worth noticing that, while (i) is present in the HNN as well, (ii) represents an additional source of noise specific to our model. For this reason, we wish to better understand its impact on the retrieval ability of the dynamics, reporting the result in Fig. 2(d). After having excluded the rightmost points η > 6, we compare the data used for the curve in panel (c), corresponding to s = 0.25 (red) to two additional data sets analogously obtained for s = 0.125 (magenta) and 0.063 (blue). Interestingly, reducing s does not seem to affect the standard deviation (i.e. the error bars); instead, one can spot a decrease in the average values ⟨M⟩ g . This suggests that small amounts of noise over the patterns can very slightly improve the retrieval of stored information. Although the cause is not entirely clear, it is worth mentioning that a similar effect can be observed in the HNN too (see Ref. [32] for further details), where, at finite size, the introduction of disorder over the patterns induces a broadening of the order parameter profile. This broadening results in an effective increase of the overlaps close to the critical point; as one moves deeper in the retrieval phase, however, the effect is reversed and the average overlap decreases at larger noise levels.
Conclusions-We have investigated the dynamics and the stationary phases of a dissipative multi-mode spinboson model, highlighting analogies between this many-body system and the HNN, such as a pattern retrieval dynamics and a crossover between a retrieval and a paramagnetic phase. Beyond mimicking the HNN, physical realizations of this spin-boson Hamiltonian, e.g. atomcavity setups, have the potential to systematically probe the quantum regime, i.e. go beyond the perturbative limit considered here. Increasing the degree of "quantumness", i.e. adding the possibility to host superpositions and entanglement offers intriguing potential for exploring quantum effects in the context of machine learning.
Acknowledgments Supplemental Material for "Signature of associative memory behavior in multi-mode spin-boson models"

DERIVATION OF THE RATES
In this section we provide details on how Eq.(2) of the main text is derived. The state of the system, ρ(t), evolves according to the Lindblad equation In order to obtain an effective description of the spin dynamics at strong coupling, we consider the evolution of the state aṡ ρ = L 0 (ρ) + L 1 (ρ), with L 0 = L − L 1 and We firstly focus on the dynamics generated by L 0 . It is worth noticing that each spin operatorσ z i is a conserved quantity under the action of the superoperator L 0 . As a result, the eigenstates ofσ z i operators, ⃗ σ⟩ = (σ 1 , ..., σ N ) T ⟩, withσ z i ⃗ σ⟩ = σ i ⃗ σ⟩ and th σ i = ±1, identify 2 N spin sectors. These are spanned by the combination ⃗ σ⟩ ⟨⃗ σ and represent 2 N degenerate stationary subspaces. Thus, we can assume the stationary state of L 0 of the form σ are a set of classical probabilities and ρ σ is the corresponding bosonic state, that we assume to be a Gaussian state.
Secondly, we consider the superoperator L 1 , which mixes different spin sectors, as a perturbation with respect to L 0 , i.e. we consider Ω ≪ g. Indeed, projecting the dynamics onto the stationary manifold of L 0 , we exploit the Nakajima-Zwanzig formalism, up to the order O(Ω 2 ), to write the evolution of the spin as where Pρ spin (τ ) = Tr B [ρ stat (τ )] = ∑ ⃗ σṗ⃗ σ (τ ) ⃗ σ⟩ ⟨⃗ σ , and Tr B is the partial trace over the bosonic modes. Furthermore, we have defined the spin-configuration dependent superoperators with M lµ = ∑ n≠l g nµ σ n , and D γµ , D κµ the dissipative terms representing loss and gain, respectively. By projecting Eq.(S1) on a state ⃗ σ ′ ⟩, the dynamics reduces to the evolution of the classical probabilities ruled by a master equation asṗ , is the transition rate for the switching ⃗ σ → ⃗ σ ′ and it allows only spin flip processes.