Photonic resource state generation from a minimal number of quantum emitters

Multi-photon entangled graph states are a fundamental resource in quantum communication networks, distributed quantum computing, and sensing. These states can in principle be created deterministically from quantum emitters such as optically active quantum dots or defects, atomic systems, or superconducting qubits. However, finding efficient schemes to produce such states has been a long-standing challenge. Here, we present an algorithm that, given a desired multi-photon graph state, determines the minimum number of quantum emitters and precise operation sequences that can produce it. The algorithm itself and the resulting operation sequence both scale polynomially in the size of the photonic graph state, allowing one to obtain efficient schemes to generate graph states containing hundreds or thousands of photons.


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). One possibility would be to prepare a graph state on emitters [20,21] and transduce it to photons, but this requires a number of emitters equal to the size of the target photonic graph state. This daunting resource overhead can be avoided by instead using sequential generation schemes. Refs. [22,23] put forward such an approach that works well for one-dimensional (1D) graph states [24] and has led to experimental demonstrations [25,26]. However, in the general case where the entanglement structure is more complicated, this method scales exponentially in the size of the target state and can lead to long generation circuits, motivating the search for more efficient approaches. Refs. [27,28] * libk@vt.edu † economou@vt.edu ‡ efbarnes@vt.edu put forward 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 [29][30][31][32]-tailored to color centers in Refs. [33,34]-and one-way computing [35,36].
Refs. [32,35] allowed for the re-interference of photons with emitters to further enhance flexibility in entanglement creation.
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 [36,37]. This is extremely resource-intensive, especially in light of the schemes for generating repeater graph states presented in Refs. [29,31], 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.

II. RESULTS AND DISCUSSION
A. Overview of the algorithm Determining how to efficiently generate an arbitrary photonic graph state from a set of quantum emitters is highly nontrivial and markedly distinct from the problem of finding an efficient quantum circuit that creates a target state on a register of qubits [38]. 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 [32,35]). 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 [39,40]. 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 (see below), which depends on the photon emission order. Optimizing this order is NP-hard in general, although we show that heuristic approaches exist for more structured graphs. 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 [41,42], which were also exclusively used in the protocols of Refs. [24,25,[27][28][29][30][31]36], 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 [43].
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. This is somewhat reminiscent of disentangling circuits used for quantum state tomography of 1D systems [44]. 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 [24]. 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 [22,23], 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 [45]. In Ref. [22], 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 [46], 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 [47], 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 [47] 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. [24], 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 exam-ple that illustrates the algorithm in detail can be found in Supplementary Notes 1-4.
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. [29] 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. [36], 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. [29] requires 5 two-qubit gates and 5 intermediate measurements.
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 [48]. 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 1024 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.

C. Photon emission ordering
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 or-dering 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 NPhard, 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. [47], 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. [47]. 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 [49], 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) re-main 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 twophoton 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-hard by mapping this to a known graph theory problem. Define Γ ij to be the adjacency matrix of the graph representing the target state |Ψ .
Ref. [45] 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 optimized maximum value of this alternative expression for h(x) (i.e., max x h(x)) is precisely equal to a graph theoretic property known as linear rank-width (LRW) [50]. The task of finding an optimal photon emission ordering is therefore equivalent to finding LRW through the graph isomorphism, which has long been studied in coding theory in the context of optimizing block code trellises [51]. Unfortunately, determining whether a simple connected graph has an LRW bounded from above by a positive integer k has been shown to be NP-complete [52,54]. 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 [53]. If the parameter k is set to larger values, a recent work [54] 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.

DATA AVAILABILITY
The data that support the findings of this study are available from the authors upon request.

CODE AVAILABILITY
A custom MATLAB code to reproduce our results is available on GitHub and archived in Zenodo (https://doi.org/10.5281/zenodo.5652105).

COMPETING INTERESTS
The authors declare that there are no competing interests.

SUPPLEMENTARY NOTE 1
This supplementary material contains several additional examples of generation protocols produced using our algorithm. We begin with a simple example that illustrates in detail how our algorithm works. We then provide solutions for more complicated examples of practical importance, including error correcting codes and repeater graph states of arbitrary size. These more complicated examples are solved numerically using MATLAB codes that are available on GitHub [S1]. This code expresses the generation sequence in terms of an MPS [S2] for bookkeeping purposes: in which the initial and final states of emitter qubits are simply product states of |0 : |ψ f = |ψ 0 = |0 ⊗ne . We denote η j as the emitter qubit that emits the j-th photon and µ j as the emitter qubit that is measured after emitting the j-th photon.Ê ηj is the emission tensor,Ê ηj = |0 j |0 ηj 0| ηj + |1 j |1 ηj 1| ηj , that describes emission of photon j from emitter η j , which can be represented as CNOT ηj ,j . U e,j is the unitary operation obtained from the j-th photon absorption step, which transforms g a as explained in the main text (Sec.II.A).M j is identity if no measurement happens (µ j is not assigned), otherwise,M j = W j H µj X sj µjπµ j , with projectionπ µj ≡ 1 2 [1 + (−1) sj Z µj ] and its random outcome s j ∈ {0, 1}. Here, W j is the unitary operation that is obtained from time-reversed measurement (Sec.IV.B), and W 0 is the unitary operation that disentangles all emitters at the final stage of the time-reversed sequence. Finally, U p,tot = j X sj j U p,j is the local Clifford operation that acts on photons with conditional X sj j flipping. The profile of the solution is stored as {U e,j , U p,j , µ j , η j , W j , W 0 }. We note that such a solution is usually not unique due to there being multiple choices for how to choose the emitter gates and emitter sites in each photon absorption and time-reversed measurement.
As discussed in the main text, the height function plays a central role in determining the number of emitters and the operation sequence needed to generate a target photonic graph state. As shown in Eq. (1) of the main text, when the stabilizers g m are in the echelon gauge, the height function can be expressed as where l(g m ) is the index of the left-most (smallest index) site on which g m acts nontrivially. In the main text, we showed that the difference in the height function across adjacent sites determines whether we perform a photon absorption or a time-reversed measurement at each step of the algorithm. Therefore, we define from which it is apparent that this difference only depends on the number of stabilizers (in the echelon gauge) that have a left-ending on site x. We begin by demonstrating our protocol solver algorithm in the case of the simple 4-photon graph state displayed in Supplementary Fig. S1(a). The stabilizers are given by We can switch to the echelon gauge by redefining g 3 → g 2 g 3 . We then calculate the height function using Supplementary Eq. (S2), finding that the maximum is 2. Therefore, at least n e = 2 emitter qubits are needed, and so we assemble a 6-qubit lattice. We can depict the complete set of 6 stabilizers as a tableau, as shown in Supplementary  Fig. S1(b). In Supplementary Fig. S1(c), we first obtain inset (1) by transforming Supplementary Fig. S1(b) to the echelon gauge. The upper left sub-block of the tableau is exactly Supplementary Eq. (S4) with g 3 → g 2 g 3 . Next we describe in detail how the generator set is updated from inset (1) to inset (17) step by step. The column label j indicates which photon we are currently focusing on, and the labels (i),...,(iv) indicate the specific step of our algorithm. For each photon, we do the following steps: • j = 4: (i) Obtain inset (1) by transforming to echelon gauge: g 3 → g 2 g 3 . (ii) Supplementary Eq. (S3) gives δh(4) = −1. Perform a time-reversed measurement on emitter site 5 by applying a Hadamard H 5 followed  S1. Step-by-step illustration of the protocol solver. Perform a time-reversed measurement on emitter site 6 by applying H 6 followed by CNOT 63 . Inset (6) is then obtained by redefining g 5 ↔ g 6 . (iii) Let g a = g 5 = σ x 3 σ x 6 in inset (6). One gets inset (7) by applying Hadamards on sites 3 and 6. Then the 3rd photon is absorbed by applying CNOT 63 . Replace g 3 → g 3 g 5 to eliminate the redundant σ z 3 . Thus, inset (7) becomes (8). (µ j = 6, η 3 = 6, U p,3 = H 3 , U e,3 = H 6 , W 3 = H 6 .) • j = 2: (i) Skip this step since inset (8) is already in echelon gauge. (ii) Skip this step since Supplementary Eq. (S3) gives δh(2) = 1. (iii) Choose g a = g 4 = σ z 2 σ z 5 σ x 6 in inset (10). One gets inset (11) by applying H 6 followed by CNOT 65 on the emitters, so that g a → σ z 2 σ z 5 . Then the 2nd photon is absorbed into emitter 5 by applying CNOT 52 . Redefine g k → g k g 4 for k = 1, 3 to eliminate the redundant σ z 's. Thus, inset (11) becomes (12). (η 2 = 5, U p,2 = 1, U e,2 = H 6 CNOT 65 ,M 2 = 1.) • j = 1: (i) Obtain inset (13) from (12) by transforming to echelon gauge: (g 3 , g 4 , g 5 , g 6 ) → (g 4 , g 5 , g 6 , g 3 ). (ii) Skip this step since Supplementary Eq. (S3) gives δh(1) = 1. (iii) Choose g a = g 2 = σ z 1 σ x 5 σ z 6 in inset (14). One gets inset (15) by applying H 5 and then CNOT 65 to transform g a → σ z 1 σ z 5 . Then the 1st photon is absorbed into emitter site 5 by applying CNOT 51 . Thus, inset (15) becomes (16). (η 1 = 5, U p,1 = 1, U e,1 = H 5 CNOT 65 , M 1 = 1.) • (iv) Finally, to recover the state |0 ⊗n , one needs to disentangle the emitter qubits. This can be done with the following gate sequence: H 5 CNOT 56 H 5 . In the last step, we permute the g m to obtain inset (17). (W 0 = H 6 CNOT 56 H 5 .) Now that the algorithm is complete, we reverse all the operations to obtain the final generation sequence. This circuit is shown in Supplementary Fig. S2(a). It is worth noting that, in this example, the emission sequence can be further optimized by swapping the 1st and 3rd photons in the emission order, such that the maximum of h(x) is reduced to 1. Thus, only one emitter qubit is needed in this case, and the corresponding generation circuit is displayed in Supplementary Fig. S2(b).

SUPPLEMENTARY NOTE 2
In this subsection, we demonstrate how to generate a useful quantum error correction code, with some continuous logical rotation. In particular, we present an emission sequence for the Shor code [S3] with 9 photonic qubits, which is able to protect a qubit from single bit-flip and phase-flip errors. The stabilizer generators of this code are well