Hitting statistics from quantum jumps

We define the hitting time for a model of continuous-time open quantum walks in terms of quantum jumps. Our starting point is a master equation in Lindblad form, which can be taken as the quantum analogue of the rate equation for a classical continuous-time Markov chain. The quantum jump method is well known in the quantum optics community and has also been applied to simulate open quantum walks in discrete time. This method however, is well-suited to continuous-time problems. It is shown here that a continuous-time hitting problem is amenable to analysis via quantum jumps: The hitting time can be defined as the time of the first jump. Using this fact, we derive the distribution of hitting times and explicit exressions for its statistical moments. Simple examples are considered to illustrate the final results. We then show that the hitting statistics obtained via quantum jumps is consistent with a previous definition for a measured walk in discrete time [Phys.~Rev.~A~{\bf 73}, 032341 (2006)] (when generalised to allow for non-unitary evolution and in the limit of small time steps). A caveat of the quantum-jump approach is that it relies on the final state (the state which we want to hit) to share only incoherent edges with other vertices in the graph. We propose a simple remedy to restore the applicability of quantum jumps when this is not the case and show that the hitting-time statistics will again converge to that obtained from the measured discrete walk in appropriate limits.


Introduction
With the advent of quantum information science and the desire to build a quantum computer, the study of quantum algorithms have become an integral part of quantum information theory [1]. Quantum walks have played a special role in quantum computing by providing a platform on which quantum algorithms may be analysed [2,3]. Moreover, they have served as a useful mechanism to describe and explain coherent transport processes in photosynthesis [4,5] and the breakdown of a driven system in an electric field [6]. This has in turn stimulated much experimental effort to realise quantum walks, see e.g. Refs. [7,8,9,10]. Theoretical quantum optics on the other hand has had fruitful applications in the analyses of quantum technologies [11,12,13]. In this paper we will apply the well-known theory of quantum jumps developed in quantum optics to define and calculate the distribution of hitting times in open (i.e. non-unitary) quantum walks.
The paper is organised as follows. We first review quantum jumps and motivate its application to hitting problems in quantum walks in Sec. 1.1. The theory of quantum walks and some existing works related to quantum jumps are then reviewed in Sec. 1.2. We then introduce the necessary background in Sec. 2 with the presentation of the quantum-jump formalism in Sec. 2.1, followed by our approach to continuous-time open quantum walks in Sec. 2.2. In Sec. 3 we explain how the quantum-jump method is to be applied to quantum walks and obtain our first result-the hitting-time distribution. Then in Sec. 4 we derive explicit expressions for the hitting-time statistics. We then relate our result to a previous definition of the hitting time devised for a discrete-time measured walk in Sec. 5. Here the hitting-time statistics of the discrete-time measured walk is shown to converge to the quantumjump approach in the continuous-time limit. A problem with the quantum-jump definition of hitting times is that it becomes inaccurate when there are coherent transitions to the final state. We overcome this problem in Sec. 6,  Figure 1: Typical scenario where quantum-jumps are applied in quantum optics. Usually an atomic system or an optical cavity with a lossy end mirror acts as a source of photons. The output light is measured by a photodetector. It is well known in quantum optics that the photodetection record contains information about the state of the light. The photodetection record can also be used to gain information about the photoemissive source. References [21,22,23,24] are some examples.
extending the quantum-jump approach to arbitrary graphs. We then conclude with a summary of our results and analyses in Sec. 7 and also comment on the relationship of our work with other studies not mentioned in the literature review of Secs. 1.1 and 1.2.

Quantum jumps-from photon counting to hitting times
The quantum-jump approach to dissipative quantum dynamics has traditionally been used for efficiently solving master equations [14,15,16,17], or calculating photon statistics in photon counting [18,19]. It also goes by the name of Monte Carlo wavefunctions or quantum trajectories due to the different contexts in which it was invented (although quantum jumps correspond only to a subset of unravellings within quantum-trajectory theory 1 ). We refer the reader to Ref. [20] for a comprehensive review of how the quantum-jump method was developed. Figure  1 illustrates a typical scenario where quantum jumps are applied e.g. Refs. [21,22,23,24]. There is usually a system (described by a master equation [25]), say a lossy cavity or a two-level atom which dissipates energy into the environment in the form of photons. The emitted photons are then measured by a photodetector. Each count, or photodetector "click" appears as a spike in the photodetection record. The quantum-jump formalism then allows one to compute the system evolution conditioned on a given photodetection record where each detector click is associated with a "jump" in the system Hilbert space. Periods of no clicks then correspond to system evolution with no jumps. The quantum-jump method states that an average over a large number of such conditioned states will reproduce the solution to the master equation. Thus in solving the master equation, one also gains access to the statistics of photon arrival times at the detector. Photon statistics are therefore often calculated using the quantum-jump method [24,26,27]. In fact, Carmichael derived the quantum-jump method [18,19] from the photon-counting distribution due to Kelley and Kleiner [28]. It is this intuition between photodetector clicks and quantum jumps that we will exploit for deriving the hitting-time statistics in quantum walks. To see this we need to first explain what hitting times are: In quantum-walk theory the system, usually called the walker, has a Hilbert space spanned by a countable set of orthonormal states {|ψ n } n , called sites, nodes, or vertices. The hitting time is defined as the time required for the walker to reach (or "hit") some "final" state |ψ N for the first time given that initially it was at |ψ m . The hitting time is a random variable since the first time that the walker arrives at |ψ N will vary from one realisation of the quantum walk to another. It therefore makes sense to speak of the average hitting time and its higher-order statistics. We will often call the hitting-time distribution simply as the hitting distribution for ease of reference. Similarly, hitting statistics (as in the title) refers to the statistics of hitting times. It should be noted that the problem we have defined here has an analogue in classical Markov chains, where the orthonormal set {|ψ n } n is replaced by a set of probabilities. The temporal evolution of the classical walker is then governed by a set of rate equations for the site probabilities [29] which plays the role of the master equation in the quantum case. In fact, one of the quantum-walk models that we will discuss two paragraphs down is based on this idea. We should also mention that hitting times are usually called first-passage times in the classical theory of Markov chains and most areas of science where this concept appears, see e.g. Refs. [30,31,32]. The quantum-jump method lends itself naturally to a hitting problem because the question one is asking, starting from the initial time, is whether the system has made a transition to the final state or not. This leads to a very natural division of the system evolution into two types-either a transition to the final state occurred, or it has PHOTODETECTION RECORD TIME JUMPS TO NO JUMPS Figure 2: An example graph to illustrate the quantum-jump approach to hitting times. We place imaginary detectors that can only detect transitions to |ψN along the edges connected to |ψN . We then lump all the detection records obtained from these photodetectors into one record. A click in the superimposed record therefore corresponds to a transition to |ψN . This can be described by a jump in the Hilbert space of the walker as prescribed by the quantum-jump method. When viewed as such, the hitting time, labelled T , is therefore the time at which the first click occurs.
not, and this corresponds to a system evolution conditioned on "jumps" and "no jumps". We illustrate this idea on an arbitrarily chosen quantum walk in Fig. 2. The walker can move between any two vertices connected by a line, usually called an edge (edges and vertices together form a graph as shown). We are not defining the edges precisely for the moment but one may think of each edge as given by a completely-positive trace-preserving map, i.e. a Kraus channel of some sort [1,33]. Consider the hitting problem defined by the initial and final states shown as |ψ 1 and |ψ N respectively on the graph. Imagine now that each transition on the graph emits a photon. Then we can find out when the quantum walker will arrive at |ψ N by associating photodetectors with the edges connected to |ψ N (assuming our photodetectors only measure photons emitted from transitions to |ψ N ). If we superimpose the photodetection records from these detectors then we will observe an overall record shown on the right of Fig. 2. In this case a photodetector click signifies that our walker is at |ψ N , and the first such click provides the hitting time.
We will not be able to distinguish which channel the quantum walker took to reach |ψ N but this does not matter if we only care about when the walker visits |ψ N . This analogy makes clear how the quantum-jump method can be applied to calculate the distribution of hitting times in quantum walks. The theory of quantum walks is more than two decades old so let us briefly review some of the key developments in this field and then discuss where our work stands in relation to the literature on this topic.

Quantum walks and quantum trajectories
The concept of a quantum walk was first introduced by Aharonov and coworkers in 1993 [34]. Quantum walks can be categorised into one of two classes-either a walk in discrete time, or a walk in continuous time [35]. Quantum walks in discrete time follow the idea originally prosposed by Aharonov and colleagues where a quantum coin is introduced. Classically, a random walk on a line involves a coin toss, the result of which determines whether the walker moves one step to the left, or one step to the right. In the quantum version the coin is simply a two-state system. The movement of the walker is then effected by a (quantum) coin toss implemented as a Hadamard gatê C, followed by applying a shift operatorŜ to change the walker's position conditioned on the coin state [36]. Each application of the unitary operatorÛ ≡ŜĈ then evolves the walker in one time step. For a unitary quantum walk, the notion of hitting becomes fuzzy due to the walker's ability to be in a superposition of sites. As a result, multiple definitions of hitting have been proposed [37,38]. Among these, the definition of hitting based on repeated measurements proposed by Krovi and Brun in Ref. [38] will be useful for our analysis. Applications of discrete-time quantum walks include the design of search algorithms [39,40], and implementing universal quantum computation [41]. Continuous-time quantum walks on the other hand do not use a coin. Instead, the quantum walker evolves according to a unitary operator generated by a Hamiltonian, i.e.Û (t) = exp(−iĤt) withĤ parametrised by hopping rates to adjacent nodes. Quantum walks in continuous time were first introduced by Farhi and Gutmann [42] in close analogy to continuous-time Markov chains [42,43]. Their hitting properties were then studied by Varbanov and colleagues [44] who defined hitting times through jumplike weak measurements. The quantum-jump treatment in the present paper can also be seen as a weak measurement of the walker position in that most of the time one gets a null result (no jumps) which does not modify the walker state very much but once in a while a collapse (a jump) happens which changes the walker state drastically [45]. As with the discrete-time case, continuous-time quantum walks have also been shown to give rise to exponential speedups over classical algorithms for searching [46] and solving black-box problems [47]. It has also been shown that universal quantum computation can be achieved using continuous-time quantum walks [48]. Discrete-time quantum walks were then generalised to allow for non-unitary evolution. These are known as open quantum walks and were first introduced by Attal and collaborators [49,50].
Here the coin changes its state according to a completely-positive trace-preserving map and the evolution operator U for the unitary walk is also replaced by a completely-positive trace-preserving map. A continuous-time model was subsequently proposed by Pellegrini as the continuous-time limit of the discrete-time open quantum walk [51]. Most recently, Liu and Balu proposed yet another continuous-time model of open quantum walks by starting with a master equation in the Lindblad form [52]. This approach to defining quantum walks is in the same spirit as Farhi and Gutmann's for unitary walks as already alluded to above. They (Farhi and Gutmann) viewed the Schrödinger equation as the analogue of the rate equation for probabilities in classical Markov chains. Generalising this to allow for non-unitary dynamics naturally leads to a master equation. We also consider here continuous-time open quantum walks described by a master equation.
Of particular interest for our work is the application of quantum trajectories in Refs. [49,50,51]. These works were motivated by the fact that quantum trajectories offer an efficient means of simulating quantum walks and have the advantage of providing a visualisation of the walk. Lardizabal and Souza have recently defined hitting times using quantum trajectories [53], but for the discrete-time model due to Attal and colleagues [49,50]. Given that quantum trajectories were invented in the context of Lindblad-form master equations the latter is in fact the most natural setting in which to consider them. For such non-unitary graphs one can therefore bypass the formalism of positive operator-valued measures (as in Ref. [44] for example) and give a quantum-trajectory treatment of the hitting problem. Our aim in this paper is to define the hitting time for open quantum walks as the time for the first jump to occur and use this as the basis for calculating the distribution of hitting times and its statistical moments. Our result also corrects Ref. [54] which claims to have calculated the average hitting time for continuous-time open quantum walks. However, an inspection of Ref. [54] shows their model to be purely classical. A summary of the relevant literature is given in Fig. 3.

Quantum jumps
Here we will explain the method of quantum jumps as illustrated in Fig. 1, i.e. in the context of photon counting. This makes the basic structure of the formalism intuitive. Assume for simplicity that the system in Fig. 1 is a single-mode field in an optical cavity with one slightly lossy mirror. The lossy mirror allows the single-mode light to couple to the field outside the cavity which we assume is in a vacuum state. This is described by a master equation which has only one dissipative channel, given bẏ whereĉ = √ γâ with γ being a damping coefficient andâ a bosonic annihilation operator. We are using the notation {Â,B} for the anticommutator of any two operatorsÂ andB. The terms containingĉ modify the usual unitary dynamics (which conserves the system's energy) to include a dissipative process (here being the loss of photons to the vacuum outside the optical cavity). The quantum-jump formalism provides a way to understand these terms as follows. The system state ρ(t) can evolve over an infinitesimal time interval dt conditioned on two types of detection events-a count, and no count (Fig. 1). Given an initial state ρ(t), the probability of registering a count in the interval [t, t + dt) is given by where we have defined This is a natural setting to apply quantum jumps Here we are assuming our detector in Fig. 1 is ideal. The state change conditioned on observing a count (a spike in the photodetection record in Fig. 1) is effected by Equation (4) is called a quantum jump andĉ is referred to as a jump operator. A quantum jump causes ρ(t) to change discontinuously in time.
The second type of evolution is conditioned on the event that no counts were registered in [t, t + dt). This is commonly referred to as the no-jump evolution and occurs with probability The state at time t + dt conditioned on a no-count observation is whereL Note that ℘ 0 (dt) can be written explicitly in terms of (7) as ℘ 0 (dt) = Tr eL dt ρ(t) .
It helps to consider two generalisations of the basic structure given by (2)-(7) for the application of quantum jumps to quantum walks coming up later.

Generalisation 1: Non-unit detection efficiency
The above formulation assumes that every jump in the system is observed with probability one. Therefore the first generalisation that we will consider is to allow for an non-ideal observation, or non-unit detection efficiency where only a fraction η (between zero and one) of the jumps are actually detected [13]. This means that the probablity of detecting a jump in time dt (when it occurs) should be (2) multiplied by η. This can be effected in the quantum-jump approach by letting The probability of detecting a jump in the time interval [t, t + dt) is now The remaining fraction of jumps that went unnoticed then contributes to the no-count evolution and we havē It is not difficult to show that this gives the probability of no counts in [t, t + dt) to be as one would expect.

Generalisation 2: Multiple decay channels
The second generalisation that we shall consider is to include the possibilility of multiple decay channels in the system. Assuming there are L such channels, the master equation in (1) is generalised tȯ If we now assume that each channel has a detection efficiency η k (k = 1, 2, . . . L), the new superoperator effecting a jump conditioned on a detection in [t, t + dt) is given by where we have defined η ≡ (η 1 , η 2 , . . . , η L ). The probability of observing a jump in any one of the L channels in duration dt is The no-count superoperator is now and the corresponding probability for no detections in the interval [t, t + dt) is The generalisation to multiple decay channels allows us to selectively tune the detection efficiency of channel k by choosing a value for η k . In particular we can choose to ignore jumps for selected channels by setting the detection efficiency for those channels to zero. The multichannel theory of quantum jumps with non-unit detection efficiency has been used to show that quantum jumps are more robust to measurement inefficiencies in disproving the existence of objective pure-state dynamical models [55]. The same formalism has also been used to study parameter estimation from multichannel photon counting 2 [24].

Open quantum walks
We now define the problem to which we would like to apply the quantum-jump formalism. Our problem can be set up in analogous fashion to a classical Markov chain problem in continuous time. We are given a quantum system that makes transitions between a set of discrete states The set S N forms an orthonormal basis for the system. The space in which the system resides is thus spanned by S N . In the language of quantum-walk theory, each state in S N is referred to as a vertex and the system is referred to as a quantum walker. We shall also refer to the elements of S N as states, or sites. We must therefore also define the exact manner in which the system makes transitions between different states (analogous to specifying a transtion matrix in classical Markov chains). We assume that our system is undergoing dynamics that can be composed from three basic processes. They are coherent transitions, incoherent population transfer, and dephasing. These processes are represented schematically in Fig. 4 for two vertices and are explained below.

Coherent transition
The ability to be in a superposition of states is what differentiates the random walk of a quantum system from a classical one, so coherent evolution, i.e. any evolution that adds coherences to our system, is indispensable in our theory. Coherent evolution between any two states |ψ j and |ψ k can be described by This is the familiar unitary evolution with a Hamiltonian given bŷ where ω k = ψ k |Ĥ jk |ψ k is the expectation ofĤ jk in the state |ψ k , and We are assuming so thatĤ jk is represented by a real and symmetric matrix. We can define the rate at which coherent evolution occurs by how quickly the transition probability varies in time under (19). The transition probability for a system evolving according to (19) is It can be shown, see e.g. [56], that this evaluates to where An inspection of (24) suggests that we should define the frequency of a coherent transition between any two states |ψ j and |ψ k to be ν jk . The symbol that we shall assign to coherent transitions is shown in Fig. 4(a).

Incoherent population transfer
This is a process which transfers a fraction of the population in one state to another. It is shown in Fig. 4(b) for the case of two states with the direction of the arrow indicating the direction of population transfer. For a system with N states, the incoherent population transfer between any two states, say |ψ m and |ψ n , occurring at a rate k mn is described by Note the order of the indices defines the direction of the process so D mn = D nm for m = n. Equation (26) can be understood by considering an arbitrary state and propagating it over an infinitesimal interval dt. For simplicity let us consider the N = 2 scenario. The result is then most apparent if we use the matrix representation in the basis S 2 . Assuming the process to occur with a rate k 21 and in the direction shown in Fig. 4 where the matrix elements of ρ(t) are abbreviated by From (27) we can see that a fraction of the population in |ψ 1 has been transferred to |ψ 2 and that this amount is given by k 21 ρ 11 (t) dt. We can interpret ρ 11 (t) as the prior probability of finding the system in state |ψ 1 at time t and k 21 dt as the probability that the system will transition to |ψ 2 during dt given that it is in |ψ 1 at time t. Notice that in transferring the population from |ψ 1 to |ψ 2 we also drive the system to a less coherent state. This is seen in (27) as the off-diagonal terms of ρ(t + dt) are less than the off-diagonal terms in ρ(t) by a factor of √ 1 − k 21 dt . This is the sense in which we call the process modelled by (26) incoherent population transfer. The usual rate-equation model of continuous-time Markov chains appears as a special case of incoherent population transfer with ρ 21 (t) = ρ 12 (t) = 0.

Loop (dephasing)
We saw above that incoherent population transfer leads to a reduction of the system coherences. However, a system can also lose coherences without the simultaneous loss of populations. Such processes are called dephasing. In general we can control how classical (or quantum) a state |ψ n is by tuning its ability to share coherences with other states in the graph. Assuming a dephasing rate of q n for |ψ n , the corresponding equation of motion for this process is Again, we illustrate this process for the simplest case of a two-state graph shown in Fig. 4(c). Allowing state |ψ 1 to dephase with a rate of q 1 , we can show that (29) changes only the system coherences by looking at the matrix representation of an arbitrary state in S 2 : It is clear that (30) has off-diagonal elements with the same form as the off-diagonal elements in (27), but now there is no transfer of populations. In fact, (29) is a special case of (26) when n = m, and That is, the population transfer is from state |ψ n back to itself so we expect the population of |ψ n to be unchanged. However, we saw above that incoherent population transfer also reduces the coherences in the graph so dephasing can be viewed simply as a loop shown in Fig. 4(c).

Distribution of hitting times
We are now in position to apply the quantum-jump prescription to calculate the distribution of hitting times. Consider now a graph where each state is connected to another by either one or all of the processes shown in Fig. 4. The question that we would like to answer is how long it takes the system (or a quantum walker) to reach a particular state ρ f ∈ S N for the first time assuming some initial state ρ i . This length of time is called the hitting time. However, each time we run the quantum walk the time taken to reach ρ f for the first time will be different, so the hitting time is actually a random variable. We denote the hitting time by T (ρ f , ρ i ). What we would like to know is actually the distribution of hitting times if we were to run the quantum walk on the same graph many times. Since the hitting time here is a continuous variable, its distribution is governed by a probability density. The hitting probability density is a function We are using the notation Pr{A|B} for the probability of event A occurring given that event B has occurred and an equation like Pr{ρ(t) = ρ } is to be read as the probability of finding the system in the state ρ at time t. The hitting distribution depends on the choice of ρ i and ρ f for a given graph but once chosen they are fixed so ρ i and ρ f are to be thought of as parameters in (32). Without loss of generality we shall define the final state as the N th state in S N as it is just a matter of labelling whereas the initial state will be left unspecified. That is Since ρ i and ρ f are fixed we will not write these out explicitly as arguments of h unless it is helpful to do so.

Derivation without dephasing
We now derive a closed-form expression for the hitting time distribution on an arbitrary graph using quantum-jump formalism. For simplicity we first consider a graph where there are no loops, i.e. no dephasing: Note that we have divided the sum over commutators by half to compensate for double counting. This is becausê H mn =Ĥ nm . In general not every state will be connected to every other state via all possible processes. When this is the case we can simply choose which processes to switch off in (34) by setting the relevant rates zero. One may have also noticed that the sums over the coherent transitions (the commutator terms) do not include transitions to the final state |ψ N . This is because only the incoherent transitions to |ψ N are detectable in the quantum-jump formalism. Later on (in Sec. 6) we will show how the hitting probability density can be obtained when there are coherent transitions to |ψ N by introducing an artifical dissipative channel. In order to apply the quantum-jump technique we first note that (32) can also be rewritten using Bayes's rule as Each factor on the right-hand side in (35) can now be calculated using the prescription laid out in Sec. 2.1. To do this we map (13) directly to (34) by making the following identification In terms of the multichannel theory of quantum jumps in Sec. 2.1.2, each channel is now labelled by an ordered pair (m, n) where m = 1, 2, . . . N , and n = 1, 2, . . . N , and a quantum jump is simply the detection of a transition between any two nodes, say |ψ m and |ψ n . The measurement efficiency for this quantum jump can then be labelled by η mn , and is a number which corresponds to the fraction of jumps observed in the |ψ n → |ψ m incoherent transition. The actual number of jumps from |ψ n to |ψ m will always be greater than the observed number of jumps unless η mn = 1. We can then express the probability of finding the system in the final state |ψ N using the probability of detecting a jump to |ψ N provided that all transitions to |ψ N are monitored with unit efficiency. This identification then dictates how we should set η mn for all m and n, namely that η N n = 1 for any n. Transitions to states other than |ψ N are irrelevant so we simply do not monitor these transitions. Therefore we set η mn = 0 for every m = N . We must use this set of detection efficiencies in the formalism outlined in Sec. 2.1.2 in order for it to be applicable for hitting-time calculations. This prompts us to assign a special label, η , for the detection efficiencies to be used: The jump superoperator (14) on using (36) and (38), becomes An application of J (η ) effects an instantaneous collapse of the system to the state |ψ N . From (15) and (36), the probability of detecting a jump in time dt given ρ(t) is thus If ρ(t) in (40) is a state where no transitions to |ψ N are observed in the time interval [0, t) given some initial state ρ(0), then (40) is exactly the first factor on the right-hand side of (35). Such a state can be directly obtained by solving the no-jump evolution defined by (16), (36), and (37): It helps to rewriteL(η )ρ(t) so that its meaning is more obvious. Expanding the sum containing (1 − δ mN ) we get It is simple to see that and similarly Substituting (43) and (44) into (42) gives where we have definedD As already described at the outset, the first term in (45) (45) is similar in form to D mn except that terms of the formQ N nρ (t)Q † N n have been removed. But these are the terms which describe jumps to the final state. ThereforeD N n is a superoperator which describes how the system evolves conditioned on no jumps to |ψ N . Substituting the formal solution of (41) into (40) we thus obtain Note that (41) produces a trace-decreasing state, which is why we have used an overbar for the density operator and the generator of time evolution. Following Sec. 2.1.2 we see that the norm ofρ(t) is precisely the probability of no jumps in the interval [0, t), given some initial state i.e.
The product of (47) and (48) simply cancels the norm ofρ(t) giving the final expression for the hitting probability density as This completes our derivation of the hitting distribution using the quantum-jump method.

Examples
We now illustrate (49) with a simple example. Consider the graph shown in Fig. 5(a). This is an example with N = 4 so by our convention the final state is ρ f = |ψ 4 ψ 4 |. We will let the initial state be ρ i = |ψ 1 ψ 1 |. In Fig. 5  1 that a Hamiltonian, sayĤ21 is parameterised by the three numbers ω1, ω2, and Ω21, which correspond to the diagonal and off-diagonal matrix elements ofĤ21 in the site basis. Here we set ω1 = 1, ω2 = 3, and ω3 = 5 which are kept constant throughout. We then set the transition rates to be Ω32 = k42 = k43 = 5 and change only Ω21 whose value is shown in the plot. Note that |ω1 − ω2| 2Ω21, so according to (25)  level of coherence. Classically, only incoherent transitions are permitted so hitting distributions do not exhibit oscillations. We can see from Fig. 5(b) that for Ω 21 = 5 the oscillations in the hitting distribution are dying out and the quantum walk is in transition to the classical regime. Furthermore, the oscillations in the hitting distribution can actually reach zero (or near zero) at its minimum meaning that there are times at which it is impossible to find the walker at |ψ 4 . Again, such features would not appear in a purely classical walk since incoherent transitions cannot produce any interference effects. For variation we show in Fig. 6 what happens when we change Ω 32 instead of Ω 21 while keeping every other parameter constant (see the figure caption for the actual values). Like Fig. 5, the amount of oscillations increases as Ω 32 is increased, but now this is accompanied by a lengthening of the tail of the hitting distribution. The shape of the hitting distribution depends on the exact manner in which the probability amplitudes for the different paths of the quantum walk interfere. A comparison of Fig. 5 with Fig. 6 illustrates that hitting times of a given graph can be much more sensitive to one coherent transition but not others. We will have more to say about this in Sec. 4.

Including dephasing
We obtained (49) based on (34) which does not have the dephasing processes defined in (29). The omission of dephasing connections was simply to make the derivation easier to follow. If we include dephasing in L then the dephasing terms will simply carry through toL(η ). The expression for the hitting-time probability density will still be given by (49), just with anL(η ) that includes dephasing connections. Stated more formally, we can consider a graph defined by Note that the total number of incoherent channels in this case is given by L = N 2 and we can make the following identification for the Lindblad operators in (13): The Hamiltonian in (13) is still given by (37). The detection efficiencies should now be set according to This then givesL Of the three types of connections that we consider, only two can induce transitions between states-coherent and incoherent population transfer. The third type, i.e dephasing, changes only the coherences in the system. Of course changing coherences in the graph does affect how quickly the final state is reached. That is to say, the distribution of hitting times for a state diagram with dephasing will be a different function of t compared to one without dephasing but otherwise identical in all respects.
For illustration we include in Fig. 7 hitting distributions for the graph with dephasing of |ψ 2 , see Fig. 7(a). The plots in Fig. 7 are identical to Fig. 5(d) except with various amounts of dephasing added as indicated in the plots. Here the dephasing rate is q 2 and as we increase this rate the quantum walk becomes more classical. This can be seen in passing from Fig. 7(b) to (d) where the oscillations which we attributed to the quantum-mechanical nature of the walk in Fig. 5 gradually die out. This is to be expected since by increasing q 2 we are making the walk more classical, in other words, making it more and more difficult to share coherences between |ψ 2 and other nodes in the graph. Next we calculate the moments of h t;Q N , ρ i .
where we are using the notation E[X] to mean the ensemble average of an arbitrary random variable X. Although having the hitting-time distribution is formally equivalent to knowing the statistics of T (Q N , ρ i ), deriving a closedform expression for (54) is still a nontrivial task. We will show in Sec. 4.2 that (54) evaluates to where we have defined the shorthandL andL + is the Moore-Penrose pseudoinverse ofL , defined bȳ Note the superoperator Hermitian conjugate in (57) is defined with respect to the Hilbert-Schmidt inner product Tr[X †Ŷ ] between any two bounded operatorsX andŶ . We say an arbitrary superoperator A is Hermitian, i.e. A = A † if and only if The Moore-Penrose pseudoinverse always exists, is unique, and recovers the standard inverse whenL is invertible. One may wonder howL + can actually be computed for a given graph. To find the pseudoinverse it will be most convenient to use a matrix representation forL . This can be accomplished by using {Q mn } m,n as an operator basis, in which case the Hilbert-Schmidt inner product allows us to represent an arbitrary state ρ by a vector ρ whose elements are ρ mn . Note in this representation each row of ρ is labelled by two indices so ρ is a N 2 × 1 vector. SimilarlyL can be represented by an N 2 × N 2 matrix (L ) whose elements are The superoperator pseudoinverseL + will then be faithfully represented by the Moore-Penrose pseudoinverse of (L ), i.e. (L + ) = (L ) + . Equation (59) will be useful in Sec. 4.3 when we consider simple examples. There, we will consider one example whereL −1 does not exist and another where it does. This will explicitly illustrate the use of the Moore-Penrose pseudoinverse in (55).
The two most considered hitting statistics are the average and variance. We therefore pay special attention to these two quantities. On setting n = 1 in (55) we obtain the average hitting time: Setting n = 2 in (55) and using the result for the average we arrive at an expression for the variance Our expression for the nth moment makes the evaluation of the hitting statistics a simple and efficient procedure givenL . For example, evaluating the average hitting time by setting n = 1 in (54) and evaluating the resulting integral requires the determination of an optimal upper limit to truncate the integral in order to balance the amount of simulation time and numerical precision. In particular, when h(t) has a long tail, the calculation of the average hitting time via (54) can be cumbersome. Equation (60) has the advantage that it avoids these issues.

Proof of the nth moment
To derive (55) we consider the moment-generating function M T (x) defined by Given M T (x), all higher-order moments of T can be calculated as A formula for M T (x) is simple to derive by substituting the expression for h(t;Q N , ρ i ) from (49) into (62). We obtain Here we are using 1 for the superoperator identity. We now need to evaluate the t −→ ∞ limit in (66). To do so we first note that This simply says that the steady state for the dynamics defined by the generatorL is the zero state (represented by a matrix with all elements equal to zero). We can prove this by showing that (45) produces a trace-decreasing density operator. It is simple to check from (26) that D mn is a traceless superoperator: It is also simple to see that the trace of a commutator is always zero. Therefore taking the trace of (45) we get The traces on the right-hand side are positive since they are averages of positive operators in the stateρ(t). They are in fact just the diagonal elements of the unnormalised density operatorρ(t). Since the rates k N m must also be positive for every value of m, Tr[ρ(t)] has a negative rate of change, which means thatρ(t) will eventually decay to zero. We now assume the evolution under exp(L t) dominates over exp(xt) [i.e. the approach to zero under exp(L t) is faster than the rise of exp(xt)] then the t −→ ∞ limit can be dropped to give It is easy to see that Using (63) we therefore arrive at We assumed thatL has an inverse in our proof and it is natural to ask what ifL −1 does not exist? If the inverse ofL does not exist then one can replaceL −1 byL + , the Moore-Penrose pseudoinverse. As mentioned in Sec. 4.1, L + always exists, is unique, and coincides withL −1 if it can be found. The replacement ofL −1 byL + in (73) then gives (55). Equations (60) and (61) are obtained from (55) which we have derived using the moment-generating function but they should also be derivable by directly doing the integral in (54) with n = 1 and n = 2 respectively. We show that this is indeed possible in Appendices A and B, which serve as independent proofs for the mean and variance of the hitting time. Let us now illustrate how these results work with a few simple examples below. We consider only the mean and variance of T .

Example 1: An absorbing state for N = 2
For our first example we consider the case of a simple incoherent transition between two states shown in Fig. 8 with k 12 = 0. It is clear that |ψ 2 is an absorbing state and we will see that our hitting statistics says so as well. The graph is defined by where we have usedQ † 21Q 21 =Q 1 . The average and variance of T (Q 2 , ρ i ) according to (60) and (61) are thus This then givesL asL Using (59) the matrix representingL reads: This is clearly not an invertible matrix since the last row has only zeros. However, it has a Moore-Penrose pseudoinverse given by the same matrix with the first three diagonal entries inverted. In operator form we havē Repeated application of (79) then gives Taking the average of (80) and (81) against |ψ 1 we get where we have defined ρ mn = ψ m |ρ i |ψ n . The average and variance of T (Q 2 , ρ i ) according to (75) and (76) are thus PHOTODETECTION RECORD TIME Figure 8: A simple example with N = 2. We are concerned with the hitting time for state |ψ2 so we are only detecting the |ψ1 −→ |ψ2 transition.
Note the mean and variance are independent of any initial coherences in the graph so the random walk is classical. This can be expected because incoherent transitions do not couple the site populations to its coherences [see (27) of Sec. 2.2.2]. For ρ i = |ψ 1 ψ 1 | we have ρ 11 = 1 and This says the average amount of time we would have to wait to see the walker reach |ψ 2 starting at |ψ 1 is the inverse of the incoherent transition rate and this has a spread equal to the square of the mean as expected. If on the other hand ρ i = |ψ 2 ψ 2 |, then ρ 11 = 0, and we obtain When the initial and final states coincide, i.e. ρ f = ρ i , the quantity T (ρ i , ρ i ) is known as the recurrence time of ρ i [29]. It is the time taken by a walker starting at ρ i to return to it. If ρ i has an infinite recurrence time then it is referred to as a null-recurrent state, and positive recurrent if its recurrence time is finite. An absorbing state ρ i can thus be seen as a special case of a positive-recurrent state with zero recurrence time since the walker never leaves ρ i . This is precisely what (86) says. The variance should also be zero as there should be no spread in the recurrence time of an absorbing state. Although simple, this example illustrates the fact that we actually just require the inverse of the nonzero block inL that describes the transitions between |ψ 1 and |ψ 2 . The reason why the last row (and column) inL is zero is simply because there are no transitions out of |ψ 2 .

Example 2: A recurrent state for N = 2
We contrast the example above with the graph shown in Fig. 8 where another incoherent transition from |ψ 2 to |ψ 1 has been added. In this case we modify (77) tō The matrix representation ofL can be shown to be It is clear that (L ) has linearly independent columns. ThereforeL has a standard inverse which in operator form is given byL As before, repeated application of (89) gives Substituting these into (75) and (76) and simplifying we get (93) Looking at Fig. 8, we can see that if ρ i = |ψ 1 ψ 1 | then we would expect the average hitting time to be identical to the case without the transition from |ψ 2 to |ψ 1 . This is indeed what we find on setting ρ 11 = 1 and ρ 22 = 0 in (92) and (93). If instead we now set ρ i = |ψ 2 ψ 2 | so that ρ 11 = 0 and ρ 22 = 1, we obtain the mean and variance for the recurrence time of |ψ 2 , (94) Intuitively we would say the round trip time T (Q 2 ,Q 2 ) should be the sum of the hitting times for each direction: Taking the average of (95) then reproduces the average in (94). Physically we also expect T (Q 1 ,Q 2 ) to be independent of T (Q 2 ,Q 1 ), so taking the variance of (95) we get V[T (Q 1 ,Q 2 )] + V[T (Q 2 ,Q 1 )] which is precisely the variance seen in (94).

Example 3: Effects of coherence for N = 4
So far we have only illustrated our formula for the average hitting time for the simplest of examples-graphs with only incoherent transitions and with the minimum number of nodes. Although they are oversimplified, these examples allow us to get a handle on (60) by using our intuition from classical random walks. Now that we are somewhat more comfortable with (60) we consider a more nontrivial example. We return to the graph discussed in Fig. 5(a) and examine how different amounts of coherence in the graph can change the average of T (Q 4 ,Q 1 ). We have deliberately chosen the graph parameters so they match those used in Fig. 5. The parameters which are held constant have their values shown in the caption of Fig. 9. In Fig. 9(a) we plot the average hitting time as a function of Ω 21 for three different values of Ω 32 [corresponding to the three values chosen in Fig. 5(b), (c), and (d)]. It can be seen from Fig. 9(a) that the average hitting times all behave, qualitatively the same as Ω 21 is varied for the three values of Ω 32 considered: They all start at infinity at Ω 21 = 0, decay to some minimum, and then increase again to some constant value. From Fig. 5(a) we see that if Ω 21 = 0 then |ψ 1 becomes an absorbing state and the quantum walker stays at its initial position forever. This is why E[T (Q 4 ,Q 1 )] −→ ∞ as Ω 21 −→ 0. Thus as we increase Ω 21 from zero, we allow the quantum walker to access the rest of the graph and the average hitting time must drop from infinity. In particular it reaches the minimum when there is non-negligible population in both states |ψ 2 and |ψ 3 . Increasing Ω 21 further makes the |ψ 1 ←→ |ψ 2 transition dominant which leads to a small probability of finding the walker in |ψ 3 . Hence when Ω 21 Ω 32 , essentially all first transitions to the final state are from |ψ 2 and the hitting time becomes independent of the actual values of Ω 21 and Ω 32 .
A much less trivial feature is present in Fig. 9(b). This time we fix the value of Ω 21 and plot the hitting time as a function of Ω 32 . Again, all the curves are qualitatively the same and show that Ω 32 Ω 21 also increases the hitting time and can even make it diverge. We emphasise that diverging hitting times are a purely quantum phenomenon as a classical walker always reaches the final state in a finite time (assuming of course that there are no absorbing states). We can shed some light on this quantum effect by solving the eigenvalue problem of the Hamiltonian driving the |ψ 1 ←→ |ψ 2 and |ψ 2 ←→ |ψ 3 transitions. Since ω 1 , ω 2 , and ω 3 are of the same order, let us assume they are all equal to zero for simplicity. In this case one can easily show that for Ω 32 Ω 21 the state |ψ 1 becomes the eigenstate of the Hamiltonian and hence it is a stationary state. In other words, by increasing the rate of coherent transitions |ψ 2 ←→ |ψ 3 we are trapping the population in the initial state |ψ 1 , and thereby increasing the hitting time. This new mechanism leading to diverging hitting times should be contrasted with two other mechanisms known previously in the context of measured quantum walks [44]: (i) When the measurement frequency tends to zero (since the walk is essentially not measured at all the hitting time tends to infinity); (ii) when the measurement frequency tends to infinity (quantum Zeno effect freezes the population in the initial state).

Generalised Krovi-Brun definition
There are a few noteworthy differences between what we have considered thus far and what is usually considered in the standard theory of quantum walks: First, much of quantum-walk theory considers only closed systems i.e. systems undergoing unitary evolution whereas the theory presented here allows for non-unitary evolution (although open quantum walks are starting to attract some attention recently as already mentioned in Sec. 1). Second, time is often discrete. In this case the hitting time is simply measured in the number of steps that the system takes to reach the final state. The hitting probability is then defined analogously to (32) to be Lastly, as we said already for quantum walks with only unitary dynamics, the notion of a quantum walker actually being in a particular state is not well defined. One way around this issue is to introduce a measurement of the finalstate population at the end of every step. Given that our hitting distribution is founded on measurement theory, it is most natural to compare our calculation with this version of the hitting distribution from the quantum-walk literature. In particular, we use the definition stated in Krovi and Brun's paper [38] but replace, in their definition of the hitting distribution, the unitary time-evolution operator by a non-unitary map K(t). Taking ρ f =Q N as before, the expression for f (n; ρ f , ρ i ) from Ref. [38] gives where the quantum-walk is defined by K(δt) = exp(L δt). The superoperators Q N and P N are defined by The effect of Q N is thus to project the system into |ψ N while P N projects the system into the subspacē Their effects on an arbitrary state ρ can most easily be seen when (98) are in their matrix representations where we have made use of (28) for the matrix elements of ρ (sometimes with a comma in the subscript when the indices of ρ are otherwise difficult to differentiate). Note that we have written the time step in (97) as a small but finite number δt to remind ourselves that (97) was derived in discrete time. Therefore the comparison of h(t; ρ f , ρ i ) to f (n; ρ f , ρ i ) has to be made in the limit of δt −→ 0. As a matter of fact, we will show next that the two distributions converge in the continuous-time limit, expressed as

Convergence of the generalised Krovi-Brun and quantum-jump distributions
To show the equivalence between the two hitting distributions in the δt −→ 0 limit it helps to separate the full generator for the Markov chain into jump and no-jump terms: whereL is given by (53) and we have also defined the shorthand Recall that we defined J (η ) in (39). We will assume that initially there is no population in the final state, This condition then implies that |ψ N cannot share coherences with the remaining states inS N . The initial state is thus confined to be inH N , which can be formally expressed as The idea of writing L as the sum ofL and J is that the no-jump and jump superoperators have the following useful properties which can be used to show the consistency between the discrete-time and continuous-time hitting distributions. The first thing to note is that the no-jump evolution of an ρ i defined by (106) is confined toH N . This automatically means that any future evolution of ρ i underL has zero projection onto |ψ N . We can express these observations formally as Remember that these are superoperator identities assuming they act on states satisfying (106). Jumps, on the other hand, take the system state out of its confinement inH N and into |ψ N . This means that We then have, in the limit of δt −→ 0, On rewriting the expansion in δt back in exponential form we arrive at This then permits us to writeρ This is again a state which satisfies (106) [with ρ i replaced byρ(t n−1 )] so the superoperator identities (107) and (108) still apply and we have where we have also noted that Q Nρ (t n−1 ) = 0. Taking the trace then gives which we recognise as the jump probability (40) with an unnormalised state, but we showed in (47) and (48) that this is precisely the hitting distribution. The only difference is that time is discretised. Substituting the definition of J from (39) into (117) we get Note that f (n) corresponds to the quantum-jump distribution evaluated at time t n−1 . This makes sense since there are only n time steps so what we want is the probability that a jump occurs in the interval from t n−1 to t n and that none occurred from t 0 to t n−1 -which is precisely h(t n−1 ) δt. 6 Including coherent transitions to the final state 6.1 The N + 1 model The theory above does not allow for coherent transitions to |ψ N in the graph. The quantum-jump approach works by detecting transitions to the final state but only incoherent transitions can be detected. That is to say, if there are coherent transitions fromH N to |ψ N , then the hitting statistics calculated from the quantum-jump approach will not reflect the true statistics, which contains a contribution from the coherent transitions to the final state.
On the other hand the Krovi-Brun formula (97) works by detecting the on-site populations and therefore applies to arbitrarily complicated graphs. We will now remedy this problem so that the quantum-jump approach can be applied to arbitrarily complicated graphs as well.
Let us introduce one extra state, a fictitious state, to the quantum walk. We will denote this extra state as |ψ N +1 . By introducing an incoherent transition from |ψ N to |ψ N +1 we can expect to recover the original hitting time (i.e. the time to reach |ψ N for the first time) by making the transition rate k N +1,N large in the extended graph.
Since the rate k N +1,N is an important quantity in this section, and is also a bit cumbersome to write repeatedly, we shall set v ≡ k N +1,N . (120) It is then intuitive to see that since in the limit of large v, the time taken to go from |ψ N to |ψ N +1 becomes negligible. In fact, in the limit of large v, we can also expect T (Q N +1 , ρ i ) and T (Q N , ρ i ) to be identically distributed. This is because as we make v larger and larger, we are also making the probability of reaching |ψ N +1 from |ψ N tend to one. When this probability is one, the quantum walker is as likely to be in |ψ N +1 as |ψ N , or more formally, Note the hitting distribution to reach |ψ N is written using the Krovi-Brun formula since the existence of coherent transitions to |ψ N limits the use of the quantum-jump method. We refer to the Krovi-Brun formula f (n;Q N , ρ i ) simply as the N model and the quantum-jump formula h(t;Q N +1 , ρ i ) dt in the limit of large v as the N + 1 model. We prove (122) below.

Convergence of the N and N + 1 models
By adding one extra state and an incoherent connection from |ψ N to |ψ N +1 we may decompose the dynamics of the N + 1 graph into jump and no-jump terms as we have been doing, Note that (123) now acts on density operators in an N + 1-dimensional space whereL describes the evolution of the quantum walk constrained to the original N -dimensional subspace. The superoperator J now describes a jump from |ψ N to |ψ N +1 , defined by In order to make the connection with the N model, it helps to introduce The proof of (122) can be made easier if we discretise time into steps of size δt and express h(t) dt in the limit of δt −→ 0. In this case an arbitrary t is equivalent to t n ≡ n δt for an arbitrary n. The hitting distribution calculated using quantum jumps is thus where P N +1 is defined in the same way as (98), (99), and (101) with N −→ N + 1 [recall also that we have defined Q n = |ψ n ψ n |]. The last line simply makes use of the proof in Sec. 5 to express the conditional evolution in terms of P N +1 and L N +1 . As we will be taking the limit of large v it makes sense to plug in the definition of J to make v explicit in the hitting distribution: There are in fact two places where v appears in (129): One in the product v δt, and the other implicitly in exp(L N +1 δt). The first limit is simple to see. All we have to do is note that (129) is already in the δt −→ 0 limit and that v −→ ∞ must be consistent with the interpretation of v δt as the probability to jump from |ψ N to |ψ N +1 in δt. This means that v δt −→ 1 as we make v ever so large but δt ever so small. This immediately gives The only nontrivial part now resides in what happens to the evolution of the N + 1 graph for large v: The larger we make v, the more difficult it is for the population of |ψ N to build up. Then, as v tends to infinity there is simply no appreciable population in |ψ N for all time. Since ψ N |ρ(t)|ψ N = 0 implies no coherences can be shared between |ψ N and any other state, we can write, in the v −→ ∞ limit, This simply says that the N th row and N th column of the density operator at all times is zero when v −→ ∞. Now we note that the only sensible initial state is one satisfying since the original hitting problem is defined only for the N -state graph. From this initial condition it is easy to see that The properties in (133) are simple to see since the action of P N +1 just says that we find the system to be in a state in the subspace spanned by S N . The first condition in (133) then follows since L N only evolves the system within the subspace spanned by S N . The second condition in (133) is true because J puts the quantum walker in the state |ψ N +1 so that one finds the population of every other state to be zero. These are the same conditions as those appearing in (107) and (108), just extended to N + 1 states. We can use (133) to simplify the limit in (130) by noting that We therefore have, upon using (133) and (134), Since P N +1 can be commuted through P N and L N , it is simple to see that Substituting this back into (130) then gives, for v −→ ∞ and δt −→ 0, As this equation already assumes δt −→ 0, this is precisely the right-hand side of (122). This therefore shows the equivalence between the N and N + 1 models. We now illustrate (122) for the graph in Fig. 10(a) where we are interested in the hitting time for vertex |ψ 4 starting at |ψ 1 . All the graph parameters except for v are fixed and have values stated in the figure caption. In this case we consider the hitting-time distribution for |ψ 5 under the quantum-jump approach and compare it to the generalised Krovi-Brun distribution taking the final state to be |ψ 4 . The hitting distribution calculated from the quantum-jump approach is shown as the (blue) solid line while the generalised Krovi-Brun approach is shown as (red) dots. The value of v is gradually increased as we pass from Fig. 10(b) to (d), and as can be seen, the hitting distribution from the quantum-jump approach (solid curve) changes until it eventually coincides exactly with the generalised Krovi-Brun hitting distribution (dots). The parameter values are Ω43 = Ω32 = k42 = k43 = 5, Ω21 = 50, ω1 = 1, ω2 = 3, and ω3 = 5. The convergence between the N + 1 and N models is clearly seen as v is increased from v = 5 in (b) to v = 500 in (d).

Discussion
The most essential result of this paper is the explicit expression for the statistics of hitting times for a model of continuous-time open quantum walks given in (55). This is a finite sum and in principle can be entered into a computer to extract meaningful quantities about the quantum walk such as the average hitting time and its variance. Below we summarise our results in conjunction with other findings obtained in this paper. We also explain the relationship between our work and some other literature which we have not mentioned so far but are nevertheless worthy of a discussion.

Summary
We have shown here how the quantum-jump method lends itself naturally to the hitting problem in continuous open quantum walks. In this language the hitting time for the walker to reach a prescribed state ρ f from some fixed ρ i in an N -state graph is defined by the time of the first jump. Based on this, the hitting distribution is just the distribution of times for the first jump, given in (49). Simple examples of this result were considered. We then derived expressions for the statistical moments of the hitting time in terms of the graph dynamics. The results are given by (60), (61), and (55). Again, some simple examples were used to illustrate the final results. In the process we showed analytically how the quantum expression simplifies to the expected classical result when the graph has only incoherent transitions. We then generalised the hitting distribution for the discrete-time unitary quantum walk in Ref. [38] to allow for non-unitary evolution and showed that this is equivalent to the quantum-jump distribution if the time step of the discrete walk approached zero. A caveat of the quantum-jump method is that it assumes the edges connected to ρ f are incoherent. The hitting time becomes much less straightforward to define when coherent transitions to ρ f are present. Thus open quantum walks have an intrinsic advantage over unitary walks in that there is a set of graphs for which hitting times can be unambiguously defined (i.e. graphs where ρ f only shares incoherent edges with other nodes). We then proposed a solution to this problem by adding one fictitious state to the graph giving us a model with N + 1 states. We then proved that our N + 1 model predicts the same hitting statistics as the N model according to the generalised definition of Ref. [38]. This is illustrated for a simple graph where the convergence of the N + 1 model to the N model can be seen "explicitly" (Fig. 10). It is interesting to point out that the N and N + 1 models differ qualitatively in that the former decides if the final state has been hit or not by detecting the walker's presence at the final state, whereas the N + 1 model works by detecting if the walker is in transit to the final state.

Relation to other work and possible future explorations
Our appeal to photon counting as a way of thinking about quantum jumps for the hitting problem raises another interesting question related to the use of measurements in quantum walks. This is related to our allusion to other unravellings made earlier in Sec. 1.1. In general, an unravelling is a specific decomposition of the master equation into stochastic trajectories such that the ensemble average of them recovers the dynamics of the master equation. Two classes of unravellings have special significance in quantum optics-the jump unravellings used in our paper, and the so-called diffusive unravellings-because they turn out to have concrete interpretations in terms of quantum optical measurements (as already seen with the jump unravellings in our work) [13]. Canonical examples of quantum-optical measurements that give rise to (or essentially realise) diffusive unravellings are the homodyne and heterodyne detection schemes [57,58]. Since we are defining the graph of a continuous-time open quantum walk by a master equation, and we know that a master equation can be unravelled in different ways, one might wonder if another type of unravelling can be used to study a hitting problem. Given that diffusive unravellings form the second major class of unravellings in quantum optics we might then ask if they can be used for the hitting problem in this paper. The short answer is that there is not much motivation to do so because they are unhelpful for solving hitting problems with discrete-state quantum walks. To understand this we recall that in our quantum-jump approach we imagined the incoherent transitions in the graph gave off photons which were measured by photodetectors. Using a diffusive unravelling is equivalent to superimposing the stream of photons coming from an incoherent transition with a local-oscillator field (a laser in a coherent state) on a beam splitter before it is detected by the photodetector. In this case the photodetector actually measures a quadrature of the stream of photons determined by the phase of the local oscillator. A diffusive unravelling thus gives us information about a continuous variable that does not directly reveal to us if the walker has reached a prescribed state or not. To get around this one can then try to introduce the fidelity between ρ f (here taken to be a parameter) and the walker's state to define when ρ f is reached or not. As one can see, this already introduces an extra step into the analysis, not to mention that the fidelity for a fixed ρ f is a nonlinear function of the density operator. The jump unravelling on the other hand does give us direct information about whether the walker has reached ρ f , and is much more natural for our quantum-walk problem. If one wants to think about applying diffusive unravellings one could start with a completely different quantum-walk model altogether where diffusive unravellings would appear to be natural. A suitable quantum-walk model would be one in which the state space is continuous. To this end, it might be interesting to use quantum Brownian motion as a model. A diffusive unravelling of this can then be expected to reveal the walker's position directly.
There is in fact an example of a hitting-time calculation for a diffusively unravelled quantum master equation in the context of rapid purification [59,60,61]. The basic idea is to purify qubits as quickly as possible by using a continuous diffusive measurement and feedback. There are two approaches to rapid purification-one can either maximize the average purity of a qubit in a fixed time, or minimize the average time for a qubit to reach a fixed purity [62]. It is the latter approach that is a hitting problem where the hitting time is simply the time to reach a predefined purity. In this case, the average time to reach a predefined purity (i.e. the average hitting time) is calculated by first converting the qubit dynamics under a diffusive measurement into a set of stochastic Bloch equations, and then from these, the equivalent classical Fokker-Planck equation. Once the classical Fokker-Planck equation is known, standard techniques for calculating the average time to reach a given purity can be applied [63] (where the hitting time is known as the first-passage time). Aside from the difference in the unravelling that is employed, there is another difference between the calculation of the average first-passage time in the rapid purification work and our results here: Our average hitting-time calculation (and in fact all higher statistics) uses a quantum method-quantum jumps-and hence our results are expressed explicitly in terms of the quantum evolution equation of the random process. The rapid-purification work on the other hand re-expresses the evolution of the random process in terms of the equivalent classical equation first (the Fokker-Planck equation), from which  the average hitting/first-passage time can then be obtained using a method invented for classical diffusive processes. In light of this, a noteworthy point of our paper is that the nature of our problem, and the ensuing method used, allows us to bypass the extra step of making a quantum-classical correspondence as in rapid purification before the hitting problem can be solved.
We have seen that in the language of photon counting, the hitting time is the analogue of the time of the first count. In fact, our calculation of the hitting-time distribution follows much the same strategy as what a quantum optician would do to derive the distribution of waiting times (the time interval between successive photon counts). Is the hitting-time distribution in our quantum-walk problem therefore equivalent to the waiting-time distribution of photon counting? The answer is no, and to see this we just have to refer back to Fig. 2. If we look at the record of jumps shown on the right in Fig. 2, the waiting time is the interval between successive jumps. However, from the graph on the left we can see that this is the amount of time the walker spends being away from the final state before returning to it again. The time taken for a random walker to start in a given node and return to it turns out to be an interesting quantity as well. It is known as the recurrence time (recall Examples 1 and 2 in Sec. 4.3). The recurrence time in a hitting problem is thus the quantity that one should regard as analogous to the waiting time of quantum optics and it is not difficult to see that it is different to the hitting time. It is in fact simple to construct a graph for which the hitting time is very different to the recurrence time. We summarise this in Table 1.
It is then interesting to ask whether there is additional knowledge from photodetection statistics that one might be able to borrow for quantum walks. One potential application is the use of the waiting-time distribution for estimating an unknown parameter in a quantum emitter that one is interested in characterising [24,64]. The simplest model being a driven two-level atom whose Rabi frequency is not known to an experimenter. Kiilerich and Mølmer have shown that the waiting-time distribution of the photodetection record obtained from monitoring the atomic fluorescence is useful for estimating the Rabi frequency [64]. This idea may have applications to quantum walks where one has only a partially characterised graph-e.g. the frequency of one of the coherent transitions in the graph is not known. It might then be possible to identify a recurrent state whose distribution of recurrence times can be used for identifying the unknown transition frequency. Clearly, one will need the multichannel extension of this idea described in Ref. [24]. A successful application of this theory would then extend the quantum-walkquantum-optics analogy further.

A Independent derivation of the average hitting time
The goal here is to provide a derivation of the average of hitting times that is independent of the moment-generating function used in Sec. 4.2. This will add further confidence in our result. The ensemble average of T (ρ f , ρ i ) with ρ f =Q N is defined by We will again suppress the dependence of T and h on the initial and final states in the following. We integrate (139) by parts directly, in which case it will be convenient to define the indefinite integral We may then write The first term is given by We can get some insight into this function by actually doing the indefinite integral in (140). The function H(t) is Assuming thatL is invertible we can see that which in general has a value between zero and some finite number. Therefore we obtain lim t→0 t H(t) = 0 .
As for the first term in (142), we need to see how H(t) behaves as t −→ ∞ . This is given by where we have used (67), stated here again for ease of reference, This means that H(t) −→ 0 as t −→ ∞, which in turn implies that the first term in (142) has an indeterminate form of the type ∞ × 0. We can therefore use L'Hôpital's rule to calculate the limit, giving where the third equality is obtained from (140). For the final equality we assume that h(t) is decaying faster than t 2 in the limit t −→ ∞, which is the case for most physically relevant probability densities. We note here that it is not always the case that h(t) is such a "nice" distribution as this depends on the actual graph defining the quantum walk (e.g. a graph containing absorbing states may give rise to infinite hitting times). We will not consider such cases in this paper. Equation (142) thus vanishes completely and the average hitting time becomes From (145) we have The definite integral in (154) is thus given by Using again the fact thatL produces a zero steady state we have B Independent derivation of the variance of the hitting time The variance of the hitting time reads All we have to do is to calculate E T 2 (ρ f , ρ i ) which is given by We can proceed in a similar manner as we have done with the average. Integrating (162) by parts we get where we defined H(t) in (145). Repeated application of the L'Hôpital rule gives We have assumed again that h(t) is a normalisable function. Using again the fact that H(t) is a continuous function and finite at t = 0 we get This means that we are left with only the integral in (163).
where we have introduced The finiteness of F (0) and L'Hôpital's rule give once again We thus have E T 2 (Q N , ρ i ) = 2 (172) As with the average, we can replaceL −1 byL + to obtain (61).