text:
# Nondiffusive heat transfer in muscle tissue. Preliminary results

## I Introduction

During the last two decades the description of heat transfer in living tissue (bioheat transfer) and mathematical modeling of temperature distribution in it under local strong heating have attracted attention due to several reasons. First, this problem is an important part of planning and optimizing hyperthermia treatment and thermotherapy of tumors where the precise control over the tissue temperature is essential (see, e.g., (I1, ; I2, ; I3, )). Second, what is interesting from the standpoint of fundamental physics, living tissue is a complex medium in which highly branching blood vessels form a fractal network of fast heat transport I4 , endowing heat propagation in it with anomalous properties BLG . Third, living tissue is an active medium, in particular, it responds to strong heating, which causes the blood perfusion rate to grow locally by tenfold Song84 . In this case one meets also the basic problem of how a natural hierarchical system comprising a huge amount of elements can respond properly to changes in its state when none of the elements possesses the total information required to govern the system as a whole LGenv95 ; gaf1 .

By now a number of significant steps towards developing the bioheat transfer theory have been made. In particular, Chen & Holmes CH80  actually stated the microscopic description of the problem. Baish et al. Bi86 , Lagendijk et al. LD1 ; LD2 , Roemer et al. RD1 ; RD2 , and Gafiychuk & Lubashevsky BLG  studied the effect of the vessel discreteness on the temperature field. Weinbaum et al. WJ85 ; WXZE97  and Lubashevsky & Gafiychuk BLG  took into account sophisticated details of heat interaction between different vessels, namely, the counter-current effect, Chato C80  and Weinbaum & Jiji WJ85 , that renormalizes the kinetic coefficients of heat propagation in the tissue. Lubashevsky & Gafiychuk BLG ; LGenv95 ; gaf1  proposed a hierarchical model for the tissue response to local heating. There is also an approach by Lagendijk et al. (see, e.g., Lag  and references therein) that dealing with individual vessels tackles the bioheat transfer problem numerically.

All these models and approaches are based on the assumption that the cellular tissue is a passive medium and heat propagation in it obeys Fick’s law, i.e. the heat flux $`𝐐`$ in the cellular tissue is proportional to the temperature gradient $`T`$. Namely, $`𝐐=\kappa T`$ and the thermal conductivity $`\kappa `$ is isotropic or, at least, is slightly anisotropic. Indeed, according to the experimental data summarized in GMLM95 , for example, the value of the temperature diffusivity $`D=\kappa /(c\rho )`$ ($`c`$ and $`\rho `$ are the specific heat and density of the tissue) is approximately the same for different tissues, $`D10^31.8`$, 1.3, 1.4, 1.5, 1.2 cm2/s for human muscle, kidney, spleen, liver, and blood, respectively. We note that models developed in food engineering also regard the cellular tissue as a passive medium (see, e.g., FE ). This assumption is justified by a number of self-consistent results obtained in a variety of experimental works devoted to measuring thermal characteristics of perfused and nonperfused biological tissues (see, e.g., GMLM95  and references therein, and, also, ExC ; V1 ; V2 ; V3 ; GrE ). However, measuring thermal characteristics of living tissue is a difficult task because of its complex structure. So, interpretation of the obtained experimental data and verification of the corresponding mathematical models are related problems which are far from being completed (see, e.g., LX ; Var ), concerning especially strong heating. Nevertheless, up to now all the peculiarities of heat transfer in living tissue are assumed to be due to the convective transport with blood.

During hyperthermia treatment the tissue temperature grows up to the upper boundary of the tissue surviving, i.e. up to $`4446^{}`$C, moreover, during a thermotherapy course based on the thermal tissue coagulation the temperature can attain values about 100 C. So the question of whether the cellular tissue undergoing so strong heating can be described in the framework of the passive medium model is not evident. Therefore the thermal properties of cellular tissue under local strong heating deserve individual investigations to clarify if the physical background of the bioheat transfer theory has to be modified. In particular, the fact that heat propagation in cellular tissue can exhibit anomalous properties was noted in HDE .

In this paper we present preliminary experimental data that have enabled us to doubt whether cellular tissue, at least some of its types, is in fact a passive medium and heat propagation in it does obey Fick’s law.

## II Temperature dynamics in the tissue bulk under strong local heating

First, we studied time variations of the temperature inside a fragment of cow haunch muscle not undergone freeze after death and cut out so that muscle fibers be parallel to its surface (Fig. 1). A small resistance $`r`$ of size $`1\times 2`$ mm was embedded into the tissue at a depth of 3 cm through this free surface oriented horizontally. Passing a current $`I`$ through it we got a local heat source of controlled power. Thermocouples measuring the temperature were located at a distance $`\rho `$ of 1 cm from the resistance along the fibers and in the transverse direction.

We studied the heating dynamics for various values of the initial tissue temperature from 0 C to 30 C and various specimens of muscle, the results were practically identical. Fig. 2 demonstrates typical time variations of the tissue temperature obtained for two values of the heat source power equal to 1.2 W ($`I=300`$ mA, $`r=13`$ O) and 2.7 W ($`I=300`$ mA, $`r=30`$ O). The time origin, $`t=0`$, is the beginning of heating and the vertical dashed line, $`t20`$ min, represents the heating tune-out. All the shown curves were obtained for tissue regions undergone heating for the first time because under strong heating irreversible transformations of the living tissue structure can occur. It should be noted that the curves corresponding to heating power of 1.2 W were obtained for the tissue with the initial temperature $`T_{\mathrm{init}}=20.2^{}`$C, so in Fig. 2 they are shifted down by 1.4 C to make the comparison more easy.

As seen in Fig. 2 the time variations of the tissue temperature along the fibers and in the transverse direction exhibit substantially different features. In fact, for the longitudinal direction the characteristic time $`\tau _l`$ of the temperature variations can be estimated as $`\tau _l68`$ min from the $`e`$-fold decrease after the heating tune-out. The same estimate will be get if we apply to the point where the curve describing the temperature increase under 1.2 W heating transfers from the concave to convex form. For the temperature increase under 2.7 W heating the latter time seems to be even shorter. The analogous curves describing the temperature variations in the transverse direction are characterized by a substantially longer time scale $`\tau _t`$. In particular, the temperature increase curves remain concave up to the heating tune-out. So, if we related the time scales $`\tau _{l,t}`$ to the temperature diffusivities $`D_{l,t}`$ along the fibers and transversely them with the expression $`\tau _{l,t}\rho ^2/D_{l,t}`$ then we would have been to set that $`D_l/D_t3`$. From our point of view the latter estimate contradicts the physical properties of heat transfer in cellular tissue because it does not contain any regions where the heat conductivity takes so different values that could give rise to this estimate.

Under strong heating the cellular tissue coagulates and changes its properties, however, on the time scales under consideration it is essential only for temperatures exceeding $`5060^{}`$C J . Therefore, this effect seems not to be substantial enough to play a key role in the heat propagation through the cellular tissue in the given case. Indeed, if we again assume the tissue temperature to be governed by the diffusion equation with, may be, an anisotropic diffusivity, then the stationary temperature field $`T`$ caused by a $`\delta `$-source will be scaled up as $`(TT_{\mathrm{init}})(TT_{\mathrm{init}})\lambda `$ when all the spatial scales are reduced by the factor $`\lambda `$, i.e. $`\rho \rho /\lambda `$. So, under 2.7 W heating the tissue temperature should attain values about 60 C at points being at a distance of 0.5 cm from the heat source along the fibers and at a substantially shorter distance in the transverse direction. Thus the region where the cellular tissue has to coagulate during the heating is of the relative volume not exceeding 10% with respect to the 1$`\times `$1$`\times 1`$ cm3 region actually tested in the experiment. Under 1.2 W heating the region of the coagulated tissue is still less and even seems not to form at all because of finite dimensions of the resistance.

Fig. 2 also shows that heat propagation in cellular tissue under strong heating must be nonlinear. At least, the change of the heating power from 1.2 W to 2.7 W, i.e. by 2.3 times, corresponds to the temperature increase by 2.7 times at the points distant from the heat source along the fibers and only by 1.4 times at the points in the transverse direction.

Fig. 3 demonstrates our test of the repeat irreversibility of the heating process showing if irreversible transformations of the tissue structure caused by strong heating are essential in the case under consideration. Namely, to answer the question of whether heat transfer in cellular tissue exhibiting such anomalous properties can go in the same way under repeated acts of heating we tested the temperature dynamics for the source power of 1.2 W. The obtained curves present time variations of the tissue temperature at the points shown in Fig. 1 during two successive processes of heating of the same tissue region that are separated in time by 40 min. We see not too significant change in the temperature variations only at the point distant from the heat source along the fibers that occurs at the stage of direct heating. During the repeated heating the temperature increase is characterized solely by a faster growth rate. So the mechanism responsible for these anomalous properties of heat transfer seems to be of reversible nature, at least, at first approximation.

We assume that these heat transfer anomalies are due to motion of the interstitial liquid induced by heating, which is demonstrated in the next section.

## III Temperature dynamics in a tissue layer under strong surface heating

Following actually a standard way of measuring the tissue heat conduction we studied heat propagation through a tissue layer of thickness of 2 cm (Fig. 4). One of its surfaces was kept up at high temperatures about 100 C whereas the temperature growth at the opposite surface was recorded by a thermocouple. The tissue layer was oriented vertically to avoid the gravitation effect. We studied four arrangements distinguished by the fiber orientation with respect to the walls bounding the layer and by the wall permeability to the interstitial liquid. In one case we used the impermeable walls to suppress the induced motion of the interstitial liquid, in the other we applied a grid to make it as free as possible. The distance between the walls was controlled by a certain small fixed force pressing the system so that there were no air-gaps between the tissue and the walls.

It should be noted that in the given situation heat propagation through the cellular tissue differs substantially from the case considered in the previous section. In fact, when the thickness of the tissue layer is much less than the other dimensions and its surface is heated uniformly, the temperature field can practically vary solely along the direction perpendicular to the layer surface. Under such conditions the interstitial liquid flow can appear only in the same direction and inevitably has to be affected by the properties of the opposite layer surface. For local heating of the cellular tissue at internal points distant far from its surface the structure anisotropy should give rise to the interstitial liquid flow of partly swirling form. In this case the effect of the surface properties seems to be weakened.

Fig. 5 exhibits the obtained results. When the tissue layer is bounded by the impermeable walls and the interstitial liquid has no way to move the evolution of tissue temperature is to be governed by the conventional heat diffusion and to depend weakly on the fiber orientation. Exactly such a behavior was observed. Then assuming the tissue temperature to obey the diffusion equation we can directly estimate the corresponding temperature diffusivity $`D`$. Namely, in the given case the tissue temperature at one boundary is fixed and equal approximately to 100 C and at the other we may impose the zero flux condition on the temperature field. The initial temperature is 20 C and according to Fig. 5 the tissue temperature attains values about 50 C at latter boundary in 10 min. Then fitting the solution of the diffusion equation subject to the appropriate conditions to these data we get the estimate $`D2\times 10^3`$ cm2/s. The latter corresponds well to the available data for the tissue temperature diffusivity within the accuracy of our measurements.

In the tissue layer bounded by the grid heat transfer exhibits another behaviour. When the muscle fibers were parallel to the layer boundaries the temperature growth turned out to be of the same form as for the impermeable walls. By contrast, when the tissue fibers were perpendicular to the boundaries the temperature growth deviated substantially from the previous dependence (Fig. 5). We can explain this effect only admitting an essential flux of the interstitial liquid to appear due to the nonuniform heating of the cellular tissue. It is likely that the fiber arrangement gives rise to long channels in the muscle tissue running along the fibers. So, when there are no substantial obstacles to motion of the interstitial liquid nonuniform heating of the tissue is able to induce the interstitial liquid flux affecting remarkably heat transfer. To partly justify the given assumption we note that for the latter fiber orientation the blood efflux through the grid was visually strong in comparison with the other cases when the measured tissue temperature attained the values about 33 C.

In order to make this phenomenon more evident the next section presents its certain visualization.

## IV Mass propagation caused by nonuniform heating

If nonuniform heating of cellular tissue causes the interstitial liquid to flow in it affecting substantially heat propagation itself then the interstitial liquid flux can be directly visualized through the induced propagation of a certain substance dying the tissue. Indeed, diffusion coefficients of simple molecules in cellular tissue (of order 10-5 cm2/s or even less) are much less than the temperature diffusivity, so mass propagation in cellular tissue is to be essentially accelerated by this interstitial liquid flux.

Figure 6 illustrates the experiment showing the effect of nonuniform heating on propagation of aniline dye in the cellular tissue. We put an aniline dye drop in a lunula on the horizontal tissue surface and pressed down the tissue with a glass plate to make the tissue surface plane and, thereby, to exclude the influence of the surface relief. A metallic disc of diameter of 1 cm was embedded into the glass plate placed so that this disc be right over the drop. An ohmic resistance was attached to the external disc surface and passing a current through it we locally heated the tissue. The colour patterns formed by the aniline propagation are represented in Fig. 7. Figure $`(a)`$ exhibits the colour pattern appeared in the tissue without heating after a lapse of 8 hours. As is seen, it is roughly of symmetrical form and the muscle fibers have no substantial effect on the aniline propagation in this case. The opposite situation is shown in figure $`(b)`$ for the aniline propagation under a strong local heating. Already after 15 min the colour pattern attained the previous one in size and got a substantially anisotropic form, the pattern dimension along the fibers is more that twice its mean size in the transverse direction. Moreover, in this figure it is clearly seen as aniline dry propagates inside small grooves on the tissue surface whereas without heating such grooves practically have no effect.

## V Discussion and hypotheses

The presented preliminary experimental data enable us to suggest that heat transfer in cellular tissue under local strong heating is a more complex phenomenon than a simple heat diffusion. On one hand, the observed strong anisotropy in the heat propagation along the muscle fibers and in the transverse direction cannot be explained keeping in mind only the heat conduction through the cellular tissue. Indeed, otherwise, we should assume that the temperature diffusivity along the fibers, $`D_l`$, is fully three times greater than the corresponding value in the transverse direction, $`D_t`$ ($`D_lD_t`$), whereas the cellular tissue contains no elements differing from one another so remarkably in the heat characteristics. On the other hand, the fiber arrangement endows the muscle tissue with a significant anisotropy for motion of the interstitial liquid through the tissue as through a porous medium. Therefore, if the nonuniform temperature field can cause the interstitial liquid to flow then heat propagation in the muscle tissue will be able to exhibit substantial anisotropy in properties. We think that in the presented experiments exactly this phenomenon, i.e., the interstitial liquid flux induced by nonuniform heating and, in turn, its effect on heat transfer was observed. Obviously, it should be nonlinear, which corresponds to the data presented in Fig. 2.

The experimental data shown in Fig. 3 enable us to assume that heat transfer in cellular tissue under such conditions is characterized by the repeat reversibility. The latter means that the time variations in the temperature of a cellular tissue undergoing repeated acts of heating go approximately in the same way. So a possible microscopic mechanism responsible for the induced flux of the interstitial liquid is also of the reversible nature, at least, it is not due to the tissue thermal coagulation or another thermal damage of the tissue. For example, it could be implemented through the liquid redistribution between the tissue cells and the intercellular space due to the change in the osmotic equilibrium induced by heating. In the given case an additional amount of the interstitial liquid appearing in the heated region has to give rise to its flow from this region.

We point out again that the interstitial liquid flux induced by heating seems to be the integral feature of heat transfer in cellular tissue. Therefore, only under special conditions heat propagation in cellular tissue is reduced solely to diffusion. Moreover, mass and heat propagation have to be related with each other through this flux, which is illustrated in Figs. 5 and 7.

Concluding the paper we would like to note that we observed anomalous properties of heat transfer in muscle tissue where the fibers endow the tissue with well developed anisotropy. There are organs, for example, kidney and liver whose cellular tissue does not possess significant anisotropy on the average regarding such an organ as a whole. For rigorously isotropic media the effect of a similar liquid flux induced by heating is likely to be depressed. However, on smaller scales these organs are certainly to exhibit local anisotropy with random orientation in space. Therefore, local strong heating again has to induce the interstitial liquid flux that, may be, changes its direction randomly in different regions of a given organ. In this case heat transfer is also to exhibit anomalous properties, which cannot be described, in general, by the conventional diffusion equation (see, e.g., I92 ).

###### Acknowledgements.

These investigations were support in part by the Russian State Programme “Integration”, Grant 0075.


format: text

text:
# References

Programming and Computer Software. Vol. 26, No. 5, 2000, pp.275–280

Original Russian publication: Programmirovanie. 2000, No. 5, pp.53–62.

Universal Numerical Algorithms and Their Software

Implementation

G.L. Litvinov and E.V. Maslova

International Sophus Lie Center,

Nagornaya ul. 27–4–72, Moscow 113186 Russia

e-mail: litvinov@islc.msk.su, glitvinov@mail.ru

Received April 10, 2000

## Abstract

The concept of a universal algorithm is discussed. Examples of this kind of algorithms are presented. Software implementations of such algorithms in $`C^{++}`$-type languages are discussed together with means that provide for computations with an arbitrary accuracy. Particular emphasis is placed on universal algorithms of linear algebra over semirings.

INTRODUCTION

Modern achievements in software development and mathematics make us consider numerical algorithms and their classification from a new point of view. Conventional numerical algorithms are oriented to software (or hardware) implementation with the use of the floating point arithmetic and fixed accuracy. However, it is often desirable to perform computations with variable (and arbitrary) accuracy. For this purpose, algorithms are required that are independent of the accuracy of computations and of a particular computer representation of numbers. In fact, many algorithms are not only independent of the computer representation of numbers, but also of concrete mathematical (algebraic) operations on data. In this case, operations may be considered as variables. Such algorithms are implemented by generic programs based on the abstract data types technique (abstract data types are defined by the user, in addition to predefined types of the language used). The corresponding program tools appeared as early as in Simula-67, but modern object-oriented languages (like $`C^{++}`$, see, e.g. ) are more convenient for generic programming.

The concept of a generic program was introduced by many authors; for example, in , such programs were called program schemes. In this paper, we discuss universal algorithms implemented as generic programs and their specific features. This paper is closely related to papers , in which the concept of a universal algorithm was defined and software and hardware implementation of such algorithms was discussed in connection with problems of idempotent mathematics . In this paper, the emphasis is placed on software implementation of universal algorithms, computations with arbitrary accuracy, universal algorithms of linear algebra over semirings, and their implementation in $`C^{++}`$.

1. UNIVERSAL ALGORITHMS

Computational algorithms are constructed on the basis of certain basic operations. Basic operations manipulate data that describe “numbers”. These “numbers” are elements of a “numerical domain”, i.e., a mathematical object like the field of real numbers, the ring of integers, or an idempotent semiring of numbers (idempotent semirings and their role in idempotent mathematics are discussed in  and below in this paper). In every particular computation, elements of the numerical domains are replaced by their computer representations, i.e., by elements of certain finite models of these domains. Examples of models that can be conveniently used for computer representation of real numbers are provided by various modifications of floating point arithmetics, approximate arithmetics of rational numbers , and interval arithmetics. The difference between mathematical objects (“ideal” numbers) and their finite models (computer representations) results in computational (e.g., rounding) errors.

An algorithm is called universal if it is independent of a particular numerical domain and (or) of its computer representation. A typical example of a universal algorithm is computation of the scalar product $`(x,y)`$ of two vectors $`x=(x_1,\mathrm{},x_n)`$ and $`y=(y_1,\mathrm{},y_n)`$ by the formula $`(x,y)=x_1y_1+\mathrm{}+x_ny_n`$. This algorithm (formula) is independent of a particular domain and its computer implementation, since the formula is defined for any semiring. It is clear that one algorithm can be more universal than another. For example, the simplest rectangular formula provides the most universal algorithm for numerical integration; indeed, this formula is valid even for idempotent integration (over any idempotent semiring ). Other quadrature formulas (e.g., combined trapezoid or Simpson formulas) are independent of computer arithmetics and can be used (e.g., in the iterative form) for computations with arbitrary accuracy. In contrast, algorithms based on Gauss–Jacobi formulas are designed for fixed accuracy computations: they include constants (coefficients and nodes of these formulas) defined with fixed accuracy. Certainly, algorithms of this type can be made more universal by including procedures for computing the constants; however, this results in an unjustified complication of the algorithms.

Computer algebra algorithms used in such systems as Mathematica, Maple, REDUCE, and others are higly universal. Standard algorithms used in linear algebra can be rewritten in such a way that they will be valid over any field and complete idempotent semiring (including semirings of intervals; see , where an interval version of the idempotent linear algebra and the corresponding universal algorithms are discussed).

As a rule, iterative algorithms (beginning with the successive approximation method) for solving differential equations (e.g., methods of Euler, Euler–Cauchy, Runge–Kutta, Adams, a number of important versions of the difference approximation method, and the like), methods for calculating elementary and some special functions based on the expansion in Taylor’s series and continuous fractions (Padé approximations) and others are independent of the computer representation of numbers.

2. UNIVERSAL ALGORITHMS AND ACCURACY

OF COMPUTATIONS

Calculations on computers usually are based on a floating-point arithmetic with a mantissa of a fixed length; i.e., computations are performed with fixed accuracy. Broadly speaking, with this approach only the relative rounding error is fixed, which can lead to a drastic loss of accuracy and invalid results (e.g., when summing series and subtracting close numbers). On the other hand, this approach provides rather a high speed of computations. Many important numerical algorithms are designed to use a floating-point arithmetic (with fixed accuracy) and ensure the maximum computation speed. However, these algorithms are not universal. The above mentioned Gauss–Jacobi quadrature formulas, computation of elementary and special functions on the basis of the best polynomial or rational approximations or Padé–Chebyshev approximations, and some others belong to this type. Such algorithms use nontrivial constants specified with fixed accuracy.

Recently, problems of accuracy, reliability, and authenticity of computations (including the effect of rounding errors) have come to the fore; in part, this fact is related to the ever-increasing performance of computer hardware. When errors in initial data and rounding errors strongly affect the computation results (ill-posed problems, analysis of stability of solutions, etc.), it is often useful to perform computations with improved and variable accuracy. In particular, the rational arithmetic, in which the rounding error is specified by the user , can be used for this purpose. This arithmetic is a useful complement to the interval analysis . The corresponding computational algorithms must be universal (in the sense that they must be independent of the computer representation of numbers).

4. MATHEMATICS OF SEMIRINGS

A broad class of universal algorithms is related to the concept of a semiring. We reiterate here the definition of a semiring (see, e.g., ). Let $`S`$ be a set on which associative binary operations $``$ and $``$, called addition and multiplication, respectively, are defined. We assume that addition is commutative and that multiplication is distributive over addition; i.e., $`x(yz)=(xy)(xz)`$ and $`(xy)z=(xz)(yz)`$ for all $`x,y,zS`$. In this case, $`S`$ is called a semiring. We assume that the semiring $`S`$ contains identity $`\mathrm{𝟏}`$ and zero $`\mathrm{𝟎}`$; i.e., $`\mathrm{𝟏}x=x\mathrm{𝟏}=x`$ and $`\mathrm{𝟎}x=x`$, $`\mathrm{𝟎}x=x\mathrm{𝟎}=\mathrm{𝟎}`$; in addition, $`\mathrm{𝟎}\mathrm{𝟏}`$. As is customary, we sometimes omit the multiplication symbol.

A semiring $`S`$ is called commutative if multiplication $``$ is commutative. A semiring $`S`$ is called idempotent if $`xx=x`$ for all $`xS`$. If a semiring is a group under addition, it is called a ring (in this case, it cannot be idempotent). If every nonzero element of a commutative ring (semiring) is invertible under multiplication, this ring (semiring) is called a field (semifield).

The best known and most important examples of semirings are “numerical” semirings consisting of real numbers. For example, the set $`𝐑`$ of all real numbers is a field under ordinary arithmetic operations; i.e., $`=+`$, $`=`$, $`\mathrm{𝟎}=0`$, $`\mathrm{𝟏}=1`$. The set $`𝐑_{\mathrm{max}}=𝐑\{\mathrm{}\}`$ equipped with operations $`=\mathrm{max}`$ and $`=+`$ provides an example of an idempotent semiring (and semifield). Here $`\mathrm{𝟎}=\mathrm{}`$ and $`\mathrm{𝟏}=0`$. This semifield is often called the Max-Plus algebra. The semiring $`𝐑_{\mathrm{min}}=𝐑\{+\mathrm{}\}`$ equipped with operations $`=\mathrm{min}`$ and $`=+`$ is isomorphic to the Max-Plus algebra. Here $`\mathrm{𝟎}=+\mathrm{}`$ and $`\mathrm{𝟏}=0`$. Another example is the set $`S_{\mathrm{max},\mathrm{min}}^{[a,b]}`$ consisting of the elements of an interval $`[a,b]`$, where $`\mathrm{}a<b+\mathrm{}`$, equipped with operations $`=\mathrm{max}`$ and $`=\mathrm{min}`$; here $`\mathrm{𝟎}=a`$ and $`\mathrm{𝟏}=b`$. This commutative semiring is not a semifield.

An important example of a noncommutative semiring is the set $`\mathrm{Mat}_n(S)`$ of all matrices of order $`n\times n`$ with elements from a commutative semiring $`S`$ with ordinary standard operations. The sum of matrices $`A=(a_{ij})`$ and $`B=(b_{ij})`$ is the matrix $`AB=(a_{ij}b_{ij})`$, and the product of these matrices is the matrix $`AB=(_{k=1}^na_{ik}b_{kj})`$, where $`i,j=1,\mathrm{},n`$. Operations on rectangular matrices can be defined similarly. Zero $`O`$ and identity $`I`$ in $`\mathrm{Mat}_n(S)`$ are defined in the conventional way. If the semiring $`S`$ is idempotent, then $`\mathrm{Mat}_n(S)`$ is also idempotent. Many other important examples can be found in , , , .

On any idempotent semiring, a canonical partial order $``$ is defined by the following rule: $`xy`$ is equivalent to $`xy=y`$. Moreover, $`xy=sup\{x,y\}`$ with respect to the canonical order. The canonical order is compatible with the semiring addition and multiplication in the common way. For the semirings $`𝐑_{\mathrm{max}}`$ and $`S_{\mathrm{max},\mathrm{min}}^{[a,b]}`$, the canonical order coincides with the standard order $``$ defined on the set of real numbers; for the semiring $`𝐑_{\mathrm{min}}`$, it is inverse to the standard order.

There exists a (heuristic) correspondence between important, useful, and interesting constructs and results of traditional mathematics over fields and similar constructs and results of idempotent mathematics (i.e., mathematics over idempotent semirings). This idempotent correspondence principle is closely related to the Bohr correspondence principle in quantum mechanics. Traditional mathematics can be considered as a “quantum” theory and idempotent mathematics as its “classical” analogue (see ). Consistent application of the idempotent correspondence principle leads to various and surprising results, including a methodology for constructing universal algorithms and patenting computer devices , .

The fundamental equations in quantum theory are linear (superposition principle). There is an idempotent version of the superposition principle ): the Hamilton–Jacobi equation, i.e., the basic (nonlinear) equation of classical mechanics, can be considered as linear over the semiring $`𝐑_{\mathrm{min}}`$; various modifications of the Bellman equation, i.e., the basic equation of optimization theory, are also linear over appropriate idempotent semirings. For example, the finite-dimensional time-independent Bellman equation can be written as

$$X=AXB,$$
(1)
where $`A`$ is a square matrix with elements from an idempotent semiring $`S`$ and $`X`$ and $`B`$ are vectors (or matrices) with elements from $`S`$. The solution $`X`$ is found from (1) when $`A`$ and $`B`$ are given.

In particular, standard problems in dynamic programming correspond to the case $`𝐑_{\mathrm{max}}`$, and the well-known shortest path problem corresponds to $`S=𝐑_{\mathrm{max}}`$. It is shown in  that the principal optimization algorithms for finite graphs correspond to standard methods for solving systems of linear equations of form (1) over semirings. The Bellman algorithm for the shortest path problem corresponds to a semiring version of the Jacobi method; the Ford algorithm corresponds to the Gauss–Seidel iterative method; and so on. These algorithms are universal and may be used for solving linear algebra problems over a broad class of semirings that includes all idempotent semirings and all fields.

Idempotent analogues of standard numerical algorithms are very important and can be used systematically for solving, for example, optimization problems. Linear algebra algorithms are of prime importance, since standard infinite-dimensional linear problems over semirings can be reduced to finite-dimensional (or finite) approximations, and nonlinear algorithms can often be approximated by linear ones.

We note that the available methods used for parallelizing linear algebra algorithms can be applied to their semiring analogues.

5. UNIVERSAL LINEAR ALGEBRA

ALGORITHMS OVER SEMIRINGS

The most important linear algebra problem is solving the system of linear equations

$$AX=B,$$
(2)
where $`A`$ is a matrix with elements from the basic field and $`X`$ and $`B`$ are vectors (or matrices) with elements from the same field. It is required to find $`X`$ if $`A`$ and $`B`$ are given. If $`A`$ in (2) is not the identity matrix $`I`$, then system (2) can be written in form (1), i.e.,

$$X=AX+B.$$
$`(1^{})`$
It is well known that form (1) or ($`1^{}`$) is convenient for using the successive approximation method. Applying this method with the initial approximation $`X_0=0`$, we obtain the solution

$$X=A^{}B,$$
(3)
where

$$A^{}=I+A+A^2+\mathrm{}+A^n+\mathrm{}$$
(4)
On the other hand, it is clear that

$$A^{}=(IA)^1,$$
(5)
if the matrix $`IA`$ is invertible. The inverse matrix $`(IA)^1`$ can be considered as a regularized sum of the formal series (4).

The above considerations can be extended to a broad class of semirings. The unary operation $`AA^{}`$ in $`\mathrm{Mat}_n(S)`$ is defined (partially) if a unary (partial) operation $`xx^{}`$, called closure, is defined on the semiring $`S`$ such that the identity

$$x^{}=\mathrm{𝟏}(x^{}x)=\mathrm{𝟏}(xx^{})$$
(6)
holds true if $`x^{}`$ is defined. It follows from (6) that

$$x^{}=\mathrm{𝟏}xx^2\mathrm{}x^nx^{}x^{n+1}$$

for any positive integer $`n`$; thus, $`x^{}`$ can be considered as a regularized sum of the formal series

$$x^{}=\mathrm{𝟏}xx^2\mathrm{}x^n\mathrm{}$$

If $`S`$ is a field, then, by definition, $`x^{}=(\mathrm{𝟏}x)^1`$ for any $`x\mathrm{𝟏}`$. If $`S`$ is an idempotent semiring, then, by definition

$$x^{}=\mathrm{𝟏}xx^2\mathrm{}=sup\{\mathrm{𝟏},x,x^2,\mathrm{},x^n\mathrm{}\},$$
(7)
if this supremum (with respect to the canonical order $``$) exists. In this case, $`x^{}=\mathrm{𝟏}`$ if $`x\mathrm{𝟏}`$. Therefore, $`x^{}=\mathrm{𝟏}`$ in the semiring $`S^{[a,b]}`$ for all $`x`$. For the semifield $`𝐑_{\mathrm{max}}`$ the closure operator $`xx^{}`$ is not defined for $`\mathrm{𝟏}x`$ (however, $`𝐑_{\mathrm{max}}`$ can be supplemented by $`+\mathrm{}`$, which turns this semifield into a semiring; in this case, $`x^{}=+\mathrm{}`$ for $`\mathrm{𝟏}x`$). It is clear that, for $`x\mathrm{𝟏}`$, $`x^{}=\mathrm{𝟏}`$ in $`𝐑_{\mathrm{max}}`$, as well as in other idempotent semirings. These examples show that the closure $`x^{}`$ of $`x`$ is often calculated very simply for idempotent semirings.

The closure operation for matrix semirings $`\mathrm{Mat}_n(S)`$ can be defined and computed in terms of the closure operation for $`S`$; some methods are described in . One such method is described below ($`LDM`$-factorization). The closure operation $`AA^{}`$ in $`\mathrm{Mat}_n(S)`$ satisfies identity (6), which implies that if $`A^{}`$ is defined, then $`X=A^{}B=A^{}B`$ is the solution to the matrix equation (1).

Consider a nontrivial universal algorithm applicable to matrices over semirings with the closure operation defined.

Example: Semiring $`LDM`$-Factorization

Factorization of a matrix into the product $`A=LDM`$, where $`L`$ and $`M`$ are lower and upper triangular matrices with a unit diagonal, respectively, and $`D`$ is a diagonal matrix, is used for solving matrix equations $`AX=B`$ . We construct a similar decomposition for the Bellman equation $`X=AXB`$.

For the case $`AX=B`$, the decomposition $`A=LDM`$ induces the following decomposition of the initial equation:

$$LZ=B,DY=Z,MX=Y.$$
(8)
Hence, we have

$$A^1=M^1D^1L^1,$$
(9)
if $`A`$ is invertible. In essence, it is sufficient to find the matrices $`L`$, $`D`$ and $`M`$, since the linear system (8) is easily solved by a combination of the forward substitution for $`Z`$, the trivial inversion of a diagonal matrix for $`Y`$, and the back substitution for $`X`$.

Using (8) as a pattern, we can write

$$Z=LZB,Y=DYZ,X=MXY.$$
(10)
Then

$$A^{}=M^{}D^{}L^{}.$$
(11)
A triple $`(L,D,M)`$ consisting of a lower triangular, diagonal, and upper triangular matrices is called an $`LDM`$-factorization of a matrix $`A`$ if relations (10) and (11) are satisfied. We note that in this case, the principal diagonals of $`L`$ and $`M`$ are zero.

The modification of the notion of $`LDM`$-factorization used in matrix analysis for the equation $`AX=B`$ is constructed by analogy with the construct suggested by Carré in  for $`LU`$-factorization.

We stress that the algorithm described below can be applied to matrix computations over any semiring under the condition that the unary operation $`aa^{}`$ is applicable every time it is encountered in the computational process. Indeed, when constructing the algorithm, we use only the basic semiring operations of addition $``$ and multiplication $``$ and the properties of associativity, commutativity of addition, and distributivity of multiplication over addition.

If $`A`$ is a symmetric matrix over a semiring with a commutative multiplication, the amount of computations can be halved, since $`M`$ and $`L`$ go into each other under transposition.

We begin with the case of a triangular matrix $`A=L`$ (or $`A=M`$). Then, finding $`X`$ is reduced to the forward (or back) substitution.

Forward substitution

We are given:

* $`L=l_j^i_{i,j=1}^n`$, where $`l_j^i=\mathrm{𝟎}`$ for $`ij`$ (a lower triangular matrix with a zero diagonal);
* $`B=b^i_{i=1}^n`$.

It is required to find the solution $`X=x^i_{i=1}^n`$ to the equation $`X=LXB`$. The program fragment solving this problem is as follows.

| for $`i=1`$ to $`n`$ do |
| --- |
| { | $`x^i:=b^i`$; |
|  | for $`j=1`$ to $`i1`$ do |
|  |  | $`x^i:=x^i(l_j^ix^j)`$; } |

Back substitution

We are given

* $`M=m_j^i_{i,j=1}^n`$, where $`m_j^i=\mathrm{𝟎}`$ for $`ij`$ (an upper triangular matrix with a zero diagonal);
* $`B=b^i_{i=1}^n`$.

It is required to find the solution $`X=x^i_{i=1}^n`$ to the equation $`X=MXB`$. The program fragment solving this problem is as follows.

| for $`i=n`$ to 1 step $`1`$ do |
| --- |
| { | $`x^i:=b^i`$; |
|  | for $`j=n`$ to $`i+1`$ step $`1`$ do |
|  |  | $`x^i:=x^i(m_j^ix^i)`$; } |

Both algorithms require $`(n^2n)/2`$ operations $``$ and $``$.

Closure of a diagonal matrix

We are given

* $`D=\mathrm{diag}(d_1,\mathrm{},d_n)`$;
* $`B=b^i_{i=1}^n`$.

It is required to find the solution $`X=x^i_{i=1}^n`$ to the equation $`X=DXB`$. The program fragment solving this problem is as follows.

| for $`i=1`$ to $`n`$ do |
| --- |
|  | $`x^i:=(d_i)^{}b^i`$; |

This algorithm requires $`n`$ operations $``$ and $`n`$ multiplications $``$.

General case

We are given

* $`L=l_j^i_{i,j=1}^n`$, where $`l_j^i=\mathrm{𝟎}`$ if $`ij`$;
* $`D=\mathrm{diag}(d_1,\mathrm{},d_n)`$;
* $`M=m_j^i_{i,j=1}^n`$, where $`m_j^i=\mathrm{𝟎}`$ if $`ij`$;
* $`B=b^i_{i=1}^n`$.

It is required to find the solution $`X=x^i_{i=1}^n`$ to the equation $`X=AXB`$, where $`L`$, $`D`$, and $`M`$ form the $`LDM`$-factorization of $`A`$. The program fragment solving this problem is as follows.

| FORWARD SUBSTITUTION |
| --- |
| for $`i=1`$ to $`n`$ do |
| { $`x^i:=b^i`$; |
| for $`j=1`$ to $`i1`$ do |
| $`x^i:=x^i(l_j^ix^j)`$; } |
| CLOSURE OF A DIAGONAL MATRIX |
| for $`i=1`$ to $`n`$ do |
| $`x^i:=(d_i)^{}b^i`$; |
| BACK SUBSTITUTION |
| for $`i=n`$ to 1 step $`1`$ do |
| { for $`j=n`$ to $`i+1`$ step $`1`$ do |
| $`x^i:=x^i(m_j^ix^j)`$; } |

Note that $`x^i`$ is not initialized in the course of the back substitution. The algorithm requires $`n^2n`$ operations $``$, $`n^2`$ operations $``$, and $`n`$ operations $``$.

LDM-factorization

We are given

* $`A=a_j^i_{i,j=1}^n`$.

It is required to find the $`LDM`$-factorization of $`A`$: $`L=l_j^i_{i,j=1}^n`$, $`D=\mathrm{diag}(d_1,\mathrm{},d_n)`$, and $`M=m_j^i_{i,j=1}^n`$, where $`l_j^i=\mathrm{𝟎}`$ if $`ij`$, and $`m_j^i=\mathrm{𝟎}`$ if $`ij`$.

The program uses the following internal variables:

* $`C=c_j^i_{i,j=1}^n`$;
* $`V=v^i_{i=1}^n`$;
* $`d`$.

| INITIALISATION |
| --- |
| for $`i=1`$ to $`n`$ do |
|  | for $`j=1`$ to $`n`$ do |
|  |  | $`c_j^i=a_j^i`$; |
| MAIN LOOP |
| for $`j=1`$ to $`n`$ do |
| { | for $`i=1`$ to $`j`$ do |
|  |  | $`v^i:=a_j^i`$; |
|  | for $`k=1`$ to $`j1`$ do |
|  |  | for $`i=k+1`$ to $`j`$ do |
|  |  |  | $`v^i:=v^i(a_k^iv^k)`$; |
|  | for $`i=1`$ to $`j1`$ do |
|  |  | $`a_j^i:=(a_i^i)^{}v^i`$; |
|  | $`a_j^j:=v^j`$; |
|  | for $`k=1`$ to $`j1`$ do |
|  |  | for $`i=j+1`$ to $`n`$ do |
|  |  |  | $`a_j^i:=a_j^i(a_k^iv^k)`$; |
|  | $`d=(v^j)^{}`$; |
|  | for $`i=j+1`$ to $`n`$ do |
|  |  | $`a_j^i:=a_j^id`$; } |

This algorithm requires $`(2n^33n^2+n)/6`$ operations $``$, $`(2n^3+3n^25n)/6`$ operations $``$, and $`n(n+1)/2`$ operations $``$. After its completion, the matrices $`L`$, $`D`$, and $`M`$ are contained, respectively, in the lower triangle, on the diagonal, and in the upper triangle of the matrix $`C`$. In the case when $`A`$ is symmetric about the principal diagonal and the semiring over which the matrix is defined is commutative, the algorithm can be modified in such a way that the number of operations is reduced approximately by a factor of two. For details see .

Other examples can be found in , .

6. SOFTWARE IMPLEMENTATION

OF UNIVERSAL ALGORITHMS

Object-oriented languages (e.g., $`C^{++}`$ and Java) and programming systems that allow abstract data types to be defined provide convenient means for the software implementation of universal algorithms. In this case, program units can operate with abstract (and variable) operations and data types. Specific values of operations are determined by the input data types, these operations (and data types) are implemented by additional program units. Recently, this type of programming technique has been dubbed generic programming (see, e.g., ). To help automate the generic programming, the so-called Standard Template Library (STL) was developed in the framework of $`C^{++}`$ . However, high-level tools, such as STL, possess both obvious advantages and some disadvantages and must be used with caution.

Using the generic programming technique, a program package was developed in $`C^{++}`$ for solving problems in linear algebra over fields and semirings (for various computer implementation of the corresponding numeric domains) and optimization problems on graphs. A hierarchy of abstract data types for basic numeric fields, rings, semifields, and semirings was developed for various computer representations. In particular, various versions of the rational arithmetic  can be used and computations can be performed with any given accuracy. Solving systems of linear Bellman equations over idempotent semirings (by various methods), standard optimization problems on graphs can be solved (the dynamic programming problem, shortest path problem, widest path problem, etc.), including interval versions of those problems , . The system provides a basis for a more powerful program package based on universal algorithms . This system will be described in detail in subsequent publications.

ACKNOWLEDGMENTS

This work was supported by the Russian Foundation for Basic Research, project no. 99–01–01 198. We are grateful to A. Ya. Rodionov and A. N. Sobolevskii for useful advice and help.


text:
# Searching a Database under Decoherence

## I  Introduction

In many computations, the search algorithm is one of the most time-consuming activity of the computer program . With the rapid proliferation of the internet and the increasing need to retrieve information from this ever expanding internet network, it is crucial to consider alternatives to current search algorithms. Recently, quantum algorithms have been shown to reduce hitherto computationally difficult problems for which there are no known polynomial time algorithm into a tractable problem involving polynomial time. This problem involves factoring a large composite integer into its prime factors. Besides prime number factorization, it has also been shown that quantum algorithms like the quantum search algorithm could substantially improve the searching process for a database despite only a square root speed-up.

Indeed, in recent years, this immense potential in quantum algorithm has spurred a fecundity of ideas on the actual implementation of a quantum computer. However, the prospect of actually building one in the next decade remains taunting and seemingly insurmountable with current technology due primarily several reasons, like decoherence and scalability. Unlike classical bits, quantum bits (qubits) are highly susceptible to collapse due to the difficulty of isolating the quantum mechanical systems from its environment. Such decoherence inevitably leads to a loss of information within the system. Thus it is necessary to consider fault tolerant computation through quantum error correction  or more recently using decoherence-free subspaces.

Recently, some interesting works have been done to consider the robustness of the quantum search algorithm under a noisy environment . In ref. , they modeled the decoherence of the search algorithm using a stochastic white noise. In this paper, by considering a specific decoherence model, we obtain analytic expression for the probability of the search item in a quantum search algorithm after a certain number of iterations. Our analytic expression can certainly facilitate the study of the behavior and efficiency of the search algorithm within the model.

This paper is organized as follows. In section II, we briefly describe the Grover’s search algorithm under a noise-free environment. In section III, we describe our decoherence model and in section IV derive an explicit analytical formula for the probability of the search item after a certain number of iterations within the model. In section V, we consider the robustness of the search algorithm under the decoherence model. Finally in section VI, we summarize our results and discuss some implications of the model.

## II Grover’s Search

In a series of seminal papers, Grover  considered a quantum algorithm that could achieve a speed-up in the computational implementation with only $`O(\sqrt{N})`$ for a large structured database with $`N`$ records. This algorithm compares favorably with the classical result which can only execute a search with $`O(N)`$ efficiency. Moreover, it has been shown that Grover’s algorithm is optimal .

Grover’s search algorithm can be summarized neatly into the following main steps: (i) Initialization of the system into a superposition of states; (ii) Subjecting the system to a hashing function, C(S), represented by a unitary operator, $`U`$, given by

$$U=\left(\begin{array}{ccccc}1& 0& 0& \mathrm{}& 0\\ 0& 1& 0& \mathrm{}& 0\\ \mathrm{}& 0& 1& \mathrm{}& \mathrm{}\\ \mathrm{}& 0& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& \mathrm{}& 1\end{array}\right)$$
(1)
(assuming that the first entry satisfies the search criteria and therefore undergoes a rotation of $`\pi `$ radians) followed by a diffusion matrix, $`D`$, defined by

$$D_{ij}=\delta _{ij}+\frac{2}{N}$$
(2)
for $`O(\sqrt{N})`$ iterations; (iii) Measuring the resulting state. The heart of the process therefore hinges significantly on the on step (ii) and the $`O(\sqrt{N})`$ iterations of the matrix $`DUS`$. In this paper, we look closely into this important step and scrutinized the ideas behind the efficiency.

Before proceeding further, we first consider the eigenvalues of the matrix $`S`$. It is not difficult to show that the eigenvalues of $`S`$ all lie on the locus $`|z|=1`$, unit circle, on the complex Argand plane and are explicitly $`\{\underset{(N2)}{\underset{}{1,\mathrm{},1}},\eta ,\eta ^{}\}`$ where the root $`\eta `$ and $`\eta ^{}`$ satisfy the equation $`z^2{\displaystyle \frac{2(N2)}{N}}z+1=0`$. Indeed, it can be further shown that if the matrix $`S`$ is diagonalized as $`P\mathrm{\Lambda }P^1`$, with $`\mathrm{\Lambda }=\text{Diag}[\underset{(N2)}{\underset{}{1,\mathrm{},1}},\eta ,\eta ^{}]`$ and

$$\eta =\frac{1}{N}\left(N22i\sqrt{N1}\right)$$
(3)
then the matrix $`P`$ assumes the form

$$\left(\begin{array}{ccccccc}0& 0& 0& \mathrm{}& 0& i\sqrt{N1}& i\sqrt{N1}\\ 1& 1& 1& \mathrm{}& 1& 1& 1\\ 0& 0& 0& \mathrm{}& 1& 1& 1\\ 0& 0& 0& \mathrm{}& 0& 1& 1\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 1& \mathrm{}& 0& 1& 1\\ 0& 1& 0& \mathrm{}& 0& 1& 1\\ 1& 0& 0& \mathrm{}& 0& 1& 1\end{array}\right).$$
(4)
Since $`|\eta |=1`$, we can rewrite the complex number $`\eta `$ as $`\mathrm{e}^{\mathrm{i}\theta }`$, where, using eq(3),

$$\mathrm{cos}\theta =\frac{N2}{N}\text{ and }\mathrm{sin}\theta =\frac{2\sqrt{N1}}{N}.$$
(5)
Finally, we note that the probability of finding the search item can be evaluated to be $`{\displaystyle \frac{1}{N}}(\mathrm{cos}(m\theta )+\sqrt{N1}\mathrm{sin}(m\theta ))^2`$ while the probability of getting one of the other entries in the database is $`{\displaystyle \frac{1}{N}}(\mathrm{cos}(m\theta ){\displaystyle \frac{1}{\sqrt{N1}}}\mathrm{sin}(m\theta ))^2`$.

## III Noise induced through decoherence

### A Superoperators

Decoherence can be studied and understood within the context of superoperators through the evolution of a bipartite quantum system. A superoperator describes a linear map, $`\$`$, from input density matrix, $`\rho _{\text{in}}`$, to output density matrix, $`\rho _{\text{out}}`$. Although the superoperator, $`\$`$, is not a unitary operator, it is a linear map that preserves hermiticity as well as the trace.

The standard procedure of understanding the behavior of one part of a bipartite quantum system is to extend the system to a larger one (in which the environment (E) is incorporated) so that the evolution of state becomes unitary under the transformation, $`𝒰`$. By assuming complete positivity of the superoperators, it is possible to study the non-unitarity evolution of the state of a subsystem using an operator sum representation. In terms of the operator sum or Kraus representation, we can express this map, $`\$`$, more succinctly as

$$\rho _{\text{out}}=\$(\rho _{\text{in}})=\underset{\mu }{}M_\mu \rho _{\text{in}}M_\mu ^{}.$$
(6)
Unitarity of the operator, $`𝒰_{CE}`$, also requires that the Kraus operators satisfy the condition

$$\underset{\mu }{}M_\mu M_\mu ^{}=1$$
(7)
### B Decoherence Model

For two-spin-1/2 system, we can consider the following Kraus representation

$`M_0`$ $`=`$ $`\left(\begin{array}{cc}\sqrt{1p}& 0\\ 0& \sqrt{1p}\end{array}\right)`$ (10)
$`M_1`$ $`=`$ $`\left(\begin{array}{cc}0& \sqrt{p/3}\\ \sqrt{p/3}& 0\end{array}\right)=\sqrt{p/3}\sigma ^x`$ (13)
$`M_2`$ $`=`$ $`\left(\begin{array}{cc}0& \sqrt{p/3}\\ \sqrt{p/3}& 0\end{array}\right)=\sqrt{p/3}i\sigma ^y`$ (16)
$`M_3`$ $`=`$ $`\left(\begin{array}{cc}\sqrt{p_3}& 0\\ 0& \sqrt{p_3}\end{array}\right)=\sqrt{p_3}\sigma ^z`$ (19)
where $`\sigma ^i,i=x,y,z`$ are the usual Pauli matrices. We can check easily that the unitarity condition for the larger space is satisfied since $`{\displaystyle \underset{\mu }{}}M_\mu ^{}M^\mu =1`$. This decoherence model is sometimes alluded to as the depolarizing channel.

This decoherence model is an idealization of a transmission and storage process with some elegant mathematical symmetry in which the quantum state in the channel can undergo a bit-flip and phase errors. The construction of a depolarizing channel arises from the interaction of a quantum state, $`|\mathrm{\Psi }`$ with an environment in which there is a probability $`p`$ that the quantum state will survive and a probability $`p/3`$ that it would execute a pure bit-flip, a pure phase error or a combination of both. Thus, schematically, we have

$`|\mathrm{\Psi }`$ $`\stackrel{1p}{}`$ $`|\mathrm{\Psi }`$ (20)
$`|\mathrm{\Psi }`$ $`\stackrel{p/3}{}`$ $`\sigma _1|\mathrm{\Psi }`$ (21)
$`|\mathrm{\Psi }`$ $`\stackrel{p/3}{}`$ $`\sigma _3|\mathrm{\Psi }`$ (22)
$`|\mathrm{\Psi }`$ $`\stackrel{p/3}{}`$ $`\sigma _2|\mathrm{\Psi }`$ (23)
Furthermore, it is possible to show that for general spin-1/2 particle in which the density matrix describing the mixed states is given by $`\rho _0=1/2(1+\stackrel{}{n}\sigma )`$, that the original density matrix $`\rho _0`$ evolves under the depolarizing channel as

$$\rho ^{}=p\frac{1}{2}𝐈_\mathrm{𝟐}+(1p)\rho _0$$
(24)
where $`𝐈_\mathrm{𝟐}`$ is the $`2\times 2`$ unit matrix. In order to apply this model to the quantum search algorithm, it is necessary to consider the extension to $`N`$ qubits. It is not hard to verify that in this case, Eq.(24) becomes

$$\rho ^{}=p\frac{1}{N}𝐈_𝐍+(1p)\rho _0$$
(25)
where $`𝐈_𝐍`$ is the $`N\times N`$ unit matrix.

## IV Search with decoherence

Let us now study the effects of the decoherence (depolarizing) model to the quantum search algorithm. As in any search algorithm, we first initialize the system into a superposition of states with equal probabilities. Thus the initial density matrix of the system is

$$(\rho _0)_{ij}=\frac{1}{N}.$$
(26)
This state is then subject to the usual inversion-diffusion process so that the density matrix at the end of the first transformation is

$$\rho _0^{}=S\rho _0S^{}.$$
(27)
In the noise-free Grover’s search algorithm, this transformation is repeated a certain number of times. In our model, before we perform the second iteration, we allow the system to evolve under a depolarizing channel so that using Eq.(25), we have

$$\rho _{i+1}=p\frac{1}{N}𝐈_𝐍+(1p)\rho _i.$$
(28)
It is easy to show inductively that for $`m`$ iterations,

$`\rho _m`$ $`=`$ $`p{\displaystyle \frac{1}{N}}𝐈_𝐍(1+(1p)+\mathrm{}`$ (30)
$`+(1p)^{m1})+(1p)^mS^m\rho _0S^m\rho _i.`$
$`=`$ $`{\displaystyle \frac{1}{N}}\left[1(1p)^m\right]`$ (32)
$`\text{ }+(1p)^mS^m\rho _0S^m`$
A schematic representation of the model is shown in Fig. 1.

After the $`m`$ iterations, the state of the system is measured using the standard Von-Neumann measurement or positive operator valued measurement through appropriate extension to a higher dimensional space, so that the state after the measurement reads

$$\rho _f=\underset{i}{}E_i\rho _mE_i^{}$$
(33)
where $`E_i`$ are orthogonal projection operators (POVM) $`|ii|`$. However, we note that $`_iE_iS^m\rho _0S^mE_i^{}`$ is just the probability of finding the desired state under the search without decoherence. Moreover, using the results in section II or ref. , we have

$$\underset{i}{}E_iS^m\rho _0S^mE_i^{}=\frac{1}{N}(\mathrm{cos}m\theta +\sqrt{N1}\mathrm{sin}m\theta )^2,$$
(34)
where $`\mathrm{cos}\theta =1{\displaystyle \frac{2}{N}}`$ and $`\mathrm{sin}\theta ={\displaystyle \frac{2\sqrt{N1}}{N}}`$. Finally, using Eq.(32) and manipulating using some simple algebra, we find that the probability of searching successfully for the desired state under decoherence is

$`P(m)`$ $`=`$ $`{\displaystyle \frac{1}{N}}(1(1p)^m[(\mathrm{cos}m\theta +`$ (36)
$`\text{ }\sqrt{N1}\mathrm{sin}m\theta )^2+1])`$
Eq. (36) is the main result in our paper.

For the case of a generalized measurement, we can consider the set of $`r`$ POVMs given by $`F_i={\displaystyle \underset{ij}{}}\lambda _{ij}E_j,i=1\mathrm{}r,j=1\mathrm{}N`$ where the parameters $`\lambda _{ij}`$ must satisfy the unitarity condition $`{\displaystyle \underset{i=1}{\overset{r}{}}}\lambda _{ij}=1`$ for each $`j`$. For simplicity, we can chose $`r=N`$ and $`\lambda _{ii}=(1ϵ),\lambda _{ij}=\frac{ϵ}{r1},ij`$. Indeed, it is not difficult to imagine an experimental setup for these POVMs since they can be used to describe possible errors in the photon detectors. In this case, we find the probability of successfully searching for the desired state is

$`P(m)`$ $`=`$ $`(1ϵ)P(m)_{\text{ortho}}+`$ (39)
$`\text{ }{\displaystyle \frac{ϵ}{N}}(1(1p)^m[(\mathrm{cos}m\theta +`$
$`\text{ }{\displaystyle \frac{1}{\sqrt{N1}}}\mathrm{sin}m\theta )^2+1])`$
where $`P(m)_{\text{ortho}}`$ is the probability of the search item under Von Neumann orthogonal measurements and given in Eq.(36).

## V Analysis

Within the search algorithm, the number of iterations needed for the search is usually fixed by maximizing the probability of getting the desired state. Without any decoherence, the number of iterations $`m_{\text{max}}`$ of the first maximum for large $`N`$ is found to be

$$m_{\text{max}}=\frac{\pi \sqrt{N}}{4}.$$
(40)
With decoherence, from Eq.(36), the condition is imposed by the following transcendental equation,

$$\frac{m}{1p}\left(1\frac{N}{2}(1\mathrm{cos}\zeta )\right)=\frac{N\zeta }{2m+1}\mathrm{sin}\zeta $$
(41)
where $`\zeta =2(2m+1)\theta `$.

It is interesting to compare the probability of the search item as a function of the number of iterations, $`m`$, for the specific case for various values of $`p`$. Such a graph for $`N=128`$ and $`p=0.01,0.04`$ and $`0.083394`$ is shown in Fig. 2. The value of $`p=p_c=0.083394`$ refers to the maximum allowed value of the decoherence parameter $`p`$ subject to $`m_{\text{max}}=\pi \sqrt{N}/4`$ such that the probability of the search item do not fall below 0.5. Fig. 2 clearly shows that the presence of decoherence tends to provide a decaying effect on the periodic probability distribution of the search item. In the noise-free situation, we have iterations in which the probability of search items reaches a periodic maximum without significant decay in the maximum probability. In the case of a noisy search, we see that there is gradual reduction in the maxima so that subsequent maxima are usually rendered useless in the search as the magnitude of decoherence gets larger.

As we have noted earlier, if we maintain the number of iterations to the value required under the noise-free situation, namely $`m_{\text{max}}`$, we can solve Eq.(36) for the value of $`p_c`$ such that the probability of the search item is above chance level (i.e. $`P(m)>0.5`$) for $`p<p_c`$. In Fig 3, we plot the critical decoherence parameter $`p`$ as a function of the database size $`n=log_2(N)`$ up to a database size of about $`2\times 10^6`$. This graph shows that as the database gets larger, it is imperative and necessary to overcome the problem of decoherence. At a database size of $`2^{21}2\times 10^6`$, the amount of decoherence allowed is only about $`0.000609`$, which is a very small value indeed.

In a general search problem, the number of iterations can always be controlled. For a noisy search algorithm, it seems from our analysis that it may be better to reduce the number of iterations so that the probability of the search items becomes a little larger. Using a database of $`N=1024`$, we plot the variation of the probability of the search items as a function of iterations $`m`$ for $`p=0.001,0.014`$ and $`p=p_c=0.0274`$ in Fig. 4. The graph clearly shows a shift in the maxima for the curves. Indeed for the near critical value of $`p=p_c`$, we see that the maxima hovers around the value of 0.5 at $`m=21`$. A three-dimensional plot of the probability as a function of $`m`$ and $`p`$ is shown in Fig 5. Indeed, by varying the decoherence parameter $`p`$, it is possible to solve for the value of $`m_{\text{max}}`$ corresponding to the probability of the search item. Fig. 6 shows the value of $`m_{\text{max}}`$ as a function of $`p`$. The step function refers to the integer part of $`m_{\text{max}}`$.

It is interesting to note that the graph of $`m_{\text{max}}`$ corresponding to the maximum probability of the search item is a linear function of $`p`$. A interpolation fit gives $`m_{\text{max}}=24.6254127.7426p`$. If we then plot the probability of the search item as a function of $`p`$ using the integer part of this linear function, $`i.e.\mathrm{Int}(m_{\text{max}})`$ and compare it with the corresponding probability for $`m_{\text{max}}=\mathrm{Int}(\pi \sqrt{N}/4)`$, we see that the difference is reasonably negligible, especially for small $`p`$. This plot is illustrated in Fig. 7.

In the context of a generalized measurement, we find that probability amplitude of the search item is given by the expression in Eq.(39) and this probability generally reduces in the first maximal iteration. The plots of the probability of the search item as a function of the number of iterations $`m`$ for various values of $`p`$ at $`N=128`$ and $`r=131`$ is plotted in Fig. 8. In the plots, we have chosen a reasonably small value for the parameter, $`ϵ=0.1`$, since under experimental conditions, one do not expect a large value for this parameter. In fact, the results show that for sufficiently small $`ϵ`$, it is possible to confine our study to orthogonal Von Neumann measurements.

## VI Discussion and Concluding Remarks

In this paper, we study the effects of the search algorithm analytically and numerically using a simple decoherence model. Several observations are in order. Firstly, we see that the search algorithm can be rendered useless for small degree of decoherence under the depolarizing channel as the database gets larger. Indeed for a database size of about 2 million entries, the degree of decoherence, $`p`$, is a factor of 1000 smaller than the corresponding value for a database size of 4. It appears that for large $`N`$, the factor varies approximately as $`\sqrt{N}`$.

Secondly, as we would expected, decoherence generally diminishes the probability of the search items. With high decoherence, it may not be useful to iterate beyond the first maxima in the probability distribution of the search item.

Thirdly, the number of iterations needed to maximize a search probability decreases with increasing decoherence. This is a somewhat unexpected result since decoherence ostensibly tends to reduces the probability distribution. Nevertheless, although the associated maximum probability has decreased, it is interesting to note that under a noisy channel, one could in general improve the efficiency by reducing the iterations prior to measurement. It is interesting to note here that the number of iterations (considered as a continuous variable, i.e. prior to taking the integer part) falls off linearly with the degree of decoherence, $`p`$. However, if we continue to use the maximum number of iterations in noise-free situation, the difference in probability appears to be less than 10%.

Finally, we note that although the decoherence model that we have employed here is very simple, it certainly provides a wealth of information regarding the behavior of the search algorithm under a noisy environment.