2.1. Dynamics in local base

We begin the discussion by considering a general framework upon which the current work is based. Specifically, the physical model under study is a finite size one-dimensional generalized DNLS type system subject to zero boundary conditions that includes a nearest-neighbor coupling, on-site potential, and a collection of cubic nonlinearities given by

(2.1)\begin{equation} i\frac{dA_n}{dz} + A_{n-1}+A_{n+1} + V_nA_n + f(A_n,A_{n\pm m}) = 0, \end{equation}

where $f(A_n,A_{n\pm \tilde n})$ is a complex cubic polynomial of the optical field amplitude $A_n$. Here, $n = 1,2,\dots, M$ indicates the lattice site, $z$ is the propagation distance, $V_n$ is the on-site potential and $\tilde n$ is an arbitrary integer. It is assumed that Eq. (2.1) remains invariant under the gauge transformation:

(2.2)\begin{equation} A_n \rightarrow A_ne^{i\theta}, \end{equation}

with real constant $\theta$. While at this point the specific structure of the cubic nonlinearities appearing in Eq. (2.1) is left arbitrary, the gauge symmetry (2.2) would eventually limit their form. In this regard, we have

(2.3)\begin{equation} f(A_n,A_{n\pm \tilde n}) \sim A_{n+n_1}A_{n+n_2}A^*_{n+n_3}, \end{equation}

where star indicates complex conjugation and $n_1$, $n_2$, $n_3$ are arbitrary integers. This means that lattices with nonlinearities of the form $|A_{n\pm n_1}||A_{n\pm n_2}|A_{n\pm n_3}$ (or their alike) are not included in our study.

2.2. Eigenvalue problem

In the linear regime (where all nonlinear cubic terms are absent) the ansatz

(2.4)\begin{equation} A_n(z) = \psi_n e^{i\epsilon z}, \end{equation}

leads to the following linear matrix eigenvalue problem:

(2.5)\begin{equation} \mathbb{M} | \psi^{(j)} \rangle = \epsilon_j | \psi^{(j)} \rangle, \end{equation}

where $\mathbb{M}$ is an $M\times M$ tridiagonal matrix whose main diagonal is the on-site potential $V_n$ and its two off diagonals are ones. The notation $| \psi^{(j)} \rangle,\,\, j = 1,2, \dots, M,$ is used to denote the $j^{\text{th}}$ supermode (also known as normal mode), defined by

(2.6)\begin{equation} | \psi^{(j)}\rangle \equiv \{\psi_n^{(j)}\}_{n=1}^M = \begin{pmatrix} \psi^{(j)}_{1}\\ \\ \psi^{(j)}_{2}\\ \vdots\\ \psi^{(j)}_{M} \end{pmatrix}, \end{equation}

and $\epsilon_j$ being the matrix eigenvalues. The supermodes constitute an orthonormal base satisfying the orthogonality condition:

(2.7)\begin{equation} \langle \psi^{(j)} | \psi^{(j')} \rangle =\delta_{j,j'}. \end{equation}

Without loss of generality, it is assumed that the eigenvalues are arranged in ascending order $\epsilon_1\leq\epsilon_2\leq \cdots \leq\epsilon_M$ where $\epsilon_1$ ( $\epsilon_M$) corresponds to the highest (lowest) order supermode. For the special case where $V_n=0$ and in the presence of zero boundary conditions, i.e. $\psi^{(j)}_0=\psi^{(j)}_{M+1}=0$, an analytical form for the supermodes $\psi_{n}^{(j)}$ and the eigenvalues $\epsilon_j$ of Eq.(2.5) can be obtained. They are given by

(2.8)\begin{equation} \epsilon_j=2 \cos\left(\frac{j\pi}{M+1}\right), \end{equation}

(2.9)\begin{equation} \psi_n^{(j)}=\sqrt{\frac{2}{M+1}}\sin\left(\frac{jn\pi}{M+1}\right). \end{equation}

2.3. General formulation in modal space

We seek solutions to Eq. (2.1) in the form of

(2.10)\begin{equation} A_n(z)=\sum\limits_{j=1}^M c_j(z)\psi_n^{(j)}, \end{equation}

where $c_j(z)$ are the so-called complex projection coefficients. Substituting the expansion (2.10) into (2.1) we obtain the following generic nonlinear evolution equation governing the dynamics of the projection coefficients (or modal amplitudes):

(2.11)\begin{equation} i\frac{dc_j}{dz}+\epsilon_j c_j+\sum_{k,l,m=1}^{M} T_{j,k,l,m}\,c_k c_l c^*_m=0, \end{equation}

for all $j=1,2,\dots, M$. In the above, $T_{j,k,l,m}$ is a rank 4 (in general) complex mixing tensor that arises due to the presence of cubic nonlinear terms. For the type of nonlinearity shown in Eq. (2.3) this tensor takes the form

(2.12)\begin{equation} T^{\text{typical}}_{j,k,l,m} =\sum_{n=1}^M\psi_n^{{(j)}^*}\psi^{(k)}_{n+n_1} \psi^{(l)}_{n+n_2} \psi^{{(m)}^*}_{n+n_3}. \end{equation}

Note that Eq. (2.11) is invariant under the phase transformation $c_j \rightarrow c_je^{i\vartheta}$ with real constant $\vartheta$. If one relaxes condition (2.2) then three other groups of nonlinear cubic terms would appear. They assume the form: $T_{j,k,l,m}\,c_k c^*_l c^*_m,\; T_{j,k,l,m}\,c^*_k c^*_l c^*_m$ and $T_{j,k,l,m}\,c_k c_l c_m.$ Such scenario is not considered in this paper and is left for a future study. System (2.11) constitute the basis of our work as it provides a unified framework for the study of the statistical properties of DNLS systems with the aforementioned type of cubic nonlinearities. When it comes to the study of equilibrium of many-body systems in the near-integrable regime, recent studies have proposed new means for diagnosing thermalization through the scaling properties of the Lyapunov spectrum [Reference Malishava and Flach33]. In our work, we adopt the more conventional approach by selecting a physically meaningful observable and follow its relaxation toward a prescribed equilibrium distribution. In this regard, the quantity of interest (observed quantity) is defined by

(2.13)\begin{equation} \langle |c_j|^2 \rangle\equiv \langle |c_j(Z)|^2 \rangle, \end{equation}

where $\langle \cdots \rangle$ denotes the ensemble average over many realizations of the random initial conditions $c_j(0)$ and $Z$ is the length of the simulation interval, i.e., $0\leq z \leq Z$. Of particular importance is the dependence of $\langle |c_j|^2 \rangle$ (at thermal equilibrium) on the eigenvalues $\epsilon_j$. We note that if the dynamical system admits a thermal equilibrium, then the associated quantity $\langle |c_j|^2 \rangle$ should be independent of the initial conditions $c_j(0)$. When this is the case, we have

(2.14)\begin{equation} \langle |c_j|^2 \rangle=h(\epsilon_j), \end{equation}

where $h$ denotes the energy distribution function. According to the theory of optical thermodynamics (see Sec. 3) this functional dependence (for lattices with two conservation laws, i.e. power and energy) take the Rayleigh-Jeans form:

(2.15)\begin{equation} h=\frac{a}{\epsilon_j+b} , \end{equation}

where $a$ and $b$ are related to the conservation laws.

Throughout the rest of the paper, equation (2.11) will be augmented with initial condition $c_j(0)$ and zero boundary conditions: $c_0(z)=c_{M+1}(z)=0$. The initial condition is prepared as follows:

(2.16)\begin{equation} c_j(0)\equiv \tilde c_j e^{2\pi i r_j}, \end{equation}

where $\tilde c_j$ is a deterministic positive real amplitude while the phase $r_j$ is a random real number sampled from a uniform distribution on the interval $(0,1)$. Furthermore, the amplitude $\tilde c_j$ is chosen to be one of the following (depending on the case at hand):

• • linear power distribution across all supermodes:

  (2.17)\begin{equation} \tilde c_j=\sqrt{\alpha(\epsilon_j-\epsilon_1)}, \quad \alpha=\frac{P}{\sum_i (\epsilon_i-\epsilon_1)}, \end{equation}

• • piecewise uniform distribution across the supermodes:

  (2.18)\begin{equation} \tilde c_j= \begin{cases} \displaystyle \sqrt{\frac{P}{j_2-j_1}} , \quad j_1\leq j\leq j_2\, ,\\ \\ 0 , \quad \quad \quad \,\,\,\, \text{otherwise}. \end{cases} \end{equation}

Both cases ensure that the total initial power (at $z=0$), and consequently at every $z$, attains the predetermined value $P$ given by

(2.19)\begin{equation} P=\sum\limits_{j=1}^M |c_j|^2. \end{equation}

In what follows, we list some important assumptions about the symmetries of the tensor $T_{j,k,l,m}$.

(2.20)\begin{equation} \text{Quasi-Hermiticity:}\,\, T_{j,k,l,m} = T^*_{l,m,j,k}\, , \end{equation}

(2.21)\begin{align} \nonumber\\ &\text{Permutation symmetry:}\,\, T_{j,k,l,m} = T_{m,l,k,j}. \end{align}

Symmetry condition (2.20) guarantees the conservation of power. The second property (2.21) allows one to derive Eq. (2.11) from Hamilton’s equations of motion for the conjugate pair of canonical variables $c_j,\,c_j^*$:

(2.22)\begin{equation} \frac{dc_j}{dz}=i\{c_j,\mathcal{H}\},\;\; \mathcal{H}=\mathcal{H}_0+\mathcal{H}_{\text{NL}}, \end{equation}

where

(2.23)\begin{equation} \mathcal{H}_0 =\sum\limits_{j=1}^{M}\epsilon_j|c_j|^2, \end{equation}

and

(2.24)\begin{equation} \mathcal{H}_{\text{NL}}= \frac{1}{4}\sum\limits_{jklm}\left(T_{j,k,l,m} c^*_k c^*_l c_m c_j + T^*_{j,k,l,m}c_k c_l c^*_m c^*_j\right). \end{equation}

In Eq. (2.22), $\{\}$ denotes the standard Poisson bracket defined by

(2.25)\begin{equation} \{D,\tilde D\}=\sum_{n=1}^M\left(\frac{\partial D}{\partial c_n}\frac{\partial \tilde D}{\partial c^*_n} - \frac{\partial D}{\partial c^*_n}\frac{\partial \tilde D}{\partial c_n}\right), \end{equation}

where $D$ and $\tilde D$ are two arbitrary functionals of the canonical variables $c_n$ and $c^*_n$. Following this definition, one can derive the Poisson brackets for the respective canonical coordinates of interest, that is, $\{c_j^*,c_n\}=\delta_{j,n}$ and $\{c_j,c_n\}=\{c_j^*,c_n^*\}=0$. We reiterate that, in order to derive Eq.(2.11) using the above Hamiltonian formulation, it is necessary to ensure that the tensor $T_{j,k,l,m}$ remains invariant under the permutation symmetries $k\leftrightarrow l$ and $j\leftrightarrow m$. A detailed derivation of these two symmetries is presented in Appendix A.

5.1. Kerr lattice

Here, we will provide a closed form expression for the tensor $T_{j,k,l,m}\equiv T^{\text{Kerr}}_{j,k,l,m}$ in the presence of Kerr nonlinearity. To do so, we start from Eq. (4.7) with $g=1$ while setting the rest of the parameters to zero. This leads to

(5.1)\begin{equation} i\frac{dA_n}{dz}+ A_{n-1}+A_{n+1} +|A_n|^2A_n=0. \end{equation}

Following similar ideas outlined in Section 2.3 (see also Appendix B for further details) we arrive at

(5.2)\begin{equation} T^{\text{Kerr}}_{j,k,l,m} = B\sum\limits_{n=1}^{M} \prod_{i=1}^4 \sin(q_i x_n)\equiv \frac{B}{16}\sum_{\iota=1}^8 (-1)^{\iota+1}\gamma_{j,k,l,m}^{(\iota)}, \end{equation}

where $B=4/(M+1)^2$ and $x_n=n\pi/(M+1)$. For the rest of the paper we will make a frequent use of the short set-type notation

(5.3)\begin{equation} \{q_i\}_{i=1}^4 \equiv \{j,k,l,m \},\, \, \,\, \{p_i\}_{i=1}^4 \equiv \{j,m,k,l \}. \end{equation}

The auxilary tensor $\gamma_{j,k,l,m}^{(\iota)}$ is given by

(5.4)\begin{equation} \gamma_{j,k,l,m}^{(\iota)} = \Bigg\{ \begin{array}{lr} (-1)^{w_\iota+1}-1, & \text{if } w_\iota \neq 2(M+1)\kappa_\iota\, ,\\ 2M, & \text{if } w_\iota=2(M+1)\kappa_\iota, \end{array} \end{equation}

where $\kappa_{\iota}$ is an arbitrary integer and $w_\iota$ assumes one of the eight distinct integer combinations of the indices $j, k, l, m \in \{1,2,..,M\}$ presented in Table 1. Before analyzing the tensor, it is convenient to list some important properties associated with the integers $w_{\iota}$ that are essential for deriving its closed form representation:

1. (1) If any member of the $\{w_{\iota}\}_{\iota=1}^8$ family is odd, then the rest are also odd.

2. (2) Correspondingly, if any element of the $\{w_{\iota}\}_{\iota=1}^8$ set is even, then all others must also be even. This includes the $w_{\iota}=0$ case.

3. (3) The maximum and minimum values of $\{w_{\iota}\}_{\iota=1}^8$ are $4M$ and $1-3M$, respectively. As a result, the only possibility that gives rise to a multiple of $2M+2$ for any $\iota$ is when $\kappa_{\iota}=0,\pm 1$.

With this at hand, one can show that the tensor now assumes the alternative and surprisingly simple form

(5.5)\begin{equation} T^{\text{Kerr}}_{j,k,l,m}=\frac{1}{2(M+1)}\sum\limits_{\iota=1}^8(-1)^{\iota+1}\delta_{0,\rho_\iota}. \end{equation}

Here, $\rho_{\iota}$ denotes the remainder of the fraction $w_{\iota}/(2M+2)$ and $\delta_{0,\rho_\iota}$ is the Kronecker delta function. Equation (5.5) gives an elegant analytic form of the tensor in terms of the number of supermodes $M$. From the first property of the above list, one can see that $T_{j,k,l,m}^{\text{Kerr}}$ vanishes when (i) one index from the set $\{j,k,l,m\}$ is odd while the rest are even, and (ii) one integer is even and the others are odd. Thus, the only non-zero entries of the tensor would arise only from cases where the elements of the set $\{w_{\iota}\}$ are even (zero included). To derive the non-zero elements of the tensor, we distinguish between two scenarios: (a) $w_1=2M+2$ and (b) $w_1$ being an arbitrary even integer not equal to $2M+2$. Note that these are the only relevant choices imposed by the upper and lower bounds on $w_1$ as well as the tensorial structure of $\gamma^{(\iota)}_{j,k,l,m}$ (see Eq. (5.4)). The alternative (a) would lead to $\gamma^{(1)}_{j,k,l,m}=2M$. This “root”, in turn, gives rise to the tree structure shown in Fig. 1 (left side). In particular, for each subsequent $w_{\iota}$ a tree branch is constructed via the following rules: Determine whether $w_2$ leads to $\gamma^{(2)}=2M$ or negative $2$ which can be established using the upper and lower bounds of $w_2$ (taking into account the constraint imposed by $w_1=2M+2$). Once this is accomplished, this process is repeated. That is to say, find all possible combinations of the rest of $\{w_{\iota}\}_{\iota=3}^8$ (conforming to all constraints listed in Table 1 and, in doing so, record the respective values of their auxiliary tensor $\gamma^{(\iota)}_{j,k,l,m}$. The aforementioned procedure is implemented for option (b) which gives rise to $\gamma^{(1)}=-2$ and the corresponding tree structure (see Fig. 1). After some algebra it can be shown that the tensor values belong to the set

(5.6)\begin{equation} S=\Big\{0,\, \frac{\pm1}{2M+2},\, \frac{1}{M+1},\, \frac{3}{2M+2}, \,\frac{2}{M+1}\Big\}, \end{equation}

Figure 1. The tree diagram associated with the Kerr tensor depicting all possible branches that lead to the set $S$. The labels appearing on the far left vertical line represent the values of the auxiliary tensor $\gamma^{(\iota)}$ defined in Eq. (5.4). For the ease of representation, the tensor indices have been suppressed. The horizontal axis corresponds to the values of the tensor for each respective tree branch as defined by Eq. (5.2). The left tree starts form the “root” $\gamma^{(1)}=2M$ while the right one from $\gamma^{(1)}=-2$. These “roots” originate from the possibilities of $w_1$ precisely equals to $2M+2$ or an arbitrary even integer not equal to $2M+2$ as dictated by its bounds (see Table 1). An example of a tree branch is given by : $2M\to-2\to2M\to-2\to2M\to-2\to-2\to-2$, which leads to tensor value $T^{\text{Kerr}}=3/(2M+2)$. For a fixed set of indices $j,k,l,m$ one can uniquely identify any specific branch on the tree. The branches are labeled from left to right with the first branch appearing on the far left and the 20th branch on the far right.

Table 1. A list of all possible combinations of $w_{\iota}$, $\iota=1,2,\dots,8$ that arise from the derivation of Eq. (5.2) along with their upper and lower bounds (which helps determine if $w_{\iota}$ is equal to an integer multiple of $2M+2$ or not). These possibilities determine the ultimate value of the auxiliary tensor $\gamma_{j,k,l,m}^{(\iota)}$ and, in turn, the tensor $T_{j,k,l,m}^{\text{Kerr}}$ as defined by Eq. (5.2)

where, as a reminder, $M$ denotes the number of supermodes. It is remarkable that the elements of the Kerr tensor appear in such a simple and elegant way, especially on their dependence on the number of supermodes. With these results at hand, equation (2.11) (for even $M$) reads

(5.7)\begin{equation} -i\frac{dc_j}{dz}=\epsilon_jc_j +\frac{1}{M+1}\left(\frac{3}{2}\!\sum_{\substack{\text{branch}\#\\2,3,5,10}}c_kc_lc_m^*+\!\!\!\!\!\!\!\sum_{\substack{\text{branch}\#4,6,7,\\11,12,15}}\!\!\!\!c_kc_lc_m^* +\frac{1}{2}\!\!\sum_{\substack{\text{branch}\#\\8,13,16,18}}\!\!c_kc_lc_m^*-\frac{1}{2}\!\!\sum_{\substack{\text{branch}\#\\9,14,17,19}}\!\!c_kc_lc_m^*\right), \end{equation}

where $\#$ labels the branch number appearing in Fig. 1 oriented from left to right (with branch $\#1$ appearing on the far left while branch $\#20$ to the far right). In what follows, we list a few examples of the above system. When $M=2$ we get

(5.8)\begin{equation} i\frac{dc_1}{dz}+\underbrace{(\epsilon_1+ P-|c_1|^2/2)}_{\text{non-interacting}}c_1+\underbrace{c_2^2c_1^*/2}_{\text{interacting}}=0, \end{equation}

(5.9)\begin{equation} i\frac{dc_2}{dz}+\underbrace{(\epsilon_2+ P-|c_2|^2/2)}_{\text{non-interacting}}c_2+\underbrace{c_1^2c_2^*/2}_{\text{interacting}} =0, \end{equation}

where again $P=|c_1|^2+|c_2|^2$ is the total power (which is a conserved quantity). If one only keeps the non-interacting terms, we find $\langle|c_j(z)|^2\rangle=|c_j(0)|^2$ which indicates lack of thermalization. Thus, the middle terms in Eq. (5.8) and Eq. (5.9), do not contribute to the wave mixing process. However, the presence of the last interacting terms produces a nonlinear wave mixing, which would prohibit $\langle|c_j(z)|^2\rangle$ from remaining constant. The equations for the projection coefficients $c_j$ become slightly involved when the number of supermodes increases. This can be seen for the $M=3$ case (see Appendix B for more details). Here, one can also identify the non-interacting ( ${\cal N}_j$) and interacting terms ( ${\cal M}_j$) for which the evolution of $c_j$ is governed by

(5.10)\begin{equation} i\frac{dc_j}{dz}+{\cal N}_jc_j+{\cal M}_j(c_j)=0. \end{equation}

As the number of supermodes increases, the equation of motion governing the dynamics of the projection coefficients becomes more involved. This is clearly demonstrated trough the tensor structure depicted in Fig. 2 for five supermodes.

Figure 2. The tensor $T_{jklm}$ appearing in Eq.(2.11) corresponds to the Kerr lattice as defined by Eq. (5.5) (for simplicity we only shown the slices associated with $j=1,3,5$) when the number of supermodes is $M=5$. The index $k$ counts the number of tensor slices oriented from left to right. For a fixed tensor slice the indices $l$ and $m$ denote the number of its rows and columns, respectively.

This example clearly illustrates the important symmetry structure of the Kerr tensor mentioned earlier along with the values it assumes. Below, we summarize several properties associated with $T_{j,k,l,m}^{\text{Kerr}}$:

• • As one can see in Fig. 2, at least half of the tensor elements are zero. This is due to the first property of the set $\{w_{\iota}\}$ i.e., if any member of the $w_{\iota}$ family is odd, then the rest are also odd. The locations where the tensor vanishes can be identified when an individual index from the set $\{j,k,l,m\}$ is odd (even) while the rest are even (odd).

• • The tensor elements are $0,\pm1/12, 1/4, 1/6, 1/3$ which coincides with the set $S$ when $M=5$. These numerical values can be found from the tree structure shown in Fig. 1.

• • For any two fixed indices, the resulting tensor slice forms a Hermitian matrix. This property is due to the high symmetry of the tensor $T_{jklm}^{\text{Kerr}}$. As we shall later see, such symmetry is absent for other cubic lattices.

• • The tensor shown in Fig. 2 obeys the quasi-Hermiticity and permutation symmetries mentioned in Sec. 2.3 i.e., invariance under the transformations $k\leftrightarrow l$, $j\leftrightarrow m$, and $k\leftrightarrow m$, $l\leftrightarrow j$.

From the above observations, we expect the modal occupancies $|c_j|^2$ to follow a Rayleigh-Jeans distribution once thermal equilibrium is reached. In fact, we have performed numerical simulations using Eq. (2.11) for $M=20$. The results are depicted in Fig. 3. This in turn confirms the theory of optical thermodynamics which also agrees with the results obtained from direct simulations using the local base [Reference Parto, Wu, Jung, Makris and Christodoulides39]. To this end, we point out that the tensor analysis presented here will lay the foundation for our study of thermalization in random tensors discussed in detail in Sec. 6.2.

Figure 3. Modal occupancies versus eigenvalues for the Kerr lattice with $M=20$, power $P=2$ and energy $H_{0}=-1.9$. The initial power distribution among the various supermodes is shown in green diamonds, see Eq. (2.17). The numerical results (red stars) indicate the modal occupancies averaged over 400 realizations of random phases evaluated at propagation distance $z=10000$. These results were obtained by simulating equation (2.11) where the (Kerr) tensor $T_{j,k,l,m} $ is given by Eq. (5.5) (see Fig. 1 that helps construct the tensor values). The theory of optical thermodynamics for those values predicts a temperature $T=0.153$ and a chemical potential $\mu=-2.48$. In this case the Rayleigh–Jeans distribution Eq. (3.2) is also shown in a solid blue line.

5.2 Nonlocal lattices

In this section we embark on the task to build an analytic form of the mixing tensor for several nonlocal optical lattices associated with Eq. (4.7). The main motivation is the identification of the tensorial symmetries which will be later tied to thermalization.

5.2.1 Case I

We first consider the case where $a_1=1,\,\,g=a_2=b_1=b_2=0$. In this situation, the evolution of the optical field $A_n$ is governed by

(5.11)\begin{equation} i\frac{dA_n}{dz}+A_{n-1}+A_{n+1}+A_n\left(|A_{n-1}|^2+|A_{n+1}|^2\right)=0. \end{equation}

The nonlocal nonlinear terms appearing in Eq. (5.11) adhere to the general form presented in Eq. (2.3), with $n_1=0, \, n_2=n_3=\pm 1$. As such, the tensor given in Eq. (2.12) will have two contributions leading to

(5.12)\begin{equation} T_{j,k,l,m} = B \sum\limits_{n=1}^{M} \left(\prod_{i=1}^2 \sin(q_i x_n) \sin(q_{i+2}x_{n+1}) \right. + \left. \prod\limits_{i=1}^2 \sin(q_i x_{n+1}) \sin(q_{i+2}x_{n})\right), \end{equation}

where $x_n=n\pi/(M+1)$. As a reminder, we refer to Eq. (5.3) for the definition of $q_i$. Interestingly enough, one can relate the tensor in Eq. (5.12) to the Kerr one as follows:

(5.13)\begin{equation} T_{j,k,l,m} = 2T_{j,k,l,m}^{\text{Kerr}}\cos(x_l)\cos(x_m)+ 2\Gamma^{(1)}_{j,k,l,m}\sin(x_l)\sin(x_m), \end{equation}

where,

(5.14)\begin{equation} \Gamma^{(1)}_{j,k,l,m} =B\sum\limits_{n=1}^M \prod_{i=1}^{2}\sin\left(q_ix_n\right)\,\cos\left(q_{i+2}x_n\right). \end{equation}

One can further simplify Eq. (5.14) to obtain an alternative closed form similar to the Kerr case given in Eq.(5.5). Doing so, we get

(5.15)\begin{equation} \Gamma^{(1)}_{j,k,l,m} = -\frac{B}{16}\sum_{\iota=1}^4 (-1)^{\iota+1}\left(\gamma_{j,k,l,m}^{(\iota)}-\gamma_{j,k,l,m}^{(\iota+4)}\right). \end{equation}

Here, the auxiliary tensor $\gamma_{j,k,l,m}^{(\iota)}$ is defined in Eq. (5.4). From Eq. (5.15) and the tree structure depicted in Fig. 1, it can be shown that the tensor elements associated with $\Gamma^{(1)}_{j,k,l,m}$ assume the values given by the set

(5.16)\begin{equation} {\cal S}=\left\{0, \frac{\pm1}{M+1},\frac{\pm1}{2M+2}\right\}. \end{equation}

With this result at hand, an example that illustrates the structure of the tensor $T_{j,k,l,m}$ is shown in Fig. 4. Scrutinizing Eq. (5.12) reveals that the tensor satisfies the two symmetry conditions given in Eqs. (2.20) and (2.21). With this fact at hand, we pose the question: Will the lattice given in Eq. (5.11) thermalize to a RJ or not? To answer this, we have numerically solved Eq. (5.11) subject to the initial random distribution given in (2.17). Our findings are shown in Fig. 5 (a). One can see that a RJ distribution has been reached validating our tensorial symmetry-thermalization connection.

Figure 4. An illustrative example of the tensor $T_{j,k,l,m}$ given in Case I of Sec. 5.2 with $M=5$. For the ease of presentation we only show the cases with $j=1$ and $j=5$.

Figure 5. Evolution of the ensemble averaged modal occupancies for various nonlocal nonlinear lattices with $M=100$ supermodes. (a) The first nonlocal lattice (Case I) with total power of $P=8$ and linear initial distribution across the modes (indicated by a solid black line at $z=0$). The lattice approaches a Rayleigh-Jeans distribution (red line at $z=10000$), agreeing with the theoretically predicted values for the temperature $T=0.19$ and a chemical potential $\mu=-2.7$. (b) and (c) correspond to the nonlocal cases II and III, respectively. Here, the initial condition corresponds to an equally distributed power of $P=8$ among the supermodes $65\leq j \leq 89$ (b) and $P=5$ among $15\leq j \leq 39$ (c) supermodes. Both simulations lead to a RJ distribution with temperatures and chemical potentials $T=0.066$, $\mu=-2.15$ (b) and $T=-0.1$, $\mu= -0.7$ (c) as the theory predicts. Lastly part (d) corresponds to the nonlocal case IV with power $P=5$ linearly distributed (at $z=0$) across the supermodes with indices $50\leq j \leq 100$. The temperature at thermal equilibrium is $T=0.07$ and the chemical potential $\mu=-2.4$. The red solid line indicates the RJ distribution at $z=10^4$. For each subfigure, a colored point represents the modal occupancy at a given eigenvalue and are connected with a gray dashed line at each $z$ cross section.

5.2.2 Case II

We next consider the DNLS model given in Eq. (4.7) with $a_2= 1$ while setting the rest of the parameters to zero. Thus, we have

(5.17)\begin{equation} i\frac{dA_n}{dz}+A_{n-1}+A_{n+1}+A_n^*\left(A_{n-1}^2+A_{n+1}^2\right)=0. \end{equation}

In this case, the mixing tensor assumes the form

(5.18)\begin{equation} T_{j,k,l,m} = B \sum\limits_{n=1}^{M}\left( \prod_{i=1}^2 \sin(p_i x_n) \sin(p_{i+2}x_{n+1}) \right. + \left. \prod_{i=1}^2 \sin(p_i x_{n+1}) \sin(p_{i+2}x_{n}) \right), \end{equation}

where $\{p_{i}\}_{i=1}^4 \equiv \{j,m,k,l\}$. In a similar fashion to the previous nonlocal case, one can express the above lattice mixing tensor in terms of the Kerr one. That is,

(5.19)\begin{equation} T_{j,k,l,m} = 2T_{j,k,l,m}^{\text{Kerr}}\cos(x_k)\cos(x_l)+2\Gamma^{(2)}_{j,k,l,m}\sin(x_k)\sin(x_l), \end{equation}

where,

(5.20)\begin{equation} \Gamma^{(2)}_{j,k,l,m} =B\sum\limits_{n=1}^M \prod_{i=1}^{2}\sin\left(p_ix_n\right)\,\cos\left(p_{i+2}x_n\right), \end{equation}

which after some simplifications becomes

(5.21)\begin{equation} \Gamma^{(2)}_{j,k,l,m} = -\frac{B}{16}\sum_{\iota=1}^2 \left(\gamma_{j,k,l,m}^{(\iota)}+\gamma_{j,k,l,m}^{(\iota+4)}\right)+\frac{B}{16}\sum_{\iota=1}^2\left(\gamma_{j,k,l,m}^{(\iota+2)}+\gamma_{j,k,l,m}^{(\iota+6)}\right). \end{equation}

The entries of the tensor $\Gamma^{(2)}_{j,k,l,m}$ also belong to the set ${\cal S}$ given in Eq. (5.16). The thermalization properties and their connection to the tensor symmetries will be discussed in detail in Sec. 6.

5.2.3 Case III

Here, we consider an alternative type of nonlocal lattice governed by Eq. (4.7) for which $b_1=1$ and the rest of the model parameters are set to zero. In such scenario, the evolution equation for the optical field envelope $A_n$ in local base reads

(5.22)\begin{equation} i\frac{dA_n}{dz}+A_{n-1}+A_{n+1} +A_n^2A_{n+1}^*+ 2|A_n|^2A_{n+1}+|A_{n-1}|^2A_{n-1}=0\, . \end{equation}

When viewed in supermode space, the resulting tensor appearing in Eq. (2.11) assumes the form

(5.23)\begin{equation} T_{j,k,l,m} = B \sum\limits_{n=1}^{M}\left( \prod_{i=1}^3 \sin(q_i x_n) \sin(m x_{n+1}) \right. + \left. \prod_{i=2}^4 \sin(q_i x_{n})\sin(j x_{n+1})\right.\left.+2\prod_{i=1}^3 \sin(p_i x_{n})\sin(l x_{n+1}) \right). \end{equation}

In order to identify the underlying symmetries of the tensor, we can split the last product in Eq. (5.23) as a sum of two contributions that would uncover the required permutation and quasi-Hermiticity symmetries. In terms of the Kerr tensor, Eq. (5.23) assumes the alternative form

(5.24)\begin{equation} T_{j,k,l,m}= T_{j,k,l,m}^{\text{Kerr}}\sum_{i=1}^4\cos(x_{q_i})+\sum_{i=1}^4\Xi_{j,k,l,m}^{(i)}\sin(x_{q_i}), \end{equation}

with

(5.25)\begin{equation} \Xi_{j,k,l,m}^{(i)}=-\frac{B}{8}\sum_{\iota = 1}^8\beta_{\iota}^{(i)}\zeta_{j,k,l,m}^{(\iota)}. \end{equation}

Here, $\beta^{(1)}_{\iota}=1$, when $\iota=1,4,5,8$ and $-1$ otherwise; $\beta^{(2)}_{\iota}=1$, for $\iota=1,2,7,8$ and $-1$ for the remaining cases; $\beta^{(3)}_{\iota}=1$, if $\iota=1,3,6,8$, else it is equal to $-1$ and finally $\beta^{(4)}_{\iota}=\pm 1$, where plus/minus sign corresponds to odd/even integer $\iota$, respectively. Furthermore, the auxiliary tensor $\zeta^{(\iota)}$ is given by

(5.26)\begin{equation} \zeta_{j,k,l,m}^{(\iota)} = \Bigg\{ \begin{array}{lr} \cot\left(\frac{\pi w_{\iota}}{2(M+1)}\right), & \text{if}\,\, w_\iota\,\,\text{odd,}\\ 0, & \text{if}\,\, w_\iota\,\,\text{even,} \end{array} \end{equation}

where, again, the $w_{\iota}$ are given in Table 1. A detailed derivation of the above results can be found in Appendix B. A few typical tensor slices for a small number of supermodes ( $M=3$) are

(5.27)\begin{equation} T_{1,1,l,m}=\frac{1}{4} \begin{pmatrix} 3\sqrt{2} & 1/2 & -\sqrt{2}/2\\ 1/2 & \sqrt{2} & 3/2\\ -\sqrt{2}/2 & 3/2 &0 \end{pmatrix} , \end{equation}

(5.28)\begin{equation} T_{1,2,l,m}=\frac{1}{4} \begin{pmatrix} 1/2 & \sqrt{2} & 3/2\\ \sqrt{2} & 1 & 0\\ 3/2 & 0 & -3/2 \end{pmatrix} , \end{equation}

(5.29)\begin{equation} T_{1,3,l,m}=\frac{1}{4} \begin{pmatrix} -\sqrt{2}/2 & 3/2 & 0\\ 3/2 & 0 & -3/2\\ 0 & -3/2 & \sqrt{2}/2 \end{pmatrix} . \end{equation}

It is remarkable that the entries of these tensor slices are rather simple looking which is unexpected given the fact that the type on nonlinearity in local space is rather non-trivial. This pattern also shows up if one slices the tensor along a different hyperplane.

5.2.4 Case IV

Lastly, we analyze the structure of the nonlinear mixing tensor for the generalized DNLS equation (4.7). Here we take $b_2=1$ and choose the rest of the model parameters to zero. Under such assumptions, the governing dynamical system is reduced to

(5.30)\begin{equation} i\frac{dA_n}{dz}+A_{n-1}+A_{n-1} +A_n^2A_{n-1}^*+ 2|A_n|^2A_{n-1}+|A_{n+1}|^2A_{n+1}=0. \end{equation}

Following similar analysis discussed in previous cases, after some algebra we arrive at

(5.31)\begin{equation} T_{j,k,l,m} = B \sum\limits_{n=1}^{M} \left(\prod_{i=1}^3 \sin(q_i x_n)\sin(m x_{n-1}) \right. + \left. \prod_{i=2}^4 \sin(q_i x_{n})\sin(j x_{n-1}) \right. + \left.2\prod_{i=1}^3 \sin(p_i x_{n}) \sin(l x_{n-1}) \right) \, . \end{equation}

Using similar arguments as discussed in Case III one can establish the validity of the quasi-Hermiticity and permutation symmetries of the tensor $T_{j,k,l,m}$. It is again interesting to note that the tensor under discussion is related to its Kerr counterpart. As such, we have

(5.32)\begin{equation} T_{j,k,l,m}= T_{j,k,l,m}^{\text{Kerr}}\sum_{i=1}^4\cos(x_{q_i})-\sum_{i=1}^4\Xi_{j,k,l,m}^{(i)}\sin(x_{q_i}). \end{equation}

It is evident that the tensor structure in Eq. (5.32) is intricate which makes its simplification a rather formidable task. Nonetheless some representative examples of tensor slices when $M=3$ that highlight its intrinsic structure are given:

(5.33)\begin{equation} T_{1,1,l,m}=\frac{1}{4} \begin{pmatrix} 3\sqrt{2} & -1/2 & -\sqrt{2}/2\\ -1/2 & \sqrt{2} & -3/2\\ -\sqrt{2}/2 & -3/2 &0 \end{pmatrix} , \end{equation}

(5.34)\begin{equation} T_{1,2,l,m}=\frac{1}{4} \begin{pmatrix} -1/2 & \sqrt{2} & -3/2\\ \sqrt{2} & -1 & 0\\ -3/2 & 0 & 3/2 \end{pmatrix} , \end{equation}

(5.2.)\begin{equation} T_{1,3,l,m}=\frac{1}{4} \begin{pmatrix} -\sqrt{2}/2 & -3/2 & 0\\ -3/2 & 0 & 3/2\\ 0 & 3/2 & \sqrt{2}/2 \end{pmatrix} \end{equation}

Similar arguments can be made about the simplicity of these tensor slices given the fact that a priori one would have expected more complicated results.

In the next section, we will address the issue of “correlation” between the above-mentioned symmetries and equilibration to a RJ distribution of all the nonlocal models discussed above.

6.1. Connections to tensorial symmetries

One of the main goals of the current study is to investigate, through systematic numerical simulations, a possible connection between the symmetries of the tensors $T_{j,k,l,m}$ that appear in Eq. (2.11) and the thermalization properties of the corresponding cubic nonlinear lattices. In other words, we would like to pose the question whether nonlinear tensors preserving the quasi-Hermiticity and permutation symmetries assumed in Eq. (2.20) and (2.21) lead to a Rayleigh-Jeans equilibrium distribution. If this unidirectional hypothesis (tensorial symmetries imply thermalization) holds, it would provide a useful tool to investigate thermalization of cubic lattices solely by scrutinizing their nonlinear tensorial symmetries.

With this in mind, we shall use the discrete model Eq. (4.7) as a testbed example for which all the underlying tensorial structures have been analytically derived. Furthermore, these tensorial symmetry conditions can be readily verified by inspecting their exact closed from provided in Sec. 5. To this goal, we next perform direct numerical simulations on Eq. (4.7) with $V_n=0$ corresponding to random initial conditions of the form:

(6.1)\begin{equation} A_n(0)=\sum\limits_{j=1}^M \tilde c_je^{2\pi i r_j}\psi_n^{(j)}, \end{equation}

where $r_j$ is a random field uniformly distributed on the interval $(0,1)$ and $\tilde c_j$ are deterministic modal amplitudes defined by either Eq. (2.17) or (2.18). Below, we report on our numerical findings for all lattices (Kerr and nonlocal) in the same order as presented in Sec. 5.

1. (1) Kerr Lattice. The Hamiltonian is given by

   (6.2)\begin{equation} H =H_0+\frac{1}{2}\sum\limits_{n=1}^M|A_n|^4, \end{equation}

   with $H_0$ being the part corresponding to the nearest neighbor coupling. Thermalization properties of the Kerr lattice have been thoroughly studied over the past several years [Reference Kottos and Shapiro25, Reference Parto, Wu, Jung, Makris and Christodoulides39]. It is well known that equilibration to a RJ distribution has been reported (see for example Fig. 3). The Kerr tensor given in Eq. (5.2) does satisfy the quasi-Hermiticity and permutation symmetries postulated. As such, the Kerr nonlinear lattice conforms with this hypothesis. It is interesting to note that the Kerr tensor exhibits high symmetry, i.e., invariance under any index permutation. Compared with the Kerr case, all other nonlocal tensors admit reduced symmetries.

2. (2) Nonlocal Lattices.

• • Case I. This is the first example, which we use to test a possible link between the constructed tensorial symmetries and their implication on relaxation to a RJ distribution. To do so, we consider the dynamical lattice given in Eq. (5.11) for which the corresponding Hamiltonian reads

  (6.3)\begin{equation} H =H_0+\sum\limits_{n=1}^M|A_n|^2|A_{n+1}|^2. \end{equation}

  Note that the resulting mixing tensor and its associated symmetries have been exactly identified in Sec. 5.2. Thus, to verify the symmetry-thermalization connection, we need to resort to numerical simulations performed on the given Hamiltonian. The results of such computation are summarized in Fig. 5 (a), where the dependence of the quantity $\langle|c_j|^2\rangle$ on the system’s linear eigenenergies is shown. The system size, total (dimensionless) propagation distance, and number of ensembles used in this and the following cases are fixed at $M=100$, $z=10^4$ and 800, respectively. These parameters were selected to ensure sufficient number of linear modes, propagation distance and the number of samples large enough to observe equilibration. As it can be seen in every case presented in Fig. 5 the relaxation to an equilibrium state start to develop approximately half the total distance used in these simulations. Additionally, the input power is always chosen so that the ratio $\mathcal{H}_{0}/\mathcal{H}$ (see Eq. (2.22) and Eq. (2.23)) is close to one. For the simulation performed in this paper we have used this ratio to be, e.g., 0.9. This in turn guarantees the dynamics remain predominantly linear with the nonlinear contribution to the total energy remaining small. As one can see, its form coincides with theoretically predicted RJ distribution which supports the above-mentioned conjecture.

• • Case II. The corresponding Hamiltonian for the system given in (5.17) is given by

  (6.4)\begin{equation} H =H_0+\frac{1}{2}\sum\limits_{n=1}^M\big(A_n^2(A_{n+1}^*)^2 + \text{c.c.}\big). \end{equation}

  We have performed direct numerical simulations using the above Hamiltonian to check if thermalization to RJ is possible. In Fig. 5 (b) we show the ensemble averaged modal occupancies in terms of the energy eigenstates. It is clear that the system indeed relaxes to a RJ distribution. Given the fact that the tensor corresponding to this Hamiltonian has been computed in Sec. 5.2 (which obeys the quasi-Hermiticity and permutation symmetries), then one can positively affirm the unidirectional connection between tensorial symmetries and thermalization.

• • Case III. In this situation, the evolution equation of the optical field $A_n$ is governed by Eq. (5.22) which can be derived from the Hamiltonian functional

  (6.5)\begin{equation} H = H_0+ \sum\limits_{n=1}^M \left(A_n^*A_{n+1}^*A_n^2+\text{c.c.}\right) \, , \end{equation}

  using the standard Poisson brackets defined in Eq. (2.25). We have studied the possibility of thermalization to a distribution that follows the Rayleigh-Jeans law by simulating the above Hamiltonian system subject to uniform modal distribution (see Fig. 5 (c)). Indeed, the system reached a RJ distribution which coincides with the theoretically predicted formula of optical thermodynamics given in Eq. (3.2). Interestingly enough, the corresponding tensorial symmetries found in Sec. 5.2 appear to conform with our symmetry-thermalization connection.

• • Case IV: Lastly, we probe the tensorial symmetries and their connection to thermalization by considering another type of nonlocal nonlinear lattice as given in Eq. (5.30) with Hamiltonian structure

  (6.6)\begin{equation} H = H_0+ \sum\limits_{n=1}^M \left(A_n^*A_{n+1}^*A_{n+1}^2+\text{c.c.} \right). \end{equation}

  As before, we performed numerical simulations on the above Hamiltonian and found that the statistically averaged modal occupancies agree with the theoretically predicted RJ distribution, see Fig. 5 (d). Furthermore, the associated tensor preserves the two postulated symmetries given in Eq. (2.20) and (2.21). This case provides another indication that supports the link between tensorial symmetries and RJ distribution.

6.2. Random tensors

Numerical simulations related to thermalization of nonlinear lattices often require less computational efforts when carried out in the local base rather than in their modal representation. One of the main reasons being the intricate structure of the tensor $T_{j,k,l,m}$ and the combinatorial sum that appears in Eq. (2.11). For example, for a system with $M$ supermodes, the number of nonlinear terms is of the order of $M^3$ and the tensor is of size $M^4$. Nonetheless, Eq. (2.11) offers several advantages:

• • Provides an elegant and unified formulation for all cubic nonlinear lattices. In other words, Eq. (2.1) with an arbitrary number of nonlinear cubic terms would look the same as in Eq. (2.11) with the proper tensor.

• • The nonlinear evolution equation in supermode base (Eq. (2.11)) allows one to probe new nonlinear systems without any direct reference to the dynamics on a local base.

• • Opens the opportunity to study thermalization for random tensors. This highlights the importance of the quasi-Hermiticity and permutation symmetries (which are not visible in the local base) in the theory of thermalization of nonlinear lattices. In other words, it provides a unique testbed for our hypothesis, which asserts that tensorial symmetries lead to thermalization.

The study of thermalization with random tensors provides a unique opportunity to test the hypothesis formulated earlier regarding the role that tensor symmetries play in thermalization. We thus, consider Eq. (2.11) where now the tensor is replaced by a random tensor i.e.,

(6.7)\begin{equation} T_{j,k,l,m}\longrightarrow T^{\text{rand}}_{j,k,l,m}. \end{equation}

The elements of the tensor $T^{\text{rand}}_{j,k,l,m}$ are generated from a uniform probability distribution on the interval $(-0.1,0.1)$ such that they adhere to the presumed symmetry conditions:

(6.8)\begin{equation} \text{quasi-Hermiticity: }\,\,\,T^{\text{rand}}_{j,k,l,m} = T^{\text{rand}}_{l,m,j,k}\, , \end{equation}

(6.9)\begin{equation} \text{permutation symmetry: }\,\,\, T^{\text{rand}}_{j,k,l,m} = T^{\text{rand}}_{m,l,k,j}\, . \end{equation}

We have performed numerical simulations using Eq. (2.11) with the random tensor given in Eq. (6.7). A summary of our results is depicted in Fig. 6 (a), where the ensemble averaged modal occupancies $\langle|c_j|^2\rangle$ over many realizations of initial random modal phases is shown to relax to a RJ distribution. These numerical findings highlight the important role that the tensorial quasi-Hermiticity and permutation symmetries play in thermalization of cubic lattices. It is interesting to mention that if one relaxes the permutation symmetry condition (while preserving the quasi-Hermiticity) then the system no longer approaches a RJ distribution. Instead, it reaches an equipartition state as is seen in Fig. 6 (b). The results corroborate our proposition that the intrinsic tensorial symmetries play a key role in determining the ultimate functional form of the statistical distribution of the modal amplitudes, as well as in the overall process of thermalization of nonlinear cubic optical lattices.

Figure 6. Modal occupancies over eigenvalues diagram for a random tensor that preserves (a) quasi-Hermiticity and permutation symmetries (b) only the quasi-Hermiticity condition. The total power $P=1$ is initially linearly distributed across $M=20$ modes in an ascending order from higher to lower order modes (purple/green line, respectively). After 400 ensembles over random phase initial conditions the distribution of the modal occupancies is depicted (black stars) at $z=5000$. In (a) the solid blue line is the theoretically predicted RJ distribution with temperature $T=0.08$ and chemical potential $\mu=-2.5$. In (b) a RJ distribution is not observed. Instead, the system tends to converge toward an equilibrium state characterized by power equipartition.

6.3. Dependence on supermodes

So far, we have studied the structure of the mixing tensor $T_{j,k,l,m}$ that appears in Eq. (2.11) for various types of nonlinear cubic lattices with nearest neighbors coupling in the absence of any external potentials. In this case, each supermode given by Eq. (2.9) occupies all lattice sites. This implies the following

\begin{equation*} \text{extended supermodes } \longrightarrow \text{extended mixing tensor.} \end{equation*}

The above statement is valid irrespective of whether the nonlinearity (in local base) is short (e.g. Kerr) or long ranged (such as the nonlinearity in Eq. (4.7)). When the tensor is dense (like the examples presented in this paper), the dimension of the mixing terms in supermode base is large, causing the system to reach a RJ distribution relatively fast (see the numerical results presented in Sec. 6.1). It is important to note here that this conclusion is based solely on numerical experiments performed on Eq. (4.7). Other dynamical systems, like the FPUT lattice, might not conform to such a behavior [Reference Berchialla, Giorgilli and Paleari13]. The situation becomes even more intricate when the supermodes are in a localized state, as is the case when the potential $V_n$ corresponds to either the Anderson disordered model [Reference Anderson7] or the Aubry-André quasi-periodic case [Reference Aubry and André8]. The intricacies can be attributed to the localization-delocalization properties of the supermodes in the presence of Kerr nonlinearity in the presence of Anderson and Aubry-André type potentials [Reference Milovanov and Iomin36, Reference Pikovsky and Shepelyansky42]. In other words, the sparseness of the tensor depends on the localization length of the supermodes as well as on the type of nonlinearity. That is to say

\begin{equation*} \text{strongly localized tensors } \longrightarrow \text{slow thermalization}\, \end{equation*}

In Fig. 7 we show a typical example of the mixing tensor for various strength of disorder lattice potentials uniformly distributed on the interval $(-W,W)$. As one can see, for small values of disorder ( $W=0.1$), the tensor slice is in an extended state, whereas it becomes progressively more localized as the disorder strength gets larger. The punchline is the following: The less (more) the supermodes are localized, the sparser (denser) the mixing tensor becomes, which in turn forces the thermalization process to evolve slower (faster). This behavior has been observed on the Anderson model with Kerr nonlinearity [Reference Pyrialakos, Wu, Jung, Ren, Makris, Musslimani, Khajavikhan, Kottos and Christodoulides46] and it has also been reported in other fields beyond nonlinear optics, such as coupled nonlinear disordered oscillators [Reference Basko10, Reference Dhar and Lebowitz18, Reference Lando and Flach30, Reference Mulansky, Ahnert, Pikovsky and Shepelyansky37, Reference Oganesyan, Pal and Huse38, Reference Senyange, Manda and Skokos53], quantum many body systems [Reference Bhattacharjee, Andreanov and Flach15, Reference Lando and Flach29, Reference Zhang, Lando, Dietz and Flach60] and physics related to superfluidity like Bose-Einstein condensates [Reference Billy, Josse, Zuo, Bernard, Hambrecht, Lugan, Clément, Sanchez-Palencia, Bouyer and Aspect16]. In the next section, we shall elaborate more on this issue through the perspective of the mixing tensor.

Figure 7. Typical tensor slices (in absolute value) for the nonlinear Anderson model with random potential $V_n$ uniformly distributed on the interval $(-W,W)$. (a) $W=0.1$, (b) $W=1$ (c) $W=5$. The number of supermodes is 60. The number of nonzero tensor elements decreases as the potential disorder increases resulting in a small number of nonlinear mixing terms (i.e., short range nonlinearity).

6.4. Optical lattices in modal space: an inverse approach

The standard approach to the study of thermalization of nonlinear lattices follows a prescribed dynamics in local base which often times is governed by a “nice” set of equations, such as the ones considered in this paper. This implies that the evolution of the projection coefficients $c_j(z)$ obeys a rather complicated coupled dynamical system. As such, the number of nonlinear terms is of the order of $M^3$, making their analysis difficult. This is the case for our current study: the nonlinear interactions in local base appear in an elegant form while “looking messy” in supermode base. Since the role of weak nonlinearity is to steer the system into an equilibrium state (whenever it exists) and the nonlinearity type is “irrelevant” (Kerr or nonlocal) this raises the question whether one can devise a method whereby the nonlinear interactions in supermode base look “simple” at the expense of dealing with “unpleasant” nonlinearities in local base. Since the quantity of interest is given by $\langle|c_j|^2\rangle$, this inverse approach would provide an advantage to probe the role of mixing tensors in reaching the equilibrium state.

The idea is schematically shown in Fig. 8. To this end, we aim at deriving a large class of lattice models that are embedded in the supermode base containing short range nonlinearities that preserve the power and Hamiltonian. The starting point is the evolution equation (2.11) which for simplicity we write again

(6.10)\begin{equation} i\frac{dc_j}{dz}+\epsilon_j c_j+\sum_{k,l,m=1}^{M} \mathcal{T}_{j,k,l,m}\,c_k c_l c^*_m=0, \end{equation}

Figure 8. A schematic presentation of the forward and inverse approaches to the study of thermalization. In the conventional case, the dynamics of the wave amplitude at site $n$ in local base contains short range nonlinearities (e.g., Kerr, and/or nonlocal terms, see Eq. (4.7)). Since the supermodes are in an extended state, under the transformation $A_n\rightarrow c_j$ the dynamics in supermode base contains long range nonlinear couplings. Contrary to this, if one begins with an evolution equation for the projection coefficients $c_j$ under the assumption of short range interactions then the inverse transform $c_j\rightarrow A_n$ produces long range nonlinear coupling in the local base.

with the lattice on-site energies $\epsilon_j$ given by either Eq. (2.8), Anderson random type, or the Aubry-André model. Furthermore, it is assumed that the tensor $\mathcal{T}_{j,k,l,m}$ is constructed from

(6.11)\begin{equation} \mathcal{T}_{j,k,l,m}=\sum_{n=1}^M\prod_{q_i=1}^4 \Phi_n^{(q_i)}\, , \, \, \{q_i\}_{i=1}^4 = \{j,k,l,m \}, \end{equation}

with

(6.12)\begin{equation} \Phi_n^{(q_i)} =\sum_{\xi=0}^{N-1}\alpha_{\xi}\delta_{n-\xi,q_i}, \end{equation}

where $N\leq M$ counts the number of peaks of the wavefunction $ \Phi_n^{(q_i)}$ with real amplitude $\alpha_{\xi}$ and $\delta_{i,j}$ is the Kronecker delta function. Substituting Eqs. (6.11), (6.12) into Eq. (6.10) we arrive at the evolution equation

(6.13)\begin{equation} i\frac{dc_j}{dz}+\epsilon_jc_j+\sum_{\zeta=0}^{N-1}\left(\alpha_{\zeta}\Bigg|\sum_{\xi=0}^{N-1}\alpha_{\xi}c_{j+\zeta-\xi}\Bigg|^2\bigg(\sum_{\xi=0}^{N-1}\alpha_{\xi}c_{j+\zeta-\xi}\bigg)\right)=0. \end{equation}

Equation (6.13) is a novel set of dynamical systems that provides an opportunity to probe the effect of nonlinear wave mixing on thermalization processes. To reconstruct the corresponding lattice equation in local base we use the inverse approach described in Fig.8. This leads to

(6.14)\begin{equation} i\frac{dA_n}{dz}+A_{n+1}+A_{n-1}+\sum_{k,l,m=1}^MT_{\text{dual}}^{n,k,l,m}A_kA_lA_m^*=0, \end{equation}

where

(6.15)\begin{equation} T_{\text{dual}}^{n,k,l,m}=\sum_{j=1}^M\sum_{\{\xi\}=0}^{N-1}\alpha_{\xi_0}\alpha_{\xi_1}\alpha_{\xi_2}\alpha_{\xi_3}\psi_n^{(j)}\psi_k^{{(j+\xi_0-\xi_1)}^*}\psi_l^{{(j+\xi_0-\xi_2)}^*}\psi_m^{(j+\xi_0-\xi_3)}, \end{equation}

and $\{\xi\}$ stands for a short notation for $\{\xi_0,\xi_1,\xi_2,\xi_3\}$. In deriving equations (6.14) and (6.15) we used $\epsilon_j$ that corresponds to the free lattice in which $\psi_n^{(j)}$ is given by Eq. (2.9). Instead, one can obtain the corresponding local equation in the presence of any $\epsilon_j$ set of eigenvalues, such as that of the Anderson model. To this end, we next present several examples in support of the hypothesis that long range nonlinear couplings in supermode base facilitate nonlinear wave mixing, which eventually leads to faster thermalization.

• • $N=1$: This case corresponds to $\Phi_n^{(q_i)}$ with a single peak (extreme localization regime) for $\alpha_0=1$. As a result Eq. (6.13) becomes

  (6.16)\begin{equation} i\frac{dc_j}{dz}+\epsilon_jc_j+|c_j|^2c_j=0. \end{equation}

  Consequently, $|c_j|^2$ is a constant of motion and thus the initial distribution among the $M$ supermodes remains invariant, i.e.,

  (6.17)\begin{equation} \langle|c_j(z)|^2\rangle= \langle|c_j(0)|^2\rangle. \end{equation}

  In other words, in the extreme localization regime there is no relaxation to a RJ distribution.

• • $N=2$: Here, the wave functions $\Phi_n^{(q_i)}$ are composed of two peaks (for sake of simplicity are taken to be of equal heights with $\alpha_0=\alpha_1=1$). As a result, Eq. (6.13) takes the surprisingly simple form:

  (6.18)\begin{equation} i\frac{dc_j}{dz}+\epsilon_jc_j+|c_{j}+c_{j-1}|^2(c_{j}+c_{j-1})+|c_{j}+c_{j+1}|^2(c_{j}+c_{j+1})=0. \end{equation}

  Unlike the previous scenario ( $N=1$), here the nonlinear mixing terms are short range nonlocal, i.e., the nonlinear coupling is between field amplitudes located at sites $j$ and $j\pm1$. An important and immediate question that arises is whether Eq. (6.18) thermalizes to a RJ distribution or not. To answer this, we simulated Eq. (6.18) using a linear initial distribution among $M=100$ supermodes with $\epsilon_j$ given by Eq. (2.8). The results can be summarized in Fig. 9(a). Clearly, a highly localized supermode with two peaks leads to a highly localized mixing tensor, which prevents the system from thermalizing.

• • $\underline{N=3}$: This is an interesting case where one starts to observe a change in the modal statistical distribution. Similar to the previous analysis one arrives at the coupled system:

  (6.19)\begin{equation} i\frac{dc_j}{dz}+\epsilon_jc_j+|C_j|^2C_j+|C_{j-1}|^2C_{j-1}+|C_{j+1}|^2C_{j+1}=0, \end{equation}
  

  Figure 9. Numerical simulation of Eq. (6.13) for various values of $N$ with a linear initial distribution among $M=100$ supermodes and eigenvalues $\epsilon_j$ given by Eq. (2.8). The simulation results for the modal occupancies are shown at $z=120000$ and are averaged over 800 ensembles of random phase initial conditions. (a) For $N=2$ and $P=2$ the averaged modal occupancies remain nearly unchanged from their initial distribution. (b) For $N=3$ and $P=2$ an almost exact Rayleigh–Jeans distribution (blue solid line) is attained, although some discrepancies appear for higher-order modes in the range $-2 \lt \epsilon_j \lt 0$. (c) For $N=4$ and $P=2$, we also get a good match between theory (red solid line) and simulation (black stars) and (d) for $N=5$ and $P=1$ the modal occupancies (black stars) match the theoretically predicted RJ distribution (solid yellow line) nearly perfectly.

  where $C_{j}\equiv c_{j+1}+c_{j}+c_{j-1}$. It is now evident that the degree of nonlocality has increased (thus inducing a long-range nonlinear wave mixing process) relative to the previous cases. That is, the nonlinear coupling goes beyond nearest neighbors and includes optical fields located at sites $j\pm 2$ (on top of $j$, $j\pm1$). As a result, one might expect a fundamental difference in the statistical behavior of the projection coefficients. We have numerically solved Eq. (6.19) under the same computational conditions as in the $N=2$ case and despite the sparsity of the nonlocal mixing terms the statistical distribution of the modal amplitudes $\langle|c_j|^2\rangle$ almost approached the theoretically predicted RJ distribution as is seen in Fig. 9(b). To further clarify the role that nonlinear nonlocal wave mixing plays in accelerating the process of reaching the thermalization state, we have simulated Eq. (6.13) for $N=4$ and $N=5$. The results were in good agreement with the theoretically predicted RJ distribution (see Fig. 9(c) and (d)) .