format: text

text:
# LIMIT THEOREMS AND ABSORPTION PROBLEMS FOR QUANTUM RANDOM WALKS IN ONE DIMENSION

## 1 Introduction

The classical random walk (CRW) in one dimension is the motion of a particle located on the set of integers. The particle moves at each step either one unit to the left with probability $`p`$ or one unit to the right with probability $`q=1p.`$ The directions of different steps are independent of each other. This CRW is often called the Bernoulli random walk. In the present paper, we consider quantum variations of the Bernoulli random walk and refer to such processes as quantum random walks (QRWs) here. Sometimes the QRW is also called the quantum walk.

Very recently, considerable work has been done on discrete-time and continuous-time QRWs by a number of groups in connection with quantum computing. Examples include Aharonov et al. (2001), Ambainis et al. (2001), Bach et al. (2002), Brun, Carteret and Ambainis (2002a, 2002b, 2002c), Childs, Farhi and Gutmann (2002), Childs et al. (2002), Dür et al. (2002), Kempe (2002), Kendon and Tregenna (2002a, 2002b), Konno (2002a, 2002b), Konno, Namiki and Soshi (2002), Konno, Namiki, Soshi and Sudbury (2003), Leroux (2002), Mackay et al. (2002), Moore and Russell (2001), Severini (2002a, 2002b), Travaglione and Milburn (2002), Yamasaki, Kobayashi and Imai (2002). For a more general setting including quantum cellular automata, see Meyer (1996). The present paper is concerned only with the discrete-time case.

QRWs behave quite differently from CRWs. It is well known that the probability distribution of a CRW is given by the form of a binomial distribution. The variance of the CRW increases linearly with the number of time steps. By contrast, the probability distribution of a QRW has a complicated oscillatry form and the variance of the QRW increases quadratically with the number of time steps (see Ambainis et al. (2001), Konno (2002a), for examples). Such behavior of the QRW is due to the interference between the separate paths of the walk.

The QRW was introduced by Gudder (see Section 7.4 in his book (1988)) to describe the motion of a quantum object in discrete space-time (very recently, I happened to know this fact from Severini (2002a)). QRWs have been re-discovered and studied intensively by Ambainis et al. (2001) in the context of quantum computation recently. So for simplicity we call the former QRW the G-type QRW, and the latter QRW the A-type QRW respectively in this paper. Both QRWs are essentially same, however they have a difference. One of the main purpose of this paper is to clarify the different behaviour between the two types QRWs, for examples, limit distribution, symmetry of distribution, absorption probability.

Here we present not only our recent results on limit theorems, symmetry of distribution, and absorption probabilities for A-type QRWs, but also the related results for G-type QRWs, by using our PQRS method based on combinatorics which is different from Fourier analysis. This paper is organized as follows. Section 2 treats the definitions of both A-type and G-type QRWs and explains the PQRS method. In Section 3, we give a new type of limit theorems for both types of QRWs. Section 4 is devoted to absorption problems for QRWs.

## 2 Definition and PQRS Method

The time evolution of the one-dimensional QRW is given by the following unitary matrix:

$`U=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]`$
where $`a,b,c,d𝐂`$ and $`𝐂`$ is the set of complex numbers. So we have

$`|a|^2+|b|^2=|c|^2+|d|^2=1,a\overline{c}+b\overline{d}=0,`$
$`c=\mathrm{}\overline{b},d=\mathrm{}\overline{a}`$
where $`\overline{z}`$ is a complex conjugate of $`z𝐂`$ and $`\mathrm{}=detU=adbc.`$ We should note that the unitarity of $`U`$ gives $`|\mathrm{}|=1.`$

For the A-type QRW, each coin performs the evolution

$`|LU|L=a|L+c|R,`$
$`|RU|R=b|L+d|R`$
at each time step for which that coin is active, where R and L can be respectively thought of as the heads and tails states of the coin, or equivalently as an internal chirality state of the particle. The value of the coin controls the direction in which the particle moves. When the coin shows L, the particle moves one unit to the left, when it shows R, it moves one unit to the right. (The G-type QRW is also interpreted in a similar way. We will give a precise definition of both types of QRWs later.) In this meaning, the QRW can be considered as a quantum version of the CRW with an additional degree of freedom called the chirality which takes values left and right. As for ways to regain the CRW from the QRW, see Brun, Carteret and Ambainis (2002a, 2002b, 2002c).

The amplitude of the location of the particle is defined by a 2-vector $`𝐂^2`$ at each location at any time $`n`$. The probability the particle is at location $`k`$ is given by the square of the modulus of the vector at $`k`$. For the $`j`$-type QRW ($`j=A,G`$), let $`|\mathrm{\Psi }_{j,k}(n)`$ denote the amplitude at time $`n`$ at location $`k`$ where

$`|\mathrm{\Psi }_{j,k}(n)=\left[\begin{array}{cc}\psi _{j,k}^L(n)& \\ \psi _{j,k}^R(n)& \end{array}\right]`$
with the chirality being left (upper component) or right (lower component). Then the dynamics for $`|\mathrm{\Psi }_{A,k}(n)`$ for the A-type QRW is given by the following transformation:

(2.3) $`|\mathrm{\Psi }_{A,k}(n+1)=P_A|\mathrm{\Psi }_{A,k+1}(n)+Q_A|\mathrm{\Psi }_{A,k1}(n)`$
where

$`P_A=\left[\begin{array}{cc}a& b\\ 0& 0\end{array}\right],Q_A=\left[\begin{array}{cc}0& 0\\ c& d\end{array}\right]`$
It is noted that $`U=P_A+Q_A.`$ The unitarity of $`U`$ ensures that the amplitude always defines a probability distribution for the location.

On the other hand, the G-type of QRW can be determined by

$`|\mathrm{\Psi }_{G,k}(n+1)=P_G|\mathrm{\Psi }_{G,k+1}(n)+Q_G|\mathrm{\Psi }_{G,k1}(n)`$
where

$`P_G=\left[\begin{array}{cc}a& 0\\ c& 0\end{array}\right],Q_G=\left[\begin{array}{cc}0& b\\ 0& d\end{array}\right]`$
In this G-type case also, we see that $`U=P_G+Q_G`$, and the unitarity of $`U`$ ensures that the amplitude always defines a probability distribution for the location. In our one-dimensional setting, Gudder considered the following simple model (see Eqs. (7.33) and (7.34) in page 279 of Gudder (1988)) and computed $`|\mathrm{\Psi }_{G,k}(n)`$ (in our notation) for this model (see Corollary 7.24 in page 285 of his book) by using a Fourier analysis which is different from our PQRS method:

$`U=\left[\begin{array}{cc}a& ib\\ ib& a\end{array}\right]`$
where $`0<a<1`$ and $`b=\sqrt{1a^2}`$.

The simplest and well-studied example of an A-type QRW is the Hadamard walk whose unitary matrix $`U`$ is defined by

$`H={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cc}1& 1\\ 1& 1\end{array}\right]`$
The dynamics of this walk corresponds to that of the symmetric CRW. In general, the following unitary matrices can also lead to symmetric QRWs:

$`U_{\eta ,\varphi ,\psi }={\displaystyle \frac{e^{i\eta }}{\sqrt{2}}}\left[\begin{array}{cc}e^{i(\varphi +\psi )}& e^{i(\varphi \psi )}\\ e^{i(\varphi \psi )}& e^{i(\varphi +\psi )}\end{array}\right]`$
where $`\eta ,\varphi ,`$ and $`\psi `$ are real numbers (see pp.175-176 in Nielsen and Chuang (2000), for example). In particular, we see $`U_{0,0,0}=H`$.

However symmetry of the Hadamard walk depends heavily on the initial qubit state, see Konno, Namiki and Soshi (2002). Another generalization of the Hadamard walk is:

$`H(\rho )=\left[\begin{array}{cc}\sqrt{\rho }& \sqrt{1\rho }\\ \sqrt{1\rho }& \sqrt{\rho }\end{array}\right]`$
where $`0\rho 1`$. Note that $`\rho =1/2`$ is the Hadamard walk, that is, $`H=H(1/2)`$.

In the present paper, the study on the dependence of some important properties and quantities (e.g., symmetry of distribution, limit distribution, absorption probability) on initial qubit state is one of the essential parts, so we define the collection of initial qubit states as follows:

$$\mathrm{\Phi }=\{\phi =\left[\begin{array}{cc}\alpha & \\ \beta & \end{array}\right]𝐂^2:|\alpha |^2+|\beta |^2=1\}$$

Let $`X_{j,n}^\phi `$ be the $`j`$-type QRW at time $`n`$ starting from initial qubit state $`\phi \mathrm{\Phi }`$ with $`X_{j,0}^\phi =0`$ for $`j=A,G`$. In our treatment of QRWs, as well as the matrices $`P_j`$ and $`Q_j`$, it is convenient to introduce

$$R_A=\left[\begin{array}{cc}c& d\\ 0& 0\end{array}\right],S_A=\left[\begin{array}{cc}0& 0\\ a& b\end{array}\right]$$

and

$$R_G=\left[\begin{array}{cc}0& a\\ 0& c\end{array}\right],S_G=\left[\begin{array}{cc}b& 0\\ d& 0\end{array}\right]$$

We should remark that both A-type and G-type $`P_j,Q_j,R_j,S_j(j=A,G)`$ form an orthonormal basis of the Hilbert space $`M_2(𝐂)`$ which is the vector space of complex $`2\times 2`$ matrices with respect to the trace inner product $`A|B=`$ tr $`(A^{}B)`$. Therefore we can express any $`2\times 2`$ matrix $`X`$ conveniently in the form,

(2.10) $`X=\text{tr}(P_j^{}X)P_j+\text{tr}(Q_j^{}X)Q_j+\text{tr}(R_j^{}X)R_j+\text{tr}(S_j^{}X)S_j`$
for each $`j=A,G`$. We call the analysis based on $`P_j,Q_j,R_j,S_j(j=A,G)`$ the PQRS method.

The $`n\times n`$ unit and zero matrices are written $`I_n`$ and $`O_n`$ respectively. For instance, if $`X=I_2`$, then Eq. (2.2) gives

(2.11) $`I_2=\overline{a}P_j+\overline{d}Q_j+\overline{c}R_j+\overline{b}S_j`$
for each $`j=A,G`$. The next table of products of $`P_j,Q_j,R_j,S_j(j=A,G)`$ is very useful in computing some quantities:

|  | $`P_j`$ | $`Q_j`$ | $`R_j`$ | $`S_j`$ |
| --- | --- | --- | --- | --- |
| $`P_j`$ | $`aP_j`$ | $`bR_j`$ | $`aR_j`$ | $`bP_j`$ |
| $`Q_j`$ | $`cS_j`$ | $`dQ_j`$ | $`cQ_j`$ | $`dS_j`$ |
| $`R_j`$ | $`cP_j`$ | $`dR_j`$ | $`cR_j`$ | $`dP_j`$ |
| $`S_j`$ | $`aS_j`$ | $`bQ_j`$ | $`aQ_j`$ | $`bS_j`$ |

where $`P_jQ_j=bR_j`$, for example. We should remark that the algebraic structure for both types is the same.

In order to consider absorption problems stated in Section 4, now we describe the evolution and measurement of the A-type QRW starting from location $`k`$ on $`\{0,1,\mathrm{},N\}`$ with absorbing boundaries (for examples, see Ambainis et al. (2001), Bach et al. (2002), and Kempe (2002) for more detailed information). As for the G-type QRW, we can define its evolution and measurement in a similar fashion.

First we consider $`N=\mathrm{}`$ case. In this case, an absorbing boundary is present at location 0. The evolution mechanism is described as follows:

Step 1. Initialize the system $`\phi \mathrm{\Phi }`$ at location $`k`$.

Step 2. (a) Apply Eq. (2.1) to one step time evolution. (b) Measure the system to see where it is at location 0 or not.

Step 3. If the result of measurement revealed that the system was at location 0, then terminate the process, otherwise repeat step 2.

In this setting, let $`\mathrm{\Xi }_{A,k}^{(\mathrm{})}(n)`$ be the sum over possible paths for which the particle first hits 0 at time $`n`$ starting from $`k`$ for the A-type QRW. For example,

$`\mathrm{\Xi }_{A,1}^{(\mathrm{})}(5)=P_A^2Q_AP_AQ_A+P_A^3Q_A^2=(ab^2c+a^2bd)R_A`$
The probability that the particle first hits 0 at time $`n`$ starting from $`k`$ is

$`P_{A,k}^{(\mathrm{})}(n;\phi )=|\mathrm{\Xi }_{A,k}^{(\mathrm{})}(n)\phi |^2`$
So the probability that the particle first hits 0 starting from $`k`$ for the A-type QRW is

$`P_{A,k}^{(\mathrm{})}(\phi )={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}P_{A,k}^{(\mathrm{})}(n;\phi )`$
Next we consider $`N<\mathrm{}`$ case. This case is similar to the $`N=\mathrm{}`$ case, except that two absorbing boundaries are presented at locations 0 and $`N`$ as follows:

Step 1. Initialize the system $`\phi \mathrm{\Phi }`$ at location $`k`$.

Step 2. (a) Apply Eq. (2.1) to one step time evolution. (b) Measure the system to see where it is at location 0 or not. (c) Measure the system to see where it is at location $`N`$ or not.

Step 3. If the result of either measurement revealed that the system was either at location 0 or location $`N`$, then terminate the process, otherwise repeat step 2.

Let $`\mathrm{\Xi }_{A,k}^{(N)}(n)`$ be the sum over possible paths for which the particle first hits 0 at time $`n`$ starting from $`k`$ before it arrives at $`N`$ for the A-type QRW. For example,

(2.12) $`\mathrm{\Xi }_{A,1}^{(3)}(5)=P_A^2Q_AP_AQ_A=ab^2cR_A`$
In a similar way, we can define $`P_{A,k}^{(N)}(n;\phi )`$ and $`P_{A,k}^{(N)}(\phi )`$.

## 3 Limit Theorem

This section treats limit theorems for both A-type and G-type QRWs $`X_{j,n}^\phi (j=A,G).`$ To study $`P(X_{j,n}^\phi =k)`$ for $`n+k=`$ even, it suffices to understand the following combinatorial expression. For fixed $`l`$ and $`m`$ with $`l+m=n`$ and $`ml=k`$, we consider

$$\mathrm{\Xi }_j(l,m)=\underset{l_i,m_i0:m_1+\mathrm{}+m_n=m,l_1+\mathrm{}+l_n=l}{}P_j^{l_1}Q_j^{m_1}P_j^{l_2}Q_j^{m_2}\mathrm{}P_j^{l_n}Q_j^{m_n}$$

since

$$P(X_{j,n}^\phi =k)=|\mathrm{\Xi }_j(l,m)\phi |^2$$

Since $`P_j,Q_j,R_j,S_j(j=A,G)`$ are a basis of $`M_2(𝐂)`$, $`\mathrm{\Xi }_j(l,m)`$ has the following form:

$$\mathrm{\Xi }_j(l,m)=p_j(l,m)P_j+q_j(l,m)Q_j+r_j(l,m)R_j+s_j(l,m)S_j$$

Next problem is to obtain explicit forms of $`p_j(l,m),q_j(l,m),r_j(l,m),`$ and $`s_j(l,m)`$. In the above example of $`n=l+m=4`$ case, we see that for $`j=A,G,`$

$`\mathrm{\Xi }_j(4,0)=a^3P_j,\mathrm{\Xi }_j(3,1)=2abcP_j+a^2bR_j+a^2cS_j,`$
$`\mathrm{\Xi }_j(2,2)=bcdP_j+abcQ_j+abdR_j+bc^2S_j,`$
$`\mathrm{\Xi }_j(1,3)=2bcdQ_j+bd^2R_j+cd^2S_j,\mathrm{\Xi }_j(2,2)=d^3Q_j`$
So, for example,

$$p_j(3,1)=2abc,q_j(3,1)=0,r_j(3,1)=a^2b,s_j(3,1)=a^2c$$

Note that $`p_j(l,m),q_j(l,m),r_j(l,m),s_j(l,m)`$ do not depend on the type $`j`$.

In $`abcd=0`$ case, the argument is much easier. So from now on we focus only on $`abcd0`$ case. In this case, the next key lemma is obtained by a combinatorial method.

###### LEMMA 1

We consider both A-type and G-type QRWs with $`abcd0`$. Suppose that $`l,m0`$ with $`l+m=n`$, then we have

(i) for $`lm(=\mathrm{min}\{l,m\})1`$,

$`\mathrm{\Xi }_j(l,m)=a^ld^m{\displaystyle \underset{\gamma =1}{\overset{lm}{}}}\left({\displaystyle \frac{|b|^2}{|a|^2}}\right)^\gamma \left({\displaystyle \genfrac{}{}{0pt}{}{l1}{\gamma 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{m1}{\gamma 1}}\right)\times \left[{\displaystyle \frac{l\gamma }{a\gamma }}P_j+{\displaystyle \frac{m\gamma }{d\gamma }}Q_j+{\displaystyle \frac{1}{c}}R_j+{\displaystyle \frac{1}{b}}S_j\right]`$
(ii) for $`l(=n)1,m=0`$,

$$\mathrm{\Xi }_j(l,0)=a^{l1}P_j$$

(iii) for $`l=0,m(=n)1`$,

$$\mathrm{\Xi }_j(0,m)=d^{m1}Q_j$$

where $`j=A,G`$.

The proofs of parts (ii) and (iii) are trivial. The proof of part (i) is based on a consequence of enumerating the paths of drift $`l+m=n`$. To do so, we consider the following 4 cases:

$`\stackrel{w_1}{\stackrel{}{P_jP_j\mathrm{}P_j}}\stackrel{w_2}{\stackrel{}{Q_jQ_j\mathrm{}Q_j}}\stackrel{w_3}{\stackrel{}{P_jP_j\mathrm{}P_j}}\mathrm{}\stackrel{w_{2\gamma }}{\stackrel{}{Q_jQ_j\mathrm{}Q_j}}\stackrel{w_{2\gamma +1}}{\stackrel{}{P_jP_j\mathrm{}P_j}}`$
$`\stackrel{w_1}{\stackrel{}{Q_jQ_j\mathrm{}Q_j}}\stackrel{w_2}{\stackrel{}{P_jP_j\mathrm{}P_j}}\stackrel{w_3}{\stackrel{}{Q_jQ_j\mathrm{}Q_j}}\mathrm{}\stackrel{w_{2\gamma }}{\stackrel{}{P_jP_j\mathrm{}P_j}}\stackrel{w_{2\gamma +1}}{\stackrel{}{Q_jQ_j\mathrm{}Q_j}}`$
$`\stackrel{w_1}{\stackrel{}{P_jP_j\mathrm{}P_j}}\stackrel{w_2}{\stackrel{}{Q_jQ_j\mathrm{}Q_j}}\stackrel{w_3}{\stackrel{}{P_jP_j\mathrm{}P_j}}\mathrm{}\stackrel{w_{2\gamma }}{\stackrel{}{Q_jQ_j\mathrm{}Q_j}}`$
$`\stackrel{w_1}{\stackrel{}{Q_jQ_j\mathrm{}Q_j}}\stackrel{w_2}{\stackrel{}{P_jP_j\mathrm{}P_j}}\stackrel{w_3}{\stackrel{}{Q_jQ_j\mathrm{}Q_j}}\mathrm{}\stackrel{w_{2\gamma }}{\stackrel{}{P_jP_j\mathrm{}P_j}}`$
where $`w_1,w_2,\mathrm{},w_{2\gamma +1}1`$ and $`\gamma 1`$. The above each type of paths of $`P_j`$ and $`Q_j`$ corresponds to the each term of $`P_j,Q_j,R_j,S_j`$ respectively in the right hand side of the equation of part (i). As for the A-type case, the details of this proof appear in Konno (2002b). A similar proof also can be seen in Appendix A of Brun, Carteret and Ambainis (2002b). The proof of the G-type case is the same.

By this lemma, the characteristic function of $`X_{j,n}^\phi (j=A,G)`$ for $`abcd0`$ case is obtained. Moreover, the $`m`$th moment of $`X_{j,n}^\phi `$ can be derived from the characteristic function in the standard fashion. Here we give only the result of the $`m`$th moment.

###### PROPOSITION 2

We consider both A-type and G-type QRWs with $`abcd0`$.

(i) When $`m`$ is odd, we have

$`E((X_{j,n}^\phi )^m)`$ $`=`$ $`|a|^{2(n1)}[n^m\mathrm{\Gamma }_j`$
$`+{\displaystyle \underset{k=1}{\overset{\left[\frac{n1}{2}\right]}{}}}{\displaystyle \underset{\gamma =1}{\overset{k}{}}}{\displaystyle \underset{\delta =1}{\overset{k}{}}}\left({\displaystyle \frac{|b|^2}{|a|^2}}\right)^{\gamma +\delta }\left({\displaystyle \genfrac{}{}{0pt}{}{k1}{\gamma 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{k1}{\delta 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{nk1}{\gamma 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{nk1}{\delta 1}}\right)`$
$`\times {\displaystyle \frac{(n2k)^{m+1}}{\gamma \delta }}\mathrm{\Delta }_j]`$
where

$`\mathrm{\Gamma }_A=\left(|a|^2|b|^2\right)\left(|\alpha |^2|\beta |^2\right)+2(a\alpha \overline{b\beta }+\overline{a\alpha }b\beta ),\mathrm{\Gamma }_G=|\alpha |^2|\beta |^2,`$
$`\mathrm{\Delta }_A=\{n(|a|^2|b|^2)+\gamma +\delta \}(|\alpha |^2|\beta |^2)+\left({\displaystyle \frac{\gamma +\delta }{|b|^2}}2n\right)\mathrm{\Theta }_A,`$
$`\mathrm{\Delta }_G=(\gamma +\delta n)(|\alpha |^2|\beta |^2){\displaystyle \frac{\gamma +\delta }{|b|^2}}\mathrm{\Theta }_G,`$
$`\mathrm{\Theta }_A=a\alpha \overline{b\beta }+\overline{a\alpha }b\beta ,\mathrm{\Theta }_G=a\beta \overline{c\alpha }+\overline{a\beta }c\alpha `$
(ii) When $`m`$ is even, we have

$`E((X_{j,n}^\phi )^m)`$ $`=`$ $`|a|^{2(n1)}[n^m`$
$`+{\displaystyle \underset{k=1}{\overset{\left[\frac{n1}{2}\right]}{}}}{\displaystyle \underset{\gamma =1}{\overset{k}{}}}{\displaystyle \underset{\delta =1}{\overset{k}{}}}\left({\displaystyle \frac{|b|^2}{|a|^2}}\right)^{\gamma +\delta }\left({\displaystyle \genfrac{}{}{0pt}{}{k1}{\gamma 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{k1}{\delta 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{nk1}{\gamma 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{nk1}{\delta 1}}\right)`$
$`\times {\displaystyle \frac{(n2k)^m}{\gamma \delta }}\{(nk)^2+k^2n(\gamma +\delta )+{\displaystyle \frac{2\gamma \delta }{|b|^2}}\}]`$
It should be noted that when $`m`$ is even, $`E((X_{j,n}^\phi )^m)`$ is independent of initial qubit state $`\phi `$ and types $`j=A,G`$. In particular, we use the above result of $`m=1`$ case in order to study symmetry of distributions for the QRWs. Moreover we have the following new type of limit theorem (as for the A-type QRW, see Konno (2002a)):

###### THEOREM 3

We consider both A-type and G-type QRWs with $`abcd0`$. Let $`\mathrm{\Theta }_A=a\alpha \overline{b\beta }+\overline{a\alpha }b\beta ,`$ and $`\mathrm{\Theta }_G=a\beta \overline{c\alpha }+\overline{a\beta }c\alpha `$. If $`n\mathrm{}`$, then

$$\frac{X_{j,n}^\phi }{n}Z_j^\phi (j=A,G)$$

where $`Z_j^\phi `$ has a density

$$f_j(x;{}_{}{}^{t}[\alpha ,\beta ])=\frac{\sqrt{1|a|^2}}{\pi (1x^2)\sqrt{|a|^2x^2}}\{1(|\alpha |^2|\beta |^2+\frac{\mathrm{\Theta }_j}{|a|^2})x\}$$

for $`x(|a|,|a|)`$ with

$`E(Z_j^\phi )=\left(|\alpha |^2|\beta |^2+{\displaystyle \frac{\mathrm{\Theta }_j}{|a|^2}}\right)\times (1\sqrt{1|a|^2})`$
$`E((Z_j^\phi )^2)=1\sqrt{1|a|^2}`$
and $`Y_nY`$ means that $`Y_n`$ converges in distribution to a limit $`Y`$.

We remark that standard deviation of $`Z_j^\phi `$ is not independent of initial qubit state $`\phi ={}_{}{}^{t}[\alpha ,\beta ]`$. The above limit theorem suggests the following result on symmetry of distribution for both A-type and G-type QRWs. This is a generalization of the result give by Konno, Namiki and Soshi (2002) for the A-type Hadamard walk and by Konno (2002a) for the A-type QRW. For $`j=A,G,`$ define

$`\mathrm{\Phi }_{j,s}`$ $`=`$ $`\{\phi \mathrm{\Phi }:P(X_{j,n}^\phi =k)=P(X_{j,n}^\phi =k)\text{for any}n𝐙_+\text{and}k𝐙\}`$
$`\mathrm{\Phi }_{j,0}`$ $`=`$ $`\{\phi \mathrm{\Phi }:E(X_{j,n}^\phi )=0\text{for any}n𝐙_+\}`$
$`\mathrm{\Phi }_{j,}`$ $`=`$ $`\{\phi ={}_{}{}^{t}[\alpha ,\beta ]\mathrm{\Phi }:|\alpha |=|\beta |,\mathrm{\Theta }_j=0\}(j=A,G)`$
and $`𝐙`$ (resp. $`𝐙_+`$) is the set of (resp. non-negative) integers. For $`\phi \mathrm{\Phi }_s`$, the probability distribution of $`X_{j,n}^\phi `$ is symmetric for any $`n𝐙_+`$. Using the explicit form of $`E(X_{j,n}^\phi )`$ given by Proposition 2 (i) ($`m=1`$ case), we have

###### THEOREM 4

We consider both A-type and G-type QRWs with $`abcd0`$. Then for $`j=A,G`$ we have

$$\mathrm{\Phi }_{j,s}=\mathrm{\Phi }_{j,0}=\mathrm{\Phi }_{j,}$$

From now on we give an outline of our proof of Theorem 3 (for more details, see Konno (2002b)). To do so, we introduce the Jacobi polynomial $`P_n^{\nu ,\mu }(x)`$, where $`P_n^{\nu ,\mu }(x)`$ is orthogonal on $`[1,1]`$ with respect to $`(1x)^\nu (1+x)^\mu `$ with $`\nu ,\mu >1`$. Then the following relation holds:

$$P_n^{\nu ,\mu }(x)=\frac{\mathrm{\Gamma }(n+\nu +1)}{\mathrm{\Gamma }(n+1)\mathrm{\Gamma }(\nu +1)}{}_{2}{}^{}F_{1}^{}(n,n+\nu +\mu +1;\nu +1;(1x)/2)$$

where $`{}_{2}{}^{}F_{1}^{}(a,b;c;z)`$ is the hypergeometric series and $`\mathrm{\Gamma }(z)`$ is the gamma function. In general, as for orthogonal polynomials, see Koekoek and Swarttouw (1996). Remark that

$`{\displaystyle \underset{\gamma =1}{\overset{k}{}}}\left({\displaystyle \frac{|b|^2}{|a|^2}}\right)^{\gamma 1}{\displaystyle \frac{1}{\gamma }}\left({\displaystyle \genfrac{}{}{0pt}{}{k1}{\gamma 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{nk1}{\gamma 1}}\right)={\displaystyle \frac{1}{k}}|a|^{2(k1)}P_{k1}^{1,n2k}(2|a|^21)`$
$`{\displaystyle \underset{\gamma =1}{\overset{k}{}}}\left({\displaystyle \frac{|b|^2}{|a|^2}}\right)^{\gamma 1}\left({\displaystyle \genfrac{}{}{0pt}{}{k1}{\gamma 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{nk1}{\gamma 1}}\right)=|a|^{2(k1)}P_{k1}^{0,n2k}(2|a|^21)`$
By using the above equations and an expression of the characteristic function for $`X_{j,n}^\phi `$, we obtain the next asymptotics of characteristic function $`E(e^{i\xi X_{j,n}^\phi /n})`$: if $`n\mathrm{}`$ with $`k/n=x((1|a|)/2,(1+|a|)/2)`$, then

$`E(e^{i\xi X_{j,n}^\phi /n}){\displaystyle \underset{k=1}{\overset{\left[\frac{n1}{2}\right]}{}}}|a|^{2n4k2}|b|^4`$
$`\times [\{{\displaystyle \frac{2x^22x+1}{x^2}}(P_{k1}^{1,n2k})^2{\displaystyle \frac{2}{x}}P_{k1}^{1,n2k}P_{k1}^{0,n2k}+{\displaystyle \frac{2}{|b|^2}}(P_{k1}^{0,n2k})^2\}\mathrm{cos}((12x)\xi )`$
$`{\displaystyle \frac{(12x)\mathrm{\Gamma }_j}{x^2}}(P_{k1}^{1,n2k})^22I_j(|\alpha |^2|\beta |^2{\displaystyle \frac{\mathrm{\Theta }_j}{|b|^2}})P_{k1}^{0,n2k}P_{k1}^{1,n2k}\left\}i\mathrm{sin}((12x)\xi )\right]`$
where $`f(n)g(n)`$ means $`f(n)/g(n)1(n\mathrm{}),P_{k1}^{i,n2k}=P_{k1}^{i,n2k}(2|a|^21)(i=0,1),`$ and $`I_A=1,I_G=1.`$

Next we prepare the following asymptotic results for the Jacobi polynomial $`P_n^{\alpha +an,\beta +bn}(x)`$ derived by Chen and Ismail (1991): if $`n\mathrm{}`$ with $`k/n=x((1|a|)/2,(1+|a|)/2)`$, then

$`P_{k1}^{0,n2k}{\displaystyle \frac{2|a|^{2kn}}{\sqrt{\pi n\sqrt{\mathrm{\Delta }}}}}\mathrm{cos}(An+B)`$
$`P_{k1}^{1,n2k}{\displaystyle \frac{2|a|^{2kn}}{\sqrt{\pi n\sqrt{\mathrm{\Delta }}}}}\sqrt{{\displaystyle \frac{x}{(1x)(1|a|^2)}}}\mathrm{cos}(An+B+\theta )`$
where $`\mathrm{\Delta }=(1|a|^2)(4x^24x+1|a|^2)`$, $`A`$ and $`B`$ are some constants (which are independent of $`n`$), and $`\theta [0,\pi /2]`$ is determined by $`\mathrm{cos}\theta =\sqrt{(1|a|^2)/4x(1x)}`$.

Combining the above observations with the Riemann-Lebesgue lemma (see p.462 in Durrett (1996)), we see that, if $`n\mathrm{}`$, then

$`E(e^{i\xi \frac{X_n^\phi }{n}})`$
$`{\displaystyle \frac{1|a|^2}{\pi }}{\displaystyle _{\frac{1|a|}{2}}^{\frac{1}{2}}}𝑑x{\displaystyle \frac{1}{x(1x)\sqrt{(|a|^21)(4x^24x+1|a|^2)}}}`$
$`\times \left[\mathrm{cos}((12x)\xi )(12x)\left(|\alpha |^2|\beta |^2+{\displaystyle \frac{\mathrm{\Theta }_j}{|a|^2}}\right)i\mathrm{sin}((12x)\xi )\right]`$
$`={\displaystyle _{|a|}^{|a|}}{\displaystyle \frac{\sqrt{1|a|^2}}{\pi (1x^2)\sqrt{|a|^2x^2}}}\left\{1\left(|\alpha |^2|\beta |^2+{\displaystyle \frac{\mathrm{\Theta }_j}{|a|^2}}\right)x\right\}e^{i\xi x}𝑑x`$
$`=\varphi (\xi )`$
Then $`\varphi (\xi )`$ is continuous at $`\xi =0`$, so the continuity theorem (see p.99 in Durrett (1996)) implies that $`X_n^\phi /n`$ converges in distribution to the $`Z^\phi `$ with characteristic function $`\varphi `$. Therefore Theorem 3 is obtained.

Muraki (2002) introduced a notion of quasi-universal product for algebraic probability spaces and showed that there exist only five quasi-universal products, that is, tensor product (case 1), free product (case 2), Boolean product (case 3), monotone product (case 4), and anti-monotone product (case 5). Algebraic central limit theorems describe limit behaviors of rescaled sum of algebraic random variables $`\sigma _1,\sigma _2,\mathrm{}`$ with mean 0 and variance 1,

$$\frac{\sigma _1+\sigma _2+\mathrm{}+\sigma _n}{\sqrt{n}}$$

converges weakly, in the limit $`n\mathrm{}`$, to the Gaussian distribution $`e^{x^2/2}/\sqrt{2\pi }`$ (case 1), the semi-circle distribution $`\chi _{[2,2]}(x)\sqrt{4x^2}/2\pi `$ (case 2), the Bernoulli distribution $`(\delta _1+\delta _1)/2`$ (case 3), and arcsine distribution $`\chi _{(\sqrt{2},\sqrt{2})}(x)/\pi \sqrt{2x^2}`$ (cases 4 and 5), where $`\chi (A)=1`$ if $`xA`$, $`=0`$ if $`xA`$ with $`A𝐑`$ (see Muraki (2001), Hashimoto (2002), for examples). So our limit theorem (Theorem 3) would belong to another category which is different from ones given by Muraki.

Now we compare our analytical results (given by Theorem 3) with the numerical ones for the A-type Haramard walk.

We see that Theorem 3 implies that if $`\sqrt{2}/2<a<b<\sqrt{2}/2`$, then as $`n\mathrm{}`$,

$`P(aX_{A,n}^\phi /nb){\displaystyle _a^b}{\displaystyle \frac{1(|\alpha |^2|\beta |^2+\mathrm{\Theta }_A)x}{\pi (1x^2)\sqrt{12x^2}}}𝑑x`$
for any initial qubit state $`\phi ={}_{}{}^{t}[\alpha ,\beta ]`$. For the symmetric CRW $`Y_n^o`$ starting from the origin, the well-known central limit theorem implies that if $`\mathrm{}<a<b<\mathrm{}`$, then as $`n\mathrm{}`$,

$`P(aY_n^o/\sqrt{n}b){\displaystyle _a^b}{\displaystyle \frac{e^{x^2/2}}{\sqrt{2\pi }}}𝑑x`$
This result is often called the de Moivre-Laplace theorem. When we take $`\phi ={}_{}{}^{t}[1/\sqrt{2},i/\sqrt{2}]`$ (symmetric case), then we have the following QRW version of the de Moivre-Laplace theorem: if $`\sqrt{2}/2<a<b<\sqrt{2}/2`$, then as $`n\mathrm{}`$,

$`P(aX_{A,n}^\phi /nb){\displaystyle _a^b}{\displaystyle \frac{1}{\pi (1x^2)\sqrt{12x^2}}}𝑑x`$
So there is a difference between the QRW $`X_{A,n}^\phi `$ and the CRW $`Y_n^o`$ even in a symmetric case for $`\phi ={}_{}{}^{t}[1/\sqrt{2},i/\sqrt{2}]`$.

Noting that $`E(X_{A,n}^\phi )=0(n0)`$ for any $`\phi \mathrm{\Phi }_{A,s}`$, we have

$`sd(X_{A,n}^\phi )/n\sqrt{(2\sqrt{2})/2}=0.54119\mathrm{}`$
where $`sd(X)`$ is the standard deviation of $`X`$. This rigorous result reveals that numerical simulation result 3/5 = 0.6 given by Travaglione and Milburn (2002) is not so accurate.

As in a similar way, when we take $`\phi ={}_{}{}^{t}[0,e^{i\theta }]`$ where $`\theta [0,2\pi )`$ (asymmetric case), we see that if $`\sqrt{2}/2<a<b<\sqrt{2}/2`$, then as $`n\mathrm{}`$,

$`P(aX_{A,n}^\phi /nb){\displaystyle _a^b}{\displaystyle \frac{1}{\pi (1x)\sqrt{12x^2}}}𝑑x`$
So we have

$`E(X_{A,n}^\phi )/n(2\sqrt{2})/2=0.29289\mathrm{},sd(X_{A,n}^\phi )/n\sqrt{(\sqrt{2}1)/2}=0.45508\mathrm{}`$
When $`\phi ={}_{}{}^{t}[0,1]`$ ($`\theta =0`$), Ambainis et al. (2001) gave the same result. In the paper, they took two approaches, that is, the Schrödinger approach and the path integral approach. However their result comes mainly from the Schrödinger approach by using a Fourier analysis.

In another asymmetric case $`\phi ={}_{}{}^{t}[e^{i\theta },0]`$ where $`\theta [0,2\pi )`$, a similar argument implies that if $`\sqrt{2}/2<a<b<\sqrt{2}/2`$, then as $`n\mathrm{}`$,

$`P(aX_{A,n}^\phi /nb){\displaystyle _a^b}{\displaystyle \frac{1}{\pi (1+x)\sqrt{12x^2}}}𝑑x`$
Note that $`f_A(x;{}_{}{}^{t}[e^{i\theta },0])=f_A(x;{}_{}{}^{t}[0,e^{i\theta }])`$ for any $`x(\sqrt{2}/2,\sqrt{2}/2)`$. Therefore concerning the $`m`$th moment of the limit distribution, we have the same result as in the previous case $`\phi ={}_{}{}^{t}[0,e^{i\theta }]`$. So the standard deviation of the limit distribution $`Z_A^\phi `$ is given by $`\sqrt{(\sqrt{2}1)/2}=0.45508\mathrm{}.`$ Simulation result 0.4544 $`\pm `$ 0.0012 in Mackay et al. (2002) (their case is $`\theta =0`$) is consistent with our rigorous result.

## 4 Absorption Problem

From now on we consider absorption problems for both A-type and G-type QRWs located on the sets $`\{0,1,\mathrm{},N\}`$ or $`\{0,1,\mathrm{}\}`$. Results in this section for the A-type QRW appear in Konno, Namiki, Soshi and Sudbury (2003).

Before we move to a quantum case, first we describe the CRW on a finite set $`\{0,1,\mathrm{},N\}`$ with two absorption barriers at locations $`0`$ and $`N`$ (see Doyle and Snell (1984), Grimmett and Stirzaker (1992), for examples). The particle moves at each step either one unit to the left with probability $`p`$ or one unit to the right with probability $`q=1p`$ until it hits one of the absorption barriers. The directions of different steps are independent of each other. The CRW starting from $`k\{0,1,\mathrm{},N\}`$ at time $`n`$ is denoted by $`Y_n^k`$ here. Let

$$T_m=\mathrm{min}\{n0:Y_n^k=m\}$$

be the time of the first visit to $`m\{0,1,\mathrm{},N\}`$. Using the subscript $`k`$ to indicate $`Y_0^k=k`$, we let

$$P_k^{(N)}=P_k(T_0<T_N)$$

be the probability that the particle hits $`0`$ starting from $`k`$ before it arrives at $`N`$. The absorption problem is also known as the Gambler’s ruin problem.

Now we review some known results and conjectures on absorption problems related to this paper for A-type QRWs.

In the case of $`U=H`$ (the Hadamard walk), when $`N=\mathrm{}`$, that is, the state space is $`\{0,1,\mathrm{}\}`$ case, Ambainis et al. (2001) proved

(4.13) $`P_{A,1}^{(\mathrm{})}({}_{}{}^{t}[0,1])=P_{A,1}^{(\mathrm{})}({}_{}{}^{t}[1,0])={\displaystyle \frac{2}{\pi }}`$
and Bach et al. (2002) showed

$`\underset{k\mathrm{}}{lim}P_{A,k}^{(\mathrm{})}(\phi )=\left({\displaystyle \frac{1}{2}}\right)|\alpha |^2+\left({\displaystyle \frac{2}{\pi }}{\displaystyle \frac{1}{2}}\right)|\beta |^2+2\left({\displaystyle \frac{1}{\pi }}{\displaystyle \frac{1}{2}}\right)\mathrm{}(\overline{\alpha }\beta )`$
for any initial qubit state $`\phi ={}_{}{}^{t}[\alpha ,\beta ]\mathrm{\Phi }.`$ Furthermore, in the case of $`U=H(\rho )`$, Bach et al. (2002) gave

$`\underset{k\mathrm{}}{lim}P_{A,k}^{(\mathrm{})}({}_{}{}^{t}[0,1])={\displaystyle \frac{\rho }{1\rho }}({\displaystyle \frac{\mathrm{cos}^1(12\rho )}{\pi }}1)+{\displaystyle \frac{2}{\pi \sqrt{1/\rho 1}}}`$
$`\underset{k\mathrm{}}{lim}P_{A,k}^{(\mathrm{})}({}_{}{}^{t}[1,0])={\displaystyle \frac{\mathrm{cos}^1(12\rho )}{\pi }}`$
The second result was conjectured by Yamasaki, Kobayashi and Imai (2002).

When $`N`$ is finite, the following conjecture by Ambainis et al. (2001) is still open for the $`U=H`$ case:

$`P_{A,1}^{(N+1)}({}_{}{}^{t}[0,1])={\displaystyle \frac{2P_{A,1}^{(N)}({}_{}{}^{t}[0,1])+1}{2P_{A,1}^{(N)}({}_{}{}^{t}[0,1])+2}}(N1),P_{A,1}^{(1)}({}_{}{}^{t}[0,1])=0.`$
Solving the above recurrence gives

(4.14) $`P_{A,1}^{(N)}({}_{}{}^{t}[0,1])={\displaystyle \frac{1}{\sqrt{2}}}\times {\displaystyle \frac{(3+2\sqrt{2})^{N1}1}{(3+2\sqrt{2})^{N1}+1}}(N1)`$
However in contrast with $`P_{A,1}^{(\mathrm{})}({}_{}{}^{t}[0,1])=2/\pi `$, Ambainis et al. (2001) proved

$$\underset{N\mathrm{}}{lim}P_{A,1}^{(N)}({}_{}{}^{t}[0,1])=1/\sqrt{2}$$

Let $`T_0`$ be the first hitting time to 0. We consider the conditional $`m`$th moment of $`T_0`$ starting from $`k=1`$ given the event $`\{T_0<\mathrm{}\}`$, that is, $`E_{j,1}^{(\mathrm{})}((T_0)^m|T_0<\mathrm{})=E_{j,1}^{(\mathrm{})}((T_0)^m;T_0<\mathrm{})/P_{j,1}^{(\mathrm{})}(T_0<\mathrm{})`$ for the $`j`$-type QRW.

From now on we begin with the classical case. In this case, to obtain $`P_k^{(N)}`$, we use the following difference equation:

(4.15) $`P_k^{(N)}=pP_{k1}^{(N)}+qP_{k+1}^{(N)}(1kN1)`$
with boundary conditions:

(4.16) $`P_0^{(N)}=1,P_N^{(N)}=0`$
To consider a similar equation even in the quantum case and to use the PQRS method are our basic strategy.

From now on we focus on $`1kN1`$ case. So we consider only $`n1`$ case. Noting that $`\{P_j,Q_j,R_j,S_j\}`$ for each $`j=A,G`$ is a basis of $`M_2(𝐂)`$, $`\mathrm{\Xi }_{j,k}^{(N)}(n)`$ can be written as

$`\mathrm{\Xi }_{j,k}^{(N)}(n)=p_{j,k}^{(N)}(n)P_j+q_{j,k}^{(N)}(n)Q_j+r_{j,k}^{(N)}(n)R_j+s_{j,k}^{(N)}(n)S_j`$
From the definition of $`\mathrm{\Xi }_{j,k}^{(N)}(n)`$, it is easily shown that there exist only two types of paths, that is, $`P_j\mathrm{}P_j`$ and $`P_j\mathrm{}Q_j`$, since we consider only a hitting time to $`0`$ before it arrives at $`N`$. Therefore we see that $`q_{j,k}^{(N)}(n)=s_{j,k}^{(N)}(n)=0(n1)`$.

We assume $`N3`$. Noting that the definition of $`\mathrm{\Xi }_{j,k}^{(N)}(n)`$, we have

$`\mathrm{\Xi }_{j,k}^{(N)}(n)=\mathrm{\Xi }_{j,k1}^{(N)}(n1)P_j+\mathrm{\Xi }_{j,k+1}^{(N)}(n1)Q_j(1kN1)`$
The above equation is a QRW version of the difference equation, i.e., Eq. (4.7) for the CRW. As an example, see Eq. (2.4). Then we have

$`p_{j,k}^{(N)}(n)=ap_{j,k1}^{(N)}(n1)+cr_{j,k1}^{(N)}(n1)`$
$`r_{j,k}^{(N)}(n)=bp_{j,k+1}^{(N)}(n1)+dr_{j,k+1}^{(N)}(n1)`$
Note that the above equations do not depend on types of QRWs. Next we consider boundary conditions related to Eq. (4.8) in the classical case. When $`k=N`$,

$$P_{j,N}^{(N)}(0;\phi )=|\mathrm{\Xi }_{j,N}^{(N)}(0)\phi |^2=0$$

for any $`\phi \mathrm{\Phi }`$. So we take $`\mathrm{\Xi }_{j,N}^{(N)}(0)=O_2`$. In this case, Eq. (2.2) gives

$$p_{j,N}^{(N)}(0)=r_{j,N}^{(N)}(0)=0$$

If $`k=0`$, then

$$P_{j,0}^{(N)}(0;\phi )=|\mathrm{\Xi }_{j,0}^{(N)}(0)\phi |^2=1$$

for any $`\phi \mathrm{\Phi }`$. So we choose $`\mathrm{\Xi }_{j,N}^{(N)}(0)=I_2`$. From Eq. (2.4), we have

$$p_{j,0}^{(N)}(0)=\overline{a},r_{j,0}^{(N)}(0)=\overline{c}$$

It is noted that the above boundary conditions also do not depend on types of QRWs. Therefore, from now on we will omit subscript $`j`$ of $`p_{j,k}^{(N)}(n)`$ and $`r_{j,k}^{(N)}(n).`$ Let

$`v_k^{(N)}(n)=\left[\begin{array}{cc}p_k^{(N)}(n)& \\ r_k^{(N)}(n)& \end{array}\right]`$
Then we see that for $`n1`$ and $`1kN1`$,

(4.22) $`v_k^{(N)}(n)=\left[\begin{array}{cc}a& c\\ 0& 0\end{array}\right]v_{k1}^{(N)}(n1)+\left[\begin{array}{cc}0& 0\\ b& d\end{array}\right]v_{k+1}^{(N)}(n1)`$
and for $`1kN`$,

$`v_0^{(N)}(0)=\left[\begin{array}{cc}\overline{a}& \\ \overline{c}& \end{array}\right],v_k^{(N)}(0)=\left[\begin{array}{cc}0& \\ 0& \end{array}\right]`$
Moreover,

$`v_0^{(N)}(n)=v_N^{(N)}(n)=\left[\begin{array}{cc}0& \\ 0& \end{array}\right](n1)`$
So the definition of $`P_{j,k}^{(N)}(\phi )`$ gives

###### LEMMA 5

$`P_{j,k}^{(N)}(\phi )={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}P_{j,k}^{(N)}(n;\phi )`$
$`P_{j,k}^{(N)}(n;\phi )=C_{j,1}(n)|\alpha |^2+C_{j,2}(n)|\beta |^2+2\mathrm{}(C_{j,3}(n)\overline{\alpha }\beta )`$
where $`\mathrm{}(z)`$ is the real part of $`z𝐂`$, $`\phi ={}_{}{}^{t}[\alpha ,\beta ]\mathrm{\Phi }`$ and

$`C_{A,1}(n)=|ap_k^{(N)}(n)+cr_k^{(N)}(n)|^2,C_{A,2}(n)=|bp_k^{(N)}(n)+dr_k^{(N)}(n)|^2`$
$`C_{A,3}(n)=\overline{(ap_k^{(N)}(n)+cr_k^{(N)}(n))}(bp_k^{(N)}(n)+dr_k^{(N)}(n))`$
$`C_{G,1}(n)=|p_k^{(N)}(n)|^2,C_{G,2}(n)=|r_k^{(N)}(n)|^2,C_{G,3}(n)=\overline{p_k^{(N)}(n)}r_k^{(N)}(n)`$
To solve $`P_{j,k}^{(N)}(\phi )`$, we introduce generating functions of $`p_k^{(N)}(n)`$ and $`r_k^{(N)}(n)`$ as follows:

$`\stackrel{~}{p}_k^{(N)}(z)={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}p_k^{(N)}(n)z^n,\stackrel{~}{r}_k^{(N)}(z)={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}r_k^{(N)}(n)z^n`$
By Eq. (4.9), we have

$`\stackrel{~}{p}_k^{(N)}(z)=az\stackrel{~}{p}_{k1}^{(N)}(z)+cz\stackrel{~}{r}_{k1}^{(N)}(z)`$
$`\stackrel{~}{r}_k^{(N)}(z)=bz\stackrel{~}{p}_{k+1}^{(N)}(z)+dz\stackrel{~}{r}_{k+1}^{(N)}(z)`$
Solving these, we see that both $`\stackrel{~}{p}_k^{(N)}(z)`$ and $`\stackrel{~}{r}_k^{(N)}(z)`$ satisfy the same recurrence:

$`d\stackrel{~}{p}_{k+2}^{(N)}(z)\left(\mathrm{}z+{\displaystyle \frac{1}{z}}\right)\stackrel{~}{p}_{k+1}^{(N)}(z)+a\stackrel{~}{p}_k^{(N)}(z)=0`$
$`d\stackrel{~}{r}_{k+2}^{(N)}(z)\left(\mathrm{}z+{\displaystyle \frac{1}{z}}\right)\stackrel{~}{r}_{k+1}^{(N)}(z)+a\stackrel{~}{r}_k^{(N)}(z)=0`$
From the characteristic equations with respect to the above recurrences, we have the same root: if $`a0`$, then

$`\lambda _\pm ={\displaystyle \frac{\mathrm{}z^2+1\sqrt{\mathrm{}^2z^4+2\mathrm{}(12|a|^2)z^2+1}}{2\mathrm{}\overline{a}z}}`$
where $`\mathrm{}=detU=adbc`$.

From now on we consider mainly $`U=H`$ (the Hadamard walk) case with $`N=\mathrm{}`$. Remark that the definition of $`\mathrm{\Xi }_{j,1}^{(\mathrm{})}(n)`$ gives $`p_1^{(\mathrm{})}(n)=0(n2)`$ and $`p_1^{(\mathrm{})}(1)=1`$. So we have $`\stackrel{~}{p}_1^{(\mathrm{})}(z)=z`$. Moreover noting $`lim_k\mathrm{}\stackrel{~}{p}_k^{(\mathrm{})}(z)<\mathrm{}`$, the following explicit form is obtained:

$`\stackrel{~}{p}_k^{(\mathrm{})}(z)=z\lambda _+^{k1},\stackrel{~}{r}_k^{(\mathrm{})}(z)={\displaystyle \frac{1+\sqrt{z^4+1}}{z}}\lambda _+^{k1}`$
where $`\lambda _\pm =(z^21\pm \sqrt{z^4+1})/\sqrt{2}z.`$ Therefore for $`k=1`$, we have $`\stackrel{~}{r}_1^{(\mathrm{})}(z)=(1+\sqrt{z^4+1})/z`$. By using these, we obtain

###### PROPOSITION 6

For each $`j`$-type QRW, we have

(4.25) $`P_{A,1}^{(\mathrm{})}(\phi )={\displaystyle \frac{2}{\pi }}+2\left(1{\displaystyle \frac{2}{\pi }}\right)\mathrm{}(\overline{\alpha }\beta ),`$
(4.26) $`P_{G,1}^{(\mathrm{})}(\phi )=|\alpha |^2+\left({\displaystyle \frac{4}{\pi }}1\right)|\beta |^2`$
for any initial qubit state $`\phi ={}_{}{}^{t}[\alpha ,\beta ]\mathrm{\Phi }`$.

Eqs. (4.10) and (4.11) are a generalization of Eq. (4.5) given by Ambainis et al. (2001). Then we get $`(4\pi )/\pi P_{j,1}^{(\mathrm{})}(\phi )1`$ for each $`j`$-type QRW. Moreover we obtain a result on the conditional $`m`$th moment of the first hitting time to 0 starting from $`k=1`$ given an event $`\{T_0<\mathrm{}\}`$, that is, $`E_{j,1}^{(\mathrm{})}((T_0)^m|T_0<\mathrm{})=E_{j,1}^{(\mathrm{})}((T_0)^m;T_0<\mathrm{})/P_{j,1}^{(\mathrm{})}(T_0<\mathrm{})`$ in a similar way.

###### PROPOSITION 7

For each $`j`$-type QRW, we have

$`E_{j,1}^{(\mathrm{})}(T_0|T_0<\mathrm{})={\displaystyle \frac{1}{P_{j,1}^{(\mathrm{})}(\phi )}}`$
$`E_{j,1}^{(\mathrm{})}((T_0)^m|T_0<\mathrm{})=\mathrm{}(m2)`$
Next we consider a finite $`N`$ case. By using boundary conditions: $`\stackrel{~}{p}_1^{(N)}(z)=z`$ and $`\stackrel{~}{r}_{N1}^{(N)}(z)=0`$, we see that $`\stackrel{~}{p}_k^{(N)}(z)`$ and $`\stackrel{~}{r}_k^{(N)}(z)`$ satisfy

(4.27) $`\stackrel{~}{p}_k^{(N)}(z)=\left({\displaystyle \frac{z}{2}}+E_z\right)\lambda _+^{k1}+\left({\displaystyle \frac{z}{2}}E_z\right)\lambda _{}^{k1}`$
(4.28) $`\stackrel{~}{r}_k^{(N)}(z)=C_z(\lambda _+^{kN+1}\lambda _{}^{kN+1})`$
where

(4.29) $`C_z={\displaystyle \frac{z^2}{\sqrt{2}}}(1)^{N2}(\lambda _+^{N3}\lambda _{}^{N3})`$
$`\times \left\{(\lambda _+^{N2}\lambda _{}^{N2})^2{\displaystyle \frac{z}{\sqrt{2}}}(\lambda _+^{N2}\lambda _{}^{N2})(\lambda _+^{N3}\lambda _{}^{N3})(1)^{N3}(\lambda _+\lambda _{})^2\right\}^1`$
(4.30) $`E_z={\displaystyle \frac{z}{2(\lambda _+^{N2}\lambda _{}^{N2})}}[2(1)^{N3}(\lambda _+\lambda _{})(\lambda _+^{N3}\lambda _{}^{N3})`$
$`\times \left\{(\lambda _+^{N2}\lambda _{}^{N2})^2{\displaystyle \frac{z}{\sqrt{2}}}(\lambda _+^{N2}\lambda _{}^{N2})(\lambda _+^{N3}\lambda _{}^{N3})(1)^{N3}(\lambda _+\lambda _{})^2\right\}^1`$
$`+(\lambda _+^{N2}+\lambda _{}^{N2})]`$
So we obtain

###### THEOREM 8

For each $`j`$-type QRW, we have

$`P_{j,k}^{(N)}(\phi )=C_{j,1}|\alpha |^2+C_{j,2}|\beta |^2+2\mathrm{}(C_{j,3}\overline{\alpha }\beta )`$
where $`\phi ={}_{}{}^{t}[\alpha ,\beta ]\mathrm{\Phi }`$ and

$`C_{A,1}={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}|a\stackrel{~}{p}_k^{(N)}(e^{i\theta })+c\stackrel{~}{r}_k^{(N)}(e^{i\theta })|^2𝑑\theta `$
$`C_{A,2}={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}|b\stackrel{~}{p}_k^{(N)}(e^{i\theta })+d\stackrel{~}{r}_k^{(N)}(e^{i\theta })|^2𝑑\theta `$
$`C_{A,3}={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}\overline{(a\stackrel{~}{p}_k^{(N)}(e^{i\theta })+c\stackrel{~}{r}_k^{(N)}(e^{i\theta }))}(b\stackrel{~}{p}_k^{(N)}(e^{i\theta })+d\stackrel{~}{r}_k^{(N)}(e^{i\theta }))𝑑\theta `$
$`C_{G,1}={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}|\stackrel{~}{p}_k^{(N)}(e^{i\theta })|^2𝑑\theta ,C_{G,2}={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}|\stackrel{~}{r}_k^{(N)}(e^{i\theta })|^2𝑑\theta `$
$`C_{G,3}={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}\overline{\stackrel{~}{p}_k^{(N)}(e^{i\theta })}\stackrel{~}{r}_k^{(N)}(e^{i\theta })𝑑\theta `$
with $`a=b=c=d=1/\sqrt{2}`$, here $`\stackrel{~}{p}_k^{(N)}(z)`$ and $`\stackrel{~}{r}_k^{(N)}(z)`$ satisfy Eqs. (4.12) and (4.13), and $`C_z`$ and $`E_z`$ satisfy Eqs. (4.14) and (4.15).

Here we consider $`U=H`$ (the Hadamard walk), $`\phi ={}_{}{}^{t}[\alpha ,\beta ]`$ and $`k=1`$. From Theorem 8, noting that $`\stackrel{~}{p}_1^{(N)}(z)=z`$ for any $`N2`$, we have

###### COROLLARY 9

For each $`j`$-type QRW, if $`N2`$, then

$`P_{A,1}^{(N)}(\phi )={\displaystyle \frac{1}{2}}\left(1+{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}|\stackrel{~}{r}_1^{(N)}(e^{i\theta })|^2𝑑\theta \right)(1+2\mathrm{}(\overline{\alpha }\beta ))`$
$`P_{G,1}^{(N)}(\phi )=|\alpha |^2+\left({\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}|\stackrel{~}{r}_1^{(N)}(e^{i\theta })|^2𝑑\theta \right)|\beta |^2`$
where $`\stackrel{~}{r}_1^{(2)}(z)=0,\stackrel{~}{r}_1^{(3)}(z)=z^3/(2z^2)`$ and for $`N4`$

$`\stackrel{~}{r}_1^{(N)}(z)={\displaystyle \frac{z^2J_{N3}(z)J_{N4}(z)}{\sqrt{2}(J_{N3}(z))^2zJ_{N3}(z)J_{N4}(z)\sqrt{2}(1)^{N3}}}`$
with

$$J_n(z)=\underset{k=0}{\overset{n}{}}\lambda _+^k\lambda _{}^{nk},\lambda _++\lambda _{}=\sqrt{2}\left(z\frac{1}{z}\right),\lambda _+\lambda _{}=1$$

In particular, when $`\phi ={}_{}{}^{t}[0,1]=|R,k=1`$ and $`N=2,\mathrm{},6`$ cases, the above corollary for the A-type QRW implies that the values $`P_{A,1}^{(N)}({}_{}{}^{t}[0,1])(N=2,\mathrm{},6)`$ satisfy the conjecture given by Eq. (4.6).

Acknowledgments. This work is partially financed by the Grant-in-Aid for Scientific Research (B) (No.12440024) of Japan Society of the Promotion of Science. I would like to thank Takao Namiki, Takahiro Soshi, Hideki Tanemura, and Makoto Katori for useful discussions.

Department of Applied Mathematics

Faculty of Engineering

Yokohama National University

Hodogaya-ku, Yokohama 240-8501, Japan

norio@mathlab.sci.ynu.ac.jp


text:
# Study of Radiation Damage in Lead Tungstate Crystals Using Intense High Energy Beams

## I Introduction

The BTeV  experiment is being readied to study beauty and charm physics at the Fermilab Tevatron collider. The goals are to make an exhaustive search for physics beyond the Standard Model (SM) and make precise measurements of the SM parameters. The important measurements to make involve CP violation, mixing, and rare decays of hadrons containing b or c quarks. Since detection of photons, mostly from $`\pi ^o`$ or $`\eta `$ decays is essential to accomplish our physics objectives, we have decided to use an electromagneticcalorimeter (EMCAL) made of lead tungstate PbWO4 (PWO) crystals. These crystals produce light proportional to the incident electromagnetic energy; this light will be sensed by photomultiplier tubes. This system is ideal for a heavy quark experiment at a hadron collider because of excellent energy and position resolution, a compact shower size that minimizes overlapping showers (due to the small Moliere radius), fast signals that minimize shower overlaps in time and expected excellent radiation hardness.

Pioneering work on PWO crystals performance was done at Protvino . These results showed the promise of such crystals. However, the technology of mass producing such crystals with high purity was not yet known. The CMS group worked with companies both in Russia and China to perfect these techniques .

In high luminosity collider experiments, PWO crystals will be irradiated by high energy particles and accumulate significant absorbed doses, up to a few Mrad. The radiation hardness of PWO crystals has been studied by the CMS group using radioactive sources and electron beams . The general conclusion is that lead tungstate crystals were radiation hard, and that the damage in crystals depends only on the dose rate . It is, however, important to measure radiation damages of PWO crystals in high energy particle environments which are more similar to that which these crystals will be exposed to.

It is important to emphasize  that in a hadron collider experiment radiation effects from hadronic interactions and neutrons could be much more serious than seen with photons or electrons of the same doses. Compared to photons or electrons, high-energy hadrons will be able to induce inelastic nuclear reactions which will locally destroy the crystal lattice. In particular, they can create nuclear fragments with very high energy transfer and lead to extended clusters of crystal lattice distortion. A simple calculation suggests that such interaction may produce significant number of additional crystal defects over the life of BTeV. Therefore it is crucial to study the radiation hardness of PWO-crystals using a hadron environment which is similar to the BTeV EMCAL expectations. Such radiation studies with lead tungstate crystals have been carried out for the first time. The results of this study are presented in this paper.

The general goal of our test beam studies was to evaluate the performance of lead tungstate crystals produced by two manufacturers in Russia, Bogoroditsk and Northern Crystal in Apatity, and one in China, Shanghai Institute of Ceramics. More specific goals were to understand how to set specifications for purchasing crystals, confirm energy and position resolution predictions, measure the radiation rate dependence of light output, and measure the correlation between light output and the LED calibration system at varying radiation loads. The 2B beam channel at the Protvino accelerator U70 has been specifically developed to provide these measurements . Results on energy and position resolutions of the PWO crystals which were obtained in these runs have been published elsewhere .

This paper is organized as follows. A general picture of radiation damage of PWO crystals as well as the results of simulations on dose rate profiles in the PWO crystals with the use of the MARS program  are described in Sec. II. These calculations are made for the BTeV experiment and for the two types of radiation studies of PWO crystals which have been carried out in Protvino for BTeV. In these studies we irradiated crystals with (a) moderate dose rates (1-60 rad/h) of high-intensity high-energy electron and pion beams in the secondary particle channel 2B and (b) super-intensive dose rates of mixed beam at a dedicated facility that was several meters away from the main ring of the U70. The test beam facility for approach (a), including phototube monitoring as well as the results of the moderate dose rates irradiation are discussed in Sec. III. Three accelerator runs, each up to a month long were devoted to these studies. The results from approach (b) are given in Sec. IV. The conclusions of the entire radiation studies are presented in Sec. V.

## II Radiation damage and absorbed dose profiles in the crystals

Radiation hardness studies of detectors and electronics are an important concern in EMCAL design  All crystal scintillators suffer from radiation damage. The most common radiation damage is due to color center formation, which results from trapping of electrons in crystal defects such as vacancies, displacements and impurities . These electrons are often in metastable states and can be excited by visible photons to higher energies. Color centers reduce light transparency of crystals, resulting in reduced light output. Additional damage may be caused by hadrons when they create crystal defects by displacing nuclei or changing nuclei to different nuclei. This kind of damage can not only reduce light transparency, but, in principle, also reduce primary scintillation light itself. It would be more difficult to monitor the latter effect. Since the trapped electrons are in metastable states of varying lifetimes and “potential barriers”, some of them may disappear very quickly, whereas others may be almost permanent.

When a PWO crystal no longer receives radiation, its color centers (semi-stable excited states) disappear, and it recovers from transmission degradation by natural room-temperature annealing. In fact, this annealing goes on even during radiation exposure. In general the rate of radiation damage decreases with the amount of damage. Therefore, when crystals are exposed to a constant radiation level, they lose light only up to the point when the rates of radiation damage and natural recovery balance. Raised temperatures accelerate the recovery process and so may ultra violet irradiation. Because the damage may recover at room temperature, it leads to a dose rate dependence of the light output.

The CMS experimental data, mainly from photon and electron irradiation, indicate that the light transmission of crystals deteriorates due to formation of color centers by radiation, while the scintillation mechanism itself seems unaffected. Besides dependence on the dose rate, the radiation damage of PWO crystals could also be sensitive to the type of radiation. In particular, the properties of crystals could be significantly degraded in hadron beams by displacement damage effects, i.e. distortions of the crystal structure. In these studies it is very useful to know the hadron fluence, the hadron spectra and the absorbed dose rate.

The BTeV EMCAL extends radially outward from the beam line. The crystals near the beam pipe receive the maximum dose. In order to ascertain the level of radiation in the crystals we performed calculations using the MARS code. Results are given in Table I.

We tried to emulate the BTeV conditions as much as possible. A 27 GeV electron beam and a 40 GeV $`\pi ^{}`$ beam have been used to irradiate the crystals with moderate dose rates. The beams were directed into the secondary beam channel from the accelerator, where primary 70-GeV protons interacted with an internal target. The MARS calculations of the absorbed dose rates in the crystals from the secondary beam channel are compared with the absorbed dose rates expected in BTeV in Fig. 1. The $`\eta `$ (pseudo-rapidity) shown here reflects the coverage of the BTeV EMCAL, where $`\eta `$ of 4.45 is at the extreme inside near the beam and $`\eta `$ of 2.27 is on the extreme outside. Electron and pion dose profiles in the crystals are different. The crystals receive damage from pions almost uniformly along their length starting from a distance of 5-7 cm from the front. For electrons an absorbed dose rate at shower maximum is two orders of magnitude higher than near the crystal ends. Because the BTeV dipole magnet sweeps particles vertically, the radiation profile at the calorimeter is different in the horizontal and vertical planes. Thus in BTeV the mix of charged hadrons and photons changes and the ratio between shower maximum and the crystal ends is only a few times in the vertical plane and an order of magnitude in the horizontal plane. That is why both electron and pion beams are used to study radiation damage of the crystals.

Two crystals, one manufactured in Bogoroditsk and the other in Shanghai were placed near the vacuum pipe of the Protvino U-70 accelerator in the first dedicated super-intensive dose rate study. These crystals were irradiated by secondary particles coming out the internal target of the accelerator. The energy spectra of neutrons, gamma-quanta and charged hadrons at the place where the crystals were irradiated are shown in Fig. 2(b). For comparison the expected particle spectra at the front face of the BTeV EMCAL are presented at the top part of the same Figure. We can see that the spectra look similar, although the dose rate in the IHEP irradiation zone is about three orders of magnitude higher than expected in BTeV. In the second dedicated intensive study, four more crystals from Bogoroditsk and Shanghai were exposed to radiation at the same facility. The intensity of the second run was reduced by two orders of magnitude. Absorbed dose rates as a function of longitudinal position for these two exposures are presented in Fig.3.

## III Moderate dose rate irradiation

In this Section, we describe the testbeam facility for moderate dose rate irradiation studies, discuss phototube gain monitoring, and present the results of irradiating crystals with electrons and pions. In our radiation studies we wanted to use radiation conditions as close to the BTeV conditions for the crystals as possible. Absorbed dose rates as a function of longitudinal profile at the BTeV EMCAL and at IHEP testbeam have been already discussed in Section II and presented in Fig. 1. We used 27-GeV electrons and 40-GeV pions to irradiate crystals in the three accelerator runs.

### A Test beam facility

The test beam setup consisted of 5x5 PWO crystal array situated inside a temperature controlled light-tight box (ECAL), a beam with a momentum tagging system and a scintillation counter trigger system .

All the crystals we used were rectangular in shape. The Bogoroditsk and Shanghai crystals were $`27\times 27`$ mm2 in cross section and 220 mm in length. The Apatity crystals were $`22\times 22`$ mm2 in cross section and 180 mm in length. Light from each crystal was collected by a 10-stage 1-inch diameter Hamamatsu R5800 photomultiplier tube (PMT). All the crystals were wrapped by a 170 $`\mu `$m thick tyvek. A radioactive source study at University of Minnesota showed that tyvek is radiation hard up to at least a few Mrad. This study as well as the Belarussian State University(Minsk) one also showed that a borosilicate glass did not lose any light at least up to 10 krad, a quartz glass up to 1 Mrad, both with an accuracy of 1$`\%`$. Six quartz PMT’s were used for a part of our test beam study, the rest were the borosilicate PMT’s.

We accumulated absorbed doses in our crystals up to a few krad. No changes inside the box or PMT HV values were made during the irradiation period. The PWO light yield strongly depends on crystal temperature . The 25 crystals were surrounded by a set of four copper plates that were water cooled, which enabled a temperature control using a Lauda cryothermostat. The temperature for the study described in this paper was fixed at $`20^0C`$ $`\pm `$ $`0.1^0C`$. To measure the temperature of the crystals, 24 temperature sensors were mounted on the front and rear faces of the crystals.

For the most of the results presented in this paper, the crystal array was monitored with the four different wave length light emitting diodes (LED). The LEDs had the following wavelengths : 660 nm (red), 580 nm (yellow), 530 nm (green), and 470 nm (blue). Transmission of red light in the crystals is not affected much by radiation damage , so the red LED monitors the PMT gain change. One LED generator with a multiplexer was placed into the light-tight box with the crystals and used for the all the LEDs. The LED temperature dependence is on average 1$`\%`$/$`1^{}`$C, and thus limited to 0.1% because of our careful temperature control inside the box. We had one bunch of fibers between the LED generator and the crystals. In each accelerator cycle 10 pulses data for one LED color were collected. Four cycles were needed to collect all the LED signals.

An $`\alpha `$-source (YAP-light pulser ) was mounted on the photocathode of a separate PMT in addition to the fiber to monitor LEDs themselves. It had 20 decays/sec with about 5,000 photons/pulse. Forty pulses were collected each spill. A signal from last dynode of this PMT was used to form an $`\alpha `$-trigger. The size of the YAP crystal was 3$`\times `$3 $`mm^2`$ with thickness of 0.1 mm. Its emission spectrum has the maximum at 360 nm. The YAP crystal temperature dependence of the light output was 0.4$`\%`$/$`1^{}`$C. The $`\alpha `$-spectrum as well as $`\alpha `$-stability is presented in Fig. 4. One can see that this stability over 85 hours is better than 0.2$`\%`$. A Hamamatsu PIN diode S6468-05 with integrated amplifiers was also used to monitor the LEDs because it has a good sensitivity in the red region as well as a gain stability. It’s temperature dependence is much less than 0.4$`\%`$/$`1^{}`$C.

We did not use an optical grease coupling between the crystals and the PMT’s in order to avoid a contribution of a possible radiation damage of the grease. The PMT’s were attached to the crystals without any optical material between them. High voltage to the tubes was supplied by a LeCroy 1440 HV system. Signals were sent to the control room patch-panel without any connection to ground inside the crystal box to avoid ground loops. A LeCroy 2285 15-bit integrating ADC was used to measure charge over 150 ns without pedestal subtraction. The ADC sensitivity was 30 fC per count. At HV values around 1000 V in the tubes we had about 2 MeV/ADC count.

### B Phototube gain change monitoring

We used high-intensity high-energy electron beam to irradiate the crystals and at the same time monitor the light output. The beam particles travel along the length of the crystals toward the PMT. We needed to take into account the possibile phototubes gain changes, for example, from varying in the beam intensity. Thus, we carried out two types of PMT gain change studies to separate the effect of PMT gain change from crystal radiation damage. We investigated the possible changes in PMT gains at a dedicated stand at IHEP after the accelerator runs. We also monitored the PMT’s continuously during one of the runs using the red LED.

Fig. 5 shows a schematics of the dedicated stand setup to study the PMT behavior, where the average anode current was adjustable by changing the intensity of DC light shining on the PMT. The setup consisted of a high quality referenced PMT(Hamamatsu R5900), a blue LED light pulser, a DC LED. Both pulsed and DC LED lights were injected into the test PMT through optical fibers. The stability of the pulsed LED itself was monitored by a $`Pu`$ radioactive source implanted in a crystal and mounted at the photocathode of the reference PMT. The read-out and control electronics were placed in a CAMAC crate which had an interface with a PC. The average anode current was chosen for each test PMT to be the same as what we had at the test beam. The anode current was measured directly by an ammeter. Fig. 6 shows a timing diagram of various measurements. Each set of measurements took 2 minutes. At the beginning of each set, we measured the pulse heights of two groups of 2000 light pulses. It took 20 sec to collect 2000 pulse data and there was a 10-sec interval between the two groups. The data from a radioactive source in a self trigger mode were collected during the remaining 70 seconds. This 2-minute set was then continuously repeated. The intensity of the DC LED to induce a finite average anode current in the tested PMT was allowed to change, if needed, in the 10 second time intervals. This system allowed us to make PMT long-term stability measurements with a precision of 0.2$`\%`$.

We used positive HV for PMT’s with grounded photocathodes for the first accelerator run and negative HV for the second and the third runs. The red LED response of PMT number 743 using negative HV during irradiation study is presented in Fig. 7. The behavior of the same PMT at the dedicated stand is shown in Fig. 8. We see that short-term loss of the signal is 3-5 $`\%`$ when the test beam intensity is at the level of $`10^4`$ e-/sec averaged over the entire accelerator cycle.

A similar signal loss was seen (Fig. 8) when the additional green LED was turned on to produce the anode current of $`5\mu `$A. Another similarity is that when the “beam” or green LED was turned off, the PMT gain rose by a few percent. We compared the behavior of each PMT at the stand and during the beam test and found a satisfactory agreement between these results.

Six phototubes with the quartz glasses (Hamamatsu R5800Q) of the same size were used to avoid the possible radiation damage to phototube windows. These PMT’s had a gain change of 5-6$`\%`$ and one of them even 10$`\%`$. The gain variation of one of these quartz phototubes is presented in Fig. 7. For the PMT’s with borosilicate glasses the signal loss has been measured not to exceed 3$`\%`$ for both positive and negative HV.

The blue LED signal amplitude over 85 hours is presented in Fig. 9. Fig. 9(a) shows the electron beam intensity over this time period in a sample crystal. Fig 9(b) shows the raw blue LED signal for the same crystal. We see that the blue LED signal fell by 5-6$`\%`$ when the beam was off. The time diagram of the blue LED corrected by the red LED is shown in Fig. 9(c). Note that most of our PMT’s lost gain when the beam was on. In our plot we selected this PMT with the opposite behaviour to show that we could correct for this big gain change even though the sign of the change was atypical.

All of our analyses included corrections using the red LED data. We corrected the signals from electrons and blue LED on the signal from red LED to subtract a PMT gain variation effect from the total signal for each PMT. When the green or yellow LED signals were used in the analysis, they also were corrected using the signal from red LED.

To check and correct the stability of the red LED, we used the $`\alpha `$-source. The instability of blue, green and yellow LEDs was corrected using the PIN diode. The ratio of the PIN to $`\alpha `$ signals was stable to an accuracy of 0.1 $`\%`$. To decrease any possible remaining LED instability left after these corrections, we kept for further analysis only accelerator spills with similar beam intensity. We conservatively estimate that the error for the blue LED signal is 0.2$`\%`$.

### C Irradiation by high-energy electrons

The crystal array was irradiated by 27 GeV electrons for one week with an accelerator efficiency of 85$`\%`$. The beam intensity at the crystal array was 6$`\times `$105 particles/spill most of the time during this period. The 80$`\%`$ of the beam entered in one of the six central crystals. About a half of the time, the beam was centered on one crystal in the array and during the rest of the time it was centered on another crystal. Coordinates of the electrons entering the crystal array were measured by the drift chambers. The events with electrons near the center of the crystals were selected for data analysis.

We now describe the analysis of the electron beam data. All the information which was accumulated during 85 beam hours (one position of the beam at the array, see above) was divided into pieces of two hours long each. This choice was made to have enough statistics to measure the average energy deposit in a crystal with an accuracy of 0.3$`\%`$, and thus we could continuously monitor the crystal signal loss. Prior to the irradiation study, the PMT gain of the each crystal in the array was adjusted to 10,000 ADC counts when 27 GeV electron hit the center of the crystal. Since this corresponds to 76$`\%`$ of the full electron energy , one ADC count corresponded to 2 MeV. The size of the beam spot was chosen 4x4 $`mm^2`$ for most irradiated crystals and 6x6 $`mm^2`$ for crystals with lower doses in order to equalize the statistics. The true coordinates of a particle at the array was calculated with the information from the last drift chamber which was close to the array. The accumulated energy peaks were fitted by a Gaussian. Then the mean values were corrected using the red LED.

For each crystal a dose rate was defined as an effective number N of electrons per second hitting this crystal multiplied by $`25.910^4`$ (see Fig. 1(d) ). The number N was calculated as energy deposit in this crystal in GeV/sec divided by 20.5 GeV (it corresponds to 76$`\%`$ of 27 GeV energy deposit when electron hits the center of the crystal  in accordance with the MARS simulation).

A typical result for an irradiated crystal is presented in Fig. 10. In Fig. 10(b) we see an intensity of the electron beam which is shown in dose rate units at the shower maximum according to the MARS simulation results presented in Fig. 1. The absorbed dose is given in Fig. 10(c). The main result is shown in Fig. 10(a) which is the normalized electron signal. We see that finally the crystal lost 12$`\%`$ of the signal under an electron beam irradiation mostly with 15 rad/h dose rate after it accumulated 1.2 krad absorbed dose. It also appears that the radiation damage is saturating. For dose rates of 10-25 rad/h under 27 GeV electron beam irradiation, eight crystals lost an average of 8$`\%`$ after a total accumulated dose of 1-2 krad.

In order to use the light monitoring system to track the effects of radiation damage, it is necessary to determine the relation between the change observed by the monitoring system and the change in the signal from beam electrons. Because of the different optical paths taken by the injected monitoring light as compared to the scintillation light this constant is not expected to be unity. Furthermore the LED system monitors the transparency of the crystal at a specific wavelengths and thus does not sample the entire spectrum of scintillation light.

The blue LED emits at 470 nm and the scintillation peak is at 430 nm. The typical blue LED and electron signal behavior under irradiation for one of the crystals is shown in Fig. LABEL:fig:apr02-3a(a). The blue LED (as well as the electron signal) is corrected by the red LED, and the red LED by the $`\alpha `$-source. The same was provided for the green and yellow LEDs. For the green LED a signal loss was smaller than for the blue LED, and for the yellow LED the signal loss was smaller still (not shown). In Fig. LABEL:fig:apr02-3a(b) we see a strong correlation between the change in the blue LED light level and the beam signal. We fit such distributions by the straight lines, ignoring some deviations from linearity. The results for a few crystals are presented in Fig. LABEL:fig:apr02-3a(c). We did not observe a significant difference in the crystals from different manufacturers. Constants of proportionality vary from 0.3 to 0.6 for these crystals. The dependence of a relative electron signal on the absorbed dose is presented in Fig. 12.

A simple model is used to describe signal loss. The signal loss $`dy`$ is proportional to the signal value $`y`$ and the number of the produced color centers, which are proportional to the absorbed dose $`dR`$. Crystal recovery is proportional to a difference between the asymptotic value $`y_0`$(after recovery) and the current signal value. Also it is proportional to the recovery time $`dt`$:

$$dy=P1ydR+P2(y_0y)dt=((P1\frac{dR}{dt}+P2)y+P2y_0)dt$$
(1)
In our case the dose rate($`\frac{dR}{dt}`$) was almost the same during the 85 hours of irradiation. Integration of this equation gives us the expression:

$$y=P0\mathrm{exp}^{(P1\frac{dR}{dt}+P2)t}+\frac{P2y_0}{P1\frac{dR}{dt}+P2}$$
(2)
We can present the signal loss behaviour function as

$$f(t)=a\mathrm{exp}^{t/\tau }+(1a),$$
(3)
The results of the fit for Fig. LABEL:fig:apr02-3a(a) are listed in Table II.

The parameter $`a`$ defines the saturated light loss value that is reached as $`t`$ goes to infinity at a constant dose rate. Close to the asymptotic value, the crystal lost 10$`\%`$ in the electron signal and 5$`\%`$ in the blue LED signal. The $`\tau `$ parameter defines the saturation time constant, which is 30 hours for our crystal and our dose rate.

The time constants for the ten studied crystals are between 20 and 30 hours. There is no significant difference in $`\tau `$ for the LED and electron signals.

We should make a note at the end of this section. When a crystal is irradiated, the red LED light is slightly absorbed. Herewith, the blue LED light is absorbed more, in 3-6 times more  compare to the red LED light in the crystals. We can estimate that the electron signal is absorbed in about two times more than the blue LED signal (see Fig. LABEL:fig:apr02-3a(c)). We assumed that red light was unchanged under crystal irradiation, and assigned the PMT gain change to the red LED change. It means that the absolute electron signal loss values might be in about 1.1 times higher than the presented ones.

### D Irradiation by high-energy pions

After the electron irradiation program was finished, we irradiated the same crystals with pions for a four day period. We used a 40 GeV $`\pi ^{}`$ beam. The size of the 40 GeV pion beam was 8 cm horizontally and 6 cm vertically, $`i.e.`$ 90% of the beam was contained within these dimensions. The beam intensity was 6 $`\times 10^6`$/sec. Six crystals were irradiated with a dose rate ranging from 10 to 30 rad/h. Five cycles of irradiation (15-20 hours each) were alternated by low intensity electron beam exposures to measure the scintillation signals in the crystals.

The radiation damage region in the crystals is different for an electron and a pion irradiation (see Fig. 1). Thus, if a crystal was irradiated first by electrons until saturation in radiation damage was reached for a given dose rate, then we expect to get an additional signal loss with pion irradiation even at the same dose rate. Fig. 13 shows the additional loss of signal for one of the crystals (from Apatity). This crystal was irradiated by 27 GeV electrons and then by 40 GeV pions. During the 85 hours of e- irradiation the dose rate was 12 rad/h. Then in the next 85 hours the dose rate was an order of magnitude less, and the crystal recovered. As a result, the first filled square point for pion irradiation data is above many open points for electrons. After that the crystal was irradiated by pions with the dose rate 12 rad/h for 100 hours. We see that the crystal lost 8$`\%`$ of the signal during the electron irradiation period and 14$`\%`$ of the signal during the pion irradiation period with the same dose rate. The constant of proportionality between the blue LED and the electron signal is 0.3 for electron irradiation (if one fits by a straight line). This is about the same for pions at the beginning of the pion irradiation, but then increases up to 1 during the further pion irradiation. The crystals (manufactured in Bogoroditsk and Apatity) lost about 14$`\%`$ on pion irradiation. The surrounding crystals which were irradiated with a dose rate of about 1 rad/h lost less than 1 $`\%`$ of their light output.

The dependence of a signal loss on dose rate was studied in a separate run using 40 GeV pion irradiation. Each beam exposure lasted for 6 continuous hours. The beam intensity started from 2$`\times 10^5`$/sec and was increased in a few steps up to 8$`\times 10^6`$/sec by the end of the study. The beam was present in 1 sec of the full accelerator cycle of 9 sec. After each 6 hour irradiation exposure we lowered intensity by a few orders of magnitude, down to 3$`\times 10^4`$/sec, so that we could avoid pile-up and see a minimum ionizing peak (MIP) for pions traversing the crystals without interacting. The crystals light output signals were monitored using the MIP peak; this procedure took 2 hours at low intensity. After that we took again high intensity beam exposure for the next 6 hours to continue irradiating the crystals. Then again switched to the low intensity MIP exposure.

To check our procedure for obtaining the change in scintillation light from time to time, we used pure muon beams and 27 GeV electrons to measure the light output changes due to pion irradiation. We continued this procedure of alternating high intensity and low intensity beams for 10 days in a row. The dependence of the normalized MIP signal on an absorbed dose for the two crystals in the array is shown in Figs. LABEL:fig:mip\_22(a) and (b). (The normalized MIP signal is defined as the ratio of the MIP signal after some absorbed dose to the one before the pion irradiation).

We have observed the dependence of light output loss on the dose rate. Like electron radiation, the light loss exhibits saturation effect when the dose was kept at a constant level. The correlation between a change in the LED signal and a change in the MIP signal under irradiation was also measured (see Fig. LABEL:fig:mip\_22(c) as an example). The constant of proportionality, if one fits by a straight line is different for different crystals and is on average 0.7 (the LED signal decreases less than the MIP one). In Fig. LABEL:fig:mip\_led\_23(a) we show the decrease in the LED signal for moderate dose rates. Different crystals received different absorbed doses during 10-day irradiation period. The open circles stand for Bogoroditsk crystals and the filled circles stand for Shanghai crystals. Six points in Fig. LABEL:fig:mip\_led\_23(a) represent the six crystals described above which accumulated absorbed doses of more than 1 krad each. Fifteen other crystals were irradiated by the beam halo and received absorbed doses less than 500 rad each. They are shown on the left side of Fig. LABEL:fig:mip\_led\_23(a) . Three Bogoroditsk crystals are close to each other in their radiation hardness, however the twelve Shanghai crystals differ among each other by an order of magnitude.

Irradiation of lead tungstate crystals creates color centers which reduce the light attenuation length. One expects that the change of attenuation length will affect the longitudinal uniformity. This can degrade the energy resolution. On the other hand, if the loss of light collected in the crystal after irradiation is relatively small, the energy resolution itself might not be degraded so that the radiation damage can be regarded as only a calibration issue. The non-uniformity of the light yield (LY) along the crystal contributes to the energy resolution. To measure changes in the LY non-uniformity the crystal array was rotated by 90 0 with respect to the beam direction, before the irradiation by pions and just after the 10 days of irradiation. The crystals were scanned with the muon beam. The position of the muon track going through the crystal was reconstructed with the drift chambers. The data were binned along the crystal lengths in 1 cm intervals. The energy deposit distribution was fitted in each bin by a convolution of Gaussian and Landau function. The non-uniformity of the light yield in the front part of crystal (3-10 radiation lengths) was about 0.5 $`\%`$/cm. The non-uniformity did not change significantly after a dose up to 4 krad at a dose rate of up to 60 rad/h, which caused the signal loss of up to 30$`\%`$. As a result, the energy resolution of the crystals did not change. The relation between the change in transparency seen by the LED light and the change seen by the scintillation light varies from crystal to crystal. A plot of such constants of proportionality for seven crystals is shown in the Fig. LABEL:fig:mip\_led\_23(b). The first four points show the Shanghai crystals, and the next three points show the Bogoroditsk crystals. Points 4 and 5 represent the super-intensive dose rates obtained by the Shanghai crystal S25 and the Bogoroditsk crystal B21 (details will be given in the next Section). We can see that the constants of proportionality for 40 GeV pion irradiation (Fig. LABEL:fig:mip\_led\_23(b) ) are larger than the constants of proportionality for 27 GeV electron irradiation (Fig. LABEL:fig:apr02-3a(c) ).

After the irradiation by electrons and then pions was finished, we kept the PMT HV on and studied crystal recovery for 15 days using the LED pulser. The results for the Apatity 1434 crystal are presented in Fig. 16. We fitted the dependences of recovery on time for the six crystals with an exponential function. The average recovery time is (200$`\pm `$40) hours, and the LED damage recovery for 400 hours is (87$`\pm `$5)$`\%`$ for these six crystals.

## IV Super-intensive beam irradiation

Six crystals from Bogoroditsk and Shanghai were irradiated by secondary particles coming out the internal target of the 27-th magnet block of the Protvino U-70 accelerator(see Fig.17). Two of them were irradiated at a dose rate of 100 krad/h, and the other four at 1 krad/h. For the latter case the intensity of the primary proton beam was lowered by two orders of magnitude. To measure the absorbed dose, thermo-luminescence dosimeters (TLD) were attached to the front face of the crystals. They were of LiF type doped by Mg, Cu and P, 5 mm in diameter and 200 $`\mu `$m in thickness. In addition, an ionization chamber(IC) filled by Xenon was installed behind the crystals. The sensitive volume of the chamber was as 18.5 mm in diameter and 36 mm in length. Both TLDs and the IC were calibrated using a Cs-137 gamma source. The accuracy of the absorbed dose measurements by TLDs and IC in this mixed radiation field was estimated to be $`30\%`$ each. These measurements were in general agreement with the results of the MARS calculations; the worst case difference was a factor of 1.5. The dominant systematic error of the calculations was due to the accuracy of the irradiation facility geometry. The IC was used to monitor the number of protons produced at the internal target for each run of the crystal irradiation. $`Al`$ activation detector in Fig.17 was used to measure a fluence of the hadrons (number of hadrons per cm2) with energy greater than 20 MeV.

Two crystals, Bogoroditsk B21 and Shanghai S25, were irradiated in the first exposure at about 100 krad/h dose rate. The longitudinal profiles of the absorbed dose rates are shown in Fig. 3(a). The maximal values of the absorbed doses accumulated in crystals during the five exposures are given in Table III. Four crystals, Bogoroditsk B17, B9 and Shanghai S22, S18, were irradiated in the second exposure at 1 krad/h dose rate. The longitudinal profiles of the absorbed dose rates are shown in Fig. 3(b). The absorbed doses accumulated in the crystals during the five exposures are given in Table IV.

The results of the irradiation of the two crystals in the first exposure are presented in Fig. LABEL:fig:bogoroditsk. The procedure was to irradiate the crystals and then measure their light output immediately thereafter using the 27 GeV electron beam. In some cases, we measured the light output again after letting the crystals sit without any radiation. The Bogoroditsk crystal (see Fig. LABEL:fig:bogoroditsk(a) ) lost $`33\%`$ of the initial signal after first 3.4 krad dose. After the second irradiation, the absorbed dose increased up to 43 krad and the signal loss was increased up to $`46\%`$. After 47 hours of recovery time, the signal rose up to $`70\%`$. After the third irradiation, the Bogoroditsk crystal accumulated 313 krad and the signal was at the level of 49$`\%`$. 32 hours of recovery time returned it to a level of 57$`\%`$. After the fourth dose the total radiation was 2300 krad and the signal level was at 37$`\%`$. After 15 hours recovery time, the signal recovered slightly to 39$`\%`$.

The Shanghai crystal (see Fig. LABEL:fig:bogoroditsk(b) ) lost $`18\%`$ of the signal after first 1.8 krad dose. After the second irradiation, the integrated dose increased up to 23 krad and the signal loss increased up to $`33\%`$ relative to the signal before the irradiation. After 47 hours of recovery time, the signal rised up to $`69\%`$. After the third irradiation, the Shanghai crystal accumulated 163 krad and the signal was at the level of 66$`\%`$. One should mention that the signal was pretty stable between the second and the third irradiations including recovery time and was at the level 66-69$`\%`$ for the absorbed doses of 23-163 krad. After getting 2800 krad, the signal dropped down to 33$`\%`$.

One of the most important conclusions of this work is that even after an integrated dose about 2.5 Mrad obtained with a super-intensive dose rate 100 krad/h both crystals remained usable, although they lost 2/3 of their light. In BTeV we expect that only $`0.1\%`$ of the crystals will receive this much dose in a year. As was expected, the constants of proportionality in the MIP-Electron correlations for both the crystals are about 1. The LED-Electron correlations for both the crystals are shown in Fig. LABEL:fig:e\_led\_rus. The constants of proportionality for the two crystals are 0.5 and 0.66. The degradation of single crystal energy resolution for 27 GeV electrons was only $`20\%`$ for Bogoroditsk crystal and $`50\%`$ for the Shanghai crystal.

The four crystals irradiated in the second exposure with a dose rate of 0.5-1 krad/h, and a total dose of 350-700 rad, lost up to 10$`\%`$ of their light output for Shanghai crystals and up to 25$`\%`$ for Bogoroditsk crystals. After each of the next three runs no signal loss was seen, within the 3$`\%`$ accuracy (the systematic error due to a PMT gain change effect).

## V Summary and Conclusions

Radiation hardness of lead tungstate crystals is an important issue for the BTeV experiment at Fermilab. Simulation of absorbed dose profiles in the crystals with the use of the MARS program has shown that the dose rates for the crystals range from 0.1 up to 700 rad/h. About 95$`\%`$ of the crystals in the BTeV electromagnetic calorimeter will get the absorbed doses from 0.1 to 30 rad/h assuming that the Tevatron luminosity is $`2\times 10^{32}cm^2sec^1`$. Almost 5$`\%`$ of the crystals will get from 30 rad/h up to 200 rad/h, and only 0.5$`\%`$ more than 200 rad/h.

A study of radiation damage in lead tungstate crystals has been carried out in Protvino in 2001-2002 for BTeV. The crystals were manufactured in Bogoroditsk (Russia) and Shanghai (China) at the very end of 2000, and in Apatity (Russia) in early 2002. There were two approaches in the study. First, crystals were irradiated by high-intensity high-energy electron and hadron beams at radiation doses ranging from 0.1 to 60 rad/h. Secondly, crystals were irradiated by charged hadrons, $`\gamma `$-quanta and neutrons from the internal target of the U70 in a wide energy spectrum from 10 eV up to 70 GeV at dose rates between 0.5 and 100 krad/h.

The dependence of light output loss on a dose rate has been measured. The light loss exhibited saturation when the dose rate was kept constant. At larger dose rates, the light output loss still saturates but at lower light output levels. Each crystal had a different percentage of light loss when it saturated. More quantitatively: no light output loss was observed for dose rates less than 1 rad/h. For dose rates of 10-25 rad/h with 27-GeV electron irradiation, eight crystals lost on the average 8$`\%`$. For 40 GeV pions this average was 12$`\%`$ at comparable irradiation dose rates. The difference between the damage due to electron and pion irradiation can be attributed entirely to their difference in the radiation profile along the length of the crystal. Much of electron energy is deposited near the shower maximum, from 4 to 10 cm from the front of the crystal. For pion beams, the radiation dose profile reaches its maximum around 5-7 cm and stretches all the way to the rear-end of the crystal. However, a possible effect due to the difference between the physical processes by which electrons and pions interact with crystals cannot be ruled out.

For dose rates of 30-60 rad/h using 40-GeV pion irradiation, five crystals lost on the average 20$`\%`$. For a dose rate of 500 rad/h using irradiation by charged hadrons, $`\gamma `$-quanta and neutrons with the average energy of 10 GeV, two crystals lost 10$`\%`$, and two other crystals lost 25$`\%`$ when they were exposed to 1 krad/h of radiation. Two crystals got extremely high dose rate of 100 krad/h and accumulated about 2.5 Mrad absorbed dose (maximum annual dose of any BTeV crystals!) also with the same mixed particle spectra irradiation. They remained useable. Their light output loss was a factor of 3. This is far from the BTeV environment, where 700 rad/h will be the highest $`0.1\%`$ crystals.

There is a correlation between a change in the LED signal and a change in the beam (electron or MIP) signal under irradiation. The constant of proportionality is different for different crystals and varies from 0.3 to 0.6 for electron irradiation and from 0.5 to 0.9 for pion irradiation.

The non-uniformity (maximum 0.5$`\%`$ per cm at one third of the crystal length) of the light yield does not change significantly when the dose rate is up to 60 rad/h. After 2.5 Mrad absorbed dose with a dose rate of 100 krad/h the uniformity became 1.5 times poorer, at least for one of our crystals.

When irradiation decreases or stops, crystals recover. The average recovery time for six crystals which lost from 7 to 20$`\%`$ of the LED signal, is (200$`\pm `$40) hours, and the damage recovery after 400 hours was (87$`\pm `$5)$`\%`$.

To summarize, lead tungstate crystals lose light from irradiation by high-intensity high-energy beams. This loss level depends on dose rate. If dose rate does not change, the light loss saturates. If the dose rate is reduced, the light output recovers. Crystals have to be calibrated continuously during the BTeV experiment. We did not see a significant difference in radiation hardness of the crystals from the three manufacturers.

## VI Acknowledgments

We would like to thank the IHEP management for providing us a beam line and accelerator time for our testbeam studies. Special thanks are due to Fermilab for providing the equipment for data acquistion. We would like to thank A.P.Bugorsky, O.A.Grachov, I.V.Kotov, V.P.Kubarovsky and R.Y.Zhu for useful discussions. This work was partially supported by the U.S. National Science Foundation and the Department of Energy.

## VII References