Quantum walks and quantum search on graphene lattices

Quantum walks have been very successful in the development of search algorithms in quantum information, in particular in the development of spatial search algorithms. However, the construction of continuous-time quantum search algorithms in two-dimensional lattices has proved difficult, requiring additional degrees of freedom. Here, we demonstrate that continuous-time quantum walk search is possible in two-dimensions by changing the search topology to a graphene lattice, utilising the Dirac point in the energy spectrum. This is made possible by making a change to standard methods of marking a particular site in the lattice. Various ways of marking a site are shown to result in successful search protocols. We further establish that the search can be adapted to transfer probability amplitude across the lattice between specific lattice sites thus establishing a line of communication between these sites.


I. INTRODUCTION
Quantum walks have been a particularly fruitful field of research in quantum information going back to ideas from Feynman [1], Meyer [2], and Aharonov [3]. One focus of recent work on quantum walks is their application to spatial search algorithms. It is well-known that Grover's algorithm [4,5] for searching on unstructured databases offers a quadratic improvement in search time over classical models. Grover's algorithm is not designed to search physical systems where only local operations are possible and quantum walk algorithms have been employed for this purpose. In developing such algorithms, the major consideration is the search time attempting to match the quadratic improvement over the classical case offered by Grover's algorithm. The first spatial search algorithm to do this, using the standard model of quantum walks, was developed by Shenvi, Kempe and Whaley [6] for searches over a hypercube.
While there exist discrete-time searches on ddimensional square lattices which are faster than classical searches for d ≥ 2 [7] , effective continuous-time quantum searches only exist for d ≥ 4 [8] or else they require additional memory (in the form of spin degrees of freedom) in order to improve their search time in lower dimensions [9]. In [10], we demonstrated that effective searches over two-dimensional lattices may be achieved in an arguably simpler way which does not require extra degrees of freedom, and could, therefore, be viewed as more efficient. This is achieved through the choice of a different lattice, specifically, a honeycomb lattice which is the underlying lattice structure of carbon atoms in the material graphene.
The association with graphene is important as, although we first study a purely theoretical problem in quantum information, the use of a graphene lattice also offers a potential physical realisation. We thus envisage using the quantum walk and quantum search algorithm framework to investigate the effect of perturbations on the dynamics on graphene and other carbon struc-tures. This offers not only the possibility for demonstrating two-dimensional continuous-time quantum searching, but also paves the way for looking for novel effects in the material graphene.
In this paper, we offer a detailed account of the theory and numerical results thus expanding on the findings as presented in [10]. The paper is structured as follows: In Section II we give an introduction to the formalism of continuous-time quantum walks and the construction of quantum walk search algorithms. We also explain why previous search algorithms struggled on lowerdimensional lattices and why graphene offers a solution. Section III will detail the relevant properties of graphene and set-up the notation we shall use throughout. Sections IV & V contain our main analytical results, where we detail the specifics of our search algorithm and offer an analysis of the search running time and success probability. We shall then show in Section VI how this search protocol can be adapted to demonstrate novel communication setups. Sections VII & VIII contain numerical work demonstrating the possibility of using alternative methods of marking to create search behaviour and other carbon nanostructures. We conclude with a review and discussion of our results in Section IX.
Generalisations of the results in [10] have been given in [11], where a general approach to solving the problem has been taken through encoding extra degrees of freedom into crystal lattices. Previous work has also investigated the use of honeycomb lattices in discrete-time quantum searches [12]; however, due to the nature of the unitary evolution operator used there, no improvement over square lattices was found.

II. QUANTUM WALKS AND SEARCHING
Continuous-time quantum walks (CTQW), first defined in [13], are the quantum analogue of continuoustime Markov chains. They are defined purely on the state space, that is, the Hilbert space H p , spanned by the states |j which represent the j th site of the lattice. Thus, the time-evolution of such systems is defined by the Schödinger equation where α j = j| ψ (t) is the probability amplitude at the j th vertex of a system described by the state vector |ψ (t) , and H is the Hamiltonian describing the connectivity of the lattice. The dynamics of a quantum walk over a network is defined by the nature of the interaction between connected sites. Therefore, the Hamiltonian is generally constructed from the adjacency matrix of the underlying lattice. The adjacency matrix A is defined as if j and l are connected 0 if j and l are not connected.
Typically, the Hamiltonian is chosen as H = D I + vA.
The parameter v determines the coupling strength between connected sites and the parameter D is an on-site energy that only enters the dynamics in a trivial way and thus can be set to a desired value. If v = −1 and D is equal to the valency of the lattice H is the discrete Laplacian. This form of Hamiltonian is closely related to the tight-binding model for condensed matter systems [14].
As first explained in [15,16], a quantum walk is transformed into a search protocol by introducing a localised perturber state, forming an avoided crossing in the spectrum of the search Hamiltonian between an unperturbed eigenstate and the localised perturber. Thus, initialising the system in the unperturbed eigenstate involved in the crossing and allowing the system to evolve in time, one finds that the system rotates into the localised perturber state.
The first CTQW search over d-dimensional square lattices [8] introduced the localised perturber state as a projector onto a single site w, resulting in the search Hamiltonian Here, γ is a free parameter governing the strength of interactions between sites in the lattice. This parameter is chosen carefully [8] such that the perturber state |w is brought into resonance with the ground state of the unperturbed Hamiltonian, the uniform superposition |s . Let us assume for the moment that all other unperturbed eigenstates are energetically sufficiently separated from the ground state -we will come back to this assumption later. Then there will be an avoided crossing of two eigenvalues corresponding to the perturber state and the ground state. Perturbation theory estimates that the energy splitting of these resonant states is proportional to the overlap of the localised perturber state |w and the uniform ground state |s , which scales as and that the corresponding eigenstates are of the form (|s ± |w )/ √ 2. By preparing the system initially in |s and allowing it to evolve for a set period of time T = π /∆E ∝ √ N , a measurement of the system will result in the state |w being measured with high probability. This explains the speed-up of the quantum walk search. The Grover algorithm and its speed-up can be understood in analogous terms and one can show that all other states are indeed energetically very well separated.
In the present case this separation only holds above a critical dimension d c = 4. A simple argument for the critical dimension can be obtained by comparing the scaling behaviour of the energy splitting ∆E ∼ N −1/2 with the energy separation of the first excited state from the ground state in the unperturbed lattice. For a square lattice we have a quadratic dispersion relation which allows us to estimate the energy separation For d > 4 one then has ∆E/(E(k) − E 0 ) → 0 as N → ∞ and a detailed analysis indeed proves that the quantum walk search works with optimal speed-up [8]. At the critical dimension d = 4 the two energy scales scale in the same way -the detailed analysis shows that the search still works but the dynamics is more complicated due to the interference of excited states. While there is still a speedup it is only almost optimal; the optimal search time T ∝ N 1/2 gets multiplied with a logarithmically increasing factor. For the experimentally relevant regimes of 2and 3-dimensional square lattices, all states participate in the dynamics as any avoided crossing gets dissolved completely as the number of sites grows. As a result no speed-up over classical searches is found.
The estimate of the critical dimension above offers a simple way how to reduce the critical dimension. If one can construct a search around a uniformly distributed state at an energy E 0 where the dispersion relation is conic (linear in |k|), then in d dimensions E − E 0 = c |k| ∼ N 1/d which results in a critical dimension d c = 2. In [9] this was implemented using a modified Dirac Hamiltonian. However, this requires the addition of a spin degree of freedom, essentially a doubling of memory, with the added complication that it is not immediately clear how such a system would be physically realised.
Instead, our solution here is to change the lattice from a square to a graphene lattice. This change of topology automatically applies the first step in our solution, as one of the important electronic features of graphene is the conic dispersion relation around the Dirac energy. This arises naturally from the tight-binding descriptions of graphene. Note that graphene has the critical dimension. Indeed our construction in [10] has an almost optimal speed-up with logarithmic corrections that need to be evaluated in a detailed analysis that goes beyond the simple perturbative description given above. (Color online) Dispersion relation for infinite graphene sheet ( D = 0).

III. RELEVANT PROPERTIES OF GRAPHENE
Graphene is a single layer of carbon atoms arranged in a honeycomb lattice. The lattice is bipartite with two sublattices, labelled A and B, and a unit cell containing two carbon atoms. The spatial and reciprocal lattices are shown in Figure 1. The primitive vectors describing the lattice are a 1 (2) , such that the position of a unit cell in the lattice is given by R (α, β) = αa 1 + βa 2 . The reciprocal lattice shows two important points, the Dirac points K and K at the two inequivalent corners of the Brillouin zone.
The energy spectrum of electrons in graphene was first derived by Wallace [17] when considering the band theory of graphite using a tight-binding Hamiltonian. The tight-binding model for graphene and the derivation of the solution are well-known [18,19] and give rise to the dispersion relation shown in Figure 2 for an infinite graphene lattice. As there are two atoms per unit cell the spectrum has two branches, the upper branch being the conduction band and the lower the valence band, which meet at the corners of the Brillouin zone, the K-points. The energy at the K-points is D which we name the Dirac energy. It is around these points, that the behaviour of the spectrum is conical, that is, with a reduced density of states, a necessary feature for the creation of the search dynamics.
As the lattice possesses a translational symmetry the Hamiltonian can be solved using linear superpositions of Bloch functions over both sublattices. As a basis we use the orthonormal states {|α, β A , |α, β B } to denote states on either the A or B sublattice in the cell at position R (α, β). For the majority of what follows (except in Section VIII), we will focus on finite-sized lattices with assumed periodic boundary conditions along the axes of both primitive vectors so that the topology of our lattice is a torus; that is, our wavefunction is of the general form |ψ = m α=1 n β=1 ψ A α,β |α, β A + ψ B α,β |α, β B . Our boundary conditions imply that the state vector must satisfy ψ α,β+n where m, n denote the period of the lattice. Thus, the wavefunctions on the torus take the Bloch function form where N is the number of sites in the lattice, k is the momentum, and C (k) is a relative phase contribution dependent upon whether the state belongs to the conduction or valence bands; it can be calculated by explicitly working through the tight-binding model. Application of the periodic boundary conditions results in the following quantised momenta where p ∈ {0, 1, . . . m − 1} and q ∈ {0, 1, . . . n − 1}. In the following and whenever we consider quantum walk dynamics on a torus, the number of cells in each direction is generally chosen to be the same, that is, m = n = N 2 . This choice is purely for simplifying the notation; alternative torus dimensions are possible as are other choices of other boundary conditions corresponding to alternative carbon structures (e.g. nanotubes or a graphene sheet) as will be demonstrated in Section VIII. For our choice of torus dimensions, we find there are momenta equal to the K-points and, consequently, eigenstates with energies equal to the Dirac energy when both m and n are some multiple of 3.
In fact, using the quantised momenta in Eq. (8) obtained for periodic boundary conditions, we find that there are four degenerate eigenstates with an energy that coincides exactly with the Dirac energy when m and n are both multiples of 3. These four states, known as the Dirac states, can be constructed to live only on one sublattice, and are given by where σ = 0 for states on the A-sublattice or σ = 1 for states on the B-sublattice.

IV. QUANTUM SEARCH ON GRAPHENE -REDUCED MODEL
In this section we will describe how a site is marked and will derive our optimal search starting state. Our approach will be to analyse the system's spectrum and its dynamics in a reduced Hamiltonian model involving only the relevant states from the avoided crossing. In the next section we will then validate this approach with a more detailed analysis, using the results obtained here as an initial guide.
As already established, we introduce the localised perturber state in a region of the spectrum with a conic dispersion relation and thus a low density of states. Simply altering the on-site energy in Eq. (3) as done in [8] does, however, not work here. As the on-site energy D and the Dirac energy are equal, one finds that the perturbation only interacts with the Dirac states in the limit of zero perturbation strength, returning the unperturbed lattice. Therefore, we choose an alternative perturbation method: namely, we modify the coupling strength between the site we wish to mark and its neighbouring vertices. We focus in this section on changing the coupling to all three neighbouring vertices of a particular site equally. Our choice of perturbation matrix W to mark the A-type vertex (α o , β o ) A is then where the state | is the symmetric superposition over the three neighbours of the perturbed site where γ is a free parameter. In what follows, we always set the on-site energy of sites in the lattice to D = 0. Considering our search Hamiltonian, we can see that setting γ = 1 corresponds to a coupling strength of v = 0 from the perturbed vertex and its nearest-neighbours; A from the lattice. Note that vacancies are a common, naturally occurring defects in graphene lattices [21].
In order to establish the critical value of γ, we numerically calculate the spectrum of H γ as a function of γ, plotted in Figure 3 for a torus of dimensions m = n = 12. As W is a rank-2 perturbation, we see in Figure 3 two perturber states working their way through the spectrum to an avoided crossing around D = 0 when γ = 1, that is, when the perturbed vertex is removed from the underlying lattice.
As well as establishing a critical value for γ, we also find from this figure the states involved in the avoided crossing: the four Dirac states at the Dirac energy and the two perturber states which form our perturbation matrix W. However, with further consideration, we may reduce the number of states involved further. As we have removed the site (α 0 , β 0 ) A from the lattice it can no longer interact with the rest of the lattice and the corresponding state drops out. Also, by direct calculation one finds that the B-type Dirac states do not interact with the perturbation, that is, Thus, they remain an eigenstate of the search Hamiltonian and do not interact with the perturbation. We are left with three states taking part in the avoided crossing: We reduce our search Hamiltonian in this three state basis at the critical point γ = 1 to obtain the following reduced Hamiltonian describing the local dynamics at the avoided crossing Using the eigenvectors of the reduced Hamiltonian, we can construct a search starting state which is a superposition of Dirac states Allowing our search starting state |s to evolve under the reduced Hamiltonian we find = cos Ẽ + t |s − i sin Ẽ + t | , so that our system rotates from our Dirac superposition |s to a state which is localised on the neighbours of the Thus, we find a O √ N search time, a polynomial improvement over the classical search time.
We plot in Figure 4 a numerically calculated search for a lattice with N = 288 sites. The system has been prepared in |s and allowed to evolve under the full search Hamiltonian from Eq. (12) at the critical value γ = 1. One finds that the system localises on the neighbours of the marked vertex and, for this particular system size, the probability for a measurement of the system to return one of the neighbours is around 45% at its peak. This is orders of magnitude higher than the average probability to measure other vertices in the lattice, 100/N , which for this system size is around 0.35%. The localised state interacting with the spectrum under the full Hamiltonian has a tail into the rest of the lattice; this leads to a loss of probability to be found at the neighbouring vertices below 100%.
It is worth emphasising that the marked vertex plays no role in the search -it is removed from the lattice. Rather, we force the system to localise on the nearestneighbours, essentially making the marked vertex conspicuous by its absence. This feature differentiates our algorithm from previous searches, which have generally concentrated on localising directly on the marked vertex. An approach considering the neighbouring vertices was detailed in [26], however, the focus was here still on localising on a single site but with classical post-processing steps considering the tail of the localised state extending into the neighbourhood around the marked vertex.
At this point, there are a number of issues which need to be mentioned. The first is that our starting state in Eq. (16) contains information about the marked vertex in the form of the relative phase between the Dirac states. However, this phase can only take three values. Therefore, we have three optimal starting states for an A-type perturbation and the same is true for B-type perturbations; there are thus in total six possible optimal states. As not all of these states are orthogonal, we only find an increase of necessary runs for a successful search by a factor of 4; this additional overhead is independent of N and, therefore, does not affect the overall time complexity of the search. The particular representation of the Dirac states in Eq. (9) was chosen in such a way as to make the calculations and conceptualisation of the system dynamics easier. In reality, constructing starting states experimentally which exist on one sublattice only is extremely difficulty; rather, in an experiment, one is likely to excite a superposition of Dirac states, reducing the success probability, on average, by a factor of 1/4.
Our search is based on the conic dispersion relation in the spectrum and the O 1/ √ N scaling between successive energy levels in the linear regime, giving rise to the O √ N search time found in our reduced model. However, previous searches in continuous-time using a modified Dirac Hamiltonian [9] or in discrete-time [7] have found logarithmic corrections to the search time. As we have seen in Eq. (17) the search time is related to the energy gap at the avoided crossing. In the inset of Figure 3 we have numerically calculated the scaling of the gap in the spectrum of the full search Hamiltonian and we indeed find evidence for a logarithmic correction. The lack of a logarithmic term in the search time for our reduced model is due to the neglect of contributions from the rest of the spectrum.
In what follows we establish a more accurate estimate of the running time and success probability by working with the full Hamiltonian. We also justify the findings of our reduced model: while this model is insufficient to estimate the finer details, such as logarithmic correction terms, the accurate calculations show that our search algorithm indeed takes place mainly in a two-dimensional subspace spanned by the Dirac states and a state localised on the neighbours of the marked vertex.

V. QUANTUM SEARCH ON GRAPHENE -DETAILED ANALYSIS
We follow here closely our derivation as presented in [10] which builds on ideas given in [8] for a similar a similar calculation for regular rectangular lattices. We focus first on the search success amplitude, that is, where T is the search time, that is, the time our search probability reaches a maximum, and our search Hamiltonian H is given by Eq. (12) with γ = 1, that is, here, |ψ a , E a are the eigenstates and eigenenergies of H. We assume, without loss of generality, a specific starting state |start where the marked vertex is chosen such that e i 2π 3 (αo+2βo) = 1 leading to In the following, we suppress the sublattice superscript as the analysis is the same regardless on which sublattice the perturbation lives. We also denote (k) the positive eigenenergies of −A from Eq. (5) and set D = 0. For eigenstates |ψ a with eigenenergies E a that are not in the unperturbed graphene spectrum (E a = (k) for all points k in the dual lattice), we may rewrite Eq. (19) in the form where √ R a = | ψ a and the phase of |ψ a is chosen such that | ψ a ≥ 0.
At this point, we can remove several states from the summation in Eq. (18). Note that the basis state associated with the marked vertex, |α o , β o , is itself an eigenstate of the search Hamiltonian with H |α o , β o = 0, and also, | α o , β o = 0 so that |α o , β o does not contribute to the time-evolution. We may also remove eigenstates of H which are at the same time eigenstates of −A with the same energy.
This can be seen in the following way. We first consider an unperturbed eigenstate |ψ o a such that −A |ψ o a = E a |ψ o a . Let us assume that there is an eigenvector |ψ a of the search Hamiltonian with the same eigenenergy E a , that is, H |ψ a = E a |ψ a . Considering the It is then clear that eigenstates of the search Hamiltonian whose eigenenergies remain in the spectrum of −A do not play a role in the time-evolution of the search.
We now derive an eigenvalue condition for those perturbed eigenvalues which are in the spectrum of H. Using the orthogonality of eigenstates By expressing α o , β o | in terms of the eigenstates of the unperturbed walk Hamiltonian −A, we may write this as a quantization condition were N is the total number of sites in the lattice. Note this is written as two summations to incorporate the fact that the spectrum of −A as well as H is symmetric around E = 0.
Choosing |ψ a to be normalised ψ a | ψ a = 1, Eq. (21) also implies which allows R a to be rewritten as We may now rewrite the amplitude in Eq. (18) in the form where we have used the adjoint of Eq. (21). Note again that eigenstates of H, which are also in the spectrum of −A, do not contribute to the time-evolution of the search; it may be seen from the definition of As we saw in the previous section, the avoided crossing is formed by two perturbers approaching symmetrically either side of the Dirac states with energy D = 0. Thus, we concentrate on evaluating the contribution to the time-evolution from these perturbed eigenstates of H either side of the Dirac point and we label these states |ψ ± . In order to evaluate these contributions we calculate the perturbed eigenenergy E + and also derive a leading-order expression for F (E + ) (since the spectrum is symmetric it is true that Using the definition of F (E) in Eq. (23), we estimate F (E + ) by separating out the Dirac points, where (K) = K = 0, from the summations and then Taylor expanding the remaining terms at E = 0, to find where the sums I n are given by As the unperturbed spectrum is symmetric only those I n with even n are non-zero, thus from now we focus only on I 2k where k ≥ 1.
We stated earlier that the spectrum of graphene is wellapproximated by a conic dispersion relation around the Dirac points. Thus, it is from around these points that the major contributions to the I 2k sums arise. By applying the linear approximation from Eq. (6) and the momenta quantum numbers from Eq. (8) we may approximate the I 2k summations as Here Z 2 (S, x) is the Epstein zeta-function [20] Z 2 (S, x) = 1 2 (p,q)∈Z 2 \(0,0) for a real positive definite real symmetric 2 × 2 matrix S. For our purposes we use which describes the spectrum close to the Dirac points. The linear approximation around both Dirac points K and K is the same and, thus, the matrices S K and S K are equal. Our summation in Eq. (30) is over the rectangular regions L and L of the lattice Z 2 . Both are centered on (0, 0) and have side lengths proportional to √ N , however the center, (0, 0), corresponding to the relevant Dirac point, is omitted from the summation.
Convergence of the Epstein zeta function is well-known for k ≥ 2 [20], and leads to the bounds lim N →∞ for k ≥ 2. A sharp estimate for is given in Appendix A. In order to calculate an estimate for E + , we truncate the Taylor expansion of F (E) in Eq. (28) at the I 2 term (the first term in the summation), and apply the eigenvalue condition F (E + ) = 0. That is, we solve 4 √ 3 N E+ − I 2 E + = 0 and obtain the ap- . This solution also leads to the estimate F (E + ) ≈ −2I 2 .
We consider next whether our solution for E + lies inside the radius of convergence of the Taylor expansion. We note that each term in Z 2 S K , k is smaller than the corresponding term in Z 2 S K , 2 for k > 2, and so it follows that Z 2 (S K , k) < Z 2 (S K , 2) for k > 2. This property of the Epstein zeta function and the estimate from Eq. (33) imply i.e. the infinite sum in Eq. (28) converges for E + < 1/ √ N . Thus, for large N , our solution lies within the radius of convergence of our Taylor expansion in Eq. (28).
In order to show that the dominant contributions to the search come from |ψ ± , we need establish the leading order error term. All the I 2k sums are positive so that, given the sign of all the terms in the Taylor expansion in Eq. (28), the true value of E + > 0 has to be smaller than the estimate we have obtained. Thus, we write the true value of E + as with ∆ > 0. One may rewrite F (E + ) = 0, using the true value of E + , to give We follow the same arguments as used for the calculation of the radius of convergence to obtain an upper bound for the summation. This together with the already established fact that E + is inside this radius for sufficiently large N , we get the following inequality: Thus, we obtain It also follows that This shows that our perturbed eigenstates |ψ ± have an O (1)-overlap with the starting state and are, therefore, the relevant states to be considered in the timeevolution of the algorithm. Using the definitions of |ψ a and R a in Eqs. (21), (26), the inner product of the starting state and the perturbed eigenvectors can be expressed as where start| α o , β o is the overlap of the starting state with the marked vertex state. Applying our previous approximations for E + and F (E + ), that is, for our perturbed eigenstates closest to the Dirac point, we find Thus, our starting state is indeed a superposition of the perturbed eigenstates, |ψ ± as assumed in the previous section. This makes it possible to investigate the running time and success amplitude of the algorithm in more detail, that is, It is clear from our earlier results for E ± and I 2 that our algorithm localises on the neighbour state | in time This confirms the logarithmic correction observed numerically and displayed in the inset of Figure  3.

VI. COMMUNICATION
It has been demonstrated in [16] that discrete-time search algorithms can be modified to create a communication protocol. We show here that a communication setup can be established also in the continuous-time search algorithm with minor changes due to the subtleties of our search.
We use the same unperturbed walk Hamiltonian as before, that is, H o = −A and the same type of perturbation matrix as in Eq. (10), but now at two different sites in the lattice. Explicitly, our communication Hamiltonian is where the perturbation matrices W s/t mark the source and target sites, respectively.
As prescribed in [16], the communication protocol operates by preparing a state localised on the source perturbation, in our case on the neighbours of the source site, and allowing the system to evolve under our communication Hamiltonian. The system then localises on the neighbouring vertices of the target site. The effect of the two perturbation matrices W s/t amounts to disconnecting the two vertices from the underlying lattice.
As noted before, our search algorithm has several different optimal starting states depending on the position of the marked vertices. Therefore, we divide our communication analysis into three different cases: the two sites i) are on the same sublattice and have the same optimal search starting state; ii) are on the same sublattice but do not share the same optimal search starting state; iii) are on different sublattices. Both of the first two cases can be treated by applying the reduced Hamiltonian method demonstrated earlier, but communication between different sublattices is not tractable by this method and so we focus on numerics in this case.
For signal transfer between two sites on the same sublattice, we will assume, without loss of generality, that our communication takes place between two sites on the A-sublattice, α s/t , β s/t A ; here, the subscript s or t denotes the source or target vertices respectively. Using arguments as employed in Sec. IV, we can reduce the number of relevant states. Ultimately, the search dynamics takes place in the subspace spanned by the basis |K where µ s/t ≡ 2π 3 α s/t + 2β s/t and ν s/t ≡ 2π 3 2α s/t + β s/t . In the first case considered we assume that there are two perturbations on the graphene lattice located at the points (α s , β s ) A and (α t , β t ) A , chosen such that e i 2π 3 (αs+2βs) = e i 2π 3 (αt+2βt) . This implies that a search for either vertex, using our search algorithm from Secs. IV & V, would use the same optimal starting state. As the phases are equal, we drop the subscript on the perturbation coordinates in what follows. Fixing our phases and diagonalising the reduced Hamiltonian, we find it has eigenvalues λ ± 2 = ±2 6 N and λ 1,2 0 = 0 with eigenvectors Using these eigenstates, we may rewrite the source neighbour state as Placing the system in the source state | s and allowing the system to evolve under the reduced Hamiltonian, one finds We note that, in the last line above, the term in the brackets involving only the Dirac states is actually the optimal search starting state for both perturbed vertices, defined in Eq. (16). The system thus oscillates between the states localised on the neighbours of the perturbed vertices, | s and | t , in a time T = π 2 N 6 , via their optimal search starting state. Figure 5 shows the system evolving under the full communication Hamiltonian, Eq. (44). The initial state used for the time evolution shown in Figure 5 is the true localised state, that is, we run the quantum search with a single perturbation located at vertex (α s , β s ) A until it reaches maximum success probability, and then apply the second perturbation to the vertex (α t , β t ) A . The figure confirms the behaviour expected from the reduced model calculation.
The communication mechanism essentially works in the same way as the quantum search algorithm, as it can be viewed as one marked vertex 'finding' another. The initial localised source state decays back towards the search starting state, and the system then searches for the target state.
In our second case of signal transfer we analyse the behaviour of a communication system where we have two perturbed vertices on the same sublattice, but with the restriction that e i 2π 3 (αs+2βs) = e i 2π 3 (αt+2βt) . As such, the two marked sites cannot interact via the same search starting state. However, the optimal search starting states are not orthogonal so that signal transfer is still possible, but the resulting interference effects make the analysis slightly more complicated.
We rewrite the coordinates of the target in terms of the source, α t = α s + x and β t = β s + y. Again, fixing our phases and diagonalising the reduced Hamiltonian in Eq. (45), we find it has eigenvalues λ ± √ 3 = ± √ 3 6 N and Using these eigenstates, one may write the source perturbation as The full expression for the time-evolution is rather cumbersome; therefore, we only show the terms and prefactors we are interested in, namely | s and | t |ψ (t) = 1 2 cos λ + √ 3 t + cos λ + 1 t | s (56) We can see here that the prefactors do not depend upon the coordinates of either the source or the target sites; the transport signal between sites on the same sublattice but with different optimal search starting state is thus independent of the position of the source or target site. In Figure 6 we show the system evolved under the full communication Hamiltonian. Again, the initial state is the true localised state on the nearest-neighbours of the source vertex, obtained by running the search algorithm using one marked vertex until it reaches its peak success probability. The time-evolution is radically different to the previous communication case; the behaviour here is erratic with uneven peaks of probability at the two perturbations involved in the protocol. However, there are still significant probability revivals.
The transport behaviour from Eq. (56), calculated using the reduced Hamiltonian, is also shown in Figure 6. The probability at time t = 0 has been scaled to match that shown in Figure 6. Our calculated behaviour has the same signal pattern as the numerically calculated behaviour from the full Hamiltonian, although over a shorter period of time. As our reduced model only makes use of the Dirac states and the perturber states, we lose the contribution to the time-evolution from the rest of the spectrum giving rise to logarithmic corrections as discussed in Sec. V. This leads to differences in the overall time scales of for the reduced model and the full Hamiltonian, however, the reduced model captures the signal pattern well. This also supports our findings that signal transfer between all non-equivalent vertices on the same sublattice is the same. The communication protocol is essentially again set up by a the search mechanism in reverse, where one vertex finds another. The slightly erratic behaviour that emerges here is due to interference between the two separate search mechanisms which interact due to the non-zero inner product of the three possible optimal search starting states for vertices on the A-sublattice.
In our final case of signal transfer, we consider communication between the two sublattices which we can not treat in the reduced Hamiltonian model. It has been demonstrated previously that perturbations to one sublattice do not interact with the Dirac states which live on the other sublattice. Therefore, when attempting to reduce the communication Hamiltonian as done before, we merely find that it decouples into two non-interacting reduced Hamiltonians, each describing a search protocol on one sublattice. The interaction between the two sublattices is here facilitated due to interactions via the bulk of the spectrum. In what follows, we focus on numerics and inspect the behaviour for systems involving  Figure 7. We can see from these examples that the communication mechanism takes place over a much longer timescale with a superimposed oscillatory dynamics.
The increase in timescale can be attributed to the weak nature of the interaction between the two perturbations due the fact that the localised states live mainly on one sublattice. The signal pattern in Figure 7 thus decouples into a fast oscillation between one site, say the source site, and the delocalised lattice state and a slow time scale on which a small amount of probability amplitude escapes into the other sublattice due to the weak interaction. This process continues until the recurrence probability at the target perturbation reaches the same peak as the initial localised source state, and then the behaviour reverses.

VII. ALTERNATIVE METHODS OF MARKING
In the previous sections, we have discussed a method of marking a vertex through modifying the hopping potential from the perturbed site to all three of its nearestneighbours. Here, we will discuss alternative methods of marking a vertex which still keep the perturber interaction in the conic dispersion region of the spectrum. This can be done by altering the hopping potential in different ways as discussed above; this approach is necessarily a rank-2 perturbation and these two perturber states must meet at the Dirac energy. We will discuss several ways of applying such a perturbation in the next paragraph. Another approach is based on coupling extra sites to the lattice; this set-up will be discussed at the end of this section.
Focussing on hopping potential perturbations, we may consider perturbing the bonds to any number of nearestneighbours to any strength. We have seen previously that a symmetric three-bond perturbation successfully creates a search. Here , we demonstrate a search based on perturbing the hopping potential from a given site to only one of its nearest-neighbours. That is, we use the same search Hamiltonian as in Eq. (12) with the perturbation matrix and eigenstates Note that it no longer makes sense to speak of a single marked vertex. Rather, our perturbation marks both vertices (α o , β o ) A/B simultaneously; thus, it may be more accurate to view our perturbation matrix as marking a single cell of the lattice.
For this type of perturbation, the avoided crossing used to generate search behaviour is not necessarily at γ = 1, see Figure 8, where the spectrum of our search Hamiltonian as a function of γ for a 12 × 12 cell torus is shown. As stated before, our single bond perturbation is also a rank-2 matrix, and so again there are two perturber states approaching the spectrum from the negative and positive regions of the spectrum. Inspecting the region around the Dirac energy, see inset in Figure 8, we see that avoided crossings are formed by four states in total, the two blue curves and the two red curves. At γ = 1 3 there is an exact crossing between the red and blue curves indicating that these states are orthogonal.
We can make this avoided crossing picture clearer by breaking down the eigenstates of the search Hamiltonian in terms of the symmetries of the lattice. In particular, we use the rotation operator C 2 , which is a rotation of π about the mid-point of the perturbed bond. Specifically, considering the action of C 2 on the marked vertices and the eigenstates of the perturbation matrix, one finds C 2 |α o , β o A/B = |α o , β o B/A , C 2 |W g = − |W g and C 2 |W e = |W e . The inset of Figure 8 shows the states nearest to the Dirac energy, where the same colour indicates the same C 2 parity. It now becomes clear that the states near the Dirac point actually form two avoided crossings, one for each parity with a minimum energy gap at γ = 1 3 . Without going into details, we can again find optimal starting states using a reduced Hamiltonian approach (one for each of the avoided crossings); the evolution of the system using these optimal states is shown in Figure 9  for a 12 × 12 cell lattice. We find there is a significant localisation on the marked vertices and their nearestneighbours, with the probability of being found on either of the marked vertices peaking at around 16 − 18%. We can also see that the probability of being found on each of the nearest-neighbours peaks at around 8%, resulting in a total probability of being found on the marked vertices and their nearest-neighbours of approximately 48%. The success probability fluctuates, as the probability amplitude oscillates between the marked vertices and their nearest-neighbours. This is due to the probability amplitude being constrained in the local area by the increased hopping potential between the two marked vertices.
Although our demonstration is for a 12×12 cell lattice, numerical investigations show that the search behaviour remains the same as the lattice size increases with critical value fixed at γ = 1 3 . We also find that the gap at the two avoided crossings scales as O 1/ √ N ln N , in the same way as for the three-bond perturbation shown in Figure 3. As the search time is inversely proportional to the energy splitting at the avoided crossing, it also gives an estimate of the running time of the search T = O √ N ln N .
We now turn to coupling additional sites to the lattice as a way of introducing a perturbation; such a treatment is in many ways closer to experimental realisations as defects due to additional add-on atoms are quite common in graphene [22]. We focus here on the idealised case where an additional site is coupled to a single lattice vertex. We use the perturbation matrix for coupling the additional site |site to the A-type vertex (α o , β o ) A and γ is a free parameter related to the onsite energy of the additional vertex. The choice of the coupling terms is such that the binding energy between the additional site and the lattice vertex is the same as the internal couplings in the lattice. Thus, our search Hamiltonian is of the form In Figure 10, the spectrum of the search Hamiltonian is given as function of γ. As our free parameter γ only changes a single term, our spectrum only has a single perturber state. One finds a clear avoided crossing around the Dirac energy (E = 0) when γ = 0, that is, when the on-site energy of the additional site matches the on-site energy of the lattice vertices.
In the reduced Hamiltonian picture, we write Eq. 61 in terms of a basis consisting of the Dirac states and our perturber state, |site . Only three states are involved in the search, that is, {|K A , K A , |site }, and H reduces, up to a different prefactor, to the same 3×3 matrix found in Eq. (13) for the three-bound perturbation. Thus, our previous reduced Hamiltonian analysis holds for this case and the optimal search starting states are the same; the time-evolution for the additional site search is shown in the inset in Figure 10. The only differences between this case and the three-bond marking are that the system localises on the additional site, and also the change in prefactor leads to a search time of T = π 4 √ N . We note that other, more realistic types of perturbations involving additional sites coupling to a lattice site and its neighbours, can be shown to result in effective search protocols. Also, the single additional site perturbation described here can, like the three-bond perturbation, be used to setup a communication protocol, with the same signal patterns and behaviours as found previously. The initial state is in this case completely localised on the additional site.

VIII. ALTERNATIVE NANOSTRUCTURES
So far we have described the dynamics of searches on graphene lattices with periodic boundary conditions, that is, graphene on a torus. Here we consider more realistic boundary conditions such as nanotubes and graphene sheets.
For the dynamics on nanotubes, we move to periodic boundary condition along one axis only and impose Dirichlet boundary conditions along the other directions. The properties of nanotubes are well-known [23], and it has been shown that the band structure of armchair nanotubes, that is, nanotubes with armchair boundaries, always allows for an energy at the Dirac energy regardless of the nanotube diameter; we will focus on these types of nanotubes. We are interested in searching on finite length nanotubes, where the band structure becomes discrete; the length of the nanotube are in addition chosen such that there exists an eigenenergy at the Dirac energy. An example of the cell we use to construct finite armchair nanotubes is shown in Figure 11. We choose the finite length of the nanotube to be along the horizontal axis and we close the underlying graphene lattice into a nanotube along the vertical axis.
We denote the basis states of sites in our nanotube as |m, A/B, l , where m indicates the m th A/B-type vertex in the horizontal direction in the l th cell. It is simple to see that, due to the periodic boundary conditions along the circumference of the tube, the vertical component of the eigenstates must be Bloch states and the horizontal component of the amplitudes must be sinusoidal in nature. That is, the eigenstates of the finite nanotube are standing waves along its length.
By working through the tight-binding model for this system, one can show that (k x , k y ) = 2π 3 , 0 is the only point where there exists an eigenenergy equal to the Dirac energy. It is also possible to show that the spectrum includes Dirac points when the number of sites along the length of the nanotube N x = 3r−1, where r is an integer.
As there is only one potential Dirac point for finite armchair nanotubes, it follows that there are only two Dirac states (one from the bonding and the anti-bonding regions of the spectrum). As we have k y = 0 at the Dirac point, the Bloch wave around the circumference of the nanotube is simply a uniform superposition. Another important feature of the Dirac states on the nanotube, which we have not encountered previously, is the existence of nodal points where the amplitude of the eigenstate is 0. We find these nodal points occur at every third site along the horizontal axis.
We focus on applying our three-bond perturbation from Sections IV & V to the armchair nanotube. In our new labelling our perturbation takes the form We have assumed that we are perturbing an A-type vertex on an even horizontal coordinate, so that the perturbed site and its nearest-neighbours remain within one nanotube cell. While it is easy to re-write this perturbation matrix for other sites we restrict ourselves to this form for simplicity. Note even in this form it can still be expressed in terms of a marked site,  (12). Inspecting the spectrum of the search Hamiltonian as a function of the free parameter γ where the perturbation is not located on a nodal point of a Dirac state, one finds an avoided crossing around the Dirac energy when γ = 1. This spectrum is very similar to the one shown in Figure 3 for the search on the torus and will be omitted. If the perturbation is located on a nodal site the picture is different, however; one finds that there is an exact crossing at the Dirac energy when γ = 1 and searching is not possible. Recall that the effect of the perturbation at the critical point is to completely remove the site from the lattice, if the site already has zero amplitude then the perturbation will not interact with the Dirac state.
Similar to the previous section, we numerically reduce the full search Hamiltonian in a basis consisting of the two Dirac states and the state living purely on the neighbours, | . The perturbed site, |m o , A, l o , is decoupled from the lattice. Through this process, we find two possible starting states, one for each sublattice. The starting states are weighted superpositions of the Dirac states over the sublattice containing the marked vertex. The nodal points are excluded from this treatment. Figure 12 shows the system evolution for two searches using the numerically found optimal starting states. One search has the marked site located in the centre of the nanotube, the other has the marked site positioned near the edge of the nanotube. The dimensions of the nanotube have been chosen so that there are 320 sites, comparable to the searches shown in earlier sections. Comparing these searches to the three-bond perturbation search on the torus in Figure 4, we see that there is a marked difference in behaviour and success probability induced by relaxing the periodic boundary conditions along one axis.
Using a perturbation located near the edge, the lower image in Figure 12, we see a reduction in success probability by a factor of 2-3 and strong fluctuations in peak height when compared to the search on the torus. Searching in the bulk, shown in the lower image of Figure 12, displays behaviour closer to searches on a graphene torus, but again with significant fluctuations at each peak. We propose that this effect is due to the reflection of probability amplitude from the edges of the nanotube. This is supported by the changes in the interference pattern in the signal as the perturbation is moved across the lattice. Numerics demonstrate that additional site perturbations and communication protocols can also be used on nanotube lattices.
We now move to working on graphene sheets, that is, removing the periodic boundary conditions along both axes. A detailed account of a specific version of this setup together with experimental results has been presented in [30] for graphene sheets consisting of armchair boundaries only. Note that such a configuration can not be achieved on rectangular sheets; the simplest configuration has the form of a parallelogram, see [30] for further details. The advantage of such a configuration is that the form of the boundary does not admit so-called 'edgestates'. In the following, we will look at the influence of these edge-states in more detail by considering rectangular graphene sheets, such as in Figure 13. In addition to armchair edges, here along the vertical axis, these sheets have boundaries formed by bearded edges (left) and zigzag edges (right) along the horizontal axis. These boundaries support edge states, i.e. states localised along these edges, with an eigenenergy close to the Dirac energy [24,25]. In the following, we will investigate, how the existence of these edge states influences the search.
As in the case of the finite armchair nanotube, the imposition of Dirichlet boundary conditions at an edge generates sinusoidal eigenstates. As a result, there are no extended (bulk) eigenstates at the Dirac energy due to the inability to equate the quantised momenta with the necessary points in k-space. We thus need to find other extended eigenstates in the sea of edge states near the Dirac point to undertake a search in this set-up.
We mark sites using the triple-bond perturbation such as in Eq. (12). Throughout this section we choose the dimensions in terms of primitive cells of the graphene sheets, (N x , N y ) = (10, 10); a bearded edge sheet thus consists of N = 200 sites and a zigzag sheet is formed of N = 218 sites. The spectra of the search Hamiltonians for the bearded lattice is shown in Figure 14 constrained to the energy region of interest. The spectrum for zigzag edges is very similar and is not show here. One can clearly see an avoided crossing around the Dirac energy at γ = 1, the critical value for this type of perturbation. There are several states very close to the Dirac energy; these are the edge states which are all non-degenerate. Therefore, it is not possible for us to construct a superposition of degenerate eigenstates which is optimal for searching. Rather, our initial starting state must be a single eigenstates of the unperturbed Hamiltonian H o = −A.
Using the search Hamiltonian from Eq. (12) and fix- ing γ = 1, we proceed by allowing the system to evolve after being prepared in one of the unperturbed eigenstates. One finds for both types of lattice, that the edge states near the Dirac energy fail to produce localisation behaviour; the probability at the neighbouring vertices of the marked site do not rise above noise levels. Rather, search behaviour only begins to emerge when we use the first delocalised, non-edge state as our initial state. We also note that we cannot search for sites near the zigzag or bearded edges, where the edge states exist. Only as we move further into the bulk of the lattice or along an armchair-type edge, the localisation behaviour returns. Note that the edge states can be viewed as a kind-of one-dimensional system and, as we saw from the scaling argument towards the end of Sec. II, one-dimensional systems imply an energy spacing between successive energy levels of E n+1 − E n = O N −1 . Thus, the perturber state interacts with many states in a dense part of the spectrum and the search fails.
In Figure 15 we show the search behaviour arising from marked sites placed at different positions on a bearded graphene lattice (zigzag lattice types display similar behaviour). The results are similar to those found for the nanotube searches. There is a slight increase in variation of signal pattern with position, and an increase in success probability maxima as we move towards the centre of the lattice from either of the bearded edges.

IX. DISCUSSION
We have shown that continuous-time quantum search can be done effectively on a two-dimensional graphene lattice without the use of internal degrees of freedom. This is achieved by making use of the conical (linear) dispersion relation in the graphene spectrum. The N . Amplification methods [27][28][29] may be used to reduce the total search time further.
Our main result focusses on perturbations which involve altering the hopping potential from a marked site to all three of its nearest-neighbours equally. We have also demonstrated other types searches based on perturbing the hopping potential in a cell and the adding extra sites. We have shown that search mechanisms can be utilised for the purposes of signal transfer.
Our findings point towards applications in directed signal transfer, state reconstruction, or sensitive switching. This opens up the possibility of a completely new type of electronic engineering using single atoms as building blocks of electronic devices. Our results demonstrate that a range of nanostructures constructed from graphene could be used to this end. Recent experiments on artificial microwave graphene [30] have shown that additional site perturbations can be used to create a switching protocol.

(A1)
Here the sums over integers p and q is over a rectangular region L of the lattice Z 2 which is centered at (0, 0) and has side lengths proportional to √ N -the center (0, 0), corresponding to the relevant Dirac point, is omitted from the sum. As stated in the main text, for k > 1 the corresponding sums converge which proves Eq. (33). For k = 1 we will establish constant C 1 and C 2 such that C 1 ln N < (p,q)∈L 1 S K,11 p 2 + 2S K,12 pq + S K,22 q 2 < C 2 ln N (A2) which then directly leads to Eq. (34). To establish C 1 note that because each term in the sum (A2) is positive its value decreases by restricting it to a square region −a 1 √ N ≤ p ≤ a 1 √ N , −a 1 √ N ≤ q ≤ a 1 √ N which is completely contained in L. Up to an error of order one the sum over a square region can in turn be written as a sum over eight terms of the form a1 √ N p=1 p q=1 1 S K,11 p 2 ± 2S K,12 pq + S K,22 q 2 . (A3) For fixed p we can find q max such that p q=1 1 S K,11 p 2 ± 2S K,12 pq + S K,22 q 2 > p S K,11 p 2 ± 2S K,12 pq max + S K,22 q 2 max . (A4) We may choose q max = b 1 p for some constant b 1 ≥ 0, so a1 √ N p=1 p q=1 1 S K,11 p 2 ± 2S K,12 pq + S K,22 q 2 > c √ N p=1 1 p (A5) which diverges as ln N .
Establishing C 2 and the corresponding logarithmic bound from above, and thus proving Eq. (34), follows the same line by first extending the sum to a square of side length 2a 2 √ N that completely contains L and then establishing p q=1 1 S K,11 p 2 ± 2S K,12 pq + S K,22 q 2 < p S K,11 p 2 ± 2S K,12 pq min + S K,22 q 2 min (A6) with q min = b 2 p.