text:
# Robust generation of entanglement between two cavities mediated by short interactions with an atom

## I Introduction

In recent years entanglement, a fundamental feature of many-body quantum mechanical systems, has come to be seen as a useful resource for many tasks in quantum information processing Nielsen and Chuang (2000); Alber et al. (2001); Plenio and Vedral (1998). In particular, certain states have been identified, which can be directly used for many such tasks. Foremost among these are the so-called Bell states, also known in the literature as EPR states, which exhibit the maximum entanglement which a two qubit system can possess. Such states are the fundamental ingredient of in a variety of quantum information processing protocols (see Nielsen and Chuang (2000); Alber et al. (2001); Plenio and Vedral (1998) for details). Bell states are also of relevance to the study of the foundations of quantum mechanics. It was shown by Bell Bell (1964) that these states possess correlations which are incompatible with a local realistic theory (see Hardy (1998) for a review). Thus, a loophole-free demonstration of such correlations would have profound implications for our understanding of the physical world. For these reasons, then, it is important to be able to generate Bell states in the laboratory, and many experimental proposals have been put forward and realized in both optical Wu and Shaknov (1950); Kwiat et al. (1995); Lamas-Linares et al. (2001) and atomic systems Rowe et al. (2001); Fry et al. (1995).

This paper contains a proposal for a scheme utilizing the techniques of cavity QED to generate Bell states between the modes of two spatially separated cavities. Rauschenbeutel and coworkers have generated Bell states between two modes in a single cavity using a Rydberg atom Rauschenbeutel et al. (2001) coherently interacting with each mode in turn. This scheme could be adapted in a straightforward way to generate such states in spatially separated cavity modes. Cabrillo et al have proposed a non-deterministic scheme to generate Bell entangled states between atoms conditionally, by driving them with a weak laser pulse, and subsequently detecting a photon spontaneously emitted by one of the atoms Cabrillo et al. (1999). Protsenko and coworkers have then shown that this scheme can be adapted to implement quantum logic gates Protsenko et al. (2002). At the same time Plenio et al have developed schemes to entangle atoms inside optical cavities and between atoms in different cavities via spontaneous decay Plenio et al. (1999); Bose et al. (1999); Jané et al. (2002). Various other conditional schemes and deterministic schemes have been proposed which employ the detection of photons to generate entanglement and perform quantum computation Knill et al. (2001); Kok et al. (2002); Fiurasek (2002). In schemes such as these, however, a very precise control of the experimental parameters is required. Such a precise control of experimental parameters is often difficult to achieve. Generally it is desirable to devise schemes in which the requirements on the experimental control are as small as possible. Such schemes would then be far more robust to errors and would generally lead to much improved fidelities of the generated states.

In this paper, we present a novel proposal for the robust experimental generation of a single photon entangled Bell state $`|\mathrm{\Psi }^+=\sqrt{\frac{1}{2}}\left(|0_A|1_B+|1_A|0_B\right)`$ in the modes of two spatially separated cavities, labeled $`A`$ and $`B`$. The cavity modes do not interact directly, but the generation of the entangled state is mediated by a resonant atom. This atom, prepared in its excited state, passes through the two cavities, prepared in their vacuum state, and interacts with each for the same short effective interaction time. The internal state of the atom is then measured. Depending on the outcome of this measurement, either a high fidelity Bell state has been generated, or the cavities have returned to the vacuum state, and are ready for the process to be repeated with a fresh atom. The scheme has the advantage that provided the interaction times between the atom and both cavities are the same, which can be attained if a well-collimated atomic beam is utilized, the actual value of this interaction time can vary within a large range of values, without greatly affecting the fidelity of the Bell state produced. This would reduce the need for velocity selection of the atoms and for the sophisticated timing and feedback required to tune the interaction times in other experiments. Additionally, the interaction times, which this scheme would utilize, can be much shorter than in non-conditional schemes, which would reduce the effects of decoherence during the entanglement generation process.

In section II, we introduce the proposed scheme in detail. We begin by discussing the scheme in the limit where the effective duration of the atom-cavity interactions is very small, and the same in both cavities. We then discuss the fidelity of the state generated by a successful run of the scheme with finite interaction times and calculate the probability that a single run of the scheme is successful. We find that, in a wide range of interaction times, the probability of success is suitably high, and the fidelity of the generated state is close to unity. In section III we discuss practical aspects of implementing our scheme in three subsections. In the subsection III.1, we discuss specific physical systems with which the scheme could be implemented. In the two subsections which follow this we consider two practical aspects of the scheme which would pose particular requirements to any implementation of the scheme. In subsection III.2 we look at the effect on the scheme of inaccurate control of the path of the atom through the cavities due to poor collimation of the atomic beam. This would lead to non-equal interaction times in the two cavities, potentially reducing the fidelity of the state generated. In subsection III.3, we discuss how employing detectors with less then perfect efficiency would affect the scheme.

## II The proposed scheme

In this section, we will describe the proposal which is the main focus of this paper, a scheme to generate the entangled Bell state $`|\mathrm{\Psi }^+`$ in two spatially separated microwave cavities. Two identical cavities, labeled A and B, are positioned as in figure 1, such that their axes are parallel, and their centers in alignment. The cavities must be prepared in their vacuum state. This can be achieved, for example, in microwave cavities by cooling to low temperatures. Two-level atoms, whose transition between the ground state $`|g`$ and excited state $`|e`$ is resonant with modes of the cavities, are used to mediate the generation of the entangled state in the cavities. The atoms are prepared in their excited state, $`|e`$. They are then passed, one at a time through the center of both cavities, such that the effective interaction time with each cavity is much less than the period of a one-photon Rabi oscillation. Upon leaving the second cavity and before the next atom enters the cavities, the internal state of the atom leaving the cavities is measured. For brevity in this paper, let us refer to the process of a single atom passing through both cavities and being measured as a “run” of the scheme. The outcome of the measurement is non-deterministic. If the excited state $`|e`$ is detected, then no entangled state has been generated in the cavities, we thus refer to a run which ends in such a measurement as a “null run” After a null run, the cavity is reset to its initial state, so a further run can be implemented immediately. If the ground state $`|g`$ of the atom is detected at the end of a run, a state is generated in the cavities, which, so long as certain conditions are fulfilled, as described below, is very close to the Bell state $`|\mathrm{\Psi }^+`$. The process is complete at this point, and we shall thus label this a “successful run”.

Before we describe the scheme in detail, let us summarize how the resonant interaction between atom and cavity is described. Consider a two-level atom with energy levels $`|g`$ and $`|e`$, separated by a transition energy resonant with a mode of the cavity. The atom is situated at position $`𝒓`$ within the cavity, whose modes are described by field operators $`a`$ and $`a^{}`$. Since the atomic transition is resonant with the cavity mode, the internal atomic state interacts with this mode via the following Hamiltonian, in the interaction picture, where the rotating wave approximation has been employed,

$$H=\mathrm{}g_0u(𝒓)\left(a^{}|ge|+a|eg|\right).$$
(1)
Here $`g_0`$ describes the strength of the coupling. When the atom is situated in the center of the cavity, the one-photon Rabi frequency is $`g_0/\pi `$. The spatial variation of the cavity field is contained in the mode function $`u(𝒓)`$. If we consider an atom moving through the cavity on a classical path $`𝒓(t)`$, the mode function can then be written as a function of time $`u(t)`$. Since all the time-dependence in $`H`$ is contained in $`u(t)`$, the unitary evolution of the system can be written as follows

$$U=\mathrm{exp}\left[\frac{\mathrm{i}\left(_{t_0}^tu(t^{})𝑑t^{}\right)\widehat{H}}{\mathrm{}}\right],$$
(2)
where $`\widehat{H}=H/u(𝒓)`$ is the interaction Hamiltonian for an atom situated in the center of the cavity. In this article, we will consider cases where the atom travels on a straight path through the cavity, and can thus introduce the total effective interaction time $`\tau =_{\mathrm{}}^{\mathrm{}}u(t^{})𝑑t^{}`$ to characterize the interaction. The complete unitary evolution of of the system due to the interaction $`U_{\text{total}}`$ is then

$$U_{\text{total}}=\mathrm{exp}\left[\frac{\mathrm{i}\tau \widehat{H}}{\mathrm{}}\right].$$
(3)
In our scheme, the atom passes through and interacts with two cavities in turn, labeled $`A`$ and $`B`$. Let us label the interaction Hamiltonians for the interaction between the atom and each cavity $`H_A`$ and $`H_B`$ respectively. These take the same form as $`H`$ in equation (1) for each respective cavity. Analogously let us label the interaction Hamiltonians for the atom in the center of each cavity $`\widehat{H}_A`$ and $`\widehat{H}_B`$. The cavities are assumed to be identical and thus the coupling constant $`g_0`$ is the same for both. Let us consider how the state of the system will evolve if a resonant atom, initially in its excited state, passes through and interacts with the cavities, one after the other. If the cavities are initially in some pure state $`|\psi _{\text{cav}}`$. The initial state of the atom-cavities system is thus $`|e|\psi _{\text{cav}}`$. The atom is passed through both cavities such that the effective interaction times are $`\tau _A`$ and $`\tau _B`$. After these interactions, the quantum state of the system has undergone a unitary evolution described by the operator $`U_{AB}`$

$$U_{AB}=e^{\mathrm{i}g_0\tau _B\widehat{H}_B}e^{\mathrm{i}g_0\tau _A\widehat{H}_A}.$$
(4)
In the scheme we propose, the cavities are aligned as shown in figure 1 and the atom passes along the straight line through the center of both cavities at constant velocity. This means that the effective interaction times between the atom and each cavity will be equal. It may, of course, be difficult to control the path of the atom with sufficient accuracy that the interaction times are exactly equal, and the effect of this is discussed in section III.2. For now, however, let us assume that $`\tau _A`$ and $`\tau _B`$ are equal, and write them both $`\tau `$. In the limit when the effective interaction times are very small, i.e. when $`g_0\tau 1`$, $`U_{AB}`$ can be expanded to the first order in $`g_0\tau `$, and takes the following form

$$\begin{array}{cc}\hfill U_{AB}& \mathbf{1}\mathrm{i}\frac{H_A\tau }{\mathrm{}}\mathrm{i}\frac{H_B\tau }{\mathrm{}}\hfill \\ & =\mathrm{𝟏}\mathrm{i}g_0\tau \left[(a_A+a_B)|eg|+(a_A^{}+a_B^{})|eg|\right].\hfill \end{array}$$
(5)
The state of the system after the atom has left the second cavity, is $`|\psi =U_{AB}|\psi _{\text{init}}`$, in this limit

$$|\psi |e|\psi _{\text{cav}}\mathrm{i}g_0\tau |g(a_A^{}+a_B^{})|\psi _{\text{cav}}.$$
(6)
When the state of the atom is now measured, the cavity modes are projected into one of two states, depending on the measurement outcome. If $`|e`$ is detected, the cavity returns to its initial state. This is important for a non-deterministic process, because it means that it can be repeated immediately from the same starting conditions. If the ground state $`|g`$ is detected, the cavity modes are now in the state, neglecting normalization, $`(a_A^{}+a_B^{})|\psi _{\text{cav}}`$. Thus, if the cavities are initially in the vacuum state, the state generated in the cavities would be $`(a_A^{}+a_B^{})|0_A|0_B=|1_A|0_B+|0_A|1_B`$, which, when normalized is the Bell state $`|\mathrm{\Psi }^+`$ introduced above.

If, following a successful run, one were to repeat the scheme immediately and carry on until $`n`$ atoms had been detected in the ground state, in the limit that $`g_0\tau `$ is small, the state generated would have the form $`(a_A^{}+a_B^{})^n|\psi _{\text{cav}}`$, or, if the cavities are initially in their vacuum state, $`(a_A^{}+a_B^{})^n|0_A|0_B`$. This is equivalent to the state produced when a $`n`$-photon Fock state and a vacuum state are incident together on a 50:50 beam splitter. However, numerical results have suggested that, the fidelity of states produced via this method would decrease swiftly with increasing $`n`$. This is due to two reasons, firstly, the short interaction time approximation becomes worse when higher photon numbers are present in the cavities since the timescale of the interactions is faster (the Rabi frequency scales with $`\sqrt{n+1}`$), secondly, when more than one photon is in the cavities, a null run, the measurement of the atom leaving the cavities to be in its excited state, does not reset the state of the cavities to the state before the run, so, as more repetitions are made, the fidelity of the final state gets worse and worse. For these reasons, this does not appear to be a good scheme for the generation of such states.

However, as we will show below, the generation of single-photon Bell states is not affected by these problems. Firstly, we find that the fidelity of the generated state remains close to unity for values of $`g_0\tau `$ much greater than the above approximation would be valid. Secondly, in this case, null runs reset the state of the cavities to the vacuum state exactly, so the fidelity of the state generated is unaffected by the number of runs which were required. In the limit that $`g_0\tau `$ tends to zero, the probability of the detector measuring a ground state is approximately $`2(g_0\tau )^2`$, which, in this limit, is vanishingly small. In a practical scheme, one would need to work in a parameter range, where the probability of success was high enough that few repetitions would be required to achieve a successful run. This is the case for higher values of $`g_0\tau `$, where the above approximation would no longer be valid. Fortunately, starting with the simple pure state $`|e|0_A|0_B`$, it is straightforward to solve the Schrödinger equation and calculate exactly the state of the system after the interactions have taken place. This state $`|\psi `$, under the assumption that the interaction times $`\tau `$ are exactly equal is,

$$\begin{array}{cc}\hfill |\psi =& \mathrm{cos}^2(g_0\tau )|e|0_A|0_B\mathrm{i}\mathrm{cos}(g_0\tau )\mathrm{sin}(g_0\tau )|g|0_A|1_B\hfill \\ & \mathrm{i}\mathrm{sin}(g_0\tau )|g|1_A|0_B.\hfill \end{array}$$
(7)
Let us consider a measurement of the atom’s state. If the excited state is detected the state of the cavities is projected back to the vacuum state, independent of $`\tau `$, as in the approximate case. This resets the cavities to their initial state, so the process can be immediately repeated with a fresh atom.

If the ground state is detected, the following entangled state is generated in the cavity,

$$|\psi _{\text{cav}}=\mathrm{cos}(g_0\tau )|0_A|1_B+|1_A|0_B,$$
(8)
where normalization has been omitted. As $`g_0\tau `$ approaches zero, this tends to the desired state $`|\mathrm{\Psi }^+`$. One can quantify how close this state generated is to $`|\mathrm{\Psi }^+`$ in terms of the fidelity $`F=|\mathrm{\Psi }^+|\psi _{\text{cav}}|^2`$.

$$F=\frac{1}{2}+\frac{\mathrm{cos}(g_0\tau )}{\mathrm{cos}^2(g_0\tau )+1}=1\frac{(g_0\tau )^4}{16}+O\left((g_0\tau )^6\right)$$
(9)
This is plotted in figure 2. We see that the fidelity remains very close to unity for a surprisingly large range of $`g_0\tau `$. For example, the fidelity remains above $`14\times 10^3`$ for $`g_0\tau =0.5`$. Additionally, the function is extremely flat in the wide range of values between $`g_0\tau =0`$ and $`g_0\tau =0.5`$. Therefore the scheme would be insensitive to variations in the interaction times within this range.

For the scheme to be useful, the probability of a single run leading to the successful generation of an entangled state, $`P_{\text{success}}`$, needs to be high enough that prohibitively many repetitions are not required.

$$P_{\text{success}}=1\mathrm{cos}^4(g_0\tau )=2(g_0\tau )^2+O\left((g_0\tau )^4\right)$$
(10)
The success probability is plotted in figure 3. For $`g_0\tau =0.5`$, this probability is approximately 0.4, so a successful run would probably be achieved in 2 or 3 repetitions. The optimal parameter range for the scheme depends upon the fidelity of the state required. The higher the value of $`g_0\tau `$ chosen the higher the success probability will be but the lower the fidelity of the state generated. If, for example, fidelities of $`0.95`$ were acceptable, and a minimum success probability of 0.5 were desired, the scheme could operate between approximately $`g_0\tau =0.6`$ and $`g_0\tau =0.9`$. The exact value of $`g_0\tau `$, however, can lie anywhere in this range, so no fine tuning of the interaction time is required. Note that even if extremely high fidelities such as 0.999 are required, the necessary parameters (up to $`g_0\tau =0.35`$), still allow a success probability of up to 0.22, meaning that a successful run would probably be reached after 4 or 5 repetitions.

## III Practical Considerations

### III.1 Implementation

The above description is quite general and no particular type of atom or cavity has been specified . A number of aspects must be considered in the choice of a physical system to implement the scheme. Firstly, it would be desirable that the entangled states, once generated, would be as long-lived as possible. This favors microwave cavities over optical cavities, since the lifetime of a photon in an optical cavity is currently at the very most one microsecond Lange (2002), whereas microwave cavities with a photon lifetime of a millisecond have been made Raimond et al. (2001). It takes an atom traveling at 500ms<sup>-1</sup> around 20$`\mu `$s to traverse a typical microwave cavity, which is much smaller than this lifetime, so there would be sufficient time for further atoms to probe and interact with the cavity mode before the entangled state has dissipated. If a microwave cavity is used, an atom with a microwave transition is then needed. Microwave transitions occur in atomic fine and hyperfine structure, however dipole transitions between these states are forbidden, and their interaction with the cavity mode would be much too weak to implement this scheme – the interaction times that would be required would be much greater than the photon lifetime of the cavity. Rydberg atoms, on the other hand, although more difficult to prepare, have large dipole moments and thus would interact strongly with cavity modes.

Experiments have been carried out, for example by the Paris group of Haroche and coworkers Raimond et al. (2001), which have parameters close to that required in this scheme. Recall that our scheme requires that the product of parameters $`g_0`$ and $`\tau `$ is at the minimum 0.2 and maximally 0.5 to 0.8, depending on the fidelity of Bell state one wants to generate. In the Paris experiments, a Rydberg atom interacts with a microwave cavity Rydberg atoms interact resonantly with microwave cavity modes with a Rabi frequency of approximately 47kHz. This means that $`g_0=47000\pi s^1=1.48\times 10^5s^1`$. Atoms in an atomic beam from an oven source travel at speeds of the order of hundreds of meters per second. In the Paris experiment atoms with a speed of 500ms<sup>-1</sup> are selected. This means that the effective interaction time is such that a single Rabi oscillation is performed, i.e. $`g_0\tau =\pi `$. This is a factor of 4 to 8 lower than the parameter range for our scheme. It would be difficult to lower $`\tau `$ by using faster atoms, since the velocity of the atoms scales with the square root of the atom oven temperature, so a lower $`g_0`$ would be required. This could be obtained be using a larger cavity, and since $`g_0`$ scales with $`1/\sqrt{V}=1/L^{\frac{3}{2}}`$, so a cavity of mirror separation 3 or 4 times as great as in the Paris experiments would lead to a effective interaction times in the required range. Thus Rydberg atoms and microwave cavities could be employed to implement the scheme.

The disadvantage in using Rydberg atoms is that, at present, the efficiency of state detection schemes is low. In Raimond et al. (2001), for example, they report a detection efficiency of 40%. We will discuss the implications this has for the scheme in subsection III.3. First, however, we consider the effect on the fidelity of the entangled states produced in the scheme if the path of the atom through the cavities is not well controlled and deviates from the line through the center of the cavities.

### III.2 The Atomic Path

If the effective interaction times with both cavities are not exactly the same, this can reduce the fidelity of the entangled state produced by a successful run of the scheme. In our discussion above, we assumed that the two interaction times were exactly equal. In practice, however, it could be difficult to control the path of the atom so precisely. Let us consider first the effect that differing interaction times would have in the fidelity of the Bell state generated by a successful run of the scheme. Let the effective interaction between the atom and cavity A and the atom and cavity B be $`\tau `$ and $`\tau (1ϵ)`$. We can rewrite equation (9), to take the differing interaction times into account, and find the following expression for the fidelity of the entangled state generated by a successful run of the scheme.

$$F=\frac{1}{2}+\frac{\mathrm{cos}(g_0\tau )\mathrm{sin}(g_0\tau )\mathrm{sin}\left(g_0\tau (1ϵ)\right)}{\mathrm{cos}^2(g_0\tau )\mathrm{sin}^2\left(g_0\tau (1ϵ)\right)+\mathrm{sin}^2(g_0\tau )}$$
(11)
The fidelity is plotted for $`g_0\tau =0.5`$ as a function of $`ϵ`$ in figure 4. The asymmetry of the plot is partly an artifact of the choice of parameterization, but, if we take this into account, by plotting $`F`$ again $`\mathrm{ln}(1ϵ)`$, as in figure 5, we see that the asymmetry remains. This is due to the asymmetry in the scheme itself regarding the interactions with the two cavities. When the atom enters the first cavity it is always in the product state $`|e`$, whereas, when entering the second it is always entangled with the first cavity. This makes the interaction with cavity $`B`$ in some sense, slightly weaker, and is the reason that the maximal value of $`F`$ occurs when $`ϵ`$ has a small negative value. The slightly longer interaction time compensates for the interaction being slightly weaker.

Since it is reasonably easy to position the cavities in the desired place to very high precision, let us assume that they are perfectly aligned in the layout illustrated in figure 1.

The factor which would be harder to control would be the atom’s path through the apparatus, since it would be traveling on a ballistic path after being ejected from a heat oven, and although such atomic beams can be highly collimated, there will generally be a small amount of spread in the transverse direction, meaning that the atom’s path will diverge slightly from the central path.

To calculate how much this would affect the interaction times we consider the cavity geometry. The spherical mirrors commonly used in cavity QED experiments support Gaussian cavity modes, which have the following mode function

$$u(x,y,z)=e^{\frac{x^2+y^2}{w_0^2}}\mathrm{cos}\frac{2\pi z}{\lambda },$$
(12)
where $`w_0`$ is the mode waist. The z-axis lies along the line connecting the centers of the two mirrors and the origin is in the center of the cavity. The waist of a Gaussian mode, $`w_0`$, is a function of the cavity geometry and the field wavelength,

$$w_0=\left[\frac{\lambda \sqrt{L}}{2\pi }\sqrt{RL}\right]^{\frac{1}{2}},$$
(13)
where $`L`$ is the separation between and $`R`$ the radius of curvature of the mirrors. If the atomic travels along the central axis of the cavity with constant speed $`v`$, the effective interaction time is $`\sqrt{\pi }w_0/v`$. Let us consider a general straight path through the system, which can be defined in terms of four parameters, $`y_0`$ and $`z_0`$, the initial displacements from the central line in the $`y`$ and $`z`$ directions, and $`\varphi `$ and $`\theta `$ the angles between the atomic path and the central line in the $`y`$ and $`z`$ directions, illustrated for $`z_0`$ and $`\theta `$ in figure 6. We can calculate the effective interaction times between the atom and cavity as it travels along this path with speed $`v`$ and find, that the interaction time with each cavity, in terms of displacements $`\delta _y`$ and $`\delta _z`$ which differ for each cavity, is

$$\begin{array}{cc}\hfill \tau _{\mathrm{eff}}=& \frac{\sqrt{\pi }\omega _0}{v\mathrm{cos}\theta }\mathrm{exp}\left[\frac{\delta _y^2}{\omega _0^2}\left(1\frac{\mathrm{sin}^2\theta }{\mathrm{cos}^2\varphi }\right)\right]\mathrm{exp}\left[\frac{k^2\omega _0^2\mathrm{tan}^2\theta }{4}\right]\hfill \\ & \mathrm{cos}\left[k\delta _zk\delta _y\left(\frac{\mathrm{sin}\theta \mathrm{sin}\varphi }{\mathrm{cos}^2\theta }\right)\right].\hfill \end{array}$$
(14)
For the interaction with cavity $`A`$, $`\delta _y=y_0+\mathrm{cos}(\varphi )D_0`$ and $`\delta _z=z_0+\mathrm{cos}(\theta )D_0`$, and for the interaction with cavity $`B`$ $`\delta _y=y_0+\mathrm{cos}(\varphi )(D_0+D_1)`$ and $`\delta _z=z_0+\mathrm{cos}(\theta )(D_0+D_1)`$. We can use these expressions to calculate $`ϵ`$ in terms of these parameters, written here to the second order in $`y_0`$, $`z_0`$, $`\varphi `$ and $`\theta `$ and their products,

$$\begin{array}{cc}\hfill ϵ& \frac{1}{w_0^2}\left[(D_1\varphi )^2+(D_1\varphi )(2D_0\varphi )+(2y_0)(D_1\varphi )\right]\hfill \\ & +\frac{2\pi ^2}{\lambda ^2}\left[(D_1\theta )^2+(D_1\theta )(2D_0\theta )+(2z_0)(D_1\theta )\right].\hfill \end{array}$$
(15)
Let us discuss the constraints this would have on the collimation of a typical experiment. In the cavity QED experiments of Haroche and coworkers Brune et al. (1996); Raimond et al. (2001), Rydberg atoms interact resonantly with microwave cavities. In a typical experiment, a cavity mode with waist $`w_0=5.97`$mm and wavelength $`\lambda =5.87`$mm is employed. For $`ϵ`$ to be small, the quantities in parentheses in equation (15) must be much smaller than $`w_0`$ and $`\lambda /\sqrt{2}\pi =1.32`$mm. This means that the atomic beam must be collimated so that the effective beam radius is much smaller than this distance. In Brune et al. (1996) an effective beam radius of $`0.25`$mm is reported. If we assume from this that, in the worst case, this would mean that $`y_0,z_0`$0.25mm and $`D_1\varphi ,D_1\theta `$0.25mm, we can estimate that $`ϵ`$ would be less than $`0.2`$. If $`g_0\tau =0.8`$ and $`ϵ=0.2`$, this would corresponds to a reduction of the fidelity of the entangled state produced by a successful run from 0.96 if both interaction times are exactly equal to 0.93. Thus, with the beam collimation currently available in the laboratory, Bell states with a fidelity high enough, for example, to exhibit significant violations of the Bell inequalities Hardy (1998), could be generated.

Equation (15) also provides another reason why optical cavities would be unsuitable for the scheme. The typical waist of an optical cavity tends to be much smaller than that of a microwave cavity , so the demands on the atomic beam collimation required if optical cavities were used would be extremely high.

### III.3 Detector Efficiency

Our analysis in the previous section assumes that the atomic state detector has perfect efficiency. In practice, the detection efficiency will be less than unity. Indeed, as mentioned above, current state detection methods for Rydberg atoms have an efficiency of just 40% Raimond et al. (2001). The detection process ionizes the atom, destroying the state, so increased efficiency could not be obtained by placing detectors in series. A single detection failure will disrupt the scheme, since it will cause a mixed state to be created in the cavities. The scheme must then be halted, and one would then have to wait until this state would have dissipated from the cavities, and the cavities have returned to the vacuum state. Otherwise, if further atoms are sent through the cavities immediately, they will interact with the mixed state, and any “successful” run, will have generated a mixed state with much reduced fidelity. Rather than halting the flow of atoms through the cavities, their interaction with the cavities could be prevented for the cavity dissipation time by the application of an electric field to the system, to create a Stark shift in the atoms such that they are no longer resonant with the cavities. Therefore, one would like a detection efficiency high enough that the probability of a detection failure, during the typical number of runs needed before a successful ground state measurement is made, is low. The mean number of runs to generate the entangled state in the cavity, if the detector were ideal, is simply the inverse of the success probability. The probability that the detector works every time during the process is therefore simply $`D^{1/P_{\text{success}}}`$, which for all values of $`P_{\text{success}}`$ is less than or equal to $`D`$. Therefore, for a reliable scheme, a high detection efficiency would be desirable.

For example, let us consider an implementation of this scheme with Rydberg atoms traveling through the cavities such that $`g_0\tau =0.5`$. The success probability $`P_{\text{success}}`$ is 40.7% and with current detectors with efficiency 40%, $`P_{\text{det}}`$ would be around 10%. This would mean, typically one would have to repeat the whole process, including preparation of the cavities, 10 times before it could reach its successful conclusion. However, in light of the fact that the cavity dissipation time is of the order of milliseconds, and the time taken for each run much less than this, even in this case, the time needed to repeat the scheme enough times to generate the Bell state would be a fraction of a second.

## IV Conclusion

We have proposed a scheme for the generation of high-fidelity Bell states between two spatially separated cavity modes. The scheme is non-deterministic, but we have shown that within the range of parameters $`g_0\tau `$ between approximately 0.3 and 0.9, fidelities higher than 0.95 are obtained, with success probabilities for a single run greater than 1/5 for this entire range.

The most appropriate physical system to implement this scheme would be a combination of microwave cavities and Rydberg atoms. The low detection efficiency for Rydberg states would increase the number of times which the scheme would need to be repeated before a successful run, and require extra time after each detection failure to allow the mixed state produced in the cavities to dissipate. Nevertheless, the scheme would still be successful within a reasonable number of repetitions. The scheme requires that the atomic beam used is highly collimated otherwise the fidelity of the states produced may b degraded, but the collimation which has already been achieved in current experiments is high enough that this effect would be small.

The principle of using brief interactions and measurement to generate entanglement, could be adapted to many other physical systems, and may be especially useful in systems where interaction times are hard to control.

###### Acknowledgements.

This work was supported by EPSRC, the EQUIP project of the European Union, European Science Foundation program on “Quantum Information theory and Quantum Computing”, Hewlett-Packard Ltd. and US-Army grant DAAD19-02-1-0161.

Note added in proof: After this work was completed, we became aware of Messina (2002), in which the author presents a different scheme for the generation of Bell states in two cavities via the passage of a single atom through the cavities.


format: text

text:
# 1 Abstract

## 1 Abstract

It is shown how various ideas that are well established for the solution of Poisson’s equation using plane wave and multigrid methods can be combined with wavelet concepts. The combination of wavelet concepts and multigrid techniques turns out to be particularly fruitful. We propose a modified multigrid V cycle scheme that is not only much simpler, but also more efficient than the standard V cycle. Whereas in the traditional V cycle the residue is passed to the coarser grid levels, this new scheme does not require the calculation of a residue. Instead it works with copies of the charge density on the different grid levels that were obtained from the underlying charge density on the finest grid by wavelet transformations. This scheme is not limited to the pure wavelet setting, where it is faster than the preconditioned conjugate gradient method, but equally well applicable for finite difference discretizations.

## 2 Introduction

Poisson’s equation and Schrödinger’s equation are the central equations for atomistic simulations. In case force fields  are used, Poisson’s equation handles the long range electrostatic interactions. If the forces are calculated quantum mechanically, electronic structure calculations have to be performed. Selfconsistent electronic structure calculations require the solution of a system where Schrödinger’s equation is coupled with Poisson’s equation. Multiscale approaches for the solution of these two central equations are widely used, since they are much more efficient for big system sizes than traditional approaches.

The multigrid method has been used for the solution of Poisson’s equation in the context of classical molecular dynamics simulations  and it has been used by various groups for electronic structure calculations. Several flavors of multigrid for electronic structure calculations have been proposed: Its direct solution by the Full Approximation scheme , the Rayleigh Quotient Multigrid method  and a scheme where the the solution of the linear system of equations arising from the preconditioning step is performed by multigrid . The linear system solved in the preconditioning steps is not Schrödinger’s equation, but Poisson’s equation. The reason for this is that it is too difficult to find coarse grained representations of the Hamiltonian operator if nonlocal pseudopotentials are used and the deferred defect correction scheme  justifies the replacement of the Hamiltonian by the Laplacian. Thus, it turns out that in all cases it is Poisson’s equation that has to be solved using multigrid.

Another very promising multiscale approach to electronic structure calculations is the use of wavelets as basis sets . As with any large systematic basis set it is also in the wavelet context very important to use efficient preconditioning schemes. Diagonal preconditioning  is most widely used. However it would be desirable to have at our disposal more powerful preconditioning schemes. Using multigrid ideas for preconditioning has already been proposed in the context of interpolating wavelets . Our discussion of prolongation and restriction schemes based on wavelet theory will show that the class of schemes based on interpolating wavelets is not optimal.

Because of its central importance for atomistic simulations and because it is a prototype equation for the study of new methods we will from now on consider only Poisson’s equation

$$^2V(𝐫)=4\pi \rho (𝐫)$$
(1)
The solution of the differential equation Eq. 1 can be written as an integral equation

$$V(𝐫)=\frac{\rho (𝐫^{})}{|𝐫𝐫^{}|}$$
(2)
Gravitational problems are based on exactly the same equations as the electrostatic problem, but we will use in this article the language of electrostatics, i.e. we will refer to $`\rho (𝐫)`$ as a charge density. The most efficient numerical approaches for the solution of electrostatic problems are based on Eq 1 rather than Eq. 2. However preconditioning steps found in these methods can be considered as approximate solutions of Eq. 2. The fact that the Green’s function $`\frac{1}{|𝐫𝐫^{}|}`$ is of long range makes the numerical solution of Poisson’s equation difficult, since it implies that a charge density at a point $`𝐫^{}`$ will have an non-negligible influence on the potential $`V(𝐫)`$ at a point $`𝐫`$ far away. A naive implementation of Eq. 2 would therefore have a quadratic scaling. It comes however to our help, that the potential arising from a charge distribution far away is slowly varying and does not depend on the details of the charge distribution. All efficient algorithms for solving electrostatic problems are therefore based on a hierarchical multiscale treatment. On the short length scales the rapid variations of the potential due to the exact charge distribution of close by sources of charge are treated, on the large length scales the slow variation due to some smoothed charge distribution of far sources is accounted for. Since the number of degrees of freedom decreases rapidly with increasing length scales, one can obtain algorithms with linear or nearly linear scaling. In the following, we will briefly summarize how this hierarchical treatment is implemented in the standard algorithms

* Fourier Analysis:
If the charge density is written in its Fourier representation

$$\rho (𝐫)=\underset{𝐊}{}\rho _𝐊e^{i\mathrm{𝐊𝐫}}$$

the different length scales that are in this case given by $`\lambda =\frac{2\pi }{K}`$ decouple entirely and the Fourier representation of the potential is given by

$$V(𝐫)=\underset{𝐊}{}\frac{\rho _𝐊}{K^2}e^{i\mathrm{𝐊𝐫}}$$
(3)
The Fourier analysis of the real space charge density necessary to obtain its Fourier components $`\rho _𝐊`$ and the synthesis of the potential in real space from its Fourier components given by Eq. 3 can be done with Fast Fourier methods at a cost of $`Nlog_2(N)`$ where N is the number of grid points. The solution of Poisson’s equation in a plane wave is thus a divide and conquer approach where the division is into the single Fourier components.
* Multigrid methods (MG):
Trying to solve Poisson’s equation by any relaxation or iterative method (such as conjugate gradient) on the fine grid on which one finally wants to have the solution leads to a number of iterations that increases strongly with the size of the grid. The reason for this is that on a grid with a given spacing $`h`$ one can efficiently treat Fourier components with a wavelength $`\lambda =\frac{2\pi }{K}`$ that is comparable to the the grid spacing $`h`$, but the longer wavelength Fourier components converge very slowly. This increase in the number of iterations prevents a straightforward linear scaling solution of Eq. 1. In the multigrid method, pioneered by A. Brandt , one is therefore introducing a hierarchy of grids with a grid spacing that is increasing by a factor of two on each hierarchic level. In contrast to the Fourier method where the charge and the potential are directly decomposed into components characterized by a certain length scale, it is the residue that is passed from the fine grids to the coarse grids in the MG method. The residue corresponds to the charge that would give rise to a potential that is the difference between the exact potential and the approximate potential at the current stage of the iteration.

The solution of partial differential equations in a wavelet basis is typically done by preconditioned iterative techniques . The diagonal preconditioning approach, that is based on well established plane wave techniques, will be presented in the next section. The section after the next will introduce multigrid for Poisson’s equation in a wavelet basis. Even though the fundamental similarities between wavelet and multigrid schemes have been recognized by many workers (such as in ref. ) this sections contains to the best of our knowledge the first thorough discussion of how both methods can profit from each other.

Within wavelet theory  one has two possible representations of a function $`f(x)`$, a scaling function representation

$$f(x)=\underset{j}{}s_j^{Lmax}\varphi _j^{Lmax}(x)$$
(4)
and a wavelet representation.

$$f(x)=\underset{j}{}s_j^{Lmin}\varphi _j^{Lmin}(x)+\underset{l=Lmin}{\overset{Lmax}{}}\underset{j}{}d_j^l\psi _j^l(x).$$
(5)
In contrast to the scaling function representation, the wavelet representation is a hierarchic representation. The wavelet at the hierarchic level $`l`$ is related to the mother wavelet $`\psi `$ by

$$\psi _i^l(x)=\sqrt{2}^l\psi (2^lxi)$$
(6)
The characteristic length scale of a wavelet at resolution level $`l`$ is therefore proportional to $`2^l`$. A wavelet on a certain level $`l`$ is a linear combination of scaling functions at the higher resolution level $`l+1`$

$$\psi _i^l(x)=\underset{j=m}{\overset{m}{}}g_j\varphi _{2i+j}^{l+1}(x)$$
(7)
Scaling functions at adjacent resolution levels are related by a similar refinement relation

$$\varphi _i^l(x)=\underset{j=m}{\overset{m}{}}h_j\varphi _{2i+j}^{l+1}(x)$$
(8)
and hence also any wavelet at a resolution level $`l`$ is a linear combination of the highest resolution scaling functions. The so-called fast wavelet transform allows us to transform back and forth between a scaling function and a wavelet representation.

Let us now introduce wavelet representations of the potential and the charge density

$`V(x)`$ $`=`$ $`{\displaystyle \underset{j}{}}V_j^{Lmin}\varphi _j^{Lmin}(x)+{\displaystyle \underset{l=Lmin}{\overset{Lmax}{}}}{\displaystyle \underset{j}{}}V_j^l\psi _j^l(x)`$ (9)
$`\rho (x)`$ $`=`$ $`{\displaystyle \underset{j}{}}\rho _j^{Lmin}\varphi _j^{Lmin}(x)+{\displaystyle \underset{l=Lmin}{\overset{Lmax}{}}}{\displaystyle \underset{j}{}}\rho _j^l\psi _j^l(x)`$ (10)
Different levels do not completely decouple, i.e the components on level $`l`$, $`V_j^l`$, of the exact overall solution do not satisfy the single level Poisson equation

$$^2\left(\underset{j}{}V_j^l\psi _j^l(x)\right)4\pi \left(\underset{j}{}\rho _j^l\psi _j^l(x)\right)$$
(11)
within the chosen discretization scheme. This is due to the fact that the wavelets are not perfectly localized in Fourier space, i.e. many frequencies are necessary to synthesize a wavelet. However the amplitude of all these frequencies is clearly peaked at a nonzero characteristic frequency for any wavelet with at least one vanishing moment. From the scaling property (Eq. 6) it follows, that the frequency at which the peak occurs changes by a factor of two on neighboring resolution grids. This suggest that the coupling between the different resolution levels is weak.

In the preceding paragraph we presented the mathematical framework only for the one-dimensional case. The generalization to the 3-dim case is straightforward by using tensor products . Also in the rest of the paper only the one-dimensional form of the mathematical formulas will presented for reasons of simplicity. It has to be stressed however that all the numerical results were obtained for the three-dimensional case and with periodic boundary conditions.

## 3 The diagonal preconditioning approach for wavelets

Preconditioning requires finding a matrix with a simple structure that has eigenvalues and eigenvectors that are similar to the ones of the matrix in question . The structure has to be simple in the sense that it allows us to calculate the inverse easily. The simplest and most widely used structure in this respect is the structure of a diagonal matrix. As will be shown a diagonal preconditioning matrix can be found in a wavelet basis set and preconditioned conjugate gradient type methods are then a possible method for the solution of Poisson’s equation expressed in differential form (Eq.1).

As one adds successive levels of wavelets to the basis set the largest eigenvalue grows by a factor of 4 for each level. This can easily be understood from Fourier analysis. As one increases the resolution by a factor of 2 (i.e. increases the largest Fourier vector $`k_{max}`$ by a factor of 2) the largest eigenvalue increases by a factor of 4. This basic scaling property of the eigenvalue spectrum can easily be modeled by a diagonal matrix, where all the diagonal elements are all equal on one resolution level, but increase by a factor of 4 as one goes to a higher resolution level. It is of course clear that all the details of the true spectrum are not reproduced by this approximations. The true spectrum consists of a large number of moderately degenerate eigenvalues. The spectrum of the approximate matrix consists of a few highly degenerate eigenvalues where each eigenvalue has all the scaling functions of one resolution level as its eigenfunctions. The true eigenfunctions are of course mixtures of scaling functions on different resolution levels, but if the wavelet family is well localized in Fourier space the contributions from neighboring resolution levels are weak. The localization in Fourier space increases with the number of vanishing moments  and therefore this diagonal preconditioning method works for instance much better for lifted interpolating wavelets  than for ordinary interpolating  wavelets .

Another line of arguments, that shows the weakness of the diagonal preconditioning, is the following. The preconditioning matrix can also be considered to be the diagonal part of the matrix representing the Laplacian . Since the diagonal elements increase by a factor of 4 on each higher resolution level,

$`{\displaystyle \stackrel{~}{\psi }_0^{l+1}(x)\frac{^2}{x^2}\psi _0^{l+1}(x)𝑑x}`$ $`=`$ $`2{\displaystyle \stackrel{~}{\psi }_0^l(2x)\frac{^2}{x^2}\psi _0^l(2x)𝑑x}`$
$`=`$ $`4{\displaystyle \stackrel{~}{\psi }_j^l(x)\frac{^2}{x^2}\psi _j^l(x)𝑑x}`$
the spectrum of the matrix has again the correct scaling properties. Eq. 3 involves both the wavelets $`\psi `$ and their duals $`\stackrel{~}{\psi }`$ since it was written down for the most general context of a Petrov-Galerkin approach in a biorthogonal basis. In a pure Galerkin context or for orthogonal wavelet families the $`\stackrel{~}{\psi }`$ have to be replaced by the $`\psi `$’s. In the 3-dimensional case there are three different types of wavelets (products of 2 scaling functions and 1 wavelet, products of 1 scaling function and 2 wavelets and products of 3 wavelets). Each type of wavelet gives rise to a different diagonal element, but again all these diagonal elements differ by a factor of 4 on different resolution levels. Because of the weak coupling between different resolution levels discussed above, we expect the matrix elements of the Laplacian involving wavelets at two different resolutions levels to be small. The numerical examination of the matrix elements (Fig. 1) confirms this guess. It also shows that within one resolution level the amplitude of the matrix elements decays rapidly with respect to the distance of the two wavelets and is zero as soon as they do not any more overlap. Nevertheless, some matrix elements coupling nearest neighbor wavelets are not much smaller than the diagonal elements. One also finds a few matrix elements between different resolution levels that are less than one order of magnitude smaller than the one within a single resolution level.

The fact that all off-diagonal matrix elements are neglected in current precondition schemes explains their relatively slow convergence. It amounts to finding an approximate Greens function that is diagonal in a wavelet representation. This is obviously a rather poor approximation. Let us nevertheless stress that this diagonal matrix obtained by inverting a diagonal approximation to the Laplacian is a much more reasonable approximation for preconditioning purposes than the diagonal part of the Greens function. The diagonal part of the Greens function has actually completely different scaling properties. The elements increase by a factor of 2 as one goes to higher resolution levels instead of decreasing by a factor of 4. The multigrid methods to be discussed later include also in an approximative way through Gauss-Seidel relaxations this off-diagonal coupling within each resolution block as well as the coupling between the different resolutions levels.

In the following we will present some numerical results for the solution of the 3-dimensional Poisson equation in a wavelet basis using the diagonal preconditioning approach. All the methods presented in this paper will have the property that the convergence rate is independent of the grid size. We have chosen $`64^3`$ grids for all the numerical examples. The fact that the number of iterations necessary to reach a certain target accuracy is independent of the system size together with the fact that a single iteration involves a cost that is linear with respect to the number of grid points ensures that the Poisson’s equation can be solved with overall linear scaling. Whereas we use here only simple equidistant grids, this linear scaling has already been demonstrated with highly adaptive grids in problems that involve many different length scales .

The preconditioning step using simply the diagonal is given by

$$\mathrm{\Delta }V_j^l=const\mathrm{\hspace{0.25em}\hspace{0.25em}4}^l\mathrm{\Delta }\rho _j^l$$
(13)
In analogy to Eq. 9,10, the $`\mathrm{\Delta }\rho _j^l`$’s are the wavelet coefficients on the $`l`$-th resolution level of the residue $`\mathrm{\Delta }\rho (𝐫)=^2\stackrel{~}{V}(𝐫)+4\pi \rho (𝐫)`$ in a wavelet representation. $`\stackrel{~}{V}(𝐫)`$ is the approximate solution at a certain iteration of the solution process. The preconditioned residue $`\mathrm{\Delta }V`$ is then used to update the approximate potential $`\stackrel{~}{V}`$. In the case of the preconditioned steepest descent method used here this update simply reads

$$\stackrel{~}{V}\stackrel{~}{V}+\alpha \mathrm{\Delta }V$$
(14)
where $`\alpha `$ is an appropriate step size. As discussed above the constant in Eq. 13 depends in the three dimensional case upon which type of wavelet is implied since it is the inverse of the Laplace matrix element between two wavelets of this type.

Fig. 2 shows numerical results for several wavelet families. The slow convergence of the interpolating wavelets is due to the fact that they have a non-vanishing average and therefore a non-vanishing zero Fourier component . Hence they are all localized in Fourier space at the origin instead of being localized around a non-zero frequency. This deficiency can be eliminated by lifting. The Fourier power spectrum of the lifted wavelets tends to zero at the origin with zero slope for the family with two vanishing moments considered here. Lifting the wavelet twice leads to 4 vanishing moments and a even better localization in Fourier space. The improvement in the convergence rate is however only marginal. The higher 8-th order lifted interpolating wavelet is smoother than its 6-th order counterpart and hence better localized in the high frequency part. This also leads to a slightly faster convergence.

Combining the diagonal preconditioning (Eq. 13) with a FGMRES convergence accelerator  instead of using it within a steepest descent minimization gives a significantly faster convergence. The number of iterations can nearly be cut into half as shown in Fig 2.

Up to now we have only considered the case where the elements of the matrix representing the Laplacian were calculated within the same wavelet family that was used to analyze the residue by wavelet transformations to do the preconditioning step. More general schemes can however be implemented. It is not even necessary that the calculation of the Laplacian matrix elements is done in a wavelet basis. One can instead use simple second order finite differences, which in the one-dimensional case are given by

$$\frac{1}{h^2}(V_{i1}+2V_iV_{i+1}),$$
(15)
or some higher order finite differences for the calculation of the matrix elements. The scaling relation Eq. 3 does not any more hold exactly, but it is fulfilled approximately and the schemes works as well as in the pure wavelet case as is shown in Fig. 3.

## 4 The MG approach for wavelets

The aim of this part of the article is twofold. One aspect is how to speed up the convergence of the solution process for Poisson’s equation expressed in a wavelet basis set compared to the diagonal preconditioning approach. The other aspect is how to accelerate multigrid schemes by incorporating wavelet concepts. The part therefore begins with a brief review of the multigrid method.

Fig. 4 schematically shows the algorithm of a standard multigrid V cycle . Even though the scheme is valid in any dimension, a two dimensional setting is suggested by the figure, since the data is represented as squares. Since less data is available on the coarse grids, the squares holding the coarse grid data are increasingly smaller. It has to be stressed that the remarks of the end of the first part remain valid and that in particular all the numerical calculations are 3-dim calculations. The upper half of the figure shows the first part of the V cycle where one goes from the finest grid to the coarsest grid and the lower half the second part where one goes back to the finest grid.

In the first part of the V cycle the potential on all hierarchic grids is improved by a standard red-black Gauss-Seidel relaxation denoted by GS. The GS relaxation reduces the error components of wavelengths $`\lambda `$ that are comparable to the grid spacing $`h`$ very efficiently. In the 3-dimensional case we are considering here, the smoothing factor is .445 (page 74 of ref ). Since we use 2 GS relaxations roughly one quarter of the error around the wavelength $`h`$ survives the relaxations on each level. As a consequence the residue $`\mathrm{\Delta }\rho `$ contains mainly longer wavelengths which then in turn are again efficiently eliminated by the GS relaxations on the coarser grids. Nevertheless, the remaining quarter of the shorter wavelengths surviving the relaxations on the finer grid pollutes the residue on the coarser grid through aliasing effects. Aliasing pollution means that even if the residue on the finer grid would contain only wavelengths around $`h`$ (and in particular no wavelength around $`2h`$) the restricted quantity would not be identically zero.

In the second part of the V cycle the solutions obtained by relaxation on the coarse grid are prolongated to the fine grids and added to the existing solutions on each level. Aliasing pollution is again present in the prolongation procedure. Due to the accumulated aliasing errors 2 GS relaxations are again done on each level before proceeding to the next finer level.

To a first approximation the different representations of $`\rho `$ at the top of Fig. 4 represent Fourier filtered versions of the real space data set $`\rho `$ on the finest grid. The large data set contains all the Fourier components, while the smaller data sets contain only lower and lower frequency parts of $`\rho `$. In the 1-dimensional case only the lower half of the spectrum is still dominating when going to the coarse grid, in the 3-dimensional case it is only one eight of the spectrum. Because of the various aliasing errors described above the Fourier decomposition is however not perfect. Obviously it would be desirable to make this Fourier decomposition as perfect as possible. In the absence of aliasing errors, the GS relaxations would not have to deal with any Fourier components spilling over from higher or lower resolution grids.

Let us now postulate ideal restriction and prolongation operators and discuss their properties. As follows from the discussion above, they should provide for a dyadic decomposition of the Fourier space. Consequently the restriction operator would have to be a perfect low pass filter for the lower half of the spectrum (in the 1-dim case, in the 3-dim case only 1/8 would survive). We will refer to this property in following as frequency separation property. The degree of perfectness can be quantified by the $`k`$ dependent function

$$S^l(k)=\sqrt{\underset{i}{}(s_i^l)^2/N}$$
(16)
where $`N`$ is the number of grid points on resolution level $`l`$. The $`k`$ dependence enters through the requirement that the signal $`s_i^{Lmax}`$ on the finest resolution level is a pure harmonic,

$$s_i^{Lmax}=\mathrm{exp}(I2\pi ki/N)$$
(17)
The function $`S^l(k)`$ for a perfect restriction operator is shown in Fig 5. Such ideal grid transfer operators have to satisfy a second property, that will be baptized the identity property. The prolongation operator has to bring back exactly onto the finer grid the long wavelength associated with the coarser grid. This can only be true if prolongation followed by a restriction gives the identity. A third desirable property would be that the coarse charge density represents as faithfully as possible the significant features of the original charge density. In particular the coarse charge density should have the same multipoles and most importantly the same monopole. The conservation of the monopole just means that the total charge is identical on all grid levels. This third property will be called the multipole conservation property in the following.

With the ideal grid transfer operators, the solution of Poisson’s equation would be a divide and conquer approach in Fourier space and it could be done with a single V cycle with a moderate number of GS relaxations on each resolution level. In contrast to a solution in a plane wave basis the division would not be into single Fourier components but into dyadic parts of the Fourier spectrum. For the case of our postulated ideal grid transfer operators it also would not matter whether the GS relaxations are applied when going up or going down, only the total number of GS relaxations on each level would count.

To establish the relation between multigrid grid transfer operators and wavelet theory, let us point out a formal similarity. For vanishing $`d`$ coefficients, the wavelet analysis step is given by (Eq. 26 of ref. )

$$s_i^{2h}=\underset{j=m}{\overset{m}{}}\stackrel{~}{h}_js_{j+2i}^h$$
(18)
and is formally identical to a restriction operation. A wavelet synthesis step for the $`s`$ coefficients is given by (formula 27 of ref. )

$`s_{2i}^h`$ $`=`$ $`{\displaystyle \underset{j=m/2}{\overset{m/2}{}}}h_{2j}s_{ij}^{2h}`$ (19)
$`s_{2i+1}^h`$ $`=`$ $`{\displaystyle \underset{j=m/2}{\overset{m/2}{}}}h_{2j+1}s_{ij}^{2h}.`$
and is formally identical to a prolongation operation. One can now for example easily see that the injection scheme for the restriction and linear interpolation for the prolongation part corresponds to a wavelet analysis and synthesis steps for 1st order interpolating wavelets. Using the values of the filters $`\stackrel{~}{h}`$ and $`h`$ for interpolating wavelets we obtain

$$s_i^{2h}=s_{2i}^h$$
(20)
and

$`s_{2i}^h`$ $`=`$ $`s_i^{2h}`$ (21)
$`s_{2i+1}^h`$ $`=`$ $`{\displaystyle \frac{1}{2}}s_i^{2h}+{\displaystyle \frac{1}{2}}s_{i+1}^{2h}.`$
which is the standard injection and interpolation. As a consequence of the fact that it can be considered as a wavelet forward and backward transformation, the combination of injection and interpolation satisfy the identity property of our ideal grid transfer operator pair, namely that applying a restriction onto a prolongation gives the identity.

Usually injection is replaced by the full weighting scheme,

$$s_i^{2h}=\frac{1}{4}s_{2i1}^h+\frac{1}{2}s_{2i}^h+\frac{1}{4}s_{2i+1}^h.$$
(22)
This scheme has the advantage that it conserves averages, i.e it satisfies the monopole part of the multipole conservation property of an ideal restriction operator. Applying it to a charge density thus ensures that the total charge is the same on any grid level. Trying to put the full weighting scheme into the wavelet theory framework gives a filter $`\stackrel{~}{h}`$ with nonzero values of $`\stackrel{~}{h}_1=\frac{1}{4}`$, $`\stackrel{~}{h}_0=\frac{1}{2}`$, $`\stackrel{~}{h}_1=\frac{1}{4}`$ This filter $`\stackrel{~}{h}`$ does not satisfy the orthogonality relations of wavelet theory (formula 8 of ref. ) with the $`h`$ filter corresponding to linear interpolation. Hence a prolongation followed by a restriction does not give the identity.

A pair of similar restriction and prolongation operators that conserve averages can however be derived from wavelet theory. Instead of using interpolating wavelets we have to use lifted interpolating wavelets . In this way we can obtain both properties, average conservation and the identity for a prolongation restriction sequence. Using the filters derived in ref.  we obtain

$$s_i^{2h}=\frac{1}{8}s_{2i2}^h+\frac{1}{4}s_{2i1}^h+\frac{3}{4}s_{2i}^h+\frac{1}{4}s_{2i+1}^h\frac{1}{8}s_{2i+2}^h$$
(23)
$`s_{2i}^h`$ $`=`$ $`s_i^{2h}`$ (24)
$`s_{2i+1}^h`$ $`=`$ $`{\displaystyle \frac{1}{2}}s_i^{2h}+{\displaystyle \frac{1}{2}}s_{i+1}^{2h}.`$
Let us finally discuss the first property of our postulated ideal restriction operator, namely that it is a perfect low pass filter. Obviously any finite length filter can only be an approximate perfect low pass filter. Fig 6 shows the $`S`$ function for several grid transfer operators. One clearly sees that the Full Weighting operator is a poor low pass filter, the filter derived from second degree lifted wavelets is already better and the filters obtained from 10th degree Daubechies wavelets and twofold lifted 6th order interpolating wavelets are best. The Daubechies 6th degree filter is intermediate and of nearly identical quality as the one that is a mixture of the Full Weighting and second degree lifted wavelet filters. In contrast to the former, the later does however not fulfill the identity property.

The degree of perfectness for the frequency separation is also related to the multipole conservation property of our postulated ideal grid transfer operators. As one sees from Fig. 6 filters which correspond to wavelet families with many vanishing are closer to being ideal for frequency separation than those with few vanishing moments. At the same time the number of vanishing moments determines how many multipoles are conserved when the charge density is brought onto the coarser grids.

The right panel of Fig. 7 shows the convergence rate of a sequence of V cycles for the full weighting/interpolation (Eq. 22,21) scheme and various wavelet based schemes, namely the scheme obtained from second order lifted wavelets (Eq. 23,24), the corresponding scheme, but obtained from twofold lifted 6-th order wavelets as well as schemes obtained from 6th and 10th order Daubechies wavelets. The numerical values for the filters are listed in the Appendix. One can observe a clear correlation between the convergence rate and the the degree of perfectness of the $`S`$ function. A high degree of perfectness is particularly useful in connection with high order discretizations of the Laplacian. Most of the filters of the grid transfer operators are longer than the standard Full Weighting filter, which just has 3 elements. The lifted 2nd order interpolating wavelet restriction filter has for instance 5 elements and the 6-th degree Daubechies filter 6 elements. This does however not lead to a substantial increase of the CPU time. This comes from the fact that on modern computers the transfer of the data into the cache is the most time consuming part. How many numerical operations are then performed on these data residing in cache has only a minor influence on the timing. The new wavelet based schemes for restriction and prolongation are therefore more efficient than the Full Weighting scheme, both for finite difference discretizations and scaling function basis sets. It is also obvious that the multigrid approach for scaling/wavelet function basis sets is more efficient than the diagonal preconditioning approach.

The identity property for a restriction prolongation operator pair was only necessary for the case of operators where the restriction part is a perfect low pass filter. One might therefore wonder how useful it is in the context of the only nearly perfect filters. The numerical experience suggests that it is nevertheless a useful property. One can for example compare the convergence rates using either the 6-th order Daubechies filters or the filter that is the average of Full Weighting and lifted 2nd order wavelet filters. Fig. 6 shows that their restriction parts have very similar $`S`$ functions. Nevertheless we always found a better convergence rate with the Daubechies filter which satisfies the identity property.

For the 2-dimensional Poisson equation it has been shown, that the convergence rate compared to the standard Full Weighting scheme can be improved by tailoring grid transfer operators for the relaxation scheme used . The theoretical foundations for this is furnished by local Fourier analysis . The same approach could certainly also be used for the 3-dimensional case considered here. It is to be expected that the grid transfer operators found by such an optimization would be very close to the ones that we have obtained from wavelet theory.

The main justification for the relaxations in the upper part of the traditional multigrid algorithm shown in Fig. 4 is to eliminate the high frequencies. This can however be done directly by fast wavelet transformations based on wavelets that have good localization properties in frequency space such as lifted interpolating wavelets. As a consequence the traditional multigrid algorithms can be simplified considerably as shown in Fig. 8. Using wavelet based restriction and prolongation operators we can completely eliminate the GS relaxation in the first part of the V cycle where we go from the fine grid to the coarsest grid. We baptize such a simplified V cycle a halfway V cycle. The numerical results, obtained with the halfway V cycle, shown in the right hand plots of Fig. 7, demonstrate that the convergence is slightly faster than for the traditional multigrid algorithm based on the same restriction and prolongation scheme. In addition one step is faster. It is not necessary to calculate the residue after the GS relaxations. Otherwise the number of GS relaxations and restrictions/prolongations is identical in the full and halfway V cycle. On purpose no CPU times are given in this context because optimization of certain routines  can entirely change these timings. Because the residue is never calculated in the halfway V cycle, the memory requirements are also reduced.

The number of GS relaxations in the halfway V cycle was chosen to be 4 in order to allow for an unbiased comparison with the traditional V cycle where also 4 GS relaxations were done on the finest grid level. For optimal overall efficiency putting the number of GS relaxation to 3 is usually best, with the values of 2 and 4 leading to a modest increase in the computing time. The convergence rate of halfway V cycles as a function of the number of GS relaxations on the finest grid level is shown in Fig 9.

In all the previous examples we specified the number of GS relaxations on the finest grid level. On the coarser grid levels the number of iterations was allowed to increase by a factor of two per grid level. In this way it was practically always possible to find the exact solution on the most coarse grid. In addition we found that this trick slightly reduces the number of iterations and the total CPU time. The overall behavior of all the different methods is however identical when the number of GS relaxation is constant on each grid level.

## 5 Conclusions

Our results demonstrate that halfway V cycles with the restriction and prolongation steps based on wavelet theory are the most efficient approach for the solution of the 3-dimensional Poisson’s equation. It is most efficient both for finite difference discretizations and for the case where scaling functions or wavelets are used as basis functions. We expect that the approach should also be the most efficient one in connection with finite elements. It is essential that the wavelet family used for the derivation of the restriction and prolongation schemes has at least one vanishing moment and conserves thus average quantities on the various grid levels. Wavelet families with more vanishing moments do not lead an appreciable increase of the convergence rate compared to the case of one vanishing moment for low order discretizations of Poisson’s equation, but lead a modest further increase for high order discretizations. In the case where a wavelet family was used to discretize the Laplace operator, it is best to use the same wavelet family to construct the grid transfer operators. In addition to increased efficiency of the proposed halfway V cycle in terms of the CPU time, it is also simpler than the standard V cycle. This makes not only programming easier, but also reduces the memory requirements.

## 6 Acknowledgments

I was fortunate to have several discussion with Achi Brandt about this work. His great insight on multigrid methods, that he was willing to share with me, helped a lot to improve the manuscript. I thank him very much for his interest and advice.

## 7 Appendix

Filter for twofold lifted 6th order interpolating wavelets :
$`h_1`$=75/128,$`h_3`$=-25/256,$`h_5`$=3/256,
$`\stackrel{~}{h}_0`$=2721/4096, $`\stackrel{~}{h}_1`$=9/32, $`\stackrel{~}{h}_2`$=-243/2048, $`\stackrel{~}{h}_3`$=-1/32,
$`\stackrel{~}{h}_4`$=87/2048, $`\stackrel{~}{h}_6`$=-13/2048, $`\stackrel{~}{h}_8`$=3/8192 .
The values for negative indices follow from the symmetry $`h_i=h_i`$ and $`\stackrel{~}{h}_i=\stackrel{~}{h}_i`$.

Filters for 6-th order Daubechies wavelets :
$`h_2`$=0.3326705529500826159985d0, $`h_1`$=0.8068915093110925764944d0,
$`h_0`$=0.4598775021184915700951d0, $`h_1`$=-0.1350110200102545886963d0,
$`h_2`$=-0.0854412738820266616928d0, $`h_3`$= 0.0352262918857095366027d0 .
$`\stackrel{~}{h}_i=h_i`$.

Filters for 10-th order extremal Daubechies wavelets :
$`h_4`$=.1601023979741929d0, $`h_3`$=.6038292697971897d0,
$`h_2`$=.7243085284377729d0, $`h_1`$=.1384281459013207d0,
$`h_0`$=-.2422948870663820d0, $`h_1`$=-.0322448695846384d0,
$`h_2`$=.0775714938400457d0, $`h_3`$= -.0062414902127983d0,
$`h_4`$=-.0125807519990820d0, $`h_5`$=.0033357252854738d0 .
$`\stackrel{~}{h}_i=h_i`$.