Thermalization dynamics in photonic lattices of different geometries

The statistical mechanical behavior of weakly nonlinear multimoded optical settings is attracting increased interest during the last few years. The main purpose of this work is to numerically investigate the main factors that affect the thermalization process in photonic lattices. In particular, we find that lattices with identically selected properties (such as temperature, coupling coefficient, lattice size, and excitation conditions) can exhibit very different thermalization dynamics and thus thermalization distances. Our investigation is focused on two different two-dimensional lattices: the honeycomb lattice and the triangular lattice. Our numerical results show that, independently of the excitation conditions, the honeycomb lattice always thermalizes faster than the triangular lattice. We mainly explain this behavior to the quasilinear spectrum that promotes wave-mixing in the honeycomb lattice in comparison to the power-like spectrum of the triangular lattice. In addition, we investigate the combined effects of temperature as well as the sign and magnitude of the nonlinearity. Switching either the sign of the Kerr nonlinear coefficient or the sign of the temperature can lead to significant differences in the thermalization dynamics, a phenomenon that can be physically explained in terms of wave instabilities. Larger absolute values of the temperature |T| result in more uniform distributions for the power occupation numbers and faster thermalization speeds. Finally, as expected, increasing the magnitude of the nonlinearity results in accelerated thermalization. Our findings provide valuable insights into optical thermalization in discrete systems where experimental realization may bring about new possibilities for light manipulation and applications.


Introduction
Over the last few years, optical thermodynamics, the formulation of statistical mechanics for optical systems, has been established through extensive research efforts [1][2][3].This theory was first put forward to explain the peculiar phenomenon of beam self-cleaning [4][5][6].An intense light injected into a highly multi-mode optical parabolic graded-index fiber exhibits a bell-shaped output intensity pattern after nonlinear propagation.The beam self-cleaning phenomenon has been initially attributed to multi-wave mixing processes as an explanatory framework [7].However, recent studies have introduced novel theoretical approaches based on thermodynamics, offering a more concise method for addressing and describing it, leading to the experimental observation of Rayleigh-Jeans (R-J) distribution in graded-index multimode fibers [8][9][10].The emergence of the R-J distribution is directly linked to principles of statistical mechanics, in particular, the maximization of entropy within a multimode system with conserved quantities [1,[11][12][13].Remarkably, it has been demonstrated that thermal equilibrium manifests across multiform multimode optical systems that obey certain conservation laws (i.e., the total optical power, and the internal energy), each of them characterized by distinct nonlinearities even beyond the multi-wave mixing paradigm [14].In different optical systems, multiple theoretical studies have predicted various interesting thermal behaviors, such as entropic processes [15,16], optical pressure in a multimoded optical setting [17], and also thermodynamics in topological and non-Hermitian structures [18,19].Further advances reporting the thermal equilibrium state with the additional conserved quantity of optical angular momentum (OAM) in fiber have been also investigated theoretically and experimentally [20,21].Also, the entropic behaviors under negative temperatures in both fiber and photonic lattices have been experimentally investigated [22][23][24].
In this work, we analyze the thermalization process and the associated thermalization distances in multimoded weakly nonlinear photonic lattices.In doing so, we utilize the Kullback-Leibler divergence (KLD), a measure that is additive for independent distributions [25][26][27].In particular, we show that the thermalization processes between different lattice geometries can lead to thermalization distances that are significantly different.Specifically, here we select to compare honeycomb lattices with triangular lattices.To ensure a fair comparison, both lattices are chosen to have the same number of waveguides and the same nearest neighbor coupling coefficients.Furthermore, in both lattices, we excite the same number of consecutive modes with the same amount of power using a uniform distribution.These sequences of modes are selected to lead to approximately the same value of the optical temperature.We observe that even though both lattices have the same parameter set, there is a large difference in the thermalization dynamics with the honeycomb lattice being significantly faster.We attribute this difference to the spectrum of the honeycomb lattice that has a quasi-linear structure and thus promotes wave mixing phenomena.In addition, we examine the influence of different magnitudes and signs of optical temperature and nonlinearity on the thermalization speed of these two lattice geometries.Such a temperature depends on the spectral feature of lattice geometry and the excitation conditions, in particular power and internal energy, but is not related to the nonlinear strength.In general, larger absolute values of the temperature lead to more uniform distributions of the power occupation numbers and faster thermalization dynamics.Thus, it is significantly more difficult to observe thermalization close to the condensation limits || → 0. In addition, wave instabilities affect the thermalization process leading to different thermalization distances, depending on the sign of  • , where  is the Kerr nonlinear coefficient.Finally, as expected, the thermalization speed is increased by increasing the magnitude of the Kerr nonlinear coefficient.Our findings not only provide valuable insights into the intricate details of the thermalization process with a concise method but also offer practical guidance for the experimental realization of thermodynamic principles in lattice-based systems.

Thermodynamics of photonic lattices
In the presence of Kerr nonlinearity, paraxial light propagation in a photonic lattice structure is governed by the following nonlinear Schrödinger type equation where (, , ) is the electric field envelope,  and y are the transversal coordinates and  the propagation distance, (, ) is the lattice potential that is proportional to the refractive index,  is the Kerr nonlinear coefficient, and  is the wavenumber.In particular, using coupled-mode theory, if next-nearest-neighbour couplings are negligible with respect to the coupling between the nearest neighbors, Eq. ( 1) can be simplified into the following discrete nonlinear Schrödinger equation where, in Eq. ( 2),   is the complex field amplitude at the th lattice site, Z represents the dimensionless propagation distance, and  is the coupling coefficient between nearest neighbors ⟨, ⟩.We expand the wavefunction |⟩ = ∑   |  ⟩  into the eigenfunctions (or supermodes) |  ⟩ of the system, which satisfy where  is the linear part of the Hamiltonian,   is the associated propagation constant of the th mode, and, since each waveguide in isolation is a single mode, the total number of nodes or lattice sites  is equal to the total number of modes.The power occupation number on each super mode is and the total power  = ⟨|⟩ = ∑     is considered to be the first conservation law with the second being the internal energy approximately equal to its linear part  = ⟨||⟩ = ∑       .Following the calculations, the optical entropy is found to satisfy [1]  = ∑ ln    =1 + . ( In the case that the lattice is relatively large along both directions, the entropy  = (, , ) can be approximated to be extensive with respect to (, , ) [1,16].The ensemble average value of the entropy increases in a spontaneous manner with the propagation distance (d⟨⟩/d ≥ 0) until asymptotically approaching its maximum value.In thermal equilibrium, the ensemble average power occupation on the eigenmodes obeys the Rayleigh-Jeans (R-J) distribution where  = −/ and  = 1/ denote the Lagrange multipliers of the constraints imposed by the conservation laws.This process of reaching thermal equilibrium is called thermalization, and  and  are defined as the optical temperature and chemical potential, respectively, being the conjugate variables to the internal energy and total power given by  −1 = / and  = −/.The extensive variables , ,  are related to ,  through the following equation of the state [1]  −  = .
For a given lattice system with eigenspectrum   , which is excited with power  and energy , the optical temperature is uniquely determined by the expression while the chemical potential is then obtained from Eq. ( 7).

Thermalization dynamics
In the rest of this paper, we examine and compare the thermalization process of two photonic lattices under different conditions of the optical temperature and Kerr nonlinearity.Specifically, we focus our attention on investigating the thermalization processes associated with the honeycomb and the triangular lattices [see Figs.1(a1,a2)].Both of them belong to the same hexagonal Bravais group but have different features.In particular, the honeycomb lattice is a compound lattice structure with two atoms per unit cell, largely studied in optics due to a host of interesting topological phenomena it is associated with.Interestingly, it possesses a pair of energy bands that touch each other forming a Dirac-like cone, and such a singular band touching is demonstrated to induce the property of pseudospin around the Dirac points [28,29].On the other hand, the triangular lattice is one of the simplest two-dimensional lattice structures, with the unit cell containing a single atom.Figure 1(a1,a2) illustrates the schematic representations of the finite-sized honeycomb and triangular lattices examined here, where blue and red circles mark their respective lattice-site locations, while black lines connecting adjacent sites show the direction of the nearest neighbor couplings.The number of lattice sites for both lattices is set to  = 144, the nearest neighbor hopping amplitudes  are taken to be 1, and zero boundary conditions are assumed throughout the analysis in order to consider a realistic scenario and account for boundary effects.Under these conditions, the band structures are directly calculated by solving the eigenvalue problems of Eq. ( 3) and they are plotted in Fig. 1(b1,b2) as a function of the mode number.One can see that the energy spectrum associated with the finite-sized honeycomb lattice displays a quasi-straight-line profile as the mode number increases [Fig.1b1], with the band-structure profile being symmetric with respect to zero energy level due to the chiral symmetry, which is preserved in this model.On the other hand, the band structure of the triangular counterpart [Fig.1b2] exhibits a power-law-like relation with an exponent greater than that one for increasing mode-number values.Indeed, the band-structure profile of the triangular model is asymmetric with respect to the zero-energy level due to the lack of chiral symmetry [30].In addition, numerical results in Fig. 1(c1,c2) compare the ensemble average thermalization process when light propagates in the two lattice structures in the presence of self-focusing nonlinearity.The initial conditions consist of a narrow band of eigenmodes that are uniformly excited (see yellow curves in Fig. 1(c1,c2)).The numerical simulations are carried out by solving Eq. ( 2) using a Runge-Kutta numerical method.In particular, the total optical power is selected to be  = 12.9 , and it is uniformly distributed among 21 consecutive eigenmodes, with a spectral bandwidth of ∆  = 0.01 − (−0.91) = 0.92 for the honeycomb lattice and ∆  = −0.24− (−1.16) = 0.92 for the triangular lattice (see yellow lines in Fig. 1(c1,c2)) so that the resulting optical temperature at thermal equilibrium is the same for both lattices:  = 0.5.The power occupation numbers at each propagation distance  are averaged over 200 realizations.As shown in Fig. 1(c1,c2), during propagation, the ensemble average power occupations progressively evolve from the initial state towards the respective theoretically predicted R-J distributions (which are shown with the red curves in Fig. 1(c1,c2)).After a propagation distance  = 1000 both systems are visually observed to reach their thermal equilibrium state.Looking at the ensemble average entropy curves, which are calculated from Eq. ( 5), as a function of the propagation distance , they increase and asymptotically approach the theoretically predicted thermal equilibrium value, as shown in Fig. 1(d1,d2).Furthermore, shaded gray areas enclosed between the entropy curves and the vertical dashed lines mark the locations where entropy reaches 98% of its maximum at the thermal equilibrium.Remarkably, they highlight the difference in the rates of changes of the entropies for both systems, thus indicating qualitatively that the thermalization process in the honeycomb lattice develops faster than the one occurring in the triangular system.In summary, our numerical results of Fig. 1 show that the honeycomb lattice tends to reach thermal equilibrium at shorter distances than the triangular lattice in comparable lattice and excitation conditions.A physical explanation can be found in the fact that the thermalization process is the outcome of different supermodes experiencing a nonlinear four-wave-mixing process, the latter of which is initiated by the presence of a small Kerr nonlinearity.The cause of different thermalization speeds in these two systems is attributable to the energy-spectral profile of their lattice geometries.In particular, honeycomb lattice features a quasi-straight-line profile of its energy spectrum or better a quasiconstant energy spacing, thereby optimizing the phase-matching conditions  1 +  2 =  3 +  4 and promoting the nonlinear wave-mixing phenomenon when compared with the triangular counterpart.

Kullback-Leibler Divergence (KLD)
In analyzing the thermalization dynamics, it is important to find a measure that directly compares the numerical ensemble average power occupation numbers with the theoretical predictions.In principle, we could use different normed metrics to measure the distance between these two distributions.However, we prefer to use the Kullback-Leibler divergence (KLD) (also known as the relative entropy), which is not formally a metric since it does not satisfy the symmetry relation (, ) = (, ) or the triangle inequality.In particular, the relative entropy of the ensemble average power occupation numbers ⟨  ⟩ at distance , from the theoretically predicted thermal equilibrium distribution    is defined as The KLD divergence is a useful measure in statistical mechanics that defines the distance between two distributions.It is a direct outcome of the Gibbs inequality, and as a result, it can only take non-negative values   (|) ≥ 0 and becomes zero   (|) = 0 if and only if  = .A property of the KLD, important in statistical mechanics, is that it is additive for independent distributions.An illustrative example where the   is plotted as a function of the propagation distance for the honeycomb lattice with  = 0.5 and  = 1 is presented in Fig. 2. One can observe that the KLD has a decreasing trend as a function of the propagation distance, even though, especially for large propagation distances, it exhibits fluctuations.We have selected to set the thermalization threshold to the value   = 0.05, which is supported by our numerical results.For example, in Fig. 2(b1)-(b4) we show the trend for the system to reach thermal equilibrium as the propagation distance increases.In particular, in Fig. 2(b1) we see the initial spectrum.At  = 300, the system is still far from equilibrium [Fig.2(b2)].Although the higher eigenvalues are approaching the theoretically predicted distribution, the first eigenstates exhibit large deviations.In Fig. 2(b3) we show the power occupation numbers at the selected thermalization onset where   = 0.05.We actually observe very good agreement between theory and numerical results.Finally, at  = 900 (see Fig. 2(b4)) we have an even better agreement between the two distributions.For the purpose of confirming the high degree of similarity between the theoretical R-J distribution and the numerically calculated ensemble average power occupation numbers observed when   = 0.05, we have statistically checked via a two-sample Kolmogorov-Smirnov test [31] that the confidence coefficient between the two curves is 0.99.Accordingly, we conclude that due to the high degree of similarity between the R-J and numerical distributions, the value   = 0.05 can be considered as a thermal equilibrium threshold.

Thermalization dynamics in honeycomb and triangular lattices
Using as a main measure the mutual entropy, we will investigate the thermalization dynamics associated with photonic lattices under different temperature conditions and signs of the Kerr nonlinearity.As discussed in the previous section, it is reasonable to define the thermalization distance as the distance at which the mutual entropy   is equal to 0.05.First, we would like to examine whether different lattices under the same conditions have similar thermalization distances or not.In this respect, for simplicity, we are examining the two different types of lattice structures that are shown in Fig. 1(a): the honeycomb and the triangular lattices.To equalize their properties, both lattices have  = 144 nodes, coupling coefficient between first neighbors  = 1, and magnitude of the nonlinearity equal to || = 1.In each configuration, the initial power is  = 12.9, which is uniformly distributed in a band of consecutive eigenmodes.This band is selected so that the predicted value of the temperature under thermal equilibrium is almost the same for both types of lattices.In Fig. 3(a), we present the thermalization distances for both triangular and honeycomb lattices, computed for five different positive and negative temperatures in the presence of self-focusing and self-defocusing nonlinearities (|| = 1).By comparing the thermalization speeds between the honeycomb and the triangular lattice, we see that under the same conditions (at least those presented in Fig. 3) the honeycomb lattice is always significantly faster than the triangular (note that the scale in the figure is logarithmic).We render this behavior to the quasi-linear type of spectrum of the honeycomb lattice that promotes power exchange between the modes.It is apparent from Fig. 3(a) that, independently of the sign and value of the optical temperature and the sign of nonlinearity, the thermalization process is accelerated as the absolute value of the optical temperature || increases.Thus, it is significantly more difficult to achieve thermalization close to the condensation limits || → 0. On the other hand, the configurations that are easier to achieve thermalization are associated with large absolute values of the temperature ||, where the Rayleigh-Jeans distribution tends to become more uniform.The honeycomb lattice has a chiral symmetry and, thus, for each eigenmode |  ⟩ with eigenvalue   , there exists an eigenmode |  ⟩ with eigenvalue −  , where  is the chiral operator that changes the sign of one sublattice.The initial excitation conditions for opposite values of the temperature  → − are related through the application of the chiral operator |⟩ → |⟩.As a result, reversing both the temperature and the sign of the nonlinearity  → − leaves the discrete nonlinear Schrödinger equation invariant.Thus, the ensemble average results obtained via the change of variable (, ) → (−, −) should actually be identical.This is what we observe in Fig. 3(a), where the small fluctuations are the outcome of using a finite ensemble set.Note that such a symmetry is not present in the triangular lattice.As a result, the thermalization curves for positive and negative temperatures are different.In the case of the honeycomb lattice, it is interesting to observe that for the same absolute values of the temperature and nonlinearity, systems that satisfy the condition  •  < 0 will always thermalize faster.We can qualitatively understand this behavior in terms of modal instabilities.For example, in the case of positive temperature  > 0, the thermal equilibrium condition as described via the Rayleigh-Jeans distribution has larger power occupation numbers for the lower order modes.However, for  > 0, the lower order modes (which are more "in-phase") are associated with instabilities and, thus, it is far more difficult to excite them in comparison to the  < 0 regime where the lower order modes are more stable.The triangular lattice does not possess a chiral symmetry and its spectrum is power-like.We see that, on average, we achieve faster thermalization when the condition  •  < 0 is satisfied.However, approaching the condensation limit  → 0 − leads to far slower thermalization speeds as compared to the limit  → 0 + .This can be physically explained by the form of the spectrum.Due to the power-law relation, it is significantly more difficult to excite the higher-order modes, which have larger spacing between them.
Finally, in Fig. 3(b1,b2) we depict the dynamics of the relative entropy for light propagating in both honeycomb and triangular lattices under self-focusing and self-defocusing Kerr nonlinearities (|| = 1), for both positive and negative temperatures with || = 1.One can clearly see that under the same conditions, the   curves of the honeycomb lattice always cross the   = 0.05 threshold lines earlier than the triangular counterpart, indicating faster thermalization.We see that in the honeycomb lattice with  •  < 1 thermalization is significantly faster than the rest of the combinations.

Effect of nonlinearity in the thermalization process
In this section, we investigate the effect of the strength of the nonlinearity in the thermalization process.We select to keep the same initial excitation and vary only the sign and value of the Kerr nonlinear coefficient  .In our numerical simulations, the amount of nonlinearity is maintained small enough to avoid the formation of optical solitons during propagation, which will affect the thermalization process.Precisely, numerical simulations are performed for temperatures  = 0.5 by maintaining the same parameter set as in Fig. 3a, except for changing the strength of  .As seen in Fig. 4(a), the thermalization process in the honeycomb lattice is observed to occur faster as compared to the triangular lattice for any sign and strength of the nonlinear coefficient .The main outcome is that increasing the strength of the nonlinearity || leads to a faster thermalization process.This is an expected result since, by using a modal decomposition, the wave mixing coefficients directly depend on .In Fig. 4(b1, b2), typical   dynamics of the relative entropy are shown for  = 1.5 and −1.5, respectively.The   curve of the honeycomb lattice for  = −1.5 decreases faster and crosses the 0.05 line at an earlier distance as compared to the honeycomb lattice for the same temperature with  = 1.5.
At the end of this section, we would like briefly to mention some ways to achieve thermalization dynamics in realistic physical realizations of triangular and honeycomb lattice structures.The main requirements for this purpose are a larger number of available supermodes and waveguide length necessary to reach the thermal equilibrium.In photonic platforms, for example, lattice structures fabricated via femtosecond writing laser technique in glass or appropriate arrangements of coupled nonlinear single-mode fiber (SMF) could be employed.However, clear experimental observations of the thermalization phenomenon require an order of hundred coupling lengths.

Discussion and Conclusions
In conclusion, we have thoroughly investigated the thermalization dynamics occurring in honeycomb and triangular photonic lattices.Specifically, we have examined how the use of different lattice structures can affect the thermalization process, as well as the combined effects arising from the use of different temperatures, and signs and magnitudes of the Kerr nonlinearity.We would like to mention that the results of this work are focused in the case of honeycomb and triangular lattices.However, we expect that the results from this investigation may be applicable to other types of lattice structures, which certainly merit further investigation.However, we expect that our results will be particularly useful in this newly developed field of nonlinear optical thermodynamics, providing valuable insights into optical thermalization in discrete systems where experimental realization may bring about new possibilities for light manipulation and applications.

Fig. 1 .
Fig. 1. (a1) Schematic illustration of a finite-sized honeycomb lattice with coupling coefficient  = 1, formed by  = 144 lattice sites.(b1) Corresponding real-space propagation constants.(c1) Ensemble average dynamics of the power occupation numbers illustrating the thermalization process occurring for Kerr-type nonlinearity with  = 1 under uniform spectrum excitation of a band of supermodes (yellow line), and reaching thermal equilibrium as predicted by the Rayleigh-Jeans distribution (red line).(d1) Ensemble average value of the optical entropy (black curve) asymptotically approaching the predicted theoretical maximum value (red line), where the vertical dashed line marks the propagation distance where 98% of the maximum is reached.(a2-d2) Same as in (a1-d1) for a finite-sized triangular lattice with 144 nodes.

Fig. 2 .
Fig. 2. (a) Relative entropy in a honeycomb lattice as a function of propagation distance Z. (b1-b4) Power occupations at selected Z under uniform initial excitation (yellow line) of the power  = 12.9 over 21 consecutive modes, compared with theoretical R-J distribution (red lines).The predicted value of the temperature is  = 0.5: (b1)  = 0; (b2)  = 300; (b3)  = 628; and (b4)  = 900.Power occupation curves approach the R-J distribution as   decreases below the 0.05 threshold value (dotted black line) at b3, displaying a confidence coefficient of 0.99 when estimated by a two-sample Kolmogorov-Smirnov test.

Fig. 3 .
Fig. 3. (a) Thermalization distances  (in logarithmic scale) at different optical temperatures T for self-focusing (solid) or self-defocusing (dotted) Kerr nonlinearity with || = 1 in triangular (blue curves) and honeycomb (red curves) lattices with  = 144 modes and coupling coefficient  = 1.The initial excitation consists of 21 consecutive supermodes that are selected to match the same value of the optical temperature with total power  = 12.9 .(b1, b2) Dynamics of the K-L divergence for (b1)  = 1 and (b2)  = −1, respectively.Dotted black lines in (b1) and (b2) mark the threshold value at 0.05.