text:
# ORTHOGONAL POLYNOMIALS IN ANALYTICAL METHOD OF SOLVING DIFFERENTIAL EQUATIONS DESCRIBING DYNAMICS OF MULTILEVEL SYSTEMS11footnote 1 Integral Transforms and Special Functions, 2000, Vol. 10, No. 3–4, pp. 299–308. Received December 12, 1999. © OPA (Overseas Publishers Association) N.V. with permission from Gordon and Breach Publishers.

## 1 INTRODUCTION

Orthogonal polynomials and functions are used as the complete bases of the solutions of differential and difference equations. In quantum physics they are of primary importance in the construction of wave functions which are orthogonal by definition. In present paper orthogonal polynomials and functions are used to describe excitation of quantum multilevel systems in high-power laser field.

The Schrödinger equation for a multilevel system being exciting in a monochromatic laser field looks like

$$ida_n(t)/dt=f_{n+1}\mathrm{exp}\left(i\epsilon _{n+1}t\right)a_{n+1}(t)+f_n\mathrm{exp}\left(i\epsilon _nt\right)a_{n1}(t),n=0,1,2,\mathrm{}$$
(1)
where $`f_n`$ is a dipole moment function for the transitions $`n1n`$ (that is $`\mu _{n1,n}=f_n\mu _{0,1}`$), $`\epsilon _n`$ is a frequency detuning. Probability amplitudes $`a_n(t)`$ allow to determine levels populations $`\rho _n(t)=\left|a_n(t)\right|^2`$, an average number of absorbed photons $`<n>=_{n=1}^{n_{max}}n\rho _n(t)`$ etc.

To obtain the solution of (1) we introduce an analytical method  based on orthogonal functions $`p_n(z)`$.

## 2 FROM ORTHOGONAL POLYNOMIALS TO QUANTUM SYSTEMS

We seek the solution of (1) in form

$$a_n(t)=_A^B\sigma (x)\frac{p_0(x)}{d_0}\frac{p_n(x)}{d_n}\mathrm{exp}\left[it\left(r\theta (x)+s_n\right)\right]𝑑x,$$
(2)
where $`d_n`$ is the squared norm and $`s_n`$ determines the frequency detuning: $`\epsilon _n=s_ns_{n1}`$.

At $`t=0`$ the solution (2) passes into the orthogonality relation for the functions $`p_n(x)`$. Hence we obtain the initial condition for Eq. (1): $`a_n(t)=\delta _{n,0}`$. The functions $`p_n(x)`$ should satisfy the recurrence relation

$$\left(r\theta (x)+s_n\right)\frac{p_n(x)}{d_n}=f_n\frac{p_{n1}(x)}{d_{n1}}+f_{n+1}\frac{p_{n+1}(x)}{d_{n+1}}.$$
(3)
Initially we select the sequence of orthogonal functions which determines the quantum system to be investigated. The takeoff criterion is the relevant recurrence relation which allows to find both the constant $`r`$ and the functions $`f_n`$ and $`s_n`$. After that we obtain the analytical solution of (1) using the integral transform like (2).

The application of analytical method (13) makes possible describing dynamics of manifold quantum models including: $`N`$-level and infinite-level, equidistant and non-equidistant, excited resonantly and nonresonantly, closed and nonclosed, excited by monochromatic field and laser pulse. Gegenbauer, Jacobi, Pollaczek, Laguerre, Meixner, Krawtchouk, Hahn polynomials; Szegö-Jacobi polynomials built using the Szegö method; Christoffel-Legendre orthogonal polynomials built by using the Christoffel formula; orthogonal Gauss hypergeometric functions are used by us to obtain the exact solutions for these systems.

This quantum systems manifold includes the systems with various characteristics, first of all with qualitatively different dipole moment functions $`f_n`$. In particular Jacobi, Legendre and Chebyshev of the $`1^{\text{st}}`$ kind polynomials lead to the quantum systems with $`f_n1`$, Jacobi, Gegenbauer and Chebyshev of the $`2^{\text{nd}}`$ kind polynomials — to the systems with $`1f_n<\sqrt{n}`$, Pollaczek, Laguerre and Meixner polynomials give rise to two sets of the systems with $`\sqrt{n}<f_nn`$ and $`f_n>n`$. It is known that Hermite polynomials correspond to the harmonic oscillator ($`f_n=\sqrt{n}`$).

Thus the existing diversity of orthogonal polynomials makes possible to obtain the solutions of (1) describing the dynamics of various quantum systems with different properties. Each orthogonal polynomials sequence (OPS) can be put in unique correspondence with some quantum system. Below we shall discuss some examples.

## 3 JACOBI POLYNOMIALS

Let’s consider the quantum systems which have been built by means of Jacobi polynomials. The dipole moment function following from a recurrence relation (3) looks like

$`f_n`$ $`=`$ $`{\displaystyle \frac{2r}{2n+\alpha +\beta }}\left[{\displaystyle \frac{n\left(n+\alpha \right)\left(n+\beta \right)\left(n+\alpha +\beta \right)}{\left(2n+\alpha +\beta 1\right)\left(2n+\alpha +\beta +1\right)}}\right]^{{}_{}{}^{1}/_{2}^{}},`$ (4)
$`r`$ $`=`$ $`{\displaystyle \frac{\alpha +\beta +2}{2}}\left[{\displaystyle \frac{\alpha +\beta +3}{\left(\alpha +1\right)\left(\beta +1\right)}}\right]^{{}_{}{}^{1}/_{2}^{}},\alpha >1,\beta >1.`$ (5)
Thus the dipole moments of higher transitions tend to the constant value $`r/2`$. Energy levels of the Jacobi systems are nonequidistant but the frequency detuning decreases with $`n`$ approximately as $`n^3`$, i.e. the highest levels are practically equidistant and they are being excited resonantly.

Jacobi OPS contains two parameters $`\alpha `$ and $`\beta `$ which determine the energy levels positions. Therefore we obtain two-parameter continual family of multilevel quantum systems. Each pair of $`(\alpha ,\beta )`$ corresponds to the specific Jacobi quantum system which is described by dipole moment $`f_n`$ and frequency detuning $`\epsilon _n`$ functions.

The solution of (1) for the Jacobi systems dynamics is expressed through the Kummer hypergeometric functions:

$`a_n(t)=\left(2irt\right)^n\left[{\displaystyle \frac{\left(\alpha +1\right)_n\left(\beta +1\right)_n}{\left(\alpha +\beta +2\right)_{2n}\left(n+\alpha +\beta +1\right)_nn!}}\right]^{{}_{}{}^{1}/_{2}^{}}\times `$
(6)
$`e_1^{it\left(r+s_n\right)}F_1(n+\alpha +1;2n+\alpha +\beta +2;2irt).`$
At $`\alpha =\beta `$ Jacobi polynomials $`P^{(\alpha ,\alpha )}(x)`$ pass into Gegenbauer ones $`C^{\left(\lambda \right)}(x)`$, $`\lambda =\alpha +1/2`$. Thus Kummer function in (3) is reduced to a cylindrical function $`J_{n+\lambda }(rt)`$. The case $`\alpha =\beta =1`$ (Chebyshev polynomials of the $`2^{\text{nd}}`$ kind) gives so called equal-Rabi system $`f_n1`$, for $`\alpha =\beta =0`$ (Chebyshev polynomials of the $`1^{\text{st}}`$ kind) we obtain the system with $`f_1=1,`$$`f_n=1/\sqrt{2}`$, $`n2`$. The levels of both Gegenbauer and Chebyshev systems are equidistant, the excitation is resonant.

An interesting situation arises in a case $`\alpha =\beta `$ when

$$f_n=\left[3\left(n^2\alpha ^2\right)/\left(1\alpha ^2\right)\left(4n^21\right)\right]^{{}_{}{}^{1}/_{2}^{}},s_n=\alpha \left[3/\left(1\alpha ^2\right)\right]^{{}_{}{}^{1}/_{2}^{}}\delta _{n,0}.$$

Here the frequency detuning $`\epsilon _n=s_0\delta _{n,1}`$, the levels are equidistant everywhere except for the first transition and only this transition is nonresonant. In particular at $`\alpha =\beta =\pm 1/2`$ we obtain $`f_n1`$, $`\epsilon _n=\delta _{n,1}`$ and the expressions for levels populations $`\rho _n(t)=\left|a_n(t)\right|^2`$ reshape an especially simple form $`\rho _n\left(t\right)=J_n^2\left(2t\right)+J_{n+1}^2\left(2t\right)`$.

## 4 KRAWTCHOUK POLYNOMIALS

Using Krawtchouk polynomials of discrete variable  we can describe the nonresonant excitation of the quantum systems with finite number of equidistant levels and dipole moment function

$$f_n=\left[\frac{n\left(Nn+1\right)}{N}\right]^{{}_{}{}^{1}/_{2}^{}},n=\overline{0,N},N=1,2,\mathrm{}.$$
(7)
The analytical solution of (7) is expressed  by way of elementary functions:

$`a_n(t)`$ $`=`$ $`\left[\left({\displaystyle \genfrac{}{}{0pt}{}{N}{n}}\right)y^n(t)\right]^{{}_{}{}^{1}/_{2}^{}}\left\{\left[1(1+\sigma )y(t)\right]^{{}_{}{}^{1}/_{2}^{}}+i\left[\sigma y(t)\right]^{{}_{}{}^{1}/_{2}^{}}\right\}^{Nn}\mathrm{exp}\left[it\epsilon (nN/2)\right];`$ (8)
$`\rho _n(t)`$ $`=`$ $`\left({\displaystyle \genfrac{}{}{0pt}{}{N}{n}}\right)y^n(t)\left[1y(t)\right]^{Nn};`$ (9)
$`<n>`$ $`=`$ $`Ny(t)`$ (10)
where

$$y\left(t\right)=\frac{\mathrm{sin}^2\left\{\left[(1+\sigma )/N\right]^{{}_{}{}^{1}/_{2}^{}}t\right\}}{(1+\sigma )},\sigma =N\epsilon ^2/4.$$
(11)
The Krawtchouk systems are suitable exactly solvable models for vibrational excitation of a molecule in high-power infrared laser fields. The solution describes an excitation with arbitrary detuning of laser frequency from the system transitions frequency. The case $`p=^1/_2`$ corresponds to the resonant excitation. Krawtchouk systems family includes both two-level and three-level systems and harmonic oscillator as special and limit cases ($`N=1`$, $`N=2`$ and $`N\mathrm{}`$ correspondingly).

Thus the Krawtchouk polynomials give a possibility to build one-parametric family of quantum multilevel systems and to obtain exact analytical solutions of the equations describing the excitation dynamics. The Krawtchouk systems family is an extensive class of multilevel models for various dynamical processes in spectroscopy, nonlinear optics and the other fields.

The method under consideration gives also the possibility of solving the equations which are more complicated then (1). Let the quantum system consist of finite number of degenerate levels. The equation describing the resonant excitation of such a system is

$`i{\displaystyle \frac{da_{n,m}(t)}{dt}}=\mathrm{\Omega }\left[f_{n+1}a_{n+1,m}(t)+f_na_{n1,m}(t)\right]+`$ (12)
$`\mathrm{\Omega }^{}\left[f_{n+1}g_{m+1}a_{n+1,m+1}(t)+f_{n+1}g_ma_{n+1,m1}(t)+f_ng_{m+1}a_{n1,m+1}(t)+f_ng_ma_{n1,m1}(t)\right]`$
where $`\mathrm{\Omega }`$ and $`\mathrm{\Omega }^{}`$ are Rabi frequencies, $`f_n=\left[n(Nn+1)/N\right]^{{}_{}{}^{1}/_{2}^{}}`$ and $`g_m=\left[m(Mm+1)/M\right]^{{}_{}{}^{1}/_{2}^{}}`$ are dipole moment functions for the interlevel and intersublevel transitions correspondingly, $`n=\overline{0,N}`$, $`m=\overline{0,M}`$; $`M,N=1,2,\mathrm{}`$ Here we assume $`a_{n,m}(0)=\delta _{n,0}\delta _{m,0}`$ as an initial condition and $`\mathrm{\Delta }n=\pm 1`$, $`\mathrm{\Delta }m=0,\pm 1`$ as the selection rules.

The solution of (12) is expressed via Krawtchouk polynomials and elementary functions:

$$a_{n,m}=\frac{i^n}{2^{Mm}}\left[\frac{\left(\genfrac{}{}{0pt}{}{N}{n}\right)}{\left(\genfrac{}{}{0pt}{}{M}{m}\right)}\right]^{{}_{}{}^{1}/_{2}^{}}\underset{j=0}{\overset{M}{}}k_m^{(^1/_2)}(j,M)\left(\genfrac{}{}{0pt}{}{M}{j}\right)\mathrm{sin}^n\tau _j\mathrm{cos}^{Nn}\tau _j$$
(13)
where $`\tau _j=(t/{}_{\sqrt{N}}{}^{})[\mathrm{\Omega }+\mathrm{\Omega }^{}(2jM)/\sqrt{M}]`$.

Selecting other orthogonal polynomials we can obtain the solutions of Eq. (12) for different types of $`f_n`$ and $`g_m`$.

## 5 CHRISTOFFEL-LEGENDRE POLYNOMIALS

Let’s consider the solution of (1) using Christoffel-Legendre orthogonal polynomials

$$p_n(x)=\frac{b}{bx}\left[P_n(x)\kappa _n^1(b)P_{n+1}(x)\right],\kappa _n(b)=P_{n+1}(b)/P_n(b)$$
(14)
built with Christoffel formula  using Legendre polynomials $`P_n(x)`$. The polynomials $`p_n(x)`$ are orthogonal on an interval $`[1,1]`$ with respect to the weight

$$\sigma (x)=(bx)/b,b1$$
(15)
and squared norm

$$d_n=\left\{\left[2/\left(n+1\right)\right]b/\kappa _n\left(b\right)\right\}^{{}_{}{}^{1}/_{2}^{}}$$
(16)
They satisfy the recurrence relation (3) ($`\theta (x)x`$) where

$`f_n`$ $`=`$ $`r\left(n/(2n+1)\right)d_{n1}/d_n,`$
$`r`$ $`=`$ $`3d_1/d_0,`$ (17)
$`s_n`$ $`=`$ $`r\left({\displaystyle \frac{n+1}{2n+1}}\kappa _n(b){\displaystyle \frac{n+2}{2n+3}}\kappa _{n+1}(b)\right).`$
The application of Christoffel-Legendre polynomials in (2) gives  the analytical solution

$$a_n(t,b)=i^n\left(\frac{\pi }{2}\right)^{\frac{1}{2}}\frac{d_0}{d_n}\frac{\mathrm{exp}\left(is_nt\right)}{(rt)^{{}_{}{}^{1}/_{2}^{}}}\left[J_{n+0,5}\left(rt\right)i\kappa _n^1\left(b\right)J_{n+1,5}\left(rt\right)\right]$$
(18)
of (1). This solution describes the excitation of one-parameter family of multilevel Christoffel-Legendre systems with nonequidistant energy levels, i.e. $`\epsilon _n\text{const}`$.

## 6 LEGENDRE FUNCTIONS OF THE FIRST KIND

Another quantum system can be built by using the functions

$$p_n^\lambda (x)=(i)^\mu \mathrm{\Gamma }(\mu +1)P_\lambda ^{\lambda \mu }(x)|_{\mu =n},$$
(19)
being a special case of Legendre function of the $`1^{\text{st}}`$ kind $`P_\nu ^\mu (x)`$. The functions $`p_n^\lambda (x)`$ satisfy the recurrence relation (3) where

$`f_n`$ $`=`$ $`\left|{\displaystyle \frac{n\left(2\lambda +n1\right)(\lambda +1)}{2(\lambda +n1)(\lambda +n)}}\right|^{{}_{}{}^{1}/_{2}^{}},s_n0,`$
$`r`$ $`=`$ $`i\sqrt{2\left|\lambda +1\right|},\theta (x)={\displaystyle \frac{x}{\sqrt{x^21}}},`$
and

$$d_n=\left|\mathrm{\Gamma }(12\lambda n)\mathrm{\Gamma }(n+1)/\left(\lambda +n\right)\right|^{{}_{}{}^{1}/_{2}^{}},\lambda <1,0n\left|\lambda \right|.$$
(21)
For $`p_n^\lambda (x)`$ is correct

$$p_n^\lambda (x)=\mathrm{\hspace{0.17em}2}^\lambda i^n\frac{\mathrm{\Gamma }(2\lambda n+1)n!}{\mathrm{\Gamma }\left(\lambda +1\right)}\left(x^21\right)^{\lambda /2}C_n^\lambda \left(\frac{x}{\sqrt{x^21}}\right),\lambda <1,0n|\lambda |$$
(22)
where $`C_n^\lambda `$ are Gegenbauer functions  of a complex argument $`z`$

$$C_\alpha ^\lambda (z)=\frac{\mathrm{\Gamma }(\alpha +2\lambda )}{\mathrm{\Gamma }(\alpha +1)\mathrm{\Gamma }(2\lambda )}{}_{2}{}^{}F_{1}^{}\left(\genfrac{}{}{0pt}{}{\alpha ,\alpha +2\lambda }{\lambda +1/2}|\frac{1z}{2}\right).$$
(23)
Functions $`p_n^\lambda `$ are orthogonal on the discrete interval:

$$\underset{l=0}{\overset{N1}{}}\sigma ^{(N)}(x_l)p_m^\lambda (x_l)p_n^\lambda (x_l)=\delta _{m,n}d_n^2,N>\left|1\lambda \right|,\lambda <1$$
(24)
with respect to the weight function

$$\sigma ^{(N)}(x_l)=\frac{𝒦_N}{𝒦_{N1}}\left(\frac{2^\lambda \mathrm{\Gamma }(1\lambda )d_{N1}}{\mathrm{\Gamma }(N)\mathrm{\Gamma }(22\lambda N)}\right)^2\left(x_l^21\right)^\lambda \left(C_{N1}^\lambda (z_l)\frac{d}{dz_l}C_N^\lambda (z_l)|_{z_l=x_l/\sqrt{x_l^21}}\right)^1$$
(25)
where $`z_l`$ are the roots of a polynomial of degree $`N`$, $`C_N^\lambda (z_l)=0,x_l=z_l/\sqrt{z_l^21}`$ and $`𝒦_N`$ is the higher coefficient of this polynomial. Thus the solution of (1) for the system of $`N`$ equidistant levels in a resonant field (at $`f_N0`$ and $`\epsilon _n0`$) looks like

$$a_n^{(N)}(t)=\underset{l=0}{\overset{N1}{}}\sigma ^{(N)}(x_l)\left(\frac{p_m^\lambda (x_l)}{d_m}\right)^{}\frac{p_n^\lambda (x_l)}{d_n}\mathrm{exp}\left[\frac{itrx_l}{\sqrt{x_l^21}}\right].$$
(26)
## 7 FROM QUANTUM SYSTEMS TO ORTHOGONAL POLYNOMIALS

Previously we initially selected orthogonal functions to obtain the solution of the equation (1) for qualitatively and quantitatively different quantum systems. Having done so we get the compact solutions which are expressed via known special or even elementary functions. When OPS contains the parameters the solution describes the dynamics of continual or discrete family of quantum systems. From the other side the imperfection of this method is that it is not always possible to obtain exact analytical solutions for a target system with the given dipole moments function $`f_n`$.

Let’s consider the quantum system with an arbitrary dipole moments function being restricted with the case of resonant excitation and a finite number of equidistant levels. Now (1) reduces to the system of $`N+1`$ uniform linear differential equations

$`ida_0(t)/dt`$ $`=`$ $`f_1a_1(t);`$
$`ida_n(t)/dt`$ $`=`$ $`f_{n+1}a_{n+1}(t)+f_na_{n1}(t);n=\overline{1,N1}`$ (27)
$`ida_N(t)/dt`$ $`=`$ $`f_Na_{N1}(t)`$
(some special cases of $`f_n`$ have been studied in ).

Let us seek the solution of (7) in the form

$$a_n(t)=\underset{k=0}{\overset{N}{}}\sigma _kp_n\left(\mathrm{\Lambda }_k\right)\mathrm{exp}\left(i\mathrm{\Lambda }_kt\right)$$
(28)
where $`p_n`$, $`\mathrm{\Lambda }_k`$ and $`\sigma _k`$ are the functions to be determined. Substituting (28) into (7) we come to a three-term recurrence relation for $`p_n`$:

$$f_{n+1}p_{n+1}\left(\mathrm{\Lambda }_k\right)+f_np_{n1}\left(\mathrm{\Lambda }_k\right)=\mathrm{\Lambda }_kp_n\left(\mathrm{\Lambda }_k\right),p_1=0,p_0=\text{const}$$
(29)
Therefore $`p_n`$ are the polynomials orthogonal with respect to variable $`\mathrm{\Lambda }_k`$ (generally on a nonuniform grid).

This truncated orthogonal polynomials $`p_n\left(\mathrm{\Lambda }_k\right)`$ differ from the common polynomials listed above which satisfy the same dipole moment function $`f_n`$ with an interval of $`n`$ change. Let’s discover the basic characteristics of these polynomials.

It is known that the recurrence relation like (3) determines an orthogonal polynomial up to an arbitrary constant factor. The comparison of (3) and (29) shows that the polynomials $`p_n\left(\mathrm{\Lambda }_k\right)`$ are orthonormalized. Therefore it is natural to choose the normalization $`p_0\left(\mathrm{\Lambda }_k\right)1`$.

Suppose that initially all the particles are concentrated on the lower level, in other words $`a_n\left(0\right)=\delta _{n,0}`$. Then it follows from (28) that

$$\underset{k=0}{\overset{N}{}}\sigma _kp_n\left(\mathrm{\Lambda }_k\right)p_0\left(\mathrm{\Lambda }_k\right)=\delta _{n,0}.$$
(30)
Thus the integration constants $`\sigma _k`$ are the weight coefficients of polynomials $`p_n`$ with respect to the variable $`\mathrm{\Lambda }_k`$.

To determine $`\mathrm{\Lambda }_k`$ we shall use the last equation of (7); it gives

$$f_Np_{N1}\left(\mathrm{\Lambda }_k\right)=\mathrm{\Lambda }_kp_N\left(\mathrm{\Lambda }_k\right).$$
(31)
Comparing (31) and (29) it is possible to conclude that $`\mathrm{\Lambda }_k`$ are the roots of the polynomial $`p_{N+1}\left(\mathrm{\Lambda }_k\right)`$.

Calculation of the polynomials $`p_n\left(\mathrm{\Lambda }_k\right)`$ by the recurrence relation (29) and determination of weight coefficients $`\sigma _k`$ from the system (30) are possible for any values of $`N`$. At the same time (30) is reduced eventually to the solving in radicals of the polynomial equations of the degree $`N+1`$. Generally it is possible only for $`N+14`$. However in the case of a resonant excitation the polynomials $`p_n`$ contain only even or odd degrees of an argument (depending on $`n`$ parity). Thus the roots of the polynomials up to the ninth degree can be always found analytically. However using some computer algebra system (e.g. Mathematica) it is possible to obtain the roots for the polynomials of higher degrees for some simple $`f_n`$. In any case it is possible to find these roots numerically with an arbitrary precision.

The method allows to obtain the solutions describing the dynamics of the quantum systems with an arbitrary dipole moment dependence on level number. When the solution of (1) with the same $`f_n`$ but for infinite number of levels can be received under the formula (2) by selecting the relevant orthogonal polynomials it is interesting to compare the truncated polynomials with the common ones. From the comparison of (3) and Eq. (29) it follows that

$$p_n\left(\mathrm{\Lambda }_k\right)=p_n^{(c)}\left(\mathrm{\Lambda }_k/r\right)$$
(32)
where $`p_n^{(c)}`$ are common polynomials, $`r=b_0d_0/d_1`$ , $`d_n`$ is a squared norm, $`b_0`$ is defined from the recurrence relations for unnormalized polynomials: $`p_{n+1}(x)+b_np_{n1}(x)=`$ $`\left(rx+s_n\right)p_n(x)`$. In particular if $`R_k^{\left(N+1\right)}`$ is the $`k^{\text{th}}`$ root of $`p_{N+1}^{(c)}`$ then $`\mathrm{\Lambda }_k=rR_k^{\left(N+1\right)}`$.

Methods (13) and (28) can be generalized to the system excitation in an arbitrary laser field, to the account of non-adjacent dipole transitions and for the systems of degenerate levels.


format: text

text:
# Flip dynamics in octagonal rhombus tiling sets

\[

## Abstract

We investigate the properties of classical single flip dynamics in sets of two-dimensional random rhombus tilings. Single flips are local moves involving 3 tiles which sample the tiling sets via Monte Carlo Markov chains. We determine the ergodic times of these dynamical systems (at infinite temperature): they grow with the system size $`N_T`$ like $`Cst.N_T^2\mathrm{ln}N_T`$; these dynamics are rapidly mixing. We use an inherent symmetry of tiling sets and a powerful tool from probability theory, the coupling technique. We also point out the interesting occurrence of Gumbel distributions.

PACS numbers: 05.40.-a, 02.50.Ga, 61.44.Br

\]

After the discovery of quasicrystals , quasiperiodic tilings  as well as their randomized counterpart, random rhombus tilings , rapidly appeared to be suitable paradigmatic models for quasicrystalline alloys . Simultaneously, these systems also became an active topic in discrete mathematics (see  or  for examples). Fig. 1 displays a random tiling, which belongs to the class of plane tilings with octagonal symmetry. Beyond this case, plane tilings with larger symmetries including Penrose tilings  and space tilings with icosahedral symmetries were proposed to model every kind of quasicrystal. The present letter is devoted to dynamical properties of random rhombus tilings in terms of local dynamical rules, the so-called phason-flips, which consist of local rearrangement of tiles (Fig. 1). These dynamics are of interest for several reasons. On the one hand, it is more and more clear that phason-flips exist in real quasicrystals  and can be modeled in a first approximation by tile-flips; they are a new source of atomic mobility, as compared to usual crystalline materials. In particular, they could carry their own contribution to self-diffusion , even if the efficiency of such processes remains controversial . They are also involved in some specific mechanical properties of quasicrystals, such as plasticity via dislocation mobility . Therefore a complete understanding of flip dynamics is essential in quasicrystal physics. The present work is a first step in this direction. On the other hand, a lot of numerical work has been carried out to characterize statistical properties of tiling sets, a part of which was based on Monte Carlo techniques which rely on a faithful sampling of tiling sets (see ). So far, no systematic study of the relaxation times between two independent numerical measures has been accomplished, whereas it is an essential ingredient for a suitable control of error bars. However, there exist exact results in the simplest case of random rhombus tilings with hexagonal symmetry  and several estimates of relaxation times in larger symmetries, either numerical or in the approximate frame of Langevin dynamics .

Random rhombus tilings are made of rhombi of unitary side length. They are classified according to their global symmetries . The simplest class of hexagonal tilings – made of 60o rhombi with 3 possible orientations – has been widely explored . Tilings with octagonal symmetry are made of 6 different tiles (characterized by their shape and orientation): two squares and four 45o rhombi (Fig. 1). Beyond these two cases, one can define tilings with higher symmetries (e.g. Penrose tilings ) or of higher dimensions . For sake of technical simplicity, we focus on tilings filling a centrally symmetric polygon with integral side lengths (Fig. 1). We are interested in the large size limit where the polygon becomes large keeping a fixed shape. Such tilings can schematically be seen as frozen near their boundary, and strain-free in their central region . Dynamics of fixed boundary tilings should be the manifestation of dynamics of their strain-free center, thus relating both boundary conditions. Here we suppose that all tilings have the same probability; we work at infinite temperature.

The set of all the tilings of such a region together with the flip dynamical rule define a discrete time Markov chain: at each step, a vertex of the tiling is uniformly chosen at random and if this vertex is surrounded by 3 tiles in flippable configuration, then we flip it. Since sets of plane rhombus tilings are connected via elementary flips , this process can reach any tiling. It converges toward the uniform equilibrium distribution, since it satisfies detailed balance. All the difficulty is to characterize how many flips one needs to get close to equilibrium. Generally speaking, let us consider a Markov chain on a finite configuration space $`L`$, which converges toward a stationary distribution $`\pi `$. Let $`x_0`$ be any initial configuration and $`P(x,t|x_0,0)`$ be the probability that the process has reached the configuration $`x`$ after $`t`$ steps. Then

$$\mathrm{\Delta }(t,x_0)=1/2\underset{xL}{}|P(x,t|x_0,0)\pi (x)|$$
(1)
usually measures the distance between both distributions . Given $`\epsilon >0`$ we define the ergodic or mixing time $`\tau (\epsilon )`$ so that whatever $`x_0`$, after $`\tau (\epsilon )`$ steps, one is sure to stay within distance $`\epsilon `$ of equilibrium:

$$\tau (\epsilon )=\underset{x_0}{\mathrm{max}}\underset{t_0}{\mathrm{min}}\{t_0/tt_0,\mathrm{\Delta }(t,x_0)\epsilon \}.$$
(2)
In this letter, we prove that for a tiling of $`N_T`$ tiles

$$\tau (\epsilon )Cst.(N_T)^\nu \mathrm{ln}N_T\mathrm{ln}(1/\epsilon ).$$
(3)
More precisely, we establish, with the help of reduced numerical work, that such a bound holds for $`\nu =3`$, then we argue that $`\nu =2`$ should be the correct exponent. Whatever this exponent, this proves that flip dynamics are rapidly mixing at infinite temperature. As for the $`\mathrm{ln}(N_T)`$ correction in (3), as discussed in , it is a feature of our choice of distance $`\mathrm{\Delta }(t,x_0)`$: an Euclidean norm would not display this correction but it is less natural in the context of measure of probability distributions convergence.

The coupling technique  has been successfully applied to estimate mixing times of several systems, such as hexagonal tilings . It relies on the surprising idea that following the dynamics of couples of configurations instead of a single one might provide the properties of the original dynamics on single configurations. A coupling is a Markov chain on $`L\times L`$; couples of configurations are updated simultaneously and are strongly correlated, but each configuration, viewed in isolation, performs transitions of the original Markov chain. Moreover, the coupled process is designed so that when both configurations happen to be identical, then they follow the same evolution and remain identical forever. Then the central idea of the technique is that the average time the two configurations need to couple (or to coalesce) provides a good upper bound on the original mixing time $`\tau (\epsilon )`$: given an initial couple $`(x_0,y_0)`$ at time $`t=0`$, define the coalescence time $`T(x_0,y_0)`$ as the minimum time that both configurations need to coalesce, and the coupling time as

$$T=\underset{(x_0,y_0)L\times L}{\mathrm{max}}<T(x_0,y_0)>,$$
(4)
where the last mean is taken over realizations of the coupled Markov chain. The following theorem, central in the coupling approach, provides an upper bound for the ergodic times of the original (not coupled) process :

$$\tau (\epsilon )Te\mathrm{ln}(1/\epsilon )+1Te\mathrm{ln}(1/\epsilon ),$$
(5)
where $`e=\mathrm{exp}(1)`$. If the configuration set $`L`$ can be endowed with a partial order relation $``$ with unique minimum and maximum elements $`\widehat{0}`$ and $`\widehat{1}`$, the implementation of the technique is highly facilitated, provided the coupled dynamics is monotonous, i.e. if $`x(t)y(t)`$, then $`x(t+1)y(t+1)`$: let $`(x_0,y_0)`$ be any initial couple such that $`\widehat{1}x_0y_0\widehat{0}`$. After any number of steps, the 4 configurations remain in this order. When the iterates of $`\widehat{0}`$ and $`\widehat{1}`$ have coalesced, the iterates of $`x_0`$ and $`y_0`$ also have. Thus $`T(x_0,y_0)T(\widehat{0},\widehat{1})`$ and $`T=<T(\widehat{0},\widehat{1})>`$.

A convenient representation of random rhombus tilings was introduced by de Bruijn . It consists of following in a tiling lines made of adjacent tiles sharing an edge with a fixed orientation (Fig. 1). The set of lines associated with an orientation is called a de Bruijn family. In an octagonal tiling, there are 4 families. When removing a family from an octagonal tiling, one gets an hexagonal tiling. Conversely, this remark enables one to propose a convenient construction of tilings : directed paths are chosen on an hexagonal tiling, called the base tiling. They are represented by dark lines in Fig. 2. They go from left to right without crossing (but they can have contacts). When they are “opened” following a new edge orientation, they generate de Bruijn lines of the fourth family. In this one-to-one representation, a tiling flip involving tiles of the fourth family becomes a path flip: the path jumps from one side of a tile to the opposite side.

Defining couplings on the whole octagonal tiling sets seems to be an infeasible task. We instead use the paths-on-tiling point of view to decompose the configuration space into smaller subsets where couplings can be defined: let $`J_a`$ denote the set of tilings which have the same base hexagonal tiling $`a`$, called “fibers” . $`L`$ is a disjoint union of fibers. The only possible flips inside $`J_a`$ are those which involve the fourth de Bruijn family. Note that we can construct four such fibrations, one for each family.

Now, let $``$ denote the symmetric transition matrix associated with the Markov chain on the whole set $`L`$: given two configurations $`x`$ and $`y`$, the matrix entry $`(x,y)`$ is equal to the transition probability $`P(x,t+1|y,t)`$. In the same way, we define the symmetric transition matrices $`_i`$ associated with the Markov chains where only flips involving the $`i`$-th de Bruijn family are allowed. Since fibers have been disconnected, $`_i`$ is block-diagonal. The following result interconnects the four fibrations:

$$=(_1+_2+_3+_4)/3\text{Id}/3,$$
(6)
where Id is the identity. Indeed, each coefficient $`(x,y)`$ appears in all four matrices $`_i`$ but one, since the corresponding flip involves 3 de Bruijn lines.

Now we implement the above coupling technique on each fiber. To begin with, we suppose that there is only one line in the flipping family, denoted by $`\mathrm{}`$. As in reference , in order to have a monotonous coupling, we slightly modify the Markov chain: at each step, we choose uniformly at random an internal vertex of $`\mathrm{}`$, the $`n`$-th one starting form the left, and a number $`r\{0,1\}`$. If this vertex is flippable upward (resp. downward) and $`r=0`$ (resp. $`r=1`$), then we flip it. Note that this Markov chain has a time unit different from the original one.

We now define the order relation ($``$): given two lines $`\mathrm{}_1`$ and $`\mathrm{}_2`$, $`\mathrm{}_1\mathrm{}_2`$ if $`\mathrm{}_1`$ is entirely above $`\mathrm{}_2`$. The maximum (resp. minimum) configuration clearly lies on the top (resp. bottom) boundary of the hexagonal domain. If the two flips on $`\mathrm{}_1`$ and $`\mathrm{}_2`$ occur with same $`n`$ and $`r`$ and if $`\mathrm{}_1\mathrm{}_2`$ then their images satisfy the same relation. Indeed, as in reference , thanks to the introduction of $`r`$, if a flip could bring $`\mathrm{}_1`$ below $`\mathrm{}_2`$, then the same flip would also apply to $`\mathrm{}_2`$, thus preserving the order between lines: the coupling is monotonous.

In the general case with $`p`$ non-intersecting lines in each configuration (Fig. 2), let us denote by $`\mathrm{}_i^{(j)}`$, $`j=1,\mathrm{},p`$, the $`p`$ lines of each configuration $`\gamma _i`$. Then $`\gamma _1\gamma _2`$ if for any $`j`$, $`\mathrm{}_1^{(j)}\mathrm{}_2^{(j)}`$. The configuration $`\gamma `$ is maximum (resp. minimum) when each of its lines is maximum (resp. minimum). At each step, the index $`j`$ of the line to be flipped is chosen between 1 and $`p`$, the same $`j`$ for both $`\gamma _i`$.

To begin with, we numerically study the diagonal case, where the 4 sides of the octagonal tilings are equal to $`k`$. For a given base tiling $`𝒯_a`$, we run a number $`m`$ of couplings until they coalesce, and then estimate the coupling time $`T(𝒯_a)`$. We then make a second average on $`M`$ different tilings, in order to get the time $`\overline{T}`$ averaged on tilings $`𝒯_a`$. We also keep track of the standard deviation $`\mathrm{\Delta }T`$. From our numerical data (see Fig. 3), we draw the following conclusions: $`\mathrm{\Delta }T/\overline{T}`$ decreases toward a constant ($`0.07`$) as $`k\mathrm{}`$, which means that the average coupling time $`T(𝒯_a)`$ goes on depending on the base tiling $`𝒯_a`$ at the large size limit. However, most $`T(𝒯_a)`$ are of order $`\overline{T}`$, and the mixing times $`\tau (\epsilon )`$ on most fibers are controlled by $`\overline{T}`$. Nevertheless, the effect of few “slower” fibers will deserve a detailed discussion below. Moreover, the measures of $`\overline{T}`$ are compatible with a $`k^4\mathrm{ln}k`$ behavior (Fig. 3, inset). In particular, this fit with logarithmic corrections is much better than a simple power-law fit. This result is consistent with known results in the case of hexagonal tilings , where $`T`$ also grows like $`k^4\mathrm{ln}k`$.

We also have explored coupling times on fibers in non-diagonal cases and our conclusions remain identical. As a consequence, couplings in fibers behave like couplings in hexagonal tiling problems, up to different numerical prefactors: the dynamics on each fiber is rapidly mixing.

Now we return to the dynamics on the whole set of tilings. Examining Fig. 4, one remarks that two fibrations are to a certain extent “transverse”: in octagonal tiling sets, one can connect any two tilings using flips of only two fibers ; if the dynamics is rapidly mixing in each fiber, the combination of dynamics on two (and even four) fibrations will certainly also be rapidly mixing. We now establish properly this point: it is common in the field of Markov processes to relate rates of convergence to spectra of transition matrices. Generally speaking, given a transition matrix $`M`$, 1 is always the largest eigenvalue in modulus, and the difference $`g(M)`$ between 1 and the second largest eigenvalue is called the first gap of $`M`$. Then $`\tau (\epsilon )\mathrm{ln}(1/\epsilon )/g()`$ for small $`\epsilon `$ . By simple arguments from linear algebra and Euclidean geometry (since matrices are symmetric), the following central gap relation can be established , based upon (6) and the smallness of the intersection of two fibers:

$$g()\underset{i}{inf}\left(g(_i)\right).$$
(7)
This result restores the symmetry lost in the fibration process. It implies that the mixing time $`\tau (\epsilon )`$ on $`L`$ is smaller than the mixing time on the slowest fibration derived from (5). Hence $`\tau (\epsilon )`$ is smaller than the mixing time on the slowest fiber. We have seen that the average coupling time $`\overline{T}`$ grows like $`k^4\mathrm{ln}k`$ in diagonal cases. Since the number of tiles is $`N_T=6k^2`$, $`\overline{T}\kappa _4N_T^2\mathrm{ln}(N_T)`$, where $`\kappa _4=1.189\pm 0.003`$ in the original time unit. However, the coupling time $`T`$ depends on the base tiling $`𝒯_a`$ and the distribution of times $`T`$ in a given fibration have a certain width around $`\overline{T}`$ (Fig. 3). And even if the typical values of coupling times are of the order of magnitude of the previous average value, this does not exclude the existence of rare slow fibers in the upper tails of these distributions. However, we recall that coupling times $`T`$ are maxima of coalescence time distributions (4). Therefore the expected shape of the distribution of $`T(𝒯_a)`$ ought to be sought in the specific class of extreme-value distributions, namely Gumbel distributions : consider $`N`$ independent identical random variables $`T_a`$, whose probability densities decay rapidly at large $`T`$:

$$p(T)\frac{C_1}{T^\alpha }\mathrm{exp}(C_2T^\beta ),$$
(8)
where $`C_1,C_2,\beta >0`$. If $`T_{\text{max}}=\mathrm{max}_aT_a`$, then at large $`N`$, the probability density of $`T_{\text{max}}`$ satisfies

$$p(u)=\mathrm{exp}(u\mathrm{exp}(u)),$$
(9)
where $`u=(T_{\text{max}}T_0)/\delta T`$ is a suitably rescaled variable.

Now, even if coalescence times are not strictly speaking independent variables , our numerical distributions appear to be well fitted by this kind of distribution (Fig. 3). This result provides the large $`T`$ behavior of coupling time distributions: $`p(T_{\text{max}})\mathrm{exp}(u)\mathrm{exp}(T_{\text{max}}/\delta T)`$, and therefore an estimation of the largest coupling time. Let us focus on diagonal cases: $`T_0`$ as well as $`\delta T`$ behave like $`k^4\mathrm{ln}k`$. But for a fibration $`i`$, there are $`N_i`$ base tilings $`𝒯_a`$. Therefore if $`T^{}`$ is the largest coupling time on all tilings $`𝒯_a`$, it is estimated by $`N_i\mathrm{exp}(T^{}/\delta T)1`$. Now $`N_i`$ grows exponentially with the number of tiles : $`\mathrm{ln}N_iCst_1k^2`$. Thus

$$T^{}Cst_2k^6\mathrm{ln}kCst_3N_T^3\mathrm{ln}N_T$$
(10)
and $`\tau (\epsilon )Cst_4N_T^3\mathrm{ln}N_T\mathrm{ln}(1/\epsilon )`$. However these extreme values should not be relevant: because of the exponential decay of $`p(u)`$, slow fibers are rare and can be bypassed via rapid ones. More precisely, perturbation theory arguments suggest that rare slow fibers can be seen as a small perturbation of an otherwise rapid transition matrix and have a vanishing influence on its spectral gap : they do not significantly slow rapid dynamics and the typical coupling time $`\overline{T}`$ drives the dynamics on fibers, leading to $`\nu =2`$ in (3) and $`Cst.e\kappa _4`$.

What does this analysis become in the case of larger symmetry tilings or of higher dimensional tilings, such as icosahedral ones? Plane rhombus tilings with $`2D`$-fold symmetry could be addressed by our approach without significant technical complication, leading to laws similar to (3), up to different prefactors $`Cst.`$ . As for higher dimensional tilings, the fibration process remains valid, but the connectivity of fibers is not established , making impossible a naïve generalization of this approach.

I would like to thank M. Latapy, K. Frahm, R. Mosseri and D.S. Dean for fruitful discussions.