text:
# Quantum Trajectories and Quantum Measurement Theory

## 1 Introduction

The concept of stochastic Schrödinger equations (sses)  or quantum trajectories  has generated an impressive degree of activity in the quantum optics community over the past few years 222Related ideas have also been developed by workers outside quantum optics. Space limitations prevent me from discussing these.. To one outside the field, it may be hard to synthesize all of this activity, not least because there are (in my mind) three distinct interpretations of quantum trajectories. I am here classifying authors by the ontological status which they ascribe to their quantum trajectories: not real; subjectively real; or objectively real. Alternatively, these three interpretations could be stated as quantum trajectories being respectively: a mathematical tool; an application of quantum measurement theory; or a theory of microscopic state reduction. In the bulk of this paper I am concerned with quantum trajectories as a part of quantum measurement theory (that is, with what I am calling the subjective interpretation of quantum trajectories). Note that by “measurement theory” I mean the body of mathematical techniques used to describe the process of measurement within the standard Copenhagen interpretation  of quantum mechanics. I am not concerned with metaphysical arguments regarding the adequacy of the Copenhagen viewpoint, nor with physical extensions to standard quantum theory. However, interpretational issues have been raised in the context of quantum trajectories, and the potential confusion regarding the reality of sses may obscure the results I wish to present. For these reasons, I will begin with a brief review of the interpretations of quantum trajectories. I am not aware of a comparison of this sort being in print.

### 1.1 Non-Real Quantum Trajectories

The “non-real” interpretation of quantum trajectories treats the trajectories simply as numerical tools for solving the master equation. This is essentially the attitude taken by Dalibard, Castin and Mølmer  who introduced the concept, although they do mention that their “Monte Carlo Wavefunction” approach also provides “new physical pictures”. Note the emphasis on wavefunctions in their terminology. The value of quantum trajectories for pure states as a computational tool is that it only takes $`2N2`$ real numbers to store a state vector in an $`N`$dimensional Hilbert space, compared to $`N^21`$ real numbers for a state matrix. Of course, this gain is offset by the fact that many trajectories are required to obtain a reliable ensemble average. However, for a large Hilbert space, quantum trajectories may offer a real advantage for numerical computation. One very practical application is computing the centre-of-mass motion of a fluorescent atom .

### 1.2 Subjectively Real Quantum Trajectories

In the body of this paper, quantum trajectories are attributed some degree of reality, in that they accurately model the evolution of a system conditioned on continuous observation. This can be called “subjective reality”, because the nature of the trajectories depend on the measurement scheme chosen by the observer. Once a particular measurement scheme has been chosen, the conditioned evolution of the system (as modeled by quantum trajectories) is objectively real in the sense that all observers will agree on the results. The subjectivity is present in the sense that it is not permissible to make a statement such as “an excited atom in free space will spontaneously emit a photon at a random point in time” unless there is a photodetector to record the emitted photon. A different detection scheme, such as heterodyne detection, will register a continuous flow of radiation from the atom, not a jump. In other words, the stochastic system evolution is only real in so far as the system is a subject of observation; it is not real for the system as an object in an observerless universe. As well as trajectories arising from different types of measurements, trajectories which do not preserve pure states (because of imperfect detection or bath preparation) also arise in this interpretation.

It is the subjective interpretation of quantum trajectories, in addition to their computational convenience, which is the motivation for the work by Carmichael and co-workers . Most of this work was numerical, but analytic results were obtained for some particular systems and particular measurement schemes . Both quantum jump (direct detection) and quantum diffusion (homodyne detection) trajectories were used. The measurement theory interpretation was also pointed out by Gardiner, Parkins, and Zoller  for the case of direct detection (which was formulated much earlier by Srinivas and Davies ).

There are a number of reasons for wishing to consider the measurement interpretation of quantum trajectories. Firstly, the simulated photocurrents can be used to calculate such quantities as the waiting time distribution for photodetections, which are quite difficult to calculate by traditional methods . Another example is the all-or-nothing fluorescence exhibited in the electron shelving of a single atom . The technique of linear quantum trajectories recently published by Goetsch and Graham  has opened the possibility for at least partly analytical solutions of measurement theory problems using quantum trajectories. This is the subject of the body of this paper, and I will thus defer further discussion until there.

A second motivation is that subjectively real quantum trajectories give valuable physical insight into irreversible quantum processes, different from that offered by the master equation. For example, a system as simple as a stationary two-level atom shows strikingly different quantum trajectory evolution depending on the chosen method of detection, and that this evolution sheds light on the respective measurement results . Another example is that coherences which are lost in the non-selective (master equation) picture are seen to be merely distributed among different quantum trajectories . These reconstructed coherences can restore noisy oscillatory dynamics to a system which would otherwise exhibit only damped oscillations .

Thirdly, the measured results can be used to select certain individual systems from the ensemble. In this way, it may be possible to see the coherences which are lost in the non-selective ensemble evolution . This could be looked upon as a very basic sort of feedback, in which the measured result is used to eliminate some members of the ensemble. The use of quantum trajectories in more active feedback schemes is a fourth motivation for subjectively real quantum trajectories, which has been pursued by myself and Milburn and others .

### 1.3 Objectively Real Quantum Trajectories

The third attitude towards quantum trajectories is a belief in their objective existence, independent of any measurement scheme. This belief is presumably founded on mistrust of the standard interpretations of quantum mechanics, and a desire to replace subjective collapses due to measurements with objectively real collapses. The roots of this approach are the dynamical state reduction models of Gisin , Diosi , Pearle  and Ghirardi-Rimini-Weber , which in turn hark back to the original objective interpretation of the wavefunction by Schrödinger . In this desire, proponents of objectively real quantum trajectory theory are akin to hidden variable theorists, but unlike the latter, they accept the wavefunction as a representation of reality. It is only when irreversible processes occur (described by a master equation) that stochasticity enters quantum mechanics, giving rise to quantum trajectories. Because the state vector is considered a basic element of reality in this interpretation, it only makes sense to consider quantum trajectories which preserve pure states. Thus it is that stochastic Schrödinger equations (as opposed to the more general stochastic master equations of subjective quantum trajectories) are essential only if one regards them as not real or objectively real.

At least three pairs of authors have published on this theme. The first is Teich and Mahler (TM). They published a single paper  presenting a theory of quantum jumps in fluorescent atoms. It has been shown to be nonphysical by myself and Milburn in reference . The second pair, Gisin and Percival (GP), have published many papers , also with other authors , and have received considerable attention. As pointed out by myself and Milburn, the quantum trajectory which GP consider can be given an interpretation in terms of unit-efficiency heterodyne detection (or ambi-quadrature homodyne detection) of the electromagnetic field radiated by the system . However GP attribute an objective reality to their sse, saying that it “represents the evolution of an individual quantum system in interaction with its environment”. By adding ad hoc irreversible processes to a model of a three level atom, they are able to induce their equation to exhibit jump-like behaviour , as seen in electron shelving experiments . In general, however, GP cannot associate their trajectories with the detection of photons for the very good reason that they actually correspond to heterodyne, not direct detection. A third pair, Breuer and Petruccione (BP), have very recently published two papers  in which they claim to show that the quantum jump unraveling of reference  is the sse for a Markovian open quantum system. Although it is in many ways the most natural unraveling, the quantum jump sse (which under some circumstances corresponds to direct detection ) does not have any fundamental privileged status. BP derive it by choosing a particular basis in which to diagonalize the bath, namely its energy eigenstates.

The basic problem with the objective interpretation of quantum trajectories lies with the origin of irreversible evolution. In standard quantum optics, the irreversibility of the master equation for a system is an approximation to the exact reversible evolution for the system and its environment . None of the above authors try to define in what contexts the master equation evolution is meant to be truly irreversible. If there were some fundamental physical mechanism by which irreversibility entered the world, for example at some sufficiently large scale, then it is conceivable that one of the three schemes outlined above could be correct 333Of course, these are not the only three conceivable dynamical state reduction schemes; two which have gained considerable attention are those of GRW  and Penrose .. The quantum measurement “problem” would then be solved, as measurement results would be determined by the irreducibly stochastic quantum evolution. Irrespective of this speculation, the models of TM, GP and BP cannot be objectively real for the systems to which the authors apply them 444These criticism now appears to be accepted by Percival, as evidenced by reference .. Irreversible stochastic evolution, in the manner of TM, GP or BP, does not take place within an individual atom. The entanglement between an atom and its outgoing field is a fact which is not represented by the models of TM, GP and BP. This entanglement is most simply shown by the fact that one can choose different detection schemes to observe the atom, as discussed in . Gisin and Percival’s scheme only gives the correct results for heterodyne detection, Breuer and Petruccione’s for direct detection, while Teich and Mahler’s scheme will not give the correct results for any detection scheme .

### 1.4 Outline

Having discussed the non-real and objectively real interpretations of quantum trajectories, the remainder of this paper is devoted to quantum trajectories as a part of quantum measurement theory. I begin with a discussion of traditional quantum measurement theory, and why this is inadequate to describe practical measurements. This leads to an outline of the theory of operations and effects, and the different ways in which probability can be incorporated in the theory. Sec. 3 introduces quantum trajectories as an application of quantum measurement theory to continuous observation. The description of such trajectories as linear stochastic Schrödinger equations is covered in Sec. 4. Sec. 5 presents the general solution for the linear sse for the case of an undriven system with a large local oscillator used in the detection process. This allows a new derivation of the operations and effects for completed homodyne and heterodyne measurements. In Sec. 6 I apply the results of Sec. 5 to a measurement scheme which has not previously been considered and which is probably intractable without linear quantum trajectories. The measurements which I am describing here are adaptive measurements, in which the past photocurrent is used to modify the present conditions of the measurement. Here I discuss the general approach, and present results for the simple but non-trivial example of adaptive phase measurements of a field containing at most one photon. Sec. 7 concludes.

## 2 Quantum Measurement Theory

Traditional quantum measurement theory is founded on projective measurements. Consider a measurement of fixed duration $`T`$ with a set of $`N`$ discrete possible measurement results labeled by $`r`$. With each result $`r`$ there is an associated projection operator $`\mathrm{\Pi }_r(T)`$, such that $`\mathrm{\Pi }_r(T)\mathrm{\Pi }_s(T)=\mathrm{\Pi }_r(T)\delta _{r,s}`$. If the initial state of the system (assumed pure for simplicity) is $`|\psi (t)`$, then the probability for obtaining the result $`r`$ is

$$P_r(T)=\psi (t)|\mathrm{\Pi }_r(T)|\psi (t).$$
(1)
This can also be written

$$P_r(T)=\stackrel{~}{\psi }_r(t+T)|\stackrel{~}{\psi }_r(t+T),$$
(2)
where the tilde denotes an unnormalized state vector

$$|\stackrel{~}{\psi }_r(t+T)=\mathrm{\Pi }_r(T)|\psi (t).$$
(3)
Furthermore, the final state of the system is simply this state, normalized if desired

$$|\psi _r(t+T)=\frac{|\stackrel{~}{\psi }_r(t+T)}{\sqrt{P_r(T)}}.$$
(4)
Although this theory may be fundamentally correct in some sense, it is inadequate to describe practical measurements such as those in quantum optics. That is because the experimenter never makes a direct measurement on the system of interest. Rather, the system of interest (such as an atom) interacts with its environment (the continuum of electromagnetic field modes), and the experimenter observes the radiated field. Of course one could argue that the experimenter does not observe the radiated field, rather that the field interacts with a photodetector, which triggers a current in a circuit, which is coupled to a CRO, which radiates more photons, which interact with the experimenter’s retina, and so on up the von-Neumann chain . The point is that at some stage one has to consider the measurement as being made. If one attempts to apply a projection postulate directly to the atom, one will obtain nonsensical predictions. However, it can be shown that assuming a projective measurement of the field will yield results negligibly different from those assuming a projective measurement at any later stage, because of the rapid decoherence of macroscopic objects . For this reason, it is sufficient to consider the field as being measured. Because the field has interacted with the system, their quantum states are entangled. If the initial state of the field is assumed known (which is not unreasonable in practice, because it is often close to the vacuum state), then the projective measurement of the field results in a measurement of the atom. Such a measurement however is not projective 555This is why the original quantum Zeno paradox of Misra and Sudarshan  can be avoided in practice.. Instead, we need a more general formalism to describe such measurements. This is called the theory of operations and effects , which I will now present in brief.

Say the initial system state vector is $`|\psi (t)`$ as above, and that the measurement results are again labeled by $`r`$. Each $`r`$ is again associated with an operator $`\mathrm{\Omega }_r(T)`$, but it need not be a projection operator. I will call these measurement operators. The effect  for the result $`r`$ is defined by

$$F_r(T)=\mathrm{\Omega }_r(T)^{}\mathrm{\Omega }_r(T).$$
(5)
This operator is used to generate the probability for the result via

$$P_r(T)=\psi (t)|F_r(T)|\psi (t).$$
(6)
Conservation of probability yields the single restriction which exists for the measurement operators , namely

$$\underset{r}{}F_r(T)=1.$$
(7)
The set of all effects $`\{F(r)\}`$ constitutes a positive-operator-valued measure (povm) on the space of results $`r`$ . The state of the system given the result $`r`$ is

$$|\psi _r(t+T)=\frac{\mathrm{\Omega }_r(T)|\psi (t)}{\sqrt{P_r(T)}}.$$
(8)
If one were doing only a single measurement, then the conditioned state $`|\psi _r`$ would be irrelevant. However one often wishes to consider a sequence of measurements, in which case the conditioned system state is vital. In terms of the state matrix $`\rho `$, which allows the possibility of mixed initial states, the conditioned state is

$$\rho _r(t+T)=\frac{𝒥[\mathrm{\Omega }_r(T)]\rho (t)}{P_r(T)},$$
(9)
where $`P_r(T)=\mathrm{Tr}[\rho (t)F_r(T)]`$ and for an arbitrary operator $`a`$,

$$𝒥[a]\rho =a\rho a^{}.$$
(10)
The superoperator $`𝒥[\mathrm{\Omega }_r(T)]`$ is known as the operation for $`r`$ . If the measurement were performed but the result ignored, the final state of the system would be

$$\rho (t+T)=\underset{r}{}P_r(T)\rho _r(t+T).$$
(11)
Consider a hypothetical experiment involving a sequence of measurements for which one knows the set of measurement operators $`\{\mathrm{\Omega }_r(T)\}`$. As a theorist, one may be interested in knowing a) the probability for obtaining the possible measurement sequences, and b) the state of the system during and at the end of the measurement, conditioned on the results obtained. If initially $`\rho (t)=|\psi (t)\psi (t)|`$ then the above formalism shows that one could obtain this information as follows. Start with $`|\psi (t)`$, and generate the probabilities for the results using equation (6). One would then choose results at random using these probabilities, and see what the new conditioned state of the system is using equation (8). This is the input into the next measurement. The probability for each sequence of results to actually occur is simply the probability that it was generated by the simulation procedure. The simulation is repeated as may times as necessary to obtain good statistics. Call this method A.

As in the case for projective measurements, this theory can be alternatively stated as follows. The final state of the system conditioned on the result $`r`$ is represented by an unnormalized state vector

$$|\stackrel{~}{\psi }_r(t+T)=\mathrm{\Omega }_r(T)|\psi (t).$$
(12)
The probability for this result having been obtained is then

$$P_r(T)=\stackrel{~}{\psi }_r(t+T)|\stackrel{~}{\psi }_r(t+T).$$
(13)
The spirit of this formulation is that one could simulate the measurement process by randomly generating a result $`r`$ with probability $`1/N`$, then calculating the appropriate unnormalized conditioned state. This is the input into the next measurement step. The actual probability assigned to each sequence of measurement results is the final norm squared of the state vector (13), where here $`T`$ would be the total time for the sequence of measurements. This simulation would be less efficient because all measurement results are given equal computational time, even those with vanishingly small probabilities. However it would have the advantage that probabilities need to be calculated only once (at the end). Call this method B.

In fact, it is possible to formulate the theory in any number of ways which include these two as special cases. Define an (in general unnormalized) conditioned state vector by

$$|\overline{\psi }_r(t+T)=\frac{\mathrm{\Omega }_r(T)|\psi (t)}{\sqrt{\mathrm{\Lambda }_r(T)}}.$$
(14)
Here the $`\mathrm{\Lambda }_r(T)`$ are positive numbers chosen such that

$$\underset{r}{}\mathrm{\Lambda }_r(T)=1.$$
(15)
They can thus be interpreted as probabilities in some sense, and I will call $`\mathrm{\Lambda }_r`$ the ostensible probability for the result $`r`$. The actual probability for getting result $`r`$ is

$$P_r(T)=\mathrm{\Lambda }_r(T)\overline{\psi }_r(t+T)|\overline{\psi }_r(t+T).$$
(16)
If the $`\mathrm{\Lambda }_r(T)`$ are chosen such that $`\mathrm{\Lambda }_r(T)`$ happens to equal $`\psi (t)|F_r(T)|\psi (t)`$, then this formulation is equivalent to the method A given above, with $`|\psi _r(t+T)=|\overline{\psi }_r(t+T)`$. If they are chosen such that $`\mathrm{\Lambda }_r(T)=1/N`$, then this formulation is equivalent to the method B, with $`|\stackrel{~}{\psi }_r(t+T)=|\overline{\psi }_r(t+T)\sqrt{N}`$. In general, one could do the simulation of the measurement by choosing the results $`r`$ with probability $`\mathrm{\Lambda }_r(T)`$, then calculating the conditioned state vector $`|\overline{\psi }_r(t+T)`$. Each result would then be weighted by $`\overline{\psi }_r(t+T)|\overline{\psi }_r(t+T)`$, where again $`T`$ would be the total time. Because of the probability of generation, this weighting gives the correct overall weighting to each member of the ensemble as given in equation (16). Call this method C.

This flexibility of formulation has been known for some time (see Goetsch and Graham  and references therein). As presented here it must appear rather pointless. However we are about to consider quantum trajectories, as an application of measurement theory for continuously (in time) monitored systems. In this context, it turns out that the ability to separate out the probability of each possible result into two components (a probability of generation $`\mathrm{\Lambda }_r(T)`$ and a weighting $`\overline{\psi }_r(t+T)|\overline{\psi }_r(t+T)`$), will be of enormous use in solving the quantum trajectories.

Before turning to the topic of quantum trajectories, there is one more useful fact which can be explained in terms of the general formalism. Recall that the final state of the system if all measurement results are ignored is given by equation (11). This can also be written as

$$\rho (t+T)=\underset{r}{}𝒥[\mathrm{\Omega }_r]\rho (t).$$
(17)
Now if we define a unitary rearrangement of of the measurement operators $`\mathrm{\Omega }_r`$ by

$$\mathrm{\Omega }_r^{}=\underset{s}{}U_{r,s}\mathrm{\Omega }_s,$$
(18)
where $`U`$ is a $`c`$number matrix satisfying $`_rU_{r,s}U_{r,q}^{}=\delta _{s,q}.`$ Then the unconditioned final state under the new measurement operators $`\{\mathrm{\Omega }_r^{}\}`$ is

$`{\displaystyle \underset{r}{}}\mathrm{\Omega }_r^{}\rho (t)\mathrm{\Omega }_{r}^{}{}_{}{}^{}`$ $`=`$ $`{\displaystyle \underset{r,s,q}{}}U_{r,s}\mathrm{\Omega }_s\rho (t)\mathrm{\Omega }_q^{}U_{r,q}^{}`$ (19)
$`=`$ $`{\displaystyle \underset{s,q}{}}\delta _{s,q}\mathrm{\Omega }_s\rho (t)\mathrm{\Omega }_q^{}={\displaystyle \underset{s}{}}\mathrm{\Omega }_s\rho (t)\mathrm{\Omega }_s^{}.`$ (20)
which is the same as that under the old measurement operators $`\{\mathrm{\Omega }_r\}`$ . That is to say, one can unitarily rearrange the measurement operators without changing the non-selective system evolution.

## 3 Quantum trajectories

One case of the above measurement theory which is of great importance is continuous observation, with an infinitesimal measurement interval $`T=dt`$. The form of such a measurement theory is restricted by the requirement that the non-selective evolution generated by the measurements should give a valid evolution equation for the state of the system. Explicitly, this non-selective evolution is

$$\rho (t+dt)=\underset{r}{}\mathrm{\Omega }_r(dt)\rho (t)\mathrm{\Omega }_r^{}(dt).$$
(21)
This is obviously Markovian , with the increment in the system state from time $`t`$ to time $`t+dt`$ depending only on the state of the system at time $`t`$. A Markovian evolution equation for the state matrix of a system is called a master equation. Physically, it can be derived from the interaction of the system with its environment, if the environment can be treated as a bath. That is to say, it is necessary for the environment to irreversibly carry away information from the system. Alternatively, one can view the environment as a continuous stream of independently prepared apparatuses which are wheeled up to the system, interact for an infinitesimal time $`dt`$, and are wheeled away again to be measured. Although this sounds unrealistic, the electric-dipole interaction between a system and the electromagnetic vacuum behaves like this to a good approximation .

The most general form of master equation is the so-called Lindblad form

$$\dot{\rho }=i[H,\rho ]+\underset{\mu }{}𝒟[c_\mu ]\rho .$$
(22)
Here $`H`$ is an Hermitian operator and $`𝒟`$ is a superoperator taking an operator argument (enclosed in square brackets). It is defined for an arbitrary operator $`r`$ by

$$𝒟[r]\rho r\rho r^{}\frac{1}{2}(r^{}r\rho +\rho r^{}r).$$
(23)
The operators $`c_\mu `$ are completely arbitrary.

Consider now a special case where there is but one Lindblad operator $`c`$ for the system, so that

$$\dot{\rho }=i[H,\rho ]+𝒟[c]\rho \rho .$$
(24)
Then it is possible to write the master equation as the nonselective evolution due to a measurement \[as in equation (21)\] with just two measurement operators $`\mathrm{\Omega }_r(dt)`$ by choosing

$`\mathrm{\Omega }_1(dt)`$ $`=`$ $`\sqrt{dt}c,`$ (25)
$`\mathrm{\Omega }_0(dt)`$ $`=`$ $`1\left(iH+\frac{1}{2}c^{}c\right)dt.`$ (26)
This is because

$$\rho (t+dt)=\underset{r=0,1}{}𝒥[\mathrm{\Omega }_r(dt)]\rho (t)=(1+dt)\rho (t),$$
(27)
where $``$ is as given in equation (24).

With the measurement operators $`\mathrm{\Omega }_0(dt),\mathrm{\Omega }_1(dt)`$, it is evident that the measurement record will be a point process , as I will now explain. For almost all infinitesimal time intervals, the measurement result is $`r=0`$, because $`P_0(dt)=1O(dt)`$. The result $`r=0`$ is thus regarded as a null result. In the case of no result, the system state changes infinitesimally, but not unitarily, via the operator $`\mathrm{\Omega }_0`$. At randomly determined (but not necessarily Poisson distributed) times, there is a result $`r=1`$, which I will call a detection. When this occurs, the system undergoes a finite evolution induced by the operator $`\mathrm{\Omega }_1`$. This change can validly be called a quantum jump, although it must be remembered that it represents a sudden change in the observer’s knowledge, not an objective physical event as in Bohr’s original conception . Real measurements which correspond approximately to this ideal measurement model are made routinely in experimental quantum optics. If $`c`$ is the lowering operator for the dipole of the quantum system, multiplied by the square root of the damping rate, then this theory describes the system evolution in terms of photodetections .

Now these measurement operators can be unitarily rearranged by the matrix

$$\left(\begin{array}{cc}U_{00}& U_{01}\\ U_{10}& U_{11}\end{array}\right)=\left(\begin{array}{cc}1\frac{1}{2}|\gamma |^2dt& \gamma \sqrt{dt}\\ \gamma ^{}\sqrt{dt}& 1\frac{1}{2}|\gamma |^2dt\end{array}\right).$$
(28)
This is equivalent to the transformation

$$cc+\gamma ;HHi\frac{1}{2}(\gamma ^{}c\gamma c^{}),$$
(29)
under which the master equation (24) is of course invariant. This transformation preserves the property that the collapse operator $`\mathrm{\Omega }_1`$ is infinitesimal, and the smooth evolution operator $`\mathrm{\Omega }_0`$ is close to unity. Explicitly,

$`\mathrm{\Omega }_1(dt)`$ $`=`$ $`\sqrt{dt}(c+\gamma ),`$ (30)
$`\mathrm{\Omega }_0(dt)`$ $`=`$ $`1dt\left[iH+\frac{1}{2}(c\gamma ^{}c^{}\gamma )+\frac{1}{2}(c^{}+\gamma ^{})(c+\gamma )\right].`$ (31)
Physically, it amounts to adding a coherent field to the output radiation field before it is detected by the system. This can be done using a low-reflectance beam-splitter. In the limit $`|\gamma |^2c^{}c`$, this procedure is known as homodyne detection.

The specification of the measurement operators is all that is necessary to define the measurement. However, it is useful to give a more explicit description of the evolution of the system under this scheme. This is done by introducing a stochastic increment $`dN(t)`$ which equals one if the result 1 is obtained in the time interval $`(t,t+dt)`$ and 0 otherwise. Formally, $`dN(t)`$ is defined by

$`dN(t)^2`$ $`=`$ $`dN(t),`$ (32)
$`\mathrm{E}[dN(t)]`$ $`=`$ $`\mathrm{\Omega }_1^{}(dt)\mathrm{\Omega }_1(dt)=dt\psi (t)|(c+\gamma )^{}(c+\gamma )|\psi (t),`$ (33)
where a classical expectation value is denoted by $`\mathrm{E}`$ and the quantum expectation value by angle brackets. From the measurement operators (30,31), the stochastic evolution equation for the conditioned state vector is

$`d|\psi (t)`$ $`=`$ $`[dN(t)({\displaystyle \frac{c+\gamma }{\sqrt{(c^{}+\gamma ^{})(c+\gamma )(t)}}}1)`$ (34)
$`+dt({\displaystyle \frac{c^{}c(t)}{2}}{\displaystyle \frac{c^{}c}{2}}+{\displaystyle \frac{c^{}\gamma +\gamma ^{}c(t)}{2}}\gamma ^{}ciH)]|\psi (t).`$
This gives the quantum trajectory for the system explicitly as a non-linear stochastic Schrödinger equation. Putting $`\gamma =0`$ gives the form generally used for numerical solutions of the master equation.

## 4 Linear Quantum Trajectories

The evolution described by equation (34) corresponds to the method A of Sec. 2. That is to say, the results $`dN(t)`$ are generated with the distribution they would have in actuality. No other weighting is necessary. In this case, method B would be outstandingly inefficient. It would imply either generating or not generating a photodetection with equal probability for every instant of time ($`\mathrm{E}[dN]=\frac{1}{2}`$). The norms of the vast majority of states would rapidly become vanishingly small, and their contribution to the ensemble would be negligible. The only states which would significantly contribute to the ensemble would be those in which photodetections were sparse, as they would be in actuality. However, there is a half-way house, corresponding to method C of Sec. 2. This would be achieved by keeping $`\mathrm{E}[dN]=O(dt)`$. Specifically, a convenient choice is

$$\mathrm{\Lambda }_1(dt)=\mathrm{E}[dN]=|\gamma |^2dt,$$
(35)
so that by necessity $`\mathrm{\Lambda }_0(dt)=1|\gamma |^2dt`$. Then the unnormalized conditioned system state would evolve via

$`|\overline{\psi }_1(t+dt)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_1(dt)|\overline{\psi }(t)}{\sqrt{\mathrm{\Lambda }_1(dt)}}}=\left(1+{\displaystyle \frac{c}{\gamma }}\right)|\psi (t)`$ (36)
$`|\overline{\psi }_0(t+dt)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_0(dt)|\overline{\psi }(t)}{\sqrt{\mathrm{\Lambda }_0(dt)}}}=\left[1dt\left(\frac{1}{2}c^{}c+c\gamma ^{}+iH\right)\right]|\psi (t)`$ (37)
The ostensible statistics for $`N(t)`$ are those of a Poisson process. The actual probability for a detection varies in time according to

$$P_1(dt)=\mathrm{\Lambda }_1(dt)\overline{\psi }_1(t+dt)|\overline{\psi }_1(t+dt)=dt\psi (t)|(c+\gamma )^{}(c+\gamma )|\psi (t),$$
(38)
as before.

Writing out this evolution explicitly, we have

$$d|\overline{\psi }(t)=\left[dN(t)\frac{c}{\gamma }dt\left(iH+\frac{1}{2}c^{}c+\gamma ^{}c\right)\right]|\overline{\psi }(t).$$
(39)
This is considerably simpler than equation (34), and, more importantly, it is linear in $`|\overline{\psi }(t)`$. This has been achieved by transferring some of the information about the probability of each measurement record into the norm of $`|\overline{\psi }(t)`$. The choice of $`|\gamma |^2dt`$ for $`\mathrm{\Lambda }_1(dt)`$ is obviously not a good one if the local oscillator amplitude $`\gamma `$ vanishes. In the remainder of this paper I am concerned with the opposite limit, $`|\gamma |\mathrm{}`$. Consider a short time $`\delta t1`$ in which the system changes negligibly, but the number of detections $`\delta N(t)|\gamma |^2\delta t`$ is very large. This allows the Poisson process $`\delta N(t)`$ to be approximated by a Gaussian process:

$$\delta N(t)|\gamma |^2\delta t+|\gamma |\delta W(t),$$
(40)
where $`\delta W(t)`$ is a Gaussian random variable of zero mean and variance $`\delta t`$. Since this time is infinitesimal as far as the system is concerned, we can substitute this into equation (39) to obtain

$$d|\overline{\psi }(t)=\left[dW(t)e^{i\mathrm{\Phi }}cdt\left(iH+\frac{1}{2}c^{}c\right)\right]|\overline{\psi }(t),$$
(41)
where $`\mathrm{\Phi }=\mathrm{arg}\gamma `$ and $`dW(t)`$ can be treated as an infinitesimal Wiener increment  satisfying

$$\mathrm{E}[dW(t)]=0,dW(t)^2=dt.$$
(42)
This equation has been previously derived by Goetsch and Graham, but in a different way. Their method involved considering ideal measurements of one quadrature of the radiated field, without involving photodetection or a local oscillator. As there is no known physical mechanism for undertaking such an ideal measurement, I believe that it is valuable to be able to derive the linear sse (41) as a limit of a physical detection scheme.

## 5 Solving the Linear Quantum Trajectory

It is now convenient to introduce a new notation. Let the recorded photocurrent be defined

$$I(t)=\underset{|\gamma |\mathrm{}}{lim}\frac{\delta N(t)|\gamma |^2\delta t}{|\gamma |\delta t}=\frac{dW(t)}{dt},$$
(43)
and let the complete record of the photocurrent from time 0 to just before time $`t`$ be denoted $`𝐈_{[0,t)}`$. Then as a reminder that $`|\overline{\psi }(t)`$ is conditioned on the photocurrent, denote it by $`|\overline{\psi }_t^𝐈`$. Goetsch and Graham  have solved the linear SSE (41) for the case of homodyne detection, where the local oscillator phase $`\mathrm{\Phi }`$ is held constant, and for the case of heterodyne detection where it varies rapidly: $`\mathrm{\Phi }(t)=\mathrm{\Phi }_0+\mathrm{\Delta }t`$. This was done for some particularly simple internal Hamiltonians $`H`$. Here I show that a solution may be found in a similar manner for any variation of $`\mathrm{\Phi }`$ with time. I consider only the case $`H=0`$, which amounts to a freely damped cavity. That is because I am interested in describing a single measurement of some property of the initial state of the cavity, in which all of the light is allowed to leak out and be detected.

Consider an initial coherent state $`|\psi _0=|\alpha `$ where $`a|\alpha =\alpha |\alpha `$, where $`a=c`$ is the annihilation operator for the damped mode, obeying $`[a,a^{}]=1`$. Now coherent states are unique with respect to the damping master equation (24) with $`H=0`$, in that they are the only states which remain pure under that evolution equation. Specifically, they have the solution $`|\psi _t=|\alpha e^{t/2}`$. This means that the solution to the linear SSE (39) must be of the form

$$|\overline{\psi }_t^𝐈=\overline{\psi }_t^𝐈|\alpha e^{t/2}.$$
(44)
Now using the Itô calculus, equation (41) can also be written as

$$|\overline{\psi }_{t+dt}^𝐈=\mathrm{exp}(\frac{1}{2}a^{}adt)\mathrm{exp}\left(dW(t)e^{i\mathrm{\Phi }(t)}a\frac{1}{2}e^{2i\mathrm{\Phi }(t)}a^2dt\right)|\overline{\psi }_t^𝐈.$$
(45)
Now using the fact that

$$\mathrm{exp}(\frac{1}{2}a^{}at)|\alpha =\mathrm{exp}\left(\frac{1}{2}|\alpha |^2(1e^t)\right)|\alpha e^t,$$
(46)
one has

$$\overline{\psi }_{t+dt}^𝐈=\mathrm{exp}\left(\frac{1}{2}|\alpha |^2e^tdt\frac{1}{2}e^{2i\mathrm{\Phi }(t)}\alpha ^2e^tdt+e^{i\mathrm{\Phi }(t)}\alpha e^{t/2}I(t)dt\right)\overline{\psi }_t^𝐈,$$
(47)
where it is to be remembered that $`dW(t)=I(t)dt`$. The total solution satisfying $`\overline{\psi }_0=1`$ is thus

$$\overline{\psi }_t^𝐈=\mathrm{exp}\left(\frac{1}{2}|\alpha |^2(1e^t)\frac{1}{2}\alpha ^2_0^te^{2i\mathrm{\Phi }(s)}e^s𝑑s+\alpha _0^te^{i\mathrm{\Phi }(s)}I(s)𝑑s\right).$$
(48)
Now it follows from this that the solution for the initial coherent state $`|\psi _0=|\alpha `$ can also be written

$$|\overline{\psi }_t^𝐈=\mathrm{exp}(\frac{1}{2}a^{}at)\mathrm{exp}\left(\frac{1}{2}S^{}(t)a^2+R^{}(t)a\right)|\psi _0,$$
(49)
where

$$R(t)=_0^te^{i\mathrm{\Phi }(s)}e^{s/2}I(s)𝑑s;S(t)=_0^te^{2i\mathrm{\Phi }(s)}e^s𝑑s.$$
(50)
That is to say, the solution at time $`t`$ depends not on the full photocurrent record $`𝐈_{[0,t)}`$ (which is a continuous infinity of real numbers), but rather on that record only through the two complex integrals (50) 666It might be thought that $`S(t)`$ does not even depend on the record $`𝐈_{[0,t)}`$, but it might if the local oscillator phase $`\mathrm{\Phi }(t)`$ is adjusted according to $`𝐈_{[0,t)}`$, so as to make an adaptive measurement as will be considered below.. Furthermore, because equation (39) is a linear equation by construction, the solution (49) is true for arbitrary initial conditions, not just for an initially coherent state. Thus any state conditioned on photodetection in the large local oscillator limit depends on the photocurrent only by the two integrals $`R(t)`$ and $`S(t)`$, regardless of how the local oscillator phase varies.

The measurement is complete in the limit $`t\mathrm{}`$, when there is no light in the cavity and hence no information left to be obtained. In that limit, $`lim_t\mathrm{}\mathrm{exp}(a^{}at)=|00|,`$ so that the final conditioned state $`|0`$ is obviously of no interest. What is of interest is the probability for having obtained the photocurrent record $`𝐈_{[0,t)}`$. As we have seen above, the only relevant aspects of this record are $`R(\mathrm{})`$, which I will denote $`A`$, and $`S(\mathrm{})`$, which I will denote $`B`$. From Sec. 3, the probability for obtaining the results $`A`$ and $`B`$ is

$$P(A,B)=P_0(A,B)\overline{\psi }_{\mathrm{}}^{A,B}|\overline{\psi }_{\mathrm{}}^{A,B},$$
(51)
where $`|\overline{\psi }_{\mathrm{}}^{A,B}`$ is the unnormalized conditioned state of equation (49) and $`P_0(A,B)`$ is the ostensible distribution for $`A`$ and $`B`$, given by

$$P_0(A,B)=𝑑𝐈_{[0,\mathrm{})}P_0\left(𝐈_{[0,\mathrm{})}\right)\delta ^{(2)}\left(AR(\mathrm{})\right)\delta ^{(2)}\left(BS(\mathrm{})\right).$$
(52)
Here $`R(\mathrm{})`$ and $`S(\mathrm{})`$ are the functionals of $`𝐈_{[0,\mathrm{})}`$ defined above, and $`P_0\left(𝐈_{[0,\mathrm{})}\right)`$ is the ostensible distribution for $`𝐈_{[0,\mathrm{})}`$, which as seen above, is simply that of Gaussian white noise.

From equation (49) we can rewrite the probability for obtaining results $`A,B`$ as

$$P(A,B)=\mathrm{Tr}[F(A,B)\rho _0],$$
(53)
where $`\rho _0`$ is the initial system state (which is allowed to be mixed), and $`F(A,B)`$ is the effect for the results $`A,B`$ defined by

$$F(A,B)=P_0(A,B)\overline{\mathrm{\Omega }}^{}(A,B)\overline{\mathrm{\Omega }}(A,B),$$
(54)
where the (unnormalized) measurement operators are

$$\overline{\mathrm{\Omega }}(A,B)=|00|\mathrm{exp}\left(\frac{1}{2}B^{}a^2+A^{}a\right).$$
(55)
This is a completely general result describing the sorts of measurements which can be made on an optical mode by photodetection in the large local oscillator limit. From it can be derived the standard results for homodyne and heterodyne detection (as shown below), and also some more interesting results for adaptive measurements (as I will discuss in Sec. 6). At finite times one has less than maximal information about the system and as a result the effect cannot be factorized in terms of measurement operators. This is a minor generalization of the theory of operations and effects as stated in Sec. 2 ; the same condition (7) for the effects applies. Explicitly, equation (49) shows that the povm for the measurement up to time $`t`$ is

$$F_t(R(t),S(t))=P_0(R(t),S(t))\mathrm{exp}\left(\frac{1}{2}S(t)a_{}^{}{}_{}{}^{2}+R(t)a^{}\right)\mathrm{exp}(a^{}at)\mathrm{exp}\left(\frac{1}{2}S^{}(t)a^2+R^{}(t)a\right),$$
(56)
where $`P_0(R(t),S(t))`$ is defined analogously to (52).

### 5.1 Homodyne Detection

For homodyne detection of a given quadrature of the system, one desires the local oscillator phase $`\mathrm{\Phi }`$ to be constant. This means that $`B`$ is simply given by $`e^{2i\mathrm{\Phi }}_0^{\mathrm{}}e^t𝑑t=e^{2i\mathrm{\Phi }}`$. Thus it can be ignored as a measurement result, and we need consider only

$$A=e^{i\mathrm{\Phi }}_0^{\mathrm{}}e^{t/2}𝑑W(t).$$
(57)
Since $`dW(t)`$ is ostensibly white noise, it follows from equation (50) that we can write $`A=e^{i\mathrm{\Phi }}X`$ where the ostensible distribution for $`X`$ is

$$P_0(X)dX=\frac{1}{\sqrt{2\pi }}\mathrm{exp}(X^2/2)dX.$$
(58)
Thus the effect for such a measurement is

$$F(X)=\frac{1}{\sqrt{2\pi }}\mathrm{exp}(X^2/2)\overline{E}(X)|00|\overline{E}^{}(X),$$
(59)
where $`\overline{E}(X)`$ is an operator defined by

$$\overline{E}(X)=\mathrm{exp}(\frac{1}{2}e^{2i\mathrm{\Phi }}a^{}{}_{}{}^{2}+Xe^{i\mathrm{\Phi }}a^{})$$
(60)
Now using some simple boson operator algebra , it is easy to show that

$$(ae^{i\mathrm{\Phi }}+a^{}e^{i\mathrm{\Phi }})\overline{E}(X)|0=\overline{E}(X)(ae^{i\mathrm{\Phi }}+X)|0=X\overline{E}(X)|0.$$
(61)
That is to say, $`\overline{E}(X)|0`$ is an unnormalized eigenstate of of $`ae^{i\mathrm{\Phi }}+a^{}e^{i\mathrm{\Phi }}`$ with eigenvalue $`X`$. The other factor in equation (59) supplies the normalization, so we can write

$$F(X)=|XX|,$$
(62)
where $`|X`$ is an eigenstate of $`ae^{i\mathrm{\Phi }}+a^{}e^{i\mathrm{\Phi }}`$ satisfying the delta-function normalization $`X|X^{}=\delta (XX^{})`$. This is as has been known for a long time; that homodyne detection can measure a quadrature of the field arbitrarily accurately (see reference  for a “traditional” derivation using quantum Langevin equations). Specific cases of this result have been derived using the theory of quantum trajectories (see for example reference ). However, as far as I know, this is the first time it has been derived as a completely general result using this theory.

### 5.2 Heterodyne Detection

Under heterodyne detection, the local oscillator phase $`\mathrm{\Phi }(t)`$ varies rapidly at a rate $`\mathrm{\Delta }1`$. This means that $`B`$ averages over time to zero. The remaining measurement result $`A`$ is defined by

$$A=_0^{\mathrm{}}e^{i\mathrm{\Phi }(0)+i\mathrm{\Delta }tt/2}𝑑W(t).$$
(63)
It is easy to verify that in the limit $`\mathrm{\Delta }1`$, this has the ostensible distribution

$$P_0(A)d^2A=\frac{1}{\pi }\mathrm{exp}(|A|^2)d^2A$$
(64)
Thus the effect for a completed heterodyne measurement is

$$F(A)=\frac{1}{\pi }\mathrm{exp}(|A|^2)\mathrm{exp}(a^{}A)|00|\mathrm{exp}(aA^{}).$$
(65)
Using the techniques of reference , it is again simple to show

$$\mathrm{exp}(|A|^2/2+a^{}A)=\mathrm{exp}(a^{}AaA^{})\mathrm{exp}(aA^{}),$$
(66)
and thus that

$$\mathrm{exp}(|A|^2/2+a^{}A)|0=\mathrm{exp}(a^{}AaA^{})|0=|A,$$
(67)
where this last state is the normalized coherent state $`a|A=A|A`$. Thus we have

$$F(A)=\frac{1}{\pi }|AA|.$$
(68)
In this case the factor of $`\pi ^1`$ is necessary because the coherent states are overcomplete. This result, that the distribution function for completed heterodyne measurements is the $`Q`$ function $`\pi ^1A|\rho |A`$, has also been known for some time (and again there is a succinct derivation in reference ). It has not previously been derived using quantum trajectories as here.

## 6 Adaptive Measurements

In the two examples considered above, homodyne and heterodyne measurement, the quantum trajectory method for finding the finite time effect simply confirmed already known results. A type of optical measurement for which the effects have never before been calculated, and for which the quantum trajectory method is perhaps the only tractable approach, is continuously adaptive measurements. These are measurements in which the results so far (the photocurrent record $`𝐈_{[0,t)}`$ up to time $`t`$) are used to alter the unitary matrix which determines the measurement operators (28) at time $`t`$. Note that, as expressed in equation (20), this does not change the overall evolution of the system as given by the master equation (24), and so is distinct from feedback onto the system dynamics as explored in Refs. . In the infinite local oscillator amplitude limit being considered here, there is only one parameter which determines (28); that is the local oscillator phase $`\mathrm{\Phi }`$. The feedback-control of this phase is just what is required for a measurement which is of particular interest: that of the optical phase of the system.

A standard measurement of the phase of an optical system consists of making a completed heterodyne measurement yielding the result $`A`$ as above, and then discarding the amplitude information in $`A`$. That is to say, the effect for such a standard phase measurement is

$$F_{\mathrm{std}}(\varphi )=|A|d|A|\frac{1}{\pi }|AA|,$$
(69)
where $`A=|A|e^{i\varphi }`$. Because it discards information, this effect cannot be factorized into measurement operators as in equation (5). In the language of reference , it is an unsharp measurement of phase. By contrast, an ideal measurement of phase  has the effect

$$F_{\mathrm{ideal}}(\varphi )=\frac{1}{2\pi }|\varphi \varphi |,|\varphi =\underset{n=0}{\overset{\mathrm{}}{}}e^{in\varphi }|n,$$
(70)
defined in terms of the unnormalized phase eigenstates $`|\varphi `$. Ideal measurements are the optimal way to extract information which has been encoded in the phase of the field .

Standard measurements of the phase are suboptimal because they measure both the phase and amplitude quadratures of the field, but use only the former information. It is intuitively obvious that one could make a better measurement (i.e. closer to ideal) by concentrating on measuring the phase quadrature. However assuming that one knows the phase quadrature before making the measurement is not in the spirit of a true phase measurement. Rather, one would have to estimate the phase quadrature during the course of the measurement, from the photocurrent so far. This is the essence of the adaptive measurement of phase. Precisely what algorithm one used to estimate the phase quadrature would depend on the knowledge one had about the initial state of the system. Nevertheless, the analysis given above shows that whatever algorithm one used, the estimated phase at time $`t`$, $`\widehat{\phi }_t`$ should depend on the photocurrent $`𝐈_{[0,t)}`$ only through the two complex numbers $`S(t)`$ and $`R(t)`$. In order to measure the phase quadrature, the local oscillator phase $`\mathrm{\Phi }(t)`$ should therefore be set equal to

$$\mathrm{\Phi }(t)=\widehat{\phi }_t[R(t),S(t)]+\pi /2.$$
(71)
When the measurement is over (at $`t=\mathrm{}`$), one’s estimate for the original phase of the system is $`\varphi =\widehat{\phi }_{\mathrm{}}`$.

Exactly how $`\widehat{\phi }_t`$ should be chosen as a function of $`R(t)`$ and $`S(t)`$ is a complicated issue. As stated above, it depends on the knowledge the observer has regarding the preparation of the system. This is evident in figure 1, which shows the Wigner function of the effect $`F_t(R(t),S(t))`$ (56) for some typical values of $`R`$ and $`S`$ which might occur in an adaptive phase measurement. The position and shape of this positive operator in phase space ($`q=a+a^{}`$ and $`p=ia+ia^{}`$) indicates the inherent uncertainties in the results of the single measurement up to time $`t`$ yielding results $`R`$ and $`S`$. As this figure shows, the state $`F_t(R(t),S(t))`$ is a Gaussian state. As $`t\mathrm{}`$ it becomes a minimum uncertainty squeezed state. For heterodyne detection this reduces to a coherent state of complex amplitude $`A`$. For homodyne detection it is an eigenstate of $`p\mathrm{cos}\theta +q\mathrm{sin}\theta `$, where $`\theta `$ is in this case equal to the chosen local oscillator phase $`\mathrm{\Phi }`$. In general $`\theta =\mathrm{arg}(S)/2`$ is a measurement result which changes over the course of the adaptive measurement. The displacements $`x`$ and $`y`$ are defined in terms of the measurement results $`R`$ and $`S`$ by

$$x=\frac{2\mathrm{}(Re^{i\theta })}{1e^t|S|};y=\frac{2\mathrm{}(Re^{i\theta })}{1e^t+|S|}.$$
(72)
The variances in the $`x`$ and $`y`$ quadrature are

$$V_x=\frac{1e^t+|S|}{1e^t|S|};V_y=\frac{1e^t|S|}{1e^t+|S|}.$$
(73)
Evidently, the phase of this state is not unambiguously defined in general. One reasonable estimate would be the phase of the centre of the Wigner function, which would be

$$\widehat{\phi }_t=\theta _t+\mathrm{arctan}(y_t/x_t)=\mathrm{arg}\left[R(1e^t)+SR^{}\right].$$
(74)
However, depending on one’s initial knowledge of the state of the system, a different estimate may be more appropriate.

There is one example for which a different estimate is definitely more appropriate, and for which there is no ambiguity about the estimated phase of the system. This is when there is at most one photon in the system (and the observer knows this). First, I wish to present the results in this case for ideal and standard phase measurements. Projecting the optimal mode effects into the subspace spanned by $`\{|0,|1\}`$ gives the ideal phase measurement effect

$$F_{\mathrm{ideal}}(\varphi )=\frac{1}{2\pi }|\varphi \varphi |,|\varphi =|0+e^{i\varphi }|1,$$
(75)
and the marginal heterodyne effect

$$F_{\mathrm{std}}(\varphi )=\frac{\sqrt{\pi }}{2}F_{\mathrm{ideal}}(\varphi )+\left(1\frac{\sqrt{\pi }}{2}\right)\frac{1}{2\pi }.$$
(76)
That is to say, the standard technique has an efficiency of $`\sqrt{\pi }/288\%`$, in the sense that the same povm would arise from a ideal measurement which worked $`88\%`$ of the time, and which gave a random answer on the interval $`[0,2\pi )`$ the other $`22\%`$ of the time.

I have considered the adaptive phase measurement of a field with at most one photon in detail in reference . The lack of ambiguity in the estimated phase quadrature $`\widehat{\phi }_t`$ is due to the fact that with at most one photon in the field, the measurement result $`S(t)`$ is irrelevant to the estimation of its state, which is obvious from equation (49). The estimated phase quadrature in this case is $`\widehat{\phi }_t=\mathrm{arg}R(t)`$. Substituting this into equation (71) leads to an analytic solution for $`P_0(A)`$:

$$P_0(A)d^2A=\pi ^1\delta (|A|^21)d^2A.$$
(77)
In other words, $`A`$ can contain only phase information, not amplitude information. This is precisely what is desired for a phase measurement, and indeed the povm for this adaptive measurement is

$$F_{\mathrm{adapt}}(\varphi )=_0^{\mathrm{}}|A|d|A|P_0(A)(1+|10|A)|00|(1+|01|A^{})=F_{\mathrm{ideal}}(\varphi ),$$
(78)
Here I have again used the result (54) for an optical mode, projected into the lowest two levels $`\{|0,|1\}`$. Thus for any field with at most one photon, a simple adaptive measurement scheme can produce an optimal measurement of the phase, which is significantly better than the result with no feedback (76).

## 7 Conclusion

In this paper I have investigated the reciprocal links between quantum trajectories and quantum measurement theory. The first link is that one interpretation of quantum trajectories is as a special case of quantum measurement theory, that of continuous observation. Here each measurement has a duration which is infinitesimal on the time scale of the non-trivial system evolution, and is described in terms of operations and effects. This leads naturally to discrete measurement results (photodetections). The equations are actually easy to solve in the limit of a very large local oscillator, in which case the discrete photodetections are replaced by a noisy photocurrent. This is the limit which was considered in detail here.

The second link is that by solving the quantum trajectories over a finite time one can obtain the operations and effects for measurements of finite duration. These results would be difficult to obtain were it not for a formulation of quantum trajectories as linear stochastic Schrödinger equations due to Goetsch and Graham which I have rederived here in a new way. For the case where the the system is not driven, but simply allowed to radiatively damp, the finite measurement time is naturally extended to infinity in order to detect all of the radiated light. Then it is possible to rederive the standard results for homodyne and heterodyne detection, and also to derive new results which would be very difficult to derive by other means. These new results are those pertaining to adaptive measurements, in which, for example, the photocurrent up to time $`t`$ is used to control the local oscillator phase at time $`t`$. Here I presented the results for one case which can be completely solved analytically: adaptive phase measurements on a system with an upper photon number bound of one.

The topic of real-time adaptive measurements begins rather than ends with phase measurements. Feedback-control of measurements can only improve the quality with which properties are measurable compared with standard non-adaptive techniques. For measuring quantities other than phase, one would not wish to work in the infinite local oscillator amplitude limit. Rather, one would have a finite local oscillator, in which case there would be two parameters to be controlled: its amplitude and phase. This is obviously a more difficult problem, and it may be that only a numerical solution using stochastic Schrödinger equations would be possible. Even using numerical techniques, it is not obvious what the general algorithm would be for optimizing such measurements, nor even how optimization should be defined. These are problems for future work, continuing the fruitful interaction between formal quantum measurement theory and its principal realm of application in quantum optics.

## Acknowledgments

I would like to acknowledge valuable discussions with S.M. Tan. This project was supported by the New Zealand Foundation for Research, Science and Technology.


format: text

text:
# Kinetic and Dynamic Delaunay tetrahedralizations in three dimensions

## 1 The Delaunay Triangulation

### 1.1 Conventions

For the sake of clarity, the illustrations in this article will be two-dimensional, unless noted otherwise. Following the notation in the literature  we denote by the term vertex a position111If one extends the algorithms towards weighted triangulations a vertex in addition contains a weight. in three-dimensional space. By an $`n`$-simplex in $`^d`$ ($`nd`$) we understand the convex hull of a set $`T`$ of $`n+1`$ affinely independent vertices, which reduces in the three-dimensional case to tetrahedra (3-simplices), triangles (2-simplices), edges (1-simplices) and vertices (0-simplices). Every $`n`$-simplex has a uniquely defined $`n`$-circumsphere. Recall that a tetrahedron is bound by four triangles, six edges and four points in three dimensions. These $`(nd)`$-simplices $`\sigma _U`$ – formed by the convex hull of a subset $`UT`$ – are also called faces of $`T`$. Since we will work in three dimensions, we will shortly denote 3-simplices by the term simplex. A collection of these simplices $`𝒦`$ is called a simplicial complex if:

* The faces of every simplex in $`𝒦`$ are also in $`𝒦`$
* If $`\sigma _T𝒦`$ and $`\sigma _T^{}𝒦`$, then $`\sigma _T\sigma _T^{}=\sigma _{TT^{}}`$.
(the intersection of two simplices is at most a face of both, the simplices are ’disjoint’)

In numerical calculations with kinetic vertices the above criterion can be destroyed: A vertex might move inside another simplex thus yielding two $`n`$-simplices whose intersection is again an $`n`$-simplex. We will refer to this situation as an invalid triangulation. To be more exact, a triangulation is defined as follows. If $`S`$ is a finite set of points in $`^d`$, then a simplicial complex $`𝒦`$ is called a triangulation of $`S`$ if

* each vertex of $`𝒦`$ is in $`S`$
* the underlying space of $`𝒦`$ is $`conv(S)`$

By the degree of a vertex in a triangulation we will denote the number of simplices in the triangulation containing the vertex as endpoints. Furthermore, we will use the terms tetrahedralization and triangulation in three dimensions synonymously, unless noted otherwise.

### 1.2 Elementary Topological Transformations

To an existing triangulation in $`^3`$ several topological transformations can be applied. We will briefly remind the main ideas. For a more detailed discussion see e.g. . The discussion basically relies on Radon’s theorem (see e.g. ):
Let $`X`$ be a set of $`d+2`$ points in $`^d`$. Then a partition $`X=X_1X_2`$ with $`X_1X_2=\mathrm{}`$ exists such that $`conv(X_1)conv(X_2)\mathrm{}`$.
If $`X`$ is in general position – meaning that every subset of $`X`$ with at most $`d+1`$ elements is affinely independent – then this partition is also unique. In our case this simply means that

* no four points lie on a common plane
* no five points lie on a common sphere

Figure 1 illustrates the idea of Radon’s theorem in three dimensions.

From the Radon partition in $`^3`$ one finds that there exist four possible flips in three dimensions, two for every partition in figure 1. For the case of figure 1a the two possible flips are shown in figure 2.

The flip changing the triangulation from $`1`$ to $`4`$ simplices corresponds to adding a vertex to an existing triangulation. Note however, that in practice the inverse transformation may not always be applicable, since the configuration of one vertex ($`E`$ in figure 2) being the endpoint of exactly four simplices is rarely ever present in a triangulation. This fact – in combination with the existence of non-tetrahedralizable polyhedra in three dimensions – complicates the deletion of vertices from triangulations .

The second partition in figure 1b following from Radon’s theorem requires a more careful evaluation, see figure 3.

The flips $`23`$ (and $`32`$) can only be performed if the polyhedron formed by the five points in $`^3`$ is convex, otherwise the operation would yield overlapping simplices in the triangulation. The convexity of $`A,B,C,D,E`$ in figure 3 can be tested by checking if for every edge $`A,B`$ and $`B,C`$ and $`C,A`$ there exists a hyperplane which has the remaining three points ($`D,E,A/B/C`$) on the same side .

### 1.3 The Delaunay Criterion

Every tetrahedron in $`^3`$ has a uniquely defined circumsphere, if the four vertices do not lie on a common plane (i.e. if the tetrahedron is not flat). The Delaunay triangulation is a triangulation where all the simplices satisfy the Empty-Circumsphere-Criterion: No vertex of the triangulation may lie inside the circumspheres of the triangulation simplices. Thus, the Delaunay triangulation is uniquely defined if the vertices are in general position, i.e. if no five vertices must lie on a common sphere and no four vertices may lie on a common plane .

The simplest method to determine, whether a vertex $`V`$ lies outside or inside the circumsphere of a simplex $`(A,B,C,D)`$ is to solve the associated four sphere equations. However, this problem can be solved more efficiently by adding one more dimension .

Suppose we would like to know whether the vertex $`E`$ lies in- or outside the circumsphere of the simplex $`(A,B,C,D)`$, which we will – without loss of generality – assume to be positively oriented. Then one can proceed as follows (e.g. ): Project the coordinates in $`^3`$ onto a paraboloid in $`^4`$ via

$`A=(A_x,A_y,A_z)A^+=(A_x,A_y,A_z,A_x^2+A_y^2+A_z^2).`$ (1)
The four points $`A^+,B^+,C^+,D^+`$ define a hyperplane in $`^4`$. If $`E`$ is within the circumsphere of $`(A,B,C,D)`$, then $`E^+`$ will be below this hyperplane in $`^4`$ and above otherwise. Consequently, the in-circumsphere-criterion in $`^3`$ reduces to a simple orientation computation in $`^4`$, i.e. by virtue of this lifting transformation one finds

$`\mathrm{in}\mathrm{\_}\mathrm{circumsphere}((A,B,C,D),E)`$ (7)
$`=`$ $`\mathrm{oriented}(A^+,B^+,C^+,D^+,E^+)`$
$`=`$ $`\mathrm{sign}\left|\begin{array}{ccccc}A_x& A_y& A_z& A_x^2+A_y^2+A_z^2& 1\\ B_x& B_y& B_z& B_x^2+B_y^2+B_z^2& 1\\ C_x& C_y& C_z& C_x^2+C_y^2+C_z^2& 1\\ D_x& D_y& D_z& D_x^2+D_y^2+D_z^2& 1\\ E_x& E_y& E_z& E_x^2+E_y^2+E_z^2& 1\end{array}\right|`$
$`=`$ $`\mathrm{sign}\left|\begin{array}{cccc}A_xE_x& A_yE_y& A_zE_z& (A_x^2+A_y^2+A_z^2)(E_x^2+E_y^2+E_z^2)\\ B_xE_x& B_yE_y& B_zE_z& (B_x^2+B_y^2+B_z^2)(E_x^2+E_y^2+E_z^2)\\ C_xE_x& C_yE_y& C_zE_z& (C_x^2+C_y^2+C_z^2)(E_x^2+E_y^2+E_z^2)\\ D_xE_x& D_yE_y& D_zE_z& (D_x^2+D_y^2+D_z^2)(E_x^2+E_y^2+E_z^2)\end{array}\right|,`$ (12)
where a positive sign is to be taken as an affirmative answer222In the general case one will have to multiply by the orientation of $`(A,B,C,D)`$..

Several algorithms have been developed for the construction of static triangulations  as well as for the maintenance of kinetic triangulations , some of which will be discussed in section 2.

### 1.4 The Geometric Dual: Voronoi Tessellation

The most general Voronoi tessellation (sometimes also called Dirichlet tessellation) of a set of generators $`\{c_i\}`$ in $`^d`$ is defined as a partition of space into regions $`V_i`$:

$`V_i`$ $`=`$ $`\{x^n:𝒫(x,c_i)𝒫(x,c_j)ji\},`$ (13)
where $`𝒫(x,c_i)`$ can be an arbitrary function, which reduces in the standard case of the simple Voronoi tessellation to the normal euclidian distance $`𝒫(x,c_i)=\left|xc_i\right|`$. In other words, the Voronoi cell around the generator $`c_i`$ contains all points in $`^d`$ that are closer to $`c_i`$ than to any other generator $`c_j`$. Note that this partition is – unlike the Delaunay triangulation – uniquely defined also for point sets that do not fulfill the general position assumption. Voronoi tessellations have many interesting applications in practice – for a survey see e.g.  – since they do describe influence regions.

In 2 dimensions Voronoi cells are convex polygons completely covering the plane, see e.g. figure 4. This finding generalizes to arbitrary dimensions: The boundaries between two $`d`$-dimensional Voronoi regions $`V_i`$ and $`V_j`$ as defined in (13) reduce to the equation for a $`d1`$ hyperplane. Therefore per definition the Voronoi cells around generators $`Z_i`$ situated on the convex hull of the point set $`Z=\{Z_1,Z_2,\mathrm{},Z_n\}`$ will extend to infinity and thus will have an infinite volume.

Voronoi tessellations can be constructed like the well-known Wigner-Seitz cell in solid state physics , but fortunately there are much more efficient ways to construct the Voronoi tessellation. In this work, we will exploit the geometric duality with the Delaunay triangulation: In any dimension, the corners of the Voronoi polyhedra are the centers of the circumspheres of the $`n`$-simplices contained in the Delaunay triangulation of the Voronoi generators, see figure 4.

The introduction of influence regions also enables the definition of proximity between vertices: We understand two vertices to be direct neighbors (in the sense that their influence regions touch) if they share a common face in their Voronoi diagram or – equivalently – if they are direct neighbors in the dual Delaunay triangulation, see figure 4.

The concept of the Voronoi cell can easily be extended towards generators with a varying strength, the weighted Voronoi tessellation . In such extensions, every generator is assigned a weight, i.e. the functions $`𝒫(x,c_i)`$ in (13) are then given by a function describing the influence strength of the generator $`i`$ at $`c_i`$ on $`x`$. Obviously, the weighted Dirichlet regions can – in contrast to the standard unweighted case – be empty, e.g. if a vertex with a weak influence is surrounded by strong vertices. Among many possible choices for weight functions  we will explicitly mention here the case of power-weighted Voronoi diagrams, also often called the Laguerre complex . It is obtained by assigning a weight $`\omega _i`$ to every generator $`c_i`$, i.e. by using $`𝒫(x,c_i)𝒫(x,c_i,\omega _i)=(c_ix)^2\omega _i^2`$ in (13). In the power-weighted case one can still show that in three dimensions the Laguerre cells are convex polyhedra, whose corners can be obtained from the weighted centers of the corresponding weighted Delaunay triangulation tetrahedra – where the empty circumsphere criterion is simply replaced by its weighted counterpart. Therefore, the Laguerre tessellation or its geometric dual – the weighted Delaunay triangulation – is suitable for collision detection between differently sized spheres. We would like to stress that the algorithms in this article can be generalized to the power-weighted case in a straightforward way.

### 1.5 Voronoi surfaces and volumes

Within the framework of growth models , tissue simulations  and the solution of partial differential equations on irregular grids , not only the neighborship relations in the Delaunay triangulation but also the corresponding Voronoi cell volumes as well as the contact surface between two Voronoi cells may become important. Obviously, to the surface of a Voronoi cell around a vertex $`Z`$ every incident simplex $`\sigma `$ with $`Z\sigma `$ contributes, see figure 4. In two dimensions every triangle $`(Z,N_i,N_{i+1})`$ contributes two surfaces, spanned by the half distances between the vertices $`A=1/2(N_iZ)`$ and $`B=1/2(N_{i+1}Z)`$, and the center of the circumcircle of $`R=COC(Z,N_i,N_{i+1})Z`$. If $`(A,B)`$ are positively oriented (which we will further on assume without loss of generality), then the oriented 2-volume contribution of the simplex $`(Z,N_i,N_{i+1})`$ to the Voronoi cell volume of the vertex $`Z`$ is given by

$`V_Z={\displaystyle \frac{1}{2}}\left(\left|\begin{array}{cc}A_x& R_x\\ A_y& R_y\end{array}\right|+\left|\begin{array}{cc}R_x& B_x\\ R_y& B_y\end{array}\right|\right),`$ (18)
see also figure 5 for illustration. Obviously, the sum of the three volumes has to equal the total simplex volume $`V_Z+V_{N_i}+V_{N_{i+1}}=V_S`$ . One can show algebraically, that this identity holds true for any $`R`$, i.e. also in the interesting case where $`R`$ being the center of the circumcircle lies outside the triangle as in figure 5. In this case one of the two area contributions in (18) will be negative. Now when considering figure 5 it becomes clear that by adding the oriented volume contributions of all simplices containing the vertex $`Z`$ one obtains the correct volume for the Voronoi cell around $`Z`$.

This finding generalizes to three-dimensional volume computation, and the three-dimensional contact surface contribution is in fact nothing but a two-dimensional volume computation, since the contact surfaces between two Voronoi polyhedra are plane polygons.

The algorithm leads to a wrong Voronoi cell volume at the boundary, where the Voronoi cells are per definition infinite. However, if the centers of the circumspheres of the simplices at the boundary do not lie outside the convex hull of the triangulation, then the volume summation yields the part of the Voronoi cells which is inside the convex hull. This is desirable for some configurations, such as the solution of partial differential equations .

The numerical complexity of the volume computation is linear with the number of simplices surrounding the vertex, whereas the complexity of contact surface calculation between two generators grows linear with the number of simplices containing both generators as endpoints. The algorithm can be checked by using that the sum over all Voronoi cell volumes computed this way has to equal the sum of the simplex volumes in the dual Delaunay triangulation.

## 2 Algorithms and Implementation

### 2.1 The Data Structure

As has already been mentioned in the introduction, three-dimensional Delaunay triangulations are much more complicated than in two dimensions. The first difference is that the number of triangulation simplices may vary for a constant number of mobile vertices. Our triangulation basically consists of two data structures:

* a list of the vertices in the triangulation
* a list of connected 3-simplices (tetrahedra) contained in the triangulation

Both structures are organized in a list to enable for the kinetic movement, insertion, and deletion of vertices and the corresponding dynamic update of the simplex list.

A vertex consists of the $`x`$, $`y`$ and $`z`$ coordinates, a weight (which will be assumed to be equal for all points throughout this article) and – to compute the Voronoi cell data efficiently – a vector of all incident tetrahedra. A simplex consists of four pointers on vertices and of four pointers on the neighboring simplices. The latter is required by the fact that we perform a walk in the triangulation (see section 2.3).

The construction of the Delaunay triangulation basically relies on two basic predicates: The determination whether two points lie on the same side of a plane defined by three others and the question whether a point lies in- or outside the circumsphere circumscribing the simplex of four others (enforcement of the Delaunay criterion). By using these two simple predicates the whole triangulation can be built up.

### 2.2 Incremental Insertion Algorithms

Unfortunately in Delaunay triangulations the insertion of one new vertex can change the whole triangulation, but this only holds true for some extreme vertex configurations, for some examples see . In practice, the effect of adding a new vertex to a Delaunay triangulation will nearly always be local. Anyhow, the algorithms we describe here can of course also cope with these worst case scenarios.

So let us assume we have a valid Delaunay triangulation with $`n`$ vertices. Let us furthermore assume that the new vertex lies within the convex hull of the $`n`$ vertices. Then the updated Delaunay triangulation can be constructed as follows (see figure 6):

* Identify all invalid simplices in the triangulation, i.e. all those containing the new vertex within their circumsphere.
* Collect the external facets of the invalid simplices. (Those are the triangles facing valid simplices.)
* Replace the invalid simplices by new ones formed via combining the external facets with the new vertex.

Sometimes this incremental algorithm is also called Bowyer-Watson Algorithm . Once all the invalid simplices have been found, its computational cost is very low (linear with the total number of invalid simplices). At first it actually suffices to find the one simplex which contains the new vertex within its convex hull – the remaining simplices can be found by iteratively checking all neighbors for invalidity, see also figure 6.

The result of this procedure is a Delaunay triangulation with $`n+1`$ vertices.

The algorithm shown in figure 6 is only slightly different than the Green-Sibson Algorithm , which needs the simplex containing the new vertex as an input. Then the elementary topological transformation $`14`$ is performed with this simplex and the resulting (Delaunay-invalid) triangulation is transformed to a Delaunay triangulation by performing $`23`$ and $`32`$ flips until all simplices fulfill the Delaunay property.

Note that for weighted triangulations the simplex containing the new vertex within its convex hull is not necessarily invalid, since the weighted circumsphere does not generally contain the complete simplex. This corresponds to the case of an empty Laguerre cell – the vertex therefore has to be rejected.

The initial triangulation can be given by an artificial large simplex which contains all the data to be triangulated within its convex hull. Therefore, the convex hull of the points to be triangulated is contained within the artificial simplex and is not reproduced by the triangulation. However, in the framework of kinetic proximity structures this has the advantage that one does not have the problem of maintaining the convex hull of moving points, since the artificial simplex does not move. The initial simplex must therefore be large enough to contain the data within its insphere throughout the full time evolution of the simulation. One choice for such an initial simplex is a $`\mathrm{CH}_4`$ configuration.

### 2.3 Location of Simplices

The incremental insertion algorithm requires a first initial simplex containing the new vertex. Many all implementations of Delaunay triangulations perform a walk in the triangulation, for an overview of different walking strategies see e.g. . Note that points can also be located by using the history the triangulation has been constructed (e.g. the so-called Delaunay tree  or history dag ). However, we are aiming at kinetic triangulations, where the length of a history stack could not be controlled.

Here we will use a stochastic visibility walk  to locate a simplex containing a point. Starting with an arbitrary initial simplex $`A`$ and a new vertex $`v`$ to be inserted in the triangulation, in the normal visibility walk one of the four neighbor simplices of $`A`$ is chosen using the following criterion:

* For all four vertices $`a_{i=1,2,3,4}`$ of the simplex $`A`$ check with the new vertex $`v`$:
Are the vertices $`a_i`$ and $`v`$ on different sides of the plane defined by the other three vertices $`a_{ji}`$, i.e. is it visible from the plane defined by the $`a_{ji}`$?
yes : $``$ Jump to the simplex opposite to $`a_i`$.
* If no neighbor simplex is found, the vertex $`v`$ is contained within the simplex $`A`$ and the destination is thus reached.

Obviously, the algorithm can take different pathways (see figure 7) since there may be more than one neighbor fulfilling this criterion. Note also that due to numerical errors the normal visibility walk may loop when triangulating regular lattices (such as cubic, …) that violate the general position assumption. Such situations can be easily avoided by using the stochastic visibility walk, where the order of the vertices to be checked is randomized. Such a stochastic visibility walk terminates with unit probability .

The complexity is directly proportional to the length of the path to be walked – measured in units of traversed simplices. For $`n`$ uniformly distributed vertices for example, the average total number of simplices will grow linearly ($`n`$) with the number of vertices, whereas the average distance between two arbitrarily selected simplices will grow like $`n^{1/3}`$. Once the invalid simplex has been found, the average remaining complexity for the incremental insertion will be in average constant (in $`n`$). Therefore one would expect the overall theoretical complexity to behave like $`\alpha n^{4/3}+\beta n`$ for uniformly distributed points and in higher dimensions $`d`$ as $`\alpha n^{1+1/d}+\beta n`$, see also .

Obviously, the algorithm heavily depends on a good choice of the starting simplex. The method could therefore be improved by checking whether the new vertices lies within a certain subregion which means preprocessing, or it can be speeded up by initially using larger step sizes, e.g. by using several triangulations of subsets of vertices . Alternatively, one can choose the closest vertex out of a random subset of the triangulation to find a good starting simplex . In some sense these two methods are similar, since finding a close vertex is nothing but finding the Voronoi-cell the point is situated in. The second method however does not require the maintenance of an additional triangulation in the case of kinetic vertices.

In many practical simulations, some neighborship relations may already be known when building the initial triangulation. Our implementation of the incremental algorithm therefore expects the vertices to be included in some order, such that successive vertices are also very close to each other in the final triangulation and therefore chooses the starting simplex in the walk algorithm as being the last simplex created if no other guess is given. Note that for processes as cell proliferation, the situation is even better: New cells can be created by cell division which corresponds to the insertion of a new vertex close to an existing one. Consequently, one always has a perfect guess for the starting simplex in these cases.

### 2.4 Updating the triangulation

In three-dimensional kinetic triangulations where the vertices are moving and neighborship relations can change, one has to deal with a changing number of simplices. Once deletion and insertion of points has been implemented, a simple method handling kinetic vertices would be to delete them at their old position and perform an insertion at the new position . However, since these operations involve many simplices, there exist more efficient algorithms.

It is evident that in the case of moving vertices the Delaunay criterion may be violated, i.e. after the vertices have moved one may end up with a triangulation that violates the Delaunay criterion. Even worse, if the vertices move too far, e.g. if one vertex moves inside another simplex, the triangulation will become invalid (contain overlapping simplices). This must be avoided by either computing a maximum step size (see e.g. section 2.5) or by simply keeping the displacements safely small. So let us assume here that after vertex movement one is left with a valid triangulation which possibly violates the Delaunay criterion. Recomputing the whole triangulation is usually not an option for large data sets. The elementary topological transformations in subsection 1.2 however can be exploited to restore the Delaunay criterion. Since we will neither add nor delete vertices in this subsection it is evident that the flips $`14`$ and $`41`$ are not necessary333This will change for weighted triangulations, as vertex movement may lead to trivial vertices that have no associated Laguerre cell volume and vanish from the triangulation.. Consequently, the transformations $`23`$ and $`32`$ will suffice to transform the given triangulation into a Delaunay triangulation, which has been shown to work . With a glance at figure 3 one can see that indeed the flip $`23`$ effectively creates a neighborship connection, whereas the flip $`32`$ destroys it. Therefore we have implemented routines to check either the complete list of simplices or only a small subset for Delaunay invalidity. Our simple data structure enables a convenient calculation of the flip criteria in three dimensions in subsection 1.2. The main advantage of the flip algorithm is that it is – in average – linear in the number of simplices which is also linear with the number of vertices in most practical applications.

We iterate through the list of simplices and check for flipping-possibilities among every simplex (the active simplex) and its neighbors (the passive simplices): Given a simplex $`S`$ and its neighbor $`N_i`$, the flip $`23`$ is performed if the following two conditions are met:

* The opposing vertex of the neighbor $`N_i`$ lies within the circumsphere of $`S`$.
* The five points in the union of the two simplices form a convex polyhedron.

Due to Radon’s theorem it suffices in the last criterion to check whether the edges of the common triangle (without loss of generality $`(S_A,S_B,S_C)`$) are convex with respect to the other two points ($`S_D`$ and $`N_{opp}^i`$) .

The criterion for the flip $`32`$ can be written as follows: Given the simplex $`S`$ and two of its neighbors $`N_i`$ and $`N_j`$, the flip $`32`$ is performed if the following conditions are met:

* The simplex $`N_i`$ is also a neighbor of $`N_j`$.
* The neighbor pairs $`(S,N_i)`$,$`(S,N_j)`$ and $`(N_i,N_j)`$ violate the Delaunay criterion.

If any flip is performed, the new simplices must be inserted at the end of the list of simplices to be checked again. The algorithm terminates as the end of this list is reached.

Note however, that for these flips to be possible, all simplices must be disjoint (the intersection of two simplices may at most be a triangle), i.e. the triangulation must be valid – flips cannot be used to recover from an invalid triangulation. This becomes an issue when computing a maximum step size for the vertex kinetics, compare subsection 2.5.

### 2.5 Deletion of Vertices

In many problems (e.g. mesh generation) the deletion of vertices from a Delaunay triangulation is not of great importance, since there is no great advantage other than a negligible gain in efficiency. However, if the triangulation is used for example for proximity structures or data interpolation, vertex deletion may become important.

Several algorithms have been developed to manage the deletion of vertices in two dimensions, see e.g. . As has already been argued, there exist some fundamental differences between the two-dimensional and the higher-dimensional case. Simply removing a vertex together with its incident simplices leaves a star-shaped hole in the triangulation, which is not necessarily convex. Unlike in two-dimensional case, where a star-shaped polygon always admits a triangulation which can be transformed by flips into the Delaunay triangulation  in three dimensions a star-shaped polyhedron may not admit a tetrahedralization. The simplest example for such a polyhedron is Schönhardt’s polyhedron , reported among others in . The star-shaped holes emerging in Delaunay triangulations however, will always possess a tetrahedralization, which has been proven in . Note however, that this does not generally hold true for constrained Delaunay triangulations .

Another approach for deletion is given in , where the history is used to reconstruct the triangulation such that the vertex has never been inserted. Again, in our approach we did not want to use some kind of history, since for kinetic triangulations there is no way to control the size of the history.

The basic idea of our approach to delete a vertex is to move it towards its nearest neighbor in several steps, each followed by a sequence of flips $`23`$ and $`32`$ restoring the Delaunay property, until the simplices between the two vertices are very flat and can be clipped out of the triangulation without harming its validity. In some sense we project the problem of vertex deletion on the already presented algorithm for vertex movement. Figure 8 illustrates the idea of the algorithm.

The main questions to be answered all reduce to the problem of the step size. How far can a vertex $`v_i`$ be moved into a certain direction without invalidating the triangulation, i.e. without creating overlapping simplices? If the vertex $`v_i`$ penetrates another simplex, the orientation of at least one of its surrounding simplices will change. Therefore one can derive a step size criterion by demanding that the orientation of the simplices incident to $`v_i`$ may not change sign. We define the pseudo-orientation of a simplex $`S_i=(A^{(i)},B^{(i)},C^{(i)},D^{(i)})`$ as follows:

$`𝒱_0^{(i)}`$ $`=`$ $`\left|\begin{array}{cccc}A_x^{(i)}& A_y^{(i)}& A_z^{(i)}& 1\\ B_x^{(i)}& B_y^{(i)}& B_z^{(i)}& 1\\ C_x^{(i)}& C_y^{(i)}& C_z^{(i)}& 1\\ D_x^{(i)}& D_y^{(i)}& D_z^{(i)}& 1\end{array}\right|`$ (23)
$`=`$ $`\left|\begin{array}{ccc}A_x^{(i)}B_x^{(i)}& B_x^{(i)}C_x^{(i)}& B_x^{(i)}D_x^{(i)}\\ A_y^{(i)}B_y^{(i)}& B_y^{(i)}C_y^{(i)}& B_y^{(i)}D_y^{(i)}\\ A_z^{(i)}B_z^{(i)}& B_z^{(i)}C_z^{(i)}& B_z^{(i)}D_z^{(i)}\end{array}\right|.`$ (27)
In the second line we have reordered the terms such that the vertex to be moved is in the first column. In fact, this is – up to a factor of $`1/6`$ – the signed volume of the simplex $`S_i`$. Now suppose that one of the vertices – without loss of generality we have chosen $`A`$ – is moved along the direction of $`\mathrm{\Delta }`$, i.e. $`AA^{}=A+\lambda _i\mathrm{\Delta }`$ with $`\lambda `$ and $`\mathrm{\Delta }=(\mathrm{\Delta }_x,\mathrm{\Delta }_y,\mathrm{\Delta }_z)`$. Then the new pseudo-orientation is obtained via

$`𝒱_1^{(i)}=𝒱_0^{(i)}+\lambda _i\left|\begin{array}{ccc}\mathrm{\Delta }_x& B_x^{(i)}C_x^{(i)}& B_x^{(i)}D_x^{(i)}\\ \mathrm{\Delta }_y& B_y^{(i)}C_y^{(i)}& B_y^{(i)}D_y^{(i)}\\ \mathrm{\Delta }_z& B_z^{(i)}C_z^{(i)}& B_z^{(i)}D_z^{(i)}\end{array}\right|.`$ (31)
If the orientation of the simplex $`S_i=(A_i,B_i,C_i,D_i)`$ is not allowed to change one has found an upper bound for $`\lambda _i`$ via

$`\lambda _i(S_i)={\displaystyle \frac{\left|\left(𝒱_0^{(i)}\right)\right|}{\mathrm{abs}\left|\begin{array}{ccc}\mathrm{\Delta }_x& B_x^{(i)}C_x^{(i)}& B_x^{(i)}D_x^{(i)}\\ \mathrm{\Delta }_y& B_y^{(i)}C_y^{(i)}& B_y^{(i)}D_y^{(i)}\\ \mathrm{\Delta }_z& B_z^{(i)}C_z^{(i)}& B_z^{(i)}D_z^{(i)}\end{array}\right|}}.`$ (35)
Of course this check has to be done for all simplices incident to the moving vertex $`A`$, i.e. with

$`\lambda =\underset{S_k:AS_k}{\mathrm{min}}\lambda _k`$ (36)
one has an overall measure of the maximum step size of $`A`$ in the direction of $`\mathrm{\Delta }`$. If $`\lambda >1`$, then the vertex can simply be moved along the complete path $`(\mathrm{\Delta }_x,\mathrm{\Delta }_y,\mathrm{\Delta }_z)`$, whereas if $`\lambda <1`$ the vertex $`A`$ can only be moved by a fraction $`\alpha \mathrm{\Delta }:\alpha <\lambda `$. Let us furthermore define $`A^{}`$ to be the nearest neighbor of $`A`$. These vertices will have a certain number of simplices in common. For the remaining simplices we define the quantity $`\lambda _{\mathrm{REST}}`$ in analogy to $`\lambda `$ via

$`\lambda _{\mathrm{REST}}=\underset{S_k:AS_k,A^{}S_k}{\mathrm{min}}\lambda _k`$ (37)
Thus, our algorithm for deleting a vertex $`A`$ can be summarized as follows:

1. Find the nearest neighbor vertex $`A^{}`$.
2. repeat

* set $`\mathrm{\Delta }=A^{}A`$
    * determine $`\lambda =\mathrm{min}_{S_k:AS_k}\lambda _k`$
    * determine $`\lambda _{\mathrm{REST}}=\mathrm{min}_{S_k:AS_k,A^{}S_k}\lambda _k`$
    * if $`\lambda _{\mathrm{REST}}1.0`$ move $`AA+\alpha \lambda \mathrm{\Delta }`$ with $`\alpha <1`$ and update the simplices surrounding $`A`$ with flips to restore Delaunay property

until $`\lambda _{\mathrm{REST}}>1.0`$
3. * delete the simplices containing both $`A`$ and $`A^{}`$
    * replace $`A`$ by $`A^{}`$ in all simplices surrounding $`A`$
    * set the correct neighborship relations in these simplices
    * update the simplices incident to $`A^{}`$ with flips

The simplices containing both $`A`$ and $`A^{}`$ will change their orientation in the last step, since their volume vanishes when $`A`$ and $`A^{}`$ merge. However, since these simplices are deleted anyway, their orientation does not need to be maintained within this last step. The orientation of the simplices containing $`A`$ but not $`A^{}`$ (described by $`\lambda _{\mathrm{REST}}`$) however, needs to be maintained, since these simplices will not be deleted afterwards. Therefore, the quantity $`\lambda _{\mathrm{REST}}`$ should be the criterion for the last vertex step, whereas $`\lambda `$ accounts for the maximum length of the previous steps.

Another problem is posed by rounding errors in (35): If the numerator becomes very small – i.e. if one has simplices with an extremely small volume or very skinny simplices, then $`\lambda `$ may tend to assume very small values. Rounding errors are then likely to happen as well. This problem can be weakened by using exact arithmetics  when computing (35) or – when working with random data – by distributing the data over a larger region of space (larger simplices). In our experiences with triangulations of reasonable data such situations actually never occurred.

We have run several tests on the deletion algorithm by first triangulating a number of points and then deleting all the points one by one. Before deleting a point we collected all simplices not containing the point (roughly speaking: the triangulation with a star-shaped hole) and checked for their existence in the resulting triangulation after point deletion. In accordance with the theorem proved in  they were all recovered in the final triangulation – which has also been compared with one constructed by complete re-triangulation.

## 3 Performance

To test our implementation, we performed calculations on a $`1.533`$ GHz AMD Athlon processor with 1 GByte of RAM. The code has been compiled using the GNU g$`++`$ 3.3 compiler with compiler optimization set. The times were then obtained using the clock() command. The seed values of the random number generator have been determined using the system time. In all test runs, the data consisted of 64-bit double variables.

### 3.1 Incremental Insertion Algorithm

The complexities of the walk algorithm and of the Green-Sibson algorithm have been extensively studied . Here we have studied the computation time in dependence of the number of points to be triangulated. Test runs were performed for different configurations of points ranging from $`10^3`$ to $`10^6`$. To avoid the handling of vertex rejection we sticked to the case of equally-weighted vertices (standard Voronoi tessellation).

In a first series of runs, we considered a slightly perturbed cubic lattice with the average lattice constant $`a=1.0`$ (diamonds in figure 9). As a starting simplex for the simplex walk we always took a simplex in the center of the cubus. The expected algorithmic complexity of $`\alpha _cN^{4/3}+\beta _cN`$ is in complete agreement with the simulation. In a second test run, we took the same lattice configuration but gave an imperfect guess for the walk algorithm (squares in fig. 9). This guess was the last simplex created and therefore good within the cubus and bad at the surface of the cubus. For large numbers of points – where the ratio between vertices at the boundary of the cubus and the total number of vertices in the cubus becomes small – we find a linear behavior, as is expected if the cost for simplex location becomes constant. To test for robustness of our code we also fed an unperturbed cubic lattice (data not shown), in this case one finds the worst-case quadratic scaling of $`T_{\mathrm{run}}8.210^8N^2`$. However, in our calculations such extreme triangulations will never occur. For comparison with a uniform random distribution we triangulated different numbers of points within the cubus $`[10.0,+10.0]\times [10.0,+10.0]\times [10.0,+10.0]`$. A much better behavior of the algorithmic complexity is found (spheres in fig. 9). Since for random data nothing is known about the final neighborships, no good guess can be given without some sort of preprocessing. However, by choosing some simplex associated with the closest vertex (out of the last 100 inserted) in analogy to  one still finds a considerable gain in efficiency and a nearly linear scaling in the observed range (triangles in fig. 9). Furthermore, figure 9 shows that the running times for random data with our simple data structure are comparable with the more sophisticated three-dimensional DCFL data structure , and other code , where the used algorithms code scale similarly on random data.

It is evident that the incremental insertion of data points depends on the choice of the initial simplex. Figure 10 shows the increase in the average number of steps necessary for the location of the simplex containing a point using the visibility walk. The expected average $`n^{1/3}`$ relation is found.

### 3.2 Incremental Deletion Algorithm

In simulation of growth models it will often be necessary to delete vertices from a Delaunay triangulation. It turns out, that vertex deletion is more efficient than vertex insertion, since there is no cost associated with simplex location, as the vertices provide $`𝒪\left(1\right)`$ access to the incident simplices. Furthermore, one would expect the average algorithmic complexity of vertex deletion to be constant (i.e. not to depend on the total number of points). In this experiment we have first created a Delaunay triangulation and deleted it afterwards by removing point by point. Again, the mean out of ten test runs has been calculated. Figure 11 gives an impression of the expected linear relationship.

### 3.3 Restoring the Delaunay property

A simulation hosting kinetic vertices will be especially sensitive on the cost of checking all simplices for Delaunay invalidity and restoring the Delaunay property. We have triangulated $`N[10^4,10^6]`$ random points in a vectorizable random lattice  (equivalent to a strongly perturbed cubic lattice) with the additional condition of a minimum distance of 0.1 between all the points. This seemed more realistic to us, since in proximity structures there will always be some minimum distance defined by the object sizes. Afterwards all points are moved by a random small step towards a hypothetical new position $`\stackrel{}{X}\stackrel{}{X}+[RND(m,m),RND(m,m),RND(m,m)]`$ with $`m`$ chosen constant. These hypothetical new positions are put on a list, which is being iterated: Possible steps that do not invalidate the triangulation are performed instantly, the others are divided in smaller sub-steps where the criterion for the maximum step size is the same as that in section 2.5. Then the Delaunay criterion is restored (with $`23`$ or $`32`$ flips). The algorithm terminates when all new positions have been reached.

In practice, this procedure would correspond to one timestep of the application and in the ideal case (where the step sizes are small enough) just one iteration should suffice. However, the additional cost of performing a step size check before actually moving the vertices produces only a factor of roughly 2 in the restoration time and should therefore be preferred to increase the robustness of the algorithm considerably.

As expected, the complexity behaves linear in the number of points, see e.g. figure 12.

By looking at figure 12 one finds that restoring the Delaunay criterion in a slightly perturbed Delaunay triangulation being about 20 times as fast as recomputing the whole triangulation as long as the vertex displacements are small compared to the average vertex distance. Locally updating a triangulation is also faster than using a combination of delete and insert operations , since much less simplices have to be flipped.

### 3.4 Mixed algorithms

To check whether a simulation can cope with a varying number of kinetic vertices, we combined the algorithms on vertex insertion, vertex deletion and vertex movement. For different numbers of uniformly distributed vertices 100 time steps have been performed. In each time step, with probability $`p=0.5`$ an arbitrary vertex was deleted from the triangulation and with probability $`p=0.5`$ a random vertex was inserted. Afterwards all the vertices were moved by a small deviation followed by the restoration of the Delaunay criterion. Therefore, if in average a constant number of vertices are deleted or inserted per timestep, we can expect an overall linear behavior as in table 1.

## 4 Summary

In this article we have shown that it is possible to construct fully dynamic and kinetic three-dimensional Delaunay triangulations by using a very simple data structure. This data structure is obtained by adding neighborship entries to every simplex and by storing the tetrahedra within a list. The performance of our data structure is comparable to that of more sophisticated kinetic data structures , which may pose an advantage for parallelization.

We have proposed a new incremental method of vertex deletion which uses a maximum step size criterion. This criterion also solves the problem of maintaining a valid three-dimensional Delaunay triangulation during vertex movement.

Note that the code allows to be generalized towards power-weighted Delaunay triangulations in a mostly straightforward way by replacing the normal circumsphere criterion by its weighted counterpart. In addition, the code provides functionality to compute volumes and contact surfaces of the associated Voronoi cells which are of importance in some simulations of interacting particle systems. The resulting tessellation of space in Voronoi/Laguerre cells can be used to model growth/shrinking processes or for the numerical solution of differential equations on irregular grids. This implementation of a fully dynamic and kinetic Delaunay triangulation thus makes our code suitable for the simulation of dynamically interacting complex systems with variable particle numbers as e.g. cell tissues.

## 5 Acknowledgments

G. S. is indebted to T. Beyer for discussing many aspects of the algorithms and for testing the code and to U. Brehm for discussing many theoretical problems. G. S. has been financially supported by the SMWK.