Entropic thermodynamics of nonlinear photonic chain networks

The complex nonlinear behaviors of heavily multimode lightwave structures have been recently the focus of considerable attention. Here we develop an optical thermodynamic approach capable of describing the thermalization dynamics in large scale nonlinear photonic chain networks - a problem that has remained unresolved so far. A Sackur-Tetrode equation is obtained that explicitly provides the entropy of these systems in terms of all the extensive variables involved. This entropic description leads to a self-consistent set of equations of state from where the optical temperature and chemical potential of the photon gas can be uniquely determined. Processes like isentropic expansion/compression, Joule expansion, as well as aspects associated with beam cleaning/cooling and thermal conduction effects in such chain networks are discussed. Our results open new vistas through which one can describe in an effortless manner the exceedingly complex dynamics of highly multimoded nonlinear bosonic systems that are nowadays of importance in areas ranging from high-power photonic settings to cold atom condensates and optomechanics.

In recent years, there has been considerable interest in understanding the rich and complex dynamics of large-scale nonlinear multimoded optical systems [1][2][3][4][5] . To a great extent, what has fueled this effort is the ever increasing quest for high information capacity photonic networks [6][7][8][9][10] and high-power optical sources 11 . Similarly, the nonlinear physics of high-Q multimode microresonators has been intensively investigated in a number of areas, like comb-generation and optomechanics [12][13][14][15][16] . Yet, what kept us from fully expanding these classical settings, is the sheer complexity behind the nonlinear interactions taking place in such vastly many-moded environments. To appreciate this difficulty, consider for example a multimode optical fiber supporting thousands of modes. In this case, to correctly account for all the nonlinear processes unfolding in this structure, one has to first estimate at least a trillion or so four-wave mixing coefficients, before even attempting to numerically solve the evolution problem -a formidable task which is by itself computationally intensive. Clearly, to overcome these hurdles, an altogether new approach has to be deployed that can effectively deal with "many-body" configurations. In essence, this calls for the development of a theoretical framework, akin to that of statistical mechanics, that has so far allowed us to understand the macroscopic properties of matter, in spite of the fact that on most occasions one has to deal with a great multitude of atoms or moleculesmost often exceeding Avogadro's number 17 .
Here, by means of an optical Sackur-Tetrode equation, we investigate the entropic response of large-scale weakly nonlinear photonic chain networks. These classical networks may for example involve a large number of evanescently coupled waveguide lattices 18 that exchange power in both crystalline (Fig. 1a) and other more complex arrangements (Fig. 1b), or could be comprised of arrays of optical resonators (CROWs) 19 that allow for energy transfer in time -used as building blocks as shown in Fig. 1c. While in this work, we will be primarily dealing with coupled waveguide topologies in order to exemplify these notions, a similar discussion holds for heavily multimode chain cavity configurations and other bosonic systems 20 . In this case, the underlying Sackur-Tetrode equation can be expressed as a function of the electrodynamic momentum flow , the number of available modes , and the input optical power . This analytical formalism developed here is in turn used to study and predict several distinct phenomena that could be of practical importance. These include for example, the prospect for beam cooling (leading to an improvement in beam quality) as well as Joule expansion and aspects associated with the photocapacity and thermal conductivity of such chain networks. We begin our discussion by considering a lossless, nonlinear array network of singlemode elements that are coupled via nearest-neighbor interactions. In this nonlinear optical chain, the equation of motion are given by represents the local mode amplitude at site and stands for the coupling coefficient between two successive elements. These evolution equations can be obtained from the classical The linear eigenvalues corresponding to the optical supermodes | � of the arrangement shown in Fig. 1a are given by = 2 cos[ /( + 1)] (ref. 18,19 ). In general, the evolution of this weakly nonlinear system can be described in terms of its supermodes, i.e., where the square of the moduli ) of the complex coefficients denotes modal occupancies. In this case, the norm in this is preserved, thus implying power/photon conservation in these coupled waveguide/cavity arrangements. In addition, in the weakly nonlinear regime, the expectation value of the total Hamiltonian 〈 〉 = 〈 〉 + 〈 〉 is dominated by its linear component 〈 〉 and as a result the internal energy remains invariant during evolution, thus establishing the second conservation law associated with this class of structures. In fact, in multimode waveguide settings, the conserved internal energy so happens to be the longitudinal electrodynamic momentum flow in these systems 21 . In this respect, the two normalized constants of the motion and are uniquely determined by the initial excitation . At this point, it is important to stress that the sole role of the weak nonlinearity involved ( (2) , (3) , etc.) is to chaotically reshuffle the optical power among modes through multi-wave mixing processes 22 -thus causing the complex amplitudes to vary randomly during propagation. This in turn, allows the system to ergodically explore in a fair manner, all its accessible microstates (in its phase space) that lie on the constant energy ( ) and power ( ) manifolds. It is this ergodicity that provides the foundation for establishing an entropic thermodynamic theory for the aforementioned multimode chain networks. Moreover, in heavily multimode arrangements, the thermodynamically extensive variables ( , , ) are related to each other through the optical temperature and chemical potential associated with the system. This relation is given through a global equation of state, given by − = that explicitly involves the total number of modes 23,24 . In general, any system is expected to reach thermal equilibrium by maximizing its entropy. In this case, the thermalized average power levels conveyed by each mode are found to obey a Rayleigh-Jeans distribution 25-28 , i.e. � � 2 = − /( + ).
Based on the above premises, one can show that the relative entropy associated with a photonic "monoatomic" chain network ( Fig. 1) is given by the following equation 24 , i.e., Equation (1) is analogous to the entropic Sackur-Tetrode equation -developed more than a century ago (based on an appropriate quantization of the phase space) in order to correctly account for the properties of monoatomic gases 17,29 . Like its counterpart in statistical mechanics, the optical entropy of these nonlinear monoatomic chains is extensive with respect to the three other extensive variables ( , , ). In other words, if we let ( , , ) → ( , , ), then from Eq. (1), one directly obtains → . We note that this extensivity in entropy is crucial in developing a selfconsistent thermodynamic formulation -free of Gibbs paradoxes. On the other hand, the optical entropy given above is by nature different from that of an ideal monoatomic gas in standard thermodynamics. This is apparent, since the prefactor in Eq. (1) is now given by the number of modes ("optical volume") as opposed to the optical power which here plays the role ofthe number of gas particles. Figure 2a depicts the entropy of a chain system as a function of ( , ) as obtained from Eq. (1), when = 100. The three pertinent equations of state can in turn be derived from (1) by employing the fundamental equation of thermodynamics 24 , i.e., Equations (2), relate the three intensive thermodynamic variables ( , ,) to their respective conjugate quantities ( , , ) . As in statistical mechanics, in this nonlinear multimoded setting, the optical temperature = ( 2 − 4 2 2 )/(2 ) represents a thermodynamic force responsible for energy transfer Δ from a hot to a colder subsystem whereas the chemical potential dictates the power exchange Δ -all aiming to maximize the entropy in accord with the second law of thermodynamics. Finally, ̂ represents the relative optical pressure in this configuration. From Eq. (2a), one can readily deduce that as the system attains its lowest possible internal energy → −2 (as imposed by the band structure of the monoatomic lattice), its temperature tends to zero ( → 0 + ) and as a result the photon gas becomes a condensate that exclusively occupies the ground state of the lattice. Conversely, when the internal energy is maximized → 2 , the opposite is true (i.e., → 0 − ), and hence the highest order mode is now solely occupied (anti-condensate) 30 . The first equation of state (Eq. (2a)) directly indicates that the temperature of these lattices is positive for −2 < < 0 whereas is negative if the internal energy lies in the domain 0 < < 2 . On the other hand, when → 0 ∓ , the temperature tends to → ±∞ (as schematically shown in the inset of Fig 2b), and as a result, equipartition of power takes place among modes, i.e., � � 2 = − / = / . In general, negative temperatures correspond to bodies "hotter" than hot 17 -a relation that regulates the direction of the energy flow ∆ . While for positive temperatures the lower group of modes is mostly occupied (since � � 2 = − /( + )), the opposite is true for negative temperatures where the higher-order modes are populated (Fig 2b). We note that, the set of the three equations in (2)  at various power levels, that as we will see dictate the photo-capacity of these photonic chain networks.
The above results can be used to predict the outcome of more complex processes like, for example, that associated with Joule expansion of the photon gas in such nonlinear heavily multimoded environments. This prospect is shown schematically in Fig. 3, where as an example, . What facilitates this "beam self-cleaning" effect has to do with the fact that during thermalization, the power is reshuffled in such a way so as to favor the lower order modes (for > 0) in the Rayleigh-Jeans distribution (Fig.4b) -a process that removes the initial speckle in the beam. Yet, while the output beam at ( = ) seems to improve in the nearfield, its multimode far-field is still highly divergent and thus is of poor quality as shown in Fig.   4d. Interestingly, the thermodynamic formulation developed here suggests that this situation can be greatly improved if in this same array structure, a "cooler" beam is launched -having a perpendicular polarization � (Fig. 4e). In this latter arrangement, each waveguide site is assumed to be substantially birefringent, thus prohibiting any power exchange between the two beams.
Instead, the two wavefronts only interact through cross-phase modulation 23 . In the example provided in Fig. 4e, the internal energy and power corresponding to each polarization is = −3.95, = 2, = −0.77, = 2. Based on these initial conditions, Eqs. 2 indicate that on its own, each polarization would have settled to an equilibrium temperature of = 1.7 × 10 −3 , = 0.33. On the other hand, once interacting together, the two beams reach the same temperature = 0.075 which can be exactly predicted following the procedure of ref. 24 . In other words, the � polarized beam can be considerably cooled after exchanging energy Δ with its � polarized counterpart. This in turn leads to a ~3-fold improvement in the far-field of the original beam (Fig.   4f).
The formalism developed here, allows one to define a photo-capacity = / associated with such nonlinear multimode photonic systems -a property that is in every respect analogous to that of the heat capacity pertaining to various phases of matter. Like its counterpart (heat capacity), this quantity critically depends on the density of states and hence it is characteristic of the multimode optical system itself. In this respect, it is formally possible to assign a capacity to any possible many-mode arrangement, irrespective of whether it is discrete or continuous. For the chain networks considered here, the photo-capacity can be explicitly obtained from Eq. (2a) and it is given by = − | | 2 ( 2 2 + 4 2 2 ) −1/2 . Figure 5a shows the function (blue curve) for such a chain configuration, involving = 200 sites while conveying a power of = 20. In this case, the photo-capacity function is an even function of temperature -as one should expect from the symmetry in the density of states. The photo-capacity curves corresponding to two different types of Lieb lattices are also displayed, as an example, in Fig. 5a for the same parameters ( = 200, = 20). Note that the ( ) function for the Lieb-2 lattice with cross-couplings is in this case highly asymmetric -a feature attributed to the band structure of this specific system.
We next explore the possibility for energy Δ transport phenomena occurring when a hot and a cold optical nonlinear multimode system are linked together under non-equilibrium conditions. To some extent, these processes have much in common with heat conduction problems in solids. Figure 5b shows such an arrangement where a square multicore photonic lattice 1 with 1 = 900 sites is connected to a graphene-like system 2 ( 2 = 896) via a chain bridge 3 involving 3 = 140 evanescently coupled waveguides. The bridge 3 is appropriately designed using birefringent elements 24 so as to allow for energy exchange Δ between the two lattices 1 and 2 while prohibiting any power transfer Δ . System 1 together with the link array 3 are initially brought to thermal equilibrium where both share a temperature of 1 = 0.3 when conveying a total power of 1+3 = 65 at a total internal energy 1+3 = −60. Similarly, lattice 2 (before is connected to the bridge) is kept at 2 = 0.1 when 2 = 40 and 2 = −46. Once 2 is put in contact to the 3 bridge, a non-equilibrium thermodynamical process ensues during which the entropy of the combined system starts to increase as expected from the second law of thermodynamics (Fig. 5c). During this non-equilibrium stage, internal energy Δ starts to flow through the bridge from hot to cold ( 1 → 2 ). In this case, the local temperature along this bridge (obtained after projecting on the local eigenfunctions) is found to linearly decrease from 1 to 2 in full accord to Fick's law of heat diffusion 35 , ⃗ = − ∇ . This in turn allows one to express this energy transfer through an effective "photo-conduction coefficient" as Δ /Δ = − ( 1 − 2 ).
From our simulations, we estimate that here the photo-conduction coefficient of the 3 chain is ≈ 3.2 × 10 −5 .
By harnessing concepts from statistical mechanics, we have developed an optical thermodynamic theory that explicitly provides the Sackur-Tetrode entropy associated with large nonlinear photonic chain networks. Based on these premises, the temperature and chemical potential of such systems can be obtained in closed form as a function of the initial excitation conditions. The formalism developed here is general and can be readily deployed to describe a number of thermodynamic processes, like isentropic expansion/compression, Joule expansion, as well as aspects pertaining to beam self-cleaning and cooling in such chain networks. Our results not only provide a platform to understand and predict in a systematic way the utterly complex dynamics of heavily multimoded nonlinear optical systems but could also be of practical importance in designing new classes of high-brightness multimode optical sources.  depicts the modal occupancies corresponding to different temperatures. While for positive temperatures the lowest group of modes is occupied, the converse is true for negative temperatures. When → ±∞, equipartition of power takes place among modes. c, Energy-temperature diagrams pertaining to this family of lattices corresponding to the same power levels in (a and b), shown as colored curves. In all cases, the number of modes was taken to be = 100.

Fig. 3 | Joule expansion in a nonlinear photonic multimoded chain system.
Light conveying power is allowed to suddenly expand in a larger chain network having four times as many waveguide elements ( → 4 ). In this case, if the temperature of the photon gas before this transition is , after this abrupt and thermodynamically irreversible expansion, the system reaches equilibrium at /4 -a direct consequence of the way Joule expansion manifests itself in these heavily multimoded nonlinear optical environments. Fig. 4 | Thermodynamic beam cleaning and cooling. a, Propagation of a � polarized optical beam in a nonlinear chain of waveguides involving 30 elements. The projected far-field after the beam exits the system at = is also shown. b, Thermalization dynamics of the modal occupancies as a function of distance when = 2 and the modes are uniformly excited in the range of −0.5 ≤ ≤ 1.2 . At the output, the system settles into a Rayleigh-Jeans distribution with a temperature of = 0.33. c, The intensity distribution corresponding to the input used in (b) displays a strong speckle (upper panel). The beam selfcleans after thermalization takes place (lower panel). d, Far-field pattern associated with the thermalized � polarized output field for the same conditions used in (b). e, Cooling the � polarized wavefront using crossphase modulation with an � polarized beam of equal power. f, After cooling, the far-field of the � polarized wave experiences considerable improvement. This cooling interaction can be analytically predicted 24 from Eqs. 2. This energy transfer takes place through a bridge chain composed of birefringent elements under nonequilibrium conditions. As required by the second law of thermodynamics, Δ always flows from a hot to a cold subsystem. In this case, the temperature in the linking chain drops linearly, in accord with Fick's law. c, The internal energy transfer ceases when thermal equilibrium is reached in which case both subsystems attain the same temperature and the total entropy is maximized. In all cases, the nonlinearity was assumed to be of the Kerr type.

Global equation of state
Using the expressions for the two conserved quantities associated with the optical power and internal energy , one can write A direct substitution of the Rayleigh-Jeans distribution, � � 2 = − /( + ) [24-26] into the above expression, leads to a global equation of state: − relation associated with a multimode network chain of nonlinear coupled elements Let us consider a 1D chain of weakly nonlinear elements. These elements are evanescently coupled via nearest neighbor interactions with an exchange strength . Starting from the Rayleigh-Jeans distribution, the power conservation law in a nonlinear multimode optical system can be written as: On the other hand, the global equation of state requires that Given that the eigenspectrum of a 1-D discrete lattice is given by = 2 cos � +1 �, we obtain . Note that one can rigorously show that for positive temperatures > 0 , > + 2 whereas for negative temperatures ( < 0) the following is true < − 2 . As a result, −1 < < 1 or 2 < 1. Since these structures support a large number of modes , the summation above can now be approximated with an integral.
After applying the limits tan � 2 � = +∞, tan −1 (+∞) = 2 , we find Because is large, + 1 → , and given that = The last three expressions in this section provide the − relation associated with a multimode network chain of nonlinear coupled elements.

Optical Sackur-Tetrode equation
From the fundamental equation of thermodynamics, the temperature can be obtained from the entropy according to Finally the chemical potential in every subsystem can be extracted from the global equation of state ( − = ). From here, the Rayleigh-Jeans distribution in each subsystem follows.

Normalizations used in multimode nonlinear lattice systems
In a 1-D discrete nonlinear photonic arrangement, the evolution of the optical modal complex amplitude at the ℎ site is governed by: where is the propagation constant at the corresponding site, is the nearest-neighbor hopping coefficient, 0 is the wavenumber in vacuum and 2 represents the effective nonlinear Kerr coefficient. When the photonic lattice is truncated to have sites, the boundary condition is 0 = +1 = 0. By assuming that all sites are identical ( = 0 ), and by setting = 0 , = / 0 , = − 0 / , = 2 0 2 / 0 = 1, one obtains the normalized evolution equations: The above equation of motion (in the local mode representation) can be derived from the system total Hamiltonian that is composed by a linear and a nonlinear part: The linear part of Hamiltonian can be readily written in the supermode representation by into Eq. (S7), i.e., The arrays are supposed to be lossless and having no leakage to radiation modes. The narrow band structure formed by the fundamental mode is expected to lie away from the radiation modes of the system.

Lieb lattice
Each unit cell in the Lieb lattice contains three sites, indicated by red, green and blue. When only the vertical and horizontal nearest neighbour couplings ( 1 = 1) are taken into account, the Lieb lattice (Lieb-1) displays a highly degenerate flat band at the centre of the eigen-spectrum. When the cross-coupling ( 2 = 0.5) is present, the band structure of the Lieb lattice becomes asymmetric (Lieb-2).
FIG. S1 | Geometric structure of a Lieb lattice.

Hexagonal, chain and square lattices
The exchange coupling mechanisms between nearest neighbors in discrete hexagonal, chain and square lattices are presented schematically in Fig. S2.

Governing nonlinear equations for the two orthogonal polarizations
In a chain network, the discrete coupled nonlinear Schrödinger equations describing the evolution of the two orthogonal polarizations and associated with the normalized optical field at site are given by: Here denotes the site number, ( ) the coupling coefficients, , the propagation constants of these waveguides, and the last three terms correspond to self-phase modulation, cross-phase modulation and four-wave mixing effects. The nonlinear coefficients , and result from the (3) tensor. If the waveguide birefringence ( − ) in every element is strong with respect to the nonlinearly induced index change, the four-wave mixing process can be ignored ( = 0) in the above equations. Throughout this study we have used = 1, = 2/3 and = 0 ., values corresponding to silica glass. In the same vein, the governing equations for the two optical polarizations can be readily written for hexagonal and square lattices.

Thermal conductivity simulations
The simulations associated with thermal conductivity were performed for an arrangement where a photonic square lattice 1 with 1 sites is connected to a graphene-like system 2 with 2 sites via a chain bridge 3 involving 3 evanescently coupled waveguides. The chain bridge consists of 7 sub-chains ( 3 (1) , 3 (2) , … , 3 (7) ) with 3 ( ) number of sites (see Fig. S3). Each optical array conveys a different polarization (blue and red sides represents � and � polarizations). On the other hand, the polarizations are allowed to overlap in the thermal bridge layer (purple) in between. This bridge "diathermal" layer allows only for energy transfer from one to the other optical array via cross phase modulation. In such arrangements the field evolution is governed by Eqs. S8-S9.