Statistical mechanics of general discrete nonlinear Schrödinger models: Localization transition and its relevance for Klein-Gordon lattices

We extend earlier work [Phys. Rev. Lett. 84, 3740 (2000)] on the statistical mechanics of the cubic one-dimensional discrete nonlinear Schrödinger (DNLS) equation to a more general class of models, including higher dimensionalities and nonlinearities of arbitrary degree. These extensions are physically motivated by the desire to describe situations with an excitation threshold for creation of localized excitations, as well as by recent work suggesting noncubic DNLS models to describe Bose-Einstein condensates in deep optical lattices, taking into account the effective condensate dimensionality. Considering ensembles of initial conditions with given values of the two conserved quantities, norm and Hamiltonian, we calculate analytically the boundary of the “normal” Gibbsian regime corresponding to infinite temperature, and perform numerical simulations to illuminate the nature of the localization dynamics outside this regime for various cases. Furthermore, we show quantitatively how this DNLS localization transition manifests itself for small-amplitude oscillations in generic Klein-Gordon lattices of weakly coupled anharmonic oscillators (in which energy is the only conserved quantity), and determine conditions for the existence of persistent energy localization over large time scales.

Figure 1

Maximum value of norm density $a$ for statistical localization to occur, according to Eq. (12), from initial condition being a traveling wave of wave vector $q$, for 1D DNLS with various $\sigma$. From top to bottom at $q=0$: $\sigma =2∕5,1∕2,2∕3,1,2,3$. Inset is a blow-up for small $a$.

Figure 2

Numerically obtained distribution functions $p\left(A_m\right)$ resulting from long-time integration $\left(t=1.1\times {10}^6\right)$ of (unstable) constant-amplitude initial states ${\psi }_m\left(0\right)=\sqrt{a}$, for a 1D DNLS chain with $N=\mathrm{10\hspace{0.17em}000}$ $\left(C=1\right)$. (a) $\sigma =0.4$, $a=100$ (squares), $a=455$ (triangles), and $a=1000$ (circles). (b) $\sigma =3$, $a=0.6$ (circles), and $a=0.8$ (squares). The dashed lines represent best fits to Eq. (10).

Figure 3

(a) Time-averaged (non-normalized) distribution functions $p\left(A_m\right)$ for weakly perturbed $\pi ∕2$ type-$\mathrm{H}$ SWs with $a=1$ $(\dots ,-\sqrt{2},0,\sqrt{2},0,\dots )$ as initial conditions $\left(C=1\right)$. $(*)$: $\sigma =\frac{2}{3}$; $(+)$: $\sigma =1$; $(\times )$: $\sigma =3$. The points are obtained by averaging over 96 time instants for $\mathrm{10\hspace{0.17em}000}⩽t⩽\mathrm{200\hspace{0.17em}000}$ and $N=\mathrm{10\hspace{0.17em}000}$ $\left(\sigma =\frac{2}{3}\right)$, 177 time instants for $\mathrm{500\hspace{0.17em}000}⩽t⩽\mathrm{1\hspace{0.17em}380}000$ and $N=1000$ $\left(\sigma =1\right)$, and 91 time instants for $\mathrm{10\hspace{0.17em}000}⩽t⩽\mathrm{55\hspace{0.17em}000}$ and $N=1000$ $\left(\sigma =3\right)$. Straight lines for $\sigma =\frac{2}{3}$ and $\sigma =1$ are predictions from Eq. (10) with $\beta =0$ and $\gamma =1∕a=1$, while curve for $\sigma =3$ is a prediction from Eq. (10) with fitted values of $\beta =0.042$ and $\gamma =0.65$. (b) Average (over space and time) of the coupling part of $h$ [i.e., $⟨2C\sqrt{A_m{A}_{m+1}}\cos ({\varphi }_m-{\varphi }_{m+1})⟩$] vs time for the simulations in (a) with, from top to bottom at $t=\mathrm{10\hspace{0.17em}000}$, $\sigma =\frac{2}{3}$, $\sigma =1$, and $\sigma =3$, respectively. Magnification in the inset illustrates the slow long-time decrease for $\sigma =\frac{2}{3}$.

Figure 4

(Color online) Resulting distribution functions $p\left(A_m\right)$ after long-time integrations of initial conditions consisting of weakly perturbed 2D constant-amplitude $\left(q_x=q_y=0\right)$ unstable solutions with (a) $a=5$, (b) $a=7$, and (c) $a=8$. Curves in (a)–(c) have been obtained by averaging over a number of different random initial perturbations [8 in (a) c 14 in (b), and 100 in (c)]; system size $128\times 128$ $\left(N=\mathrm{16\hspace{0.17em}384}\right)$; integration times $t=\mathrm{500\hspace{0.17em}000}$ (a), $\mathrm{200\hspace{0.17em}000}$ (b), respectively $\mathrm{50\hspace{0.17em}000}$ (c). Curves in (d) have been obtained from one single realization for a $50\times 50$ system with $a=7$, by averaging over 20 different time instants in the intervals $500<t<\mathrm{10\hspace{0.17em}000}$ (squares) and $\mathrm{110\hspace{0.17em}000}<t<\mathrm{300\hspace{0.17em}000}$ (circles), respectively. The scales are such that, in (a)–(c), the dots with smallest probability correspond to one site in one realization, and in (d) to one site at one time instant. Straight lines are predictions according to Eq. (10) with $\beta =0$. ($\sigma =C=1$.)

Figure 5

(Color) (a) Distribution functions $p\left(A_m\right)$ after integrations over long but finite times $\left(t=5\times {10}^6\right)$ of initial conditions consisting of weakly perturbed 3D constant-amplitude $\left(q_x=q_y=q_z=0\right)$ unstable solutions with $a=9$ (blue), $a=12$ (red), and $a=15$ (black). System size $64\times 64\times 64$ $\left(N=\mathrm{262\hspace{0.17em}144}\right)$, $\sigma =C=1$. (b) Intensities, in a representative $15\times 15\times 15$ subbox of the simulation box, at the end of the simulation for $a=9$ in (a). Red and yellow patches are localized breathers.

Figure 6

(Color online) Numerical integration of KG Morse chain with $C_K=0.005$, $N=200$, and randomly perturbed constant-amplitude initial condition $u_n\left(0\right)=0.05$. (a) Time evolution of local energy density (note logarithmic time scale). (b) Main figure: $\mathcal{H}∕N$ vs $\mathcal{A}∕N$ for the simulation in (a), with $\mathcal{A}$ calculated from Eq. (36) and, from top to bottom in the left part of the figure, $\mathcal{H}$ calculated from Eqs. (39, 37, 42), respectively. Time runs from right to left (i.e., $\mathcal{A}∕N$ decreases). Lowest curve is the localization transition line (9). Inset in (b) shows the ratio of time-averaged cubic (33) to quartic (34) energies vs time, compared to the DNLS prediction (44) (lower line).

Figure 7

Thermalization of a quartic KG chain $\left(\alpha =0,{\beta }^{\prime }=1\right)$ with $C_K=0.01$, $N=800$, coupled to a thermal bath at temperature $T^{\prime }=0.005$ with dissipation constant $\eta =0.1$. (a) Time-averaged total kinetic energy $⟨{\sum }_n\left({\dot{u}}_n^2∕2\right)⟩$. (b) $\mathcal{H}∕N$ vs $\mathcal{A}∕N$ for the simulation in (a), with $\mathcal{A}$ calculated from Eq. (36), and $\mathcal{H}$ calculated from Eq. (37). Each dot represents a time average over the interval $[t-100,t]$ at $\mathrm{15\hspace{0.17em}382}$ different times $t$. Line in (b) is the localization transition line (9). Larger points in (b) show the locations of the initial conditions used in Fig. 8.

Figure 8

(a) Example of a breather appearing in the microcanonical integration of an initial condition represented by the lower large point in Fig. 7b. (b),(c) (Non-normalized) velocity distribution functions $p\left({\dot{u}}_n\right)$ obtained from long-time numerical microcanonical integrations (points) compared to Maxwellian distributions $P\left(\dot{u}\right)\sim {\left(2\pi T^{\prime }\right)}^{-1∕2}\exp (-{\dot{u}}^2∕2T^{\prime })$ (lines) at the estimated temperature. In (b) the initial condition is the same as for (a), the temperature is $T^{\prime }\approx \mathrm{0.004\hspace{0.17em}94}$, and the integration time is $1.2\times {10}^6$. (c) corresponds to the upper large point at $(0.704,-0.471)$ in Fig. 7b, with $T^{\prime }\approx \mathrm{0.004\hspace{0.17em}85}$, and integration time $0.6\times {10}^6$. In both cases, the velocities of all sites are registered in intervals of 0.6 time units. Insets in (b) and (c) show magnification of the small-velocity regime in a nonlogarithmic scale.