Tensorial symmetries and optical lattice thermalization

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 that includes a nearest-neighbor coupling, on-site potential, and a collection of cubic nonlinearities given by

$$i\frac{dA_{n}}{dz}+A_{n-1}+A_{n+1}+V_{n}A_{n}+f(A_{n},A_{n\pm m})=0\,,$$

(1)

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. (1) remains invariant under the gauge transformation:

$$A_{n}\rightarrow A_{n}e^{i\theta}\;,$$

(2)

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

$$f(A_{n},A_{n\pm\tilde{n}})\sim A_{n+n_{1}}A_{n+n_{2}}A^{*}_{n+n_{3}}\,,$$

(3)

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.

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

$$A_{n}(z)=\psi_{n}e^{i\epsilon z}\;,$$

(4)

leads to the following linear matrix eigenvalue problem:

$$\mathbb{M}|\psi^{(j)}\rangle=\epsilon_{j}|\psi^{(j)}\rangle\;,$$

(5)

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, defined by

$$|\psi^{(j)}\rangle\equiv\{\psi_{n}^{(j)}\}_{n=1}^{M}=\begin{pmatrix}\psi^{(j)}_{1}\\ \vskip 0.72229pt\\ \psi^{(j)}_{2}\\ \vdots\\ \psi^{(j)}_{M}\end{pmatrix},$$

(6)

and $\epsilon_{j}$ being the matrix eigenvalues. The supermodes constitute an orthonormal base satisfying the orthogonality condition:

$$\langle\psi^{(j)}|\psi^{(j^{\prime})}\rangle=\delta_{j,j^{\prime}}\;.$$

(7)

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.(5) can be obtained. They are given by

$$\epsilon_{j}=2\cos(\frac{j\pi}{M+1})\,,$$

(8)

$$\psi_{n}^{(j)}=\sqrt{\frac{2}{M+1}}\sin(\frac{jn\pi}{M+1})\,.$$

(9)

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

$$A_{n}(z)=\sum\limits_{j=1}^{M}c_{j}(z)\psi_{n}^{(j)}\,,$$

(10)

where $c_{j}(z)$ are the so-called complex projection coefficients. Substituting the expansion (10) into (1) we obtain the following generic nonlinear evolution equation governing the dynamics of the projection coefficients (or modal amplitudes):

$$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\,,$$

(11)

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. (3) this tensor takes the form

$$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}}\,.$$

(12)

Note that Eq. (11) is invariant under the phase transformation $c_{j}\rightarrow c_{j}e^{i\vartheta}$ with real constant $\vartheta$. If one relaxes condition (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 (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. In this regard, the quantity of interest is defined by

$$\langle|c_{j}|^{2}\rangle\equiv\langle\lim_{z\to\infty}|c_{j}(z)|^{2}\rangle\,,$$

(13)

where $\langle\cdots\rangle$ denotes the ensemble average over many realizations of the random initial conditions $c_{j}(0)$. 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

$$\langle|c_{j}|^{2}\rangle=h(\epsilon_{j})\,,$$

(14)

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

$$h=\frac{a}{\epsilon_{j}+b}\hskip 2.84544pt\,,$$

(15)

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

Throughout the rest of the paper, equation (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:

$$c_{j}(0)\equiv\tilde{c}_{j}e^{2\pi ir_{j}}\,,$$

(16)

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:

  $$\tilde{c}_{j}=\sqrt{\alpha(\epsilon_{j}-\epsilon_{1})},\quad\alpha=\frac{P}{\sum_{i}(\epsilon_{i}-\epsilon_{1})}\,,$$

  (17)

• •

  piecewise uniform distribution across the supermodes:

  $$\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}$$

  (18)

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

$$P=\sum\limits_{j=1}^{M}|c_{j}|^{2}\,.$$

(19)

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

	$\displaystyle\text{Quasi-Hermiticity:}\,\,T_{j,k,l,m}=T^{*}_{l,m,j,k}\,,$		(20)
	$\displaystyle\text{Permutation symmetry:}\,\,T_{j,k,l,m}=T_{m,l,k,j}\,.$		(21)

Symmetry condition (20) guarantees the conservation of power. The second property (21) allows one to derive Eq. (11) from Hamilton’s equations of motion for the conjugate pair of canonical variables $c_{j},\,c_{j}^{*}$:

$$\frac{dc_{j}}{dz}=i\{c_{j},\mathcal{H}\},\;\;\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{\text{NL}}\,,$$

(22)

where

$$\mathcal{H}_{0}=\sum\limits_{j=1}^{M}\epsilon_{j}|c_{j}|^{2}\,,$$

(23)

and

$\displaystyle\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)\,.$

(24)

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

$$\{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)\;,$$

(25)

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.(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.

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. (LABEL:FullLattice) with $g=1$ while setting the rest of the parameters to zero. This leads to

$$i\frac{dA_{n}}{dz}+A_{n-1}+A_{n+1}+|A_{n}|^{2}A_{n}=0\,.$$

(36)

Following similar ideas outlined in Sec. II.3 (see also Appendix B for further details) we arrive at

$$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)}\,,$$

(37)

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

$$\{q_{i}\}_{i=1}^{4}\equiv\{j,k,l,m\},\,\,\,\,\{p_{i}\}_{i=1}^{4}\equiv\{j,m,k,l\}.$$

(38)

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

$$\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}$$

(39)

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

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

(40)

Here $\rho_{\iota}$ denotes the remainder of the fraction $w_{\iota}/(2M+2)$ and $\delta_{0,\rho_{\iota}}$ is the Kronecker delta function. Equation (40) 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. (39)). 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

$$S=\Big{\{}0,\,\frac{\pm 1}{2M+2},\,\frac{1}{M+1},\,\frac{3}{2M+2},\,\frac{2}{M+1}\Big{\}}\,,$$

(41)

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 (11) (for even $M$) reads

$$i\frac{dc_{j}}{dz}+\epsilon_{j}c_{j}+\frac{1}{M+1}\left(\frac{3}{2}\sum_{\begin{subarray}{c}\text{branch}\#\\ 2,3,5,10\end{subarray}}c_{k}c_{l}c_{m}^{*}+\sum_{\begin{subarray}{c}\text{branch}\#4,6,7,\\ 11,12,15\end{subarray}}c_{k}c_{l}c_{m}^{*}+\frac{1}{2}\sum_{\begin{subarray}{c}\text{branch}\#\\ 8,13,16,18\end{subarray}}c_{k}c_{l}c_{m}^{*}-\frac{1}{2}\sum_{\begin{subarray}{c}\text{branch}\#\\ 9,14,17,19\end{subarray}}c_{k}c_{l}c_{m}^{*}\right)=0\,,$$

(42)

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

$$i\frac{dc_{1}}{dz}+\underbrace{(\epsilon_{1}+P-|c_{1}|^{2}/2)}_{\text{non-mixing}}c_{1}+\underbrace{c_{2}^{2}c_{1}^{*}/2}_{\text{mixing}}=0\,,$$

(43)

$$i\frac{dc_{2}}{dz}+\underbrace{(\epsilon_{2}+P-|c_{2}|^{2}/2)}_{\text{non-mixing}}c_{2}+\underbrace{c_{1}^{2}c_{2}^{*}/2}_{\text{mixing}}=0\,,$$

(44)

where again $P=|c_{1}|^{2}+|c_{2}|^{2}$ is the total power (which is a conserved quantity). If one only keeps the non-mixing terms, we find $\langle|c_{j}(z)|^{2}\rangle=|c_{j}(0)|^{2}$ which indicates lack of thermalization. Thus, the middle terms in Eq. (43) and Eq. (44), do not contribute to the wave mixing process. However, the presence of the last 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-mixing (${\cal N}_{j}$) and mixing terms (${\cal M}_{j}$) for which the evolution of $c_{j}$ is governed by

$\displaystyle i\frac{dc_{j}}{dz}$

$\displaystyle+$

$\displaystyle{\cal N}_{j}c_{j}+{\cal M}_{j}(c_{j})=0\,.$

(45)

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

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,\pm 1/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. II.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. (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 [13].

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. VI.

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. (LABEL:FullLattice). The main motivation is the identification of the tensorial symmetries which will be later tied to thermalization.