text:
# 1 Introduction

## 1 Introduction

In 1998, Neuberger proposed the nested conjugate gradient  for solving the propagator of the overlap-Dirac operator

$`D=m_0\left(1+\gamma _5{\displaystyle \frac{H_w}{\sqrt{H_w^2}}}\right),`$ (1)
with the sign function $`\text{sgn}(H_w)H_w(H_w^2)^{1/2}`$ approximated by the polar approximation

$`S(H_w)={\displaystyle \frac{H_w}{n}}{\displaystyle \underset{l=1}{\overset{n}{}}}{\displaystyle \frac{b_l}{H_w^2+d_l}}H_wR^{(n1,n)}(H_w^2),`$ (2)
where $`H_w=\gamma _5D_w`$, and $`D_w`$ is the standard Wilson-Dirac operator plus a negative parameter $`m_0`$ ($`0<m_0<2`$), and the coefficients $`b_l`$ and $`d_l`$ are:

$`b_l=\mathrm{sec}^2\left[{\displaystyle \frac{\pi }{2n}}\left(l{\displaystyle \frac{1}{2}}\right)\right],d_l=\mathrm{tan}^2\left[{\displaystyle \frac{\pi }{2n}}\left(l{\displaystyle \frac{1}{2}}\right)\right].`$
In principle, any column vector of $`D^1=D^{}(DD^{})^1`$ can be obtained by solving the system

$`DD^{}Y=m_0^2[2+\gamma _5S(H_w)+S(H_w)\gamma _5]Y=\text{1I}`$ (3)
with conjugate gradient, provided that the matrix-vector product $`S(H_w)Y`$ can be carried out. Writing

$`S(H_w)Y={\displaystyle \frac{H_w}{n}}{\displaystyle \underset{l=1}{\overset{n}{}}}b_lZ^{(l)},`$ (4)
one can obtain $`\{Z^{(l)}\}`$ by solving the system

$`(H_w^2+d_l)Z^{(l)}=Y`$ (5)
with multi-shift conjugate gradient . In other words, each (outer) CG iteration in (3) contains a complete (inner) CG loop (5), i.e., nested conjugate gradient.

Evidently, the overhead for the nested conjugate gradient is the execution time for the inner conjugate gradient loop (5) as well as the memory space it requires, i.e., $`(2n+3)`$ large vectors, each of $`12N_{site}`$ double complex numbers, where $`N_{site}`$ is the number of sites, and 12 = 3 (color) $`\times `$ 4 (Dirac) is the degrees of freedom at each site for QCD. The memory storage becomes prohibitive for large lattices since $`n`$ is often required to be larger than 16 in order to achieve a reliable approximation for the sign function. To minimize the memory storage for $`\{P^{(l)},Z^{(l)}\}`$, Neuberger  observed that one only needs the linear combination $`_{l=1}^nb_lZ^{(l)}`$ rather than each $`Z^{(l)}`$ individually. Since $`Z_i^{(l)}`$ and its conjugate vector $`P_i^{(l)}`$ at the $`i`$-th iteration of the inner CG are linear combinations of their preceedents $`\{P_j^{(l)},Z_j^{(l)},j=0,\mathrm{},i1\}`$ in the iteration process, thus it is possible to obtain their updating coefficients $`\{\alpha _i^{(l)},\beta _i^{(l)}\}`$ in the first pass, and then use them to update the sum $`_{l=1}^nb_lZ^{(l)}`$ successively in the second pass, with memory storage of only 5 vectors, independent of the degree $`n`$ of the rational polynomial $`R^{(n1,n)}`$.

At first sight, the double pass algorithm seems to be slower than the single pass algorithm. However, in the test run (with $`SU(2)`$ gauge field on the $`8^3`$ lattice), Neuberger found that the double pass actually ran faster by $`30\%`$ than the single pass, and remarked that the speedup most likely reflects the cache usage in the testing platform, the SGI O2000 (with four processors, each with 4MB cache memory).

In this paper, we analyze the number of floating point operations ($`F_2`$) for the double pass algorithm, and show that it is independent of the degree $`n`$ of $`R^{(n1,n)}`$, provided that the number of lattice sites ($`N_{site}`$) is much larger than the number of iterations ($`L_i`$) of the CG loop. The last condition is amply satisfied even for a small lattice (e.g., $`N_{site}=8^3\times 24=12288`$), since $`L_i`$ is usually less than 1000 (after the low-lying eigenmodes of $`H_w^2`$ are projected out). On the other hand, the number of floating point operations ($`F_1`$) for the single pass is a linearly increasing function of $`n`$. It follows that there exists a threshold $`n_F`$ such that $`F_2F_1`$ for $`nn_F`$, where the value of $`n_F`$ depends on the implementation of the algorithms ($`n_F59`$ for our codes). Corresponding to the number of floating point operations, we also obtain the formulas for the CPU times ($`T_1`$ and $`T_2`$) for the single and the double pass algorithms. Further, we show that there exists a threshold $`n_T`$ such that the double pass is faster than the single pass111 In this paper, we only consider the (faster) single pass algorithm in which the vectors $`P^{(l)}`$ and $`Z^{(l)}`$ ($`l=2,\mathrm{},n`$) are not updated after $`Z^{(l)}`$ converges. for $`n>n_T`$, where $`n_T1225`$ (for most platforms), which is quite smaller than the threshold $`n_F59`$ for the number of floating point operations. By timing the speed of each subroutine, we can account for the extra slow down in the single pass algorithm, which is unlikely to be eliminated, due to the memory bandwidth, a generic weakness of any computational system. Thus, in general (for most vector or superscalar machines), one may find that the double pass is faster than the single pass, for $`n>n_T1225`$. This explains why in Neuberger’s test run, even at $`n=32`$, the double pass is already faster than the single pass by 30%. In fact, we find that DEC alpha XP1000 and IBM SP2 SMP also have $`30\%`$ speedup at $`n=32`$ (see Table 4), which agrees with the theoretical estimate (34) using the CPU time formulas (29)-(30).

Nevertheless, the most interesting result is that the speed of the double pass algorithm is almost independent of the degree $`n`$. This implies that the matrix-vector product $`(H_w^2)^{1/2}YR^{(n1,n)}(H_w^2)Y`$ can be approximated to very high precision with sufficiently large $`n`$, without noticeably extra costs.

The outline of this paper is as follows. In section 2, we outline the single and double pass algorithms for the iteration of the (inner) CG loop (5), and analyze their major differences. In section 3, we estimate the number of (double precision) floating point operations as well as the CPU time for the single and double pass algorithms respectively, and show that there exists a threshold $`n_T`$ such that the double pass is faster than the single pass for $`n>n_T`$. In section 4, we perform some tests. In section 5, we conclude with some remarks.

## 2 The single and double pass algorithms

In the section, we outline the single and the double pass algorithms for the inner CG loop (5), and point out their major differences.

For the single pass algorithm, with the input vector $`Y`$, we initialize the vector variables $`\{Z^{(l)},P^{(l)}\},R,A,B`$ and the scalar variables $`\alpha ,\beta ,\{\gamma ^{(l)}\}`$ as

$`Z_0^{(l)}=0,P_0^{(l)}=Y,l=1,\mathrm{},n`$
$`R_0=Y,`$
$`\alpha _1=1,`$
$`\beta _0=0,`$
$`\gamma _1^{(l)}=\gamma _0^{(l)}=1,l=1,\mathrm{},n`$
Then we iterate ($`j=0,1,\mathrm{}`$) according to:

$`A_j`$ $`=`$ $`H_wP_j^{(1)}`$ (6)
$`B_j`$ $`=`$ $`H_wA_j+d_1P_j^{(1)}=(H_w^2+d_1)P_j^{(1)}`$ (7)
$`\alpha _j`$ $`=`$ $`{\displaystyle \frac{R_j|R_j}{P_j^{(1)}|B_j}}`$ (8)
$`R_{j+1}`$ $`=`$ $`R_j\alpha _jB_j`$ (9)
$`\beta _{j+1}`$ $`=`$ $`{\displaystyle \frac{R_{j+1}|R_{j+1}}{R_j|R_j}}`$ (10)
$`P_{j+1}^{(1)}`$ $`=`$ $`R_{j+1}+\beta _{j+1}P_j^{(1)}`$ (11)
$`Z_{j+1}^{(1)}`$ $`=`$ $`Z_j^{(1)}+\alpha _jP_j^{(1)}`$ (12)
together with the following updates for $`l=2,\mathrm{},n`$:

$`\gamma _{j+1}^{(l)}`$ $`=`$ $`{\displaystyle \frac{\gamma _j^{(l)}\gamma _{j1}^{(l)}\alpha _{j1}}{\alpha _j\beta _j(\gamma _{j1}^{(l)}\gamma _j^{(l)})+\gamma _{j1}^{(l)}\alpha _{j1}(1+\alpha _j(d_ld_1))}}`$ (13)
$`P_{j+1}^{(l)}`$ $`=`$ $`\gamma _{j+1}^{(l)}R_{j+1}+\beta _{j+1}\left({\displaystyle \frac{\gamma _{j+1}^{(l)}}{\gamma _j^{(l)}}}\right)^2P_j^{(l)}`$ (14)
$`Z_{j+1}^{(l)}`$ $`=`$ $`Z_j^{(l)}+\alpha _j{\displaystyle \frac{\gamma _{j+1}^{(l)}}{\gamma _j^{(l)}}}P_j^{(l)}.`$ (15)
The loop terminates at the $`i`$-th iteration if $`\sqrt{R_{i+1}|R_{i+1}/Y|Y}`$ is less than the tolerance (tol).

Since we are only interested in the linear combination $`_{l=1}^nb_lZ_{i+1}^{(l)}`$, in which each $`Z_{i+1}^{(l)}`$ can be expressed in terms of $`\{R_j,j=0,\mathrm{},i\}`$, thus we can write

$`{\displaystyle \underset{l=1}{\overset{n}{}}}b_lZ_{i+1}^{(l)}={\displaystyle \underset{j=0}{\overset{i}{}}}c_jR_j`$ (16)
where $`c_j`$ can be derived as

$`c_j={\displaystyle \underset{m=0}{\overset{ij}{}}}\left[\alpha _{j+m}\delta _m\left(b_1+{\displaystyle \underset{l=2}{\overset{n}{}}}b_l{\displaystyle \frac{\gamma _{m+j+1}^{(l)}\gamma _{m+j}^{(l)}}{\gamma _j^{(l)}}}\right)\right],`$ (17)
with

$`\delta _m=\{\begin{array}{cc}_{k=1}^m\beta _{j+k},\hfill & \text{for }m>0,\hfill \\ 1,\hfill & \text{for }m=0.\hfill \end{array}`$ (20)
Therefore, the r.h.s. of (16) can be evaluated with the CG loop (6)-(11), requiring only the storage of 5 large vectors ($`A,B,R,P^{(1)},T=c_jR_j`$), provided that the coefficients $`\{c_j,j=0,\mathrm{},i\}`$ are known. However, from (17), the determination of $`c_j`$ at any $`j`$-th iteration requires some values of $`\{\alpha \}`$, $`\{\beta \}`$ and $`\{\gamma \}`$ which can only be obtained in later iterations. Thus we have to run the first pass, i.e., the CG loop (6)-(11), to obtain all coefficients of $`\{\alpha \}`$, $`\{\beta \}`$ up to the convergence point $`i`$, and then compute all $`\{c_j,j=0,\mathrm{},i\}`$ according to (17) and (13). Finally, we run the second pass, i.e., going through (6), (7), (9), (11), and adding $`c_jR_j`$ to the r.h.s. of (16), successively from $`j=0`$ to the convergence point $`i`$.

Evidently, all operations in (6)-(12) and (14)-(15) are proportional to the number of lattice sites $`N_{site}`$ times the number of iterations $`L_i`$. On the other hand, the computations of the coefficients $`\{\gamma \}`$ (13) and $`\{c_j\}`$ (17) do not depend on $`N_{site}`$, but only on $`L_i`$ (up to a small term proportional to $`L_i^3`$). Thus, for $`N_{site}L_i`$, we can neglect the computation of $`\{c_j\}`$ (17), and focus on the major difference between the single pass and the double pass, namely, the number of operations in (14)-(15), which is proportional to $`(n1)N_{site}L_i`$, versus the number of operations in (6), (7), (9), (11) plus the vector update on the r.h.s. of (16), which is proportional to $`N_{site}L_i`$. Obviously, the number of floating-point operations in the single pass is a linearly increasing function of $`n`$, while that of the double pass is independent of $`n`$, thus it follows that the double pass must be faster than the single pass for sufficiently large $`n`$.

In the next section, we estimate the number of floating point operations as well as the CPU time, for the single pass and the double pass respectively. Even though our countings are based on our codes, they serve to illustrate the general features of the single and the double pass algorithms, which are valid for any software implementations and/or machines.

## 3 The CPU time and the number of floating point operations

For our codes, the number of floating point operations for the single pass is

$`F_1=N_{site}L_i[3552+120(n1)p]+N_{site}(288n_{ev}+48n+1776),`$ (21)
while for the double pass is

$`F_2`$ $`=`$ $`6888N_{site}L_i+N_{site}(288n_{ev}1656)`$ (22)
$`+\left[{\displaystyle \frac{L_i^3}{6}}+L_i^2(2n1)+L_i\left(13n{\displaystyle \frac{73}{6}}\right)7n+7\right]q`$
where $`N_{site}`$ is the number of sites of the lattice, $`L_i`$ is the number of iterations of the CG loop, $`n_{ev}`$ is the number of projected eigenmodes of $`H_w^2`$, and $`n`$ is the degree of the rational polynomial $`R^{(n1,n)}`$. In the single pass, Eq. (21), $`(n1)p`$ is the effective number of the $`(n1)`$ updates in (14)-(15), since $`P^{(l)}`$ and $`Z^{(l)}`$ are not updated after $`Z^{(l)}`$ converges. The value of $`p`$ depends on the convergence criteria as well as the rational polynomial $`R^{(n1,n)}`$ and its argument. Similarly, in the double pass, the sum in (17) only includes the terms which have not yet converged at the iteration $`j`$, and the reduction in the number of floating point operations can be taken into account by the factor $`q`$ in (22). (The value of $`q`$ is about 0.95 for convergence up to the zero in the IEEE double precision representation).

Taking into account of different speeds of various floating point operations, we estimate the CPU time for the single pass and the double pass as follows.

$`T_1`$ $`=`$ $`N_{site}L_i[192t_b+72t_c+3288t_d+(48t_b+72t_a)(n1)p]`$ (23)
$`+N_{site}(288n_{ev}t_e+48nt_b+24t_a+108t_c+1644t_d)`$
$`T_2`$ $`=`$ $`N_{site}L_i(240t_b+72t_c+6576t_d)`$ (24)
$`+N_{site}(288n_{ev}t_e+24t_a144t_b+108t_c1644t_d)`$
$`+q\left[{\displaystyle \frac{L_i^3}{6}}+L_i^2(2n1)+L_i\left(13n{\displaystyle \frac{73}{6}}\right)7n+7\right]t_f`$
where $`t_a,t_b,t_c`$ and $`t_d`$ denote the average CPU time per floating-point operation (FPO) for the four different kinds of vector operations (a)-(d) listed in Table 1 respectively, $`t_e`$ the average CPU time per FPO for constructing the complementary vector from the projected eigenmodes of $`H_w^2`$, and $`t_f`$ the time for computing the coefficients (17) in the double pass. Note that setting $`t_a=t_b=t_c=t_d=t_e=t_f=1`$ in (23) and (24) reproduce (21) and (22) respectively.

It should be emphasized that the numerical values of the constants and coefficients in (21)-(24) may vary slightly from one implementation to another, however, the number of different terms and their functional dependences on the variables ($`N_{site}`$, $`L_i`$, $`n`$, $`n_{ev}`$, $`p`$, $`q`$, $`t_a,t_b,t_c,t_d,t_e`$, and $`t_f`$) should be the same for any codes of the single and double pass algorithms.

For the double pass, it is clear that the first term on the r.h.s. of Eq. (22) is the most significant part, since the number of lattice site ($`N_{site}`$) is usually much larger than the number of iterations ($`L_i`$) of the CG loop such that the second and the third terms on the r.h.s. of (22) can be neglected. For example, $`N_{site}=8^3\times 24`$, $`L_i=1000`$, $`n=16`$, $`n_{ev}=32`$ and $`q=0.95`$, then the first term is $`6888N_{site}L_i8.5\times 10^{10}`$, while the sum of the second and the third terms only gives $`2.8\times 10^8`$. Thus we can single out the most significant part of $`F_2`$,

$`F_26888N_{site}L_i,`$ (25)
which comes from the first pass (6)-(11) and the second pass (6), (7), (9), (11) plus the vector update on the r.h.s. of (16). Similarly, for the single pass, the most significant part of $`F_1`$ is the first term on the r.h.s. of (21)

$`F_1N_{site}L_i[3552+120(n1)p],`$ (26)
which comes from the operations in (6)-(12), and (14)-(15).

Evidently, from (26) and (25), $`F_1`$ is a linearly increasing function of $`n`$ while $`F_2`$ is independent of $`n`$. Thus it follows that there exists a threshold $`n_F`$ such that $`F_2<F_1`$ for $`n>n_F`$. From (26) and (25), we obtain the threshold $`n_F`$,

$`n_F=1+{\displaystyle \frac{139}{5p}},`$ (27)
where the value of $`p`$ depends on the convergence criterion for removing $`\{P^{(l)},Z^{(l)}\}`$ from the updating list, as well as the rational polynomial $`R^{(n1,n)}`$ and its argument. For our codes and the tests in the next section, $`p0.48`$, thus we have

$`n_F59.`$ (28)
Assuming $`N_{site}L_i`$, we obtain the most significant parts of the CPU times (23)-(24) as

$`T_1`$ $``$ $`N_{site}L_i[192t_b+72t_c+3288t_d+(48t_b+72t_a)(n1)p],`$ (29)
$`T_2`$ $``$ $`N_{site}L_i(240t_b+72t_c+6576t_d).`$ (30)
Obviously, from (29) and (30), there exists a threshold

$`n_T=1+{\displaystyle \frac{2t_b+137t_d}{(2t_b+3t_a)p}}`$ (31)
such that $`T_2<T_1`$ (the double pass is faster than the single pass) for $`n>n_T`$.

Even though the countings in (23) and (24) are based on our codes (for $`R^{(n1,n)}`$ with argument $`H_w^2`$), the essential features of (23) and (24) should be common to all implementaions of the single and the double pass algorithms. In other words, the numerical coefficients in (29) and (30) may change from one implementaion to another, however, the existence of a threshold $`n_T`$ must hold for any implementation.

Now it is interesting to compare $`n_T`$ with $`n_F`$. From (27) and (31), one immediately sees that $`n_T<n_F`$ if

$`19t_d<11t_a+7t_b`$ (32)
is satisfied222Note that the inequality (32) is more restrictive than $`685t_d<417t_a+268t_b`$..

In practice, it turns out that $`t_a/t_d>2`$ and $`t_b/t_d>3`$ for most systems (Tables 1-2). Thus, $`n_T1225`$, which is quite smaller than $`n_F59`$.

The speedup of the double pass with respect to the single pass (for $`n>n_T`$) can be defined as

$`S={\displaystyle \frac{T_1T_2}{T_2}}`$ (33)
which is estimated to be

$`S{\displaystyle \frac{(3t_a+2t_b)p}{10t_b+3t_c+274t_d}}(nn_T)`$ (34)
where Eqs. (29)-(31) have been used.

In Table 1, we list our measurements of $`t_a,t_b,t_c`$ and $`t_d`$ for four different hardware configurations of Pentium 4, i.e., two different Rambuses of faster/slower (PC1066/PC800) speed, and with/without SSE2 (the vector processing unit of Pentium 4) codes.

Substituting the values of $`t_a`$, $`t_b`$ and $`t_d`$ into (31), we obtain the theoretical estimates for the threshold $`n_T`$,

$`n_T\{\begin{array}{cc}12,\hfill & \text{Pentium 4, PC800, with SSE2},\hfill \\ 22,\hfill & \text{Pentium 4, PC800},\hfill \\ 13,\hfill & \text{Pentium 4, PC1066, with SSE2},\hfill \\ 25,\hfill & \text{Pentium 4, PC1066},\hfill \end{array}`$ (39)
where $`p0.48`$ has been used.

Note that for each hardware configuration in Table 1, the average CPU time per FPO of the simple vector operations (a)-(c) is much longer than that of (d), Wilson matrix times vector. A simplified explanation333 It should be emphasized that the mechansim of the interactions between the CPU and the RAM is a rather complicated process, which is beyond the scope of this paper. is as follows. Since all these four vector operations involve long vectors, the CPU and its cache cannot hold all data at once. Thus it is necessary to transfer the data from/to the memory successively, every time after the CPU completes its operations on a portion of the vectors. However, for any system, the memory bandwidth is limited. Thus, there is a time interval between consecutive sets of data transferring to/from the CPU. Therefore, if the CPU finishes a computation before the next set of data is ready, then it would waste its cycles in idling. Since any one of the vector operations (a)-(c) is rather simple, the CPU finishes a computation at a speed faster than that of transferring data from/to the memory, thus the CPU ends up wasting a significant fraction of time in idling. On the other hand, for the vector operation (d), the number of FPO is much more than that of any one of (a)-(c), thus when the CPU completes its operations on a portion of the vectors, the next set of data might have been ready, so the CPU does not waste much time in the memory I/O. This explains why the average CPU time per FPO of (a)-(c) is much longer than that of (d). Further, this simple picture also explains why turning on SSE2 of Pentium 4 (see Table 1) doubles the speed of (d) but has no speedups for (a)-(c), since the bottleneck of (a)-(c) is essentially due to the memory bandwidth rather than the speed of the CPU.

If the memory bandwidth is the major cause for the inefficiency of the simple vector operations (a)-(c), then using faster memories would increase the speeds of (a)-(c) more significantly than that of (d). From Table 1, we can compare the speedups of these four vector operations as the (slower) PC800 is repalced with (faster) PC1066. We find that the speedup for (a)-(c) is $`27\%`$, but that for (d) is only $`11\%`$. Thus the speedups are consistent with above picture.

Obviously, the inefficiency of vector operations (a)-(c) should exist in any platforms, not only for the Pentium 4 systems. To check this, we measure $`t_a,t_b,t_c`$ and $`t_d`$ for IBM SP2 SMP (Power 3 at 375 Mhz) and DEC alpha XP1000 (21264A at 667 Mhz) respectively. The results are listed in Table 2, which give

$`n_T\{\begin{array}{cc}21,\hfill & \text{DEC alpha XP1000},\hfill \\ 20,\hfill & \text{IBM SP2 SMP}.\hfill \end{array}`$ (42)
Although it is impossible to go through all platforms and measure the values of $`t_a`$, $`t_b`$ and $`t_d`$ individually, it is expected that $`t_a/t_d>1`$ and $`t_b/t_d>1`$ (such that the inequality (32) is amply satisfied) is a common feature of most systems. In other words, we expect that the double pass is faster than the single pass for $`n>n_T1225`$, at least for most platforms.

Recall that in Neuberger’s test run with SGI O2000, at $`n=32`$, the double pass is faster than the single pass by $`30\%`$ . This is not a surprise at all, in view of similar speedups of other systems at $`n=32`$. For example, for IBM SP2 SMP or DEC alpha XP1000, sustituting the values of $`t_a`$, $`t_b`$, $`t_c`$ and $`t_d`$ (from Table 2) into (34), we find that $`S=T_1/T_2130\%`$ at $`n=32`$, which also agrees with the actual measurements given in the next section (see Table 4). Thus, the speedup $`S`$ of the double pass for $`n>n_T`$ with $`n_T`$ quite smaller than $`n_F`$ is a generic feature of any platform, stemming from the fact that the vector operations in the second pass is more efficient than those (14)-(15) in the single pass (i.e., $`t_a>t_d`$ and $`t_b>t_d`$).

Nevertheless, the salient feature of (25) and (30) is that the number of floating point operations and the CPU time for the double pass are almost independent of $`n`$. Thus one can choose $`n`$ as large as one wishes, with only a negligible overhead. For example, for the $`16^3\times 32`$ lattice, with $`L_i=1000`$, $`n_{ev}=20`$, and $`q=0.95`$, the increment of $`T_2`$ from $`n=16`$ to $`n=200`$ is less than $`0.05\%`$. In other words, one can approximate $`(H_w^2)^{1/2}Y`$ (i.e., preserve the chiral symmetry) to any precision as ones wishes, without noticeably extra costs. This is the virtue of Neuberger’s double pass algorithm, which may have been overlooked in the last five years.

## 4 Tests

In this section, we perform several tests on the single and the double pass algorithms, and compare the theoretical thresholds $`n_T`$ (31) and $`n_F`$ (27) with the measured values.

In Table 3, we list the number of floating point operations and the CPU time for computing one column of the inverse of

$`D(m_q)=m_q+(m_0m_q/2)(1+\gamma _5S(H_w)),`$
i.e., $`D^1(m_q)=D(m_q)^{}Y`$, where $`Y`$ is solved from

$`D(m_q)D^{}(m_q)Y=\left\{m_q^2+(m_0^2m_q^2/4)[2+(\gamma _5\pm 1)S(H_w)]\right\}Y=\text{1I}`$ (43)
with multi-mass (outer) conjugate gradient for a set of 16 bare quark masses ($`0.02m_q0.3`$), while the inner CG (5) is iterated with the single pass and the double pass respectively. The tests are performed on the $`8^3\times 24`$ lattice with SU(3) gauge configuration generated by the Wilson gauge action at $`\beta =5.8`$. Other parameters are: $`m_0=1.30`$, $`n_{ev}=32`$ (the number of projected eigenmodes), $`\lambda _{max}/\lambda _{min}=6.207/0.198`$ (after projection), and the tolerances for the outer and inner CG loops are $`1.0\times 10^{11}`$ and $`2.0\times 10^{12}`$ respectively. The total number of iterations $`L_o`$ for the outer CG loop is around $`100103`$, while the average number of iterations for the inner CG loop is $`287`$.

With the formulas (21)-(24), we can estimate the number of floating point operations and the CPU time for computing one column of $`D^1(m_q)`$ for a number $`n_q`$ of bare quark masses. For the number of floating point operations, our results are

$`G_k`$ $`=`$ $`(L_o+n_q)F_k+`$ (44)
$`N_{site}(60L_on_q+84L_o+66n_q)+16L_on_q13L_o+18n_q+2`$
where $`L_o`$ is the number of iterations of the outer CG loop (43), the subscript $`k=1(2)`$ stands for the single (double) pass respectively. Obviously, the most significant part of $`G_k`$ is the first term on the r.h.s. of (44), thus

$`G_k(L_o+n_q)F_k,k=1,2.`$ (45)
Similarly, the most significant part of the CPU time is

$`V_k(L_o+n_q)T_k,k=1,2.`$ (46)
where $`T_1`$ and $`T_2`$ are given in (23)-(24).

In Table 3, the estimated CPU times $`V_1`$ and $`V_2`$ are in good agreement with the measured CPU times (the deviation is always less than 5%). By comparing the CPU times for the single pass and the double pass, we see that the double pass becomes faster than the single pass at $`n13`$, in agreement with the theoretical estimate (39) for $`p=0.48`$, where $`p`$ is obtained by measuring the effective number of the ($`n1`$) vector pairs $`\{P^{(l)},Z^{(l)},l=2,\mathrm{},n\}`$ which are updated before $`Z^{(l)}`$ converges.

Further, comparing $`G_2`$ and $`G_1`$, we see that $`G_1G_2`$ at $`n_F59`$, in agreement with the theoretical estimate (28) for $`p=0.48`$.

Also, in Table 3, the remarkable feature of the double pass algorithm is demonstrated: the number of floating point operations ($`G_2`$) and the CPU time are almost independent of $`n`$. Thus $`n`$ can be increased to 64 or any higher value such that the chiral symmetry is preserved to any precision as one wishes. The chiral symmetry breaking or the error of the rational approximation $`R^{(n1,n)}`$ due to a finite $`n`$ can be measured by

$`\sigma =\underset{Y}{\mathrm{max}}\left|{\displaystyle \frac{W^{}W}{Y^{}Y}}1\right|,W=S(H_w)Y,`$ (47)
which is shown in the last column of Table 3.

To check the theoretical estimates for the threshold $`n_T`$ in (42), we repeat the tests of Table 3 for Pentium 4 (PC800), IBM SP2 SMP, and DEC alpha XP1000 respectively. The results are listed in Table 4. Obviously, in each case, the double pass is faster than the single pass for $`n>2022`$, in good agreement with the theoretical estimates in (42). Further, at $`n=32`$, the speed of the double pass is faster than the single pass by $`25\%`$, $`31\%`$, and $`31\%`$ for these three platforms respectively, compatible with what Neuberger found in his test run with SGI O2000 . Note that for Pentium 4, using SSE2 code increases the speedup of the double pass to $`66\%`$ at $`n=32`$ (see Table 3), thus makes the double pass alogorithm even more favorable for P4 clusters.

At this point, it may be interesting to repeat the tests of Table 3, but replacing the polar approximation (2) with the Zolotarev optimal rational approximation,

$`S_{opt}(H_w)=h_w{\displaystyle \underset{l=1}{\overset{n}{}}}{\displaystyle \frac{b_l^{}}{h_w^2+c_{2l1}^{}}}H_wR_Z^{(n1,n)}(H_w^2),h_w=H_w/\lambda _{min},`$ (48)
where

$`R_Z^{(n1,n)}(H_w^2)={\displaystyle \frac{d_0^{}}{\lambda _{min}}}{\displaystyle \frac{\underset{l=1}{\overset{n1}{}}(1+h_w^2/c_{2l}^{})}{_{l=1}^n(1+h_w^2/c_{2l1}^{})}}={\displaystyle \frac{1}{\lambda _{min}}}{\displaystyle \underset{l=1}{\overset{n}{}}}{\displaystyle \frac{b_l^{}}{h_w^2+c_{2l1}^{}}},`$ (49)
and the coefficients $`d_0^{}`$, $`b_l^{}`$ and $`c_l^{}`$ are expressed in terms of Jacobian elliptic functions  with arguments depending only on $`n`$ and $`\lambda _{max}^2/\lambda _{min}^2`$ ($`\lambda _{max}`$ and $`\lambda _{min}`$ are the maximum and the minimum of the the eigenvalues of $`|H_w|`$). The results are listed in Table 5.

Comparing Table 3 with Table 5, it is clear that for the single pass with $`n<32`$, Zolotarev optimal approximation is better than the polar approximation, in terms of the precision of the approximation ($`\sigma `$). However, for the double pass, the polar approximation seems to be as good as the Zolotarev approximation since the degree $`n`$ can be pushed to a very large value, with negligible extra CPU time. In other words, with the double pass algorithm, it does not matter which rational approximation one uses to compute $`D^1(m_q)`$ in a gauge background. This seems to be a rather unexpected result.

## 5 Concluding remarks

So far, we have restricted our discussions to the sign function with argument $`H_w`$. However, it is clear that the salient features of the double pass algorithm are invariant for other choices of the argument, e.g., improved Wilson operator. In general, the double pass algorithm is a powerful scheme for the matrix-vector product $`R(H^2)Y`$, where $`R`$ can be any rational polynomial $`R`$ with argument $`H^2`$ (positive definite Hermitian operator), not just for $`(H^2)^{1/2}`$.

The virtue of Neuberger’s double pass algorithm is its constancy in speed and memory storage for any degree $`n`$ of the rational approximation, where its constancy in speed is valid under a mild condition ($`N_{site}L_i`$) which can be fulfilled in most cases. Further, the double pass is faster than the single pass even for $`n`$ as small as 12 (Pentium 4), and it is expected that the threshold $`n_T1225`$ for most systems. Thus, it seems that there is not much room left for the single pass algorithm with Zolotarev approximation, unless the number of inner CG iterations is exceptionally large, which could happen if the low-lying eigenmodes of $`H_w^2`$ are not projected out and treated exactly.

Note that $`H_w^2`$ can be tridiagonalized by the conjugate gradient (6)-(11), with the unitary transformation matrix $`U`$ formed by the normalized residue vectors $`\{\widehat{R}_j,j=0,\mathrm{},i\}`$, and the elements of the tridiagonal matrix expressed in terms of the coefficients $`\{\alpha _j,\beta _j,j=0,\mathrm{},i\}`$  (up to the tolerance of the conjugate gradient), i.e.,

$`U^{}H_w^2U𝒯`$ (50)
where

$`U_{kj}={\displaystyle \frac{(R_j)_k}{\sqrt{R_j|R_j}}},`$ (51)
and $`𝒯`$ is a symmetric tridiagonal matrix with nonzero elements

$`𝒯_{jj}`$ $`=`$ $`{\displaystyle \frac{\beta _j}{\alpha _{j1}}}+{\displaystyle \frac{1}{\alpha _j}},`$ (52)
$`𝒯_{j+1,j}=𝒯_{j,j+1}`$ $`=`$ $`{\displaystyle \frac{\sqrt{\beta _{j+1}}}{\alpha _j}},j=0,\mathrm{},i.`$ (53)
Thus, after running the first pass of the CG loop (6)-(11), $`𝒯`$ can be constructed from the coefficients $`\{\alpha _j,\beta _j\}`$, and diagonalized by an orthogonal transformation

$`𝒯=O\mathrm{\Lambda }\stackrel{~}{O}.`$ (54)
Then the matrix-vector product $`(H_w^2)^{1/2}Y`$ can be evaluated as

$`{\displaystyle \frac{1}{\sqrt{H_w^2}}}YUO{\displaystyle \frac{1}{\sqrt{\mathrm{\Lambda }}}}\stackrel{~}{O}U^{}Y={\displaystyle \underset{j=0}{\overset{i}{}}}l_jR_j`$ (55)
where

$`l_j={\displaystyle \underset{m=0}{\overset{i}{}}}O_{jm}{\displaystyle \frac{1}{\sqrt{\lambda _m}}}O_{0m}\sqrt{{\displaystyle \frac{R_0|R_0}{R_j|R_j}}}.`$ (56)
Here the summation on the r.h.s. of (55) is obtained by running the second pass of the CG loop $`\{(\text{6}),(\text{7}),(\text{9}),(\text{11})\}`$, and adding $`l_jR_j`$ to the sum successively from $`j=0`$ to $`i`$.

It is well known that (any positive definite Hermitian matrix) $`H_w^2`$ can be tridiagonalized by Lanczos iteration  as well as the conjugate gradient. The connection between the Lanczos iteration and the conjugate gradient for the tridiagonalization of a positive definite Hermitian matrix has been well established , and both have almost the same performance in practice. In Ref. , the Lanczos approach was proposed for the matrix-vector product $`(H_w^2)^{1/2}Y`$, and its variant (replacing Lanczos iteration with the conjugate gradient) was used in Ref. .

The only difference between the Lanczos (CG) algorithm and Neuberger’s double pass algorithm is the diagonalization of the tridiagonal matrix $`𝒯`$ and the computation of the coefficients $`\{l_j\}`$ (56) in the former versus the computation of the coefficients $`\{c_j\}`$ (17) in the latter. Since the number of floating point operations for the diagonalization of a symmetric tridiagonal matrix $`𝒯`$ is $`3L_i^3`$ (where $`L_i`$ is the number of iterations of the inner CG loop, or the size of $`𝒯`$), it is compatible with that of computing the coefficients $`\{c_j\}`$, i.e., the last term on the r.h.s. of (22). Thus we expect that the performance (speed and accuracy) of Lanczos (CG) algorithm and Neuberger’s double pass algorithm are compatible.

In Table 6, we compare the Lanczos (CG) algorithm with Neuberger’s double pass algorithm, by computing one column of $`D^1(m_q)`$ (for 16 bare quark masses) on the $`16^3\times 32`$ lattice with SU(3) gauge configuration generated by the Wilson gauge action at $`\beta =6.0`$. Other parameters are: $`m_0=1.30`$, $`n_{ev}=20`$ (the number of projected eigenmodes), $`\lambda _{max}/\lambda _{min}=6.260/0.215`$ (after projection), and the tolerances for the outer and inner CG (Lanczos) loops are $`1.0\times 10^{11}`$ and $`2.0\times 10^{12}`$ respectively. The number of iterations for the outer CG loop is $`L_o=347`$, while the average number of iterations for the inner CG loop is $`300`$. Evidently, these seemingly different algorithms have almost the same speed as well as the accuracy ($`\sigma `$).

Thus, for quenched lattice QCD, one has several compatible options to compute the quenched quark propagator

$`(D_c+m_q)^1=(1rm_q)^1[D^1(m_q)r],r={\displaystyle \frac{1}{2m_0}},`$ (57)
even though we have chosen Neuberger’s double pass algorithm to solve $`D^1(m_q)`$ in our recent investigation . Nevertheless, for lattice QCD with dynamical quarks, the quark determinant $`\text{det}D(m_q)`$ could not be computed directly with existing algorithms and computers. If $`\text{det}D(m_q)`$ is incorporated through the dynamics of $`2n`$ pseudofermion fields (where $`n`$ can be regarded as the degree $`n`$ in the rational polynomial $`R^{(n1,n)}`$), then an additional degree of freedom (or the fifth dimension with $`N_s=2n`$ sites) has to be introduced. Thus a relevant question is how to reproduce $`\text{det}D(m_q)`$ accurately with the minimal $`N_s`$. A solution has been presented in Ref. . On the other hand, it would be interesting to see whether there is an algorithm to drive the dynamics of these $`N_s`$ pseudofermion fields such that the cost is almost independent of $`N_s=2n`$.

Acknowledgement

T.W.C. would like to thank Herbert Neuberger and Rajamani Narayanan for a brief but stimulating discussion at the Institute for Advanced Study, during a visit in November 2002. This work was supported in part by the National Science Council, ROC, under the grant number NSC91-2112-M002-025.


format: text

text:
# EPR and ferromagnetism in diluted magnetic semiconductor quantum wells

## Abstract

Motivated by recent measurements of electron paramagnetic resonance (EPR) spectra in modulation-doped CdMnTe quantum wells, \[F.J. Teran et al., Phys. Rev. Lett. 91, 077201 (2003)\], we develop a theory of collective spin excitations in quasi-two-dimensional diluted magnetic semiconductors (DMSs). Our theory explains the anomalously large Knight shift found in these experiments as a consequence of collective coupling between Mn-ion local moments and itinerant-electron spins. We use this theory to discuss the physics of ferromagnetism in (II,Mn)VI quantum wells, and to speculate on the temperature at which it is likely to be observed in n-type modulation doped systems.

Substitution of transition metal elements in a semiconductor lattice often adds local moments zunger  to the system’s low-energy degrees of freedom and can lead to qualitatively new physics. Important progress sstreviews  has recently been achieved in understanding the materials science and physics of (II,Mn)VI and (III,Mn)V ternary compound semiconductors in which Mn has been substituted on a relatively small fraction of the cation sites. When these systems are doped p-type, the Mn ions spontaneously align at low temperatures in both bulk and quantum-well systems with typical ferromagnetic transition temperatures $`1\mathrm{K}`$ in the quantum-well case and $`100\mathrm{K}`$ in the bulk. Although there is broad agreement that ferromagnetism in these systems is due to carrier-mediated interactions between Mn local moments, consensus bookchapter  on the details of this picture is still building as the body of experimental studies on well characterized samples grows. In this connection it is intriguing that ferromagnetism has never been observed when these semiconductors are n-doped caveat1 . Recent inelastic Raman and resistively detected EPR studies of n-doped (Cd,Mn)Te quantum wells by Teran et al. teran  provide important information about the role played by quantum-well electrons in correlating Mn ion local-moment-spin orientations. In particular, these authors have discovered an avoided crossing between well-defined electron spin and Mn-ion spin resonance modes, an effect which demonstrates that the two subsystems can couple collectively.

Quantum-well DMS systems are unusual in that they consist of a quasi-3D Mn ion system coupled to quantum-well electrons that have only two-dimensional (2D) translational degrees of freedom. In this Letter we present a theoretical analysis which completely accounts for the observations of Teran et al., including the size of the avoided crossing gap they see. The theory sheds light on the physics that controls ferromagnetism in DMS quantum wells and on the essential differences between electron- and hole-doped cases. It follows from our theory that the gap discovered experimentally occurs because quantum-well electrons and Mn-ion spins are coupled ferromagnetically in n-doped systems. According to our theory the ferromagnetic transition temperature is controlled by spin-orbit coupling strength, a quantity that can be adjusted in situ by biasing the quantum well rashba  and is much larger in p-doped systems. We conclude that ferromagnetism will occur in n-doped quantum wells, but only at temperatures well below those predicted by mean-field theory and below those which have been studied experimentally. In the following we first describe our theory of collective excitations of the coupled local-moment and quantum-well-electron spin subsystems, and then discuss ferromagnetism in n-doped (Cd,Mn)Te systems.

We consider an n-type quantum well with width $`d`$ (in $`\widehat{z}`$ direction) and one occupied subband. The quantum-well geometry makes it convenient to split the three-dimensional spatial coordinates into $`(𝐫,z)`$ where $`𝐫`$ is the $`x`$-$`y`$-plane projection. When subband mixing is neglected leewidewells  the electron wavefunction separates, $`\mathrm{\Psi }(𝐫,z)=\psi (𝐫)\chi (z)`$, where $`\chi (z)`$ can be chosen to be real. We take advantage of the large ratio between the Mn density and the electron density in the samples studied by Teran et al. and in typical dietlhandbook94  modulation-doped (II,Mn)VI quantum-well systems, replacing the random distribution of Mn local moments by a continuum density $`N_{\mathrm{Mn}}(z)`$. We are interested in the collective excitations in the presence of an external magnetic field $`𝐁=(0,0,B)`$. The Hamiltonian, $`H=H_{\mathrm{kin}}+H_{\mathrm{Zeeman}}+H_{\mathrm{sd}}`$ is the sum of three terms. The kinetic energy of the electrons is given by

$`H_{\mathrm{kin}}`$ $`=`$ $`{\displaystyle d^2r_0^d𝑑z\chi ^2(z)}`$ (1)
$`{\displaystyle \underset{\sigma }{}}\widehat{\psi }_\sigma ^{}(𝐫)\left({\displaystyle \frac{(i\mathrm{}+e𝐀/c)^2}{2m^{}}}\mu \right)\widehat{\psi }_\sigma (𝐫)`$
where $`𝐀`$ is the vector potential. The Zeeman term reads

$$H_{\mathrm{Zeeman}}=\mu _B𝐁d^2r_0^d𝑑z\left[g_\mathrm{e}\widehat{𝐬}(𝐫,z)+g_{\mathrm{Mn}}𝐒(𝐫,z)\right]$$
(2)
where $`\mu _B>0`$ is the electron Bohr magneton. Here,

$$\widehat{𝐬}(𝐫,z)=\frac{1}{2}\chi ^2(z)\underset{\sigma \sigma ^{}}{}\widehat{\psi }_\sigma ^{}(𝐫)\tau _{\sigma \sigma ^{}}\widehat{\psi }_\sigma ^{}(𝐫),$$
(3)
is the electron spin density, $`\tau _{\sigma \sigma ^{}}`$ the vector of Pauli matrices, and $`𝐒(𝐫,z)`$ is the spin density of the Mn system. The ferromagnetic ($`J_{\mathrm{sd}}<0`$) coupling between the electron and local moment sstreviews  spins is described by

$$H_{\mathrm{sd}}=J_{\mathrm{sd}}d^2r_0^d𝑑z𝐒(𝐫,z)\widehat{𝐬}(𝐫,z).$$
(4)
As in our earlier work Koenig00  on bulk DMS ferromagnets, we develop our theory of elementary spin excitations in a language where the conduction band degrees of freedom are integrated out. The local moments with spin $`S=5/2`$ are represented by Holstein-Primakoff (HP) boson with retarded conduction-band-mediated interactions. Because $`g_{\mathrm{Mn}}>0`$, the external magnetic field tends to align the Mn spins along the negative $`\widehat{z}`$ direction. For small fluctuations around this state we can approximate the spin fields in a coherent-state functional-integral representation of the HP-boson partition function by $`S^+\overline{w}\sqrt{2N_{\mathrm{Mn}}(z)S}`$, $`S^{}w\sqrt{2N_{\mathrm{Mn}}(z)S}`$, and $`S^z=N_{\mathrm{Mn}}(z)S+\overline{w}w`$ where the complex variables $`\overline{w},w`$ label boson coherent states. The partition function is formally expressed as $`Z=𝒟(\overline{w}w)\mathrm{exp}(S_{\mathrm{eff}}[\overline{w}w])`$ with the effective action

$`S_{\mathrm{eff}}[\overline{w}w]`$ $`=`$ $`{\displaystyle _0^\beta }𝑑\tau {\displaystyle d^2r_0^d𝑑z\left[\overline{w}_\tau w+g_{\mathrm{Mn}}\mu _BBS^z\right]}`$ (5)
$`\mathrm{ln}detG^1(\overline{w}w).`$
The electron Greens function in this equation may be partitioned into mean-field and fluctuating terms, $`G^1(\overline{w}w)=(G^{\mathrm{MF}})^1+\delta G^1(\overline{w}w)`$, where

$$(G^{\mathrm{MF}})^1=\left(_\tau \frac{\mathrm{}^2\stackrel{~}{}^2}{2m^{}}\mu \right)+\frac{\mathrm{\Delta }+g_\mathrm{e}\mu _BB}{2}\tau ^z,$$
(6)
with $`\stackrel{~}{}=+(ie/\mathrm{}c)𝐀`$, and the fluctuating part is

$`\delta G^1(\overline{w}w)`$ $`=`$ $`{\displaystyle \frac{|J_{\mathrm{sd}}|}{2}}{\displaystyle _0^d}𝑑z\chi ^2(z)`$ (7)
$`\left[\sqrt{2N_{\mathrm{Mn}}(z)S}\left(\overline{w}\tau ^{}+w\tau ^+\right)+\overline{w}w\tau ^z\right].`$
The exchange contribution to the conduction band spin-splitting that appears in $`G^{\mathrm{MF}}`$ is given by $`\mathrm{\Delta }=|J_{\mathrm{sd}}|\overline{N}_{\mathrm{Mn}}S`$ where $`\overline{N}_{\mathrm{Mn}}=_0^d𝑑z\chi ^2(z)N_{\mathrm{Mn}}(z)`$. Expanding to second order in the boson fields and Fourier transforming we can write the quadratic part of the spin-wave action as

$`S_{\mathrm{eff}}[\overline{w}w]`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{m}{}}{\displaystyle \frac{d^2k}{(2\pi )^2}_0^d𝑑z_0^d𝑑z^{}}`$ (8)
$`\overline{w}(𝐤,z,\nu _m)D^1(𝐤,z^{},z,\nu _m)w(𝐤,z^{},\nu _m)`$
where $`\nu _m`$ are the bosonic Matsubara frequencies. The kernel $`D^1(𝐤,z^{},z,\nu _m)`$ in Eq. (8) is obtained from a straightforward calculation which leads to a fermion loop diagram specified by an unwieldy expression that we do not reproduce here. Instead, we immediately specialize to the case of $`𝐤=0`$ probed by the EPR and Raman experiments of Teran et al.. In this limit the kernel is given by the following physically transparent expression:

$`D^1(𝐤=0,z^{},z,\nu _m)`$ $`=`$ $`\left[i\nu _m+g_{\mathrm{Mn}}\mu _BB+{\displaystyle \frac{|J_{\mathrm{sd}}|}{2}}(n_{}n_{})\chi ^2(z)\right]\delta (zz^{})`$ (9)
$`+{\displaystyle \frac{SJ_{\mathrm{sd}}^2\chi ^2(z)\chi ^2(z^{})\sqrt{N_{\mathrm{Mn}}(z)N_{\mathrm{Mn}}(z^{})}}{2}}{\displaystyle \frac{n_{}n_{}}{i\nu _m\mathrm{\Delta }g_\mathrm{e}\mu _BB}},`$
where $`n_\sigma `$ is the (2D) mean-field conduction-electron spin density for spin $`\sigma `$. At $`𝐤=0`$ this expression is valid at any Landau level filling factor. The first term on the right hand side of Eq. (9) represents the energy cost of flipping an individual Mn spin and includes the Knight shift $`K`$ contribution (see below) due to exchange-coupling with band electrons whose spin density depends on the position of a Mn ion within the quantum well. The second term is the correction to the energy cost (at $`𝐤=0`$) which occurs because the quantum-well electron system responds to Mn spin reorientations. Collective excitations are determined by locating zeroes of the kernel determinant. Since uniform spin orientation fluctuations correspond in our continuum theory to $`w_0(z)\sqrt{N_{\mathrm{Mn}}(z)}`$, the energies $`E=\mathrm{}\mathrm{\Omega }`$ of the collective excitations we seek details  are obtained by solving

$$_0^d𝑑z_0^d𝑑z^{}w_0(z)D^1(0,z^{},z,i\nu _m=\mathrm{\Omega })w_0(z^{})=0,$$
(10)
which implies collective excitations at

$$E=\frac{E_{\mathrm{Mn}}+E_\mathrm{e}}{2}\pm \sqrt{\frac{(E_{\mathrm{Mn}}E_\mathrm{e})^2}{4}+K\mathrm{\Delta }}.$$
(11)
The quantities in Eq. (11) all have simple physical interpretations. The expressions $`E_{\mathrm{Mn}}`$ and $`E_\mathrm{e}`$ denote mean-field transition energies in which the magnetic quantum number is increased by one, for Mn and electron spins respectively. Each mean-field excitation energy has two contributions, the Zeeman energy and the exchange energy due to coupling between the two spin subsystems. For the Mn spin, $`E_{\mathrm{Mn}}=E_{\mathrm{Mn}}^Z+K`$ where the Zeeman term is $`E_{\mathrm{Mn}}^Z=g_{\mathrm{Mn}}\mu _BB`$ and the Knight shift is

$$K=\frac{|J_{\mathrm{sd}}|(n_{}n_{})}{2}\frac{\overline{N}_{\mathrm{Mn}}}{N_{\mathrm{Mn}}d}.$$
(12)
$`K/(g_{\mathrm{Mn}}\mu _B)`$ is the mean exchange field experienced by the local moments because of the electron spin polarization. It is analogous to the Knight shift experienced by nuclei in a spin-polarized electron gas in NMR experiments. The average Mn density in this expression, $`N_{\mathrm{Mn}}`$, defined by $`N_{\mathrm{Mn}}d=_0^d𝑑zN_{\mathrm{Mn}}(z)`$, is not equal to the quantity $`\overline{N}_{\mathrm{Mn}}`$ which appeared previously in the mean-field quantum-well electron spin gap; the effective 3D electron density in Eq. (12) involves the quantum-well width and the Mn distribution within the quantum well in a non-obvious way that follows from our analysis.

For the quantum-well electrons, $`E_\mathrm{e}=E_\mathrm{e}^Z+\mathrm{\Delta }`$, the Zeeman term is $`E_\mathrm{e}^Z=g_\mathrm{e}\mu _BB`$. Following the NMR analogy, the mean-field electron spin-splitting $`\mathrm{\Delta }=|J_{\mathrm{sd}}|\overline{N}_{\mathrm{Mn}}S`$ corresponds to the nuclear-polarization induced Overhauser shift in electron spin-resonance frequencies. Because $`g_{\mathrm{Mn}}>0`$, and therefore $`E_{\mathrm{Mn}}^Z>0`$, the mean-field Mn-spin configuration points in the negative $`\widehat{z}`$ direction. For the conduction-band electrons, there is a competition between the Overhauser energy $`\mathrm{\Delta }>0`$ and the Zeeman term $`E_\mathrm{e}^Z<0`$ ($`g_\mathrm{e}<0`$). The Overhauser shift dominates except at very strong magnetic fields, i.e., the mean-field electron spin polarization is along the negative $`z`$ axis.

Eq. (11) is able to account quantitatively for the experiments of Teran et al.. In their sample the avoided crossing is seen in a strong magnetic field near Landau level filling factor $`\nu =3`$, where the electronic state has two occupied majority-spin Landau levels and one minority spin, i.e., $`n_{}n_{}=n/32\times 10^{11}\mathrm{cm}^2`$. For a Mn fraction that is constant across the quantum well we obtain $`N_{\mathrm{Mn}}=\overline{N}_{\mathrm{Mn}}4.4\times 10^{19}\mathrm{cm}^3`$ for the studied sample teran . The width of the quantum well is $`d10\mathrm{n}\mathrm{m}`$. At low temperatures, the Mn spins are fully polarized by the external magnetic field. Using the experimental teran  low-temperature value for $`\mathrm{\Delta }1.65\mathrm{meV}`$, we conclude that $`|J_{\mathrm{sd}}|`$ in (Cd,Mn)Te is $`0.015\mathrm{eV}\mathrm{nm}^3`$, in agreement with Ref. furdyna, . Using these values we conclude that $`K1.5\mu \mathrm{eV}`$ and hence the avoided crossing gap $`\sqrt{K\mathrm{\Delta }}0.05\mathrm{meV}`$, in close agreement with the experimental estimate $`0.03\mathrm{meV}`$. In the quantum Hall regime, the gap’s temperature dependence may be attributed to the temperature-dependent Mn spin-polarization and the consequent temperature dependence of $`\mathrm{\Delta }`$ (see Fig. 1).

At these strong fields the Mn spins polarization is described by the Brillouin function so that $`\mathrm{\Delta }(T)/\mathrm{\Delta }(0)=B_{5/2}(5g_{\mathrm{Mn}}\mu _BB/2k_B(T+T_{\mathrm{AF}}))`$ where the phenomenological parameter $`T_{\mathrm{AF}}=0.18K`$ accounts for the direct antiferromagnetic Mn-Mn coupling teran .

The close agreement between experimentally observed and calculated anticrossing gaps supports our effective model for DMSs in which the low-energy degrees of freedom are local $`S=5/2`$ Mn spins that are exchange coupled to band-electron spins. The avoided crossing establishes that both Mn and band spins behave collectively, and that they are coupled, conditions under which ferromagnetism is expected. The ferromagnetic transition temperature of this quantum-well system can be estimated by applying the mean-field approach that appears to be generally successful when applied tcmft  to bulk (III,Mn)V materials. Assuming uniform Mn doping and particle-in-a-box subband wavefunctions, the mean-field critical temperature for quantum wells tcqw  can be written as

$$k_BT_c^{\mathrm{MF}}=\frac{S+1}{4}\frac{K^{\mathrm{max}}\mathrm{\Delta }}{ϵ_\mathrm{F}}$$
(13)
where $`K^{\mathrm{max}}=n|J_{\mathrm{sd}}|/(2d)`$ is the value of the ‘Knight shift’ when the quantum-well electrons are completely spin polarized and $`ϵ_\mathrm{F}`$ is the paramagnetic 2D electron gas Fermi energy. Using experimental values for $`\mathrm{\Delta }1.65\mathrm{meV}`$, $`K^{\mathrm{max}}4.5\mu \mathrm{eV}`$, and the effective mass $`m^{}=0.107m_\mathrm{e}`$ we find that $`T_c^{\mathrm{MF}}5.6\mathrm{mK}`$.

In three-dimensional DMS ferromagnets, we expect schliemann  mean-field theory to be reliable when the carrier density is smaller than the Mn density and $`\mathrm{\Delta }<ϵ_\mathrm{F}`$. These conditions are both satisfied here. In the quasi-2D case, however, long-wavelength collective excitations have a significant negative impact on tendencies toward long-range magnetic order. In fact, for the model studied here which has no spin-orbit interactions, continuous spin-rotational invariance implies that long-range magnetic order is impossible at finite temperatures in 2D mermin-wagner . Thermal and quantum fluctuations at low temperatures in these quasi-2D ferromagnets are even more important than usual because 2D electron gas properties lead to vanishing spin stiffness kossacki00 ; details . Long-range magnetic order is possible in (II,Mn)VI quantum wells only because spin-orbit coupling favors magnetization orientations perpendicular to the quantum wells.

In symmetric quantum wells spin-orbit interactions are described by the Dresselhaus Hamiltonian, $`H_D=\gamma (\sigma _xk_x+\sigma _yk_y)k_z^2`$. Evaluating the magnetization-orientation dependence of the energy correction due to this term by second-order perturbation theory we find that the easy axis is perpendicular to the quantum well. The collective excitation energy gap is twice the anisotropy energy per Mn spin. We find that

$$E_{\mathrm{so}}=\frac{\gamma ^2k_z^2^2k_\mathrm{F}^2n}{SN_{\mathrm{Mn}}d\mathrm{max}\{\mathrm{\Delta },2ϵ_\mathrm{F}\}},$$
(14)
taking $`k_z^2=(\pi /d)^2`$ for the lowest quantum-well subband, and $`k_\mathrm{F}^2=2m^{}ϵ_\mathrm{F}/\mathrm{}^2`$. When $`E_{\mathrm{so}}`$ is small the Curie temperature will be limited by collective fluctuations schliemann . Taking account of the expected dispersionless 2D spin-wave bands and estimating the Curie temperature as the temperature at which $`S`$ magnons per Mn spin are thermally excited, we predict that $`k_BT_c^{\mathrm{coll}}(S+1/2)E_{\mathrm{so}}`$. Plugging in $`\gamma 12\mathrm{meV}\mathrm{nm}^3`$ for the Dresselhaus coefficient in the CdTe conduction band cardona , this implies that $`T_c^{\mathrm{coll}}0.4\mathrm{mK}`$ in the sample studied by Teran et al., considerably smaller than the mean-field estimate. It is interesting to note that with increasing Mn-doping concentration $`N_{\mathrm{Mn}}`$ the mean-field Curie temperature increases \[Eq. (13)\] while collective-fluctuation energy decreases \[Eq. (14)\]. For asymmetric quantum wells the spin-orbit interactions are dominated by Rashba coupling, and, depending on the degree of asymmetry, the collective-excitation energy gap will become larger until, eventually, the system can again reach the mean-field regime. Stronger spin-orbit interactions in the valence band holerefs  may explain the apparent success of mean-field theory in predicting ferromagnetic transition temperatures there. The strong antiferromagnetic interactions that occur between neighbouring Mn ions will reduce holerefs ; sstreviews  the number of Mn spins that can be spontaneously aligned by the relatively weak but longer ranged carrier mediated interactions, but should not preclude ferromagnetism when the Mn mole fraction is low.

In conclusion, we presented a theory of collective spin excitations in DMS quantum wells which explains the anomalously large Knight shift observed in recent experiments. We use our theory to discuss the physics of ferromagnetism in n-type DMS quantum-well systems, which has not yet been observed. We point out that spin-orbit interactions are necessary to support ferromagnetism, predict ferromagnetic transition temperatures in symmetric quantum wells that can be well below those implied by mean-field theory, and suggest that a crossover between collective-fluctuation-limited and mean-field-interaction-limited ferromagnetism can be observed in these systems by using a bias voltage to adjust the spin-orbit interaction strength.

We acknowledge helpful discussions with T. Dietl, D. Frustaglia and M. Governale, and thank F.J. Teran and M. Potemski for showing us their results prior to publication. This work was supported by the Deutsche Forschungsgemeinschaft under the Center for Functional Nanostructures and the Emmy-Noether program, by the Department of Energy under grant DE-FG03-02ER45958 and the Welch Foundation.


text:
# On extreme zeros of classical orthogonal polynomials

## 1. Introduction

Study of extreme zeros of the Hermite, Laguerre and Jacobi polynomials has a long history and most of the classical results are collected in . But only recently attention has been shifted to the case when the parameters may vary with the degree $`k`$ of a polynomial . Most of these results are of the asymptotic nature (with  and  being a remarkable exception) and hold under certain restrictions on the parameters. Recently the author obtained the following explicit uniform bounds  (similar inequalities for the Laguerre case were given earlier in ).

###### Theorem 1

Let $`x_1`$ and $`x_k`$ be the least and the largest zero of the Laguerre polynomial $`L_k^{(\alpha )}(x)`$ respectively, $`\alpha >1.`$ Then

$$x_1>V^2+3V^{4/3}(U^2V^2)^{1/3},$$
(1)
$$x_k<U^23U^{4/3}(U^2V^2)^{1/3}+2,$$
(2)
where $`V=\sqrt{k+\alpha +1}\sqrt{k},`$ $`U=\sqrt{k+\alpha +1}+\sqrt{k}.`$

###### Theorem 2

Let $`x_1`$ and $`x_k`$ be the least and the largest zero of the Jacobi polynomial $`P_k^{(\alpha ,\beta )}(x)`$ respectively, $`\alpha \beta >1.`$ Then

$$x_1>A+3(1A^2)^{2/3}(2R)^{1/3},$$
(3)
$$x_k<B3(1B^2)^{2/3}(2R)^{1/3}+\frac{4q(s+1)}{(r^2+2s+1)^{3/2}},$$
(4)
where

$$s=\alpha +\beta +1,q=\alpha \beta ,r=2k+\alpha +\beta +1,R=\sqrt{(r^2q^2+2s+1)(r^2s^2)},$$

and

$$A=\frac{R+q(s+1)}{r^2+2s+1},B=\frac{Rq(s+1)}{r^2+2s+1}.$$

As the zeros of the Hermite polynomials can be easily expressed through the zeros of the corresponding Laguerre polynomials we will not consider them in this paper.
Previously known results give, roughly speaking, $`V^2<x_1<x_k<U^2,`$  (see also a survey article ) for Laguerre polynomials, and $`A<x_1<x_k<B,`$  for the Jacobi case. It is also known that these bounds are asymptotically correct under certain assumptions on the parameters. On the other hand one can expect that much sharper results similar to these of Theorems 1 and 2 hold in a more general situation. In particular, analogous inequality analogous to (1) - (4) are known for the zeros of Charlier  and binary Krawtchouk polynomials .

The aim of this paper is to show that the bounds given by Theorems 1 and 2 are essentially sharp, thus locating the extreme zeros of the classical orthogonal polynomials with a high precision. Namely we shall establish (in a rather elementary way) two following theorems giving similar inequalities in the opposite direction. Our method is based on so-called Bethe ansatz equations, having some important applications to orthogonal polynomials . It is worth also noticing that the above bounds $`V^2<x_1<x_k<U^2,`$ and $`A<x_1<x_k<B,`$ for the Laguerre and Jacobi polynomials respectively, are an immediate corollary of the Bethe ansatz equation we use here (Lemma 1 below).

###### Theorem 3

Let $`\delta =\frac{1}{k}+\frac{1}{\alpha +1}<\frac{1}{50},`$ then in the notation of Theorem 1,

$$x_1<V^2+\frac{9V^{4/3}}{(U^2V^2)^{1/3}(227\delta ^{2/3})}.$$
(5)
Let $`k30,`$ then

$$x_k>U^2\frac{9U^{4/3}}{2(U^2V^2)^{1/3}}$$
(6)
provided $`\alpha 2(3+2\sqrt{3})k1,`$ and

$$x_k>U^2\frac{9U^{4/3}}{(U^2V^2)^{1/3}(23k^{2/3})},$$
(7)
otherwise.

###### Theorem 4

Let $`\alpha \beta >1,`$ then in the notation of Theorem 2, for $`k5,`$

$$x_1<A+9(1A^2)^{2/3}(2R)^{1/3},$$
(8)
and for $`k56,`$

$$x_k>B9(1B^2)^{2/3}(2R)^{1/3}.$$
(9)
It seems that the bounds in this direction received much less attention. We will use here some rather weak classical inequalities (, sec.6.2).

Theorems 1 \- 4 yield the asymptotics for the extreme zeros given in the next theorem (in the Jacobi case $`x_k`$ and $`B`$ may vanish what leads to more complicate expressions). The meaning of $`O`$-terms here is that for sufficiently large $`k,`$ say $`k>100,`$ one can replace them by absolute constants.

###### Theorem 5

(i) In the notation of Theorem 1, for sufficiently large $`k`$ and $`\alpha >50,`$ the extreme zeros of the Laguerre polynomial $`L_k^{(\alpha )}(x)`$ satisfy

$$\frac{x_1}{V^2}=1+O\left((\alpha +1)^{1/2}(\frac{1}{\alpha +1}+\frac{1}{k})^{1/6}\right),$$
(10)
$$\frac{x_k}{U^2}=1O\left(k^{1/6}(k+\alpha )^{1/2}\right).$$
(11)
(ii) In the notation of Theorem 2, for sufficiently large $`k`$ and $`\alpha \beta >1,`$ the extreme zeros of the Jacobi polynomial $`P_k^{(\alpha ,\beta )}(x)`$ satisfy

$$\frac{x_1}{A}=1+O\left(\left(\frac{(\beta +1)^2}{k(k+\alpha )(k+\beta )}\right)^{2/3}\right);r^2q^2+r^2,$$
(12)
$$\frac{x_1}{A}=1+O\left(\frac{(\beta +1)^{4/3}}{k^{2/3}(k+\beta )^{5/6}\sqrt{k+\alpha }}\right);r^2<q^2+s^2,$$
(13)
Let $`r^2=q^2+s^2+\gamma (s+1)^{2/3}(r^2s^2)^{1/3},`$ then

$$\frac{x_k}{B}=1O\left(\gamma ^1+\gamma ^{2/3}k^{2/9}\right),\gamma >0;$$
(14)
$$\frac{x_k}{B}=1O\left((\alpha k)^{1/3}\right),\gamma <\frac{3(s+1)^{4/3}}{4(r^2s^2)^{1/3}};$$
(15)
$$\frac{x_k}{B}=1O\left(|\gamma |^1+|\gamma |^{1/2}k^{1/3}\right),\frac{3(s+1)^{4/3}}{4(r^2s^2)^{1/3}}\gamma <0;$$
(16)
$$|x_k|=O\left(\frac{1}{k^{1/6}\sqrt{k+\alpha }}\right),|\gamma |1.$$
(17)
It is worth to compare the obtained inequalities with the classical results for the fixed values of the parameters. In particular, in the Laguerre case one has (, Theorem 6.32, see also  for a far-reaching generalization)

$$x_k<\left(\sqrt{4k+2\alpha +2}6^{1/3}(4k+2\alpha +2)^{1/6}i_{11}\right)^2,$$

where $`6^{1/3}i_{11}=1.85575\mathrm{},`$ and $`i_{11}`$ stands for the least positive zero of the Airy function. One can check that for a fixed $`\alpha `$ this differs from (2) only by the better factor $`c=26^{1/3}i_{11},`$ instead of $`3,`$ before the second terms of (2). It is tempting to conjecture that asymptotically for $`k\mathrm{},`$ and uniformly in all the parameters involved, one should get the same constant $`c`$ instead of $`3`$ before the second terms in all the expressions (1)-(4).

The paper is organized as follows. In the next section we establish rather general inequalities being our main tool in the sequel. In sections 3 and 4 we will prove Theorems 3 and 4, dealing with Laguerre and Jacobi polynomials respectively. Section 4 also contains a proof of Theorem 5.

## 2. Bethe Ansatz Inequalities

In this section we will consider real polynomials $`f=f(x)`$ with only real simple zeros $`x_1<x_2<\mathrm{}<x_k,`$ satisfying a differential equation

$$f^{\prime \prime }2af^{}+bf=0$$
(18)
We suppose here that $`a=a(x)`$ and $`b=b(x)`$ are meromorphic functions and none of $`x_i`$ coincides with singularities of $`a`$ or $`b.`$ For such an $`f`$ we define the discriminant $`\mathrm{\Delta }(x)=b(x)a^2(x),`$ and consider the second negative moments of $`f`$ at its zeros

$$S(f,x_i)=\underset{ji}{}\frac{1}{(x_ix_j)^2}.$$

###### Lemma 1

$$S(f,x_i)=\underset{ji}{}\frac{1}{(x_ix_j)^2}=\frac{\mathrm{\Delta }(x_i)2a^{}(x_i)}{3}.$$
(19)
Proof.  Using the logarithmic derivative and (18) we get

$$\frac{1}{(xx_j)^2}=\left(\frac{f^{}}{f}\right)^{}=\frac{f^2ff^{\prime \prime }}{f^2}=\frac{f^22af^{}f+bf^2}{f^2}.$$
(20)
Thus

$$S(f,x_i)=\underset{xx_i}{lim}\left(\frac{f^22af^{}f+bf^2}{f^2}\frac{1}{(xx_i)^2}\right).$$

The result follows on applying four times L’Hôpital’s rule and substituting $`f^{\prime \prime }`$ from (18) at each step. $`\mathrm{}`$

###### Remark 1

Results of this type are called Bethe ansatz equations and are known (or can be routinely established) in a more general situation and weaker smoothness assumptions. We refer to  and the references therein for a more detailed discussion.

###### Lemma 2

$$D(f,x_i,x)=1+(xx_i)^2\left(\frac{\mathrm{\Delta }(x_i)2a^{}(x_i)}{3}\mathrm{\Delta }(x)\right)>0,$$
(21)
provided $`x[x_1,x_k].`$ In particular, if $`a^{}(x_i)0,`$ then

$$32(xx_i)^2\mathrm{\Delta }(x_i)+3(xx_i)^2\left(\mathrm{\Delta }(x_i)\mathrm{\Delta }(x)\right)>0.$$
(22)
Proof.  From (20) we have

$$\frac{1}{(xx_i)^2}+\underset{ji}{}\frac{1}{(xx_j)^2}=(\frac{f^{}(x)}{f(x)}a(x))^2+b(x)a^2(x)\mathrm{\Delta }(x).$$

Since

$$\underset{ji}{}\frac{1}{(xx_j)^2}<\underset{ji}{}\frac{1}{(x_ix_j)^2}=S(f,x_i),$$

for $`x[x_1,x_k],`$ we obtain

$$\frac{1}{(xx_i)^2}+S(f,x_i)>\mathrm{\Delta }(x),$$

and (21),(22) follow by Lemma 1. $`\mathrm{}`$

###### Remark 2

Similar arguments can be apply to $`x[x_1,x_k],`$ say $`x_i<x<x_{i+1},`$ giving an upper bound on $`x_{i+1}x_i.`$ Indeed,

$$\mathrm{\Delta }(x)\frac{1}{(xx_i)^2}+\frac{1}{(xx_{i+1})^2}+\underset{j<i}{}\frac{1}{(xx_j)^2}+\underset{j>i+1}{}\frac{1}{(xx_j)^2}<$$

$$\frac{1}{(xx_i)^2}+\frac{1}{(xx_{i+1})^2}+\underset{j<i}{}\frac{1}{(x_jx_i)^2}+\underset{j>i+1}{}\frac{1}{(x_jx_{i+1})^2}<$$

$$\frac{1}{(xx_i)^2}+\frac{1}{(xx_{i+1})^2}\frac{2}{(x_{i+1}x_i)^2}+S_2(f,x_i)+S_2(f,x_{i+1}).$$

By substituting here $`x=\frac{x_i+x_{i+1}}{2},`$ one obtains

$$(x_{i+1}x_i)^2<\frac{18}{3\mathrm{\Delta }(\frac{x_i+x_{i+1}}{2})\mathrm{\Delta }(x_i)\mathrm{\Delta }(x_{i+1})+2a^{}(x_i)+2a^{}(x_{i+1})},$$

provided the denominator is positive.

We will solve inequality (21) for the Laguerre and Jacobi polynomials in the next section. This will require rather involved calculations but the following simple heuristic arguments show what type of bounds may be expected.
Suppose that $`\mathrm{\Delta }(x)`$ has only two real zeros $`y_1<y_2.`$ Neglecting the term $`2a^{}(x),`$ we obtain that all the zeros of $`f`$ are in the interval $`(y_1,y_2).`$ Let $`x_k`$ be, say, the largest zero of $`f,`$ we put $`x_k=y_2ϵ,`$ and choose $`x=y_2\frac{5ϵ}{9}.`$ Now, on omitting higher derivatives of $`\mathrm{\Delta },`$ that is putting $`\mathrm{\Delta }(y_2\delta )=\mathrm{\Delta }(y_2)\delta \mathrm{\Delta }^{}(y_2)=\delta \mathrm{\Delta }^{}(y_2),`$ (21) can be rewritten as

$$0<1+\frac{16ϵ^2}{81}\left(\frac{\mathrm{\Delta }(x_k)}{3}\mathrm{\Delta }(x)\right)1+\frac{32ϵ^3\mathrm{\Delta }^{}(y_2)}{729}.$$

Thus we obtain $`x_k>y_2+\frac{9}{2}(4\mathrm{\Delta }^{}(y_2))^{1/3}.`$ Notice that similar heuristic considerations given in  yield in the opposite direction $`x_k<y_2+3(4\mathrm{\Delta }^{}(y_2))^{1/3},`$ ( $`\mathrm{\Delta }^{}(y_2)`$ is negative as $`\mathrm{\Delta }(y_2)=0`$).

## 3. Laguerre Polynomials

The Laguerre Polynomials $`L_k^{(\alpha )}(x)`$ are polynomials orthogonal on $`[0,\mathrm{})`$ for $`\alpha >1,`$ with respect to the weight function $`x^\alpha e^x.`$ The corresponding ODE is

$$u^{\prime \prime }(1(\alpha +1)x^1)u^{}+kx^1u=0,u=L_k^{(\alpha )}(x).$$

We also need the explicit representation

$$L_k^{(\alpha )}(x)=\underset{i0}{\overset{k}{}}\left(\genfrac{}{}{0pt}{}{k+\alpha }{ki}\right)\frac{(x)^i}{i!}.$$
(23)
Using the notation of Theorem 1 we get $`k=\frac{(UV)^2}{4},\alpha =VU1,`$ and the condition $`\alpha >1,`$ means $`V>0.`$

We have $`a(x)=\frac{xVU}{2x},a^{}(x)=\frac{VU}{2x^2}>0,`$ and also

$$\mathrm{\Delta }(x)=\frac{(U^2x)(xV^2)}{4x^2},$$
(24)
Let $`x_1`$ and $`x_k`$ be the least and the largest zeros of $`L_k^{(\alpha )}(x)`$ respectively. We need the following (rather weak) bounds on $`x_1`$ and $`x_k,`$ (, sec. 6.2).

$$x_1\frac{(\alpha +1)(\alpha +3)}{2k+\alpha +1}=\frac{2VU(VU+2)}{V^2+U^2}.$$
(25)
By (23) we have $`_{i=0}^kx_i=k(k+\alpha ),`$ implying $`x_1<k+\alpha =\frac{(U+V)^2}{4}<x_k.`$ Moreover, as $`0<S(L_k^{(\alpha )},x_i)<\mathrm{\Delta }(x_i),`$ we get that all the zeros satisfy $`V^2<x_i<U^2,`$ and thus

$$V^2<x_1<\frac{(\alpha +1)(\alpha +3)}{2k+\alpha +1}<x_k<U^2.$$
(26)
###### Lemma 3

For $`V^2<x<x_1,`$

$$\mathrm{\Delta }(x_1)\mathrm{\Delta }(x)<\frac{U^2V^2}{4V^4}(x_1x).$$
(27)
For $`x_k<x<U^2,`$

$$\mathrm{\Delta }(x_k)\mathrm{\Delta }(x)<\frac{U^2V^2}{4x_k^2}(xx_k).$$
(28)
Proof.  Using that $`\frac{(V^2+U^2)xyV^2U^2(x+y)}{xy},`$ is an increasing function in $`x`$ and $`y`$ we obtain

$$\frac{\mathrm{\Delta }(x_1)\mathrm{\Delta }(x)}{x_1x}=\frac{V^2U^2(x+x_1)(V^2+U^2)xx_1}{4x^2x_1^2}<\frac{U^2V^2}{4xx_1}<\frac{U^2V^2}{4V^4};$$

$$\frac{\mathrm{\Delta }(x_k)\mathrm{\Delta }(x)}{xx_k}=\frac{(V^2+U^2)xx_kV^2U^2(x+x_k)}{4x^2x_k^2}<\frac{U^2V^2}{4xx_k}<\frac{U^2V^2}{4x_k^2}.$$

and the result follows. $`\mathrm{}`$

$`Proof`$ $`of`$ $`Theorem`$ 3.
(i) We choose $`x=x_1ϵ,`$ where $`ϵ=\frac{2V^{4/3}}{(U^2V^2)^{1/3}}.`$ Then (22) and (3) gives

$$0<3\frac{ϵ^2(U^2x_1)(x_1V^2)}{2x_1^2}+\frac{3ϵ^3(U^2V^2)}{4V^4}=9\frac{ϵ^2(U^2x_1)(x_1V^2)}{2x_1^2}:=F(x_1).$$

We claim that under our assumptions $`F(x)`$ has two zeros $`y_1<y_2,`$ and $`x_1<y_1.`$ As $`x_1<x_0=\frac{(\alpha +1)(\alpha +3)}{2k+\alpha +1},`$ it is enough to show that $`F(x_0)<0.`$ Putting $`b=\alpha +1,`$ we have

$$F(x_0)=9+\frac{2ϵ^2}{(b+2)^2}+\frac{8ϵ^2k(k+b)}{b^2(b+2)^2}\frac{2ϵ^2k(k+b)}{b^2}.$$

Here

$$\frac{2ϵ^2}{(b+2)^2}+\frac{8ϵ^2k(k+b)}{b^2(b+2)^2}<\frac{8ϵ^2(k+b)^2}{b^4}<\frac{16}{b}\delta ^{1/3}<16\delta ^{4/3},$$

and

$$\frac{2ϵ^2k(k+b)}{b^2}=\left(\frac{16bk(b+k)}{U^4}\right)^{2/3}>\left(\frac{bk}{k+b}\right)^{2/3}=\delta ^{2/3}.$$

Now it is left to check that $`9+16\delta ^{4/3}\delta ^{2/3}<0,`$ for $`\delta <\frac{1}{50},`$ proving the claim. For $`y_1`$ we get

$$y_1=V^2+\frac{9V^2}{h(1+\sqrt{118V^2U^2h^4})9},$$

where $`h=V^{2/3}(U^2V^2)^{1/3},`$ and

$$h(1+\sqrt{118V^2U^2h^4})9>2h\frac{18V^2U^2}{h^3}9<2h\frac{27U^2}{U^2V^2}=$$

$$2h\left(1\frac{27U^{8/3}}{2b^{2/3}(U^2V^2)^{4/3}}\right)<h(227\delta ^{2/3}).$$

As $`22750^{2/3}>0,`$ the result follows.
(ii) We choose $`x=x_kϵ,`$ where $`ϵ=\frac{2U^{4/3}}{(U^2V^2)^{1/3}}.`$ By (22) and (3) we have

$$0<3\frac{ϵ^2(U^2x_k)(x_kV^2)}{2x_k^2}+\frac{3ϵ^3(U^2V^2)}{4x_k^4}=3+\frac{6U^4}{x_k^2}\frac{ϵ^2(U^2x_k)(x_kV^2)}{2x_k^2}<$$

$$\frac{9U^4}{x_k^2}\frac{ϵ^2(U^2x_k)(x_kV^2)}{2x_k^2}.$$

Thus

$$F(x_k):=18x_k^2ϵ^2(U^2x_k)(x_kV^2)>0.$$

The equation $`F(x)=0,`$ has two zeros, $`y_1<y_2,`$ and $`x_k>y_2.`$ Indeed, as $`x_k>x_0=\frac{(V+U)^2}{4},`$ it is enough to check $`F(x_0)<0.`$ We have

$$4F(x_0)=72U^4\frac{ϵ^2(UV)^2(3U^2+10VU+3V^2)}{4}72U^43U^{8/3}(U^2V^2)^{4/3}$$

$$3U^{8/3}(U+V)^{4/3}(24(UV)^{4/3})=3U^{8/3}(U+V)^{4/3}(24(4k)^{2/3})<0,$$

for $`k30.`$ Thus,

$$x_k>y_2=U^2\frac{9U^2}{9+U^{2/3}(U^2V^2)^{1/3}\left(1+\sqrt{118V^2U^{2/3}(V^2+U^2)^{4/3}}\right)}>$$

$$U^2\frac{9U^2}{2U^{2/3}(U^2V^2)^{1/3}+9\frac{18V^2}{U^2V^2}}.$$

Finally, $`9\frac{18V^2}{U^2V^2}0,`$ if $`\alpha 2(3+2\sqrt{3})k1,`$ proving (6). Otherwise,

$$2U^{2/3}(U^2V^2)^{1/3}\frac{18V^2}{U^2V^2}=U^{2/3}(U^2V^2)^{1/3}\left(2\frac{18V^2}{U^{2/3}(U^2V^2)^{4/3}}\right)>$$

$$U^{2/3}(U^2V^2)^{1/3}\left(2\frac{9b^2}{2^{5/3}k^{2/3}(k+b)^2}\right)>U^{2/3}(U^2V^2)^{1/3}(23k^{2/3}),$$

and (7) follows. $`\mathrm{}`$

## 4. Jacobi Polynomials

The Jacobi Polynomials $`P_k^{(\alpha ,\beta )}(x)`$ are polynomials orthogonal on $`[1,1]`$ for $`\alpha ,\beta >1,`$ with respect to the weight function $`(1x)^\alpha (1+x)^\beta .`$ The corresponding ODE is

$$u^{\prime \prime }\frac{(\alpha +\beta +2)x+\alpha \beta }{1x^2}u^{}+\frac{k(k+\alpha +\beta +1)}{1x^2}u=0,u=P_k^{(\alpha ,\beta )}(x).$$

We will use the notation of Theorem 2 and put $`p=r^2+2s+1`$ throughout this section. We have

$$\mathrm{\Delta }(x)=\frac{px^2+2q(s+1)x+s^2+q^2r^2}{4(1x^2)^2}=\frac{p(xA)(Bx)}{4(1x^2)^2}.$$
(29)
As

$$a^{}(x)=\frac{((\alpha +\beta +2)x+\alpha \beta )^2+4(\alpha +1)(\beta +1)}{2(\alpha +\beta +2)(1x^2)^2}>0,$$

we can use (22) and moreover, as $`\mathrm{\Delta }(x_i)>0,`$ we obtain

$$A<x_i<B$$
(30)
In the opposite direction it is known (, sec. 6.2)

$$x_1<\frac{2k+\alpha \beta 2}{2k+\alpha +\beta }<\frac{2k+\beta \alpha 2}{2k+\alpha +\beta }<x_k.$$
(31)
It is easy to show (see e.g. ) that $`x_1<0,`$ for $`\alpha \beta .`$

###### Lemma 4

For $`A<x<x_1,`$

$$\mathrm{\Delta }(x_1)\mathrm{\Delta }(x)<\frac{R(x_1x)}{2(1x_1^2)^2}.$$

For $`x_k<x<B,`$

$$\mathrm{\Delta }(x_k)\mathrm{\Delta }(x)<\frac{R(xx_k)}{2(1x_k^2)^2}.$$

Proof.  We have

$$\mathrm{\Delta }(x_1)\mathrm{\Delta }(x)<\frac{p}{4(1x_1^2)^2}\left((x_1A)(Bx_1)(xA)(Bx)\right)<$$

$$\frac{p}{4(1x_1^2)^2}(BA)(x_1x)=\frac{R(x_1x)}{2(1x_1^2)^2}.$$

$$\mathrm{\Delta }(x_k)\mathrm{\Delta }(x)<\frac{p}{4(1x_k^2)^2}\left((x_kA)(Bx_k)(xA)(Bx)\right)<$$

$$\frac{p}{4(1x_k^2)^2}(BA)(xx_k)=\frac{R(xx_k)}{2(1x_k^2)^2}.$$

$`\mathrm{}`$

$`Proof`$ $`of`$ $`Theorem`$ 4.
(i) Choose $`ϵ=\frac{(22A^2)^{2/3}}{R^{1/3}},`$ and put $`x=x_1ϵ.`$ Then $`x>A,`$ otherwise there is nothing to prove. Using the previous lemma and (22) we obtain

$$0<3\frac{ϵ^2p(x_1A)(Bx_1)}{2(1x_1^2)^2}+\frac{3ϵ^3R}{2(1x_1^2)^2}<$$

$$\frac{3(1A^2)^2}{(1x_1^2)^2}\frac{ϵ^2p(x_1A)(Bx_1)}{2(1x_1^2)^2}+\frac{3ϵ^3R}{2(1x_1^2)^2}.$$

Thus, we get

$$18(1A^2)^2ϵ^2p(Bx_1)(x_1A):=F(x_1)>0.$$
(32)
We shall show that this quadratic has two real zeros $`z_1<z_2,`$ and $`x_1<z_1.`$ For, it is enough to prove $`F(\frac{A+B}{2})<0,`$ and $`\frac{A+B}{2}.`$ The last claim follows from (31), as

$$x_1<\frac{2k+\alpha \beta 2}{2k+\alpha +\beta }<\frac{A+B}{2}.$$

Indeed, $`\alpha ,\beta >1,`$ and we obtain

$$\frac{A+B}{2}x_1>\frac{A+B}{2}+\frac{2k+\alpha \beta 2}{2k+\alpha +\beta }=$$

$$\frac{4(2k^3+(3\alpha +\beta +4)k^2+(\alpha ^2+\alpha \beta +4\alpha +4\beta +4)k+(\alpha +1)(\alpha +\beta +2))}{(r1)p}>0.$$

Now we have

$$F(\frac{A+B}{2})=72(1A^2)^2ϵ^2p(BA)^2,$$

and it is negative whenever

$$2R^4>729p^3(1A^2)^2.$$
(33)
As

$$\frac{d}{dq}\left(\frac{1A^2}{R^2}\right)=\frac{2((s+1)^2q^2)}{(pq^2)R(qR+(s+1)(pq^2))}<0,$$

and for $`q=0,`$

$$\frac{1A^2}{R^2}=\frac{(s+1)^2}{p^2(r^2s^2)}.$$

We have

$$\frac{p^3(1A^2)^2}{R^4}<\frac{(s+1)^4}{p(r^2s^2)^2}<\frac{(s+1)^4}{16k^2(k+s)^2(2k+s)^2}<\frac{1}{16k^2}<\frac{2}{729},$$

provided $`k5.`$ This proves (33) and, thus, $`x_1<z_1.`$
Finally, solving $`F(x)=0,`$ we obtain

$$x_1<A+\frac{18(1A^2)^2}{ϵ^2R\left(1+\sqrt{1\frac{18p(1A^2)^2}{ϵ^2R^2}}\right)}<A+\frac{18(1A^2)^2}{ϵ^2R}=A+\frac{9(1A^2)^{2/3}}{(2R)^{1/3}}.$$

(ii) Choose $`ϵ=\frac{(22B^2)^{2/3}}{R^{1/3}},`$ and put $`x=x_k+ϵ.`$ Similarly to the previous case we get

$$18(1B^2)^2ϵ^2p(Bx_k)(x_kA):=F(x_k)>0.$$
(34)
We shall show that $`x_k`$ is greater than the largest zero of $`F(x)=0.`$ To prove this we establish $`F(x_0)<0,`$ where $`x_0=\frac{2k+\beta \alpha 2}{2k+\alpha +\beta }<x_k,`$ by Lemma 31. For it is enough to show

$$G=\left(\frac{18(1B^2)^2}{ϵ^2p(Bx_k)(x_kA)}\right)^3=\frac{729R^2(1B^2)^2}{2(p(Bx_k)(x_kA))^3}<1.$$

We have

$$\frac{d}{dq}\left(\frac{1B^2}{R^2}\right)=\frac{2(r^2s^2)(q(r^2s^2)+(s+1)R)}{pR^4}>0.$$

As $`q=\alpha \beta <\alpha +\beta +2=s+1,`$ we obtain

$$\frac{1B^2}{R^2}<\frac{4(s+1)^2}{p^2(r^2s^2)},$$

that is

$$1B^2<\frac{4(s+1)^2R^2}{p^2(r^2s^2)}.$$

We also have

$$p(Bx_0)(x_0A)=\frac{16(\alpha +1)((k1)(\alpha +1)+k(k+\beta )(2k+\alpha +\beta ))}{(2k+\alpha +\beta )^2}$$

$$\frac{16k(\alpha +1)(k+\beta )}{2k+\alpha +\beta +1}>\frac{2(rq1)(s+q+1)(rs)}{r}.$$

Therefore we obtain

$$G<\frac{729r^3(s+1)^4(r+s)(pq^2)^3}{p^4(rq1)^3(rs)^2(s+q+1)^3}<\frac{729(s+1)(r+s)(pq^2)^3}{r^5(rq1)^3(rs)^2}.$$

The last expression is an increasing function in $`q,`$ and substituting $`q=s+1,`$ we get

$$G<\frac{729(s+1)(rs)(r+s)^4}{r^5(rs2)^3}=\frac{2916k(s+1)(k+s)^4}{(k1)^3(2k+s)^5}<\frac{2916k}{(k1)^3}<1,$$

for $`k56.`$ Finally, solving $`F(x_0)=0,`$ we obtain

$$x_k>B\frac{18(1B^2)^2}{ϵ^2R\left(1+\sqrt{1\frac{18p(1B^2)^2}{ϵ^2R^2}}\right)}>B\frac{18(1B^2)^2}{ϵ^2R}=B\frac{9(1B^2)^{2/3}}{(2R)^{1/3}}.$$

$`\mathrm{}`$

###### Remark 3

More accurate calculations show that in fact (9) holds for $`k20,`$ instead of $`56.`$ It is also easy to improve the constant $`9`$ in (8), (9) to $`\frac{9}{2o(1)},`$ similarly to the Laguerre case.

$`Proof`$ $`of`$ $`Theorem`$ 5. The asymptotics for the Laguerre case is an easy exercise, here we will establish (12)-(17).
Notice that the inequality $`r^2s^2+q^2`$ is equivalent to $`Rq(s+1).`$ We also observe that the last term in (4) may be ignored. Indeed,

$$1B^2=\frac{(q+s+1)^2(R+pq(s+1))}{p(R+p+q(s+1))},$$

and this is an increasing function on $`R.`$ As $`q<s+1,`$ we get $`R>r^2s^2,`$ what implies

$$1B^2>\frac{(q+s+1)^2(2r^2s^2q(s+1))}{p(2r^2s^2+q(s+1))}>\frac{2(\alpha +1)^2(2r^2s^2q(s+1))}{p^2}>$$

$$\frac{(\alpha +1)^2(k+\alpha )(k+\beta )}{p^2}.$$

Now calculations yield

$$\left(\frac{q(s+1)R^{1/3}}{p^{3/2}(1B^2)^{2/3}}\right)^6<c\frac{q^6(s+1)^6k}{(\alpha +1)^8(k+\alpha )^4(k+\beta )^3},$$

for some positive constant $`c`$. This expression is a decreasing function in $`\beta `$ and for $`\beta =1,`$ is $`O\left(\frac{(\alpha +1)^4}{k^2(k+\alpha )^4}\right).`$ Thus the last term in (4) is negligible whenever $`k\mathrm{}.`$
$`Proof`$ $`of`$ (12). As $`Rq(s+1),`$ we have $`|A|>\frac{R}{2r^2},`$ and

$$1A^2<2(1+A)=\frac{2(s+1q)^2}{R+pq(s+1)}<\frac{8(\beta +1)^2}{p}<\frac{16(\beta +1)^2}{r^2}.$$

Therefore,

$$\left(\frac{(1A^2)^{2/3}}{|A|R^{1/3}}\right)^{3/2}<\frac{128(\beta +1)^2r}{R^2}<\frac{256(\beta +1)^2r}{(r^2q^2)(r^2s^2)}<\frac{32(\beta +1)^2}{k(k+\alpha )(k+\beta )},$$

and (12) follows.
$`Proof`$ $`of`$ (13). As $`R<q(s+1),`$ we get $`q^2>r^2s^2.`$ This yields $`1<\beta <2k+\alpha 2\sqrt{k(2k+2\alpha +1)},`$ $`\alpha >2k1+2\sqrt{k(2k1)},`$ and $`k<\alpha /2.`$ Thus, $`s`$ is a large positive number and $`|A|>\frac{qs}{r^2}.`$ Now, using $`R>r^2s^2,`$ we obtain

$$1A^2<\frac{2(s+1q)^2}{R+pq(s+1)}<\frac{8(\beta +1)^2}{2r^2s^2q(s+1)}<\frac{4(\beta +1)^2}{\alpha (k+\beta )}.$$

This yields

$$\left(\frac{(1A^2)^{2/3}}{|A|R^{1/3}}\right)^6<\frac{256r^{12}(\beta +1)^8}{\alpha ^4(k+\beta )^4q^6s^6(r^2q^2)(r^2s^2)}<\frac{10^8(\beta +1)^8}{k^4(k+\alpha )^3(k+\beta )^5},$$

and the result follows.
$`Proof`$ $`of`$ (14). The condition $`q^2=r^2s^2\gamma (s+1)^{2/3}(r^2s^2)^{1/3},\gamma >0,`$ implies that $`R>q(s+1),`$ and $`B>0.`$ Rewriting $`B`$ as $`\frac{r^2s^2q^2}{R+q(s+1)}`$ we obtain $`B>\frac{r^2s^2q^2}{2R}.`$ We also have

$$1B^2<2(1B)=\frac{8(\alpha +1)^2}{R+p+q(s+1)}<\frac{8(\alpha +1)^2}{r^2}.$$

Hence

$$\left(\frac{(1B^2)^{2/3}}{BR^{1/3}}\right)^3=\frac{512(\alpha +1)^4}{r^4}\left(\gamma ^3+\frac{(r^2s^2)^{1/3}}{\gamma ^2(s+1)^{4/3}}\right)<$$

$$512\gamma ^3+\frac{900k^{1/3}(\alpha +1)^4(k+\alpha +\beta +1)^{1/3}}{\gamma ^2(\alpha +\beta +2)^{4/3}(2k+\alpha +\beta +1)^4}.$$

The second term here is a decreasing function in $`\beta >1,`$ and does not exceed

$$\frac{900k^{1/3}(\alpha +1)^{8/3}(k+\alpha )^{1/3}}{\gamma ^2(2k+\alpha )^4}<900k^{2/3}\gamma ^2,$$

and the result follows.
$`Proof`$ $`of`$ (15), (16). In those case $`k<\frac{\sqrt{2\alpha ^2+2\alpha +1}\alpha }{2},`$ $`\alpha >2k1+2\sqrt{k(2k1)},`$ and so $`\alpha `$ is large. Therefore,

$$0>r^2q^2s^2>r^22(s+1)^2>r^24s^2,$$

and hence $`s<r<2s.`$ By $`q<s+1,`$ it follows

$$\gamma <\frac{2s^2+2s+1r^2}{(s+1)^{2/3}(r^2s^2)^{1/3}}<\frac{(s+1)^{4/3}}{(r^2s^2)^{1/3}}.$$

Rewriting $`B`$ as $`\frac{r^2q^2s^2}{R+q(s+1)},`$ and using $`R<q(s+1),`$ we have $`B>\frac{r^2q^2s^2}{2q(s+1)}.`$ Now,

$$\left(\frac{(1B^2)^{2/3}}{BR^{1/3}}\right)^6B^6R^2<\frac{64q^2(s+1)^6}{R^2(r^2q^2s^2)^6}=$$

$$\frac{64(s+1)^{4/3}\left((r^2s^2)^{2/3}\gamma (s+1)^{2/3}\right)^3}{\gamma ^6(r^2s^2)^2\left((s+1)^{4/3}+\gamma (r^2s^2)^{1/3}\right)}<$$

$$\frac{256(s+1)^{10/3}}{\gamma ^3(r^2s^2)^2\left((s+1)^{4/3}+\gamma (r^2s^2)^{1/3}\right)}+\frac{256(s+1)^{4/3}}{\gamma ^6\left((s+1)^{4/3}+\gamma (r^2s^2)^{1/3}\right)}=I_1+I_2.$$

Now we shall consider two cases corresponding to the restrictions in (16) and (15). If $`\gamma \frac{3(s+1)^{4/3}}{4(r^2s^2)^{1/3}},`$ that is $`q^2<r^2\frac{s^26s3}{4},`$ then

$$I_1<\frac{1024(s+1)^2}{\gamma ^3(r^2s^2)^2}<\frac{64}{\gamma ^3k^2}.$$

Otherwise, using $`\gamma <\frac{2s^2+2s+1r^2}{(s+1)^{2/3}(r^2s^2)^{1/3}},`$ and $`k<\frac{\alpha }{4},`$ for large $`\alpha ,`$ we get

$$I_1<\frac{256(s+1)^6}{(r^2s^2)^2(2s^2+2s+1r^2)}<\frac{128(s+1)^6}{\alpha ^6k^2(k+s)^2}=O\left(\alpha ^2k^2\right).$$

Similarly, $`I_2=O\left(\gamma ^6\right),`$ if $`\gamma \frac{3(s+1)^{4/3}}{4(r^2s^2)^{1/3}},`$

$$I_2=O\left(\frac{(r^2s^2)^2}{(s+1)^8}\right)=O\left(k^2a^6\right),$$

if $`\frac{3}{4}<\frac{\gamma (r^2s^2)^{1/3}}{(s+1)^{4/3}}\frac{6}{7},`$ and $`I_2=\left(k\alpha ^5\right),`$ otherwise. These readily yield (15), (16).
$`Proof`$ $`of`$ (17) In this case

$$|B|=|\frac{r^2q^2s^2}{R+q(s+1)}|<\frac{|\gamma |(s+1)^{1/3}}{(r^2s^2)^{1/6}\sqrt{(s+1)^{4/3}+\gamma (r^2s^2)^{1/3}}}=$$

$$O\left(\frac{|\gamma |}{(s+1)^{1/3}(r^2s^2)^{1/6}}\right)=O\left(\frac{1}{k^{1/6}\sqrt{k+\alpha }}\right),$$

and $`R=s\sqrt{r^2s^2}\left(1+o(1)\right).`$ Thus,

$$(1B^2)^{2/3}R^{1/3}=\frac{\left(1o(1)\right)}{s^{1/3}(r^2s^2)^{1/6}}=O\left(\frac{1}{k^{1/6}\sqrt{k+\alpha }}\right),$$

and (17) follows. $`\mathrm{}`$


text:
# Influence of moving breathers on vacancies migration

## 1 Introduction

The interaction of moving localized excitations with defects is presently a subject of great interest and can be connected with certain phenomena observed in crystals and biomolecules. Recently, Sen *et al*  have observed that, when a silicon crystal is irradiated with an ion beam, the defects are pushed towards the edges of the sample. The authors suggest that mobile localized excitations called quodons, which are created in atomic collisions, are responsible for this phenomenon. The interpretation is that the quodons are moving discrete breathers that can appear in 2D and 3D lattices and move following a quasi-one-dimensional path . The interaction of moving breathers with defects is currently of much interest  (for a review on the concept of discrete breather see, e.g. ).

In this paper, we consider a simple one-dimensional model in order to study how a moving breather can cause a lattice defect to move. This study is new in the sense that most studies that consider the interaction of moving discrete breathers with defects, assume that the position of the latter are fixed and cannot move through the lattice. The defect that we consider is a lattice vacancy, which is represented by an empty well or anti-kink in a Frenkel–Kontorova model . The aim of this paper is to determine in which conditions the vacancy moves towards the ends of the chain. This is a previous step to reproduce the phenomenon observed in  for higher dimensional lattices.

We have observed, as it will be explained in detail in Section 3, that different phenomena can occur: the vacancy can move forwards or backwards or remain at rest, and the breather can be reflected, refracted or trapped. This is quite a different scenario from the continuous case, in which the vacancy (or anti-kink) only moves backwards and the breather is always refracted .

## 2 The model

In order to study the migration of vacancies, we consider a Hamiltonian Frenkel–Kontorova model with anharmonic interaction potential :

$$H=\underset{n}{}\frac{1}{2}\dot{x}_n^2+V(x_n)+CW(x_{n+1}x_n).$$
(1)
The dynamical equations are:

$$F(\{x_n\})\ddot{x}_n+V^{}(x_n)+C[W^{}(x_nx_{n1})W^{}(n_{n+1}x_n)]=0,$$
(2)
where $`\{x_n\}`$ are the absolute coordinates of the particles; $`V(x)`$ is the on–site potential, which is chosen of the sine-Gordon type:

$$V(x)=\frac{L^2}{4\pi ^2}\left(1\mathrm{cos}\frac{2\pi x}{L}\right),$$
(3)
with $`L`$ being the period of the lattice. The choice of a periodic potential allows us to represent a vacancy easily. Thus, if we denote the vacancy site as $`n_\mathrm{v}`$ (see figure 1), the displacements of the particles with respect to their equilibrium position are:

$$\{\begin{array}{cc}u_n=x_nnL\hfill & n<n_\mathrm{v}\hfill \\ & \\ u_n=x_n(n+1)L\hfill & n>n_\mathrm{v}.\hfill \end{array}$$
(4)
The interaction potential $`W(x)`$ is of the Morse type:

$$W(x)=\frac{1}{2b^2}[\mathrm{exp}(b(xa))1]^2,$$
(5)
where $`a`$ is the distance between neighboring minima of the interaction potential. In order to avoid discommensurations, we have chosen $`L=a=1`$. The parameter $`b`$ is a measure of the inverse of the width of the interaction potential well. The interaction between particles is stronger when $`b`$ decreases. The $`1/b^2`$ factor allows a Taylor expansion of $`W(x)`$ at $`x=a`$ independent on $`b`$ up to second order and, in consequence, the curvature at the bottom of the interaction potential $`CW(x)`$ depends only on $`C`$.

The reason for the choice of this potential is twofold. On the one hand, it represents a way of modelling the interaction between atoms in a lattice so as the lager the distance between particles, the weaker the interaction between them becomes. On the other hand, if a harmonic interaction potential were chosen, apart from being unphysical in this model, the movement of the breather would involve a great amount of phonon radiation, making it impossible to perform the study developed in this paper.

Throughout this paper, the results correspond to a breather frequency $`\omega _\mathrm{b}=0.9`$. Values of $`\omega _\mathrm{b}[0.9,1)`$ lead to qualitatively similar results. Values of $`\omega _\mathrm{b}0.9`$ have not been chosen as moving breathers do not exist .

## 3 Numerical Results

### 3.1 Preliminaries

In order to investigate the migration of vacancies in our model, we have launched a moving breather towards the vacancy located at the site $`n_\mathrm{v}`$. This moving breather has been generated using a simplified form of the marginal mode method , which consists basically in adding to the velocity of a stationary breather a perturbation which breaks its translational symmetry, and letting it evolve in time. In these simulations, a damping term for the particles at the edges has been introduced in order that the effects of the phonon radiation be minimized.

The initial perturbation, $`\{\stackrel{}{V}_n\}`$ has been chosen as $`\stackrel{}{V}=\lambda (\mathrm{},0,1/\sqrt{2},0,1/\sqrt{2},0,\mathrm{})`$, where the nonzero values correspond to the neighboring sites of the initial center of the breather. This choice of the perturbation allows it to be independent on the parameters of the system $`b`$ or $`C`$. If the pinning mode were chosen as an initial perturbation, it would depend on the parameters of the system.

### 3.2 Breather–vacancy interaction

When a moving breather reaches the site occupied by the particle adjacent to the vacancy, i.e., the location $`n_\mathrm{v}1`$, it can jump to the vacancy site or remain at rest. If the former takes place, the vacancy moves backwards. However, if the interaction potential is wide enough, the particle at the $`n_\mathrm{v}+1`$ site, can feel the effect of the moving breather at the $`n_\mathrm{v}1`$ site and it can also move towards the vacancy site. In the last case, the vacancy moves forwards. Figures 2 and 3 illustrate both phenomena.

It is interesting that the vacancy can migrate along several sites before stopping if the interaction between particles is strong enough (see Figure 4). The largest jumps we have detected are of eleven sites.

There is no apparent correlation between the characteristics of the moving breather, e.g. its kinetic energy and its phase (which has no obvious definition but depends on the initial distance between the breather and the vacancy and the initial velocity of the breather). As an example, Figure 5 shows the vacancy jumps corresponding to different values of the translational kinetic energy of the breather. We have not been able to detect any pattern. The same plot with respect to the breather distance to the vacancy has a similar appeareance.

Numerical simulations show that the occurrence of the three different cases depends highly on the relative phase of the incoming breather and the particles adjacent to the vacancy. However, some conclusions can be extracted: 1) The incident breather always losses energy; 2) The breather can be reflected, trapped (with emission of energy) or refracted by the vacancy, in analogy to the interaction moving breather-mass defect ; 3) the refraction of the breather (i.e. the breather passes through the vacancy) can only take place if the vacancy moves backwards, i.e. the particle to the left jumps one site in the direction of the breather. A explanation of this fact is that the particles to the right of the vacancy, in order to support a moving breather, need a strong interaction which cannot be provided by the interaction across a vacancy site, because the distance correspond to the soft part of the Morse potential.

### 3.3 Numerical simulations

As mentioned earlier, the moving breather–vacancy interaction is highly phase-dependent in a non obvious way. That is, the interaction depends on the velocity of the breather and the distance between the breather and the vacancy. Consequently, a systematic study of the state of the moving breather and the vacancy after the interaction cannot be performed.

Therefore, we have performed a great number of simulations each one consisting in launching a single breather towards the vacancy site. In particular, we have chosen 1000 breathers following a Gaussian distribution of the perturbation parameter $`\lambda `$ with mean value $`0.13`$ and variance $`0.03`$ for different values of the parameters $`b`$ and $`C`$. Figure 6 shows the probabilities that the vacancy remains at its original site, or that it jumps backwards or forwards, for $`C=0.5`$ and $`C=0.4`$. Figure 7 shows the mean values of the number of vacancy jumps for the forward and backward movement as a function of the inverse potential width $`b`$.

An important consequence can be extracted from this figure. There are two different regions of values for the parameter $`b`$, separated by a critical value $`b_0(C)`$. For $`b>b_0(C)`$, the probability that the vacancy moves forwards is almost zero, whereas for $`b<b_0(C)`$, this probability is significant. For example, $`b_0(C=0.5)0.70`$ and $`b_0(C=0.4)0.55`$.

Figure 8 represents this dependence for an uniform distribution of $`\lambda (0.10,0.16)`$, and shows the occurrence of the same phenomenon. Thus, this result seems to be independent on how $`\lambda `$ is distributed.

### 3.4 Analysis of some results. Vacancy breather bifurcation.

The non-existence of forwards vacancy migration can be explained through a bifurcation. If we analyze the spectrum of the Jacobian of the dynamical equations (2) defined by $`𝒥_xF(\{x_n\})`$, bifurcations can be detected. A necessary condition for the occurrence of a bifurcation is that an eigenvalue of $`𝒥`$ becomes zero. Figure 9 shows the dependence of the eigenvalues closest to zero with respect to $`b`$ for $`C=0.5`$ and $`C=0.4`$. It can be observed that, in both cases, there is an eigenvalue that crosses zero in $`b(0.65,0.70)`$ for $`C=0.5`$ and in $`b(0.50,0.55)`$ for $`C=0.4`$. These values agree with the points where the probability of the jump forward vanishes.

These bifurcations are related to the disappearance of the entities we call *vacancy breathers*. They are defined as breathers centered at the site neighboring to the vacancy, e.g. the $`n_\mathrm{v}1`$ or $`n_\mathrm{v}+1`$ sites. It can be observed (figure 10) that, for $`b`$ below the bifurcation value, vacancy breathers do not exist.

## 4 Conclusions

In this paper, we have observed that a moving breather can force a vacancy defect to move forwards, backwards or let it at its site. We have also analyzed the influence of the width of the coupling potential and the coupling strength on the possibility of movement of a vacancy after the collision with a moving breather. We have observed that the width of the potential must be higher than a threshold value in order that the vacancy can move forwards. This behaviour is relevant because experiments developed in crystals show that the defects are pushed towards the edges. We have also established that the non–existence of a breather centered at the sites adjacent to the vacancy is a necessary condition for the forward vacancy movement.

The incident breathers can be trapped, in the sense that the energy becomes localized at the vacancy next–neighbors, which radiate and eventually the energy spreads through the lattice. It can also be transmitted or reflected. The transmission can only occur if the vacancy moves backwards. The moving breather always losses energy but there is not a clear correlation between the vacancy and breather behaviours.

The authors acknowledge Prof. F Palmero, from the GFNL of the University of Sevilla, for valuable suggestions. They also acknowledge partial support under the European Commission RTN project LOCNET, HPRN-CT-1999-00163. J Cuevas also acknowledges an FPDI grant from ‘La Junta de Andalucía’.