format: text

text:
# Quantum Critical Behavior and Possible Triplet Superconductivity in Electron Doped CoO2 Sheets

## Abstract

Density functional calculations are used to investigate the doping dependence of the electronic structure and magnetic properties in hexagonal NaxCoO2. The electronic structure is found to be highly two dimensional, even without accounting for the structural changes associated with hydration. At the local spin density approximation level, a weak itinerant ferromagnetic state is predicted for all doping levels in the range $`x=0.3`$ to $`x=0.7`$, with competing but weaker itinerant antiferromagnetic solutions. The Fermi surface, as expected, consists of simple rounded hexagonal cylinders, with additional small pockets depending on the $`c`$ lattice parameter. Comparison with experiment implies substantial magnetic quantum fluctuations. Based on the Fermi surface size and the ferromagnetic tendency of this material, it is speculated that a triplet superconducting state analogous to that in Sr2RuO4 may exist here.

The last few years have seen the discovery of a number of novel unconventional superconductors associated with magnetic phases. Focusing on triplet (or likely triplet) superconductors these include UGe2 (Tc $``$ 1K),  URhGe (Tc $``$ 0.25K),  and ZrZn2 (Tc $``$ 0.3K),  where ferromagnetism coexists with superconductivity, and Sr2RuO4 (Tc $``$ 1.5K),  which has a paramagnetic Fermi liquid normal state, but is “near” magnetic phases. Although the exact pairing mechanism has not been established in these materials, it is presumed that spin fluctuations are involved, most probably the quantum critical fluctuations in the materials with co-existing ferromagnetism and superconductivity.  In Sr2RuO4, strong nesting related antiferromagnetic spin fluctuations are found in local density approximation (LDA) calculations and experiment.  In addition ferromagnetic fluctuations may also be present, and if so these would favor a triplet superconducting state.

Recently, Takada and co-workers reported the discovery of superconductivity with Tc $``$ 5K, in the CoO2 layer material, NaxCoO$`{}_{2}{}^{}y`$H2O.  Based on the two dimensional 3$`d`$ transition metal oxide structural motif, they speculate that the superconductivity may be related to that of the cuprate high-Tc superconductors. Here an alternate possibility is discussed, that is the connection with the triplet superconductors mentioned above. The parent compound of NaxCoO$`{}_{2}{}^{}y`$H2O is Na0.5CoO2, which was synthesized and studied in single crystal form by Terasaki and co-workers.  The compound consists of CoO2 layers (Co is in at the center of distorted O octahedra formed by layers of O above and below the Co planes) separated by layers of of Na ions.  The result is a highly anisotropic metal with a large thermopower at room temperature and a substantial renormalization of the Fermi liquid properties at low temperature.  Modifications of this compound consisting of substitutions for the Na layer have been made to obtain technologically useful thermoelectric materials,  and, in fact, Na0.5CoO2 is itself a good thermoelectric at high temperature.  Significantly, the Seebeck coefficient of Na0.5CoO2 as calculated from the LDA band structure using standard kinetic transport theory is in quantitative agreement with experiment, indicating that the Fermi liquid description at higher temperature is robust.  Band structure calculations for the higher performance so-called misfit compounds formed by replacement of the Na layer also reasonably account for room temperature experimental data.

Turning to the low temperature properties, as mentioned, there are substantial renormalizations in Na0.5CoO2, for example in specific heat and susceptibility. LDA calculations indicate a ferromagnetic ground state for this parent compound, while experimentally ferromagnetism is not observed, one possible explanation being that quantum fluctuations suppress the ordered ferromagnetic state.  Interestingly, the related misfit compound, Ca3Co4O9 shows a large negative magnetoresistance at low temperature reaching 35% at 7T, which could be interpreted as a sign of proximity to a ferromagnetic quantum critical point. Other evidence for a proximity to magnetism is the formation of a spin density wave when Cu doped,  and the production of a weak magnetic ground state with strong magnetoresistance upon increasing the Na concentration to $`x=0.75`$.

The reported crystal structure of superconducting NaxCoO$`{}_{2}{}^{}y`$H2O differs from that of the parent compound by a reduction in the Na content to $`x0.35`$ and the intercalation of water molecules around the Na sites. The structure of the CoO2 sheets is very similar to that of the parent compound Na0.5CoO2 – the in-plane $`a`$ lattice parameter is slightly lower (2.823 Å vs. 2.84 Å) as is the height of the O (0.885 Å vs. an LDA value of 0.908 Å in the parent). It should be noted that a slight contraction of the CoO6 octahedra, as implied by these structural differences, is consistent with the lower Na content, in particular a lower electron count on the CoO2 sheets would be expected to lead to a small contraction.

Here, well converged LDA calculations are reported for NaxCo2O4 for $`x`$=0.3,0.5,0.7. In addition, calculations are reported for a strained lattice corresponding to the structure reported for superconducting NaxCoO$`{}_{2}{}^{}y`$H2O, but neglecting the intercalating water. The calculations were done using the general potential linearized augmented planewave method with local orbitals,  as described in Ref. , except that better Brillouin zone samples, corresponding to a minimum special k-points mesh of $`16\times 16\times 2`$ in the hexagonal zone is used here to obtain convergence of the magnetic energies. As discussed in Ref. , a virtual crystal method is used to account for the partially occupied Na site.

Although there is some hybridization, the valence band structure of Na0.5CoO2 consists of three manifolds of bands separated by gaps – a lower lying occupied O 2$`p`$ derived manifold, followed by Co $`t_{2g}`$ and $`e_g`$ manifolds.  As expected from ionic considerations the Fermi energy, $`E_F`$ lies near the top of the $`t_{2g}`$ manifold, which contains 0.5 holes per Co ion. Because of the actual axial site symmetry, the $`t_{2g}`$ manifold can be regarded as consisting of two-fold (also labeled $`e_g`$) and one-fold (labeled $`a_g`$) crystal field states. These overlap, but the top of the $`t_{2g}`$ manifold is primarily of $`a_g`$ character in Na0.5CoO2, with the result that the band structure near $`E_F`$ can be roughly viewed as consisting of one band per Co ion with an filling of 3/4 (1/2 hole per Co). The present LDA calculations for the electronic structure of Na0.7CoO2 and Na0.3CoO2 follow this picture. These were done holding the crystal structure fixed at that of Na0.5CoO2 and varying the Na occupancy in the virtual crystal. The dominant $`a_g`$ character at the top of the $`t_{2g}`$ manifold is, however, lost for the strained lattice, discussed below.

The magnetic properties are similar for the various Na contents in the range $`x=0.3`$ to $`x=0.7`$. Results of fixed spin moment constrained LDA calculations are shown in Fig. 1. In particular, itinerant ferromagnetism is found. In each case, the energy decreases with magnetization until a magnetization at which the band edge is reached in the majority spin. Then the energy increases rapidly reflecting the crystal field induced gap between the $`t_{2g}`$ and $`e_g`$ manifolds.

Thus, independent of $`x`$ in this range, the LDA predicts a ferromagnetic ground state, with a spin moment per Co equal to the number of holes ($`p=1x`$) and a half metallic band structure (here we refer to the hole concentration as the concentration of holes in the $`t_{2g}`$ manifold; without any Na $`p`$=1; Na electron dopes the sheets, which leads to a reduction in $`p`$). The fixed spin moment curves show a shape crossing over from parabolic at low moment to more linear as the band edge is approached. As may be seen, the shapes and initial curvatures for the different doping levels are roughly similar. The trend towards slightly weaker initial curvatures at higher Na concentration is possibly an artifact due to the fixed crystal structure used in the present calculations. It reflects increasing hybridization (increasing band width and decreasing density of states) as charge is added to the CoO2 planes. In reality, the lattice would be expected to expand, perhaps compensating this trend. In any case, for this range of $`x`$, ferromagnetism with a magnetic energy of approximately, $`E(FM)50p`$ in meV/Co and spin moment $`M(FM)=p`$ in $`\mu _B`$/ Co is found. As mentioned, calculations were also done for a strained cell with the structure of the superconducting sample, but without H2O. These calculations were done for $`x`$=0.5 and $`x`$=0.35, the latter corresponding to the experimentally determined doping level. In both cases, the LDA predicted a ferromagnetic state. The magnetic energies were $`E(FM)=`$20 and 27 meV/Co for $`x`$=0.5 and $`x`$=0.35, respectively. Thus the behavior is similar, but the magnetic energies are somewhat larger.

LDA calculations were also done for an antiferromagnetic configuration with the unit cell doubled along one of the in-plane lattice vectors. Thus, within a Co plane, each Co ion has four opposite spin nearest neighbors and two like spin nearest neighbors. At all three doping levels investigated an antiferromagnetic instability was found, but this instability is weaker than the ferromagnetic one. Details of the LDA moments and energies are given in Table I. Essentially, the energy of the antiferromagnetic configuration examined tracks the ferromagnetic energy at a value $``$ 1/4 as large.

The LDA generally provides a good description of itinerant ferromagnetic materials. It is known to fail for strongly correlated oxides where on-site Coulomb (Hubbard) repulsions play an important role in the physics. In such cases, the LDA underestimates the tendency of the material towards local moment formation and magnetism. Here, the LDA is found to predict ferromagnetic ground states for materials that are paramagnetic metals in experiment. While materials in which the LDA substantially overestimates the tendency towards magnetism are rare, a number of such cases have been recently found. These are generically materials that are close to quantum critical points, and include Sc3In,  ZrZn2,  and Sr3Ru2O7 (Ref. ). Sr3Ru2O7 displays a novel metamagnetic quantum critical point,  while, as mentioned, ZrZn2 shows coexistence of ferromagnetism and superconductivity.

Density functional theory is in principle an exact ground state theory. It should, therefore, correctly describe the spin density of magnetic systems. However, common approximations to the exact density functional theory, such as the LDA, neglect Hubbard correlations beyond the mean field level, yielding the underestimated magnetic tendency of strongly Hubbard correlated systems. Overestimates of magnetic tendencies, especially in the LDA are very much less common. Another type of correlations that is missed in these approximations are quantum spin fluctuations. This is because the LDA is parameterized based on electron gases with densities typical for atoms and solids. However, the uniform electron gas is very far from magnetism in this density range. In solids near quantum critical points, the result is an overestimate of the magnetic moments and tendency toward magnetism (i.e. misplacement of the position of the critical point) due to neglect of the quantum critical fluctuations.

The present results for NaxCoO2 show a weak ferromagnetic instability that is robust with respect to doping and structure (note the instability for the strained lattice). Based on this, and the experimentally observed renormalized paramagnetic state, it seems likely that NaxCoO2 is subject to strong ferromagnetic quantum fluctuations of this type, and that these are the reason for the disagreement between the LDA and experimental ground states.

The effects of such quantum fluctuations can be described on a phenomenalogical level using a Ginzburg-Landau theory in which the magnetic properties defined by the LDA fixed spin moment curve are renormalized by averaging with an assumed (usually Gaussian) function describing the beyond LDA critical fluctuations.  Although a quantitative theory allowing extraction of this function from first principles calculations has yet to be established, one can make an estimate based on the LDA fixed spin moment curves as compared with experiment. In particular, NaxCoO2 shows a disagreement between the LDA moment and experiment equal to $`p=1x`$, and has a very steeply rising LDA energy for moments larger than $`p`$. Thus one may estimate an r.m.s. amplitude of the quantum fluctuations of $`\xi \alpha p`$ in $`\mu _B`$, with $`1/2<\alpha <1`$, and most likely closer to 1. These are large values c.f. ZrZn2. It is therefore tempting to associate the superconductivity of NaxCoO$`{}_{2}{}^{}y`$H2O with ferromagnetic quantum critical fluctuations. Considering the simple 2D Fermi surface, which consists of rounded hexagonal cylinders plus small sections,  and the ferromagnetic fluctuations a triplet state like that originally discussed for Sr2RuO4 (Ref. ) seems plausible. I now speculate about the ingredients in a spin fluctuation mediated triplet superconducting state.

Within a spin fluctuation induced pairing approach analogous to that employed for Sr2RuO4 the key ingredient is the integral over the Fermi surface of the k-dependent susceptibility with a function of the assumed triplet symmetry,  i.e. in the simplest case, $`𝐤𝐤^{}/kk^{}`$. For a Fermi surface in the shape of a circular cylinder, radius $`k_F`$, the needed integral is proportional to $`_0^{2\pi }d\theta cos(\theta )V(2k_Fsin(\theta /2))`$, where $`V(k)`$ is the assumed pairing interaction. In any case, for a smooth variation of the spin fluctuations with $`k`$ and a maximum at $`k`$=0 (ferromagnetic), the integral is roughly proportional to $`k_F`$ times the variation of $`V`$ from $`k=0`$ to $`k=2k_F`$. This latter variation depends on the detailed shape of $`V(k)`$, but may be expected to cross over from being proportional to $`k_F^2`$ for small $`k_F`$ to proportional to $`k_F`$ for larger $`k_F`$. Neglecting small Fermi surface sections, $`k_F`$ varies as $`p^{1/2}`$. One possibility for $`V(k)`$ is a function smoothly going from a finite value at $`k=0`$ to near zero at the zone boundary (reflecting the rather weak antiferromagnetic instability relative to the ferromagnetic), with a size at $`k=0`$ given by the LDA ferromagnetic energy ($`p^2`$) or alternately a Hund’s coupling ($`p`$ independent) times $`\xi `$ ($`p`$).

Within this $`p`$-wave scenario it would be quite interesting to measure the variation of the superconducting properties of NaxCoO$`{}_{2}{}^{}y`$H2O as a function of doping level. The above arguments imply a substantial model dependent variation up to the level where proximity to the critical point suppresses $`T_c`$, with the implication that still higher values of $`T_c`$ may be obtained. However, it should be emphasized that the mechanism and superconducting symmetry of NaxCoO$`{}_{2}{}^{}y`$H2O have yet to be established, and in fact, even conventional electron-phonon superconductivity competing with spin fluctuations has not been excluded.

I am grateful for helpful discussions with R. Asahi, I.I. Mazin, A.J. Millis, W.E. Pickett, S.S. Saxena and I. Terasaki. Some computations were performed using facilities of the DoD HPCMO ASC and ARL centers. The DoD-AE code was used for some of this work.


text:
# Effects of Random Noisy Oracle on Search Algorithm Complexity

## I Introduction

Issues of fault tolerance and error correction are of both theoretical and practical interest in designing and implementing quantum algorithms. One approach to diminishing the effects of error focuses on error correcting codes, which can be used to correct errors once they have occurred Shor (1995); Preskill (1998). Another approach which has also proven successful is that of decoherence-free subspaces, which involves mapping a computation onto a subspace that is relatively free from error Bacon et al. (2000). This second approach is of interest because it is an example of passive, rather than active, error correction. The success of the passive approach leads to the question of whether existing quantum algorithms are inherently robust to errors, and, if they are not, whether there might exist modified implementations of these algorithms that are robust to certain errors, i.e. without requiring active error correction.

In this paper, we analyze the robustness of Grover’s algorithm to error. Grover’s search algorithm is the basis for a number of quantum algorithms Abrams and Williams (1999); Novak (2001); Brassard et al. (2000). Most of these algorithms involve the use of an oracle, a black box device which takes as input a quantum state in the computational basis and returns as output some function of that state. Given a perfect, noise-less oracle, Grover’s search algorithm attains a quadratic speedup over classical search algorithms. However, this speed-up is predicated upon the perfect implementation of the oracle. Although, for the purposes of analysis, the oracle is often treated as a “black box” whose inner workings are unknown, any physical implementation of Grover’s algorithm must also include a physical implementation of the oracle and this may be imperfect. Therefore it is of interest to ascertain what effect the accuracy of the oracle implementation has on the overall complexity of the algorithm. We make use of this term here as it is employed in computer science terminology, namely, the complexity of an algorithm is defined as the number of computational steps required in order to achieve a pre-determined fixed probability of success.

Several previous papers have studied the effects of oracle noise on Grover’s algorithm, using various models Long et al. (1999); Hsieh et al. (2002); Pablo-Norman and Ruiz-Altaba (1999); Bae and Kwon (2002); Shapira, Mozes, and Biham (1998). We consider here the random phase error model, addressing it within both discrete- and continuous-time implementations of the search algorithm. The effect of random phase errors on the discrete-time Grover algorithm was already studied numerically in Long et al. (1999). In this paper we derive analytic results for this model and present numerical evidence verifying the validity of these results. We analyze the complexity of the search algorithm as a function of the scaling of the errors, and arrive at bounds on the error that must be satisfied for a constant probability of success. In particular, we find that in order to achieve a constant success probability independent of the library size that is being searched, the oracle error must scale at most as $`N^{1/4}`$ where $`N`$ is the library size. We also analyze the effect of phase errors on a continuous-time analogue of the search algorithm. This is found to be relatively straightforward to study within a density matrix approach. We find that this leads to similar results for the complexity as a function of the scaling of the errors, namely that there is an upper bound of $`O(N^{1/2})`$ on algorithmic complexity for $`\delta 1/4`$ and $`O(N^{12\delta })`$ for $`\delta 1/4`$, when the error scaling goes as $`N^\delta `$. We further show that these bounds are tight; in other words, that the algorithmic complexity has a lower bound of $`\mathrm{\Omega }(N^{1/2})`$ for $`\delta 1/4`$ and $`\mathrm{\Omega }(N^{12\delta })`$ for $`\delta 1/4`$. (Following standard computer science notation, we will use the following to characterize the growth of certain functions: We will say that $`f(n)`$ is bounded above by $`g(n)`$, i.e. that $`f(n)=O(g(n))`$, if there are positive constants $`c`$ and $`k`$ such that $`0f(n)cg(n)`$ for $`nk`$. Similarly $`f(n)`$ is bounded below by $`g(n)`$, i.e. $`f(n)=\mathrm{\Omega }(g(n))`$, if $`0cg(n)f(n)`$ for constants $`c,k0`$ and $`nk`$.)

This result has practical importance in determining the amount of oracle error that is allowable as the size of the library is increased. We will show that our complexity results imply that if the library size is increased by a factor of $`k`$, the oracle error must decrease by a factor of $`k^{1/4}`$ in order to attain a similar probability of success.

The remainder of this paper is organized as follows. Sec. II describes the random phase error model for a noisy oracle, summarizes the Grover search algorithm in a discrete-time implementation, and then derives the effect of the oracle noise on this implementation. Sec. III derives the effect of oracle noise on the continuous-time formulation of Grover’s algorithm proposed by Farhi et al. in Farhi and Gutmann (1998); Farhi et al. (2000). We summarize and present conclusions in Sec. IV.

## II Discrete-Time Search

### II.1 Quantum Search with a Phase Oracle

A phase oracle is a quantum oracle that “marks” one or more computational basis states with a specified phase (usually $`1`$). For a function $`f(𝐱)\{0,1\}`$, where $`𝐱`$ denotes an n-bit binary string, a phase oracle implements the following operation:

$$|𝐱\stackrel{\mathrm{O}}{}(1)^{\mathrm{f}(𝐱)}|𝐱$$
(1)
The search problem is phrased as follows: Given an oracle, $`O`$, which implements the function $`f(𝐱)\{0,1\}`$, find a state in the computational basis $`𝐱`$ for which $`f(𝐱)=1`$. Assuming a perfect, noiseless oracle (the concept of a noisy oracle is introduced below), the well-known result for the complexity of Grover’s algorithm, i.e., for the quantum search problem, is $`O(\sqrt{\frac{N}{M}})`$, where $`M`$ is the number of states for which $`f(𝐱)=1`$ L.K. Grover (1996); Grover (1997); C.H. Bennett, E. Bernstein, G. Brassard and U.V. Vazirani (1997); Boyer et al. (1998). In contrast, the classical complexity of the search problem is $`O(\frac{N}{M})`$ D. E. Knuth (1973). Thus, the quantum algorithm provides a quadratic speed-up over the classical algorithm.

### II.2 Grover’s algorithm in absence of noise

The description of Grover’s algorithm in this paper follows the discussion in Nielsen and Chuang (2000). The search is made on the set of $`N=2^n`$ states represented by all n-bit binary strings $`x\{0,1\}^n`$. The states are denoted by $`|x`$. Within a discrete-time computation model, Grover’s algorithm searches for marked or “target” states via repeated application of the Grover iteration operator, $`G`$, which can be written as:

$$G=\left(2|\eta \eta |I\right)O$$
(2)
where $`|\eta =\frac{1}{\sqrt{N}}_{𝐱=0}^{N1}|𝐱`$ is the equal superposition over all states. From this point on, we will assume that $`M=1`$ to simplify the discussion. The arguments can easily be extended to the case where $`M>1`$. Let the state $`|\tau `$ be the “marked” state. Then we can write out the explicit form of the oracle as:

$$O=I+\left(e^{i\pi }1\right)|\tau \tau |$$
(3)
Assuming that our registers are initialized to $`|\eta `$, it can be shown that after $`O(\sqrt{N})`$ applications of $`G`$, the quantum computer is approximately in the marked state $`|\tau `$ L.K. Grover (1996); Grover (1997). This result is demonstrated by noting that the search algorithm can be mapped onto a two-state subspace spanned by the two basis vectors $`|1`$ and $`|2`$, where

$`|1`$ $`=`$ $`|\tau ,`$ (4)
$`|2`$ $`=`$ $`{\displaystyle \frac{|\eta |\tau \tau |\eta }{\sqrt{1\left|\tau |\eta \right|^2}}}.`$ (5)
In other words, $`|1`$ is the marked state and $`|2`$ is the equal superposition over all unmarked states.

When our initial state, $`|\eta `$, and the operator $`G`$ are rewritten in the $`|1`$, $`|2`$ basis,

$$|\eta =\left(\begin{array}{c}\sqrt{\frac{1}{N}}\\ \sqrt{\frac{N1}{N}}\end{array}\right)$$
(6)
$$G=\left(\begin{array}{cc}1\frac{2}{N}& \frac{2\sqrt{N1}}{N}\\ \frac{2\sqrt{N1}}{N}& 1\frac{2}{N}\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}(\mathrm{\Theta }/2)& \mathrm{sin}(\mathrm{\Theta }/2)\\ \mathrm{sin}(\mathrm{\Theta }/2)& \mathrm{cos}(\mathrm{\Theta }/2)\end{array}\right),$$
(7)
where

$$\mathrm{sin}(\mathrm{\Theta }/2)=\frac{2\sqrt{N1}}{N},$$
(8)
we see that the effect of repeated applications of $`G`$ is to make successive rotations by $`\mathrm{\Theta }/2`$ on state $`|\eta `$. This convention for the definition of $`\mathrm{\Theta }`$ is useful because $`\mathrm{\Theta }`$ represents the angle of rotation applied by the operator $`G`$ on the Bloch sphere defined for the Grover subspace spanned by the basis states $`|1`$ and $`|2`$. We can therefore view Grover’s algorithm as the process of rotating our initial state, $`|\eta `$, onto the target state, $`|\tau |1`$ by a discrete sequence of small rotations. For small $`\mathrm{\Theta }`$ (i.e. large $`N`$), $`\mathrm{sin}\mathrm{\Theta }\mathrm{\Theta }`$, yielding an incremental rotation angle of approximately $`\frac{4\sqrt{N1}}{N}\frac{4}{\sqrt{N}}`$. Then a rotation of $`\pi `$ radians on the Bloch sphere requires $`O(\sqrt{N})`$ applications of $`G`$. Hence, Grover’s search algorithm requires $`O(\sqrt{N})`$ calls to the oracle. The well-known lower bound for quantum search has been established as $`\mathrm{\Omega }(\sqrt{N})`$ Bennett et al. (1997). Thus, Grover’s algorithm is an optimal search.

### II.3 Grover’s algorithm in presence of oracle noise

This $`O(\sqrt{N})`$ bound on the running time of the search assumes that the oracle exactly implements the phase operation specified by Eq. (1). Instead, we could envision a more realistic oracle which implements Eq. (1) with some finite degree of precision. Specifically, we adopt here the concept of a noisy phase oracle introduced in Long et al. (1999) which implements the following operation:

$$|𝐱\stackrel{\mathrm{O}}{}\left(\mathrm{e}^{\mathrm{i}(\pi +ϵ)}\right)^{\mathrm{f}(𝐱)}|𝐱,$$
(9)
where $`ϵ`$ is a small, random phase error ($`ϵ\pi `$) with probability distribution $`p(ϵ)`$ 111For other models of oracle noise, see Long et al. (1999); Pablo-Norman and Ruiz-Altaba (1999); Hsieh et al. (2002).. We will make the assumption that the error is zero-mean, i.e. that $`ϵ=0`$. For $`p(ϵ)=\delta (ϵ)`$, the noiseless oracle is recovered. In practice, the root mean-square magnitude of $`ϵ`$ can be made small through a careful physical implementation of the oracle. However, due to the finite precision of any experimental implementation, the average magnitude of $`ϵ`$ will never be zero, i.e., $`ϵ_{rms}=\sqrt{ϵ^2}>0`$.

Given this more realistic oracle model, we now investigate whether Grover’s search algorithm is still $`O(\sqrt{N})`$. The explicit form of the noisy oracle can be written as,

$$O=I+(e^{i(\pi +ϵ)}1)|\tau \tau |$$
(10)
Then Grover’s iteration operator, $`G`$, becomes:

$$G=\left(\begin{array}{cc}(1+\frac{2}{N})e^{i(\pi +ϵ)}& \frac{2\sqrt{N1}}{N}\\ \frac{2\sqrt{N1}}{N}e^{i(\pi +ϵ)}& 1\frac{2}{N}\end{array}\right),$$
(11)
which is clearly still unitary. To see the effects of the random error, we first note that any single qubit unitary operator can be written as

$$U=\mathrm{exp}(i\alpha )R_{\widehat{n}}(\mathrm{\Phi })=\mathrm{exp}(i\alpha )\left(\mathrm{cos}\left(\frac{\mathrm{\Phi }}{2}\right)Ii\mathrm{sin}\left(\frac{\mathrm{\Phi }}{2}\right)\widehat{n}\stackrel{}{\sigma }\right),$$
(12)
where $`\stackrel{}{\sigma }`$ are the Pauli operators. Thus, the Grover’s iteration operator using the noisy oracle is equivalent, up to an overall global phase factor, to a rotation, $`R_{\widehat{n}}(\mathrm{\Phi })`$, on the Bloch sphere by some angle $`\mathrm{\Phi }`$ about some direction $`\widehat{n}`$. We can determine the value of $`\mathrm{\Phi }`$ and $`\widehat{n}`$ by using Eq. (11) and Eq. (12),

$`\mathrm{cos}\left({\displaystyle \frac{\mathrm{\Phi }}{2}}\right)`$ $`=`$ $`\left(1{\displaystyle \frac{2}{N}}\right)\mathrm{cos}\left({\displaystyle \frac{ϵ}{2}}\right)`$ (13)
$`\mathrm{sin}\left({\displaystyle \frac{\mathrm{\Phi }}{2}}\right)n_x`$ $`=`$ $`{\displaystyle \frac{2\sqrt{N1}}{N}}\mathrm{sin}\left({\displaystyle \frac{ϵ}{2}}\right)`$ (14)
$`\mathrm{sin}\left({\displaystyle \frac{\mathrm{\Phi }}{2}}\right)n_y`$ $`=`$ $`{\displaystyle \frac{2\sqrt{N1}}{N}}\mathrm{cos}\left({\displaystyle \frac{ϵ}{2}}\right)`$ (15)
$`\mathrm{sin}\left({\displaystyle \frac{\mathrm{\Phi }}{2}}\right)n_z`$ $`=`$ $`\left(1{\displaystyle \frac{2}{N}}\right)\mathrm{sin}\left({\displaystyle \frac{ϵ}{2}}\right).`$ (16)
A noiseless oracle can be recovered using these equations by setting $`ϵ`$ to zero. In that case we obtain the large $`N`$ Grover rotation angle, $`\mathrm{\Phi }\frac{4}{\sqrt{N}}`$, with $`n_x=n_z=0`$, and $`n_y=1`$. Hence, the noiseless Grover’s algorithm corresponds to a rotation on the Bloch sphere about the $`y`$-axis. The effect of a noisy oracle is to add small, random perturbations to this rotation axis, thereby changing the pure $`y`$-axis rotation to a rotation about an axis $`\widehat{n}`$ which also contains non-zero $`x`$\- and $`z`$\- components.

To analyze the running time of the noisy search algorithm, we will take the limit of large $`N`$. This is a useful assumption for our current purposes since we are interested in the complexity of the search algorithm for large $`N`$. From Eqs. (1316), we can express $`G`$ up to a global phase factor as,

$$G=\mathrm{exp}\left(i\mathrm{\Phi }\left(\left(1+\frac{2}{N}\right)\mathrm{sin}\left(\frac{ϵ}{2}\right)\sigma _z+\frac{2\sqrt{N1}}{N}\left(\mathrm{sin}\left(\frac{ϵ}{2}\right)\sigma _x\mathrm{cos}\left(\frac{ϵ}{2}\right)\sigma _y\right)\right)/\mathrm{sin}\left(\frac{\mathrm{\Phi }}{2}\right)\right),$$
(17)
where $`\mathrm{sin}(\frac{\mathrm{\Phi }}{2})=\sqrt{1(1\frac{2}{N})^2\mathrm{cos}^2(\frac{ϵ}{2})}`$. Since both $`ϵ`$ and $`N^{1/2}`$ are small parameters, we can use a double Taylor expansion of the terms in the exponent with respect to $`ϵ`$ and $`N^{1/2}`$. Excluding second-order terms and higher, we obtain

$$G=\mathrm{exp}\left(i\left(ϵ\sigma _z+\frac{4}{\sqrt{N}}\sigma _y\right)\right).$$
(18)
Again using the fact that $`ϵ`$ and $`N^{1/2}`$ are both small parameters, we can make use of the Baker-Campbell-Hausdorff formula Messiah (1961) to arrive at

$$\begin{array}{ccc}\hfill G& =& \mathrm{exp}\left(iϵ\sigma _z\right)\mathrm{exp}\left(i\frac{4}{\sqrt{N}}\sigma _y\right)+O(\frac{ϵ}{\sqrt{N}})\hfill \\ & & R_{\widehat{z}}(ϵ)R_{\widehat{y}}(\frac{4}{\sqrt{N}}).\hfill \end{array}$$
(19)
Examining Eq. (19), we observe that there are two important timescales that will be relevant to any physical implementation. The first is the timescale of the $`y`$-axis rotation, $`\sqrt{N}`$. Taking the limit as $`ϵ0`$, we see that the $`y`$-axis rotation, $`R_{\widehat{y}}`$ is the “driving force” of the search algorithm, in that it rotates the initial state, $`|\eta `$, towards the target state, $`|\tau `$. The second important timescale, $`1/ϵ`$, is the timescale of the random rotation about the z-axis of the Bloch sphere, $`R_{\widehat{z}}`$. This rotation is the effect of noise and does not advance the search process. Because $`ϵ`$ is a random variable, this rotation about $`z`$ will in general be different for each application of $`G`$. However, $`ϵ`$ does have a characteristic root-mean-square magnitude, $`ϵ_{rms}`$, which is determined by the probability distribution $`p(ϵ)`$. Then the characteristic timescale of the $`z`$-axis rotation is determined by $`1/ϵ_{rms}`$.

In general, $`ϵ_{rms}`$ and $`\sqrt{N}`$ are independent parameters: $`ϵ_{rms}`$ is the phase error that is specific to the experimental implementation, while $`N`$ is the size of the database. However, in reality, these two parameters will be linked for a given experimental setup because the size of the database (i.e., the number of states $`N=2^n`$ accessed by the $`n`$ qubits of the computer) will affect the accuracy of the oracle implementation. Note that the oracle acts on all $`n`$ qubits, (see Eq. (9)). It seems very unlikely that $`ϵ_{rms}`$ would decrease as $`N`$ increases, since a greater number of qubits generally introduces more potential for noise. In the best case scenario, it might be possible to find a physical implementation for which $`ϵ_{rms}`$ is essentially constant over a large range of values for $`N`$. To ensure that our analysis is as general as possible, we will assume here that $`ϵ_{rms}`$ scales as $`N^\delta `$, where $`\delta `$ is some constant that may take both positive and negative values, or zero. We can then determine what effect the scaling constant $`\delta `$ has on algorithm performance. In particular, we shall determine the effect of $`\delta `$ upon the complexity of the algorithm, i.e., upon the number of Grover iterations required to achieve a given probability of successful measurement of the target state $`|1`$. This will allow us to further determine the maximum allowable oracle error scaling that ensures a constant probability of success independent of the library size, for a given rate of library growth.

### II.4 Dynamics and scaling of the noisy Grover search algorithm

To analyze the dynamics of the Grover iteration, we observe that any pure state on the Bloch sphere can be represented in spherical coordinates by two angles, $`\theta `$ and $`\varphi `$, where we define $`\theta `$ to be the angle from the $`z`$-axis and $`\varphi `$ to be the azimuthal angle. Here we have chosen our coordinates such that the target state $`|1`$ is the south pole ($`\theta =\pi `$) and the state $`|2`$ is the north pole ($`\theta =0`$) of our Bloch sphere. Because the operator $`G`$ depends upon the random variable $`ϵ`$, each application of $`G`$ induces a new probability distribution over the Bloch sphere which we will denote by $`f(\theta ,\varphi )`$. In this notation, our initial state at time $`t=0`$ is given by a delta function (i.e. a pure state) centered on $`\varphi =0`$, with $`\theta =\mathrm{sin}^1(1/\sqrt{N})1/\sqrt{N}`$.

Using this notation, we can write the probability of obtaining the state $`|1`$ upon measurement after $`t`$ iterations as

$$\begin{array}{ccc}\hfill P(t)& =& 1\frac{1}{2}\mathrm{cos}\theta f_t(\theta ,\varphi )\mathrm{sin}\theta d\theta d\varphi \hfill \\ & =& \frac{1}{2}\left(1\mathrm{cos}\theta \right)\hfill \\ & =& \frac{1}{2}\left(1z\right).\hfill \end{array}$$
(20)
Let us define the complexity of our search algorithm to be the number of iterations required to obtain a success probability of $`2/3`$. We immediately see that if the probability of obtaining state $`|1`$ upon measurement after $`t`$ iterations is $`P(t)`$, then by repeating this procedure approximately $`\frac{2}{3P(t)}`$ times, we can boost the overall success probability to $`2/3`$. It should be noted that the choice of the constant $`2/3`$ is arbitrary; in general, the choice of constant will not affect the complexity of the algorithm. Thus, the overall complexity of our algorithm is $`O(t/P(t))`$. Furthermore, Eq. (20) states that in order to succeed with some desired probability $`P^{}`$, our probability distribution function $`f(\theta ,\varphi )`$ must be non-negligible when $`\theta \theta ^{}`$, where

$$\theta ^{}=\mathrm{cos}^1\left(12P^{}\right).$$
(21)
In other words, to obtain some desired probability of success, $`P^{}`$, there must be a high probability of reaching points on the Bloch sphere with polar angle greater than $`\theta ^{}`$.

Having established this terminology, we will now give a phenomenological description of the evolution of the probability distribution over the Bloch sphere. Specifically, we would like to know what dependence the magnitude and scaling of the error $`ϵ_{rms}`$ has on the maximum attainable polar angle.

We first rewrite the effects of the Grover iterate as a function of polar coordinates. In the polar coordinates defined above, the $`z`$-rotation can be written as

$`\varphi `$ $`\stackrel{\mathrm{R}_{\widehat{\mathrm{z}}}}{}`$ $`\varphi +ϵ,`$ (22)
$`\theta `$ $`\stackrel{\mathrm{R}_{\widehat{\mathrm{z}}}}{}`$ $`\theta .`$ (23)
If we momentarily neglect the $`y`$-axis rotations, our $`z`$-axis rotation dynamics correspond to a random walk on the variable $`\varphi `$ with periodic boundary conditions. To analyze the effects of the $`y`$-rotation, we then take advantage of the fact that $`1/\sqrt{N}`$ is a small quantity and expand in powers of $`1/\sqrt{N}`$ to obtain

$`\varphi `$ $`\stackrel{\mathrm{R}_{\widehat{\mathrm{y}}}}{}`$ $`\varphi \mathrm{sin}\varphi {\displaystyle \frac{\mathrm{cos}\theta }{\mathrm{sin}\theta }}{\displaystyle \frac{4}{\sqrt{N}}}+O(1/N),`$ (24)
$`\theta `$ $`\stackrel{\mathrm{R}_{\widehat{\mathrm{y}}}}{}`$ $`\theta +\mathrm{cos}\varphi {\displaystyle \frac{4}{\sqrt{N}}}+O(1/N).`$ (25)
Finally, we can write the combined effects of our noisy Grover operator $`G`$ using Eqs. (2225), as

$`\varphi `$ $`\stackrel{\mathrm{G}}{}`$ $`\varphi \mathrm{sin}\varphi {\displaystyle \frac{\mathrm{cos}\theta }{\mathrm{sin}\theta }}{\displaystyle \frac{4}{\sqrt{N}}}+ϵ,`$ (26)
$`\theta `$ $`\stackrel{\mathrm{G}}{}`$ $`\theta +\mathrm{cos}\varphi {\displaystyle \frac{4}{\sqrt{N}}},`$ (27)
where we have dropped terms of $`O(1/N)`$.

Having written the dynamics of the Grover operator in terms of polar coordinates on the Bloch sphere, we now consider the probability distribution of the quantum state over the Bloch sphere after $`T`$ Grover iterations. We will analyze the dynamics for $`\varphi 1`$, and then consider when this approximation is valid. In this regime, we can approximate the Grover operator dynamics as

$`\varphi `$ $`\stackrel{\mathrm{G}}{}`$ $`\varphi +ϵ,`$ (28)
$`\theta `$ $`\stackrel{\mathrm{G}}{}`$ $`\theta +{\displaystyle \frac{4}{\sqrt{N}}}.`$ (29)
It is evident that in this small $`\varphi `$ limit, the $`\varphi `$ dynamics are completely determined by $`z`$-axis rotation, which results in a random walk on the variable $`\varphi `$ that is controlled by the random variable $`ϵ`$. The central limit theorem tells us that regardless of the probability distribution $`p(ϵ)`$ from which the random variable $`ϵ`$ is sampled, after an adequate number of iterations the distribution of $`\varphi `$ will converge to a Gaussian with width $`ϵ_{rms}\sqrt{T}`$:

$$\varphi _{rms}ϵ_{rms}\sqrt{T},$$
(30)
where $`T`$ is the number of iterations.

Let us now consider when this approximation is valid. Clearly, the condition that $`\varphi 1`$ is satisfied for the initial state $`\varphi _0=0`$. We recall that $`ϵ_{rms}`$ is assumed to scale as $`N^\delta `$, where $`\delta `$ is some constant. The validity of our small $`\varphi `$ approximation is found to be highly dependent on the scaling exponent $`\delta `$. We consider two cases: $`\delta >1/4`$ and $`\delta 1/4`$.

If $`\delta >1/4`$, we examine the probability distribution after $`T=\lambda \sqrt{N}`$ iterations, where $`\lambda `$ is some small constant. Inspection of Eq. (30) shows that the probability distribution of $`\varphi `$ after $`T`$ steps will then have width $`\sqrt{\lambda }N^{1/4\delta }`$. Since $`N`$ can certainly be made arbitrarily large, this justifies our assumption that $`\varphi `$ is small. Within this regime, Eq. (28) shows that the dynamics of the variable $`\theta `$ are simply those of the deterministic $`y`$-rotation with constant increment $`\frac{4}{\sqrt{N}}`$. After $`T`$ iterations of this, we obtain

$$\begin{array}{ccc}\hfill \theta _T& =& \theta _0+\frac{4T}{\sqrt{N}},\hfill \\ & & \frac{4T}{\sqrt{N}},\hfill \end{array}$$
(31)
which is the same as the large $`N`$ limit of the noiseless Grover search, Eq. (9). Thus when $`\delta >1/4`$, the error has no effect on algorithmic complexity for significantly large $`N`$, yielding a complexity of $`O(N^{1/2})`$. This bound is trivially tight since the search problem is well-known to have a lower bound of $`\mathrm{\Omega }(N^{1/2})`$Bennett et al. (1997). For the remainder of this section, we will therefore analyze the complementary case, i.e., when $`ϵ_{rms}`$ scales as $`N^\delta `$ with $`\delta 1/4`$, and determine the effect of the error on algorithmic complexity. We note that this latter case includes the best-case physical situation of $`ϵ_{rms}`$ independent of $`N`$, i.e., $`\delta =0`$.

In order to analyze the system dynamics for $`\delta 1/4`$, let us examine the probability distribution after $`T`$ iterations such that

$$T=\lambda /ϵ_{rms}^2,$$
(32)
where $`\lambda `$ is a small constant. By Eq. (30), the distribution width $`\varphi _{rms}`$ is proportional to $`\sqrt{\lambda }`$. Thus, by selecting a small enough constant $`\lambda `$, our approximation that $`\varphi `$ is small is again justified. Using Eq. (28) and Eq. (31), we find that the $`\theta `$-rotation is again essentially deterministic and that

$`\theta _T`$ $`=`$ $`{\displaystyle \frac{4T}{\sqrt{N}}},`$ (33)
$`=`$ $`{\displaystyle \frac{4\lambda }{ϵ_{rms}^2\sqrt{N}}}`$
$`=`$ $`O\left({\displaystyle \frac{1}{ϵ_{rms}^2\sqrt{N}}}\right).`$
Using Eq. (20), we obtain

$$\begin{array}{ccc}\hfill P(T)& =& \frac{1}{2}\left(1\mathrm{cos}\left(\frac{\lambda }{ϵ_{rms}^2\sqrt{N}}\right)\right)\hfill \\ & & \frac{\lambda ^2}{4ϵ_{rms}^4N},\hfill \end{array}$$
(34)
where we have used the fact that $`\frac{\lambda }{ϵ_{rms}^2\sqrt{N}}1`$ when $`\delta 1/4`$. By proving that we can attain (essentially deterministically, as described above) a polar rotation of *at least* $`\frac{4\lambda }{ϵ_{rms}^2\sqrt{N}}`$ in $`T=\lambda /ϵ_{rms}^2`$ iterations, we have shown that the complexity of the noisy search algorithm is $`O(\frac{T}{P(T)})=O(ϵ_{rms}^2N)`$. Taking the scaling of $`ϵ_{rms}`$ into account ($`ϵ_{rms}N^\delta `$), leads to the overall algorithmic complexity $`O(N^{12\delta })`$ 222It should be noted that this analysis applies only in the limit of large $`N`$, which is what we mean when we talk of algorithmic complexity in the first place.. So for $`\delta =1/4`$ we obtain an optimal solution having the same speed-up as the noiseless quantum search. For $`0<\delta <1/4`$, we find a scaling intermediate between the noiseless quantum and classical search algorithms. The quantum speed-up factor decreases as $`\delta `$ approaches zero and is completely lost when $`\delta =0`$. (Note that the quantum search formally becomes slower than the classical search in the worst case scenario when the phase errors increase with $`N`$, i.e., $`\delta <0`$.)

In order to show that this complexity bound is tight, we now will show that we can attain *at most* a polar rotation of $`O\left(\frac{1}{ϵ_{rms}^2\sqrt{N}}\right)`$ in sub-classical time. Repeating the argument above following Eq. (33) will then lead to the identification of $`\mathrm{\Omega }(N^{12\delta })`$ as a lower bound. Thus, we find that $`O(N^{12\delta })`$ is a tight bound on the complexity in presence of noise. Obtaining a lower bound on the algorithmic complexity is important because without it, it would be unclear whether we could further reduce our search complexity by running the noisy search algorithm for more than $`O(1/ϵ_{rms}^2)`$ iterations. The lower bound will demonstrate that after $`O(1/ϵ_{rms}^2)`$ iterations, it is impossible to achieve a super-classical enhancement in success probability.

To achieve a polar rotation greater than $`O\left(\frac{1}{ϵ_{rms}^2\sqrt{N}}\right)`$, we will need to obtain a significant probability density on the region of the Bloch sphere described by $`\theta \frac{c}{ϵ_{rms}^2\sqrt{N}}`$, where $`c`$ is a constant. Let us therefore consider the dynamics of a state on this region of the Bloch sphere. Given a pure state with $`\varphi _0=0`$ and $`\frac{\pi }{2}\theta _0\frac{c}{ϵ_{rms}^2\sqrt{N}}`$, we examine the action of $`T`$ iterations of the Grover algorithm on this state, where $`T=\lambda /ϵ_{rms}^2`$. Now, instead of choosing $`\lambda `$ small as previously, we select $`\lambda 1`$, which ensures a large width for the long time distribution of the azimuthal angle $`\varphi `$, Eq. (30). Thus after $`T`$ iterations with this large width, we obtain a completely uniform distribution on the variable $`\varphi `$. We now consider what happens to the polar variable, $`\theta `$. Referring to Eq. (23) and Eq. (25), we find that after $`T`$ iterations we obtain

$$\begin{array}{ccc}\hfill \left|\theta _T\theta _0\right|& & \frac{T}{\sqrt{N}}\hfill \\ & & \frac{\lambda }{ϵ_{rms}^2\sqrt{N}}\hfill \\ & =& O\left(N^{1/2+2\delta }\right).\hfill \end{array}$$
(35)
Again, provided that $`\delta 1/4`$, this difference can be made arbitrarily small by selecting a large enough $`N`$. Thus, after $`T`$ iterations, our probability distribution of the polar coordinate $`\theta `$ is confined within an arbitrarily small window around $`\theta _0`$, but is uniformly distributed on the azimuthal coordinate $`\varphi `$. Figure 1 shows a graphical representation of these dynamics.

Let us now consider the effects when the Grover operator is subsequently applied to such a probability distribution that is uniform in $`\varphi `$. We will show that there exists a symmetry in the associated $`\theta `$ transformations in this regime, such that one half of all azimuthal angles, $`\varphi `$, are associated with a $`\theta `$ rotation in one direction, and the other half are associated with $`\theta `$ rotations of equal magnitude in the opposite direction. To reveal this symmetry, we show that the transformation $`\varphi \varphi +\pi `$ on Eqs. (26) and (27) results in an exactly opposite rotation in the $`\theta `$ direction to that associated with the azimuthal angle, $`\varphi `$. From Eq. (35), we know that the distribution of $`\theta _T`$ will be arbitrarily close to $`\theta _0`$. Since we have assumed that our initial state was located in the region $`\pi /2\theta _0\frac{c}{ϵ_{rms}^2\sqrt{N}}`$, we can rewrite Eq. (26) as follows:

$`\varphi `$ $`\stackrel{\mathrm{G}}{}`$ $`\varphi \mathrm{sin}\varphi {\displaystyle \frac{\mathrm{cos}\theta }{\mathrm{sin}\theta }}{\displaystyle \frac{1}{\sqrt{N}}}+ϵ`$ (36)
$``$ $`\varphi \mathrm{sin}\varphi {\displaystyle \frac{ϵ_{rms}^2\sqrt{N}}{c}}{\displaystyle \frac{1}{\sqrt{N}}}+ϵ`$
$`=`$ $`\varphi \mathrm{sin}\varphi {\displaystyle \frac{ϵ_{rms}^2}{c}}+ϵ.`$
We now make use of the arbitrariness in choice of the constant $`c`$, choosing $`c1`$ so that $`c\mathrm{sin}\varphi `$, making the second term arbitrarily small. Hence we obtain

$`\varphi `$ $`\stackrel{\mathrm{G}}{}`$ $`\varphi +ϵ,`$ (37)
which demonstrates that the transformation $`\varphi \varphi +\pi `$ has no net effect on the azimuthal dynamics. On the other hand, our $`\theta `$ dynamics become

$`\theta `$ $`\stackrel{\mathrm{G}}{}`$ $`\theta +\mathrm{cos}(\varphi +\pi ){\displaystyle \frac{4}{\sqrt{N}}}+O(1/N)`$ (38)
$`=`$ $`\theta \mathrm{cos}\varphi {\displaystyle \frac{4}{\sqrt{N}}}+O(1/N).`$
Eq. (38) shows that the transformation $`\varphi \varphi +\pi `$ causes the polar angle $`\theta `$ to be rotated in the opposite direction, as a consequence of the change in sign of the $`\mathrm{cos}\varphi `$ factor. It follows that once $`\varphi `$ has reached a uniform distribution, there will be an equal probability of rotating $`\theta `$ by $`4/\sqrt{N}`$ in the positive and negative directions. This behavior is valid on the entire region of the Bloch sphere specified by our initial condition $`\pi /2\theta _0\frac{c}{ϵ_{rms}^2\sqrt{N}}`$. As a result, the $`\theta `$ dynamics over this entire region no longer resemble a deterministic rotation towards $`\theta =\pi `$, but instead resemble a random walk on $`\theta `$ with step size of approximately $`4/\sqrt{N}`$. Figure 1 shows a graphical representation of these dynamics.

Using the central limit theorem, we find that

$$\theta _{rms}\sqrt{t/N},$$
(39)
where $`t`$ is the number of iterations after the uniform distribution in $`\varphi `$ has been attained, i.e., after the first $`T`$ iterations. Consequently, to “move” the polar angle of the probability distribution by some small angle $`\mathrm{\Delta }\theta `$ now takes time proportional to $`\mathrm{\Delta }\theta ^2N`$. From Eq. (20), the corresponding change in the probability of a successful measurement of the target state is

$$\mathrm{\Delta }P\frac{1}{2}\mathrm{sin}\mathrm{\Delta }\theta \mathrm{\Delta }\theta \frac{1}{2}\mathrm{\Delta }\theta ^2\mathrm{\Delta }t/N$$
(40)
Thus, once we have entered the region of the Bloch sphere characterized by $`\pi /2\theta \frac{c}{ϵ_{rms}^2\sqrt{N}}`$, the probability of a successful measurement further increases only at the rate $`t/N`$. This is the same as the classical result and consequently, we can attain no further speed-up over the classical algorithm once we have entered this regime. The polar angle $`O(\frac{1}{ϵ_{rms}^2\sqrt{N}})`$ thus constitutes an upper bound on the rotation that can be achieved in sub-classical time. Hence, by the arguments given above, our bound $`O(N^{12\delta })`$ on the complexity of the noisy quantum search algorithm is a tight bound.

### II.5 Summary of discrete time noise Grover search

Using the above results, we can now give an accurate phenomenological description of the dynamics of the noisy search algorithm. We identify the critical timescale of the algorithm as $`t_{mixing}=1/ϵ_{rms}^2`$, which is the mixing time of the $`z`$-rotation induced by the random error. For times $`tt_{mixing}`$, the algorithm proceeds more or less deterministically, and the initial state is rotated towards the target state $`|1`$, attaining a polar angle of $`\theta =\frac{1}{ϵ_{rms}^2\sqrt{N}}`$. However, for times $`tt_{mixing}`$, the $`\varphi `$ variable becomes completely mixed and the $`\theta `$-rotation, which is the driving force of the search process, has an equal probability of increasing or decreasing the polar angle of the probability distribution. As a result, the $`\theta `$ dynamics also become those of a random walk, and the searching rate becomes classical. Therefore, in order to gain the maximum speed-up over the classical algorithm, the best methodology we can employ is to run the algorithm for $`tt_{mixing}`$, measure, and then repeat this process to boost our overall success probability. This protocol will allow the sub-classical scaling to be retained.

The effects of the scaling of $`ϵ_{rms}`$ on the complexity of the search algorithm are seen to be different according to the value of the error scaling index $`\delta `$. We saw that if $`\delta >1/4`$, the search algorithm was essentially unaffected by the presence of error and the complexity is identical to the noiseless quantum search result of $`O(\sqrt{N})`$. On the other hand, if $`\delta 1/4`$, then the optimal speed-up is obtained by running the algorithm for $`t=O(1/ϵ_{rms}^2)`$ iterations, achieving a maximum polar angle of $`\theta =O(1/ϵ_{rms}^2\sqrt{N})`$ and then measuring. Expanding Eq. (20) yields the measurement success rate, $`P(t)=O(1/ϵ_{rms}^4N)`$ and a resulting algorithmic complexity of $`O(t/P(t))=O(ϵ_{rms}^2N)=O(N^{12\delta })`$. So the optimal error scaling is given by $`\delta 1/4`$, i.e., the errors scale as $`O(N^{1/4})`$. In contrast, the most physically realistic constant error scaling $`\delta =0`$, corresponding to a constant error over a range of $`N`$ values, yields no speed up over the classical search.

It is interesting to compare these results to the analytic results for a constant phase error oracle given in Long et al. (1999). In the constant phase error model, the oracle applies the same phase $`\mathrm{exp}(i(\pi +ϵ))`$ to the marked state at each iteration. It can be shown that for a constant phase error of magnitude $`ϵ`$, the error must scale as $`ϵ=O(N^{1/2})`$ in order to obtain a quadratic speed-up over the classical search. The differences in error accumulation between the constant and random phase error processes can be compared (for example) to the ballistic and diffusive regimes of Brownian motion, respectively. Constant phase errors tend to accumulate quickly since they are all in the same “direction”. On the other hand, subsequent random phase errors can cancel each other out and thus accumulate more slowly.

Let us now consider a question of physical importance. Given an oracle with error magnitude $`ϵ_{rms}`$ and a library of size $`N`$, let us assume that we can attain some constant success rate $`P`$ after $`T`$ iterations. Given a larger library of size $`N^{}=kN`$, we wish to obtain the same success rate $`P`$ after $`T^{}=\sqrt{k}T`$ iterations. In other words, we wish to obtain the ideal quantum search complexity of $`O\left(\sqrt{N^{}}\right)`$ also for a larger library size. What is the maximum allowable oracle error, $`ϵ_{rms}^{}`$, for the larger library? The answer to this question follows immediately from our previous analysis. We know that $`ϵ_{rms}`$ can scale as $`N^\delta `$ with $`\delta 1/4`$ without affecting algorithmic complexity. Thus the maximum allowable error in the oracle for the larger library is $`ϵ_{rms}^{}=k^{1/4}ϵ_{rms}`$. As an example, if the library size is doubled, the oracle error for the larger library is required to be at most 0.84 the corresponding error for the smaller library.

To verify the validity of our analytic results, we can simulate the effects of a noisy oracle numerically. Given a library of size $`N`$, we simulated the Grover’s search algorithm using the procedure described in Section II.3. For each iteration of the search algorithm, a phase error $`ϵ`$ was selected from a zero-mean Gaussian distribution with standard deviation $`ϵ_{rms}`$. The algorithm was run for $`t=\pi \sqrt{N}/4`$ iterations and the maximal probability of success attained was recorded. This process was repeated $`100`$ times and the average success probability was calculated. Figure 2 plots the average probability of success, $`P`$, versus the library size, $`n=\mathrm{log}N`$, for various values of $`ϵ_{rms}`$. Figure 3 shows the value of $`\mathrm{log}ϵ_{rms}`$ that yields a constant success probability of $`P=1/2`$ for a given library size. As predicted by our analysis above, we find that in order to obtain a constant success rate as a function of $`n`$ (or $`N`$), the error must scale as $`ϵ_{rms}N^{1/4}`$.

## III Effect of Noise on Continuous Time Analogue of Quantum Search

In this section, we analyze the effect of noise on a slightly different model for quantum search that has been proposed by Farhi and Gutmann Farhi and Gutmann (1998). In this model, one again starts in the symmetric superposition of all states, $`|\eta `$, and then applies the following Hamiltonian

$$H_0=|\eta \eta |+|\tau \tau |=H_\eta +H_\tau ,$$
(41)
where $`|\tau `$ is the marked state. Note that this Hamiltonian is directly related to the Grover’s iterate, Eqs. (2) and (3). Action of Eq. (2) amounts to simply applying the Hamiltonian $`H_\eta `$ for a time $`\pi `$, followed by applying the Hamiltonian $`H_\tau `$ for a time $`\pi .`$ Clearly the two operators Eq. (2) and Eq. (41) would be equivalent if $`H_\eta `$ and $`H_\tau `$ commuted. However, $`[H_\eta ,H_\tau ]=\frac{1}{\sqrt{N}}(|\eta \tau |+|\tau \eta |)`$. Thus the two methods become similar for large $`N`$ and are formally equivalent as $`N`$ approaches infinity.

Farhi and Gutmann calculated the time evolution of the system when one starts in the equal superposition state $`|\eta `$. They found that

$$e^{iHt}|\eta =e^{it}\left(\left(\frac{1}{\sqrt{N}}\mathrm{cos}\left(\frac{t}{\sqrt{N}}\right)i\mathrm{sin}\left(\frac{t}{\sqrt{N}}\right)\right)|\tau +\sqrt{\left(1\frac{1}{N}\right)}\mathrm{cos}\left(\frac{t}{\sqrt{N}}\right)|\eta \right).$$
(42)
They noticed that at time $`t=\frac{\pi \sqrt{N}}{2}`$, the initial state, $`|\eta `$, has evolved to the marked state, $`|\tau =|1`$. The time required to evolve to the marked state scales as $`O(\sqrt{N})`$, matching the complexity of Grover’s algorithm with respect to an oracle. As a result, we take the time that it takes to reach the marked state as a function of $`N`$ to be the measure of the complexity of a continuous time algorithm.

As noted above, in the limit of large $`N`$ the continuous time and discrete Grover’s algorithm are formally equivalent. Therefore, it seems useful to also evaluate the effect of a fluctuating ”oracle” in the continuous time picture. The continuous time Hamiltonian that models a discrete quantum search noisy oracle with phase fluctuation $`ϵ`$ is given by

$`H`$ $`=`$ $`H_\eta +(1+\xi )H_\tau `$
$`=`$ $`|\eta \eta |+(1+\xi )|\tau \tau |.`$
Here $`\xi `$ is a time-dependent, Markovian stochastic variable that satisfies $`_0^\pi \xi 𝑑t=ϵ`$. We shall assume that $`ϵ`$ fluctuates and that it can be described by a Gaussian distribution.

In order to evaluate the Hamiltonian and its effect on the initial state, it is simpler to transform to our orthonormal two state basis $`|1`$ and $`|2`$ defined by Eqs. (4) and (5). One can then decompose the transformed Hamiltonian into the corresponding spin operators to find that

$$H=\left(1+\frac{\xi }{2}\right)𝐈+\left(\frac{1}{N}+\frac{\xi }{2}\right)\sigma _z+\frac{1}{\sqrt{N}}\sqrt{\left(1\frac{1}{N}\right)}\sigma _x.$$
(44)
When the fluctuations in the time evolution operator $`e^{iHt}`$ are now considered, i.e., the resulting phase $`ϵ0`$, it is useful to approach the problem by examining the evolution of the density matrix. Consider an initial density matrix,

$$\rho (0)=\left[\begin{array}{cc}a& b\\ b^{}& d\end{array}\right].$$
(45)
It is well established that a fluctuating $`\sigma _z`$ component of the Hamiltonian leads to dephasing. For example, applying the fluctuating perturbation $`\frac{\xi }{2}\sigma _z`$ for a time $`\pi `$ yields,

$$\rho (t)=\mathrm{exp}\left(i\sigma _z\pi /2\right)\rho (0)\mathrm{exp}\left(i\sigma _z\pi /2\right).$$
(46)
Subsequently averaging over all possible values of the resulting total phase, $`ϵ`$, leads to the following evolution:

$$\frac{1}{\sqrt{2\pi }s}e^{ϵ^2/2s^2}e^{iϵ\sigma _z/2}\left[\begin{array}{cc}a& b\\ b^{}& d\end{array}\right]e^{iϵ^2/\sigma _z}𝑑ϵ=\left[\begin{array}{cc}a& be^{s^2/2}\\ b^{}e^{s^2/2}& d\end{array}\right].$$
(47)
Here $`s^2`$ is the variance of the fluctuating phase $`ϵ`$. Since it has been assumed that the fluctuating field $`\xi `$ is Markovian, one can treat this as a dephasing term described by a decay constant $`\mathrm{\Gamma }=s^2/2\pi `$ and then use the corresponding Bloch equations K. Blum (1981) to calculate the evolution of the system. This Markovian approximation is valid when the fluctuation of $`\xi `$ is much faster than $`\frac{1}{\sqrt{N}}`$. To analyze the system evolution, we therefore decompose our density matrix using the Bloch representation,

$$\rho =\frac{1}{2}𝐈+\frac{1}{2}\widehat{n}\sigma .$$
(48)
The time evolution is calculated by solving the Bloch equations for the components of $`\widehat{n}`$ K. Blum (1981):

$`\underset{x}{\overset{.}{n}}`$ $`=`$ $`{\displaystyle \frac{2}{N}}n_y\mathrm{\Gamma }n_x,`$ (49)
$`\underset{y}{\overset{.}{n}}`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{N}}}\sqrt{\left(1{\displaystyle \frac{1}{N}}\right)}n_z{\displaystyle \frac{2}{N}}n_x\mathrm{\Gamma }n_y,`$ (50)
$`\underset{z}{\overset{.}{n}}`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{N}}}\sqrt{\left(1{\displaystyle \frac{1}{N}}\right)}n_y.`$ (51)
In order to understand the effect of the dephasing on the algorithmic complexity, we examine these equations in the limit of large $`N`$ and keep only terms that are of order $`\frac{1}{\sqrt{N}}.`$ This yields only two coupled equations, in $`n_y`$ and in $`n_z`$:

$`\underset{y}{\overset{.}{n}}`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{N}}}n_z\mathrm{\Gamma }n_y`$ (52)
$`\underset{z}{\overset{.}{n}}`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{N}}}n_y.`$ (53)
Our initial density matrix is given by

$`\rho `$ $`=`$ $`|\eta \eta |`$ (54)
$`=`$ $`{\displaystyle \frac{1}{2}}𝐈+{\displaystyle \frac{1}{2}}\left(1+{\displaystyle \frac{2}{N}}\right)\sigma _z+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{2\sqrt{N1}}{N}}\right)\sigma _x.`$
We note that the quantity $`n_z`$ is the projection onto the state $`|1`$ and thus provides a measure of how well the computation is proceeding. Initially we have $`n_z=1+\frac{2}{N}`$, while our computation becomes complete when $`n_z=1.`$ Solving the differential equations in the large $`N`$ limit, (Eq. (53)), with initial conditions $`n_y(0)=0`$, $`n_z(0)=1+\frac{2}{N}`$, yields the following solution:

$$\begin{array}{cccc}n_y(t)\hfill & =\hfill & \frac{2(1+\frac{2}{N})}{\sqrt{\mathrm{\Gamma }^2N16}}\hfill & \left(\mathrm{exp}\left(\frac{t}{2}\left(\sqrt{\mathrm{\Gamma }^2}\sqrt{\mathrm{\Gamma }^216/N}\right)\right)+\mathrm{exp}\left(\frac{t}{2}\left(\sqrt{\mathrm{\Gamma }^2}\sqrt{\mathrm{\Gamma }^216/N}\right)\right)\right)\hfill \\ n_z(t)\hfill & =\hfill & \frac{(1+\frac{2}{N})}{2\sqrt{\mathrm{\Gamma }^2N16}}\hfill & ((\sqrt{\mathrm{\Gamma }^2N}+\sqrt{\mathrm{\Gamma }^2N16})\mathrm{exp}\left(\frac{t}{2}(\sqrt{\mathrm{\Gamma }^2}\sqrt{\mathrm{\Gamma }^216/N})\right)\hfill \\ & & & +(\sqrt{\mathrm{\Gamma }^2N}+\sqrt{\mathrm{\Gamma }^2N16})\mathrm{exp}\left(\frac{t}{2}(\sqrt{\mathrm{\Gamma }^2}+\sqrt{\mathrm{\Gamma }^216/N})\right)).\hfill \end{array}$$
(55)
Ideally, we would like to find the time at which $`n_z(t)=1`$, or equivalently, the time when the probability of reaching the marked state is unity,

$$P(t)=\tau |\rho (t)|\tau =\frac{1}{2}\left(1+n_z(t)\right)=1.$$
(56)
When measuring the complexity of Grover’s algorithm however, we need only to find the time required such that the probability $`P(t)`$ of being in the marked state $`\tau `$ is greater than some constant. For concreteness, we choose here the minimum time satisfying $`P(t)=\tau |\rho (t)|\tau 1/4`$.

To determine this time explicitly, we calculate $`P(t)`$ as a function of various values of $`\mathrm{\Gamma }`$ in the limit of large $`N`$. We observe that in Eq. (55) the term $`\sqrt{\mathrm{\Gamma }^216/N}`$ is imaginary when $`\mathrm{\Gamma }<\frac{4}{\sqrt{N}}`$. Consequently, we expect that the behavior will be drastically different for the two regimes a) $`\mathrm{\Gamma }<\frac{4}{\sqrt{N}}`$, and b) $`\mathrm{\Gamma }>\frac{4}{\sqrt{N}}`$. In order to make a direct comparison between the continuous time behavior and the discrete time results, we choose the dephasing constant $`\mathrm{\Gamma }`$ to scale with $`N`$ in the same way as in Section II, namely

$$\mathrm{\Gamma }=\alpha N^{2\delta }ϵ_{rms}^2,$$
(57)
where $`\alpha `$ is a proportionality constant. This provides the contact point of the continuous time search algorithm with the discrete time algorithm of Section II. We then see that regimes a) and b) correspond to the two regimes already established in Section II, i.e., in region a), $`\delta 1/4`$, and in region b), $`\delta 1/4`$. We note that our regions b) and a) both include $`\delta =1/4`$. The boundary line between the regions occurs when $`\delta =1/4`$ and $`\alpha =4`$.

For $`\mathrm{\Gamma }=\alpha N^{2\delta }<\frac{4}{\sqrt{N}}`$, i.e., in region a), we calculate from Eq. (55) that

$$P(t)=\frac{1}{2}+\frac{1}{2}\left(1+\frac{2}{N}\right)e^{t/2\mathrm{\Gamma }}\left(\mathrm{cos}(\sqrt{16/N\mathrm{\Gamma }^2}t/2)+i\sqrt{\frac{\mathrm{\Gamma }^2N}{\mathrm{\Gamma }^2N16}}\mathrm{sin}(\sqrt{16/N\mathrm{\Gamma }^2}t/2)\right).$$
(58)
We now pick a time

$$t^{}=\frac{2\pi }{\sqrt{16/N\mathrm{\Gamma }^2}}$$
(59)
and find that

$`P(t^{})`$ $`=`$ $`{\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{2}}\left(1{\displaystyle \frac{2}{N}}\right)\mathrm{exp}\left({\displaystyle \frac{\pi \mathrm{\Gamma }}{\sqrt{16/N\mathrm{\Gamma }^2}}}\right).`$ (60)
Hence,

$`P(t^{})`$ $`>`$ $`{\displaystyle \frac{1}{2}}O(1/N).`$ (61)
Since this is larger than the value corresponding to our definition of minimum time, it implies that the complexity of the search algorithm is bounded from above by $`O(t^{})`$. Inspection of Eq. (59) shows that $`t^{}`$ is an increasing function of $`\mathrm{\Gamma }`$. Therefore, in order to evaluate an upper bound for $`t^{}`$ in the regime a) where $`\delta 1/4`$, we need to evaluate $`t^{}`$ for the largest possible value of $`\mathrm{\Gamma }`$. The largest value of $`\mathrm{\Gamma }`$ in this regime lies on the boundary with regime b), namely where $`\delta =1/4`$ and $`\alpha =4`$. Hence, we evaluate $`t^{}`$ as one asymptotically approaches the boundary between regions a) and b), for a given error scaling $`\delta `$. We choose $`\delta =1/4`$ and $`\alpha =4m`$, where m is a small constant greater than zero. Substituting in Eq. (60) yields $`t^{}=\frac{2\pi \sqrt{N}}{m}`$. Thus the minimum time will be on the order of $`O(\sqrt{N})`$, corresponding to an upper bound on the algorithmic complexity of $`O(\sqrt{N})`$. Hence for $`\delta 1/4`$, in regime a), the continuous time search algorithm achieves its maximal algorithmic quantum speed up. This is in agreement with the results from Section II, which showed that the discrete time quantum search algorithm yields the maximum algorithmic speedup when the root mean square error $`ϵ_{rms}`$ is smaller than $`N^{1/4}`$, corresponding to $`\delta 1/4`$.

To complete the continuous time analysis, we solve for $`P(t)`$ in the regime b) where $`\delta 1/4`$. Here evaluation of $`n_z`$ in the limit of large $`N`$ leads to

$`P(t)={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2}}e^{4t/N\mathrm{\Gamma }}.`$ (62)
Setting $`P(t^{})=1/4`$ and solving for $`t^{}`$ yields

$`t^{}`$ $`=`$ $`{\displaystyle \frac{N\mathrm{\Gamma }\mathrm{ln}(2)}{4}}`$ (63)
$`=`$ $`O(N^{12\delta }).`$ (64)
We have verified our conclusions by numerically simulating Eq. (55) for various values of $`N`$ and $`\mathrm{\Gamma }`$. Figure 4 shows a logarithmic plot of the minimum time $`t^{}`$ required to obtain a success probability of $`P(t^{})=1/4`$, as a function of the error scaling parameter $`\delta `$, for various values of $`N`$. It is evident that the algorithmic complexity shows a marked transition at $`\delta =1/4`$, from scaling $`O(N^{1/2})`$ for $`\delta 1/4`$, to $`O(N^{12\delta })`$ for $`\delta 1/4`$, as predicted by Eq. (59) and Eq. (64) respectively. Equivalently, we can state that the continuous time quantum search algorithm with randomized phase error achieves minimal complexity and hence maximum algorithmic speedup when $`ϵ_{rms}N^{1/4}`$ ($`\delta 1/4`$). This agrees with the results for the discrete-time algorithm.

## IV Conclusions

The analysis in this paper has provided a phenomonological description of the dependence of the algorithmic complexity of Grover’s algorithm on the scaling of oracle phase error for a discrete quantum search, and on stochastic Hamiltonian errors leading to phase error in a continuous time quantum search algorithm. In both the discrete- and continuous-time versions of the algorithm, it was found that if the phase error scaled with size as $`N^\delta `$, then for $`\delta 1/4`$ the effect on the complexity of the algorithm for large $`N`$ was negligible. However, if the size scaling of the error lies in the regime $`\delta 1/4`$, then for large $`N`$ it was determined that there is tight bound of $`O\left(N^{12\delta }\right)`$ on the complexity of the search algorithm. In particular, this implies that in the presence of any constant (non-zero) amount of phase error in the oracle ($`\delta =0`$), there exists some library size $`N_{max}`$ above which the quantum search algorithm no longer provides a quadratic speed-up. In this case, for databases of size $`N>N_{max}`$, the search time is $`\mathrm{\Omega }(N)`$, which is equivalent to the classical result and there is therefore no quantum speedup. Intermediate error scaling, $`0<\delta <1/4`$, provides speedup intermediate between the classical and quantum limits, respectively. These results hold for both the discrete-time and continuous-time quantum search algorithms, and assume very little about the specific form of the underlying error processes.

The complexity analysis we have made here is also important for determining the precision needed in scaling up a quantum search. For instance, let us assume that we are able to implement a quantum search for a library of size $`N`$ with an oracle error of magnitude $`ϵ_{rms}`$. Then to perform a quantum search on a library of size $`kN`$ with equivalent accuracy, our results imply that we need to implement an oracle with an error of at most $`ϵ_{rms}/k^{1/4}`$. Since this must lie in the regime $`\delta 1/4`$, physically, this requires a system where the phase error decreases as a function of database size. Consequently, the precision must increase exponentially as a function of the number of qubits, putting severe demands on the physical realization. In contrast, if the error and hence the precision is constant in the system size, (e.g., a system where the natural line width is independent of the number of states), then there always exists a database size such that the quantum approach offers no speed up over the classical search algorithm.

The main consequence of the non-robustness of these forms of quantum searches to oracle noise that was demonstrated here is a practical limitation on the size of the library on which may be searched with a quadratic speed-up using quantum search algorithms without any explicit error correction. This result has significant consequences for physical implementation of quantum search algorithms, since although quantum error correction can be used to reduce the error present in the oracle, such error correction procedures can require significant resources Preskill (1998). In practice it will therefore be necessary to balance the cost of error correction (in both spatial and temporal resources) with the extent of speed-up attained by a noisy quantum search. The analytic results presented in this paper provide a useful bound on the maximum oracle error permissible if a quadratic speedup is to be retained. Above this maximum allowable error, we must use error correction in any physical implementation. Conversely, below this maximum error, we can be confident that error correction will not be necessary, provided that oracle phase error is the primary source of error.

As a final comment we point out that although Grover’s algorithm and its continuous time analogue are not inherently robust to phase error in the oracle, it is not clear whether other implementations of quantum search may be inherently robust. Exploration of both active and passive error correction schemes for Grover’s algorithm will therefore be a valuable direction for future work.

###### Acknowledgements.

NS thanks the University of California, Berkeley, for a Berkeleyan Fellowship. The work of KRB was supported by the Fannie and John Hertz Foundation. KBW thanks the Miller Foundation for Basic Research for a Miller Research Professorship 2002-2003. This effort is sponsored by the Defense Advanced Research Projects Agency (DARPA) and the Air Force Laboratory, Air Force Material Command, USAF, under agreement number F30602-01-2- 0524. We also thank NSF ITR/SY award 0121555.