Optical computation of a spin glass dynamics with tunable complexity
OOptical computation of a spin glass dynamics with tunable complexity
M. Leonetti,
1, 2
E. H¨ormann, L. Leuzzi,
2, 3
G. Parisi,
3, 2, 4 and G. Ruocco
1, 3 Center for Life Nano Science@Sapienza, Istituto Italiano di Tecnologia, Viale Regina Elena, 291 00161 Rome, Italia. CNR NANOTEC-Institute of Nanotechnology, Soft and Living Matter Lab, P.le Aldo Moro 5, I-00185 Rome, Italy Department of Physics, University Sapienza, P.le Aldo Moro 5, I-00185 Roma, Italy INFN, Sezione di Roma I, P.le A. Moro 2, 00185 Roma, Italy (Dated: July 2, 2020)Spin Glasses (SG) are paradigmatic models for physical, computer science, biological and socialsystems. The problem of studying the dynamics for SG models is NP-hard, i.e., no algorithm solvesit in polynomial time. Here we implement the optical simulation of a SG, exploiting the N segmentsof a wavefront shaping device to play the role of the spin variables, combining the interference atdownstream of a scattering material to implement the random couplings between the spins (the J ij matrix) and measuring the light intensity on a number P of targets to retrieve the energyof the system. By implementing a plain Metropolis algorithm, we are able to simulate the spinmodel dynamics, while the degree of complexity of the potential energy landscape and the region ofphase diagram explored is user-defined acting on the ratio the P/N = α . We study experimentally,numerically and analytically this peculiar system displaying a paramagnetic, a ferromagnetic anda SG phase, and we demonstrate that the transition temperature T g to the glassy phase from theparamagnetic phase grows with α . With respect to standard “in-silico” approach, in the opticalSG interaction terms are realized simultaneously when the independent light rays interferes at thetarget screen, enabling inherently parallel measurements of the energy, rather than computationsscaling with N as in purely in − silico simulations. The solution of large combinatorial problems demands for novel hardware architectures enabling for faster andinherently parallel calculation. An emerging trend is that of pairing an optical layer into a specific digital or analogcomputation scheme, in order to improve performance while reducing computational costs and processing times.Optical computing promises parallel processing and high bandwidth which may be eventually performed in free space,with limited power consumption (e.g., Fourier transform performed by a lens). Optical computation is an emergingscheme in quantum transport [1], quantum simulation [2], and machine learning [3], and can be implemented ondifferent platforms including free space, photonic chips [4] and optical fibers [5]. One of the advantages brought byoptics is that certain operations can be performed at the “speed of light”. Indeed, the evaluation of a matrix productcan be estimated in the time needed for a properly shaped light beam to pass through a diffractive pattern opportunelytailored to mimic the requested transfer matrix [6]. By exploiting last generation optical modulation devices, millionsof light rays can be driven simultaneously between several states and within a microsecond time frame, thus potentiallyproviding a scalable optical platform that only need to be properly projected on to the relevant and computationallyhard problem.Spin glasses [7] serve as prototype models, capable to provide nontrivial equilibrium and off-equilibrium phenomenol-ogy [8, 9]. In particular, the dynamics in an energy landscape with many equilibrium states and the origin of (multiple)relaxation times in finite dimensional systems, are open questions in modern statistical mechanics [10–13]. Complexsystems from diverse fields fall into the spin-glass universality class, like, e.g., brain functions [14], random lasers[15, 16], and quantum chromo-dynamics [17]. Indeed, novel methods for the calculation of the equilibrium states andof the dynamics of a spin glass system are highly desiderate.Here we propose an optical system able to compute the energy of a given spin-glass state. We integrated such opticallayer onto a standard digital computation layer to realize an optical spin-glass (OSG) dynamics simulation. Our ideastems from the observation that the overall intensity I = (cid:80) Pν =1 I ( ν ) at P given points ( ν ) on a screen placed at thedownstream of a strongly scattering medium shone with N coherent light rays from a single laser, can be formallywritten as a spin glass Hamiltonian. Thus scattering, coupled with an adaptive optical element, has been proposed asan instrument to access the spin-glass dynamics and employed for low complexity (small P values) simulations [18].Successively a similar approach has been employed to find ground states of large transmission matrices [19].As a starting point for our experiment it is easy to observe that, in the simplified case in which the i th light rayfield has a complex amplitude a i = A i e ıφ i , the single target contribution I ( ν ) reads a r X i v : . [ c ond - m a t . d i s - nn ] J u l FIG. 1: A sketch of the experimental setup. Laser light (Azure Light 532 nm), reflected by the de DMD, is then scattered byan opaque medium. The far facet of the scattering medium, is then imaged on a camera, after passing on a linear polarizer.The inset shows the measured P ( I off diag ) and P ( I diag ). They have bene fitted with the function P ( I ∗ ) = A ∗ exp( − I ∗ /B ∗ ).We retrieve B diag = 0 .
08 and B off diag = 0 .
04. The values of the B ∗ are consistent with the predicted behavior of I diag = v = | ξ | and I off diag = v + v = 2 v R = 2 (cid:60) [ ξ ¯ ξ ] (see methods), that is : P ( v Rij ) = 1 / (2 σ ) exp −| v Rij | /σ and P ( v ii ) =1 / (2 σ ) exp − v ii / (2 σ ) . I ( ν ) = E ( ν ) E ( ν ) † = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N (cid:88) i =1 ξ νi a i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = N (cid:88) i,j ξ ( ν ) i ¯ ξ ( ν ) j a i ¯ a j (1)= N (cid:88) i,j | ξ ( ν ) i || ξ ( ν ) j | A i A j e ı (arg ξ ( ν ) i − arg ξ ( ν ) j + φ i − φ j ) where ξ ( ν ) i = | ξ νi | e ı arg ξ ( ν ) i are the complex transmission matrix elements from the ith incoming beam to the target ν .Input illumination is controlled by a DMD with the superpixel method (see methods) . This approach enables toseparate the input laser wavefront in many segments (up to 50 ) each one composed by 4 DMD mirrors. Each segmentcan be then programmed into one between two states each characterized by a phase factor φ i = ± π , equivalent to anamplitude factor of S i ∈ {− , } , so that the single target intensity (1) can be rewritten as I ( ν ) = ,N (cid:88) ij v νi,j S i S j , (2)where all transmission matrix elements and input field amplitudes have been included in the coefficient v ( ν ) ij ≡ A i A j ξ ( ν ) i ¯ ξ ( ν ) j . (3)Let us stress that though v ( ν ) ij are complex-valued, the intensity I ( ν ) is always a real number because v ( ν ) ij = ¯ v ( ν ) ji andthe sum in Eq. (2) runs on all i, j = 1 , . . . , N . Amplitudes A i are defined by the laser intensity. By using a laserwith a Gaussian beam and expanding it to retrieve an homogeneous distribution of the intensity on the active DMDarea, it is possible to approximate all the A i , for any i , to a constant over the DMD segments. Maximizing the overallintensity with respect to DMD spin S configurations, finally corresponds to minimize the following Hamitonian H [ S ] = − ,N (cid:88) ij J ij S i S j , (4)where we have introduced the interaction matrix J ij ≡ N P (cid:88) ν =1 v ( ν ) ij , (5)properly rescaled with N in order to provide thermodynamic convergence in the large N limit also when the numberof targets grows like the number of spins: P = αN , with α = O (1). Note that the matrix J is a Hermitian matrix, J ij = J † ji , and H [ S ] is, therefore, real.Equation 4 is, by all means, a spin glass Hamiltonian, more precisely it is a generalization to complex continuousvalued patterns ξ of the Hopfield model [20, 21], for which the study of Amit, Gutfreund and Sompolinsky [22–24]predicts the existence of a low temperature/large P spin-glass phase. That is, a phase with multiple degeneratefrustrated thermodynamic states, and critical slowing down dynamics approaching the glassy transition from theparamagnetic phase. To check this, we realized experimentally our optical simulation of a spin glass. As shown inFig. 1 we exploited laser light form a stabilized Nd-YAG laser (Azure light 0.5 W), to shine a Digital MicromirrorDevice (DMD, Vialux V-7000) controlled by a CPU based workstation.The same computer simultaneously monitors the Camera (Ximea CB013MG-LX-X8G3 64 Gbit/s speed on thePCIe Gen3 lane) which collects light at the downhill of a scattering medium (Thorlabs ground glass diffuser, ThorlabsDG05-220). A two lens system, enables to monitor in the camera plane the intensity corresponding to the output facetof the diffuser with a factor 8 magnification (Focal of Lens1 is 25 , I ( ν ) changes on the targets (asingle target is a set of neighboring camera pixels) ν . If I ( ν ) increases, then the change is accepted, otherwise thechange is accepted only with a probability p = e ∆ I/T where ∆ I is the measured variation in intensity following thespin/micromirror flip and T is an user defined system temperature. The spin flipping is performed acting on therelative segment on the DMD (flipping time is 18 µs ), while the spin coupling is realized optically thanks to thenearly instantaneous light propagation into the disordered medium, the intensity reading is performed through thecamera, while the move acceptance is performed digitally by the computer.By storing spin configurations S ( t ) at each time step t we are able to extract the temporal behavior of the connectedautocorrelation function: F self ( τ ) = 1 N N (cid:88) i (cid:104) S i ( t ) S i ( t + τ ) (cid:105) c ; (6)where (cid:104) . . . (cid:105) indicates the averaging over Monte Carlo steps t . As an instance, we show the case of an optical spinmodel with N = 225 spins whose dynamics we simulated for t max = 2 × Monte Carlo steps (each step consists of N micromirror flips, with micromirrors selected uniformly and in random fashion) from which we extract the behaviorof the correlation F self ( τ ) on a more limited temporal window for τ in order to perform a proper temporal averaging.The results for a single target (number of targets P = 1) is reported in Fig. 2a, displaying F self ( τ ). In this case,the correlation function decays rapidly to zero at high temperature and lowering the temperature its behavior crossesover to a power-law decay. This slowing down of the dynamics is called critical relaxation[26–30]. Eventually, at lowenough T , F self ( τ ) tends towards a non-zero plateau. To check further the nature of the energy landscape, we run 80( p = 1 , . . . ,
80 ) very short simulations ( τ max (cid:39)
30) with fixed target ν at zero T . In Fig. 2b, we report the final states S final (the DMD square matrix has been reshaped to a 1D array and visualized as black or white squares respectivelyrepresenting S = − S = 1) for each one of the simulations. We notice that, apart for a small fluctuation dueto measurement noise on a limited set of i ’s, only two configurations S final are appearing. One, ↑ S , with an overallpositive magnetization 1 /N (cid:80) i S i , and one, ↓ S , with negative magnetization, as graphically reported reported inFig. 2e ). They are strongly anti correlated as shown by the q parameter reported in Fig. 2c and 2d, for all the S final ( p ) obtained from 80 replicas of the dynamics at T = 0. The parameter q ( p , p ) shown there, represents thedegree of similarity between the S final of two replicas of the system and is calculated as the normalized scalar product: q ( p , p ) = S final ( p ) · S final ( p ) /N . The q ( p , p ) matrix, appearing symmetric and with ones on the diagonal, hasbeen visualized with p and p sorted in order to clusterize similar S final . In this favorable configuration, states areaappearing as yellow squares on the q matrix. The two squares appearing in Fig. 2c-d are states composed two thespin-reversed states confirming that in the case P = 1 we have a transition to a ferromagnetic phase.Indeed, the occurrence of only a pair of opposite states may be explained if we go back to the formal description ofthe J ij in Eqs. 3, 4, 5, with P = ν = 1 target: H [ S ] = − N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N (cid:88) i =1 ξ (1) i S i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = − | ξ (1) · S | N (7)in which we exploited the approximation of constant amplitudes for all the incoming light rays and we rescaled A i ξ ( ν ) i (cid:39) Aξ ( ν ) i → ξ ( ν ) i . The Hamiltonian is just proportional to the scalar product of the spin configuration arrayand the array of the transmission matrix elements ξ (1) i from the micromirror-created spins ( i ) to the target (1). Thisis maximum for two configurations: the one maximizing the field along the transmission vector ( ↑ , see Fig. 1) andthe other maximizing the field along the opposite direction. It is in fact well known that, given a certain transmission FIG. 2: Experimental results for the 1-target Regime. a) F self as a function of τ , that is, the number of Metropolis sweeps.The inset shows the first steps in logarithmic τ scale. b) Final states obtained after a series of repeated zero temperaturesimulations with the same disorder realization. The degree of complexity is low (just two opposite states are present). In c)and d) we report the overlap q ( p , p ) – defined in equation (8) – between the configurations S ( p ) and S ( p ) of two replicatedsystems and sorted by a k -means algorithm in order to be organized into two clusters. Replicas are labeled by the indexes p and p running on all the simulated systems. The appearance of only two different yellow squares indicates that the configurationsof the replicas thermalize into states that can be grouped in two clusters of very similar states one opposite to the other. d) isthe same of c) for a different realization of disorder. In e) we report the two (opposite) states found in c). matrix, it is possible to find the best input configuration producing maximum intensity [31] and that this configurationis “unique”. FIG. 3: Results for P = 36, N = 225. In a) the F self as a function of τ indicating the sweeps. In b) we report the degree ofsimilarity obtained for clusterized states resulting many optimizations at zero temperature. The phase space results clusterizedin many states. In c) d) and e) we report respectively A , τ and a as a function fo 1/T, comparing the results form the digitalsimulation (red full dots) with that of the optical simulation (open blue symbols) (see methods) . In panel d) τ is reportedonly in the region where the exponential contribution results a reasonable interpolation while in panel e) a is reported only inthe critical region. The interaction matrix J ij = ξ (1) i ¯ ξ (1) j /N is a diadic matrix constructed on a single vector ξ (1) with the Hebb rule[32]. Resulting from a diadic matrix, our scattering system is, thus, behaving as a sort of trained optical memory inwhich the pattern ξ (1) is the memory learned by the network. In this neural network contest “learning” means varyingthe couplings. The network is designed in such a way that the learnt patterns ξ are retrieved as stable configurations,that is, in this example, the spin configurations ↑ S and ↓ S . In optical words, instead, when P = 1 the experiment isidentical to a standard wavefront shaping experiment [33] in which there exist a single spin configuration maximizingintensity [34]. The case where the ξ are Boolean, rather than complex and Gaussian distributed, is also called theMattis model [35].If P is larger and, in particular, if it is so large to scale with the number of variables N , the situation is similar tothe Hopfield model of Amit, Gutfreund and Sompolinksy [22]. In that seminal paper the authors demonstrate that itis possible to store states into a neural network memory by exploiting a sum of diadic matrices generated from thatmany different vectors. However, there is a limit to the number of independent states that can be stored in such a kindof memory[36] When this limit is surpassed, the multiple states start to combine, thus generating additional states tothose related to the transmission patterns. All these states are minima of a complex energy landscape displaying aspin-glass phase in a given region of the phase diagram. To experimentally engineer the complexity of the spin-glassphase in our optical system, then, the only thing that we need to do is to increment the number of targets. We,therefore, perform optical Montecarlo simulation with the Hamiltonian (4) and increasing P .We perform a OSG simulations with N = 225 and P = 4 , , , , ,
225 targets, corresponding to α (cid:39) . .
04, 0 . .
11, 0 .
16, 1 (where α is the ratio between the number of targets P and the number of degrees of freedom N ). We report the results for P = 36 in Fig. 3. In panel 3a we show the autocorrelation functions of the σ ’s as afunction of τ . It is readily noticed that the F self moves from an exponential to a non-exponential relaxation in timeas temperature decreases, that is a typical aspect of many complex systems [27, 28, 37]. In particular the F self ofthe OSG with many targets can be typically described by the generic interpolating function approaching the criticalpoint from high temperature: FIG. 4: F Self ( τ ) for various values of N , P and α . From a (cid:48) to a (cid:48)(cid:48)(cid:48)(cid:48) : each panel has a value of α (the same ) four temperatures.Panels from b (cid:48) to b (cid:48)(cid:48)(cid:48)(cid:48) report the same data organized with a single temperature and four values of α per panel. Panels fromc (cid:48) to c (cid:48)(cid:48)(cid:48)(cid:48) instead report α fixed while N (number of spins) is varied. F self ( τ ) = (1 − A ) e − ττ + A Ct a in which the first term represent the exponential relaxation time while the second term the power law decay appearingapproaching the glassy phase, i. e., crossing over for large τ . The A parameter represents the relative weight of thetwo terms: the larger it is, the more the power-law-like is relevant. The behavior of the fitting parameters is reportedin Fig. 3, estimated both by means of the data obtained in standard “digital” numerical simulations and with theoptical Monte Carlo methods introduced in the present work. As expected, as temperature is decreased from theparamagnetic phase, we observe (i) an increase of the correlation time τ , (Fig. 3d), (ii) a very sharp increase of theweight A of the power-law decay contribution in the critical region, (Fig. 3c), and (iii) a lowering of the power-lawdecay exponent a in the critical region, (Fig. 3e).The critical slowing down of the relaxation is a signature of a multi-state (free-)energy landscape. Further evidenceof complexity is also shown in an experiment in which multiple short relaxation dynamics at T = 0 have beenretrieved. In Fig.3b) we report the parameter q ( p , p ) for a thousands different simulations. Here the pattern showsmany different clusters indicating that many equilibrium states are populating the potential energy landscape. Wehave been using k -means to organize data into 36 clusters.In Fig. 4 we compare the F self ( τ ) for various values of N , the ratio α = P/N and the inverse temperature β = 1 /T . FIG. 5: Probability distribution of the Parisi parameter q as a function of β , for N=225. Optical simulations have beenperformed for few values of P . a): P = 4; α =0.02 b): P = 16; α =0.07 c): P = 25; α =0.11. P ( q ) in a given β , α configuration,has been obtained by performing 50 fast (25 sweeps) simulations (see methods). In panel d) we report simulations retrieved byperforming 2 very long (8000 sweeps) simulations (see methods) for P = 225; α =1. In panel e) we report the estimated T G bymeans of optical simulations, together with the curve T G = A (1 + √ α ) predicted from Replica theory for N → ∞ (eq. 38).Optical simulation provides a value of A of 1 . ± . A = 1). The inset of panel e) shows theKurtosis of the P ( q ) as a function of β . Continuous line is to guide the eye. In the first set of panels, from a (cid:48) to a (cid:48)(cid:48)(cid:48)(cid:48) , we report the F self ( τ ) for four values of β at fixed N = 256 ( a different valueof α is shown in for each panel). It is possible to note that at high α , the decorrelation at any given temperature isstrongly hindered, if not avoided at all, at the time-scales considered.Panels from b (cid:48) to b (cid:48)(cid:48)(cid:48)(cid:48) report the same data (fixed N = 225) organized at a single temperature in each panel andfour values of α per panel. In panels from c (cid:48) to c (cid:48)(cid:48)(cid:48)(cid:48) , instead, we show the correlation function at fixed α and varying N in each panel, demonstrating that the behavior of F self ( τ ), for the time window considered, has no strong finitesize effects already at N = 256 spins (c (cid:48)(cid:48)(cid:48) and c (cid:48)(cid:48)(cid:48)(cid:48) ).As we already mentioned, in Eqs. (4-5) the number of targets P plays the same role of the number of memoriesto be stored into a Hebbian matrix in the AGS theory [22, 24] predicting the impossibility of memory retrieval, i. e.,the inability to retrieve the encoded ξ states, for α larger than a given critical value α c . This qualitative change inthe degree of complexity of our system can be demonstrated in Fig. 5, where we report the temperature behavior ofthe probability density function P ( q ) of the overlap q ab ≡ N N (cid:88) k =1 S ( a ) k S ( b ) k (8)between couples of replicas S ( a ) and S ( b ) . At high temperature each equilibrium configuration, will be uncorrelatedfrom the others, because of strong thermal fluctuations. For low T , at α (cid:39)
0, that is
P/N → N → ∞ , themodel dynamics will tend to align, or counter-align, to the ξ patterns, as mentioned in Eq. 7 for the generalizedMattis model. The values of the overlap between such states will, then, converge to only two possible values as N increases, one the inverse of the other, and the distribution will display two distinct symmetric peaks. As the numberof patterns to be satisfied increases with the number of micromirrors, though, frustration arises as it becomes moreand more unfeasible to satisfy any of them with a configurations of S minimizing the Hamiltonian H [ S ], cf. Eq.(4). The change of the P ( q ) from a low temperature two peak distribution into a multi peaks distribution is a clearevidence of the onset of such frustration and the relative complexity in the equilibrium state organization. We stressthat in Fig. 5 the overlap distributions displayed are for a single realization of the random couplings. Moreover, fromthese curves we extracted the values of the T G (see methods), and these are reported in Fig. 5e as a function of α for N = 225 (green dots fitted with green curve). The curves represent the fit with the model from phase transitionequation from replica symmetry breaking calculation ( eq. 38) leaving A as a free parameters. By fitting the datafor T G with the model we retrieve a constant A of 1.44 ± A =1). FIG. 6: Self-correlation function for a system with N = 12100 = 110 ×
110 spins and P = 1369 targets. This has been obtainedwith a total of 0 . N for the optical (red dots) and digital (blue dots) simulations. In order to test our optical procedure we performed numerical experiments with a few spins, from N = 16 to N = 225, however the real advantage of the proposed method comes about when the sample size and the number oftargets is large. In Fig. 6 we report a measurement with a large fully connected spin glass ( N = 12100 micromirror-combined spins organized into a 110 ×
110 2D square lattice) with α = 0 . (cid:39) − s/step for this large optical experiment. This TPS turns out to be only barely dependent on N , as we showin the inset of Fig. 6 in which we report the TPS for the Optical simulation as red circles. The time per step for thedigital simulation is, instead, growing with N as expected (inset of Fig. 6 blue circles). The optical simulation is thusoutperforming the digital one at N ∼ N .In conclusion, we demonstrate the possibility to simulate the dynamics of a system with N spins and N couplingswith random couplings and a low temperature spin-glass phase. Our approach has the advantage of (i) being scalable:the possibility to simulate very large systems without affecting the calculation speed, because the instantaneousinterference enables to access the energy difference without requiring to individually calculate all the coupling terms.(ii) varying the constraints density so as to survey both a ferromagnetic-like ordering, retrieving the transmissionmatrix patterns, and a spin-glass freezing, i.e., loss of memory, by playing with the number of optical targets inwhich the intensity is monitored, (iii) vary the dynamic variables of the system, for Ising spins, needing each onefour micro-mirrors to be constructed, both to multiple discrete states variables and to complex continuous phasors.(iiii)The same approach may be employed with multi state spins, employing higher dimensional superpixel capableto generate many different phase and amplitude values. In this configuration, as far as the number of targets perspin is below the critical α c ( T ) value of memory retrieval, the analysis of the ferromagnetic Mattis-like states allowsto reconstruct the transmission matrix elements of the random opaque medium through the states S visited in theoptical simulation, as memories learned in a neural network. AppendixesRandom couplings and transmission matrix
In Eq. (5) the couplings have been defined as J ij = 1 N P (cid:88) ν =1 v ( ν ) ij ; v ( ν ) ij ≡ ξ ( ν ) i ¯ ξ ( ν ) j ; ξ ( ν ) i = ρ ( ν ) i + ıυ ( ν ) i ∈ C (9)where the transmission matrix properties [38] show that real and imaginary parts of ξ are both Gaussian distributedas ( x = ρ, υ ), P ( x ) = 1 √ πσ e − x σ . (10)The J ij matrix is Hermitian because ξ ( ν ) i ¯ ξ ( ν ) j = ξ ( ν ) j ¯ ξ ( ν ) i . The intensity per mirror I/N is then, expressed as: IN = N (cid:88) ij J ij S i S j = N (cid:88) j>i ( J ij + J ji ) S i S j + N (cid:88) i J ii = N (cid:88) i (cid:54) = j J Rij S i S j + const . (11)Note that even if J ij , i (cid:54) = j , is in general complex valued, each term of the sum J ij + J ji = 2 J Rij is real and symmetricin i ↔ j because it results form ξ i ¯ ξ j + ξ j ¯ ξ i . For a large number of targets P , as in the case P = αN with an α = O (1),the J ij element is the sum of many random numbers (each one the product of two Gaussian random numbers) and,for the central limit theorem, will be Gaussian distributed.For a single target P = 1 to derive the probability distribution of the values of J ij = ξ i ¯ ξ j /N we can, instead,integrate the expression (we discard the index ν ) P ( J ii ) = (cid:90) dρ dυ πσ e − ρ υ σ δ (cid:0) J ii N − ρ − υ (cid:1) = N σ e − JiiN σ , (12)for the diagonal entries J ii = | ξ i | , while for i (cid:54) = j the distribution of J Rij = ρ i ρ j + υ i υ j reads P ( J Rij ) = (cid:90) dρ i dυ i πσ dρ j dυ j πσ e − ρ i + υ i + ρ j + υ j σ δ (cid:0) J Rij N − ρ i ρ j − υ i υ j (cid:1) = N √ πσ (cid:90) dz √ z exp (cid:40) − z − N (cid:0) J Rij (cid:1) σ z (cid:41) = N σ e − | JRij | Nσ , (13)where we have used the relationship I ( γ ) ≡ (cid:90) ∞ dz √ z e − z − γ z = √ πe | γ | . To verify that the distributions are, indeed, exponential, thus indirectly confirming the Gaussian distribution ofthe transmission matrix elements ξ , and to estimate their variance by interpolation, we can experimentally measure0the distribution of the J (cid:48) s , or, equivalently, of the v ij = J ij N , cf. Eq. (5) with P = 1. We resorted to a pair ofmeasurements of the propagation of the signal from two superpixels S , S on a single target. In a first run we turnon a single superpixel and measure the diagonal contributions independently: the intensity at the target k = 1 , I diag ,k ≡ v kk = ξ k ¯ ξ k = ρ k + υ k , (14)In the second run we turn on the couple of superpixels S = S = 1 and measure the overall intensity at the target 1 I [ S = 1 , S = 1] = v + v + v + v = ξ ¯ ξ + ξ ¯ ξ + ξ ¯ ξ + ξ ¯ ξ = ρ + υ + ρ + υ + 2 ρ ρ + 2 υ υ whose extra, off-diagonal contribution, is I off diag ≡ I − I diag , − I diag , = ξ ¯ ξ + ¯ ξ ξ = 2 ρ ρ + 2 υ υ = v + v = 2 v R = 2 N J R . (15)Measuring several values of I and I diag ,k we find that both I off diag and I diag , k have an exponential distribution asshown in the inset of Fig. 1. Looking in semi-log scale and fitting with the interpolating function A e − BJ the slopeof the distribution P ( J Rij ) turns out to be twice the slope of P ( J ii ), as predicted by Eqs. (12,13). Superpixel method
Light reflected by the DMD produces a diffraction pattern composed by a square array of intensity maxima. Thesemaxima correspond to maxima of interference due to in-phase summation of all the contribution from the each lightray reflected by pixels. By placing a lens (L1 with focal lenght f = 200 mm) to collect DMD light we have the zerothorder fringe appears at coordinates [ x f , y f ] = [0 ,
0] and the first order at coordinates [ x f , y f ] = [ a, a ] with a the latticeprimitive vector a = λfd , where λ is the wavelength of light (0.532 nm), and d is the DMD pixel Pitch ( 13.68 µm ).As shown in [Opt Expr 22 1364] it is possible to exploit a DMD device to control both amplitude and phase ofcoherent light beam. This technique is based on the phase delays accumulated by light rays corresponding to differentpixels, if light is collected only from special positions in Fourier space. The first step is to collect light thought aspatial filter (iris I ) which is placed off at coordinates [ x f , y f ] = [ a/ , a/
2] in the Fourier space. In such case eachDMD pixel acquires a relative phase delay wick depends on its position in the real space: in particular, if pixel areidentified by a couple of integers [ x i , y i ], pixels with even T = x i + y i contribute with a +1 pre-fator, while , pixels withodd T contribute with a -1 pre-factor. The we organize the DMD into superpixels composed by 2x2 DMD pixels. Inorder to merge each pixel into a single contribution we reduce the optical resolution in order to “blur them together”this is performed thanks to the iris I which is tuned on a diameter of 2.66 mm. This diameter produces an opticalresolution of 32 µm thus blurring together four pixels.In order to ensure the validity of the proposed model we also realized a system in which the size of the superpixelmatches the size of the disorder grain in the scattering medium. This is obtained employing a ground glass diffuser(“DG10-220” thorlabs with grain size of 63 µm ) and an optical system with magnification 2, producing a superpixelimage of 64 µm side. Overlap q , P ( q ) and calculation of T G In the fast protocol (Fig. 5a-b-c) P ( q ) for a given [ β , α ] configuration, has been obtained by performing 50 fast(25 sweeps) simulations of replicas of the system (each starting from a random initial configuration) and calculatinga value of q as the scalar product between the final states for each pair of replicas. P ( q ) is obtained by extractingthe occurrence probability for each q value. Even if this approach provides states not completely thermalized, thecomparison with longer numerical simulations (Fig. 5d) demonstrates a good agreement about the onset of the spin-glass phase. In the long protocol (Fig. 5d) P ( q ) in a given β and α = 1 configuration, has been obtained by monitoringeach replica for a long time (8000 Monte Carlo sweeps in the Optical Spin Glass simulation). Then P ( q ) has beenobtained as the normalized histogram of the scalar product of two replicas at the same time.A first raw estimate of T G has been extracted by calculating at the Kurtosis K ( β ) of the P ( q ) distribution forvarious values of α as the temperature is decreased: K ( T ) = (cid:104) ( q − (cid:104) q (cid:105) ) (cid:105)(cid:104) ( q − (cid:104) q (cid:105) ) (cid:105) (16)1where the average is taken over the distribution P ( q ). We underline that the P ( q ) = P ξ ( q ) considered here are fora given realization of the disorder (transmission patterns ξ , or, equivalently, spin couplings J ). We estimate T G ( N ),at a given system size, as the temperature at which K ( T ) is sensitively less than 3, i. e., P ( q ) is no more a hightemprerature Gaussian. Numerical simulations
To check the performances of the OSG we also have performed standard numerical simulations. The implementeddynamics is a standard Metropolis Monte Carlo, as in the optical spin-glass simulation. We built the random generatorfor the coupling factors ξ , rather than for J , in order to ensure correspondence with the optical simulations, generating M N -components complex vectors ξ (1) , ξ (2) , . . . , ξ ( M ) , which represent the transmission from the superpixel on to thetargets in the optical setup. Subsequently, the J ij element of the coupling matrix is calculated from those vectors asfollows: J ij = 12 M (cid:88) ν =1 (cid:16) ξ ( ν ) i ¯ ξ ( ν ) j + c.c. (cid:17) = M (cid:88) ν =1 (cid:16) ρ ( ν ) i ρ ( ν ) j + υ ( ν ) i υ ( ν ) j (cid:17) (17)During the dynamics, for the energy difference at each step, we only perform the sum over the couples which includethe spin that have been flipped at that step; in such a way that we only need to sum O ( N ) terms instead of all possible O ( N ) terms. The final formula for the energy difference is therefore, assuming the flipped spin to be the ¯ı-th one:∆ H = σ ¯ı (cid:88) j J ¯ı ,j σ j . (18)The simulation have been performed by a multi core cluster composed by • • Replica theory of the complex Gaussian Hopfield model
Defining the extensive “energy” of a superpixel configuration s as H [ s ] = − I/ /N = O ( N ), the statistical mechan-ical model for our optical system displays the Hopfield-like Hamiltonian H [ s ] = − (cid:88) i (cid:54) = j J ij s i s j − N (cid:88) i =1 s i H i (19)with J ij ≡ N p (cid:88) µ =1 (cid:60) [ ξ ( µ ) i ¯ ξ ( µ ) j ] (20) H i ≡ √ N p (cid:88) µ =1 (cid:60) [¯ h ( µ ) ξ ( µ ) i ] = p (cid:88) µ =1 (cid:60) h ( µ )0 √ N + 1 N N (cid:88) j =1 ¯ ξ ( µ ) j ξ ( µ ) i (21)where the local fields 1 /N (cid:80) Nj =1 ¯ ξ ( µ ) j arise in the superpixel composition but their overall contribution to the Hamilto-nian is, actually, of O (1) and, thus negligible with respect to the extensive contribution of the pairwise spin-interaction.The fields h ( µ )0 can be introduced to identify Mattis states, i.e., those states most aligned along a given transmissionpattern ξ ( µ ) , and are to be sent to zero eventually, in order to be able to study also the ferromagnetic transition atlow enough α . We focus in this work to present the spin-glass behavior and we will present the analysis of Mattisstates elsewhere. Therefore, we set h ( µ )0 = 0, ∀ µ in the following.2To compute the average of the free energy over the distribution of the transmission matrix complex-valued elementswe need to replicate the system, i. e., β Φ = − lim N →∞ N E [ln Z ] ξ (22)= lim n → lim N →∞ E [ Z n ] ξ − nN where the replicated partition function for a given random transmission (given transmission matrix) reads Z n = ( N β ) p/ (cid:90) (cid:89) µa dr µa √ π dt µa √ π e − Nβ (cid:80) µa (cid:104) ( r ( µ ) a ) + ( t ( µ ) a ) (cid:105) (cid:88) { s } e F ( { ξ } , { (cid:126)s } ,(cid:126)m )with F ≡ β (cid:88) aµ (cid:34) r ( µ ) a (cid:88) i ξ ( µ ) R,i s ( a ) i + t ( µ ) a (cid:88) i ξ ( µ ) I,i s ( a ) i (cid:35) + O (1) . (23)After having averaged over the transmission matrix elements distribution, one obtains E [ Z n ] ξ = ( N β ) p/ (cid:90) (cid:89) µa dr µa √ π dt µa √ π e − Nβ (cid:80) µa (cid:104) ( r ( µ ) a ) + ( t ( µ ) a ) (cid:105) (cid:88) { (cid:126)s } E (cid:104) e F ( { ξ } , { (cid:126)s } ,(cid:126)m ) (cid:105) ξ with E (cid:2) e F (cid:3) ξ = (cid:89) µi e β σ (cid:80) ab (cid:16) r ( µ ) a r ( µ ) b + t ( µ ) a t ( µ ) b (cid:17) s ( a ) i s ( b ) i . Let us now define the vector of the spin state through the n replicas with an arrow (cid:126)s ≡ { s , . . . , s n } , to distinguishit from the configuration of all spins in a replica that we represent by a bold s . Properly rescaling some integrationvariables, defining the replica overlap q ab ≡ N (cid:88) i s ( a ) i s ( b ) i , (24)introducing the matrix M ab ≡ δ ab (1 − γ ) − γq ab (25)and defining the rescaled inverse temperature γ ≡ βσ we obtain the average replicated partition function E [ Z n ] ξ = ( βN ) p (1 − n ) (cid:90) (cid:89) a
1, yielding E [ Z n ] ξ (cid:39) e nN S [ q sp ,λ sp ] + O (cid:18) N (cid:19) β Φ = − lim N →∞ N E [ln Z ] ξ (cid:39) lim n → lim N →∞ E [ Z n ] ξ − nN = S [ q sp , λ sp ] (27)The saddle point equations, i. e., the self-consistency equations for the matrix order parameters q ab and λ ab , read λ ab = 2 αγ (cid:0) M − (cid:1) ab (28) q ab = (cid:80) { (cid:126)s } s a s b e H [ λ,(cid:126)s ] (cid:80) { (cid:126)s } e H [ λ,(cid:126)s ] = (cid:104) s a s b (cid:105) (29)3 Replica Symmetry
A first solution, certainly valid in the paramagnetic phase, can be put forward assuming replica symmetry, that is,assuming the overlap matrix form q ab = (1 − δ ab ) q (30)and, consequently, M ab = (1 − γ + γq ) δ ab − γq, a (cid:54) = b (31)The saddle point equations of the RS solution read λ = 2 αγ q (1 − γχ ) (32) q = (cid:90) dz √ π e − z / tanh (cid:16) √ λz (cid:17) (33)with χ ≡ − q , the replica symmetric free energy is β Φ = λ − q ) − (cid:90) dz √ π e − z / ln 2 cosh (cid:16) √ λz (cid:17) + α ln(1 − γχ ) − αγ q − γχ (34)In particular, the paramagnetic free energy will be ( q = 0) β Φ PM = − ln 2 + α ln(1 − γ ) (35) Stability and phase transition
The stability eigenvalues of the replicated free energy in the RS Ansatz are three: the so-called replicon eigenvalueΛ = = 12 α (cid:18) − γχγ (cid:19) − (cid:90) dz √ π e − z / (cid:104) − tanh (cid:16) √ λz (cid:17)(cid:105) (36)and the longitudinal and the anomalous eigenvalues, that coincide for n → , = 12 αγ (cid:104) (1 − γ ) + 4 qγ (1 − γ ) + 3 q γ (cid:105) − (1 − q + 3 r ) (37)On the paramagnetic phase ( q = 0) the stability conditions become one Λ = Λ , ≥ Tσ > √ α. (38)that is the critical temperature of the transition from the paramagnetic to the spin-glass phase, that depends on theratio α between targets and spins. The phase transition turns out to be continuous in the order paramater q . Aknowledgements
M.L. thanks “Fondazione CON IL SUD”, “Grant Brains2south, Project “Localitis”. Manuscript submitted toPNAS. [1] N. C. Harris, G. R. Steinbrecher, M. Prabhu, Y. Lahini, J. Mower, D. Bunandar, C. Chen, F. N. Wong, T. Baehr-Jones,M. Hochberg, et al., Nature Photonics , 447 (2017). [2] C. Sparrow, E. Martin-Lopez, N. Maraviglia, A. Neville, C. Harrold, J. Carolan, Y. N. Joglekar, T. Hashimoto, N. Matsuda,J. L. O Brien, et al., Nature , 660 (2018).[3] X. Lin, Y. Rivenson, N. T. Yardimci, M. Veli, Y. Luo, M. Jarrahi, and A. Ozcan, Science , 1004 (2018).[4] A. Peruzzo, J. McClean, S. Alan, and J. L. O Brien, Nature communications , 4213 (2014).[5] C. Roques-Carmes, Y. Shen, C. Zanoci, M. Prabhu, F. Atieh, L. Jing, T. Dubˇcek, V. ˇCeperi´c, J. D. Joannopoulos,D. Englund, et al., Conference on Lasers and Electro-Optics p. FTu4C.2 (2019).[6] N. H. Farhat, D. Psaltis, A. Prata, and E. Paek, Applied optics , 1469 (1985).[7] S. Edwards and P. Anderson, J. Phys. Lett. , 965 (1975).[8] M. M´ezard, G. Parisi, and M. Virasoro, Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications,vol. 9 (World Scientific Publishing Company, 1987).[9] A. P. Young, Spin glasses and random fields, vol. 12 (World Scientific, 1998).[10] D. S. Fisher and D. A. Huse, Phys. Rev. B , 386 (1988).[11] L. Leuzzi, G. Parisi, F. Ricci-Tersenghi, and J. J. Ruiz-Lorenzo, Phys. Rev. Lett. , 267201 (2009).[12] T. Temesvari, Nucl. Phys. B , 534 (2010).[13] M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, D. I˜niguez, A. Maiorano,E. Marinari, V. Martin-Mayor, et al. (Janus Collaboration), Phys. Rev. Lett. , 267203 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.267203 .[14] D. J. Amit and D. J. Amit, Modeling brain function: The world of attractor neural networks (Cambridge university press,1992).[15] N. Ghofraniha, I. Viola, F. Di Maria, G. Barbarella, G. Gigli, L. Leuzzi, and C. Conti, Nature communications , 6058(2015).[16] F. Antenucci, Statistical physics of wave interactions (Springer, 2016).[17] M. Halasz, A. Jackson, R. Shrock, M. A. Stephanov, and J. Verbaarschot, Physical Review D , 096007 (1998).[18] H. Erik, Optical spin glasses : a new model for glassy systems (Master Thesis, Univ. Sapienza, Oct. 2018).[19] D. Pierangeli, M. Rafayelyan, C. Conti, and S. Gigan, arXiv preprint arXiv:2006.00828 (2020).[20] J. J. Hopfield, Proceedings of the National Academy of Sciences , 2554 (1982), ISSN 0027-8424.[21] J. J. Hopfield, Proceedings of the National Academy of Sciences , 3088 (1984), ISSN 0027-8424.[22] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Physical Review Letters , 1530 (1985).[23] D. Amit, H. Gutfreund, and H. Sompolinsky, Phys. Rev. A , 1007 (1985).[24] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Annals of Physics , 30 (1987).[25] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, The journal of chemical physics , 1087(1953).[26] A. Crisanti and L. Leuzzi, Phys. Rev. B , 144301 (2007).[27] W. G¨otze, Complex Dynamics of Glass-Forming Liquids, A Mode-Couplig Theory (Oxford University Press, 2009).[28] F. Caltagirone, U. Ferrari, L. Leuzzi, G. Parisi, F. Ricci-Tersenghi, and T. Rizzo, Phys. Rev. Lett. , 085702 (2012).[29] F. Caltagirone, U. Ferrari, L. Leuzzi, G. Parisi, and T. Rizzo, Phys. Rev. B , 064204 (2012).[30] U. Ferrari, L. Leuzzi, G. Parisi, and T. Rizzo, Phys. Rev. B , 014204 (2012).[31] S. Popoff, G. Lerosey, R. Carminati, M. Fink, A. Boccara, and S. Gigan, Physical review letters , 100601 (2010).[32] D. O. Hebb, The organization of behavior: a neuropsychological theory (Science Editions, 1962).[33] I. M. Vellekoop and A. Mosk, Optics letters , 2309 (2007).[34] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, Phys. Rev. Lett. , 100601 (2010).[35] D. Mattis, Phys. Lett. , 421 (1976).[36] J. Buhmann and K. Schulten, Biological cybernetics , 313 (1987).[37] W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and S. C. Glotzer, Physical review letters , 2827 (1997).[38] C. W. Beenakker, Reviews of modern physics69