Thermalization dynamics of macroscopic weakly nonintegrable maps

We take inspiration from quantum Unitary Circuits (UCs) systems. Unitary Circuits serve as basis models for quantum computing [34, 35] and as a platform for studies of quantum chaos [36] and operator spreading [37]. Typically a state in quantum UCs is represented in terms of a one dimensional wave function with multiple components per site (usually in terms of qubits) which are then coupled by unitary operations.

Our classical version of Unitary Circuits is a map acting on a state represented by a one dimensional complex valued vector of size $N$ consisting of $N/2$ unit cells with sites $A$ and $B$:

The time evolution of the system is governed by a discrete unitary map consisting of several transformations of $\vec{\Psi}$:

where maps $\hat{C}_{A,B}$ and $\hat{C}_{B,A}$ are given by unitary matrices acting on the neighboring sites $(\psi_{n}^{A},\psi_{n}^{B})^{T}$:

and $\hat{G}_{n}$ induces a nonlinearity proportional to the strength coefficient $g$:

The local transformations $\hat{C}$ can be in general represented as arbitrary $2\times 2$ unitary matrices. Our particular choice is parametrized with a single angle $\theta$ which plays the role of a hopping parameter strength in Hamiltonian systems. At the same time the nonlinearity $g|\psi_{n}|^{2}$ is an analog of an effective mean-field potential. The maps of this type can be constructed experimentally with periodic application of magnetic field pulses [38, 39, 40, 41].

The evolution (Eq. (2)) results in the following equations of motion:

where $\varphi^{A,B}_{n}(t)$ are the components of state vector $\vec{\Psi}$ after the application of mixing maps $\hat{C}$ (see Fig. 1).

We consider two integrable limits in the parameter space $\{g,\;\theta\}$ of the classical Unitary Circuits map: vanishing coupling: $\{\theta=0,\;g\neq 0\}$, and vanishing nonlinearity: $\{g=0,\;\theta\neq 0\}$ (see Fig. 1). In each case the integrable limit is characterized by a corresponding set of $N$ observables (actions) which are decoupled and thus are conserved in time during the evolution. In what follows we will investigate each of the integrable limits in detail by deviating from it slightly and thus inducing a coupling between the actions.

Our goal is to study the ergodization of local and extended observables in two distinct integrable limits. Ergodization is a process of time averages of observables approaching their phase space averages. In line with this definition we will define sets of finite time averages of observables and study the statistical properties of these sets, in particular their variance. In both SRN and LRN regimes we study the statistical properties of sets of finite time averages of observables. For a time-dependent observable $o(t)$ we define a finite time average as:

We construct a set of finite time averages $\{\overline{o}_{T}\}_{M}$ by following $M$ trajectories. This set is characterized by a distribution with probability density function $\rho(\bar{o}_{T})$. From the ergodization hypothesis we expect each $\bar{o}_{T\rightarrow\infty}=\langle o\rangle$, where $\langle o\rangle$ is an average taken over the phase space. The distribution $\rho$ is therefore expected to peak around the phase space average, reducing its variance to zero and thus approaching a delta-function for infinite averaging times: $\rho(\bar{o}_{T\rightarrow\infty})=\delta(\bar{o}_{T\rightarrow\infty}-\langle o\rangle)$. We study the convergence by following the variance $\sigma^{2}(T)$ of the distribution $\rho(\bar{o}_{T})$. We perform the analysis for two distinct observables – local observables $|\psi_{n}|$, and extended observables $|c_{k}^{\alpha}|$ which are the amplitudes of normal modes coefficients.

To characterize chaoticity related timescales of non-integrable dynamics we compute the full spectrum of Lyapunov characteristic exponents (LCEs). An advantage of LCE spectra is that they are independent of any coordinate basis choice. Nonzero Lyapunov spectra necessarily imply chaotic mixing dynamics (while the converse is not necessarily the case). On the downside we note that while the largest LCE can be experimentally assessed, the entire Lyapunov spectrum appears to be not experimentally measurable.

We compute the full LCE spectra using methods described in [42] in both SRN and LRN limits upon variation of $\theta$ and $g$ respectively. We evolve a set of deviation vectors in the phase space and extract the Lyapunov exponents corresponding to their growth, see Appendix A for more details. The number of Lyapunov characteristic exponents equals the dimensionality of the phase space, in our case $4N$. Lyapunov exponents come ordered from largest to smallest value upon incrementing the index $i$. Due to the symplectic nature of the unitary circuits map the spectrum is symmetric with LCEs coming in pairs $\Lambda_{i}=-\Lambda_{2N-i+1}$. Norm conservation ensures two vanishing LCEs $\Lambda_{N}=\Lambda_{N+1}=0$.

Without loss of generality we consider only positive LCEs. We renormalize Lyapunov spectra $\overline{\Lambda}_{i}=\Lambda_{i}/\Lambda_{max}$ and rescale the index index $\rho=i/2N$ so that all positive LCEs $\overline{\Lambda}(\rho)$ correspond to $\rho\in[0,1]$.

An in depth investigation of Lyapunov spectra scaling of weakly nonintegrable dynamics has been performed in Ref.  [1].

Further, we perform a characterization of the tangent vector corresponding to the largest LCE (for definition of tangent vectors see Appendix A). At any given time the normalized tangent vector points in the direction of the strongest chaotization. We compute the time average of the participation number $\overline{PN}$ of the normalized tangent vector $\vec{w}(t)$:

The tangent vectors depend on the coordinate basis choice. Once the basis is fixed, a delocalized tangent vector results in $PN\sim N$ while a localized tangent vector yields a system size independent value of $PN$.

We perform the evolution of the state vector $\vec{\Psi}(t)$ in coordinate space using the unitary map defined in Eq.(2). Before measuring the observables the system is prerun to ensure thermalization. The initial conditions for each component of the state vector are set as $\psi_{n}^{A,B}=r_{n}e^{i\gamma_{n}}$. For each of the $M$ initial conditions (trajectories) we generate the amplitudes $r_{n}$ as uncorrelated random numbers to be distributed according to an exponential distribution with probability density function $f(x)\sim xe^{-x^{2}}$ in accord with the Gibbs distribution for the norm densities. We then renormalize the state vector such that the norm density is set to unity. The phases $\gamma_{n}$ are distributed uniformly on the interval $[0,\;2\pi]$.

LOs: Due to the translational invariance in the system, we assume the local observables (LOs) to be statistically identical and independent. This leads to the possibility of obtaining a set of finite time averages of LOs from a single trajectory. This way an output of a single run will generate the variance $\sigma^{2}(T)$ of a set of finite time averages of $N$ observables.

EOs: In contrast to LOs the normal mode coefficients $c_{k}^{\alpha}$ - which are the extended observables (EOs) - are not statistically identical. Thus we choose a specific value of the wave vector $k=\pi/2$ and perform $M$ trajectory runs. We extract the time average $\overline{c}_{k}^{m}(T)$ to obtain the set of $M$ finite time averages $\{\overline{c}_{k}^{m}(T)\}_{m=1}^{M}$ whose variance $\sigma^{2}(T)$ is then computed.

For the sake of clarity we briefly list and define the time scales involved in our studies. The ergodization time $T_{E}$ is the time scale up to which the actions which turn integrals of motion at the very integrable limit stay essentially constant for a weakly nonintegrable system. It follows that $T_{E}$ must diverge upon approaching the very integrable limit. We study systems with short range coupling in real space, which co-exists with any of the nonintegrable network ranges (short and long), as the latter are defined in the corresponding action space of actions which turn integrals of motion at the very integrable limit. Therefore, regardless of the type of nonintegrable network, once the time scale $T>T_{E}$ is reached, diffusion in real space will set in. A second diffusion time $T_{D}\sim N^{2}$ marks the time at which the diffusion in real space reaches the boundaries of the finite system. Note that for an infinite system size $T_{D}$ diverges for any weakly nonintegrable system which still has a finite (though potentially very large) ergodization time $T_{E}$. Finally, we will compare the above time scales with the Lyapunov time $T_{\Lambda}=1/\Lambda_{max}$ which is given by the inverse of the largest Lyapunov exponent $\Lambda_{max}$. In all cases $T_{\Lambda}\leq T_{E}$ which is an expected result - there can be no ergodization and thermalization before any chaos sets in. However, we find that the ratio $T_{E}/T_{\Lambda}$ diverges much faster with increasing $T_{E}$ for SRNs, while stays almost constant for LRNs.

In the integrable limit $\theta=0$ the LOs are decoupled and the variance of the set of finite time averages of LOs will stay constant over time $\sigma^{2}(T)=const$. Once the small deviation $\theta\neq 0$ has been introduced the LOs are coupled and the dynamics shows nonintegrable behavior. The weak coupling of LOs will manifest in nearly frozen actions with rare resonant spots in the system, where chaotic dynamics takes place [24]. The mean distance between the chaotic spots grows upon approaching the integrable limit. The strength of chaotic dynamics (largest Lyapunov exponent) diminishes upon approaching the integrable limit. This interplay between chaotic and non-chaotic parts of the system will result in the ergodization time $T_{E}$ - a time scale on which the chaotic spots diffuse over a distance of the order of the average spacing between the resonances. It follows that $\sigma^{2}(T<T_{E})$ stays approximately constant up to $T_{E}$. That time scale $T_{E}$ will diverge upon approaching the integrable limit. At finite distance from the integrable limit the resonances continue to diffuse through the system resulting in $\sigma^{2}(T)\sim T^{-1/2}$ for $T_{E}<T<T_{D}$ [25]. Once the excitations diffuse across the entire system all correlations vanish and we expect $\sigma^{2}(T)\sim T^{-1}$ due to finite size effects. The time scale of the transition from $\sigma^{2}(T)\sim T^{-1/2}$ to $\sigma^{2}(T)\sim T^{-1}$ is denoted as $T_{D}\sim N^{2}$ [25]. The schematic representation of expected behavior of $\sigma^{2}(T)$ is presented in Fig. 2.

The EOs are not conserved even at the very integrable limit. They will show fast fluctuations and quick pseudo-thermalization in analogy with Ref. [21]. Thus we expect an immediate $\sigma^{2}(T)\sim T^{-1}$ decay starting from the shortest time scales, see Fig. 2.

In the integrable limit $g=0$ the system dynamics is linear, and we expect $\sigma^{2}(T)=const$ for EOs. Upon the deviation from the limit the variance is expected to decay $\sigma(T)^{2}\sim T^{-1}$ after some ergodization time $T_{E}$ required for spread of chaos into the network (see Fig. 3). We expect time $T_{E}$ to be of the order of Lyapunov time $T_{\Lambda}$ as there appears to be no other time scale governing the dynamics, see Fig. 3.

LOs will show a more involved outcome. First, even at the integrable limit they are not conserved, and will show fast fluctuation and quick pseudo-thermalization in analogy with Ref. [21]. Thus we expect an immediate $\sigma^{2}(T)\sim T^{-1}$ decay starting from the shortest time scales. However, this only holds up to $T_{E}$ if the system is large enough. Indeed, the EOs are preserved up to $T_{E}$ and correspond to ballistically propagating modes (waves) in real space having a largest finite group velocity $v_{g}$. If the system size $N\gg v_{g}T_{E}$ the modes will start to interact and ballistic propagation is replaced by diffusive propagation in real space for $T>T_{E}$. Correspondingly the LO dynamics results in a crossover from $\sigma^{2}(T)\sim T^{-1}$ to $\sigma^{2}(T)\sim T^{-1/2}$ for $T>T_{E}$. At a time scale $T_{D}\sim N^{2}$ the diffusion reaches the system boundaries, and the LO dynamics crosses over back to a final asymptotic $\sigma^{2}(T)\sim T^{-1}$ decay. We show the expected behavior in Fig. 3.

In the SRN case we fix the nonlinearity strength $g=1$ without loss of generality. The system size $N=10^{4}$. At $\theta=0$ the LOs are frozen – the system is at the integrable limit. Upon increasing the parameter $\theta$, the LOs get intercoupled with nearest neighbors into a SRN. We study the statistics of local and extended observables in the SRN in close proximity to the corresponding integrable limit.

In the LRN case we fix $\theta=0.33\pi$ without loss of generality. The system size $N=10^{5}$. For $g=0$ the normal mode coefficients $|c_{k}|$ are constant in time. The EOs are frozen, indicating another integrable limit. Upon increasing $g\neq 0$ the EOs get intercoupled with an all-to-all coupling into a LRN. We study the statistics of local and extended observables in the LRN in close proximity to the corresponding integrable limit.

In both SRN and LRN cases we compare the ergodization times $T_{E}$ with the characteristic Lyapunov time scales $T_{\Lambda}$ of chaotic dynamics. As the variance of finite time averages of actions $\sigma^{2}(T)$ is expected to transition from nearly constant behavior to $\sigma^{2}(T)\sim T^{-1/2}$ in SRN and $\sigma^{2}(T)\sim T^{-1}$ in LRN we follow the local derivatives $d_{\sigma^{2}}$ of the variance curves $\sigma^{2}(T)$ , see insets of Fig. 4 and Fig. 5. The ergodization time $T_{E}$ is extracted as the time when $d_{\sigma^{2}}=-0.25$ in SRN and $d_{\sigma^{2}}=-0.75$ in LRN for the first time. The Lyapunov times $T_{\Lambda}$ are computed for the same parameter values and plotted against $T_{E}$ in Fig. 6 for SRN and Fig. 7 for LRN. In the SRN we notice $T_{E}\gg T_{\Lambda}$ such that the ratio $T_{E}/T_{\Lambda}$ grows quickly as the system approaches integrable limit $\theta\rightarrow 0$. In the LRN case $T_{\Lambda}\sim T_{E}$ and their ratio is practically constant upon approaching the integrable limit $g\rightarrow 0$.

In Fig. 8 and Fig. 9 we show the rescaled Lyapunov spectra $\bar{\Lambda}(\rho)=\Lambda(\rho)/\Lambda_{max}$ corresponding to the SRN and the LRN respectively.

The rescaled spectrum tends to a non-analytic function in the SRN upon approaching the integrable limit in Fig. 8. Close to the limit the spectrum shows exponential decay with a (length scale) exponent which diverges at the very integrable limit [1]. The slowing down of thermalization is characterized by two diverging scales - a time scale and a length scale.

The rescaled spectrum tends to an analytic function in the LRN upon approaching the integrable limit in Fig. 9. The slowing down of thermalization is characterized by one diverging time scale only.

In the insets of Fig. 8 and Fig. 9 we plot the participation $PN$ number of the tangent vector corresponding to the largest Lyapunov exponent as a function of system size for. For the SRN case the parameters are $g=1,\;\theta=0.001$. For the LRN case the parameters are $\theta=0.33\pi,\;g=0.001$. The participation numbers depend on the basis choice of the tangent vector. Thus we compute the $PN$ in both LO and EO representations.

In the SRN case (inset in Fig. 8) the tangent vectors in LO representation are localized due to rare chaotic resonances in real space. Therefore the $PN$ number in LO representation is predicted to be system size independent. Since a localized distribution in real space turns delocalized in reciprocal space, the EO representation results in a $PN$ which grows and scales with the system size.

In the LRN case (inset in Fig. 9) resonances are expected to appear almost everywhere in the EO representation.Thus we expect the $PN$ number in EO representation to grow and scale with the system size. This is in accord with our numerical observation. Interestingly the delocalized nature of the tangent vector in the EO representation does not imply that the LO representation will result in localized structures. Indeed the numerical computation show that the $PN$ number in LO representation also grows and scales with the system size.