Optical Phase Transitions in Photonic Networks: A Spin System Formulation

We investigate the collective dynamics of nonlinearly interacting modes in multimode photonic settings with long-range couplings. To this end, we have established a connection with the theory of spin networks. The emerging"photonic spins"are complex, soft (their size is not fixed) and their dynamics has two constants of motion. Our analysis shed light on the nature of the thermal equilibrium states and reveals the existence of optical phase-transitions which resemble a paramagnetic to a ferromagnetic and to a spin-glass phase transitions occurring in spin networks. We show that, for fixed average optical power, these transitions are driven by the type (constant or random couplings) of the network connectivity and by the total energy of the optical signal.

An emerging framework is in photonics, where light propagation in non-linear multimode optical structures have recently attracted a lot of attention [9][10][11][12][13][14][15][16][17][18][19] . On the fundamental side there are many unanswered questions associated with the energy exchange between the modes and the role of the underlying spatio-temporal complexity, originating from the disorder, the network topology and the complex intermodal interactions. Bruteforce computational attempts to answer these questions are either impossible (due to the large number of degrees of freedom involved) or unsatisfactory as far as the understanding of the underlying physics that dictates the energy redistribution. At the same time, there is a pressing need from modern technologies to develop theoretical tools that will allow us to tailor the intermodal energy exchange and harvest it to our advantage. If this endeavor is successful, it will give rise to a next generation of high power light sources 17 , high-resolution imaging schemes [20][21][22] , and high-speed telecommunication systems [23][24][25] .
In this regard, a number of recent papers have promoted well established equilibrium 9,11,19,[26][27][28] and nonequilibrium 15,29-32 thermodynamics techniques as a viable theoretical toolkit which can be used to address some of the above challenges. A decisive step along these lines has been achieved recently by the authors of Ref. 19,26 which, assuming thermal equilibrium conditions and weak non-linearity, have established a comprehensive * corresponding author:tkottos@wesleyan.edu Figure 1: (Color online) A variety of photonic nonlinear multimode networks, with long-range modal couplings, whose field dynamics is modeled by an effective coupled mode theory of the type given by Eqs. (1,2). (a) A network of coupled micro-resonators 40 ; (b) A photonic network of single-mode fibers coupled together via couplers (optical splitters) 41,42 ; (c) A multicore fiber 25 ; (d) A multimode fiber 25,43,44 ; (e) Deformed multimode (micro-) resonator with underlying chaotic dynamics [45][46][47] ; (f) A network of coupled "soft" (size-modulated) spins. The corresponding Hamiltonian has similarities but also crucial differences with the Hamiltonian that describes a photonic multimode network.
optical thermodynamics formalism allowing us to design potentially novel classes of high-power multimode optical structures or efficient cooling schemes. For weak nonlinearity, the specific nature of the nonlinear mode interaction (e.g. Kerr or saturable or thermal nonlinearities) is irrelevant. Its role is important for the thermalization process but not for the properties of the equilibrium state.
The implementation of statistical thermodynamics methods in modern photonics opens up a new arena where ideas and concepts from statistical mechanics can be transfered to optics and utilized for light control.
A prominent example is the notion of phase-transitions characterized by changes in the properties of a system as an external parameter varies 2 . Maybe, the most celebrated area of physics where such phenomena have been extensively studied is the theory of spin network models [33][34][35] . These studies indicated that the phase transition can be understood as a competition between the interactions among spins, which facilitates order, and the thermal fluctuations (the entropy) causing random disturbances. In many occasions, these phase transitions are associated with symmetry breaking phenomena. To understand them, one has to analyze the nature of the thermal equilibrium state of the interacting system. In this paper, we will show that there are certain analogies between the statistical thermodynamics of spin models and multimode photonic networks. Inspired by these analogies we will ask questions like: are the phases of the electromagnetic field correlated or are they entirely random? Is the optical power distributed more or less evenly between all modes of the entire system, or some finite fraction of the total optical power can reside in a single mode? How the topology of the inter-mode connectivity and the randomness in the coupling affect the nature of the thermal state?
It turns out that these questions can be related. In previous works, for example, it was shown that macroscopic occupation of a single photonic mode does occur and, moreover, it can happen even in linear systems (we stress again that weak non-linearity can be neglected only in the equilibrium state, while it is crucial in the thermalization process) 27,29,36 . It is quite remarkable, thus, that a purely classical system exhibits a phenomenon alike BEC transition in a quantum Bose gas. Actually, it has been a number of experiments demonstrating that, in the course of propagation along the fiber, optical power "flows" towards the lower modes 18,[37][38][39] .
In this paper, we demonstrate that the connectivity of an optical structure is an important factor in its thermalization process: it affects the type of the optical phase transitions and the nature of the thermal equilibrium state. The question is not only of fundamental importance; it pertains also to recent photonic developments where networks with complex connectivities, can be realized, see Figs. 1a-e. Specifically, we show that in the case of long-range couplings, the nonlinearity is instrumental for achieving optical phase transitions. In fact, by identifying an order parameter that is equivalent to the magnetization in spin-network models, we are able to show both theoretically and numerically the existence of a ferromagnetic-to-paramagnetic phase transition, analogous to the one occurring in spin systems. When disorder is introduced into the couplings, the system might undergo another type of a transition, namely to a spinglass phase; much as in the case of frustrated coupled spins. Although these analogies between photonics multimode networks (Figs. 1 a-e) and spin-networks ( Fig.  1f) is useful, one needs to keep in mind that the two problems have important differences. Specifically our "photonic spins" are complex dynamical variables (amplitudes of the electric field). Moreover, they fluctuate not only in their direction but also in their size. We expect that the analogies drawn from our study will bring together two seemingly different areas: statistical mechanics of spin networks and light transport in nonlinear multimode settings. This cross-fertilization will, hopefully, allow the development of better design strategies for the control of light transport in multimode photonic networks.
The structure of this paper is as follows. In the next section II we discuss the general statistical thermodynamics formalism associated with the analysis of optical thermal equilibrium states. Special attention is given to the case of weak nonlinearities where we derive the occupation number statistics. In the next section III we analyze a class of complex multimode photonic networks with long-range connectivity. Two cases are discussed in detail: the case of constant couplings and the opposite case of random couplings. We show that under specific conditions these systems demonstrate optical phase transitions from ferromagnetic to paramagnetic and spinglass phases. Finally, our conclusions are discussed at the last section IV.

II. GENERAL FORMALISM
The dynamics of nonlinear multimode photonic networks in Fig. 1 can be modeled using the framework of time dependent coupled mode theory. The associated equations are where ψ l is the degree of freedom (the complex amplitude) at node l of the "network", and J lj = J * jl is the connectivity matrix that dictates the couplings among the nodes. We will typically assume zero self-coupling terms J ll = 0. Finally, the last term in Eq. (1) describes the nonlinearity due, for instance, to Kerr effect. On a formal and general level ψ l = l|ψ are components of the electric field (with some fixed polarization) in some basis of orthonormal modes {|l } (the "basic modes"). The choice of the set of these modes depends on the problem at hand. For instance, for the case of coupled single-mode microresonators 40 (Fig.1a) the index l labels the resonators and the "basic mode" |l is the eigenmode of the l-th resonator, decoupled from the rest of the network (we treat the resonators as structureless point objects). In this framework, the coefficients J lj = J * jl (for l = j) represent evanescent couplings between different resonators and J ll = ω l is the eigenfrequency of resonator l (with nonlinearity neglected), which in case of identical resonators can be set to be zero i.e. ω l = 0(l = 1, · · · , N ). The same interpretation applies to a fiber network 41,42 Fig. 1b and to a multicore fiber 25 , Fig. 1c, where now l labels the single mode fibers. In Fig. 1c the propagation direction z plays the role of time.
It is important to clearly distinguish between the "basic modes", like an eigenmode of an isolated resonator in Fig.1a, and the eigenmodes of the entire structure, i.e. the stationary solutions of the entire system of coupled equations (1) (with the nonlinearity neglected). The latter are often referred to as "supermodes". For instance, when we are talking about condensation of the optical power in a single mode, we mean of course the supermode (with the lowest energy) and not an eigenmode of a single resonator. Sometimes, where no confusion can arise, we will use the term "mode" instead of "supermode". The "supermodes" f α (l) form a complete basis, so the field ψ l (t) can be expanded as ψ l (t) = α C α (t)f α (l), reducing Eq. (1) to with As we will see below, this mode representation of the evolution Eq. (1) is particularly useful when the nonlinearity is weak. Equations (1) and (2) also describe multimode optical fibers 25,43,44 or multimode resonators [45][46][47] , Figs. 1d,e respectively. In this case the "basic modes" |l are the eigenmodes of an ideal, undeformed fiber (or resonator) while J lj are the couplings among these modes, due to various perturbations (deformations of the ideal system, possible impurities, etc). Note that unlike the previous case , when the basic modes were localized in space (on a single resonator), now they extend over the entire structure.
The equation of motion (1) is derivable from the energy functional (the Hamiltonian) In the course of time the total energy E and the total optical power are conserved. Finally, we always assume that both the total power A, and the energy E, are extensive quantities, proportional to the number of modes N , as appropriate for thermodynamics.

A. The Problem of Thermalization
An important question is whether an isolated system of interacting modes eventually thermalizes, i.e. reaches an equilibrium state which can be described by just two parameters -the inverse temperature β and the chemical potential µ, which in turn are determined by the energy E and the total power A of the initial preparation. If such equilibrium state is reached, then the system can be analyzed using the well established methods of statistical mechanics and thermodynamics.
For example, a statistical mechanics description of the system of Eqs. (1) is achieved by calculating the classical grand -canonical partition function Z where the Lagrange multipliers β = 1/T and µ have been introduced (in analogy with the inverse temperature and the chemical potential) to ensure conservation (on average) of E and A respectively (see Eqs. (4,5)). Specifically the relation between the microcanonical quantities which describe the average optical power a and averaged energy density h per mode and the grant canonical quantities µ, β is given by Using the partition function as a starting point we can next calculate the thermodynamic potential and from there, the entire "optical thermodynamics" can be derived. For instance, the entropy is S = − ∂Ω ∂T µ . Finally we note that the problem of thermalization in non-linear lattices, under time evolution defined in Eq. (1) (primary in one spatial dimension and with J lj restricted to nearest neighbors only) has been also addressed in the framework of statistical mechanics [48][49][50][51] . In this studies, it has been pointed out that thermalization occurs only in a certain region of the (E, A)-plane. These studies indicated that for fixed total norm A, the system thermalizes only if its energy is not too large 48 . Otherwise, the equilibrium Gibbs distribution, Eq. (6), cannot be reached-the system is said to belong to the non-Gibbsian, or negative temperature region. The maximum value h = h max for a given average norm a per site corresponds to β −→ 0 (high temperature) and µ −→ −∞, and it is It turns out that for sufficiently large E the norm cannot spread uniformly over the entire system and highamplitude peaks of ψ l (breathers 52 ) emerge. Below, we will confine our analysis to the domain where the temperature is positive and thermalization can be achieved.

B. Thermal Equilibrium in the Case of Weak Non-linearities
In many optics applications the non-linearity is considered weak. It is of course essential in the mode-mixing process (see Eq. (2)), needed for reaching equilibrium. The total energy and power, however, are dominated by the linear term in the Hamiltonian, i.e. to a good approximation Here |C α | 2 is the (normalized) optical power in mode α.
Using these expressions Christodoulides et al. 19,26 were able to develop a kind of "optical thermodynamics", identifying the optical analogy of entropy, equation of states and other quantities. Below we briefly summarize this theory, using the grand-canonical formulation Eq. (6) discussed above. Assuming that the system has thermalized, its grand partition function, non-linearity being neglected, is given by the product of independent modes contributions where β = 1/T and µ can be found from the constraints in Eq. (11). It immediately follows from Eq. (12) that the average of optical power in mode α obeys the Rayleigh-Jeans distribution 19,26,27,[29][30][31] (see also Ref. 53 ) From Eq. (13) we can calculate the thermodynamic potential Eq. (9) which takes the form All other thermodynamic variables follow from Eq. (14). For example, the entropy (up to an irrelevant constant) is S = N α=1 ln (n α ) which, in equilibrium, coincides with the expression derived in Ref. 19 by counting the number of ways in which a large number of "packets of power" can be distributed over the N modes. The expression for S served in 19 as the starting point for the development of "optical thermodynamics". For instance, one can derive the following equation of state 19 E − µA = N T which connects three extensive quantities (E, A, N ) to the two intensive variables (µ, T ).

C. Fluctuations in the Case of Weak
Non-linearities Next, we briefly discuss the fluctuations of the optical power n α in the mode α. If the nonlinearity, i.e. the intermode interaction, in the thermal equilibrium state is negligibly small 19,26,27,[29][30][31] , then within the grand canonical treatment the probability density for n α is where Z α is the normalization factor. This yields Eq. (13) for the average valuen α and ∆n 2 α = (n α ) 2 for the variance. The same results can be obtained, in even simpler way, if one uses the expression Ω α = −T ln πT εα−µ for the contribution of mode α to the grand potential Ω (see Eq. (14)) and the standard formulas 54n α = −∂Ω α /∂µ, ∆n 2 α = T ∂n α /∂µ. Thus, the standard deviation ∆n 2 α 1/2 ≡ σ α comes out to be equal to the average optical powern α . For a mode with macroscopic occupation (n α 1) this result looks paradoxical. The "paradox", however, is well known in the theory of Bose-Einstein condensation and it is resolved by observing that we have here one of the rare cases when the canonical and the grand canonical ensembles yield different results 54 . Indeed, in the experiment, as well as in our numerical simulations, the total optical power α n α = A is strictly conserved while in the grand canonical treatment it is conserved only on the average. This is perfectly fine for calculating various average quantities but not for the fluctuations. When the conservation law α n α = A is strictly enforced (canonical ensemble), the large unphysical fluctuations in a macroscopically occupied mode disappear (for instance, at T = 0, when the entire power A is located on a single mode, there are no fluctuations at all).
Note, however, that at high temperatures, when there are many modes populated withn α 1, the result σ α =n α does hold for such modes. This is because the constraint α n α = A ∼ N on the total power does not significantly affect fluctuations in a single mode with n α 1 (the other modes serve as an "environment" for the mode α). One should be aware of these large fluctuations when interpreting the numerical or the experimental data.

III. MULTIMODE OPTICAL SYSTEMS WITH LONG RANGE COUPLING
We consider the connectivity matrix J lj (see Eq. (1)) of the following form: where the first term describes a fully connected network, with equal couplings, while the second term introduces some randomness into the couplings. In the case of chaotic or disordered networks 55 the couplings B lj = B jl are given by a Gaussian distribution with zero mean and a unit standard deviation. We point out that, unless stated otherwise, in all simulations below the random matrix elements B i,j remain the same (fixed) for a specific set of parameters N, J 0 , χ dictating the Hamiltonian of the photonic network. For H in Eq. (4) to remain extensive, in the large N limit, one has to scale the coupling strengths with N , as written in Eq. (15). The coupling matrix Eq. (15) gives rise to two distinct terms in the Hamiltonian Eq. (4). The first term resembles the Curie-Weiss (CW) model 56 , where N Ising spins are coupled to each other by constant distanceindependent interactions. The second term resembles the Sherrington-Kirkpatrick (SK) model for a spin glass 57,58 , where the couplings are completely random. We write "resembles" because our model differs from the well studied CW and SK models in three important respects: First, our dynamical variables ψ l are complex, as appropriate for the complex amplitudes of electric field, and they can be treated as real two-component "spins". Second, since |ψ l | 2 is not restricted to some fixed value, our spins fluctuate not only in their direction but also in their size (note though, that the non-linearity does not allow for too wild fluctuations in size). And, third, the dynamics of our "photonic spins" conserve not only the energy (as in standard spin systems) but also the optical power, see Eq. (5). This second conservation law introduces novel features into the characteristics of the thermal equilibrium state; for instance a condensation of optical power in a single mode.
Below we will distinguish between the two limiting cases corresponding to a photonic network with equal couplings (σ = 0) and to its "random" coupling analogue (J 0 = 0). We will also briefly discuss the case where the connectivity matrix Eq. (15) contains both terms. We will show that our photonic network exhibits various phases depending on the disorder strength of the coupling constants, the energy and optical power (h, a) of the initial preparation and the strength of the nonlinearity parameter χ.

A. Numerical Method
The thermalization process of an initial state {ψ n }(n = 1, · · · , N ) has been investigated numerically using a high order three part split symplectic integrator scheme [59][60][61] for the integration of Eq. (1). The method conserved, up to errors O(10 −8 ), the total energy Eq. (4) and the optical power Eq. (5) of the system. These quantities have been monitored during the simulations in order to ensure the accuracy of our results.
We have focused our interest on the electric field amplitudes ψ n (t) and the supermode amplitudes C α (t) which can be evaluated from the projection of ψ n (t) on the supermode basis, see Eq. (2). Knowledge of C α (t) (or ψ n (t)) allows us to calculate various thermodynamic quantities Q by making a time-average and invoking ergodicity.
In practice, the approach to a thermal equilibrium state involves a long time propagation of an initial preparation {ψ n }(n = 1, · · · , N ). Typical integration times were as long as 4×10 8 coupling constants. After an initial and confirmed its convergence to a steady state value. An example of such simulations for the nodal powers a n (t) ≡ |ψ n (t)| 2 and the supermode power |C α | 2 are shown in Fig.  2a,b respectively. Failure to reach a steady state value indicated that the system did not reached the thermal equilibrium state.
In the case of weak non-linearity the numerical results for |C α | 2 have been compared against the theoretical predictions Eq. (13). A good agreement between them serves as a confirmation that the system reached a thermal equilibrium state. A disagreement between the numerics and the theoretical predictions of Eq. (13), indicates that the thermal equilibrium has not been reached. Instead, the system might have reached a metastable state as happens in the case of a spin-glass behavior 62,63 . Given enough time (large relaxation times), of course, the system will reach the global free energy minimum.
In all cases, the initial conditions were generated by considering the field amplitudes {|ψ n |}, and the phases φ n being random numbers in the intervals [1 − δ, 1 + δ] and [−π, π] respectively. Out of a large number of {ψ n } configurations we have chosen only the ones that satisfy the energy and normalization constraints that define the specific state, see Eqs. (4,5) respectively. Finally, in all our simulations below we have used the normalization condition A = n |ψ n | 2 = N (i.e. a = 1, see Eq. (7)) associated with the total power.

B. Equal Coupling Networks
We start our analysis with the equal-coupling photonic network. In this case each node is coupled to all other nodes by a hopping amplitude of equal strength. The Hamiltonian that describes this system is Eq. (4) with a connectivity matrix given by Eq. (15) with σ = 0. We have where l and j run over all N sites with the only constraint that l = j.
In the absence of non-linearity the statistical mechanics of the system in Eq. (17) is trivial. The connectivity matrix (Eq. (15) with σ = 0) can be easily diagonalized. In this case the Hamiltonian Eq. (17) (with χ = 0) has one non-degenerate eigenvector (supermode) with eigenfrequency ε 1 = −(N −1)J 0 /N and (N −1)-fold degenerate eigenvectors with frequency ε α = J 0 /N (α = 2, · · · , N ). Therefore, in the large-N limit, only the lowest non-degenerate mode contributes to the total energy E and, since the latter quantity is required to be extensive, it is clear that a final fraction of the total optical power A must reside in that mode. A simple calculation, based on relations Eq. (11) and the expression Eq. (13) yields the following result: Since the total optical power is A = aN and the total energy is E = hN , then, in the large-N limit, the resulting chemical potential and the temperature are Since, obviously, J 0 a ≥ |h| must hold, the linear model does not allow for either negative temperatures or for a transition: a finite fraction (|h|/J 0 a) of the total power A is condensed into the lowest mode 64 . A refined analysis, where the finite size effects are taken into consideration, leads to the following exact expressions for the chemical potential and the temperature 20) which, in the limit of N 1, are nicely matching the results in Eq. (18). In Fig. 3 we have compared these theoretical predictions with the values of (µ, T ) that we have extracted from our numerical simulations with the Hamiltonian of Eq. (17), using two multimode photonic networks with N = 8 and 64 modes. The nice agreement indicates that these systems have reached a thermal equilibrium state. At the inset of the same figure we also report |C α | 2 t (see Eq. (16)) by making use of the scaled variablesn α where in the second equation above α = 2, · · · , N . The right hand side of Eqs. (21) has been evaluated using Eq. (13) together with Eqs. (11). It is important to point out that thermalization has been achieved even for the system with relatively small number of nodes N = 8.
The presence of non-linearity changes the picture completely and provides us with an example of a (mean field) optical phase transition from an ordered to a disordered phase. In the former phase the amplitudes ψ l on different nodes are correlated while in the latter phase the nodes become essentially decoupled from each other. The ordered (disordered) phase corresponds to low (high) temperature, i.e. to small (large) internal energy E.
The ground state of the Hamiltonian (17) corresponds to a uniform field configuration {ψ l } = A /N (1, 1, . . . , 1), where the normalization is such that the total optical power is l |ψ l | 2 = A. All the "spins" in this state point in the same direction. Note that the ground state is highly degenerated: one can rotate all the spins by an angle θ, i.e. multiply the field by an overall phase factor exp(iθ). The ground state energy is corresponding to T = 0. In the opposite limit, T → ∞, the kinetic energy (hopping) can be neglected and the system reduces to a set of uncoupled nonlinear oscillators. The probability density to find a node with optical power |ψ| 2 ≡ I is p(I) ∼ exp(βµI − 1 2 βχI 2 ) and a simple calculation yields an average thermal energy equal to χa 2 , in the T → ∞ limit. Thus, the maximal energy of the system (see Eq. (10)) is E max = χa 2 N indicating that for fixed a the thermalization occurs for E min < E < E max . One can, of course, endow the system with an energy larger than E max , for instance, by putting all the norm A on a single site (in which case the energy E would become "super-extensive", proportional to N 2 ). We do not consider, however, initial preparation with E > E max since we do not address the "non-Gibbsian" region in the (E, A)-plane (see previous discussion and also Refs. 48,49 ).
Next, we proceed with the calculation of the partition function Eq. (6) which is performed using the saddlepoint method. It is convenient to write the complex amplitudes ψ l = q l + i p l where (q l , p l ) is a pair of real variables. The Hamiltonian (23) can be interpreted in terms of interacting two-component spins. We write l,j q l q j = ( l q l ) 2 and use the identity and similarly for l,j p l p j . The integrals over q l , p l in the partition function now factorize into a product of integrals, each involving only variables for one node. Furthermore, in the large N -limit, the grand-canonical partition function Z is dominated by a saddle-point, which can be interpreted as the order parameter, i.e. the average fieldψ (the magnetization in the statistical mechanics language). Actually, there is a whole family of saddle points, distinguished from each other by an overall phase. Choosing one saddle point of the family amounts to breaking the rotational symmetry in the spin space, obtaining a non-zero value for the order parameter. We choseψ ≡ m to be real. Skipping all calculation details, we only give the final equation for m where the function F (m) is given by the integral F (m) = 2π´∞ 0 rdrI 0 (2βJ 0 mr) exp βµr 2 − 1 2 χβr 4 and I 0 (x) is the modified Bessel function of order zero.
For small m, Q(m) is a linear function of m, Q(m) = s · m, and the slope s determines whether Eq. (25) has a nontrivial solution m = 0. Such solution exists only if s > 1. The slope can be calculated using the small argument expansion I 0 (x) = 1 + 1 4 x 2 , and it can be written as s = βJ 0 K 3 /K 1 where an integral K n is defined as K n = ∞ 0 drr n exp(βµr 2 − 1 2 χβr 4 ). Let us, as an example, fix µ at the value zero and study s as a function of the temperature T = 1/β. For µ = 0 a simple expression for the slope s is obtained, namely s = J 0 (2/πT χ) 1/2 . The value s = 1 yields the critical temperature T c = 2J 2 0 /πχ. For temperature T slightly below T c one finds, by keeping the term of order m 3 in the expansion of Q(m), the standard mean field result m ∼ (T c − T ) 1/2 . When the temperature decreases further, towards T = 0, the magnetization increases and approaches the maximal value, corresponding to the fully ordered ground state. Taking again µ = 0 as an example, one obtains from Eq. (25) that, in the β → ∞ limit, m → J 0 /χ. This result becomes transparent when written in terms of the optical power A. Indeed, at T = 0 the total power resides in the (fully correlated) ground state so that the magnetization per site is m = A/N ≡ √ a. To connect a to µ we have to use the expression Eq. (22) for the ground state energy which yields µ = 1 N ∂E/∂a = −J 0 + χa. For µ = 0, one obtains m = √ a = J 0 /χ. Thus, Eq. (25) is well suited for studying m, as a function of β, for a fixed value of µ. In experiment, however, one usually controls the optical power A, rather than the conjugate variable µ. Calculating analytically m as a function of β for fixed A is more involved than the above calculation for fixed µ, and we do not attempt it in the present paper. Instead, in Fig. 4 we present some numerical results, serving a double purpose: first to verify that the system, with the appropriate initial preparation, indeed thermalizes and then to study its properties in the thermal equilibrium state as a function of the averaged energy density h. To this end, we evaluate numerically the time-averaged magnetization |m| t (see Fig. 4a) defined as In our numerics, we consider moderate values of the nonlinear parameter χ = 1 and coupling constant J 0 = 1.2.
At the ground state h = h min all "spins" have the same orientation. As a result, the magnetization acquires its maximum value |m| t = 1 indicating a ferromagnetic behavior. For higher h-values the magnetization decreases and at h → h c ≈ 0.75, which is smaller than h max = 1, it acquires a constant value |m min | t = O(1/ √ N ) (see inset of Fig. 4a), associated with finite size effects. Further numerical analysis indicates that in the thermodynamic limit of N → ∞ the magnetization, as a function of h, approaches zero following a square-root behavior |m| t ≈ 0.85 √ h c − h (see bold black line in Fig. 4a). Such behavior is characteristic of a second-order phase transition, from a "ferromagnetic" to a "paramagnetic" (optical) phase.
In Fig. 4b we show the numerical results for the timeaveraged optical power in the ground state supermode |C 1 | 2 t versus the averaged energy density of the initial beam. In the simulations, we have considered the same photonic networks as above with χ = 1, J 0 = 1.2, and different numbers of nodes N . The numerical findings are reported using the scaled variablen 1 , see Eq. (21). We find that for low averaged energy densities h, the optical power of the ground state is |C 1 | 2 t ≈ N indicating a condensate. As h increases the condensate depletes and eventually at h = h c , corresponding to the ferro-paramagnetic transition, it is completely destroyed i.e.n 1 = 0. At the same figure we also report for comparison the theoretical value ofn 1 = −h/J 0 (black solid line), applicable for the case of weak non-linearities, see Eq. (13). We stress once more that the condensation transition analyzed above is due solely to the nonlinearity unlike the previous studied cases where the transition occurred already in the linear system 29,32 .

C. Random Coupling
Next, we analyze the effect of disorder in the coupling constants of the photonic network i.e. σ = 0 in Eq. (15). First, we consider the case of extreme disorder where J 0 = 0. In this case the coupling constants are entirely random, with equal probabilities to be positive or negative, i.e. the system is completely "frustrated". Despite the vast literature on spin-networks, this model with complex, "soft spins" has not been studied up to now. Of course, certain analogies with the Sherrington-Kirkpatrick model can still be instructive. In the latter case, there is a transition from a paramagnetic phase to a spin-glass phase when the temperature (the energy E of the system) decreases. One characteristic distinction between the two phases is that for the spin-glass the thermalization time is much longer than for the paramagnet. Moreover, for large N a spin-glass does not reach a full thermal equilibrium in any reasonable time, and the system gets stuck in one of the many metastable states.
Our simulations for the random coupling multimode photonic network are presented in Fig. 5. For a weak nonlinearity χ = 0.01, the equilibrium optical powers n α of the supermodes (of the linear problem) are given by Eq. (13). In the simulation we evaluated the set {|C α (t)| 2 } as a function of time, extracted their timeaverage Eq. (16), and use their comparison to Eq. (13) as a criterion for thermalization. We find that for energy h ≈ −0.685, see Fig. 5a, the system gets eventually close to thermal equilibrium at times t ≈ 5 × 10 7 (in units of standard deviation of the coupling elements). After this time the occupation numbers change only slightly.  Fig. 3. The condensation transition, corresponding ton1 = 0, occurs for the same value of h = hc ≈ 0.75 as the one that signifies the transition from a ferromagnetic to a paramagnetic behavior in (a). In all cases we have considered an initial averaged optical power a = 1, coupling constant J0 = 1.2 while the nonlinearity is χ = 1. In this case, the maximum energy density is hmax = 1.
For higher energies (not shown) the thermalization time becomes shorter (for example for h ≈ −0.15 the thermalization time for a network of N = 16 nodes was ∝ 10 4 ). On the other hand, for energy h = −1 (see Fig. 5b), the optical powers { |C α | 2 t } are far away from {n α } even after a fairly long time t = 5 × 10 7 and, moreover, they do not show any significant change with respect to the results extracted for shorter times t = 10 5 (filled black circles in Fig. 5b). We have confirmed that the lack of thermalization (for any reasonable large time) is typical for other initial preparations (with the same h). This is a typical behavior of a spin-glass.
Indeed, the most important signature of a spin-glass, from which the term itself was derived, is that at low temperatures the directions of spins at various sites get frozen in some random configuration (metastable state). For our "optical spin glass" such behavior implies that the average values of the complex amplitudes {ψ l }, and in particular the phases {θ l }, form a random set. The "average" here refers to the thermal statistical average, for a fixed realization of the disorder. One could expect that in a numerical simulation averaging over a statistical ensemble can be replaced by the time average. However, Figure 5: (Color online) Comparison between the timeaveraged supermode optical powers for a multimode photonic network whose dynamics is described by the Hamiltonian Eq. (4) with random connectivity matrix Eq. (15) with J0 = 0 and σ = 1, and the theoretical results (red squares) of Eq. (13) for weak disorder. Black circles correspond to moderate time evolutions with t = 10 5 while green diamonds correspond to larger time evolutions with t = 5×10 7 . The initial state has (a) high averaged energy density per mode h = −0.685 (corresponding to high temperatures) or (b) low averaged energy density per mode h = −1 (corresponding to moderate/low temperatures). In both cases the average optical power is a = 1. In these simulations, the number of modes is N = 16 and the nonlinear parameter is χ = 0.01. due to a dynamic overall phase in the time evolution defined in Eq. (1), the time-average of the phase, θ l , at any site l, amounts to zero. Therefore, in order to distinguish a spin glass from a paramagnet, we use the following criterion: Let us denote by superscripts α, β two initial preparations, with the same total energy E and optical power A. Their time evolution is given by {ψ α l (t)} and {ψ β l (t)} respectively. The quantity ζ(t) = 1 N l (ψ α l (t)) * ψ β l (t) is a measure of the overlap between the two evolutions, at time t. In the paramagnetic phase ζ(t) decreases with time, approaching zero, because (in the large N -limit) the two evolutions become completely uncorrelated.
It is appropriate to consider many initial preparations, i.e. many (α, β)-pairs, and treat the real and the imaginary parts of ζ(t) as statistical quantities with probability distribution P (Re(ζ)) and P (Im(ζ)). In the long time limit, and for large but finite N , these distributions are expected to have a manifestly different form in the two phases: For a paramagnet P (Re(ζ)) and P (Im(ζ)) should be narrow distributions (with a width approaching zero when N → ∞), centered around ζ = 0.  around ζ = 0. When, on the other hand, we consider the same distributions for a set of initial preparations with low value of h = −1.4, we observed an entirely different behavior for P (Re(ζ)) and P (Im(ζ)) (see brown highlighted histogram in Figs. 6a,b). Namely, they become broad and almost flat, covering the whole allowed interval i.e. −1 < Re(ζ), Im(ζ) < 1. We stress that the above simulations were performed for a given realization of disorder and for the same energy h = −1.4. Only the initial preparations have been randomly chosen. We interpret the "flatness" of P (Re(ζ)) , P (Im(ζ)) as a signature of many metastable states, typical of a spin-glass 62,63 . It is natural to ask what happens to the network at low averaged energy densities h when both terms in the connectivity matrix Eq. (15) coexist, i.e. J 0 = 0 and σ = 0. In Fig. 7 we report the dependence of the magnetization |m| t versus the control parameter x = J 0 /σ and for h ≈ 0.88 1 . In the simulations we keep σ = 1 and change the magnitude of J 0 . Following the same scheme as in section III B, we break the rotational symmetry of the spin space by preparing the system at a real-valued configuration {ψ n }. When x = 0 the connectivity matrix has only random coupling elements (i.e J 0 = 0) "forcing" the network into the spin-glass phase. In this regime, the system evolves towards a metastable state with the "spins" pointing towards random directions. As a result, the magnetization is approaching zero x ≈ 1 becomes more abrupt as N increases. In these simulations, the initial preparation was taken to have averaged energy density h ≈ 0.88 1 (low energy regime) while the nonlinearity is χ = 0.01.
as |m| t ∝ O(1/N ) in the limit of large N -values, see Fig. 7. In the other limiting case of x → ∞ the randomness in the coupling elements are suppressed and the connectivity matrix is dominated by (essentially) constant couplings J 0 . In this case, the network is in the ferromagnetic phase and the magnetization acquires a non-zero magnitude |m| t = 0 which is dictated by the value of h (e.g. for h = h min ≈ −J 0 the magnetization is |m| t = 1). Our analysis (see Fig. 7) indicates that the transition from a spin-glass to a ferromagnetic phase occurs at x ∼ 1. The transition becomes more abrupt, as expected from statistical mechanics, in the limit of large N -values.

IV. CONCLUSIONS
We unveiled a connection between nonlinear photonic networks, consisting of many coupled single modes, and spin-networks. As opposed to standard spin models, our "photonic spins" are complex, soft (i.e. their size fluctuates), and their dynamics has two constants of motion: the total energy and the total optical power. This second conservation law is responsible for the appearance of novel optical phase transitions and the emergence of new forms of thermal photonic states. We have found that these transitions are controlled by the nature of the connectivity of the network (disorder or constant), and the amount of the averaged energy density of the initial optical preparation. Another important parameter is the strength of the non-linearity that controls the thermalization process. For strong non-linearities and constant couplings, we have discovered a ferro-paramagnetic phase transition as the averaged energy density h of an initial preparation increases. This transition is associated with the destruction of a photonic condensate and its depletion to thermal states. In the other limiting case of random coupling constants the photonic network is driven from a paramagnetic to a spin-glass phase as h decreases. Finally, we have shown that the same network, when prepared at low energies, undergoes another transition from a spin-glass to a ferromagnetic phase. The control parameter that drives this transition is the degree of randomness (frustration) of the coupling constants between the photonic spins. Our results shed light on the ongoing effort of taming the flow of electromagnetic radiation in nonlinear multimode photonic networks. We also expect to sparkle the interest of the statistical physics community since the mathematical models that are typically used to describe light transport in multimode systems are falling outside the traditional spin-network framework.