Entangled photon factory: How to generate quantum resource states from a minimal number of quantum emitters

,


I. INTRODUCTION
Entanglement is widely recognized as playing a critical role in quantum computation, error correction, communication, and sensing. A family of entangled states that features prominently in these applications are graph (or cluster) states. They are key resources in one-way quantum computing paradigms [1,2] and in quantum error correction [3][4][5][6]. In addition, many quantum repeater schemes [7][8][9][10][11] and quantum sensing protocols [12,13] rely on graph states. Photonic graph states are especially important because photons are the predominant platform for measurement-and fusion-based computing, and, as flying qubits, they are the only viable choice for quantum networks [14] and quantum imaging [15,16].
Unfortunately, creating photonic resource states is fundamentally difficult. Because photons do not interact with each other, most attempts have focused on probabilistic generation schemes using linear optics and postselection [17], which are very resource-intensive, severely limiting the size of the resulting states [18,19]. This bottleneck can in principle be overcome by instead using a deterministic approach in which entangled photons are produced directly from quantum emitters (i.e., matter qubits). A general method for finding deterministic generation protocols has been known for some time [20,21], but it scales exponentially in the size of the target state and can lead to long generation circuits, motivating the search for more efficient approaches. Ref. [22] put forward one such scheme in the context of one-dimensional (1D) graph states produced from spin-tagged quantum dots; this was subsequently demonstrated experimentally [23]. Refs. [24,25] went on to propose protocols for 2D lattice graphs that leverage the principle that entangled emitters can emit entangled photons. This idea was extended further to develop protocols that deterministically generate resource states for quantum repeaters [26][27][28][29]- * economou@vt.edu † efbarnes@vt.edu tailored to color centers in Refs. [30,31]-and one-way computing [32,33]. Despite this progress and the intense interest this approach has generated among experimentalists, existing graph state generation protocols are limited to a small subset of graphs or require a number of emitters that scales linearly with the graph size [33,34]. This is extremely resource-intensive, especially in light of the schemes for generating repeater graph states presented in Refs. [26,28], which require only two emitters regardless of the number of photons. The required resources (number of emitters and entangling gates) is a critical factor that determines the practical feasibility of the protocol. For a general graph state, finding resource-efficient generation protocols in polynomial time remains an open problem.
Here, we address this challenge by presenting a general approach to generating arbitrary photonic graph states from quantum emitters. Given a target graph state, we show how to determine in polynomial time both the minimal number of emitters required to create it and an explicit generation protocol. The latter consists of a sequence of gate operations and measurements performed on the emitters. Moreover, our protocol naturally takes into account the order in which photons should be emitted, which can be an important consideration for applications, as it is generally preferable to emit photons in the order they are measured to avoid photon storage. Our method provides a recipe for doing this. The broad applicability of our method, its practical relevance, and its efficient use of resources make it ideally suited to the generation of any photonic graph state from various types of quantum emitters. highly nontrivial and markedly distinct from the problem of finding a quantum circuit that creates a target state on a register of qubits. Several additional challenges arise in the former, including the fact that qubits are both created and removed, and that different types of qubits (photons vs. emitters), with different roles and allowed gates, are involved. Depending on the experimental setup, there may also be further restrictions, e.g., emitted photons cannot interact with any other qubits following their emission (although schemes that re-interfere photons with emitters have been proposed [29,32]). Our method addresses these challenges by leveraging three main ingredients: the notion of the height function (which is related to the entanglement entropy), the stabilizer formalism, and the concept of timereversed emission events and measurements, which we introduce here.
The first insight is to utilize the so-called height function, which is the entanglement entropy of the system as a function of the partition point when the system is arranged in a 1D lattice and partitioned into two subsystems [35,36]. This function provides information about the entanglement structure of the target state as well as the number of emitters required to produce it. The latter is equal to the maximum value of the height function, as we show below. Moreover, the height function plays a crucial role in determining the sequence of operations (gates and measurements) needed to generate the target graph state from the emitters.
A second key ingredient is the use of gates from the Clifford group. Given that arbitrary graph states can be generated solely with Clifford gates [37,38], which were also exclusively used in the protocols of Refs. [22][23][24][25][26][27][28]33], restricting ourselves to this set does not affect the generality of our approach. Clifford gates enable the use of the stabilizer formalism, such that we can manipulate Pauli operators instead of keeping track of the whole state. This makes the problem of finding the emission operation sequence tractable, reducing it from exponential to polynomial scaling due to the Gottesman-Knill theorem [39].
A final key element in our algorithm is that we timereverse the emission sequence. That is, we start from a target multi-photon graph state and an appropriate number of decoupled emitters (obtained from the height function for the target state), and we determine a sequence of emitter gates, "time-reversed measurements", and "photon absorption" events such that the target state is converted to a product state. The final state is a product state because, without loss of generality, photons that have not yet been emitted can be described by qubits prepared in the computational basis state |0 . Photon emission is then modeled as a two-qubit photon-emitter gate that brings the photon from |0 into an entangled state with the emitters [22]. Because the photon absorption steps are time-reversed versions of photon emission, these too are described by photon-emitter gates.
The run time of the protocol solver algorithm scales as O(n 4 p ), where n p is the number of photons in the target graph state. This is a direct consequence of the fact that the algorithm is based on the stabilizer formalism (see Methods section). This is in contrast to previous methods [20,21], which scale exponentially in n p due to the need to perform singular value decompositions repeatedly. We also show that the number of gates in the final emission sequence scales at most as O(n 2 p ) (see Methods). However, this assumes two-qubit gates can be applied between any pair of emitters. If this is not the case, then additional SWAP operations are needed, bringing the gate count up to O(n 3 p ). Therefore, both the protocol solver and the resulting gate sequence it obtains scale polynomially in the size of the target graph state. Now we provide a more detailed description of the protocol solver algorithm. We begin with a target graph state |ψ p of n p photons and n e decoupled emitters, so that the total state is |Ψ = |ψ p ⊗ |0 ⊗ne . An n p = 4 photon example graph is shown in Fig. 1(a). This is what the state of the total system should be at the end of the generation sequence. n p is set by the size of the desired photonic graph state |ψ p , while n e remains to be determined. We assume the graph representing |ψ p is connected; if this is not the case, then the algorithm can be run separately for each connected subgraph. The state |Ψ is fully described by a set of n = n p + n e stabilizers g m , m = 1, . . . , n, defined such that g m |Ψ = |Ψ . The full set of n qubits can be arranged in a 1D lattice with site index x ∈ {0, 1, 2, . . . , n} (see Fig. 1(b)). Sites x = 1, . . . , n p correspond to the photons and are ordered according to the desired photon emission ordering, while the sites x = n p + 1, . . . , n are the emitters. The additional x = 0 site is included as a matter of convention. We can now define the height function h(x) = S A to be the bipartite entanglement entropy when the 1D lattice is divided into the subregion A = {1, 2, . . . , x} and its complement. Note that S A = 1 1−α log 2 Tr(ρ α A ) can be any of the Rényi entropies; for stabilizer states, they are all equal [40]. In Ref. [20], it was shown that the state of the emitted photons, |ψ p , can be represented by a matrix product state (MPS) with bond dimension 2 ne . Because the entanglement entropy of a MPS is given by the base-2 logarithm of the bond dimension [41], it follows that n e is equal to the maximum value of h(x). The height function for the graph in Fig. 1(a) is shown in Fig. 1(c). In this example, its maximum is 2, implying 2 emitters are needed. In general, the maximum of the height function is in fact the minimal number of emitters capable of generating the target graph state, as fewer emitters would be insufficient to match the bond dimension of any exact MPS representation.
The height function can be computed efficiently from the stabilizers. Because products of stabilizers are also stabilizers, there are many equivalent choices for the set {g m }. Here, we focus on a particular choice of the stabilizers that we refer to as the echelon gauge [42], in which the stabilizer matrix has a row-reduced echelon form (see Methods). When the g m are in this gauge, the height function can be expressed as [42] h where l(g m ) is the index of the left-most (smallest index) site on which g m acts nontrivially. The last term in Eq.
(1) counts the number of stabilizers that act nontrivially only on sites to the right of (i.e., larger than) x. Although Eq.
(1) depends on n e , this dependence cancels out for states like |Ψ in which the emitters are decoupled. Therefore we can obtain n e from the maximum of h(x) on the photonic sites, using only the stabilizers of |ψ p .
Once we have the number of emitters n e , we can run the protocol solver algorithm to determine the sequence of gates, time-reversed measurements, and photon absorption events needed to transform the target state |Ψ into the initial state |0 ⊗n , which corresponds to decoupled emitters and no photons. We first introduce a photon index j and initialize it to j = n p . The algorithm then consists of four steps: (i) Transform the stabilizers g m into echelon gauge if they are not already, then compute the height function h(x).
Otherwise apply a time-reversed measurement and update the g m accordingly.
(iii) Apply a photon absorption operation on the j-th photon and update the g m accordingly. If j > 1, then set j → j − 1 and go to step (i). Otherwise, go to step (iv).
(iv) All photons are now in state |0 . Apply a series of gates on the emitters to disentangle them, bringing the total state to |0 ⊗n .
This algorithm involves repeated applications of two basic operational primitives: time-reversed measurement and photon absorption. During the algorithm, the height function of the current state tells us which of these we need to perform next to bring the state closer to |0 ⊗n . Each photon absorption step disentangles one photon qubit from the rest, starting with the last-emitted photon, j = n p , and working down to the first photon, j = 1. For our 4-photon example, the graphs at intermediate steps of the algorithm are shown in Fig. 1(d). A step by step explanation of this example is given in the Supplementary Information. When the algorithm concludes, we can reverse the entire sequence to obtain an operation sequence that generates |ψ p starting from n e decoupled emitters. We now describe each of the two operational primitives in more detail, the precise gates they introduce into the generation sequence, and their connection to the height function.
Photon absorption of the j-th photon refers to a timereversed version of photon emission. For concreteness, we focus on the case where emission is described by a CNOT gate between the photon and its emitter (with the emitter as the control), as in Ref. [22], although our algorithm can be adapted to any Clifford gate describing photon emission. Mathematically, the task of absorbing photon j requires finding a stabilizer g a that can be transformed to σ z j by applying CNOT ij , where i is an emitter site. It is possible to find such a stabilizer when h(j) ≥ h(j − 1). From Eq. (1), we see that this condition implies there must be at least one stabilizer, g a , such that l(g a ) = j. This stabilizer has the form  where α, β k ∈ {x, y, z} label the nontrivial Pauli operators, and 1 ≤ j ≤ n p < i 1 < · · · < i s ≤ n. Note that we can assume g a acts trivially on all photons with index larger than j since these have already been decoupled at this point in the algorithm. We also assume that g a acts nontrivially on at least one emitter site; if this is not the case, then photon absorption is unnecessary since the jth photon is then already disconnected. To transform g a into σ z j , we can first apply a local Clifford operation on the j-th site and general Clifford operations on the emitters to transform g a → σ z j σ z i , where i > n p is an emitter site. This can be done for example by applying local Clifford operations to transform g a to σ z j σ z i1 · · · σ z is , and then applying CNOT gates on pairs of emitters to transform this to σ z j σ z i . Applying CNOT ij brings this to σ z j , completing the absorption of the j-th photon. Note that we can choose any emitter to absorb the photon; typically, the emitter that requires the shortest circuit to transform g a into σ z j is preferred. The resulting circuit is included in the time-reversed generation sequence.
Time-reversed measurements are applied whenever h(j) < h(j − 1), in which case photon absorption is not possible. Indeed, in this case, Eq. (1) implies #{g m |l(g m ) = j} = 0, or in other words, a suitable g a does not exist. In order to absorb the next photon, we must therefore first find a way to increase h(j) relative to h(j − 1). This can be accomplished with a timereversed measurement on an emitter. To perform this operation, we first rotate the state to |Φ ⊗ |0 i , where |Φ is a stabilizer state involving photons 1, . . . , j and emitters other than i. This can always be done using O(n e ) Clifford gates on emitters when h(j) < h(j − 1) (see Methods). Now notice that this state is obtained from the pre-measurement state CNOT ij |Φ ⊗|+ i when emitter i is measured to be in the state |0 i . Therefore, starting from |Φ ⊗ |0 i , if we perform a Hadamard gate on emitter i followed by the gate CNOT ij , we effectively reverse the measurement on the emitter. These operations transform the stabilizers g m in such a way that h(j) now satisfies h(j) ≥ h(j − 1) (see Methods), and we can proceed with the next photon absorption. The emitter gates, Hadamard on i, and CNOT ij are all included in the time-reversed generation sequence.

B. Examples
We demonstrate our algorithm with several examples. The first is the important case of repeater graph states [10], where we use our algorithm to obtain generation protocols that are more efficient than previously known ones. As a second example, we consider random graphs containing up to hundreds of photons and demonstrate the polynomial scaling of the resulting generation circuits. Additional examples, including modified repeater graph states, error correcting codes, and a simple example that illustrates the algorithm in detail can be found in the Supplementary Information.

Repeater graph states
Next, we apply our algorithm to find operation sequences that produce repeater graph states [10]. In addition to its importance in quantum network applications, this example also illustrates how different photon emission orderings impact the required number of emitters. Ref. [26] presented a generation protocol for a particular ordering that was devised essentially through guesswork. Our algorithm can be used to systematically find protocols for any ordering. An example of a 12-photon repeater graph state is shown in Fig. 2(a). The graph contains a fully connected core of 6 photons, each of which is connected to a single external photon. Bell measurements are performed on pairs of these external photons, where the two photons in each pair come from different graph states. If a Bell measurement succeeds, then the two corresponding core photons are linked by an edge, and entanglement extends across two nodes of the repeater network. Having multiple external photons provides built-in redundancy that increases the likelihood that at least one Bell measurement between two repeater graph states is successful. Upon success, core photons are then measured in the z or x basis to remove photons connected to failed measurements or to create entanglement links between successful measurements, respectively. Because the external photons are measured first, it may be advantageous to emit these first when generating the graph state to reduce photon storage requirements. This corresponds to the photon ordering shown in Fig. 2(a). The height function for this graph and photon ordering is shown in Fig. 2(d), where it is evident that 6 emitters are needed to produce the state. However, if efficient photon storage is available, then the ordering shown in Fig. 2(b) may be preferable, where now external and core photons are emitted in an alternating sequence. This ordering reduces the number of emitters down to only 2, as shown in Fig. 2(e). As we discuss further below, this illustrates our general finding that "natural" orderings in which neighboring vertices are emitted around the same time reduce the requisite number of emitters. This reduction in quantum resources becomes still more dramatic as the size of the graph increases; for orderings as shown in Fig. 2(a), the number of emitters scales linearly with photon number, while for the natural ordering of Fig. 2(b), the number of emitters remains at 2 regardless of the number of photons. This is shown explicitly in the Supplemental Information.
As discussed in Ref. [33], some of the edges in the repeater graph can be removed without affecting the functionality of the repeater. Fig. 2(c) shows an example of this in which 4 of the core edges are deleted. As shown in Fig. 2(f), the number of emitters is still 2. However, removing the redundant edges reduces the depth of the resulting generation circuit, which is shown in Fig. 2(g). This circuit contains 4 CNOTs between emitters and 1 intermediate measurement on an emitter, whereas the original protocol presented in Ref. [26] requires 5 twoqubit gates and 5 intermediate measurements.

Random graph states
To demonstrate how our algorithm scales with the number of photons in the target state, we run it for random graphs ranging in size from n p = 16 to n p = 256 photons. These graphs are produced randomly using the Erdös-Rényi model [43]. In this approach, each random graph is constructed by connecting n p vertices randomly with fixed probability p. We discard any graphs that contain disconnected vertices when sampling these realizations. The likelihood that such graphs arise becomes very small if p is chosen sufficiently close to 1. In Fig. 3, we show the maximum value, h max , of the height function averaged over 128 realizations for each value of n p . Averaged measurement and gate counts are also shown. It is evident that h max , and hence the number of emitters, scales linearly with n p as n p becomes large. The same is also true of the number of measurements. On the other hand, the number of CNOTs and the total number of gates in the resulting generation circuits scale quadratically with the number of photons in the target state. These results confirm both the polynomial scaling of our algorithm, which allows us to easily find generation protocols for graph states containing hundreds of photons, and the polynomial scaling of the resulting protocols, which makes them practical for near-term experiments. 16 32 64 128 256

III. DISCUSSION
A powerful feature of our algorithm is that it readily incorporates a desired photon emission ordering. This is encoded when we arrange the photons and emitters in a 1D lattice to define the height function. If no specific ordering is preferred, then ideally we would want to choose the ordering that minimizes the number of emitters n e . However, the task of finding this optimal ordering is NPcomplete, as we show in Methods. Nevertheless, one can still look for heuristic solutions to the problem. In fact, the expression for the height function in Eq. (1) makes it clear that this function is suppressed for orderings in which the stabilizers, when expressed in the echelon gauge, are supported predominantly on high-index sites on the right side of the 1D lattice. This tends to occur for "natural" orderings in which neighboring photons in the graph are emitted around the same time, because in this case the stabilizers are localized on the 1D lattice. This was illustrated with our repeater graph state example in the previous section. The extent to which the stabilizers can be localized in this way depends on the graph of course. For an N × M square lattice, it is inevitable that some neighboring vertices will be separated by M steps in the emission sequence (assuming M < N ), and so the number of emitters is of order M . On the other hand, for other graph structures like those of the repeater graph states, far fewer emitters may be needed, provided a natural photon ordering is used. Note that in this example, as for many graphs, edges between remote vertices cannot be avoided (see Fig. 2(b)). Despite this, we showed that optimal orderings for which the height function remains small can still be found. Thus, emitting neighboring vertices around the same time is sufficient but not always necessary to keep the number of emitters small.
In summary, we presented an efficient algorithm to construct polynomial-depth operation sequences that produce arbitrary multi-photon graph states from a minimal number of quantum emitters. By reducing both the number of photon sources and the number of quantum operations that need to be performed on them, our method brings the wide range of quantum information applications that rely on entangled photon resource states closer to experimental reality.

A. Echelon gauge
The echelon gauge was first defined in Ref. [42], where it was called row reduced echelon form. In this gauge, the stabilizer tableau has a recursive row-reduced form based on the following three types of matrices: where σ, σ 1 , and σ 2 are nontrivial Pauli matrices, and σ 1 = σ 2 . In this work, we always choose σ 2 = σ z , and σ 1 can be either σ x or σ y . The full tableau cannot have the first form shown above (with only identities in the first column), because this case does not apply to pure states. However, the submatrix M can follow any of the above three patterns, and the structure iterates recursively. The stabilizers can be transformed into this gauge starting from any other by performing a series of row reductions, as described in Ref. [42]. In the echelon gauge, the independent stabilizers acting on A = {x + 1, . . . , n} appear at the bottom right of the tableau. Therefore, starting from the formula for the entanglement entropy for subregionĀ of a stabilizer state [44], SĀ = nĀ − |GĀ|, where nĀ is the size ofĀ and |GĀ| is the number of independent stabilizers acting onĀ, and using h(x) = S A = SĀ, we obtain Eq. (1).

B. Time-reversed measurements
Above, we saw that when the total state of the system has the form |Φ ⊗ |0 i , where i is an emitter site, we can perform a time-reversed measurement to convert this to the pre-measurement state CNOT ij |Φ ⊗ |+ i . Here, we clarify two important questions regarding this process: (i) When and how can we bring the system into the state |Φ ⊗ |0 i ? (ii) How can we see that a time-reversed measurement on this state increases h(j), as needed for a subsequent photon absorption process?
Regarding question (i), when h(j) < h(j − 1), we can always find a set of Clifford gates that act purely on the emitters that will transform the state of the system into |Φ ⊗ |0 i . To see this, first note that h(j) = h(n p ), as follows from Eq. (1) when photons j + 1 through n p are in state |0 . Using that the height function is bounded from above by n e , we then have h(n p ) = h(j) < h(j − 1) ≤ n e . On the other hand, from Eq. (1) we have h(n p ) = n e − #{g m |l(g m ) > n p }. Together, these results imply #{g m |l(g m ) > n p } > 0, or in other words, there is at least one stabilizer that is supported solely on the emitter sites. We can therefore transform this stabilizer into σ z i using at most O(n e ) Clifford gates on the emitters, bringing the state to |Φ ⊗ |0 i . We can then convert this stabilizer to σ x i by applying a Hadamard gate on site i. This prepares the system for the second part of the time-reversed measurement process, which is the gate CNOT ij .
Proof We are assuming that h(j) < h(j − 1), which from Eq. (1) implies #{g m |l(g m ) = j} = 0. Now consider how the stabilizers transform under CNOT ij . If l(g m ) < j before the gate, then l(g m ) remains invariant, and the contributions of these stabilizers to h(x) remain the same after the gate. The only potential changes to h(x) come from stabilizers g m for which l(g m ) > j.
These stabilizers necessarily have 1 on the j-th site. Stabilizers among this set that have 1 or σ z i on the i-th site will be unchanged by the CNOT ij gate. However, if one or more of these stabilizers have σ x i or σ y i before the gate, then afterward, these stabilizers will contain σ x j . Consequently, h(j) increases, while h(j − 1) remains the same. In the echelon gauge, there can only be one stabilizer with σ x j as the left-most nontrivial Pauli. Therefore, h(j) → h(j) + 1 when CNOT ij is applied. Moreover, if the i-th qubit is stabilized by σ x i , then this becomes σ x j σ x i after the gate, and so the height function for all sites between j − 1 and i increases: h(x) → h(x) + 1 ∀x ∈ {j, j + 1, · · · , i − 1}.

C. Scaling analyses
Here, we determine the complexity of both the protocol solver algorithm itself and the resulting graph state generation circuit. Regarding the algorithm, the main factor that determines the complexity is the need to restore the stabilizers to the echelon gauge after each operation is applied. Transforming a n-qubit stabilizer state into the echelon gauge generally requires O(n 3 ) steps, which is the complexity of Gaussian elimination. Another important factor is the process of determining which gates need to be applied in preparation for photon absorption or time-reversed measurement. Solving for each set of gates takes no more than O(n e n) steps, which is the number of entries in the emitter part of the stabilizer tableau. Thus, the Gaussian eliminations needed to restore echelon gauge dominate the scaling. In the worst case where n e ∝ n, our algorithm will then take O(n 4 ) steps, where the additional factor of n comes from the fact that the algorithm requires O(n p ) ∼ O(n) iterations.
As for the complexity of the output generation circuit, there are at most O(n e ) operations between any two photon emissions. For example, O(n e ) gates are needed to transform g a into the appropriate form for photon absorption. Thus, the depth of the circuit acting on the emitter qubits is at most O(n p n e ). In the worst case where n e ∼ n p , the scaling is then O(n 2 p ), which is consistent with Fig. 3. Nevertheless, due to the fact that some long-range two-qubit gates may arise, and given that these are usually decomposed as O(n e ) short-ranged two-qubit gates in real devices, the overall circuit depth may become O(n p n 2 e ).

D. Complexity of finding optimal photon emission orderings
We can show that the task of finding optimal emission orderings is NP-complete by mapping this to a known graph theory problem. Define Γ ij to be the adjacency matrix of the graph representing the target state |Ψ . Ref. [40] showed that we can obtain the height function from Γ ij using the formula h(x) = rank 2 (Γ AĀ ), where Γ AĀ is the sub-matrix of Γ ij with row indices i ∈ A = {1, 2, · · · , x} and column indices j ∈Ā. Note that this expression does not simplify the computation of h(x); it can take more steps to find the maximum compared to using Eq. (1) since the former performs Gaussian eliminations for O(n p ) rounds, while the latter only takes one round. However, the maximum value of this alternative expression for h(x) is precisely equal to a graph theoretic property known as linear rank-width (LRW) [45]. The task of finding an optimal photon emission ordering is therefore equivalent to finding the graph isomorphism that minimizes the LRW, which has long been studied in coding theory in the context of optimizing block code trellises [46]. Unfortunately, determining whether a simple connected graph has an LRW bounded from above by a positive integer k (i.e., max x h(x) ≤ k) has been shown to be NP-complete [47,48]. Therefore, it is unlikely this problem can be solved efficiently for large, arbitrary photonic graph states unless P = NP. Nevertheless, if the parameter k is set to 1, this problem can be answered in polynomial time [49]. If the parameter k is set to larger values, a recent work [45] showed that this problem can be reduced to a fixed parameter tractable problem. Specifically, its answer, along with the sequence solution (if it exists), can be determined in O(f (k)n 3 p ) steps, where f (k) is an exponentially large function of k. However, the growth of f (k) is so rapid that this result is not likely to be of practical use for photonic graph state generation.   · · · · · · · · · · · ·