text:
# The Quantum-Classical Correspondence in Polygonal Billiards

## I Introduction

The relation between wave motion in the short-wavelength regime and the corresponding ray dynamics is of fundamental importance in quantum mechanics, electromagnetics and acoustics. It has been found in quantum mechanics that the properties of the stationary solutions of the Schrödinger equation, the eigenstates $`|E_j`$ of the Hamiltonian, reflect the degree of order in the classical ray dynamics. For example, energy eigenfunctions $`x,y|E_j`$ of classically integrable systems with quasiperiodic motion on invariant tori are regular, while eigenfunctions of classically chaotic systems with ergodic motion on energy surfaces are typically irregular. For generic systems, it has been conjectured \[Berry77\] (see \[Robnik98\] for a review) that in the (semi-)classical limit $`\mathrm{}0`$ the averaged Wigner transforms of typical energy eigenfunctions “condense” uniformly onto the underlying classical stationary objects in phase space, which are usually invariant tori, chaotic components or entire energy surfaces. In this article, however, we consider a class of systems with exotic classical invariant surfaces for which it is not clear a priori whether such a condensation scenario exists.

The classical free motion inside a planar domain $`𝒬`$ with elastic reflection at the boundary has a constant of motion, Hamilton’s function $`H=p_x^2+p_y^2+V(x,y)`$. The particle’s mass is $`1/2`$, and $`V(x,y)=0`$ if $`(x,y)𝒬`$ and $`\mathrm{}`$ otherwise. A billiard with polygonal boundary has a second constant of motion $`K(p_x,p_y)`$ if all angles $`\alpha _j`$ between sides are rationally related to $`\pi `$, i.e. $`\alpha _j=m_j\pi /n_j`$, where $`m_j,n_j>0`$ are relatively prime integers \[Hobson75, ZemlyakovKatok76\]. However, this does not imply integrability in the sense of Liouville-Arnol\[Arnold78\]. The Poisson bracket $`\{H,K\}`$ vanishes identically only if $`m_j=1`$ for all $`j=1,2,\mathrm{}`$, so only rectangles, equilateral triangles, $`\pi /2,\pi /4,\pi /4`$-triangles, and $`\pi /2,\pi /3,\pi /6`$-triangles are integrable. Critical corners with $`m_j>1`$ destroy the integrability in a singular way: $`\{H,K\}`$ is zero everywhere in phase space except at a measure-zero set corresponding to the critical corners; the phase space is foliated by two-dimensional invariant surfaces $`H,K=\text{constant}`$ as in integrable systems, but they do not have the topology of tori. The motion is not quasiperiodic and characterized as pseudointegrable \[RichensBerry81\].

The quantum mechanical free wave motion with Dirichlet boundary conditions on a rational polygon barely reflects the classical pseudointegrability if only individual eigenfunctions of the Hamiltonian $`𝑯=H(𝒙,𝒚,𝒑_x,𝒑_y)`$ are considered. A typical eigenfunction is hard to distinguish from eigenfunctions in classically chaotic systems \[BiswasJain90, BellomoUzer95, VUF95, TomiyaYoshinaga96\]. It is in general not an eigenfunction of the second operator $`𝑲=K(𝒑_x,𝒑_y)`$ and has nonzero uncertainty $`\mathrm{\Delta }K=(𝑲^2𝑲^2)^{1/2}`$, since the commutator $`[𝑯,𝑲]`$ does not vanish. Only a subtle signature of pseudointegrability seems to be contained in individual energy eigenfunctions, namely in the distribution of zeroes of the associated Husimi function \[BS99\]. The closeness to the quantum mechanics of chaotic systems on the one side and to the classical dynamics of integrable systems on the other side is a first hint of an unusual quantum-classical correspondence.

A further indication pointing in this direction comes from an analogy with the metal-insulator transition of the Anderson model in three dimensions. Both kinds of systems have energy level statistics close to the semi-Poisson distribution \[BGS99\]. In the Anderson model, at the transition point between extended states with Wigner statistics and localized states with Poisson distributed energy levels, this intermediate statistics has its origin in the multi-fractal character of the wave functions, which are neither extended nor localized \[SG91\]. Analogously, a typical energy state in a rational polygon is expected to be a fractal in momentum space \[BGS99\], localized around an energy surface but neither extended nor localized within this surface; condensation onto a lower-dimensional invariant surface is ruled out.

As a final indication we will present a numerical study on a particular system, the barrier billiard. Instead of Wigner transforms in four-dimensional phase space, we study distributions in the two-dimensional $`(H,K)`$-space, where each point represents a classical invariant surface. In this context, we define condensation as $`\mathrm{\Delta }H/𝑯0`$ and $`\mathrm{\Delta }K/𝑲0`$ in the semiclassical limit. This ensures that mean values of the quantum operators in highly excited states can be interpreted as well-defined values of the classical constants of motion. Our central issue is to demonstrate that even though the energy states presumably do not show such a condensation there is an alternative basis of states which do so up to classically long times.

The paper is organized as follows. In Sec.II we construct the basis for general rational polygons. Its time evolution is discussed in Sec. III. In Sec. IV we illustrate our statements for the barrier billiard. We conclude with a brief summary in Sec. V and an Appendix with details on the numerical computations.

## II Construction of the basis

In the following, we will introduce states with relatively small $`\mathrm{\Delta }K`$ at the expense of a nonzero but also small energy uncertainty $`\mathrm{\Delta }H`$. As for coherent states in the harmonic oscillator \[Glauber63\], we will minimize the product of the uncertainties involved. To do so, we consider the uncertainty relation

$$\mathrm{\Delta }H\mathrm{\Delta }K\frac{1}{2}|[𝑯,𝑲]|.$$
(1)
It is an easy matter to show that the equality in (1) holds in a state $`|L`$ satisfying

$$(𝑯𝑯)|L=ia(𝑲𝑲)|L,$$
(2)
where $`i^2=1`$, $`a`$ is an arbitrary real number and $`\mathrm{}=L|\mathrm{}|L`$ with $`L|L=1`$. $`|L`$ has to be a right eigenvector of the operator

$$𝑳=𝑯ia𝑲.$$
(3)
$`𝑳`$ does not commute with its adjoint. It therefore does not belong to the class of normal operators, which contains Hermitian and unitary operators as special cases. In general, the right eigenvectors of a nonnormal operator cannot be used as a basis. Instead, right and left eigenvectors together form a biorthogonal basis; see, e.g., \[Wilkinson65\]. In our case, however, $`E_i|𝑳|E_j`$ is a complex-symmetric matrix due to the fact that the Hermitian matrices $`E_i|𝑯|E_j`$ and $`K_{ij}=E_i|𝑲|E_j`$ can be made real. Left and right eigenvectors are therefore identical and form separately a nonorthogonal basis. Furthermore, it can be shown that $`\mathrm{\Delta }H=|a|\mathrm{\Delta }K`$. In order to have equal uncertainties we choose $`a=1`$ (states with $`a=1`$ are identical). The real and imaginary parts of the complex eigenvalues $`L`$ have a physical interpretation as mean energies $`𝑯`$ and $`𝑲`$, respectively. Note that the classical function $`HiK`$ is a constant of motion whose constant-level surfaces are the invariant surfaces.

We derive now the properties of the $`|L`$-states in the limit $`\mathrm{}0`$. Without loss of generality, we stipulate that $`𝑲`$ and $`K`$ are homogeneous functions of the momenta of degree two, like $`𝑯`$ and $`H`$. Starting from the expansion (see, e.g., \[Brenig92, HOSW84\])

$$p_x,p_y|i[𝑯,𝑲]|x,yx,y|p_x,p_y=\mathrm{}\{H,K\}+O(\mathrm{}^2),$$
(4)
we will exploit only the fact that the Poisson bracket $`\{H,K\}`$ is everywhere zero except at isolated critical points in position space. The following line of reasoning is therefore also true for integrable billiards with a finite number of magnetic flux lines or point scatterers (we have checked this for Šeba’s billiard \[Seba90\]). In the semiclassical regime, both sides of Eq. (4) are very small in the region $`𝒬𝒞`$ excluding the union $`𝒞`$ of neighborhoods of the critical corners, the area of which shrinks to zero as $`\mathrm{}0`$. Manipulating Eq. (2) and restricting the integration over an $`|L`$-state to region $`𝒬𝒞`$ gives

$$(\mathrm{\Delta }H_{𝒬𝒞})^2+(\mathrm{\Delta }K_{𝒬𝒞})^2=i[𝑯,𝑲]_{𝒬𝒞},$$
(5)
with $`(\mathrm{\Delta }H_{𝒬𝒞})^2=(𝑯𝑯)^2_{𝒬𝒞}`$, etc (Note that $`1_{𝒬𝒞}<1`$). The smallness of the l.h.s. of Eq. (4) carries over to both sides of Eq. (5). From this we conclude that a function $`x,y|L`$ in region $`𝒬𝒞`$ is locally either very small or can be approximated by a joint eigenfunction of both operators $`𝑯`$ and $`𝑲`$. Such a joint eigenfunction cannot fulfill the boundary conditions globally, otherwise the billiard would be integrable. Loosely speaking, region $`𝒬𝒞`$ must be divided into subregions in each of which a joint eigenfunction (or a vanishing function) fulfills the boundary conditions locally. Joint eigenfunctions of neighboring subregions match smoothly, so they must have roughly the same number of nodal lines. Hence, $`|L`$-states have features of energy states of integrable systems: (i) the functions $`x,y|L`$ are in some regions regular while in other regions, separated by “caustics”, vanishing; (ii) they can be labelled by two “quantum numbers” $`(n_1,n_2)`$; (iii) the eigenvalues $`L=𝑯i𝑲`$ are regularly distributed in the complex plane. The last property holds because $`𝑯`$ (all arguments are also valid for quantities derived from $`𝑲`$) is asymptotically equal to $`𝑯_{𝒬𝒞}`$, which is approximately a homogeneous function of the quantum numbers of degree two due to the homogeneity of $`𝑯`$. The non-relevance of $`𝑯_𝒞`$ can be understood with a renormalization procedure. Going from $`(n_1,n_2)`$ to $`(2n_1,2n_2)`$, we get a new wave function with slightly larger region $`𝒬𝒞`$ containing four times more nodal lines, and region $`𝒞`$ being a four times smaller copy of the old region $`𝒞`$. Hence, $`𝑯_{𝒬𝒞}`$ roughly increases by a factor of four, while $`𝑯_𝒞`$ stays constant. This leads to $`𝑯𝑯_{𝒬𝒞}`$ as more and more renormalization steps bring us towards the semiclassical limit.

The approximate homogeneity of $`𝑯_{𝒬𝒞}`$ with respect to the quantum numbers ensures that the local transitions between different joint eigenfunctions induce only a small uncertainty $`(\mathrm{\Delta }H_{𝒬𝒞})^2`$ of order $`𝑯𝑯^2`$, increasing roughly by a factor of four under the action of the renormalization. The same factor is an upper bound for the increase of $`(\mathrm{\Delta }H_𝒞)^2`$, corresponding to a $`𝑯^2`$-scaling weighted with the size of $`𝒞`$. The sum $`(\mathrm{\Delta }H)^2=(\mathrm{\Delta }H_{𝒬𝒞})^2+(\mathrm{\Delta }H_𝒞)^2`$ therefore increases by a factor of four. Hence,

$$(\mathrm{\Delta }H)^2𝑯,(\mathrm{\Delta }K)^2𝑲$$
(6)
in the semiclassical limit. Consequently,

$$\underset{\mathrm{}0}{lim}\frac{\mathrm{\Delta }H}{𝑯}0,\underset{\mathrm{}0}{lim}\frac{\mathrm{\Delta }K}{𝑲}0,$$
(7)
i.e. the $`|L`$-states condense onto the invariant surfaces.

## III Time evolution

The time dependence of $`|L`$-states is non-trivial since $`𝑳`$ does not commute with the Hamiltonian; $`|L(t)`$ is in general not an eigenstate of $`𝑳`$ for $`t>0`$. Let us define three time scales associated with a given state $`|L`$ at time $`t=0`$: the quantum mean period, the lifetime and the classical mean free time

$$\tau _\text{Q}=\frac{2\pi \mathrm{}}{𝑯},\tau _\text{L}=\frac{2\pi \mathrm{}}{\mathrm{\Delta }H},\tau _\text{C}=\sqrt{\frac{A}{2𝑯}},$$
(8)
where $`A`$ is the area of the billiard. From Eq. (6) we get $`\tau _\text{Q}\tau _\text{L}\tau _\text{C}`$, i.e. the state is quasi-stationary with a lifetime of order of the classical mean free time.

At first glance, it seems that the state fails to condense onto an invariant surface for classical long times $`t\tau _\text{C}`$. However, the requirements for condensation, small relative uncertainties and constant mean values of $`𝑯`$ and $`𝑲`$, may persist beyond the lifetime of the state. Clearly, this is true for the relative uncertainty and the mean value of $`𝑯`$ since this operator commutes with the evolution operator $`\mathrm{exp}(i𝑯t/\mathrm{})`$. Keeping in mind that $`i[𝑯,𝑲]`$ is relatively small, we may expect that the relative uncertainty and the mean value of $`𝑲`$ change slower than the state itself. To see this, we define the time scale associated with

$`𝑲(t)`$ $`=`$ $`L(t)|𝑲|L(t)`$ (9)
$`=`$ $`{\displaystyle \underset{j,k=1}{\overset{\mathrm{}}{}}}L|E_jK_{jk}e^{i(E_jE_k)t/\mathrm{}}E_k|L.`$ (10)
in the same way as the lifetime in Eq. (8) as $`\tau _\text{K}=2\pi \mathrm{}/\mathrm{\Delta }H_K`$, where

$$(\mathrm{\Delta }H_K)^2=\frac{_{j,k}|L|E_jK_{jk}(E_jE_k)E_k|L|^2}{2_{j,k}|L|E_jK_{jk}E_k|L|^2}$$
(11)
is the mean energy difference, i.e. mean frequency difference, of the energy states involved, but in contrast to $`\mathrm{\Delta }H`$ the states are weighted according to their relevance for $`𝑲(t)`$. Two extreme cases show that definition (11) is reasonable: a diagonal matrix $`K_{jk}`$ gives $`\mathrm{\Delta }H_K=0`$ whereas a uniform matrix $`K_{jk}=\text{constant}`$ gives $`\mathrm{\Delta }H_K=\mathrm{\Delta }H`$.

The series (9) can be expanded in orders of $`t`$, with $`t`$ much smaller than $`\tau _\text{K}`$ but with the possibility of being larger than the quantum mean period and even the lifetime. Clearly, the zeroth-order term is $`𝑲`$. The next terms are roughly of order $`𝑲\mathrm{\Delta }H_Kt`$, $`𝑲(\mathrm{\Delta }H_Kt)^2`$, and so on. If we compare this to

$$𝑲(t)=𝑲+i[𝑯,𝑲]t/\mathrm{}+O(t^2),$$
(12)
we find that the order of $`\tau _\text{K}`$ can be estimated as $`2\pi 𝑲/i[𝑯,𝑲]`$. It follows that in the semiclassical limit $`𝑲(t)`$ does not change relative to its initial value for classically long times $`t`$ below $`\tau _\text{K}\tau _\text{C}`$, even though $`|L(t)`$ may differ strongly from the initial state $`|L`$.

Following the same line of reasoning shows that the time below which $`𝑲^2(t)`$ relative to its initial value is constant scales as $`𝑲^2/i[𝑯,𝑲^2]`$. Starting from Eq. (2) the following equation can be derived

$$i[𝑯,𝑲^2]2𝑲i[𝑯,𝑲]=2(𝑲𝑲)^3.$$
(13)
With the renormalization procedure it can be shown that the asymptotic behavior of the r.h.s of this equation is bounded from above by $`𝑲^2`$. Together with the already known result $`i[𝑯,𝑲]𝑲`$ we finally get the upper bound $`i[𝑯,𝑲^2]𝑲^2`$. Hence, the time scale which governs $`𝑲^2(t)`$ is of the same or of larger order as $`\tau _\text{K}`$. From this we can conclude that for times smaller than $`\tau _\text{K}`$, $`\mathrm{\Delta }K(t)`$ remains small if compared with $`𝑲(t)`$. Hence, Eq. (7) holds not only for times of order $`\tau _\text{L}\tau _\text{C}`$ but also for times much larger than $`\tau _\text{C}`$ and well below $`\tau _\text{K}`$. The condensation of $`|L`$-states onto invariant surfaces outlives the lifetime of the states up to classically long times.

## IV Example: the barrier billiard

We illustrate all statements for the barrier billiard \[Zwanzig83, HM90, Wiersig00\], a rectangle with width $`l_x`$ and height $`l_y`$ and a vertical barrier of length $`l_y/2`$ placed on the symmetry line $`x=l_x/2`$; see Fig. 1(a). The billiard is not only pseudointegrable, it is also almost-integrable \[Gutkin86\], i.e. is composed of several copies of a single integrable sub-billiard, here the rectangle shown in Fig. 1(b). The function $`K=p_x^2`$ is a second constant of motion. The general formula for the genus of the invariant surfaces \[RichensBerry81\] gives 2, i.e. the surfaces have the topology of two-handled spheres and not that of tori (single-handled spheres).

The energy eigenfunctions are solutions of the Helmholtz equation with Dirichlet boundary conditions on the polygon. The functions are odd or even with respect to the symmetry line. The odd ones are trivial eigenfunctions of the integrable sub-billiard. We therefore deal only with the even ones, which fulfill mixed boundary conditions on the symmetry reduced polygon; see Fig. 1(b). We have calculated the solutions numerically with the mode-matching method for the parameters $`\mathrm{}=1`$, $`l_x=\pi \sqrt{8\pi }/3`$, and $`l_y=3\sqrt{8\pi }/\pi `$ as described in App. 1. The statistical properties of the energy levels are found to be close to the semi-Poisson distribution \[BGS99\]. The energy eigenfunctions are not eigenfunctions of the operator $`𝑲=𝒑_x^2`$. The numerical computation of the uncertainty $`\mathrm{\Delta }K`$ is explained in App. 2. Figure 2 shows that there is no trend towards vanishing relative uncertainties. This indicates that energy states of the barrier billiard do not condense onto the invariant surfaces.

The nontrivial $`|L`$-states are also calculated with the mode-matching method. As explained in App. 3, we cannot compute as many of these states as energy states. However, even in the accessible low-energy regime our theoretical results on the semiclassical behavior of these states will be well confirmed in the following.

The regular pattern of the eigenvalues $`L=𝑯i𝑲`$ can be most clearly seen when transformed into the real action variables of the integrable sub-billiard

$$(I_x,I_y)=(\frac{l_x}{2\pi }\sqrt{𝑲},\frac{l_y}{\pi }\sqrt{𝑯𝑲}).$$
(14)
Note that the topology of the invariant surfaces in the full system rules out action-angle variables \[RichensBerry81\]. As can be seen from Fig. 3 the “action space” is split into two regions A and B separated by a transition region. Away from this transition region, and the lines $`I_x=0`$ and $`I_y=0`$, the eigenvalues are approximately located on regular lattices which are given by EBK-like quantization rules:

$$(I_x,I_y)=(n_1+\frac{3}{4},n_2+1)\mathrm{}$$
(15)
in region A and

$$(I_x,I_y)=(\frac{n_1}{2}+\frac{1}{2},2n_2+\frac{7}{4})\mathrm{}$$
(16)
in region B with $`n_1,n_2=0,1,\mathrm{}`$. For example, we have $`(I_x,I_y)(6.75,37)\mathrm{}`$ for the eigenfunction of type A and $`(I_x,I_y)(21,9.7)\mathrm{}`$ for the function of type B in Fig. 4. Both types have relatively small uncertainties (see App. 4 for numerical details) and look rather regular. While type-A functions cover the billiard uniformly, apart from a localization around the critical corner, the type-B functions are restricted to the lower ($`n_1`$ odd) or upper ($`n_1`$ even) half of the billiard, bounded by a caustic-like curve.

The eigenvalue pattern in Fig. 3 resembles strongly that of integrable systems with separatrices \[RDWW96, WWD97, WWD98\]. It is therefore not surprising that the statistical properties of the mean energies of $`|L`$-states are similar to those of energy levels of integrable systems. For example, Fig. 5 confirms that the nearest-neighbor statistics is in agreement with the Poisson distribution.

The scaling laws (6) for the uncertainties are verified in the following way. First note that lines of constant mean energy $`𝑯`$ are circles in the scaled action space $`(2I_x/l_x,I_y/l_y)`$. This plane can therefore be conveniently parametrized by the radial coordinate $`𝑯`$ and the polar angle $`\varphi =\mathrm{arctan}[(I_y/l_y)/(2I_x/l_x)]`$. The first scaling law holds if $`(\mathrm{\Delta }H)^2`$ divided by $`𝑯`$ is a function of $`\varphi `$ alone. Figure 6 shows that this condition is well satisfied. This confirms also the second scaling law since $`\mathrm{\Delta }K=\mathrm{\Delta }H`$ and $`𝑲`$ scales as $`𝑯`$.

The main features of the time dependence of $`|L`$-states can be observed in Figs. 7 and 8; see App. 5 for the numerical aspects. For small times $`t`$ around $`\tau _\text{Q}\tau _\text{L}26.7\tau _\text{Q}`$, the time evolved state $`|L(t)`$ and the initial state $`|L`$ are similar; the overlap $`|L(t)|L|^2`$ is close to one. For times $`t\tau _\text{L}`$ both states differ considerably; the overlap is close to zero. The state has lost some of its regularity, but $`𝑲(t)`$ and $`\mathrm{\Delta }K(t)`$ stay constant up to order $`𝑲`$ for longer times well below $`\tau _\text{K}2736\tau _\text{Q}`$. Interestingly, the states of type B become more uniformly distributed in configuration space in the course of time, cf. Figs. 4(b), 8(a) and (b).

## V conclusion

We have formulated the quantum-classical correspondence in rational polygonal billiards in the following way: there exists a basis of quantum states with each state condensing individually onto a classical invariant surface $`H,K=\text{constant}`$ up to classically long times in the sense that the relative uncertainties $`\mathrm{\Delta }H/𝑯`$ and $`\mathrm{\Delta }K/𝑲`$ vanish in the semiclassical limit. We have presented some hints, like the analogy to the Anderson metal-insulator transition, and numerical evidence for a particular system, the barrier billiard, that the energy states do not show such a condensation. We have then introduced an alternative basis of states for which we have explicitly shown that (i) they are quasi-stationary; (ii) they condense onto the invariant surfaces up to classically long times even exceeding their lifetimes; (iii) in configuration space they show a regular nodal-line structure, possibly with caustic-like boundaries; (iv) the eigenvalues form a regular pattern in the complex plane.

Whether these states condense uniformly onto the invariant surfaces in phase space is another question. At first sight, our numerical results on the barrier billiard seem to indicate that there cannot be uniform condensation since a fraction of the states cover only one half of the configuration space. However, each such a state becomes more uniformly distributed on the invariant surface for times beyond the lifetime but remains in the neighborhood of the surface for classically long times, i.e. there can be uniform condensation after some transient time.

It is important to note that the condensation scenario holds for general rational polygonal billiards and also if the time-dependent Schrödinger equation is replaced by the wave equation used in acoustics and electromagnetics, even though the individual time scales are different.

I would like to thank M. Sieber and J. Nöckel for discussions.

## Numerical computations on the barrier billiard

### 1 Energy states

We compute the energy eigenvalues and eigenfunctions with the mode-matching method; see, e.g.,\[RichensBerry81, WSM95, ESTV96\]. Let us first set $`\mathrm{}=1`$ and then divide the symmetry-reduced barrier billiard in two regions as shown in Fig. 1(b): region 1 with $`yl`$ and region 2 with $`y>l`$. The length of the barrier, $`l`$, is in our case fixed to $`l_y/2`$ but the following derivations hold also for general $`0<l<l_y`$. In region 1 we expand the wave function as

$$\mathrm{\Phi }_1=\underset{m=1}{\overset{\mathrm{}}{}}a_m\mathrm{sin}(2m\pi x/l_x)\mathrm{sin}(g_{2m}y)$$
(17)
with

$$g_j^2(E)=E\left(\frac{j\pi }{l_x}\right)^2.$$
(18)
This kind of expansion is nontrivial since $`g_j`$ can be imaginary and $`a_m`$ complex. By construction, the function (17) fulfills the Helmholtz equation $`^2\mathrm{\Phi }_1=E\mathrm{\Phi }_1`$ with Dirichlet boundary condition on $`x=0`$, $`y=0`$, and $`x=l_x/2`$. In region 2 we take an analog function,

$$\mathrm{\Phi }_2=\underset{m=1}{\overset{\mathrm{}}{}}b_m\mathrm{sin}[(2m1)\pi x/l_x]\mathrm{sin}[g_{2m1}(l_yy)],$$
(19)
which satisfies Dirichlet boundary condition on $`x=0`$ and $`y=l_y`$ but Neumann boundary condition on $`x=l_x/2`$. We stipulate that both functions match smoothly at $`y=l`$, i.e. for $`0xl_x/2`$ we require

$$\mathrm{\Phi }_1(x,l)=\mathrm{\Phi }_2(x,l)$$
(20)
and

$$\frac{\mathrm{\Phi }_1}{y}|_{(x,l)}=\frac{\mathrm{\Phi }_2}{y}|_{(x,l)}.$$
(21)
Inserting the identity

$$\mathrm{sin}(2nz)=\underset{k=1}{\overset{\mathrm{}}{}}A_{nk}\mathrm{sin}[(2k1)z]$$
(22)
for $`0z\pi /2`$ with the orthogonal matrix

$$A_{nk}=\frac{2}{\pi }\left(\frac{\mathrm{sin}[(2n2k+1)\pi /2]}{2n2k+1}\frac{\mathrm{sin}[(2n+2k1)\pi /2]}{2n+2k1}\right)$$
(23)
into the first matching condition (20) and solving for the coefficients of $`\mathrm{sin}[(2n1)\pi x/l_x]`$, $`n=1,\mathrm{}`$, gives

$$\underset{m=1}{\overset{\mathrm{}}{}}a_m\mathrm{sin}(g_{2m}l)A_{mn}=b_n\mathrm{sin}[g_{2n1}(l_yl)].$$
(24)
Similarly, we get from the second matching condition (21)

$$\underset{m=1}{\overset{\mathrm{}}{}}a_mg_{2m}\mathrm{cos}(g_{2m}l)A_{mn}=b_ng_{2n1}\mathrm{cos}[g_{2n1}(l_yl)].$$
(25)
We now rewrite the relations (24)-(25) by using the definitions

$`\stackrel{~}{a}_m`$ $`=`$ $`a_mg_{2m}\mathrm{cos}(g_{2m}l),`$ (26)
$`\stackrel{~}{b}_m`$ $`=`$ $`b_mg_{2m1}\mathrm{cos}[g_{2m1}(l_yl)]`$ (27)
and the real function

$$f_m(E,l)=\frac{\mathrm{tan}(g_ml)}{g_m}$$
(28)
as

$$\underset{m=1}{\overset{\mathrm{}}{}}\stackrel{~}{a}_mf_{2m}(E,l)A_{mn}=\stackrel{~}{b}_nf_{2n1}(E,l_yl)$$
(29)
and

$$\underset{m=1}{\overset{\mathrm{}}{}}\stackrel{~}{a}_mA_{mn}=\stackrel{~}{b}_n.$$
(30)
The last two equations are combined to

$$\underset{m=1}{\overset{\mathrm{}}{}}\stackrel{~}{a}_mM_{mn}=0$$
(31)
with the real matrix

$$M_{mn}(E)=[f_{2m}(E,l)+f_{2n1}(E,l_yl)]A_{mn}.$$
(32)
Equation (31) has a solution provided $`detM(E)=0`$. Finding the energy eigenvalues $`E_i`$ is therefore equivalent to finding the zeroes of the determinant of $`M`$. We approximate $`M`$ by a $`500\times 500`$ matrix which is sufficient for calculating the first $`\mathrm{100\hspace{0.17em}000}`$ zeroes. Figure 9 shows the determinant as a function of the energy. It is convenient to search roots numerically only between two consecutive poles. We therefore rewrite Eq. (28) using the identity \[GradRyzh65\]

$$\frac{\pi }{4z}\mathrm{tan}\frac{\pi z}{2}=\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{z^2(2n1)^2}$$
(33)
as

$$f_m(E,l)=\frac{2}{l}\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{E(\frac{m\pi }{l_x})^2(\frac{n\pi }{2l})^2}$$
(34)
where the summation is only over odd $`n`$. From this expression the poles can be easily read off. Note that the poles are not degenerate if $`l^2/l_x^2`$ and $`(l_yl)^2/l_x^2`$ are irrational numbers. The interval between two given consecutive poles is divided into 400 subintervals, each assumed to contain at most a single change of sign of $`detM`$. The bisection method is then employed in order to find each zero with an accuracy $`10^4`$ of the mean level spacing.

Having determined the energy levels, we get for each of them the real quantities $`\stackrel{~}{a}_m`$ and $`\stackrel{~}{b}_m`$ from Eqs. (30) and (31) and then the wavefunction (17) and (19), which we finally normalize to unity for the full billiard.

### 2 Uncertainties of energy states

We here determine the uncertainty $`(\mathrm{\Delta }K)^2=𝑲^2𝑲^2`$ in a given $`|E_j`$-state. We take $`\mathrm{\Phi }_1=\mathrm{\Phi }_1(E_j)`$, $`g=g(E_j)`$, and $`\mathrm{\Phi }_2=\mathrm{\Phi }_2(E_j)`$ from Eqs. (17)-(19) to compute $`K_{jj}=E_j|𝑲|E_j=𝑲`$ from

$`K_{jj}`$ $`=`$ $`2{\displaystyle _0^{l_x/2}}{\displaystyle _0^l}\mathrm{\Phi }_1^{}\left({\displaystyle \frac{^2}{x^2}}\mathrm{\Phi }_1\right)𝑑x𝑑y`$ (35)
$`+`$ $`2{\displaystyle _0^{l_x/2}}{\displaystyle _l^{l_y}}\mathrm{\Phi }_2^{}\left({\displaystyle \frac{^2}{x^2}}\mathrm{\Phi }_2\right)𝑑x𝑑y.`$ (36)
A straightforward calculation gives

$`K_{jj}`$ $`=`$ $`{\displaystyle \frac{l_x}{4}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}a_m^2h_{2m}\left[l{\displaystyle \frac{\mathrm{sin}(2g_{2m}l)}{2g_{2m}}}\right]`$ (37)
$`+`$ $`{\displaystyle \frac{l_x}{4}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}b_m^2h_{2m1}\left[(l_yl){\displaystyle \frac{\mathrm{sin}[2g_{2m1}(l_yl)]}{2g_{2m1}}}\right]`$ (38)
with $`a_m=a_m(E_j)`$, $`b_m=b_m(E_j)`$, $`g_m=g_m(E_j)`$, and $`h_m=(m\pi /l_x)^2`$. Analogously, we get for $`𝑲^2=𝒑_x^4`$ again Eq. (37) but with $`h_m=(m\pi /l_x)^4`$.

### 3 $`|L`$-states

Eigenvalues and eigenfunctions of $`𝑳`$ can be computed with the procedure described in App. 1 if Eq. (18) is replaced by

$$g_j^2(L)=L(1i)\left(\frac{j\pi }{l_x}\right)^2$$
(39)
with complex number $`L`$. $`M`$ in Eq. (32) is then a complex matrix. Finding the complex roots of $`detM`$ is much more cumbersome than finding real roots as in the case of the energy states. Because of this, we first compute $`𝑳`$ in energy-state representation, i.e. we calculate the matrix elements $`E_j|𝑳|E_k=E_j\delta _{jk}iK_{jk}`$ with $`K_{jk}=E_j|𝑲|E_k`$ given by Eq. (37) if $`j=k`$, otherwise

$`K_{jk}`$ $`=`$ $`{\displaystyle \frac{l_x}{4}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}a_m^{}a_mh_{2m}\left[{\displaystyle \frac{\mathrm{sin}(g_{2m}^{}l)}{g_{2m}^{}}}{\displaystyle \frac{\mathrm{sin}(g_{2m}^+l)}{g_{2m}^+}}\right]`$ (40)
$`+`$ $`{\displaystyle \frac{l_x}{4}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}b_m^{}b_mh_{2m1}\times `$ (42)
$`\left[{\displaystyle \frac{\mathrm{sin}[g_{2m1}^{}(l_yl)]}{g_{2m1}^{}}}{\displaystyle \frac{\mathrm{sin}[g_{2m1}^+(l_yl)]}{g_{2m1}^+}}\right]`$
with $`a_m^{}=a_m^{}(E_j)`$, $`a_m=a_m(E_k)`$, $`b_m^{}=b_m^{}(E_j)`$, $`b_m=b_m(E_k)`$, $`g_m^+=g_m^{}(E_j)+g_m(E_k)`$, $`g_m^{}=g_m^{}(E_j)g_m(E_k)`$, and $`h_m=(m\pi /l_x)^2`$. A LAPACK routine is used to diagonalize the complex matrix $`E_j|𝑳|E_k`$ with $`j,k1000`$, giving a rough approximation to the first 800 eigenvalues of $`𝑳`$. We use them as initial guesses for Newton’s method in order to find reliable approximations to the eigenvalues. Finally, we compute for each determined eigenvalue the complex quantities $`\stackrel{~}{a}_m`$ and $`\stackrel{~}{b}_m`$ from Eqs. (30) and (31) and then the normalized wavefunction (17) and (19).

### 4 Uncertainties of $`|L`$-states

The main advantage of using the mode-matching method for the $`|L`$-states is that we get the uncertainties $`\mathrm{\Delta }K`$ and $`\mathrm{\Delta }H`$ as accurately as for the energy states in App. 2. We take $`\mathrm{\Phi }_1=\mathrm{\Phi }_1(L_j)`$, $`\mathrm{\Phi }_2=\mathrm{\Phi }_2(L_j)`$, and $`g=g(L_j)`$ from Eqs. (17), (19), and (39) to compute firstly $`𝑲`$ starting from Eq. (35). A straightforward calculation shows that $`𝑲`$ is given by the r.h.s. of Eq. (42) with $`a_m^{}=a_m^{}(L_j)`$, $`a_m=a_m(L_j)`$, $`b_m^{}=b_m^{}(L_j)`$, $`b_m=b_m(L_j)`$, $`g_m^+=g_m^{}(L_j)+g_m(L_j)`$, $`g_m^{}=g_m^{}(L_j)g_m(L_j)`$, and $`h_m=(m\pi /l_x)^2`$. Similarly, we get for $`𝑲^2=𝒑_x^4`$ the same equation but with $`h_m=(m\pi /l_x)^4`$; for $`𝑯`$ we use $`h_m=g_m^2(L_j)+(m\pi /l_x)^2`$ and for $`𝑯^2`$ we use $`h_m=|g_m^2(L_j)+(m\pi /l_x)^2|^2`$.

### 5 Time dependence of $`|L`$-states

We here compute $`x,y|L(t)`$, $`L(t)|L`$, $`𝑲(t)`$ and $`𝑲^2(t)`$. These expressions can be written in terms of energy states as

$`x,y|L(t)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}e^{iE_jt}E_j|Lx,y|E_j,`$ (43)
$`L(t)|L`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}e^{iE_jt}|E_j|L|^2,`$ (44)
and as in Eq. (9). The first 1400 energy states are incorporated in these sums. As in the previous sections it can be shown that $`E_j|L`$ is given by the r.h.s. of Eq. (42) with $`a_m^{}=a_m^{}(E_j)`$, $`a_m=a_m(L)`$, $`b_m^{}=b_m^{}(E_j)`$, $`b_m=b_m(L)`$, $`g_m^+=g_m^{}(E_j)+g_m(L)`$, $`g_m^{}=g_m^{}(E_j)g_m(L)`$, and $`h_m=1`$. $`E_j|𝑲^2|E_k`$ is also given by the r.h.s. of Eq. (42) if $`jk`$ with $`a_m^{}=a_m^{}(E_j)`$, $`a_m=a_m(E_k)`$, $`b_m^{}=b_m^{}(E_j)`$, $`b_m=b_m(E_k)`$, $`g_m^+=g_m^{}(E_j)+g_m(E_k)`$, $`g_m^{}=g_m^{}(E_j)g_m(E_k)`$, and $`h_m=(m\pi /l_x)^4`$. For $`j=k`$ one has to use Eq. (37).


format: text

text:
# Two-state system driven by imperfect 𝜋 pulses: an estimate of the error accumulation in bang-bang control methods

## I Introduction and model

A number of dynamical methods for the suppression of decoherence in quantum systems have recently been proposed . While some of these are quite sophisticated, the simplest such approach, sometimes called “bang-bang” control , is a rather straightforward idea which could, in principle, be applied to any two-state system interacting with an environment which has a finite (nonzero) correlation time $`\tau _c`$. (See  for a recent proof-of-principle experimental demonstration.) The idea is to “flip” rapidly (faster than $`\tau _c`$) the state of the system, back and forth, in such a way that the environment’s unwanted influence on the system is constantly being undone by the environment itself.

For instance, for a system to environment coupling which involves only the operator $`\sigma _z`$ (pure phase decoherence), due to a term such as

$$\underset{k}{}\mathrm{}\sigma _z\left(g_kb_k^{}+g_k^{}b_k\right)$$
(1)
in the Hamiltonian (where the boson operators $`b_k`$ represent modes of the environment), a sequence of very short $`\pi `$ pulses, which rotate the system’s pseudospin by $`180^{}`$ around the $`x`$ axis and hence change the sign of $`\sigma _z`$, would result in an evolution in which the sign of (1) changes form one instant to the next, and its effect therefore averages to zero.

The purpose of this note is to consider the constraints which imperfections in the $`\pi `$ pulses place on the successful implementation of such a dynamical decoupling strategy. The approach is simply to consider the evolution of a two-state system subject to a series of imperfect $`pi`$ pulses, and to calculate the total error (as measured by the state’s fidelity) introduced, on the average, after $`N`$ such pulses. The average considered here is over all possible initial states of the system, as well as over the distribution of the random errors in the $`\pi `$ pulses.

I work in the interaction picture and assume that the central frequency of the pulses is exactly tuned to the resonance frequency of the transition between the two states of the system. The rotating-wave approximation is also assumed to apply. The Hamiltonian describing the interaction of the system with a single pulse is then

$$H=\mathrm{}g\left((t)\sigma _++^{}(t)\sigma _{}\right)$$
(2)
where $`\sigma _+=|+|`$, $`\sigma _{}=|+|`$, and $`|+`$ and $`|`$ are the eigenstates of $`\sigma _z`$. The coupling constant $`g`$ can be taken to be real without loss of generality. In general, the Hamiltonian (2) does not commute with itself at different times. To simplify matters, I shall take the classical field amplitude to be of the form $`(t)=E(t)e^{i\phi }`$, where $`E(t)`$ is real and $`\phi `$ is a constant, determined essentially by the timing of the pulse. Then, the evolution operator corresponding to (2) is, in the {$`|+`$,$`|`$} basis,

$`U(t)`$ $`=\mathrm{exp}\left[{\displaystyle \frac{i}{\mathrm{}}}{\displaystyle _0^t}H(t^{})𝑑t^{}\right]`$ (4)
$`=\mathrm{exp}\left[i\left(\mathrm{cos}\phi \sigma _x\mathrm{sin}\phi \sigma _y\right){\displaystyle _0^t}gE(t^{})𝑑t^{}\right]`$
If the duration of the pulse is $`t_p`$, define the angle $`\mathrm{\Theta }2_0^{t_p}gE(t^{})𝑑t^{}`$. In pseudospin terms, the operator $`U(t_p)`$ represents a rotation of the spin by an angle $`\mathrm{\Theta }`$ about an axis in the $`xy`$ plane which makes an angle $`\phi `$ with the $`x`$ axis. If $`\mathrm{\Theta }=\pi `$, one has a $`\pi `$ pulse.

In this paper I consider the state resulting from a sequence of “imperfect” $`\pi `$ pulses. Imperfections may arise from either the phase $`\phi `$, which controls the axis of the rotation, and which would be sensitive to errors in the timing of the pulses, or from errors in the amplitude and/or duration of the pulse, which would make $`\theta \pi `$. Amplitude errors are considered in Section 2, phase errors in Section 3, and the results, and their possible consequences for “bang-bang” control, are discussed in Section 4.

## II Amplitude errors

In this Section I shall assume that $`\phi =0`$, in Eq. (4), for all the pulses in the sequence. A full spin cycle consists of two $`\pi `$ pulses (separated by a short time of “free” evolution), so I shall consider a total of $`2N`$ pulses, where $`N`$ is the total number of cycles. I will take the “pulse area” $`\mathrm{\Theta }`$ to be of the form $`\mathrm{\Theta }=\pi +ϵ_i`$, where $`ϵ_i`$ is a small number which varies randomly from one pulse to the next. All the $`ϵ_i`$ are assumed to be distributed in a Gaussian way, that is, according to the following probability distribution:

$$P(ϵ_i)=\frac{1}{\sqrt{2\pi }\mathrm{\Delta }}e^{ϵ_i^2/2\mathrm{\Delta }^2}$$
(5)
with standard deviation $`\mathrm{\Delta }`$.

With this notation the evolution operator (4) for the $`i`$-th pulse is simply

$$U_i=e^{i(\pi /2+ϵ_i/2)\sigma _x}$$
(6)
and the result of $`2N`$ consecutive pulses is given by

$`U_{total}`$ $`=e^{i(N\pi +ϵ_1/2+\mathrm{}ϵ_{2N}/2)\sigma _x}`$ (9)
$`e^{i(N\pi +ϵ/2)\sigma _x}`$
$`=(1)^N\left[\mathrm{cos}{\displaystyle \frac{ϵ}{2}}i\mathrm{sin}{\displaystyle \frac{ϵ}{2}}\sigma _x\right]`$
where the random variable $`ϵ=ϵ_1+\mathrm{}ϵ_{2N}`$ is also distributed in a Gaussian way, and has standard deviation $`\sqrt{2N}\mathrm{\Delta }`$.

A general initial state of the two-level system can be written as

$$|\psi _0=\mathrm{cos}\frac{\theta }{2}|++e^{i\varphi }\mathrm{sin}\frac{\theta }{2}|$$
(10)
In this state the spin is aligned along an axis of colatitude $`\theta `$ and azimuth $`\varphi `$. The effect of (9) on the state (10) is

$$U_{total}|\psi _0=(1)^N\left[\left(\mathrm{cos}\frac{ϵ}{2}\mathrm{cos}\frac{\theta }{2}ie^{i\varphi }\mathrm{sin}\frac{ϵ}{2}\mathrm{sin}\frac{\theta }{2}\right)|++\left(\mathrm{cos}\frac{ϵ}{2}\mathrm{sin}\frac{\theta }{2}e^{i\varphi }i\mathrm{sin}\frac{ϵ}{2}\mathrm{cos}\frac{\theta }{2}\right)|\right]$$
(11)
The fidelity, $``$, of the state after all the pulses is given by $`=\left|\psi _0|U_{total}|\psi _0\right|^2`$, and it equals

$$=\mathrm{cos}^2\frac{ϵ}{2}+\mathrm{sin}^2\frac{ϵ}{2}\mathrm{sin}^2\theta \mathrm{cos}^2\varphi $$
(12)
To calculate an average fidelity, Eq. (12) may be integrated over $`\theta `$ and $`\varphi `$, assuming a uniform distribution of initial states over the whole Bloch sphere (that is, with a density $`\mathrm{sin}\theta d\theta d\varphi /4\pi `$), and over $`ϵ`$ with the appropriate probability distribution. Before proceeding to do so, however, note that $`\mathrm{sin}\theta \mathrm{cos}\varphi `$ is just $`x`$ in Cartesian coordinates; hence a reparametrization of the Bloch sphere that exchanges $`x`$ and $`z`$ will turn $`\mathrm{sin}\theta \mathrm{cos}\varphi `$ into $`\mathrm{cos}\theta ^{}`$, which makes the integral over $`\varphi ^{}`$ trivial. The probability distribution for the final state fidelity can thus be written formally as

$`P()`$ $`={\displaystyle \frac{1}{2\mathrm{\Delta }}}{\displaystyle \frac{1}{\sqrt{4N\pi }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}e^{ϵ^2/4N\mathrm{\Delta }^2}𝑑ϵ{\displaystyle _0^\pi }sin\theta ^{}𝑑\theta ^{}\delta \left(\mathrm{cos}^2{\displaystyle \frac{ϵ}{2}}\mathrm{cos}^2\theta ^{}\mathrm{sin}^2{\displaystyle \frac{ϵ}{2}}\right)`$ (14)
$`={\displaystyle \frac{1}{\mathrm{\Delta }}}{\displaystyle \frac{1}{\sqrt{4N\pi }}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle _{\sqrt{}}^{\sqrt{}}}\mathrm{exp}\left[{\displaystyle \frac{1}{N\mathrm{\Delta }^2}}\left(\mathrm{sin}^1\sqrt{{\displaystyle \frac{1}{1x^2}}}n\pi \right)^2\right]{\displaystyle \frac{dx}{\sqrt{(1)(x^2)}}}`$
The last integral may be evaluated numerically for a given value of $`N\mathrm{\Delta }^2`$ and $``$ (with $`0<<1`$); for reasonable values of $`N\mathrm{\Delta }^2`$, only a few terms in the sum over $`n`$, around $`n=0`$, are necessary. The results are shown in Figure 1 for $`N\mathrm{\Delta }^2=0.1`$, $`1`$, and $`10`$, respectively. The divergence as $`1`$ is integrable; the value of $`P()`$ as $`0`$ is finite. The average of $``$ can be calculated easily, and the result is

$$=\frac{2}{3}+\frac{1}{3}e^{N\mathrm{\Delta }^2}$$
(15)
This shows the expected decrease in fidelity with the number of cycles, but it also shows that for sufficiently large $`N\mathrm{\Delta }`$ the average fidelity saturates at the value $`2/3`$. This may look like a relatively large number, but one must keep in mind that for a two-level system a totally random mixed state with density operator $`\rho =1/2`$ would still yield a fidelity of $`0.5`$ when compared to any initial state. Thus, in fact, a fidelity of $`0.67`$ for this system corresponds to an almost random final state.

This is illustrated in Figure 2, which shows three simulations of trajectories starting from three randomly picked states and going through a total of $`N_{max}=400`$ cycles, with $`\mathrm{\Delta }`$ chosen in each case so that $`N_{max}\mathrm{\Delta }^2=0.1`$, $`1`$, and $`10`$, as in Fig. 1 (the corresponding average fidelities, as given by Eq. (15), are $`0.968`$, $`0.789`$ and $`0.667`$). Clearly, the trajectory for $`\sqrt{N}_{max}\mathrm{\Delta }=10`$ is all over the place. (The trajectories shown have been chosen from a relatively large sample of randomly generated trajectories so as to be as “average” as possible, that is, neither unusually “good” nor unusually “bad” for the particular value of $`N_{max}\mathrm{\Delta }^2`$ that they illustrate.)

In this connection it may be worthwhile to consider the “worst-case fidelity” predicted by Eq. (12): it corresponds to initial states having either $`\theta =0,\pi `$, or $`\varphi =\pi /2`$, and, when averaged over $`ϵ`$, yields

$$_{min}=\frac{1}{2}+\frac{1}{2}e^{N\mathrm{\Delta }^2}$$
(16)
For large $`N\mathrm{\Delta }^2`$, this gives $`1/2`$, i.e., the fidelity of the totally random mixed state, as discussed above.

## III Phase errors

In this Section I take the pulse area $`\mathrm{\Theta }2_0^{t_p}gE(t^{})𝑑t^{}`$ to be exactly $`\pi `$, and the phase $`\phi `$ in Eq. (4) to be a small random number, varying from pulse to pulse, with zero average and a Gaussian distribution with standard deviation $`\mathrm{\Delta }`$. Then the evolution operator for the $`i`$-th pulse is

$`U_i`$ $`=e^{i(\pi /2)(\mathrm{cos}\phi _i\sigma _x\mathrm{sin}\phi _i\sigma _y)}`$ (18)
$`=i(\mathrm{cos}\phi _i\sigma _x\mathrm{sin}\phi _i\sigma _y)`$
and its effect on the general state (10) is

$$U_i|\psi _0=e^{i\phi _i}\left(\mathrm{cos}\frac{\theta }{2}|++e^{i(\varphi +2\phi _i)}\mathrm{sin}\frac{\theta }{2}|\right)$$
(19)
that is, except for an overall phase factor, all that happens is that the angle $`\varphi `$ becomes $`\varphi +\phi _i`$. Therefore, after $`N`$ cycles totalling $`2N`$ pulses, the final state will be, up to a phase factor

$$U_{total}|\psi _0=\mathrm{cos}\frac{\theta }{2}|++e^{i(\varphi +ϵ)}\mathrm{sin}\frac{\theta }{2}|$$
(20)
where $`ϵ=2(\phi _1+\mathrm{}+\phi _{2N})`$ is now a random variable with zero average and a Gaussian distribution with standard deviation $`2\sqrt{2N}\mathrm{\Delta }`$. A little trigonometric manipulation shows that the fidelity of the state (20) can be written as

$$=\mathrm{cos}^2\frac{ϵ}{2}+\mathrm{sin}^2\frac{ϵ}{2}\mathrm{cos}^2\theta $$
(21)
which is of the same form as (12) (after the reparametrization of the Bloch sphere which exchanges $`x`$ and $`z`$), and yields the same results (14) and (15) for the probability distribution and average value of the fidelity, only with $`N`$ replaced by $`4N`$ everywhere.

## IV Discussion

From the analysis in the previous sections one can conclude that, if each pulse used for “bang-bang” control has a r.m.s. error of $`\mathrm{\Delta }`$ radians in either the phase or the (integrated) amplitude, in order to avoid loss of fidelity one must (roughly speaking) restrict the total number of cycles to

$$N_{max}\frac{1}{\mathrm{\Delta }^2}$$
(22)
In practice, of course, “bang-bang” control would be applied to a system which has an internal Hamiltonian and is, as well, interacting with other systems, so the overall picture will be rather more complicated than that described by Eq. (2), but there are grounds to believe that the conclusion above should still hold true, approximately, in spite of these complications. For instance, recent calculations of decoherence due to internal interactions in a quantum computer show that the net effect does not depend much on whether the computer is actually carrying out a calculation or merely sitting idle in a predetermined state . Similarly, Miquel et al.  found, in their simulation of a quantum computation with phase drift errors, that the effect of said errors could be well approximated by a formula which did not depend on the details of the computation, only on the number of pulses.

As a simple test, I show in Fig. 3 the results of a simulation for a system with the self-Hamiltonian

$$H_0=\mathrm{}\mathrm{\Omega }\sigma _z$$
(23)
initially prepared in the state $`|\psi _0=\frac{1}{\sqrt{2}}(|++i|)`$, which, without “bang-bang” control, would evolve in time as $`|\psi _t=\frac{1}{\sqrt{2}}(e^{i\mathrm{\Omega }t}|++ie^{i\mathrm{\Omega }t}|)`$, and hence exhibit a fidelity $`=\left|\psi _0|U_{total}|\psi _0\right|^2=\mathrm{cos}^2\mathrm{\Omega }t`$ (solid curve). Low-noise “bang-bang” control (dashed curve, amplitude noise, $`N_{max}\mathrm{\Delta }^2=0.1`$) cancels out this free evolution quite effectively, so the state fidelity remains close to 1 for the times shown, but if the noise is high instead (dotted curve, $`N_{max}\mathrm{\Delta }^2=10`$), the fidelity is completely lost. This simulation assumes that the control pulses are essentially instantaneous and very strong, so that the self-Hamiltonian (23) can be ignored during a control pulse. Also, in order to keep the fidelity close to 1, the interval between pulses $`\mathrm{\Delta }t`$ needs to be short enough; specifically, in this case, shorter than $`2\pi /\mathrm{\Omega }`$. For Figure 3, $`\mathrm{\Delta }t=0.005\pi /\mathrm{\Omega }`$.

It should be pointed out that the initial state chosen for the simulation in Fig. 3 (an eigenstate of $`\sigma _y`$) is just about the “worst” possible choice, that is, it is one of the states that minimize the fidelity as given by Eq. (12) (the average fidelity for such states is given by Eq. (16)). This explains why the simulations in Fig. 3 look somewhat worse than those in Fig. 2. Nonetheless, this initial state was chosen because it gives maximum modulation of the fidelity under the self-Hamiltonian (23).

In Ref.  it was shown that, in order to combat decoherence by the “bang-bang” method, one needs to have at least one complete cycle every $`\tau _c`$, where $`\tau _c`$ is the characteristic “coherence time” of the environment. Putting this condition together with (22) yields the result that the “bang-bang” method with imperfect pulses can protect a quantum information storing (or processing) device for a time no longer than

$$T_{max}<\frac{\tau _c}{\mathrm{\Delta }^2}$$
(24)
where $`\mathrm{\Delta }`$ is the r.m.s. error per pulse. For instance, if the $`\pi `$ pulses are known to be accurate to, say, one part in $`10^5`$, one may expect $`\mathrm{\Delta }\pi \times 10^5`$, and $`T_{max}10^9\tau _c`$

It is important to keep in mind that $`\tau _c`$ is not, in general, equal to the decoherence time of the system itself under its interaction with the environment. As an extreme example, the spontaneous emission of an atom has a decoherence time of the order of $`1/\gamma `$, where $`\gamma `$ may be mega- or gigahertz for an optical transition, but the environment responsible for the decoherence—namely, the quantized electromagnetic field—has a “coherence time” of the order of an optical period or shorter, that is, about $`10^{15}`$ s. On the other hand, for other kinds of decoherence—such as, for instance, the kinds resulting from internal interactions between the qubits—one may expect the relevant $`\tau _c`$ to be of the order of a characteristic interaction frequency, which would also determine the basic decoherence timescale. In that case, “bang-bang” control with accurate pulses could substantially lengthen the decoherence time, as Eq. (24) indicates.

Discussions with C.-P. Yang are gratefully acknowledged. This work has been supported in part by the National Security Agency (NSA) and Advanced Research and Development Activity (ARDA) under Army Research Office (ARO) contract number DAAD19-99-1-0118; and by the National Science Foundation under grant PHY-9802413.