text:
# A strange recursion operator for a new integrable system of coupled Korteweg–de Vries equations

## 1 Introduction

Recently, the Painlevé classification of coupled Korteweg–de Vries (KdV) equations was made in , and the following new system possessing the Painlevé property was found there:

$`u_t`$ $`={\displaystyle \frac{\delta c_1}{12}}u_{xxx}+{\displaystyle \frac{3\delta c_1^2}{4d_1}}v_{xxx}+c_1u_xu{\displaystyle \frac{3c_1^2}{d_1}}u_xv{\displaystyle \frac{6c_1^2}{d_1}}v_xu,`$
$`v_t`$ $`={\displaystyle \frac{\delta d_1}{12}}u_{xxx}{\displaystyle \frac{7\delta c_1}{12}}v_{xxx}+d_1u_xuc_1u_xv2c_1v_xu{\displaystyle \frac{6c_1^2}{d_1}}v_xv,`$ (1)
where $`\delta `$, $`c_1`$ and $`d_1`$ are arbitrary nonzero constants.

More recently, the coupled KdV equations (1) appeared in , also as a system possessing the Painlevé property, in the following form:

$`u_t`$ $`=u_{xxx}9v_{xxx}12uu_x+18vu_x+18uv_x,`$
$`v_t`$ $`=u_{xxx}5v_{xxx}12uu_x+12vu_x+6uv_x+18vv_x,`$ (2)
which is related to (1) by the transformation (1)$``$(2):

$$u\delta u\frac{3\delta }{2}v,v\frac{\delta d_1}{2c_1}v,t\frac{6}{\delta c_1}t.$$
(3)
It was found in  that the system (2) admits the zero-curvature representation (ZCR)

$$D_tX=D_xT[X,T],$$
(4)
which is the compatibility condition for the overdetermined linear system

$$\mathrm{\Psi }_x=X\mathrm{\Psi },\mathrm{\Psi }_t=T\mathrm{\Psi },$$
(5)
with

$$X=\left(\begin{array}{ccc}\sigma & uv& \sigma \left(\frac{8}{3}u4v\right)\\ \frac{3}{2}& 2\sigma & uv\\ 0& \frac{3}{2}& \sigma \end{array}\right),$$
(6)
where $`\sigma `$ is an essential (‘spectral’) parameter, $`D_t`$ and $`D_x`$ stand for the total derivatives, the square brackets denote the commutator, $`\mathrm{\Psi }(x,t)`$ is a three-component column, and an explicit (but cumbersome) expression for the $`3\times 3`$ matrix $`T`$ can be seen in .

In the present paper, we study this new integrable system of coupled KdV equations in its equivalent form, which is characterized by more transparent structure of $`X`$ in its ZCR. Namely, making the transformation

$`u`$ $`2u{\displaystyle \frac{1}{3}}v,v{\displaystyle \frac{4}{3}}u{\displaystyle \frac{1}{3}}v,t{\displaystyle \frac{1}{2}}t,`$
$`\mathrm{\Psi }`$ $`S\mathrm{\Psi },XSXS^1,TSTS^1,S=\mathrm{diag}({\displaystyle \frac{2}{3}},1,{\displaystyle \frac{3}{2}})`$ (7)
in (2), (5) and (6), we obtain the system

$`u_t`$ $`=4u_{xxx}v_{xxx}12uu_x+vu_x+2uv_x,`$
$`v_t`$ $`=9u_{xxx}2v_{xxx}12vu_x6uv_x+4vv_x`$ (8)
which admits the ZCR (4) with

$$X=\left(\begin{array}{ccc}\sigma & u& \sigma v\\ 1& 2\sigma & u\\ 0& 1& \sigma \end{array}\right)$$
(9)
(the explicit expression for $`T`$ is omitted, since not used in what follows).

Our aim is to obtain a recursion operator for the system (8). We derive the recursion operator from the matrix $`X`$ (9) of the system’s ZCR. We do this, following the method of gauge-invariant description of ZCRs of evolution equations, introduced in  (note also a more general gauge-invariant formulation of ZCRs, developed in ).

The paper is organized as follows. In Section 2, we find the complete class of coupled evolution equations admitting ZCRs (4) with the matrix $`X`$ given by (9). This class is determined implicitly, through a three-term recurrence, which is explicitly solved in Section 3, giving rise to the recursion operator sought. In Section 4, we prove that this operator satisfies the sufficient condition of being a recursion operator in the sense of symmetries. Section 5 contains discussion of the results.

## 2 Class of systems

Let us find the complete class of coupled evolution equations

$`u_t`$ $`=f(x,t,u,\mathrm{},u_{x\mathrm{}x},v,\mathrm{},v_{x\mathrm{}x}),`$
$`v_t`$ $`=g(x,t,u,\mathrm{},u_{x\mathrm{}x},v,\mathrm{},v_{x\mathrm{}x})`$ (10)
which admit ZCRs (4) with the matrix $`X`$ given by (9) and without any restrictions imposed on the $`3\times 3`$ matrix $`T(x,t,u,\mathrm{},u_{x\mathrm{}x},v,\mathrm{},v_{x\mathrm{}x},\sigma )`$. We solve this problem algorithmically, using the technique of cyclic bases of ZCRs, developed in  (see also  for more examples).

Since the matrix $`X`$ (9) does not contain derivatives of $`u`$ and $`v`$, the characteristic matrices of the ZCR are simply $`C_u=X/u`$ and $`C_v=X/v`$:

$$C_u=\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right),C_v=\left(\begin{array}{ccc}0& 0& \sigma \\ 0& 0& 0\\ 0& 0& 0\end{array}\right).$$
(11)
Using the operator $`_x`$, defined as $`_xY=D_xY[X,Y]`$ for any (here, $`3\times 3`$) matrix $`Y`$, we compute $`_xC_u`$, $`_xC_v`$, $`_x^2C_u`$, $`_x^2C_v`$, etc., and find that the cyclic basis (i.e. the maximal sequence of linearly independent matrices $`_x^iC_u`$ and $`_x^jC_v`$, $`i,j=0,1,2,\mathrm{}`$) is

$$\{C_u,C_v,_xC_u,_xC_v,_x^2C_u,_x^2C_v,_x^3C_u,_x^4C_u\},$$
(12)
and that the closure equations for this cyclic basis are

$`_x^5C_u`$ $`=a_1C_u+a_2C_v+a_3_xC_u+a_4_xC_v+a_5_x^2C_u+a_6_x^2C_v`$
$`+a_7_x^3C_u+a_8_x^4C_u,`$
$`_x^3C_v`$ $`=b_1C_u+b_2C_v+b_3_xC_u+b_4_xC_v+b_5_x^2C_u+b_6_x^2C_v`$
$`+b_7_x^3C_u+b_8_x^4C_u,`$ (13)
where

$`a_1`$ $`=u_{xxx}11uu_x+{\displaystyle \frac{3}{2}}vu_x27\sigma ^2u_x+6\sigma ^2v_x,`$
$`a_2`$ $`=6u_{xxx}+2v_{xxx}+54uu_x3vu_x22uv_x+3vv_x`$
$`+54\sigma ^2u_x18\sigma ^2v_x,`$
$`a_3`$ $`=4u_{xx}22u^2+3uv99\sigma ^2u+{\displaystyle \frac{45}{2}}\sigma ^2v81\sigma ^4,`$
$`a_4`$ $`=24u_{xx}+7v_{xx}+54u^236uv+{\displaystyle \frac{9}{2}}v^2+54\sigma ^2u9\sigma ^2v,`$
$`a_5`$ $`=5u_x,a_6=36u_x+8v_x,`$
$`a_7`$ $`=13u{\displaystyle \frac{3}{2}}v+18\sigma ^2,a_8=0,`$
$`b_1`$ $`={\displaystyle \frac{1}{2}}u_x,b_2=u_x+v_x,b_3=u+{\displaystyle \frac{3}{2}}\sigma ^2,`$
$`b_4`$ $`=u+{\displaystyle \frac{3}{2}}v,b_7={\displaystyle \frac{1}{2}},b_5=b_6=b_8=0.`$ (14)
Thus, in the case of the matrix $`X`$ (9), the cyclic basis is eight-dimensional. Therefore the singular basis (see  for the definition) can be one-dimensional at most; in fact, it consists of the unit matrix $`I=\mathrm{diag}(1,1,1)`$ and, due to $`_xI=0`$, has no effect on the form of the represented equations (10).

Now, we rewrite the ZCR (4) in its equivalent (characteristic) form

$$fC_u+gC_v=_xT$$
(15)
and decompose the matrix $`T`$ over the cyclic basis as

$`T`$ $`=p_1C_u+p_2C_v+p_3_xC_u+p_4_xC_v+p_5_x^2C_u+p_6_x^2C_v`$
$`+p_7_x^3C_u+p_8_x^4C_u,`$ (16)
where the coefficients $`p_1,\mathrm{},p_8`$ are functions of $`x,t,u,\mathrm{},u_{x\mathrm{}x},v,\mathrm{},v_{x\mathrm{}x}`$ and $`\sigma `$. Taking into account the closure equations (13) and the linear independence of the matrices (12), we obtain from (15) and (16) the following relations:

$`f`$ $`=D_xp_1+a_1p_8+b_1p_6,g=D_xp_2+a_2p_8+b_2p_6,`$
$`p_1`$ $`=D_xp_3a_3p_8b_3p_6,p_2=D_xp_4a_4p_8b_4p_6,`$
$`p_3`$ $`=D_xp_5a_5p_8b_5p_6,p_4=D_xp_6a_6p_8b_6p_6,`$
$`p_5`$ $`=D_xp_7a_7p_8b_7p_6,p_7=D_xp_8a_8p_8b_8p_6.`$ (17)
These relations (17), together with (14), determine $`p_1,p_2,p_3,p_4,p_5,p_7,f`$ and $`g`$ in terms of some linear differential operators applied to the two (still undetermined) functions $`p_6`$ and $`p_8`$. For what follows, we need no expression for $`T`$: it is sufficient to know that $`T`$ exists. The expressions for the right-hand sides $`f`$ and $`g`$ of the represented equations (10) are

$$h=\left(M+\lambda L+\lambda ^2K\right)r,$$
(18)
where

$$h=\left(\begin{array}{c}f\\ g\end{array}\right),\lambda =9\sigma ^2,r=\left(\begin{array}{c}p_8\\ p_6\end{array}\right),$$
(19)
$$K=\left(\begin{array}{cc}D_x& 0\\ 0& 0\end{array}\right),$$
(20)
$$L=\left(\begin{array}{cc}2D_x^3+\left(11u\frac{5}{2}v\right)D_x+\left(8u_x\frac{11}{6}v_x\right)& \frac{1}{6}D_x\\ \left(6u+v\right)D_xv_x& 0\end{array}\right),$$
(21)
$$M=\left(\begin{array}{cc}M_{11}& M_{12}\\ M_{21}& M_{22}\end{array}\right)$$
(22)
with

$`M_{11}`$ $`=D_x^5+\left(13u+{\displaystyle \frac{3}{2}}v\right)D_x^3+\left(34u_x+{\displaystyle \frac{9}{2}}v_x\right)D_x^2`$
$`+\left(33u_{xx}+{\displaystyle \frac{9}{2}}v_{xx}+22u^23uv\right)D_x`$
$`+\left(11u_{xxx}+{\displaystyle \frac{3}{2}}v_{xxx}+33uu_x{\displaystyle \frac{3}{2}}vu_x3uv_x\right),`$ (23)
$$M_{12}=\frac{1}{2}D_x^3uD_x\frac{1}{2}u_x,$$
(24)
$`M_{21}`$ $`=\left(36u_x+8v_x\right)D_x^2`$
$`+\left(48u_{xx}+9v_{xx}54u^2+36uv{\displaystyle \frac{9}{2}}v^2\right)D_x`$
$`+\left(18u_{xxx}+3v_{xxx}54uu_x+33vu_x+14uv_x6vv_x\right),`$ (25)
$$M_{22}=D_x^3+\left(u\frac{3}{2}v\right)D_x\frac{1}{2}v_x.$$
(26)
If $`\sigma `$ in (9) is a fixed constant, then (18) and (19)–(26) give the solution of our problem: a continual class of coupled evolution equations, containing two arbitrary functions $`p_6`$ and $`p_8`$, is represented by (4) with (9). However, $`\sigma `$ is a free parameter in (9), and we need to take into account that the condition $`h/\lambda =0`$ must be satisfied for (18). Using the series expansion of $`r`$

$$r=r_0+\lambda r_1+\lambda ^2r_2+\lambda ^3r_3+\lambda ^4r_4+\mathrm{}$$
(27)
(if a singularity is suspected in $`r`$ at $`\lambda =0`$, one may use (27) after an infinitesimal shift of $`\lambda `$, which does not affect final results), we obtain from (18) that

$$h=Mr_0$$
(28)
and

$`0`$ $`=Mr_1+Lr_0,`$ (29a)
$`0`$ $`=Mr_2+Lr_1+Kr_0,`$ (29b)
$`0`$ $`=Mr_3+Lr_2+Kr_1,`$ (29c)
$`0`$ $`=Mr_4+Lr_3+Kr_2,`$ (29d)
$`\mathrm{}.`$
This gives us the following implicit solution of our problem: the system of coupled evolution equations (10) admits the ZCR (4) with the matrix $`X`$ given by (9) if and only if its right-hand side $`h=(f,g)^\mathrm{T}`$ is determined through (28) by the set of two-component functions $`r_0,r_1,r_2,\mathrm{}`$ satisfying the three-term recurrence (29), where the matrix differential operators $`M`$, $`L`$ and $`K`$ are defined by (20)–(26).

## 3 Recursion of systems

Let us solve the recurrence (29), using formal inversion of differential operators. We say that the operator $`P^1`$ is formally inverse to a linear differential operator $`P`$ if $`P^1P=PP^1=I`$, where $`I`$ is the unit operator. In other words, $`P^1a`$ denotes any function $`b`$ of local variables, such that $`a=Pb`$, if $`b`$ exists. We have to remind that, in the present case, differential operators contain total (not partial) derivatives $`D_x`$ and act on differential functions (in Olver’s sense , i.e. functions of independent variables, dependent variables and finite-order derivatives of dependent variables). For this reason, $`D_x^10=\varphi \left(t\right)`$ with any function $`\varphi `$, whereas $`\left(D_x^2+u\right)^10=0`$ and $`\left(3D_x^34vD_x2v_x\right)^1v_x=\frac{1}{2}`$, to give some examples.

For the operators $`K`$, $`L`$ and $`M`$, defined by (20)–(26), we notice that the inverse operators $`M^1`$ and $`L^1`$ exist, but $`K^1`$ does not exist, and therefore we cannot use inversion of $`K`$ in what follows. We obtain

$$r_0=L^1Mr_1$$
(30)
from (29a), and

$$r_1=\left(KL^1ML\right)^1Mr_2$$
(31)
from (29b), and notice that the form of (31) is different from the form of (30), because $`KL^1M0`$ for the operators defined by (20)–(26). Then we obtain

$$r_2=\left(K\left(KL^1ML\right)^1ML\right)^1Mr_3$$
(32)
from (29c), and find that, fortunately, it is possible to rewrite (32) as

$$r_2=\left(KL^1ML\right)^1Mr_3$$
(33)
owing to the condition

$$KL^1K=0$$
(34)
satisfied by $`K`$ (20) and $`L`$ (21). We see that the form of (33) coincides with the form of (31). The same is true, also owing to (34), for the relation between $`r_3`$ and $`r_4`$, which follows from (29d); and, in general,

$$r_i=\left(KL^1ML\right)^1Mr_{i+1},i=1,2,3,\mathrm{}.$$
(35)
We have to emphasize that the recursion (35) is caused by the specific property (34) of the operators $`K`$ (20) and $`L`$ (21), and that it would be impossible to solve the three-term recurrence (29) for generic operators.

For convenience in what follows, we rewrite (35) as

$$r_i=M^1LM^1RML^1Mr_{i+1},i=1,2,3,\mathrm{},$$
(36)
where

$$R=MN^1$$
(37)
with

$$N=ML^1KL.$$
(38)
The explicit expressions for $`N`$ and $`R`$ are

$$N=\left(\begin{array}{cc}D_x^3+\left(5u+\frac{5}{2}v\right)D_x+\left(5u_x+\frac{11}{6}v_x\right)& \frac{1}{6}D_x\\ 6D_x^3+8vD_x+4v_x& 0\end{array}\right)$$
(39)
and

$$R=\left(\begin{array}{cc}R_{11}& R_{12}\\ R_{21}& R_{22}\end{array}\right)$$
(40)
with

$$R_{11}=3D_x^26u3u_xD_x^1,$$
(41)
$`R_{12}`$ $`=[2D_x^5+(2u+3v)D_x^3+(4u_x+8v_x)D_x^2`$
$`+\left(6u_{xx}+7v_{xx}+4u^26uv\right)D_x`$
$`+\left(2u_{xxx}+2v_{xxx}+6uu_x3vu_x4uv_x\right)`$
$`+u_xD_x^1v_x\left]\right(3D_x^34vD_x2v_x)^1,`$ (42)
$$R_{21}=6D_x^2+\left(6u9v\right)3v_xD_x^1,$$
(43)
$`R_{22}`$ $`=[3D_x^5+(18u+12v)D_x^3+(27u_x+18v_x)D_x^2`$
$`+\left(21u_{xx}+14v_{xx}+12u^2+12uv9v^2\right)D_x`$
$`+\left(6u_{xxx}+4v_{xxx}+12uu_x+6vu_x+6uv_x9vv_x\right)`$
$`+v_xD_x^1v_x\left]\right(3D_x^34vD_x2v_x)^1.`$ (44)
Now, using (28), (30) and (36), we obtain the following expression for $`h`$:

$$h=R^nML^1Mr_{n+1},n=0,1,2,\mathrm{}.$$
(45)
Since the leading-order terms of the operator $`R`$ (40)–(44) are given by

$$R=\left(\begin{array}{cc}3D_x^2+\mathrm{}& \frac{2}{3}D_x^2+\mathrm{}\\ 6D_x^2+\mathrm{}& D_x^2+\mathrm{}\end{array}\right)$$
(46)
and the orders of derivatives of $`u`$ and $`v`$ in $`h`$ must be finite, we can choose a sufficiently large $`n`$ in (45), so that the orders of derivatives of $`u`$ and $`v`$ in $`ML^1Mr_{n+1}`$ are less than 2 for this $`n`$. In this case, we conclude, taking into account the form of $`M`$ (22)–(26) and $`L`$ (21), that $`L^1Mr_{n+1}`$ must be a function of $`x`$ and $`t`$ only, that the orders of derivatives of $`u`$ and $`v`$ in $`Mr_{n+1}`$ must be less than 2, and that $`r_{n+1}=(0,0)^\mathrm{T}`$. Then we find that

$$L^1\left(\begin{array}{c}0\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 2\varphi \left(t\right)\end{array}\right),$$
(47)
where the function $`\varphi \left(t\right)`$ is arbitrary, and the factor 2 is taken for convenience. Finally, (45) gives us

$$h=\varphi \left(t\right)R^n\left(\begin{array}{c}u_x\\ v_x\end{array}\right),\varphi \left(t\right),n=0,1,2,\mathrm{},$$
(48)
or, equivalently,

$$h=R^n\left(\begin{array}{c}0\\ 0\end{array}\right),n=1,2,3,\mathrm{},$$
(49)
because $`R(0,0)^\mathrm{T}=(\varphi \left(t\right)u_x,\varphi \left(t\right)v_x)^\mathrm{T}`$ with any $`\varphi \left(t\right)`$.

Thus, the problem, formulated in the beginning of Section 2, is solved explicitly: the relation (49) determines the right-hand sides $`h=(f,g)^\mathrm{T}`$ of all systems (10) which admit the linear problem (5) with the matrix $`X`$ (9). The constructed operator $`R`$ (40)–(44) turns out to be a recursion operator in the Lax sense: it generates a hierarchy of integrable systems which all possess Lax pairs with the predetermined spatial part. The system (8) corresponds to $`n=2`$ in (49), and the next member of the hierarchy ($`n=3`$) is the fifth-order system

$`u_t`$ $`=6u_{xxxxx}{\displaystyle \frac{5}{3}}v_{xxxxx}60uu_{xxx}+{\displaystyle \frac{40}{3}}vu_{xxx}+{\displaystyle \frac{50}{3}}uv_{xxx}`$
$`{\displaystyle \frac{10}{3}}vv_{xxx}150u_xu_{xx}+40v_xu_{xx}+{\displaystyle \frac{125}{3}}u_xv_{xx}10v_xv_{xx}`$
$`+120u^2u_x40uvu_x+{\displaystyle \frac{5}{3}}v^2u_x{\displaystyle \frac{80}{3}}u^2v_x+{\displaystyle \frac{20}{3}}uvv_x,`$
$`v_t`$ $`=15u_{xxxxx}4v_{xxxxx}120uu_{xxx}+10vu_{xxx}+30uv_{xxx}`$
$`{\displaystyle \frac{5}{3}}vv_{xxx}360u_xu_{xx}+80v_xu_{xx}+90u_xv_{xx}{\displaystyle \frac{55}{3}}v_xv_{xx}`$
$`+160uvu_x40v^2u_x+40u^2v_x60uvv_x+{\displaystyle \frac{35}{3}}v^2v_x.`$ (50)
## 4 Recursion of symmetries

The operator $`R`$, defined by (40)–(44), is also a recursion operator with respect to symmetries of the coupled KdV equations (8). In order to see this, we have to prove that $`R`$ satisfies the following condition :

$$D_tR=[H,R],$$
(51)
where $`H`$ is the Fréchet derivative of the right-hand side of the system (8),

$$H=\left(\begin{array}{cc}H_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right)$$
(52)
with

$`H_{11}`$ $`=4D_x^3+\left(12u+v\right)D_x+\left(12u_x+2v_x\right),`$
$`H_{12}`$ $`=D_x^3+2uD_x+u_x,`$
$`H_{21}`$ $`=9D_x^312vD_x6v_x,`$
$`H_{22}`$ $`=2D_x^3+\left(6u+4v\right)D_x+\left(12u_x+4v_x\right).`$ (53)
It is convenient to study the condition (51) in its equivalent form:

$$HMD_tM=R\left(HND_tN\right),$$
(54)
which follows from (51) through the definition of $`R`$ (37) and the identity $`D_tN^1=N^1\left(D_tN\right)N^1`$. Direct computations show that the condition (54) is satisfied identically by $`H`$ (52)–(53), $`M`$ (22)–(26), $`R`$ (40)–(44), $`N`$ (39) and the expressions (8) for $`u_t`$ and $`v_t`$. Since the condition (51) is satisfied, the system of coupled KdV equations (8) possesses infinitely many symmetries: the recursion operator $`R`$ generates a ‘new’ symmetry (but not necessarily a local one) from any ‘old’ symmetry.

## 5 Conclusion

The recursion operator $`R`$ (40)–(44), constructed in this paper for the new integrable system of coupled KdV equations (8), is characterized by unusual structure of its nonlocal part, which contains not only the conventional operator $`D_x^1`$ but also the strange nonlocal operator $`\left(3D_x^34vD_x2v_x\right)^1`$. Such a phenomenon has not been observed in the literature as yet. This operator $`\left(3D_x^34vD_x2v_x\right)^1`$ cannot be expressed in terms of a finite number of operators $`D_x^1`$, if one uses only local variables, as we did throughout the paper. However, if we introduce the nonlocal variable $`w:v=3w_{xx}/w`$, we can represent the obtained recursion operator in a more conventional form, using the factorization $`\left(3D_x^34vD_x2v_x\right)^1=\frac{1}{3}w^2D_x^1w^2D_x^1w^2D_x^1w^2`$.

We found the recursion operator of the system (8) by the technique of cyclic bases of ZCRs. It is interesting if this recursion operator can be obtained by different methods as well, e.g. by those developed in .

We have seen that the relation (49) generates local expressions for $`h`$ at $`n=1,2,3`$. It is easy to check that the seventh-order system of this hierarchy ($`n=4`$) is also a local expression. Nevertheless, we have not proven as yet that this is the case for any $`n`$.

Though our method produces the recursion operator in the quotient form $`MN^1`$ (37), the operators $`M`$ (22)–(26) and $`N`$ (39) are not Hamiltonian operators. It is still unknown if the obtained recursion operator is hereditary and if the system (8) admits a bi-Hamiltonian structure.

## Acknowledgments

The authors are grateful to the Scientific and Technical Research Council of Turkey (TÜBİTAK) for support. One of the authors, S. Yu. S., is also grateful to the Middle East Technical University (ODTÜ) for hospitality.


format: text

text:
# Dynamical Upper Bounds for One-Dimensional Quasicrystals

## 1. Introduction

In this paper we study quantum dynamical properties of one-dimensional quasicrystals. That is, we will be concerned with the family of operators $`H_{\lambda ,\alpha ,\theta }`$ on $`\mathrm{}^2()`$, acting on $`u\mathrm{}^2()`$ by

(1)
$$H_{\lambda ,\alpha ,\theta }u(n)=u(n+1)+u(n1)+\lambda v_{\alpha ,\theta }(n)u(n).$$

Here, $`\lambda >0`$, $`\alpha (0,1)`$ is irrational, $`\theta `$ belongs to $`[0,1)`$, and $`v_{\alpha ,\theta }(n)`$ is given by

(2)
$$v_{\alpha ,\theta }(n)=\chi _{[1\alpha ,1)}(n\alpha +\theta mod1).$$

Operators of the form (1) have been extensively studied since the early eighties. We refer the reader to  for history and known results (“last-millennium results”). The spectral type has been found to always be purely singular continuous . Moreover, in some cases the spectral measures were even found to be absolutely continuous with respect to certain Hausdorff measures  from which one may deduce quantitative quantum dynamical lower bounds (using a general theory initiated by Guarneri , Combes , and Last ).

We will be concerned with quantum dynamical upper bounds in the spirit of Killip, Kiselev, and Last . While there is no general method known to this date which bounds the “fast” part of the wavepacket time evolution from above, the authors of  propose a method to at least bound the spreading of the “slow” part from above. They apply their method to the operator of the form (1), where $`\lambda >8`$, $`\alpha =(\sqrt{5}1)/2`$, and $`\theta =0`$. In this paper we consider the case of general $`\alpha `$ and $`\theta `$.

Before stating our main theorem, we fix some notation. Let $`\delta _n`$ be the element of $`\mathrm{}^2()`$ which is supported at $`n`$ and obeys $`\delta _n(n)=1`$. Given a function $`\psi :`$ and $`L>0`$, we define

(3)
$$\psi _L^2=\underset{n=L}{\overset{L}{}}|\psi (n)|^2+(LL)\left(|\psi (L1)|^2+|\psi (L+1)|^2\right).$$

For a function $`A:[0,\mathrm{})`$, we define

$$A(t)_T=\frac{2}{T}_0^{\mathrm{}}e^{2t/T}A(t)𝑑t.$$

Finally, for $`\alpha (0,1)`$ irrational, we consider its continued fraction expansion

(4)
$$\alpha =\frac{1}{a_1+{\displaystyle \frac{1}{a_2+{\displaystyle \frac{1}{a_3+\mathrm{}}}}}}$$

with uniquely determined $`a_k`$. Let us define the associated rational approximants $`p_k/q_k`$ of $`\alpha `$ by

(5) $`p_0`$ $`=0,`$ $`p_1`$ $`=1,`$ $`p_k`$ $`=a_kp_{k1}+p_{k2},`$
(6) $`q_0`$ $`=1,`$ $`q_1`$ $`=a_1,`$ $`q_k`$ $`=a_kq_{k1}+q_{k2}.`$
The rational numbers $`p_k/q_k`$ are known to be best approximants to $`\alpha `$ and, in particular, $`p_k`$ and $`q_k`$ are relatively prime. See Khinchin  for background on continued fraction expansions.

We let

(7)
$$=\{\alpha (0,1)\text{ irrational}:B\text{ such that }q_k+1B^kk\}.$$

It is known that $``$ has full Lebesgue measure .

Our main result reads as follows:

###### Theorem 1.

There are constants $`C_1,C_2>0`$ such that for every $`\lambda >20`$, every $`\alpha `$, and every $`\theta [0,1)`$, we have

(8)
$$e^{itH_{\lambda ,\alpha ,\theta }}\delta _1_{C_1T^{p(\lambda ,\alpha )}}^2_TC_2$$

for $`T`$ large enough (i.e., for $`TT_0(\alpha ,\theta )`$), where

$$p(\lambda ,\alpha )=\frac{6\mathrm{log}B}{\mathrm{log}\left(\frac{\lambda 8}{3}\right)}$$

and $`B`$ is associated to $`\alpha `$ via (7).

Remarks. (a) The physical interpretation of (8) is the following. At any time $`TT_0(\theta )`$, the averaged probability to find the particle in a ball of radius $`C_1T^{p(\lambda ,\alpha )}`$ is uniformly bounded away from zero. This gives an upper bound on the slow part of the wavepacket time evolution.
(b) Since for fixed $`\alpha `$, $`p(\lambda ,\alpha )0`$ as $`\lambda \mathrm{}`$, we see that the slow part of the wavepacket moves arbitrarily slowly if we choose large enough coupling.
(c) This theorem was shown by Killip et al. in  for the particular case $`\alpha =(\sqrt{5}1)/2`$ and $`\theta =0`$. They only require $`\lambda >8`$. For this particular rotation number $`\alpha `$, our extension to arbitrary phase $`\theta `$ works for the same range of $`\lambda `$-values.
(d) Our extension to arbitrary phase $`\theta `$ is based on a fine analysis of the local combinatorial structure of the sequences $`v_{\alpha ,\theta }`$ which might be of independent interest.
(e) The key ingredient in our extension to almost every rotation number $`\alpha `$ is a recent result of Liu and Wen  which generalizes a paper of Raymond  (which in turn is at the heart of the proof in ) to arbitrary rotation number.
(f) As in , the theorem, while stated for initial vector $`\delta _1`$, can be recast for other initial vectors from $`\mathrm{}^2()`$.
(g) Alternatively, the bound (8) holds for every $`T`$ if one allows for $`\theta `$-dependent $`C_1`$ and adjusted $`p(\lambda ,\alpha )=p(\lambda ,\alpha ,\theta )`$.

The organization of the article is as follows. In the next section, we recall important properties of the operators $`H_{\lambda ,\alpha ,\theta }`$ and in particular describe the results of Liu and Wen. In Section 3 we establish all the necessary ingredients for an application of the general method of Killip, Kiselev, and Last. Finally, in Section 4, we put everything together and prove Theorem 1.

## 2. Spectra of Periodic Approximants and a Result of Liu and Wen

In this section we discuss the structure of the spectra of the standard periodic approximants to a given $`H_{\lambda ,\alpha ,\theta }`$ and describe recent results of Liu and Wen.

Fix $`\lambda >0`$ and some irrational $`\alpha (0,1)`$. Let $`(a_k)_k`$ be the continued fraction expansion coefficents associated with $`\alpha `$ via (4), and let $`p_k/q_k`$ be the rational approximants to $`\alpha `$. Following , we define the matrices $`M_k(E)=M_k(E;\lambda ,\alpha )`$, for $`k1`$, by

(9)
$$M_k(E)=M(E,q_k,\lambda ,\alpha ,0),$$

where, for $`n`$,

(10)
$$M(E,n,\lambda ,\alpha ,\theta )=\{\begin{array}{cc}T(E,n,\lambda ,\alpha ,\theta )\times \mathrm{}\times T(E,1,\lambda ,\alpha ,\theta )& \text{ if }n1\hfill \\ \left(\begin{array}{cc}1& \hfill 0\\ 0& \hfill 1\end{array}\right)& \text{ if }n=0\hfill \\ T(E,n+1,\lambda ,\alpha ,\theta )^1\times \mathrm{}\times T(E,0,\lambda ,\alpha ,\theta )^1& \text{ if }n1\hfill \end{array}$$

and, for $`m`$,

(11)
$$T(E,m,\lambda ,\alpha ,\theta )=\left(\begin{array}{cc}E\lambda v_{\alpha ,\theta }(m)& \hfill 1\\ 1& \hfill 0\end{array}\right).$$

Thus, $`M_k(E)`$ is the standard transfer matrix associated with the operator $`H_{\lambda ,\alpha ,\theta =0}`$ and the interval $`[1,q_k]`$.

For $`k=1`$ and $`k=0`$, we define

(12)
$$M_1(E)=\left(\begin{array}{cc}1& \hfill \lambda \\ 0& \hfill 1\end{array}\right)\text{ and }M_0(E)=\left(\begin{array}{cc}E& \hfill 1\\ 1& \hfill 0\end{array}\right).$$

Furthermore, we let

(13)
$$t_{(k,p)}(E)=\mathrm{tr}\left(M_{k1}(E)M_k(E)^p\right)$$

and

(14)
$$\sigma _{(k,p)}=\{E:|t_{(k,p)}(E)|2\}.$$

The set $`\sigma _{(k,p)}`$ is the spectrum of a periodic Schrödinger operator whose transfer matrix over one period is given by $`M_{k1}(E)M_k(E)^p`$. Consequently, it consists of a finite number ($`=pq_k+q_{k1}`$, to be precise) of closed intervals (“bands”). In particular, the set $`\sigma _{(k+1,0)}`$ is the spectrum of the periodic Schrödinger operator $`H_{\lambda ,p_k/q_k,0}`$ (i.e., its potential results from $`\lambda v_{\alpha ,0}`$ by replacing $`\alpha `$ by the rational approximant $`p_k/q_k`$); compare .

We recall two key results :

(15)
$$M_{k+1}(E)=M_{k1}(E)M_k(E)^{a_{k+1}},$$

which implies $`t_{(k+2,0)}=t_{(k,a_{k+1})}`$, and

(16)
$$t_{(k+1,0)}^2+t_{(k,p)}^2+t_{(k,p+1)}^2t_{(k+1,0)}t_{(k,p)}t_{(k,p+1)}=4+\lambda ^2.$$

We first discuss how the spectrum of $`H_{\lambda ,\alpha ,\theta }`$ is approximated by the sets $`\sigma _{(k,p)}`$. It is easy to see (and was noted in ) that the set $`\sigma (H_{\lambda ,\alpha ,\theta })`$ is independent of $`\theta `$, that is, there exists a compact set $`\mathrm{\Sigma }_{\lambda ,\alpha }`$ such that

$$\sigma (H_{\lambda ,\alpha ,\theta })=\mathrm{\Sigma }_{\lambda ,\alpha }\text{ for every }\theta [0,1).$$

Let us define the following three types of bands:

| $`(k,I)`$-type band | : | a band of $`\sigma _{(k,1)}`$ which is contained in a band of $`\sigma _{(k,0)}`$, |
| --- | --- | --- |
| $`(k,II)`$-type band | : | a band of $`\sigma _{(k+1,0)}`$ which is contained in a band of $`\sigma _{(k,1)}`$, |
| $`(k,III)`$-type band | : | a band of $`\sigma _{(k+1,0)}`$ which is contained in a band of $`\sigma _{(k,0)}`$. |

These bands are mutually disjoint and they are called spectral generating bands of order k by Liu and Wen. Their union is denoted by $`𝒢_k`$. The following was shown in :

###### Lemma 2.1.

(a) $`\mathrm{\Sigma }_{\lambda ,\alpha }=_k𝒢_k`$.
(b) For $`k>0`$, every spectral generating band of order $`k+1`$ is contained in some spectral generating band of order $`k`$.

From this lemma, we see that for every $`k`$, every energy in the spectrum lies in some spectral generating band of order $`k`$. Thus, we can assign to it a one-sided infinite sequence over the alphabet $`𝒜=\{I,II,III\}`$.

Define the matrices $`T_m`$ by

$$T_m=(t_m(i,j))_{i,j=1,2,3}=\left(\begin{array}{ccc}0& 1& 0\\ a_m+1& 0& a_m\\ a_m& 0& a_m1\end{array}\right).$$

An $`(m1,i)`$-type band generates $`t_m(i,j)`$ bands of $`(m,j)`$ type. Moreover, every spectral generating band of order $`k`$ is associated with a unique word $`i_0i_1\mathrm{}i_k𝒜^{k+1}`$, called its type index.

A central result of , whose proof makes critical use of (15) and (16), establishes lower bounds for the derivative with respect to the energy of the trace function $`t_{(k,p)}`$ for energies from spectral generating bands. To formulate this result, we need the following sequence of matrices:

(17)
$$P_m=(p_m(i,j))_{i,j=1,2,3}=\left(\begin{array}{ccc}0& t_\lambda ^{(a_m1)}& 0\\ a_m/t_\lambda & 0& a_m/t_\lambda \\ a_m/t_\lambda & 0& a_m/t_\lambda \end{array}\right),$$

where

$$t_\lambda =\frac{3}{\lambda 8}.$$

Then, Liu and Wen prove the following result (essentially Proposition 5 of ; the explicit estimate appears in the proof of this proposition).

###### Proposition 2.2.

Let $`\lambda >20`$. For every spectral generating band $`B`$ with type index $`i_0i_1\mathrm{}i_k`$, we have for every $`EB`$,

$$|T^{}(E)|p_1(i_0,i_1)p_2(i_1,i_2)\mathrm{}p_k(i_{k1},i_k),$$

where $`T`$ is the appropriate trace function (i.e., $`t_{(k,p)}`$ if $`B`$ is a band of $`\sigma _{(k,p)}`$).

We now put the above result in a form suitable for our purpose. Denote, for $`k`$,

(18)
$$A_k=\mathrm{\Pi }_{i=1}^ka_i1$$

and

$$x_k(E)=t_{(k+1,0)}(E)\text{ and }\sigma _k=\sigma _{(k+1,0)}.$$

Then we have (cf. )

(19)
$$\mathrm{\Sigma }_{\lambda ,\alpha }=\underset{k}{}(\sigma _k\sigma _{k+1})$$

and the following result for $`x^{}(E)`$:

###### Proposition 2.3.

Let $`\lambda >20`$. Then for every $`k`$ and every $`E\sigma _k`$, we have

$$|x_k^{}(E)|A_k\xi (\lambda )^{k1},$$

where

$$\xi (\lambda )=\left(\frac{\lambda 8}{3}\right)^{1/2}.$$

###### Proof.

This is an immediate consequence of the results above. We note that every band of $`\sigma _k`$ is either a $`(k,II)`$-type band or a $`(k,III)`$-type band and hence generating. Moreover, we see from the specific form of the matrices $`P_m`$ in (17) that we pick up a factor $`(\lambda 8)/3`$ in at least every other step. Finally, for $`\lambda >20`$, we have

$$t_\lambda ^{(a_m1)}a_m$$

and we therefore pick up a factor $`a_m`$ in every step. ∎

## 3. Bounds on Transfer Matrices and Variation of the Phase

In this section we establish the input to KKL theory, which will then be used to prove Theorem 1 in the next section. Namely, we will study derivatives of traces of transfer matrices with respect to the energy, and we will then prove lower bounds for averaged transfer matrix norms. Our approach is combinatorial in nature. Namely, we will use the partition approach to Sturmian sequences, as developed in , to study the local structure of the sequences $`v_{\alpha ,\theta }`$ for general $`\alpha ,\theta `$. In particular, we will exhibit plenty of occurrences of words conjugate to words whose associated trace derivatives we can control, thanks to Propositions 2.2 and 2.3, which is crucial since the trace is invariant with respect to cyclic permutation. This will then allow us to use the results of Section 2 to study the issues at hand.

Let $`M(E,n,\lambda ,\alpha ,\theta )`$ be the transfer matrix as defined in (10). We denote, for $`L1`$,

$$M(E,\lambda ,\alpha ,\theta )_L^2=\underset{n=1}{\overset{L1}{}}M(E,n,\lambda ,\alpha ,\theta )^2+(LL)M(E,L,\lambda ,\alpha ,\theta )^2,$$

where $`M(E,n,\lambda ,\alpha ,\theta )`$ is the usual operator norm of the matrix $`M(E,n,\lambda ,\alpha ,\theta )`$. While this notation conflicts somewhat with the definition in (3), it is standard (and was used in  and previous works). We hope that this does not lead to any real confusion.

As shown in  (using a formula from ), we have the following estimate of this quantity in terms of the derivative of the trace:

(20)
$$\frac{}{E}(\mathrm{tr}M(E,L,\lambda ,\alpha ,\theta ))4M(E,\lambda ,\alpha ,\theta )_{L+1}^3.$$

Thus, combining this estimate with the result given in Proposition 2.2 (resp., Proposition 2.3), we obtain lower bounds on $`M(E,\lambda ,\alpha ,\theta )_L`$ for energies in spectral generating bands (for suitable $`L`$ and, at this point, $`\theta =0`$).

The quantity $`M(E,\lambda ,\alpha ,\theta )_L`$ for $`L0`$ is defined analogously. Clearly, an estimate similar to (20) holds.

Our immediate goal is to study the trace of $`M(E,n,\lambda ,\alpha ,\theta )`$ for general $`\theta `$. To this end, we analyze the local structure of the sequences $`v_{\alpha ,\theta }`$. We show in particular that the traces of the transfer matrices over intervals of length $`q_k`$ exhibit a rather strong invariance property with respect to a variation of the phase.

Let us first recall some combinatorial notions. As general references, we want to mention . Let $`A`$ be a finite set, called the alphabet, and denote by $`A^{},A^{},A^{}`$ the set of finite, one-sided infinite, and two-sided infinite words over $`A`$, respectively. A word $`v`$ is called a subword (or factor) of some word $`u`$ if there are words $`w_1,w_2`$ (possibly empty) such that $`u=w_1vw_2`$. If $`w=w_1\mathrm{}w_m`$ with $`w_iA`$, then the word $`w^R=w_m\mathrm{}w_1`$ is called the reversal of $`w`$. A word $`w`$ with $`w=w^R`$ is called a palindrome. If a word $`w`$ can be written as $`w=uv`$ with words $`u,v`$, we call $`u`$ a prefix of $`w`$ and $`v`$ a suffix of $`w`$. Two words $`w_1,w_2`$ are said to be conjugate if there are words $`u,v`$ such that $`w_1=uv`$ and $`w_2=vu`$. This is equivalent to $`w_2`$ being equal to a cyclic permutation of $`w_1`$. Given any word $`w`$, we denote by $`P_w`$ the set of its finite subwords, and by $`P_w(n)`$ the set of its finite subwords of length $`n`$, $`n`$. Write $`p_w(n)`$ for the cardinality of $`P_w(n)`$; the function $`p_w:`$ is called the complexity function associated with $`w`$. If the word $`w`$ is infinite and uniformly recurrent (i.e., each of its finite subwords occurs infinitely often and with bounded gaps), we define the hull associated with $`w`$ by $`\mathrm{\Omega }_w=\{sA^{}:P_s=P_w\}`$.

In our concrete setting, the alphabet will be given by $`A=\{0,1\}`$ and the word $`w`$ will be given by the restriction of $`v_{\alpha ,0}`$ to $``$, that is,

$$w=v_{\alpha ,0}(1)v_{\alpha ,0}(2)v_{\alpha ,0}(3)\mathrm{}\{0,1\}^{}.$$

The word $`w`$ is called a characteristic Sturmian sequence and its combinatorial properties have been studied extensively. For example, it is well known that its complexity function is given by

(21)
$$p_w(n)=n+1\text{ for every }n.$$

This implies that for every $`n`$, $`P_w(n)`$ contains exactly one word which has two extensions to the right to form a word in $`P_w(n+1)`$; every other word extends uniquely to the right. This special word, let us denote it by $`r_n`$, is called the right-special factor of $`w`$ of length $`n`$ and, for later use, we note the following (Proposition 2.1.23 of ):

(22)
$$r_n=v_{\alpha ,0}(n)\mathrm{}v_{\alpha ,0}(1).$$

That is, the set of right-special factors of $`w`$ coincides with the set of reversals of its prefixes. Also for later use, we note that it follows from Proposition 4.5 of  that

(23)
$$\text{For every }k,v_{\alpha ,0}(1)\mathrm{}v_{\alpha ,0}(q_k2)\text{ is a palindrome.}$$

Moreover, it is also well known that $`w`$ is uniformly recurrent and that for every $`\theta `$, we have $`v_{\alpha ,\theta }\mathrm{\Omega }_w`$. In particular, we have

(24)
$$P_{v_{\alpha ,\theta }}=P_w\text{ for every }\theta .$$

Our goal is to study the sets $`P_w(q_k)`$, $`k`$. One element of $`P_w(q_k)`$ is certainly given by the prefix $`s_k`$ of $`w`$ of length $`q_k`$. As is well known, the words $`s_k`$, $`k`$ obey recursive relations. Since this fact is crucial to our proof, we recall this result briefly:

(25)
$$s_0=0,s_1=0^{a_11}1,s_k=s_{k1}^{a_k}s_{k2}\text{ for }k2.$$

Note that this is where (15) comes from.

Our first observation is very simple:

###### Lemma 3.1.

For $`k2`$, the word $`s_k`$ has suffix $`10`$ if $`k`$ is even and it has suffix $`01`$ if $`k`$ is odd.

###### Proof.

This follows immediately from the recursion (25). ∎

Define $`\overline{}:AA`$ by $`\overline{0}=1`$, $`\overline{1}=0`$. Our next goal is to list all elements of $`P_w(q_k)`$ explicitly. Write $`s_k=s_k^{(1)}\mathrm{}s_k^{(q_k)}`$ with $`s_k^{(i)}A`$, $`1iq_k`$.

###### Lemma 3.2.

For every $`k0`$, the $`q_k+1`$ elements of $`P_w(F_k)`$ are given by

* The $`q_k`$ cyclic permutations of $`s_k`$ which are mutually distinct, and
* the word $`\overline{s_k^{(q_k)}}s_k^{(1)}\mathrm{}s_k^{(q_k1)}`$.

###### Proof.

As a preliminary remark, we note that it follows from (25) that $`w`$ contains the factor $`s_{k1}s_ks_k`$, for every $`k1`$. Namely, for $`k1`$, $`w`$ has the prefix

$`s_{k+3}`$ $`=`$ $`s_{k+2}^{a_{k+3}}s_{k+1}`$
$`=`$ $`\left(s_{k+1}^{a_{k+2}}s_k\right)^{a_{k+3}}s_k^{a_{k+1}}s_{k1}`$
$`=`$ $`\left(\left(s_k^{a_{k+1}}s_{k1}\right)^{a_{k+2}}s_k\right)^{a_{k+3}}s_k^{a_{k+1}}s_{k1},`$
which contains the factor $`s_{k1}s_ks_k`$.

First of all, the claim of the proposition can be verified easily for $`k=0`$ and $`k=1`$. We therefore consider the case $`k2`$ and write $`a`$ for the rightmost symbol of $`s_k`$. Then, using Lemma 3.1, we have the following structure somewhere in $`w`$:

$$\begin{array}{ccc}s_{k1}& s_k& s_k\\ \text{ }|\overline{a}& \text{ }|a& \text{ }|a\end{array}$$

In particular, all the words listed in the assertion of the lemma belong to $`P_w(q_k)`$. Finally, to conclude the proof all we have to show is that the $`q_k`$ cyclic permutations of $`s_k`$ are mutually distinct because by (21) there are only $`q_k+1`$ words in $`P_w(q_k)`$ and the list contains $`q_k+1`$ words which are mutually distinct. Define, for $`k0`$, the height $`h(s_k)`$ of $`s_k`$ by $`h(s_k)=`$ number of $`1`$’s in $`s_k`$. It follows from the definition, (5), and (25) that $`h(s_k)=p_k`$. Since $`p_k`$ and $`q_k`$ are relatively prime for every $`k0`$, we get that for $`k0`$, $`h(s_k)`$ and $`|s_k|`$ are relatively prime. This implies that the cyclic permutations of $`s_k`$ are mutually distinct, for otherwise $`s_k`$ could be written as a power of some shorter word, contradicting the above observation. ∎

We see that there is only one word $`b_k`$ in $`P_w(q_k)`$ which is not a cyclic permutation of $`s_k`$ and it can be described explicitly. In particular, it follows from Lemma 3.1 and Lemma 3.2 that the following holds:

(26)
$$\text{The leftmost symbol of }b_k\text{ is }\{\begin{array}{cc}1\hfill & \text{if }k\text{ is even,}\hfill \\ 0\hfill & \text{if }k\text{ is odd,}\hfill \end{array}$$

and

(27)
$$\text{The rightmost symbol of }b_k\text{ is }\{\begin{array}{cc}1\hfill & \text{if }k\text{ is even,}\hfill \\ 0\hfill & \text{if }k\text{ is odd.}\hfill \end{array}$$

Denote by $`s_k^\theta `$ the word $`v_{\alpha ,\theta }(1)v_{\alpha ,\theta }(2)\mathrm{}v_{\alpha ,\theta }(q_k)`$ and by $`t_k^\theta `$ the word $`v_{\alpha ,\theta }(q_k+1)v_{\alpha ,\theta }(q_k+2)\mathrm{}v_{\alpha ,\theta }(0)`$. We can now state the following combinatorial result:

###### Proposition 3.3.

For every $`\theta `$, we have that $`s_k^\theta `$ is a cyclic permutation of $`s_k`$ for all $`k`$ odd or for all $`k`$ even. The same statement is true for $`t_k^\theta `$.

###### Proof.

Fix $`\theta `$. If $`v_{\alpha ,\theta }(1)=0`$, then by (26), $`s_k^\theta `$ is conjugate to $`s_k`$ for all $`k`$ odd. Similarly, if $`v_{\alpha ,\theta }(1)=1`$, then by (26), $`s_k^\theta `$ is conjugate to $`s_k`$ for all $`k`$ even. A completely analogous argument, using (27), yields the claim for $`t_k^\theta `$. ∎

This has the following immediate consequence for the traces of the transfer matrices. Let

$$x_k(E,\lambda ,\alpha ,\theta )=\mathrm{tr}M(E,q_k,\lambda ,\alpha ,\theta ).$$

###### Proposition 3.4.

For every $`\lambda ,\theta `$, we have that

(28)
$$x_k(E,\lambda ,\alpha ,\theta )=x_k(E,\lambda ,\alpha ,0)\text{ for every }E$$

holds for all $`k`$ odd or for all $`k`$ even. In particular,

(29)
$$\frac{}{E}x_k(E,\lambda ,\alpha ,\theta )=\frac{}{E}x_k(E,\lambda ,\alpha ,0)\text{ for every }E$$

holds for all $`k`$ odd or for all $`k`$ even.

###### Proof.

This follows from Proposition 3.3 and the invariance of the trace of a product with respect to a cyclic permutation of the factors. ∎

Similarly, we have the following result on the left half-line. Let

$$y_k(E,\lambda ,\alpha ,\theta )=\mathrm{tr}M(E,q_k,\lambda ,\alpha ,\theta ).$$

###### Proposition 3.5.

For every $`\lambda ,\alpha ,\theta `$, we have that

(30)
$$y_k(E,\lambda ,\alpha ,\theta )=x_k(E,\lambda ,\alpha ,0)\text{ for every }E$$

holds for all $`k`$ odd or for all $`k`$ even. In particular,

(31)
$$\frac{}{E}y_k(E,\lambda ,\alpha ,\theta )=\frac{}{E}x_k(E,\lambda ,\alpha ,0)\text{ for every }E$$

holds for all $`k`$ odd or for all $`k`$ even.

###### Proof.

We have

$`y_k(E,\lambda ,\alpha ,\theta )`$ $`=`$ $`\mathrm{tr}M(q_k,E,\lambda ,\alpha ,\theta )`$
$`=`$ $`\mathrm{tr}\left(T(q_k+1,E,\lambda ,\alpha ,\theta )^1\times \mathrm{}\times T(0,E,\lambda ,\alpha ,\theta )^1\right)`$
$`=`$ $`\mathrm{tr}\left(T(0,E,\lambda ,\alpha ,\theta )\times \mathrm{}\times T(q_k+1,E,\lambda ,\alpha ,\theta )\right),`$
where in the last step we have used that the determinant is one and hence the trace is invariant with respect to inverting the matrix. By Proposition 3.3 and invariance of the trace with respect to cyclic permutations, for all even $`k`$ or for all odd $`k`$, the right-hand side is equal to $`x_k(E,\lambda ,\alpha ,0)`$ for all $`E`$. ∎

For some fixed $`\theta `$, these results determine the traces $`x_k(\theta )`$ and $`y_k(\theta )`$ in terms of $`x_k`$ for one-half of the possible values of $`k`$, which leaves us with the question of how to investigate the other half. While it is certainly not true that for every $`\theta `$ and every $`k`$, $`x_k(\theta )`$ and $`x_k`$ are equal as functions of $`E`$, we will show that this is true at least for $`k`$ large enough and all but one $`\theta `$.

To this end, we will employ another combinatorial consideration which is based on the partitions of the sequences $`v_{\alpha ,\theta }`$ introduced by D. Lenz and the author in .

###### Definition 3.6.

Let $`k_0`$ be given. A $`(k,\alpha )`$-partition of a function $`f:\{0,1\}`$ is a sequence of pairs $`(I_j,z_j)`$, $`j`$ such that:

* the sets $`I_j`$ partition $``$;
* $`1I_0`$;
* each block $`z_j`$ belongs to $`\{s_k,s_{k1}\}`$; and
* the restriction of $`f`$ to $`I_j`$ is $`z_j`$. That is, $`f_{d_j}f_{d_j+1}\mathrm{}f_{d_{j+1}1}=z_j`$.

Notice that $`d_j`$ is defined implicitly to be the left-hand endpoint of the interval $`I_j`$.

We will suppress the dependence on $`\alpha `$ if it is understood to which $`\alpha `$ we refer. In particular, we will write $`k`$-partition instead of $`(k,\alpha )`$-partition. The sequences $`v_{\alpha ,\theta }`$ have a unique decomposition property which is given in the following lemma (Lemma 3.2.(b) of ).

###### Lemma 3.7.

For every $`k_0`$, every irrational $`\alpha (0,1)`$, and every $`\theta [0,1)`$, there exists a unique $`k`$-partition $`(I_j,z_j)`$ of $`v_{\alpha ,\theta }`$.

Using this partition lemma, we can now continue our study of the local structure of a given $`v_{\alpha ,\theta }`$. Recall that $`b_k`$ denotes the unique word in $`P_w(q_k)`$ which is not conjugate to $`s_k`$. As a preliminary result, we show in the following lemma that $`b_k`$ can occur in $`v_{\alpha ,\theta }`$ only at canonical positions relative to the $`k`$-partition of $`v_{\alpha ,\theta }`$. Namely, whenever the $`k`$-partition yields the occurrence of some $`s_{k1}`$ block followed by some $`s_k`$ block, we saw above that we get an occurrence of $`b_k`$. The following lemma says that these are the only occurrences of $`b_k`$ in $`v_{\alpha ,\theta }`$.

###### Lemma 3.8.

Let $`\alpha (0,1)`$ be irrational and let $`\theta [0,1)`$. If

$$v_{\alpha ,\theta }(m)\mathrm{}v_{\alpha ,\theta }(m+q_k1)=b_k$$

for some $`m`$, then $`[m+1,m+q_k]`$ is an interval $`I_j`$ belonging to the $`k`$-partition of $`v_{\alpha ,\theta }`$. In particular, for every occurrence of $`b_k`$ in $`v_{\alpha ,\theta }`$, we have the following local structure:

$`s_{k1}`$  $`s_k`$ (blocks of $`k`$-partition)
$`b_k`$ (relative position of $`b_k`$)
###### Proof.

Recall that it follows from (21) that there is exactly one factor of length $`q_k1`$ which does not have a unique extension to the right to a factor of length $`q_k`$, and it follows from (22) and (23) that this factor is different from $`v_{\alpha ,\theta }(m+1)\mathrm{}v_{\alpha ,\theta }(m+q_k1)`$. Thus this factor extends uniquely to the right, and we necessarily have $`v_{\alpha ,\theta }(m)\mathrm{}v_{\alpha ,\theta }(m+q_k)=s_k`$. Now we can apply Lemma 3.3 of  which says that this occurrence of $`s_k`$ must correspond to one from the $`k`$-partition. That this $`s_k`$ block must be preceded by an $`s_{k1}`$ block in the $`k`$-partition is forced by $`v_{\alpha ,\theta }(m)`$ being equal to the first letter of $`b_k`$. ∎

###### Proposition 3.9.

Let $`\alpha (0,1)`$ be irrational and let $`\theta [0,1)`$. If $`\theta 1\alpha `$, we have

(32)
$$x_k(E,\lambda ,\alpha ,\theta )=x_k(E,\lambda ,\alpha ,0)\text{ for every }E,kk_0(\theta )$$

and

(33)
$$y_k(E,\lambda ,\alpha ,\theta )=x_k(E,\lambda ,\alpha ,0)\text{ for every }E,kk_1(\theta ).$$

As before, this gives (29) and (31) for the respective $`k`$-ranges.

###### Proof.

Assume first that (32) fails. Our goal is to show

(34)
$$\theta =1\alpha .$$

Fix some $`k`$ and consider the $`k`$-partition of $`v_{\alpha ,\theta }`$. If $`s_k^\theta `$ is conjugate to $`s_k`$, we have $`x_k(E,\lambda ,\alpha ,\theta )=x_k(E,\lambda ,\alpha ,0)`$ for every $`E`$. Conversely, if $`s_k^\theta `$ is not conjugate to $`s_k`$, then it follows from Lemma 3.8 that we must have the following situation:

(35)
$$\begin{array}{cc}\hfill \text{}\text{ }s_{k1}\text{ }& \text{}\text{ }s_k\text{ }\hfill \\ \hfill \mathrm{}1& 2\mathrm{}\hfill \end{array}$$

That is, the site $`1`$ is the right endpoint of $`I_0`$. If (32) fails, then we have the situation depicted in (35) for infinitely many values of $`k`$. In other words, $`v_{\alpha ,\theta }`$ restricted to $`[2,\mathrm{})`$ coincides with $`v_{\alpha ,0}`$ restricted to $`[1,\mathrm{})`$ because as a one-sided infinite word, it has infinitely many $`s_k`$’s as prefixes. Since the phase $`\theta `$ can be recovered uniquely from $`v_{\alpha ,\theta }`$ restricted to a half-line (and $`v_{\alpha ,1\alpha }`$ restricted to $`[2,\mathrm{})`$ coincides with $`v_{\alpha ,0}`$ restricted to $`[1,\mathrm{})`$), we obtain (34).
Assume now that (33) fails. By Lemma 3.8, we have the following situation for infinitely many values of $`k`$:

(36)
$$\begin{array}{ccc}\hfill \text{}\text{ }s_{k1}\text{ }& \hfill \text{}\text{ }s_k\text{ }& \text{}\text{ }s_{k1}\text{ or }s_k\text{ }\hfill \\ \hfill \mathrm{}q_k+1& \hfill \mathrm{}01& 2\mathrm{}\hfill \end{array}$$

Again, the site $`1`$ is the right endpoint of $`I_0`$ and we conclude as above that (34) holds. ∎

Next, we provide power-law lower bounds for $`M(E,\lambda ,\alpha ,\theta )_L`$ for $`\lambda >20`$ and $`\alpha `$. This will be achieved using (20) and the results from the preceding section.

###### Proposition 3.10.

Let $`\lambda >20`$ and $`\alpha `$. Then there exist constants $`C,\zeta >0`$ such that for every $`\theta [0,1)`$ and every $`E\mathrm{\Sigma }_{\lambda ,\alpha }`$, we have

(37)
$$M(E,\lambda ,\alpha ,\theta )_LC|L|^\zeta \text{ for }|L|L_0(\theta ).$$

###### Proof.

For each $`\alpha `$, there is a constant $`B`$ such that the associated sequence $`(q_k)_k`$ obeys

(38)
$$q_k+1B^k\text{ for every }k.$$

Now let $`\lambda >20`$, $`\theta [0,1)\{1\alpha \}`$, and $`k`$ such that

(39)
$$k\mathrm{max}\{k_0(\theta ),k_1(\theta )\},$$

with $`k_0(\theta )`$ and $`k_1(\theta )`$ from (32) and (33). Let $`E\mathrm{\Sigma }_{\lambda ,\alpha }`$. Then, by (19), $`E\sigma _k\sigma _{k+1}`$. Hence, by Proposition 2.3, we have either

$$\frac{}{E}x_k(E,\lambda ,\alpha ,0)A_k\xi (\lambda )^{k1}$$

or

$$\frac{}{E}x_{k+1}(E,\lambda ,\alpha ,0)A_{k+1}\xi (\lambda )^k.$$

In either case, by (18), (20), and (32) (yielding (29)), we obtain

$$4M(E,\lambda ,\alpha ,\theta )_{q_{k+1}+1}^3\xi (\lambda )^{k1}.$$

Thus, if we consider $`q_{k+1}+1Lq_{k+2}+1`$ for $`k`$ obeying (39), we get

$`M(E,\lambda ,\alpha ,\theta )_L`$ $``$ $`M(E,\lambda ,\alpha ,\theta )_{q_{k+1}+1}`$
$``$ $`\left({\displaystyle \frac{1}{4}}\xi (\lambda )^{k1}\right)^{1/3}`$
$``$ $`CB^{\zeta (k+2)}`$
$``$ $`C(q_{k+2}+1)^\zeta `$
$``$ $`CL^\zeta ,`$
with $`B`$ from (38) and, essentially,

$$C=\frac{1}{4^{1/3}}\text{ and }\zeta =\frac{\mathrm{log}\xi (\lambda )}{3\mathrm{log}B}$$

(more precisely, $`\zeta =\frac{\stackrel{~}{k}1}{\stackrel{~}{k}+2}\frac{\mathrm{log}\xi (\lambda )}{3\mathrm{log}B}`$ with $`\stackrel{~}{k}=\mathrm{max}\{k_0(\theta ),k_1(\theta )\}`$, and we obtain an exponent $`\zeta `$ which is arbitrarily close to $`\frac{\mathrm{log}\xi (\lambda )}{3\mathrm{log}B}`$ if we increase $`\stackrel{~}{k}`$ suitably).

A completely analogous proof, using (33) (yielding (31)), proves the claimed bound for negative $`L`$.

We are left with the exceptional case $`\theta =1\alpha `$. Since $`v_{\alpha ,1\alpha }`$ is obtained from the sequence $`v_{\alpha ,0}`$ by a unit shift to the right, there is a constant $`C`$, depending only on $`\lambda `$, such that for every $`n1`$, every $`E\mathrm{\Sigma }_{\lambda ,\alpha }`$, we have

$`M(`$ $`E,n,\lambda ,\alpha ,1\alpha )=`$
$`=T(E,n,\lambda ,\alpha ,1\alpha )\times \mathrm{}\times T(E,1,\lambda ,\alpha ,1\alpha )`$
$`=T(E,n1,\lambda ,\alpha ,0)\times \mathrm{}\times T(E,0,\lambda ,\alpha ,0)`$
$`=T(E,n,\lambda ,\alpha ,0)^1(T(E,n,\lambda ,\alpha ,0)\times \mathrm{}\times T(E,1,\lambda ,\alpha ,1))T(E,0,\lambda ,\alpha ,0)`$
$`T(E,n,\lambda ,\alpha ,0)^1M(E,n,\lambda ,\alpha ,0)T(E,1,\lambda ,\alpha ,0)^1^1`$
$`CM(E,n,\lambda ,\alpha ,0),`$
and a similar inequality for $`n0`$. This permits us to deduce (37) for $`\theta =1\alpha `$ from (37) for $`\theta =0`$ with the constant adjusted accordingly. ∎

## 4. Application of KKL Theory

In this section we prove Theorem 1 by combining (37) with the general theory of . To this end, we first recall a result from  and we then show how it applies in our concrete situation.

Consider a discrete, one-dimensional Schrödinger operator

$$Hu(n)=u(n+1)+u(n1)+V(n)u(n)$$

in $`\mathrm{}^2()`$, define the transfer matrices $`M(n,E)`$ as usually, and define for $`\epsilon >0`$, the characteristic length scales $`\stackrel{~}{L}_\epsilon ^\pm (E)`$ (where $`\pm \stackrel{~}{L}_\epsilon ^\pm (E)>0`$) by

$$M(E)_{\stackrel{~}{L}_\epsilon ^\pm (E)}=2M(1,E)^1\epsilon ^1.$$

Denote the spectral measure of the pair $`(H,\delta _1)`$ by $`\mu `$. Then, the following was shown in  (essentially, Theorem 1.5 of that paper):

###### Theorem 2 (Killip, Kiselev, Last).

For any $`T>0`$, $`L>2`$, we have

$$e^{itH}\delta _1_L^2_T>C\mu (\{E:\stackrel{~}{L}_{T^1}^\pm (E)L\}),$$

where $`C`$ is a universal positive constant.

We are now in a position to give the

###### Proof of Theorem 1..

Consider the case of $`\stackrel{~}{L}_{T^1}^+(E)`$ (the other one is analogous). The length scale $`\stackrel{~}{L}_{T^1}^+(E)`$ is determined by

$$M(E,\lambda ,\alpha ,\theta )_{\stackrel{~}{L}_{T^1}^+(E)}=2M(E,1,\lambda ,\alpha ,\theta )^1T.$$

Under the assumptions of the theorem on $`\lambda ,\alpha `$, we get from (37) that for $`T`$ large enough, and every $`E\mathrm{\Sigma }_{\lambda ,\alpha }`$,

$$\stackrel{~}{L}_{T^1}^+(E)^{\frac{\mathrm{log}\xi (\lambda )}{3\mathrm{log}B}}\stackrel{~}{C}T$$

and hence

$$\stackrel{~}{L}_{T^1}^+(E)CT^{\frac{3\mathrm{log}B}{\mathrm{log}\xi (\lambda )}}.$$

This, together with Theorem 2, implies the statement of Theorem 1 and concludes the proof. ∎

Let us comment on some of the remarks that are listed after the statement of Theorem 1. First of all, if the rotation number $`\alpha `$ is equal to $`(\sqrt{5}1)/2`$, then an analog of Proposition 2.2, sufficient for our purpose, was shown in  for $`\lambda >8`$. Since this result is the only place in our proof where we have to put a restriction on $`\lambda `$, our extension to every phase for this particular value of $`\alpha `$ works for every $`\lambda >8`$, rather than every $`\lambda >20`$.

Secondly, if one wants to establish the bound (8) for every $`T`$, one needs (37) for every $`L`$ which can easily be achieved for every fixed $`\theta `$ by adjusting the constants $`C`$ and $`\zeta `$ suitably.


text:
# The model:

### The model:

Calculation of elastic scattering amplitudes and angular distributions was done in the muffin-tin potential approximation by the partial wave expansion technique . To perform these calculations we used programs from the Barbieri/Van Hove Phase Shift package .

For describing inelastic scatterings, we apply the Lindhard dielectric function  together with optical-data models . The response of the medium to a passing electron, giving an energy loss of $`\mathrm{}\omega `$ and a momentum change $`\mathrm{}𝐪`$, is described in those models by a complex dielectric function  $`ϵ(𝐪,\omega )`$. The probability of an energy loss $`\mathrm{}\omega `$ per unit distance travelled by a non-relativistic electron of energy $`E`$, i.e. the differential inverse mean free path, $`\tau (E,\omega )`$ , then reads :

$$\tau (E,\omega )=\frac{1}{\pi Ea_0}_q_{}^{q_+}\frac{dq}{q}Im[ϵ(q,\omega )^1],$$
(1)
where $`a_0`$ is the Bohr radius, and

$$\mathrm{}q_\pm =\sqrt{2m_e}\left(\sqrt{E}\pm \sqrt{E\mathrm{}\omega }\right)$$
(2)
are the maximal and minimal values of the allowed momentum transfer to the solid, and $`m_e`$ denotes the mass of a free electron. The expression for $`\mathrm{}q_\pm `$ assumes that the energy and momentum transfer for an electron moving in the medium is the same as for a free particle in vacuum, i.e. there is no effective mass assumed and no energy gap. Integration of the differential inverse mean free path over the allowed values of $`\omega `$ yields the total inverse inelastic mean free path . It follows from (1) that the only quantity needed to evaluate $`\tau (E,\omega )`$ and $`\lambda _{in}(E)`$ is the dielectric response function $`ϵ(q,\omega )`$. However, most existing data on dielectric response functions were obtained from photon scattering on solids, for which the momentum transfer is zero. The problem is how to predict $`ϵ(q,\omega )`$ with $`q>0`$, knowing only its optical limit ($`q=0`$) . For that purpose a phenomenological optical model was introduced, where $`Im[ϵ(q,\omega )^1]`$ is expressed via the convolution of $`Im[ϵ(q=0,\omega )^1]`$ with some profile function of $`q`$ and $`\omega `$. We used two such models in the following. The optical model by Ashley  includes exchange between the incident electron and the electron in the medium, which is modeled in analogy with the structure of the non-relativistic Møller cross-section. The Tanuma, Powell and Penn model (TPP-2)  was adopted for calculating the differential inelastic mean free path and inelastic mean free path of electrons in the solid. The energy loss function for diamond used in these calculations was taken from Ref. . Diamond shows a dominant peak in its electron-loss function, corresponding to well-defined volume plasmons . This means that the Lindhard dielectric function describes this solid satisfactorily , at least in the resonance region that is far from the minimum energy needed to make an electron-hole pair.

For the sake of simplicity, we only consider the system at low dose rates, and very early in an exposure when the medium is generally still neutral. Thus we neglect long-range Coulomb interactions, in particular electron-electron interactions, and therefore assume that electron trajectories follow straight lines between consecutive scatterings. Consequently, we do not simulate the motion of the holes and neglect any electron-hole interaction. Effects of hole propagation, i. e. possible inelastic scatterings of the holes and subsequent ionization, are neglected as well, since these processes do not contribute significantly to the total ionization rate. Due to the low impact energy of primary electrons ($`E250`$ eV), we neglect interactions of the electrons with inner shells of the atoms (cf. ). The interactions with valence electrons are modelled as interactions with the band in the solid. For electrons in the conduction band of kinetic energy lower than $`E_{elast}14`$ eV in Ashley’s, and $`E_{elast}9`$ eV in the TPP-2 model, the inelastic mean free path predicted by these models was very large. Therefore we assumed that below these energies electrons scatter elastically only. Interactions of conduction band electrons with lattice phonons, which are important for low kinetic energies ($`E<10`$ eV) were neglected . Following Ashley’s model and the TPP-2 model, we did not allow any form of electronic recombination, so electrons excited into the conduction band stayed there during the simulations.

Fig. 1 presents the band structure of diamond at $`T=300`$ K, as derived from . The zero of the energy scale lies in the middle of the band gap. The position of the bottom of the conduction band is $`E_C=E_G/2`$, and the band gap equals $`E_G=5.46`$ eV at $`T=300`$ K. The width of the valence (V) band is estimated to be $`\mathrm{\Delta }E_V23`$ eV . The dispersion relations for electrons near the bottom of the conduction band and for holes near the top of the valence band are given by the quadratic approximation . In order to keep consistency with the Lindhard dielectric formulation, we assume henceforth the effective mass of the electrons in the conduction band to be equal to the mass of a free electron in vacuum, $`m_e`$. For simplicity, we also assume that the mass of the hole is direction-independent, and equals to the effective mass, $`m_V=0.8m_e`$, where $`m_V^3`$ is the determinant of the valence band mass tensor . In contrast with the free-electron-gas model of the insulator band , the kinetic energies of free electrons are measured with respect to the bottom of the conduction band (see Fig. 1). Following Ref. , the number of electrons or holes in an infinitesimal energy interval $`(E,E+dE)`$ for an insulator is:

$$dn(T)=dE\left\{g_C(E)\overline{\mathrm{\Theta }}(E>E_C)+g_V(E)\overline{\mathrm{\Theta }}(E<E_V)\right\}f(E,\mu ,T),$$
(3)
where $`g_{C,V}(E)`$ corresponds to the density of levels in conduction or valence band,

$$g_{C,V}(E)=\sqrt{2EE_{C,V}}\frac{m_{C,V}^{3/2}}{\mathrm{}^3\pi ^2},$$
(4)
and $`f(E,\mu ,T)`$ describes the Fermi-Dirac distribution of electrons, depending on the chemical potential $`\mu `$. The function $`\overline{\mathrm{\Theta }}`$ is the step function.

### Electron trajectories:

The trajectory of electrons inside a diamond bulk was simulated in the following way: the probability that a free electron of kinetic energy $`E`$ travels a distance $`x`$ in a solid without being scattered decreases exponentially, $`P(E,x)=exp(x/\lambda _T(E))/\lambda _T(E)`$, where $`\lambda _T(E)`$ is the total electron mean free path:

$$\frac{1}{\lambda _T(E)}=\frac{1}{\lambda _{el}(E)}+\frac{1}{\lambda _{in}(E)}.$$
(5)
For a given (probability) distribution $`f(y)`$ of $`y`$, we define a normalized probability distribution, $`P(y)`$, in the form:

$$P(y)=\frac{_{y_{min}}^yf(y^{})𝑑y^{}}{_{y_{min}}^{y_{max}}f(y^{})𝑑y^{}},$$
(6)
where $`y_{min}`$ and $`y_{max}`$ limit the allowed values of $`y`$, $`y_{min}<y<y_{max}`$. It may be easily checked that $`0P(y)1`$. In the MC simulation $`y`$ is then chosen to fulfill the relation: $`r=P(y)`$, where $`r`$ is a random number, $`0<r<1`$.

The path length between two collisions may then be obtained from the normalized probability distribution, and it reads $`x=\lambda _T(E)\mathrm{ln}(r)`$, where $`r`$ is a random number, $`0<r<1`$.

After the electron has travelled a distance x, either elastic or inelastic collisions occur as a stochastic process (probability of collision $`\frac{\sigma _{el(in)}}{\sigma _{el}+\sigma _{in}}`$). In elastic collisions, the electron preserves its energy, but it is scattered at an angle $`\mathrm{\Theta }_{el}`$. The value of $`\mathrm{\Theta }_{el}`$ is obtained from the normalized probability distribution (6) obtained from the differential elastic cross section $`\frac{d\sigma _{el}}{d\theta }(E)`$. The momentum transfer to the lattice is neglected at this point. For an inelastic collision the situation gets more complicated. The scattered electron loses part of its energy $`\mathrm{}\omega `$ and changes its momentum by $`\mathrm{}q`$. From the value of $`\mathrm{}q`$ one can estimate the inelastic scattering angle $`\mathrm{\Theta }_{in}`$. MC probabilities $`P(E,\omega ,q)`$ for this process are obtained from the doubly differential cross section $`{\displaystyle \frac{d^2\sigma _{in}}{dqd\omega }}`$.

The energy $`\mathrm{}\omega `$ and momentum $`\mathrm{}q`$ losses can then be transferred to excite a secondary electron from the valence band into the conduction band. The excited electron leaves a hole with an opposite momentum in the valence band. However, the information whether a secondary electron excitation took place, or not, is not contained implicitly in the inelastic cross section. In order to estimate the number of secondary emissions, we consider two limiting cases here:

Limiting case 1: Let us introduce a probability, $`P(E_0)`$, to excite an electron of initial energy $`E_0`$ from the valence band. The probability is derived from the carrier distribution (3), using (6). In the first case we assume that the probabilities $`P(E,\omega ,q)`$ and $`P(E_0)`$ describe two independent events: (a) a scattering of the primary electron of energy $`E`$ with energy loss $`\mathrm{}\omega `$ and momentum transfer $`\mathrm{}q`$, and (b) a possible excitation of an electron of the initial energy $`E_0`$ from the valence band into the conduction band, if $`E_0+\mathrm{}\omega `$ is greater than $`E_C`$. If, however, $`E_0+\mathrm{}\omega <E_C`$, the excitation cannot take place, and the energy loss of the primary electron is transferred to the crystal lattice.

Limiting case 2: Here we do not allow energy transfer to the lattice. We assume that the energy lost by the primary electron is given to a secondary electron (and a hole). The electron is then removed from the valence band into the conduction band without any energy transfer to the lattice. Events (a) and (b) above are thus no longer independent, and a joint probability $`P(E,\omega ,q,E_0)`$ replaces the product $`P(E,\omega ,q)P(E_0)`$. This approach gives the upper limit for the number of secondary electrons released.

In the first limiting case the momentum conservation imposed on the secondary electron emission process, reduces the number of emitted secondary electrons. The electron excited from valence band has an initial kinetic energy of $`E_{k,sec}^e=\frac{m_V}{m_e}(E_VE_0)`$ obtained from the dispersion relation in the valence band. This corresponds to the kinetic energy of the hole, $`E_{k,sec}^h=E_VE_0`$ in the valence band (see Fig. 2a). The energy $`E_V`$, $`E_V=E_CE_G`$, is the energy at the top of the V band. The final kinetic energy of the excited electron is then $`E_{k,sec}^{}=E_0+\mathrm{}\omega E_C`$. The incoming electron transfers a part of its momentum, $`\mathrm{}𝐪`$, to the secondary electron and the hole. Momentum conservation requires that $`𝐩_{sec}=𝐩_h`$ and $`𝐩_{sec}^{}=𝐩_h\mathrm{}𝐪`$ (see Fig. 2b), where $`𝐩_{sec}`$ and $`𝐩_{sec}^{}`$ are the initial and final momenta of the secondary electron, and $`𝐩_h`$ is the momentum of the hole. The momenta of electron have magnitudes $`𝐩_{sec}\sqrt{2m_eE_{k,sec}}`$, $`𝐩_{sec}^{}\sqrt{2m_eE_{k,sec}^{}}`$, and their orientations are chosen to fulfill the conservation of the momentum, which requires: $`|𝐩_{sec}𝐩_{sec}^{}|<\mathrm{}𝐪<𝐩_{sec}+𝐩_{sec}^{}`$ . This is not always possible at independently chosen $`E_0`$’s and $`(\omega ,q)`$ pairs.

For both secondary emission schemes, we write the energy conservation law for the cloud of electrons in the form :

$$E_{kinet,T}^e+E_{kinet,T}^h=EN_{ion}E_GE_{loss},$$
(7)
where $`E_{kinet,T}^e`$ is the total kinetic energy of electrons in the conduction band, $`E_{kinet,T}^h`$ is the total kinetic energy of holes in the valence band, $`E`$ denotes the impact energy of the primary electron, $`N_{ion}`$ denotes the number of secondary electrons (holes) released, and $`E_{loss}`$ is the energy lost in inelastic collisions without secondary electron excitations. The ionization of the medium continues, until the energies of all excited electrons, including the primary electron, fall below $`E_{elast}`$. After that point, only elastic scatterings occur.

### Numerical results

A set of 2000 MC simulations in diamond was performed. Each cascade was initiated by a single electron. In these simulations, the kinetic energy of the primary electron, measured in respect to the bottom of the conduction band was fixed at $`E=250`$ eV. The starting position of the primary electron at $`t=0`$ fs was at the point x=(0,0,0) of an arbitrarily chosen coordinate system. Afterwards, the electron cloud was created via secondary electron emissions. Motions of the holes were neglected.

The value of the elastic and inelastic cross sections depends on the energy of the electron, and this is changing with time. The elastic and inelastic mean free paths of low energy electrons were computed as described in . Within the energy regime considered here, inelastic collisions were rare compared to elastic interactions, and most of them happened within the first $`10`$ fs following primary electron impact. Within the first $`100`$ fs, there were about $`10000`$-$`35000`$ elastic interactions and about $`10`$-$`50`$ inelastic interactions in an individual cascade simulation, from which $`5`$-$`13`$ resulted in the ionization of the sample. The data show that elastic interactions represent the dominant mechanism for electron propagation in the system.

Spatio-temporal evolution of the cascade was analysed through (a) the number of secondary ionizations, $`N_{ion}(t)`$, (b) the electron range, $`r_{max}(t)`$, i. e. the distance of the most distant electron from the position of the primary electron emission, x=(0,0,0), (c) the kinetic temperature of the free electron gas $`kT(t)`$. It should be noted that the kinetic temperature is a non-equilibrium parameter. Quantities (a)-(c) were averaged over a number of cascades. Figure 4 shows the results.

The average number of ionization events during the first femtosecond was estimated to be $`46`$ based on Ashley’s model and $`58`$ based on the TPP-2 model. This is due to the fact that Ashley’s inelastic mean free path, considered as a function of kinetic energy in the conduction band, is larger than the corresponding TPP-2 inelastic mean free path at all energies, and less ionizations occur within the same time interval. The number of secondary ionizations increased with time, and saturated within about $`100`$ fs with a total of $`510`$ electrons (Ashley), or $`613`$ (TPP-2).

It should be stressed that the ionization rate is very sensitive to the secondary electron emission scheme used in the calculations. However, all predicted ionization rates lie within the limiting values predicted by experiment and other phenomenological models . The lower limit is given by the secondary electron emission yield, $`\delta `$, defined as the ratio of the intensity of secondary electrons emitted from the surface to the incident beam intensity. Its experimental value at $`E250`$ eV is $`\delta 3`$ . This means that at least three electrons were created in the cascade. The upper limit of ionization rate, $`N_{ion}^{max}`$ is estimated  as $`N_{ion}^{max}=E/E_{eh}`$, where $`E`$ is the energy of the primary electron, and $`E_{eh}`$ is the average energy needed to create an electron-hole pair in an insulator. For diamond $`E_{eh}17`$ eV , so $`N_{ion}^{max}15`$. Our model gives: without energy transfer to the lattice: $`10`$-$`13`$ secondary electrons, with energy transfer: $`5`$-$`6`$ secondary electrons.

The average electron range after $`90`$ fs was around $`110130`$ Å. The electron range predicted by Ashley’s model was about 24% larger than in the TPP-2 model at $`90`$ fs. Both models predicted almost linear growth of $`r_{max}`$ as a function of time after $`40`$ fs. This is due to the fact, that after this time low energy electrons are dominant in the sample, and, as a consequence, the dynamics of the cloud is dominated by isotropic elastic scatterings. It should be stressed that for low energy electrons ($`E10`$ eV) phonon coupling becomes important in insulators. Here, we did not take phonon exchanges into account.

The temperature of the electrons drops as the cascade evolves. The temperature characteristics of the electron cloud in Ashley’s and the TPP-2 model are similar. In both models they are quite insensitive to the different limiting cases of secondary electron emission. The curves show rapid decrease within $`1`$ fs after the primary electron emission. This is due to the fact that most ionizations occured during the first femtosecond, and subsequent expansion cools the system effectively. The rate of cooling drops later. Temperature data were used to check that in all simulations the energy was conserved in the system.

Energy distribution of secondary electrons. The positions and velocities of electrons recorded at times, $`t=1,\mathrm{\hspace{0.17em}10},\mathrm{\hspace{0.17em}90}`$ fs were collected from all cascades generated and put into one file. Using these data, histograms of energy distributions were obtained (Fig. 5), $`N(E)/N`$, at these time points ( $`N(E)=_{i=1}^{2000}N_i(E)`$, and $`N_i(E)`$ was the number of electrons found in the energy interval $`(E,E+\mathrm{\Delta }E)`$ for the ith cascade ). The distributions were normalized to the total number of electrons, $`N=_EN(E)`$.

As expected, the energy distribution of secondary electrons shows that the number of low energy electrons increases with time. At $`10`$ fs around $`3080`$% of electrons have energies lower than $`10`$ eV. At $`90`$ fs there are mainly low energy electrons (more than $`50`$% of electrons of energy, $`E<10`$ eV) in the sample. Average total energy losses to the lattice, $`E_{loss}`$, are $`106`$ eV in Ashley’s model and $`95`$ eV in Tanuma’s model within $`100`$ fs for the first secondary electron emission scheme, i.e. when energy transfer to the lattice is allowed.

Local electron density and the Debye length. Using the cumulative results from all cascades, we calculated the average local electron density, $`\rho (r)`$, and the local electron temperature in spherical shells of the width $`1`$ Å, centered at the center of mass of the cumulative electron cloud. Using these data, we calculated the Debye length, $`\lambda _D(r)`$ , in each shell. The results recorded at times, $`t=1,\mathrm{\hspace{0.17em}10},\mathrm{\hspace{0.17em}90}`$ fs, are shown in Figs. 6-9. The results show that the Debye length was longer than the radial size of the electron cloud. This was also the case at the highest density peak. There is an approximate exponential relation between the Debye length, $`\lambda _D(r)`$, and the radial size of the electron cloud, $`r`$, $`\mathrm{ln}\lambda _D(r)=ar+b`$, which in this case implies that $`\lambda _D(r)>r`$. Therefore the long-range Coulomb field is not shielded within this cloud, and the electron gas generated does not represent a plasma in a single impact cascade triggered by an electron of $`E250`$ eV energy. This is important as it justifies the independent-electron approximation used in the model. Moreover, crude estimations of the ionization fraction show that the ionization level needed to create an Auger electron plasma in diamond will be reached with a dose of $`2\times 10^5`$ impact X-ray photons per Å2 if these photons arrive before the cascade electrons recombine.

Spatial distribution of secondary electrons. In order to analyse the spatial distribution of secondary electrons, positions of all electrons generated in $`2000`$ cascading processes were collected (at times, $`t=1,\mathrm{\hspace{0.17em}10},\mathrm{\hspace{0.17em}90}`$ fs). The $`40`$ Å wide cross sections of the collective clouds are shown in Figs. 10-13. They were obtained for different optical models and different secondary electron emission schemes. In all models, most of the electrons lay within a sphere of radius $`50`$ Å at $`1`$ fs, $`70`$ Å at $`10`$ fs and $`150`$ Å at $`90`$ fs with respect to the origin of the cascade. It was found that the center of mass of the cloud was moved $`7`$ Å from the starting point in the direction of the electron impact with $`250`$ eV electrons. The shifts in other directions were $`1`$ Å.

The sphericity tensor, $`S^{ab}={\displaystyle \frac{\underset{i=1}{\overset{N}{}}r_i^ar_i^b}{_{i=1}^Nr_i^2}}`$, calculated for the collective cloud was almost diagonal. Non-diagonal elements were $`101000`$ smaller than the diagonal ones. The sphericity tensor was diagonalized, and the sphericity parameter, $`S`$, was obtained from its eigenvalues: $`S=3/2(\lambda _2+\lambda _3)`$, where $`\lambda _1\lambda _2\lambda _3`$. For all emission schemes $`S`$ lay within the interval $`0.830.89`$ at $`1`$ fs, $`0.950.97`$ at $`10`$ fs and $`0.980.99`$ at $`90`$ fs. This implies that the spatial distribution of the collective cloud was practically isotropic, since $`S1`$.

A small anisotropy appeared only at $`1`$ fs. It manifested distinctly in the cloud generated in the TPP-2 model when no energy loss to the lattice was allowed (cf. Fig. 13). The distribution was cylindrically symmetric along the primary impact vector but it shrank along the direction of the primary impact. This was due to the fact that at that time scale the primary electron was still much faster than the secondary electrons (cf. Fig. 5a). At later times scales the sample got dominated by low energy electrons which scattered more and more isotropically. This was reflected in the progressing isotropy in the distribution.

### Conclusions

The mechanism of ionization in solids by low energy electrons ($`10`$-$`500`$ eV) is not trivial. When the electron wavelength is comparable to distances between atoms in a solid, the Born approximation breaks down and a quantum mechanical description of the electron-atom interactions is necessary. Such a model has recently been published  on the basis of a free-electron gas approximation for the electronic structure of the solid. While this approach is appropriate for the description of metallic substances, it is inadequate for insulators. Insulators have a closed valence band and the band gap between the valence band and the conduction band is large. Here we describe simulations where the band structure of the insulator (diamond) was explicitly included. Electron impact cascades were initiated by a single electron of $`250`$ eV, corresponding to the Auger energy of carbon. The results of Monte Carlo simulations show the spatio-temporal evolution of the ionization cascades with the number of secondary electrons emitted, the range of emitted electrons, their kinetic temperature, energy distribution and the plasma characteristics in terms of the Debye length. The following conclusions can be drawn from this study:

\- Ionization rates

The average number of ionization events during the first femtosecond was between $`48`$ depending on the secondary electron emission scheme used in the calculations. Saturation was reached within $`100`$ fs with a total of $`513`$ in the various models. These ionization rates are lower than rates obtained earlier with the Fermi electron gas approximation , and agree better with experiments and other predictions . At time scales much longer than $`10`$ fs, very low energy electrons dominate the sample, and a model specific for this energy regime should be worked out.

\- Energy distribution

The average temperature of the electrons dropped rapidly within the first femtosecond following the primary emission. This is due to the fact that a large part of all ionizations occurred during the first femtosecond, and this cooled the system effectively. As expected, the number of low energy electrons increased with time, and very low energy electrons ($`E<10`$ eV) became dominant in the sample after about $`10`$ fs. As a consequence, elastic interactions were frequent in the sample. When energy transfer to the lattice was allowed, cooling was faster, with an average total energy loss to the lattice of around $`100`$ eV within $`100`$ fs in an individual cascade simulation.

\- Spatial distribution

The spatial distribution of secondary electrons is anisotropic towards the primary impact at $`1`$ fs. At longer time scales, as the system cools down and energy is distributed more and more equally among electrons, the spatial distribution of the electrons becomes isotropic. Phonon coupling and inelastic interactions between very low energy electrons ($`E<10`$ eV) and matter were neglected in the model. These effects may influence the transportation of very low energy electrons within the sample.

\- Debye length and Coulomb shielding

The Debye length increases exponentially with the radial size of the electron cloud. The long-range Coulomb field is not shielded within this cloud, and the electron gas generated does not represent a plasma in a single impact cascade triggered by an electron of $`E250`$ eV energy. This is important as it justifies the independent-electron approximation used in the model. Moreover, an analysis of the ionization fraction shows that ionization level needed to create an Auger electron plasma in diamond will be reached with a dose of $`2\times 10^5`$ impact X-ray photons per Å2 if these photons arrive before the cascade electrons recombine.

\- Implications

The results can be used to estimate damage by low energy electrons in diamond and other carbon-based covalent compounds. The Monte-Carlo code may be adopted to simulate multiionization phenomena in systems, ranging from the explosion of atomic clusters to the formation of plasmas. The model could also be used to estimate ionization rates and the spatio-temporal characteristics of secondary electron cascades in biological substances. Based on numerical figures for diamond, the maximal diameter of the secondary electron cloud reaches about $`100`$ Ångstrom within the first femtosecond with about $`4`$-$`6`$ ionization events in this period. The diameter of this cloud is $`2`$-$`3`$ times larger at than the diameter of a simple protein molecule. At $`10`$ fs, ionization is almost complete with an average of $`5`$-$`11`$ electrons released per cascade. At this point the maximal diameter of the cascade is about $`150`$ Ångstrom, i.e. the size of a large multienzyme complex. At $`100`$ femtosecond, the diameter of the secondary electron cloud reaches the size of a small virus particle ($`300`$-$`400`$ Ångstrom). The overall number of ionization events does not change appreciably from this point onwards as the electrons have cooled down and interact elastically with the sample. Note that the density of proteins is about $`3`$ times smaller than that of diamond, and so the total number of interactions and ionization rates are expected to be lower in proteins than in diamond. The figures above overestimate ionization within the volumes of small biological samples at these time points. Further studies will be necessary to quantitate expected differences.

One should also point out that there are certain practical applications for diamond in which an understanding of the damage mechanisms is important. For instance, diamond and other wide band gap materials are considered as potential cold electron emitters. Due to its low or negative electron affinity, diamond irradiated by electrons shows a high secondary electron emission yield . This offers the possibility for applying diamond for signal amplification, e. g. in scanning electron microscopy. The high radiation resistivity of diamond makes it an ideal candidate to work as a semiconductor in high-radiation environments , as, for instance, in particle detectors.

## Acknowledgments

We are grateful to Gyula Faigel, Zoltan Jurek, Michel A. van Hove, Richard W. Lee, Richard London, Leszek Motyka and Edgar Weckert for discussions. This research has been supported in part by the Polish Committee for Scientific Research with grant No. 2 P03B 05119, the EU-BIOTECH Programme and the Swedish Research Councils. A. S. was supported by a STINT distinguished guest professorship. B. Z. was supported by the Wenner-Gren Foundations.