Scheme for Universal High-Dimensional Quantum Computation with Linear Optics

Photons are natural carriers of high-dimensional quantum information, and, in principle, can benefit from higher quantum information capacity and noise resilience. However, schemes to generate the resources required for high-dimensional quantum computing have so far been lacking in linear optics. Here, we show how to generate GHZ states in arbitrary dimensions and numbers of photons using linear optical circuits described by Fourier transform matrices. Combining our results with recent schemes for qudit Bell measurements, we show that universal linear optical quantum computing can be performed in arbitrary dimensions.

Photons are natural carriers of high-dimensional quantum information, and, in principle, can benefit from higher quantum information capacity and noise-resilience.However, schemes to generate the resources required for high-dimensional quantum computing have so far been lacking in linear optics.Here, we show how to generate GHZ states in arbitrary dimensions and numbers of photons using linear optical circuits described by Fourier transform matrices.Combining our results with recent schemes for qudit Bell measurements, we show that universal linear optical quantum computing can be performed in arbitrary dimensions.
Photonics is a sophisticated platform for the development of quantum technologies, from quantum processors to distributed quantum communication [1][2][3][4][5].Until now, linear optical architectures have focused on encoding photons as qubits (two-level systems).Yet using higher dimensional systems -qudits -can in principle improve the information capacity and noise tolerance of computational resources, and potentially unlock new routes to fault-tolerant quantum computing and distributed quantum networks [6].Qudits can naturally be encoded in photons using d orthogonal optical modes in a variety of degrees of freedom, e.g.spatial modes [7][8][9], orbital angular-momentum [10][11][12][13][14], optical frequencies [15,16], and time-bins [17,18].High-precision control and arbitrary operations of single qudits have been demonstrated using programmable interferometers [7,8,19,20].However, architectures for universal quantum photonic processors based on higher-dimensional systems have, so far, been absent.While there has been significant experimental progress in photonic qudit entanglement generation [7][8][9][10][11][12][13][14][15][16]18], the post-selected schemes used so far can only generate a limited set of high-dimensional entangled states and present no clear route to scalability [21,22].In fact, in contrast to the qubit case [23,24], even determining which high-dimensional entangled states can be generated with single photons and linear optics has so far remained an open problem [6].
Here, we answer this question by showing that all highdimensional entangled states with fixed numbers of photons can be generated with linear optics and, in principle, with a scalable architecture.That is, we show that universal linear optical quantum computing (LOQC) is possible in arbitrary dimensions.Key to this result are linear optical schemes for the generation of heralded N -photon GHZ entanglement in arbitrary dimension d.Combining these schemes with previous results on Bell measurements with photonic qudits to fuse d-dimensional GHZ states, and techniques from qubit LOQC architectures, we obtain a scheme for universal measurementbased quantum computing with photonic qudits.
Heralded high-dimensional GHZ entanglement.The GHZ states for N photons, each encoding a qudit of dimension d, are defined as and represent the initial resource states in an architecture that builds cluster states of qudits for universal quantum computing.Our scheme to generate these states, shown schematically in Fig. 1a, consists of single photon sources, linear optical elements, and photon-number resolving detectors.Our scheme exploits the zero-transmission law (ZTL) in discrete Fourier transform (DFT) interferometers (previously investigated e.g. for the verification of boson sampling protocols [19,25,26]).In particular, we use a result of Tichy et al. [27] which states that, if m single photons are sent in the m individual input modes of an m-mode DFT, described by the unitary U j,k = exp jk 2πi m / √ m, then the output configurations c with non-zero amplitude must satisfy where c i represents each particle's output port, with modes indexed from 0 to m − 1.For example, c = (0, 1, 1, 4) represents the configuration with one photon in each of modes 0 and 4, two photons in mode 1, and zero photons in all remaining modes.The ZTL condition is valid for any value of m.
To describe how the ZTL can be used to generate GHZ states, we introduce a notation for the encoding of individual qudits in the output modes of the DFT.Each qudit of dimension d is encoded in a single photon propagating through d optical modes (see e.g.Refs.[7,8]).We will denote as X ⊂ {0, . . ., m − 1} the total set of N d modes used to encode the N qudits, while its complement X = {0, . . ., m − 1} \ X represents all the auxiliary modes used for heralding the desired state.Each Optimised solution for the heralded generation of Bell states (N = 2) in arbitrary dimension using the DFT scheme, requiring m = 2d + 1.The correspondence between the optimal modes and the associated computational states of each qudit is shown.c.Optimised solution for the heralded generation of 3-photon 3-dimensional GHZ states using the DFT scheme, requiring m = 25.Output modes associated to the three different qudits are highlighted with different colours, and the associated qudit computational states are labelled.
of the N qudits is defined via a set of d distinct modes so that all the qudits {Q i } form a partition of X .We say that the i-th qudit is in the logical state |k i if a single photon is present in mode x (i) k , while all other modes in the set Q i are vacuum.
For simplicity, we will focus on heralding configurations where all the m − N auxiliary photons are detected in the zero-th mode, and vacuum in the remaining heralding modes X \ {0}, as shown in Fig. 1a.In this way, all auxiliary photons give zero contributions to the sum in Eq. 2, and only contributions from the N encoding photons remain.Under this condition, we can write the set of all possible remaining N -photon output states induced by the heralding on X and the ZTL as Here b k are the output configurations allowed by the ZTL, indexed by k, and y k,i represent each photon's output port for that configuration, as in Eq. 2.
We will now show that the following conditions for the set B X are sufficient (but not necessary) to obtain a GHZ state at the output: The first step is to show that, if these three conditions are satisfied, we can specify N well-defined photonic qudits, given by Q i = {x (i) k = y k,i } k=0,...,d−1 .Note that, because the sets {Q i } are given by a simple transposition of indices of B X , condition 2 immediately implies that {Q i } also forms a partition of X .Moreover, conditions 1 and 2 imply that all y k,i are different, i.e. y k,i = y k ,i if k = k or i = i .In fact, because B X contains exactly d configurations of N elements, if any two y k,i were identical then necessarily k b k < N d, and thus B X could not form a partition of X , given that |X | = N d.Finally, from their definition and the fact that the elements y k,i are all different, it immediately follows that each set Q i contains exactly one output photon.Our specification for the sets {Q i } therefore defines N disjoint sets of d different modes, each set containing exactly one photon, and thus provides a valid encoding for the N qudits.We are now left to verify that the state of the N qudits is in fact a GHZ state.
Note that, with the definition used for the Q i s, the d elements b k of B X correspond to the logical N -qudit states b k → |k, k, . . ., k , k ∈ 0, . . ., d − 1. Due to the ZTL, the total output state is therefore a superposition of the states |k, k, . . ., k .As shown in Appendix 1, condition 3 ensures that all amplitudes in such superposition are uniform and non-zero, thus providing the desired N -photon GHZ state in d dimensions.
General and optimised schemes.The task of generating GHZ states for a given photon number N and dimension d can now be reduced to a combinatorial number theory problem: finding an integer number m and a set X ⊆ {1, 2, . . ., m − 1} (0 is occupied for the heralding) so that B X satisfies conditions 1-3.Solutions can be found for any N and d.An example of such general solutions is given by the set with m = (N N d −1)/(N d −1), with details in Appendix 2. While this particular solution is highly sub-optimal in the number of resources and success probability, and likely not suitable for practical implementations, it is general and shows that arbitrary N -photon d-dimensional GHZ states can, in principle, be generated with linear optics.
For small values of N and d, optimised solutions for our heralded GHZ generation scheme can be found numerically.For example, in Fig. 1c we report the optimised solution for the case of heralded GHZ generation for (N, d) = (3, 3), which requires m = 25.The list of modes used for the encoding in this case is X = {1, 2, 3, 4, 5, 9, 13, 16, 22}, and the only three-element combinations (including repetitions) that sum up to multiples of m, given by B X = {(1, 2, 22), (3,9,13), (4,5,16)}, satisfy conditions 1-3.This therefore provides heralded GHZ generation for the three qudits defined in the modes Note that, at the output of the DFT, adjacent modes associated to different qudits can be interleaved; in such cases, a network of swaps is required to separate the different qudits and address them individually.However, without further optimisation, the three-dimensional GHZ generation success probability remains very small, approximately 10 −10 .
The main reason for low success probabilities is the use of a single heralding pattern amongst exponentially many possible outcomes.While the choice of a single heralding pattern was done to simplify the treatment in the general case of arbitrary N and d, for a given N and d many more heralding patterns are likely to generate GHZ entanglement.We show in Appendix 3 that these (combinatorially many) valid heralding patterns can be used, in conjunction with feed-forward operations and balancing circuitry, to significantly improve the success probability.For example, Monte Carlo simulations of the (N, d) = (3, 3) GHZ generation scheme of Fig. 1c show that the success probability is boosted to approximately 10 −4 .This indicates that many orders of magnitude improvements can be found through solution-specific optimisations.We estimate that, with this optimisation, heralded (3, 3) GHZ generation could be achieved at 0.7 kHz rates with state-of-the-art quantum photonic hardware.Furthermore, in Appendix 2 we report an algorithm to estimate optimized solutions for larger values of N and d.
Constructing universal cluster states of qudits.Measurement-based quantum computing (MBQC) [28] in linear optics typically proceeds by connecting small entangled resource states using probabilistic fusion gates [24], to build large cluster states.Similarly, high-dimensional photonic GHZ states can be used as building-blocks to construct large high-dimensional cluster states for MBQC.Two recent protocols for type-II fusion (destructive Bell state measurements) of arbitrary dimensional qudits have been proposed independently: Luo et al. [13]  As shown in Fig. 2, we can combine type-II highdimensional fusion operations with three-photon highdimensional GHZ generators in a modular approach to build up universal cluster states of qudits.The square lattice of qudits shown in Fig. 2b is an example of a cluster state universal for high-dimensional MBQC [30,31].Fig. 2a shows a module for the architecture, where multiple GHZ states are fused together to link a single computational photonic qudit to the rest of the lattice.Once such states are built, universal measurement-based high-dimensional quantum computing can be performed [28,30,31]: operations on the logical qudits encoded in the rows of the lattice are performed via measurement and feed-forward, with the output state encoded in the qudits of the last layer of the lattice.Because the resource is universal, any pure highdimensional multi-photon state can be prepared as a result of the computation up to arbitrary precision.This implies that the generation of arbitrary quantum states comprised of d-dimensional photons, for any d, is possible using only linear optics.
On the other hand, because both the GHZ generation and the type-II fusion gates have low success probability, the total probability to successfully generate highdimensional qudit cluster states can in general vanish quite rapidly when increasing the number of qudits or the dimensionality.Nevertheless, by adapting techniques already developed for qubit-based loss-tolerant LOQC architectures, the approach can in principle be made scal- able and the total success probability boosted to nearunity [24,[32][33][34].For example, because the GHZ generation is heralded, gate multiplexing can be used to render the production of GHZ states near-deterministic with a resource overhead that scales approximately linearly with the generation success probability [33,35].Repeat-untilsuccess proposals provide a flexible approach to correct for the limited success probability of the fusion gates, at the cost of requiring quantum memories [24].Moreover, if the fusion success probability is improved above the percolation threshold of the lattice used, ballistic architectures can be used to correct the unsuccessful fusion gates directly on the generated lattice without the use of quantum memories [32,33].However, this approach would likely require the use of lattices with increased valency to bring the percolation threshold above the current qudit fusion gates success rates [36].In addition, it may be possible to increase the success probability of the qudit fusion gates by using additional ancilla resources, as has been shown for qubits [37,38].

Robustness against photon distinguishability.
While requiring additional resources compared to qubitbased approaches, qudits can provide improved robustness to noise, with potential benefits for quantum communication and fault-tolerant quantum technologies [39][40][41].An important source of noise in quantum photonics is distinguishability between optical modes, which can arise from imperfections in solid-state or spontaneous photon emitters.Here, we numerically analyse how distinguishability affects the multi-photon high-dimensional entanglement generated using the schemes proposed in ation is simulated using an approach by Tichy describing multi-photon interference of partially distinguishable photons [42] (see Appendix 7 for details).Because the number of photons involved in the scheme increases rapidly with N and d, this simulation quickly becomes intractable when increasing the heralded state complexity.Nevertheless, for the scheme shown in Fig. 1b with N = 2 and low values of d, the simulation remains tractable on a standard laptop.
For various levels of distinguishability, we show in Fig. 3a the fidelity of the generated state to the ideal qudit Bell state, and in Fig. 3b the logarithmic negativity as a figure of merit to assess the generated entanglement.The logarithmic negativity is an entanglement monotone that upper bounds the distillable entanglement, which in turn quantifies the amount of pure state entanglement that can be extracted under local operations and classical communication [43,44].Distinguishability weakens the interference governing the ZTL, and therefore, as expected, both the state fidelity and generated entanglement decrease when the indistinguishability is reduced.However, while the fidelity is always lower for higher dimensional states, the negativity reaches zero at larger values of photon distinguishability for higher dimensions, as highlighted in the inset of Fig. 3b.This indicates that for increased dimensionality, the generated entanglement can endure higher levels of photon distinguisha-bility, even if the number of imperfect input photons is larger.
Discussion.The large (but polynomial) resource overheads of high-dimensional LOQC could be compensated to some extent by quantum error correction protocols that are more efficient due to robustness to noise that increases with qudit dimension.Furthermore, and similarly to the improvements made in efficiency for qubit-based LOQC [23,24,33], we expect our results to provide a first and important step for developing high-dimensional LOQC architectures closer to mid-term technological capabilities.For example, as detailed in Appendix 3, by considering all valid heralding events in our GHZ generation scheme for the particular case with d = 3 (see Fig. 1c), the resource overheads for multiplexing are reduced from O(10 10 ) to O(10 4 ), showing that enormous improvements are possible via solution-specific optimisation.Crucially, this work proves that such solutions can exist for any dimension.Here we provide the technical details for the derivation of the output amplitudes in the DFT-based N -photon d-dimensional heralded GHZ generator described in the main text.We start by reporting some useful formulas for the analytical calculation of bosonic evolutions in DFT interferometers, and then use these results to obtain the output amplitudes and success probabilities of the heralded GHZ generator circuit.b.Formulas Formula 1. Zero-Transmission Law (ZTL) [27]: Given a vector a = (a 0 , a 1 , . . ., a m−1 ), with Proof.Derived in Ref. [27].
Formula 2. Given a vector a = (a 0 , a 1 , . . ., a n−1 ), with a i ∈ N 0 , and n ≤ m: Proof.This corollary of the ZTL can be simply obtained by padding m − n zeros to the vector a, i.e. defining the m-element vector b = (a 0 , . . ., a n−1 , 0, . . ., 0).Obviously we have a i = 0 mod m, and We can then simply apply Formula 1 to obtain Formula 3. Given a vector a = (a 0 , a 1 , . . ., a n−1 ), with a i ∈ N 0 , such that n−1 i=0 a i = m and a i < m ∀i, we have Proof.Because n−1 i=0 a i = m and a i < m, then at least two elements in the vector are non-zero.Therefore, given also the symmetry under permutations of a of the sum in the formula, we can consider a 0 > 0 and a n−1 > 0 without loss of generality.Defining we can rewrite T by fixing the n-th value of the permutation to be b ∈ {0, . . ., m − 1} and summing over the permutations of the remaining m − 1 elements and over all possible values of b: where For all n < n, we can expand S n (b) in a convenient iterative form.The idea is to consider the sum for all permutations in P m excluding the permutations where b actually appears in the first n elements.Given that only the first n elements are useful, the contributions in this sum are exactly the same as when summing over P m \ b, but each one repeated m − n time (the possible ways to distribute b in the remaining m − n elements).Considering also a symmetry over permutations of (a 0 , a 1 , . . ., a n−1 ), we can finally write Because n−1 i=0 a i = m and a 0 > 0 and a n−1 > 0, for all 1 < n < n we have which allows us to use Formula 2 to set the first term in the expansion of S n (b) to zero: The term S 0 (b) can be easily calculated as where we used that, because 0 < a 0 /m < 1, necessarily 1 − exp (2πia 0 /m) = 0. Substituting in Eq.S9 we obtain a simple recursive relation which provides Substituting in Eq.S4 we finally obtain which is interestingly independent on the vector a.

Output amplitudes and success probabilities in the DFT-based GHZ generator
An n-photon output configuration can be described by an n-element vector x = (x 0 , x 1 , . . ., x n−1 ), with x i ∈ 0, . . ., m − 1, representing the n output modes of the photons.Each configuration can be associated to an m-element occupancy vector s( x) = (s 0 , s 1 , . . ., s m−1 ), with each element s i ∈ 0, . . ., n indicating the number of photons in the i-th optical mode.If m photons are injected into the m distinct input modes of a m × m DFT, defined via the unitary U j,k = exp jk 2πi m / √ m, the amplitude of a given m-photon configuration x at the output is given by the matrix permanent [27] Φ( x)|Ψ = 1 In the main text we show that in the N -photon d-dimensional heralded GHZ state generator shown in Fig. 1 Now, because m i=1 y k,i = m and y k,i < m for all k and i, we can estimate this sum over the permutations using Formula 3, giving the amplitude It is immediate to see that this quantity does not depend on k, which immediately implies that for any N and d the amplitudes are equal for all k, and therefore that the state is in fact an N -photon d-dimensional GHZ state.The success probability is given by where we have now written explicitly the dependence of m on both d and N .Although, as discussed in Appendix 2, for each d and N a value of m that provides the desired GHZ structure exists, it is difficult to find an analytic dependence of m from d and N in the general case.However, because the output circuit needs to encode N photons in dimension d, necessarily m > N d, meaning that P succ (d, N ) in our scheme necessarily decreases exponentially with both N and d.Such unfavourable scaling is due to the fact that, for simplicity in the calculations to show the generality, we are considering only one out of combinatorially many possible heralding combinations that can give raise to GHZ states.
Approaches to explore all such combinations to exponentially boost the success probability are discussed in section 3.
For the case of Bell states, i.e.N = 2, a simple solution with m = 2d + 1 is shown in the main text, implying a success probability: (S19) Appendix 2: Examples of solutions for GHZ generation We here investigate solutions to the scheme to generate GHZ of N photons in arbitrary dimension d via DFT interferometers.Note first that, given a set X , verifying whether it allows a partition B X that satisfies conditions 1-3 is a special case of the k-partitioning problem, and therefore is in general NP-hard [45].Nevertheless, as discussed below, general solutions for arbitrary N and d can be found.These, however, are likely to be highly not optimal in terms of auxiliary resources required and success probability.We also report more efficient solution for cases where N and d are small enough to make the optimisation task computationally feasible.

A general solution for GHZ generation
We describe here in more detail the example provided in the main text of a general solution for m and X that satisfies conditions 1-3 from the main text.The scope of this example is to show that our scheme based on linear optics and Fourier-transform interferometers provides a way to generate of N -photon d-dimensional GHZ states for any N and d, but for simplicity we do not focus here on the efficiency of the protocol.In fact, solutions from this general example are very sub-optimal compared to the special case shown in the next section.We leave the improvement of success probabilities and the reduction of resource requirements for general solutions as an interesting open question.
To find a solution for the generation of a GHZ state of N photons and d dimensions with our scheme, we need to find an integer m and and set X ⊆ {0, . . ., m − 1} such that the set satisfies conditions 1-3 from the main text.
In order to build up a general solution, it is useful to note that, if we consider a solution X and remove all the terms y k,N , then the sums must all have different values.A simple, although possibly inefficient, way to ensure this is to start with a reduced set Y = X \ {y k,N } k such that any combinations of N − 1 elements from this set, including repetitions, provide different sums.The task of finding Y is then similar to the problem of finding sets of integers with distinct subset sums, which is a well-known problem in combinatorial number theory (note however that multisets need to be considered to include sums with repeated elements) [46][47][48].A simple solution for a set such that any combinations of N − 1 elements from this set (including repetitions) provide different sums is given by Y = {N 0 , N 1 , . . ., N (N −1)d−1 }.In fact, different combinations correspond to different integer numbers in basis N , therefore providing different sums.We can then complete Y into X by setting y k,i = N id+k for all i ∈ [0, N − 1] and k ∈ [0, d − 1], and defining the last d elements as where This ensures all b k sum up to m.We are now left to choose m.A challenge here is to pick m so that no N -elements combinations of X sum up to multiples of pm with p > 1, and that no y k,N overlaps with the elements already in Y.A way to do this is to pick m > N max k S k , therefore we use In conclusion, the resulting general solution is given by the set and m = (N N d − 1)/(N d − 1).The scope of this solution is to illustrate that for any N and d the generation of N -photon d-dimensional heralded GHZ states is possible with only linear optics using our scheme.However, we remark that this solution is far from optimal.Below we describe algorithms able to numerically find significantly better solutions for computationally practicable values of N and d.

Brute-force optimised solutions for small values of N and d
For small values of N and d, optimal solutions for the heralded GHZ generation with our scheme can be found numerically using brute force approaches.The parameter to be optimised is the number of photons and modes m, which is also connected to the success probability via Eq.S18.The brute force search for solutions with minimal m was performed as follows: we start from a value of m = m 0 large enough (we started from m 0 = 40 for the search of the solution in Fig. 1c); we then test the possible subsets X ⊂ {0, . . ., m} with N d, and check whether they satisfy conditions 1-3, i.e. if they represent valid solutions to the GHZ generation scheme.Because the number of such possible subsets increases combinatorially as m N d , if m is too large for a reasonable computational time (in our case m 30), we used a limited number of randomly sampled subsets rather than all possible subsets.We empirically found that, for the values of m, N and d tested here, testing 10000 different random subsets was enough to find possible solutions with high probability, if solutions existed.If at least one solution is found for a given m, we decrease it by one (m → m − 1) and repeat the process.If we don't find any solution for a given m, we conclude that the optimised value is m = m + 1, and a solution X is also given by one of those found at the step with m = m.Using this brute force approach, in the case (N, d) = (3, 3) we found the optimised solution with m = 25 described in the main text and shown in Fig. 1c.

3.
An algorithm for near-optimal solutions with larger N and d While the brute force method described above always provides optimised solutions for the proposed GHZ scheme (in terms of number of input photons m), it is very computationally costly and we found it not practicable to find solutions for larger values of N and d.On the other hand, the general solutions from Eq. S23 readily provide analytical solutions for any N and d, although highly sub-optimal and likely not suitable for experimental implementations.An approach which achieves a trade-off between these two methods, i.e. providing near-optimal solutions but with significantly better computational cost, is therefore worth investigating.
We developed an algorithm to find solutions for the GHZ scheme which performs such a trade-off, and is readily able to find solutions for values N ≥ 3 and d ≥ 3 on a standard laptop, significantly improving the general analytical solutions of Eq.S23 in terms of resource optimisation and the brute-force approach in terms of computational cost.The pseudocode of the algorithm is provided in Algorithm 1, and a Python implementation is freely accessible on GitHub [49].The algorithm, which we call Different-Sums Iterative (DSI) algorithm, proceeds in two parts.In the first part we iteratively build up a set Y of (N − 1)d distinct positive integers such that any sum of N − 1 of its elements, possibly including repetitions, is different and such that max(Y) is minimal.This part is thus similar to the first part of the analytical approach, where Y was built analytically at the cost of large overheads in terms of max(Y).This optimised solution for Y is found iteratively: starting from the (N − 1)-element set Y = Y 0 = {1, . . ., N − 1}, at each step we find the smallest positive integer ā such that, if it is included in the current set Y, it preserves the desired condition (i.e. that any sum of N − 1 of elements, including repetitions, is different).The number ā is then added to Y and the process is repeated (N − 1)(d − 1) times to get a set Y with (N − 1)d distinct elements.This part is fast enough to find sets Y with N and d 8 on standard laptops.
To complete the set Y to a possible solution of N d qudit-encoding modes for the GHZ scheme, we can proceed similarly to the second part of the analytical approach (Eq.S21): because each partition P = {P 1 , . . ., P d } of Y into d sets of N − 1 elements (|P i | = N − 1) will contain subsets P i with different sums σ i = p∈Pi p = σ j , choosing an additional element ȳi = m − σ i we get that each subset P i ∪ {ȳ i } will now sum up to m.We thus obtain the set X = Y ∪ {ȳ 1 , . . ., ȳd } as a possible solution.The challenge is to find the smallest m such that B X for the obtained set X satisfies conditions 1-2 of the GHZ scheme, i.e. such that m and X form a valid solution.For example, some ȳi could already be in Y, in which case the obtained B X would not satisfy conditions 1 and 2 for a valid solution to the GHZ scheme.With this approach, not all sets obtained X will thus necessarily be correct solutions, and we therefore need a verification step which is repeated until a valid solution is found.In the pseudocode provided in Algorithm 1, this part is done via brute-force starting from smaller values of m and looking if any of the possible valid partitions of Y would give raise to a set B X which satisfies conditions (1-3) of the GHZ scheme.This is repeated, increasing the value of m until a correct solution X is found, and the obtained values of X and m are finally returned.
Using the DSI algorithm we can easily obtain solutions for high-dimensional heralded GHZ state generation with up to N = 5 photons and, for N = 3, up to dimensionality d = 8.Examples for sets X and the number of modes m for such solutions are reported in Supplementary Table S1.For N = 2 the algorithm provides the same optimised solutions for any d as the ones reported in the main text (see Fig. 1b).For (N, d) = (3, 3) the solution found requires m = 32, which is slightly worse than the optimised solution with m = 25 obtained with the brute-force approach (see Fig. 1c).This indicates that, while the solutions given by the DSI algorithm are likely close to the optimised ones, they are in general non-optimal in terms of m.To further investigate this, Supplementary Fig. S1 shows in log scale the values of m for the solutions given by the DSI algorithm for different values of N and d.The curve for N = 3 shows a sub-exponential increase of m with the dimension (in particular, a scaling m ∼ d 2 for N = 3), which significantly improves the exponential increment of the analytical solution in Eq.S23, potentially making these schemes significantly more suitable for practical implementations.
We finally comment on possible routes to improve the DSI algorithm.First, we note that in the first part of the DSI algorithm in obtaining the set Y we insist that any different choice of N − 1 of its elements, including repetitions, must have a different sum.Such a constraint can be relaxed.In fact, a looser condition for the set Y is that there exists a partition P = {P 1 , . . ., P d } into d sets of N − 1 elements such that each of the sums σ i = p∈Pi p cannot be obtained by any other choice of N − 1 of elements of Y (including repetitions) that is not P i .Such a condition is sufficient to complete Y to a possible solution X as in the current algorithm.We believe this relaxation in obtaining Y is likely to provide better solutions in terms of the number of required resources m.In terms of computational resources required to find a solution, we note that the brute-force method to perform the second part of the DSI algorithm is the one that requires more memory usage and computational time, in particular when exploring all the possible partitions of Y. Finding heuristics that require to explore fewer partitions of Y are likely to improve this part significantly and enable the investigation of solutions for larger values of N and d.
Find new element ā to add to Y as the minimum positive integer not in Γ nor already in Y.

8:
Append ā to Y Update Y including new element ā.

9:
Append ā to S1 10: Iteratively update sets of sums S k after ā is added to Y. for σ ∈ Σ do 16: X ← Y ∪ {m − σi|σi ∈ σ} Complete the set with d different elements such that the partitions sum to m.

17:
Calculate BX given X and m.

18:
if BX satisfies conditions (1-3) of the GHZ scheme then return X and m If the obtained BX satisfies the GHZ scheme conditions, return X and m as a solution.
X satisfied condition Eq.S24.For such successful heralding pattern, we calculated the amplitudes a k associated to the terms |k, . . ., k in the qudits superposition.Each pattern was then accepted with probability d min k |a k | 2 , corresponding to the success probability of the Procrustean distillation.The estimation of the total success probability was finally obtained as the ratio between the number of accepted configuration and the total number of samples tested.
Repeating the procedure also provides an estimation for the uncertainty of the estimate.For the (N, d) = (3, 3) scheme in Fig. 1c we obtained an estimate of the success probability of 0.8(1) × 10 −4 , a six orders of magnitude improvement with respect to the case with a single heralding pattern.
To understand the practicability on near-term hardware for the schemes after this further success probability optimisation, we investigate the expected rates for the (N, d) = (3, 3) GHZ state generator with state-of-the-art hardware.In recent work from quantum dot research groups, a single photon generation efficiency over 92% has been reported with current technology [52].The state of the art for low-loss multiport interferometers is 99% transmission [53] and state of the art number resolving detectors have 98% efficiency [54].If we optimistically assume that these numbers can be achieved in a single experiment, which runs at the repetition rate of the quantum dot (145 MHz), we find that the optimised 0.8 × 10 −4 success probability corresponds to 677 successfully heralded and detected GHZ states per second.
The rules for generating GHZ states outlined in the main text allow us to prove that linear optical GHZ generation is possible for any (N, d).However, we have shown that as we relax these requirements, huge performance improvements can be found.We therefore expect that by further relaxing these requirements (i.e.different unitary matrices, different number of input photons), further improvements can be found for all (N, d).
of modes to the two qubits.However, considering all of these outcomes as successes would be inappropriate for some applications.For this reason, we define success probabilities for cases with fixed allocations, p f , and arbitrary allocations, p a .Fixed allocation success probability describes the probability of heralding a Bell pair when we fix which modes are allocated to each of the two qudits.Arbitrary allocation success probabilities defines the probability of heralding a Bell pair where the qudits can be described by any partition of modes.For arbitrary allocated success, it is likely that an optical switch would be required to rearrange the modes based on some heralding outcome.
As in Appendix 3, in some cases, the circuits described in this section can produce states which have entanglement with full Schmidt rank, but which are not maximally entangled due to unbalanced amplitudes.Such states can be probabilistically corrected to maximally entangled qudit Bell pairs using Procrustean distillation.To include this possibility, we also define a 'corrected' success probability, p c .Corrected success describes the total success probability of the gate if such corrections are allowed.These correction operations include the deterministic correction of the mode allocation of the qudits, and therefore we always have p c ≥ p a ≥ p f .

Heralded Bell circuit: Version 1
The heralded Bell pair generation circuit of Zhang et al. [58] uses a type-II fusion gate at its core.We aimed to generalise this circuit for creating heralded qudit Bell states by using the type-II fusion gate from Luo et al. [13].The layout of this circuit is shown for arbitrary dimension in Fig. S2a.Successful generation events are heralded by one photon being detected at the output of each DFT and a further d − 2 photons being detected in both the first and last d detectors.
The other circuits we present in the main text require ≤ 2d + 1 photons, whereas this circuit requires 3d − 2 single photon inputs and has success probabilities lower than other circuits we present, as shown in Table S3.We therefore believe this circuit to be of less practical interest compared to the other circuits we present.

Heralded Bell circuit: Version 2
This circuit comes from a more direct generalisation of the Bell pair circuit of Zhang et al. [58].This circuit uses 2d photons -which we conjecture to be the minimum number of photons required to produce a d-dimensional Bell state.The unitary, H 2d , shown in Fig. S2b, can be chosen to be any 2d × 2d complex Hadamard.In Table S3 we assess this circuit for the following three cases (labelled A, B, C): where P is a permutation which interleaves the outputs of blocks in one layer to the inputs for the next i.e. it the i-th mode of the j-th DFT in the first layer to the jth mode of the i-th DFT in the second layer.The optimal choice of H d depends on the type of success probability which is most useful for the success probability.Table S3 shows these success probabilities up to d = 5.

Heralded Bell circuit: Version 3
This circuit can be considered a generalisation of the Bell pair generation circuit from Zou et al. [59].It is also closely related to Version 2c.Reversing the first two layers allows you to swap between the two circuits.The circuit design is shown in Fig. S2c.Success probabilities are summarised in Table S3.Supplementary Table S3.Success probabilities for heralded qudit Bell state generation using different possible circuits.The labels v1, v2, v3 refer to the Bell state circuits described in sections 1, 2 and 3 of Appendix 6. A, B, C refer to the different options for the unitary H 2d as described in equations S29, S30, S31.ZTL refers to the zero-transmission-law based Bell pair generators described in Fig. 1b of the main text.Here we briefly describe how we simulated the state generation for qudit Bell states of dimensions 2 − 5 when using partially distinguishable photons in the heralded high-dimensional entanglement schemes described in the main text.As described below, the best methods to simulate the evolution of partially distinguishable photons in linear-optical interferometers calculate the probabilities for the various output configurations, but do not directly provide their amplitudes.To calculate the full quantum state we therefore use qudit quantum state tomography techniques that allow us to reconstruct the quantum state from the simulated measurement probabilities.

Qudit quantum state tomography
As mentioned in the main text, only simulations of bipartite entangled states with N = 2 remain tractable on a standard laptop.Therefore we here limit discussion to quantum state tomography of systems of two qudits.In general the density matrix for such a state may be written as [60,61] where G i are the generators of SU (d).These form an orthogonal operator basis in the space of linear operators for a qudit Hilbert space and obey the normalisation Tr(G i G j ) = 2δ ij .There are d 2 − 1 generators, which are the Pauli matrices for for d = 2 and the Gell-Mann matrices for d ≥ 3. We also include the rescaled identity G 0 = 2/d × I for a total of d 2 Hermitian operators with which to decompose a qudit.The coefficients r i,j can be calculated using the expectations of these generators by

Output probabilities with partial distinguishability
Given the generation circuit in Fig. 1b for qudit Bell states, we need to calculate the expectations of these operators G i for different values of photon indistinguishability | ψ i |ψ j | 2 as mentioned in the main text.To do so, we will need to find unitaries that rotate between the operator eigenbases and the qudit basis, since the latter is the one in which we detect photons.
We first take the spectral decomposition of the Hermitian operator photons...

2 .
B X forms a partition of X .3. For all b k ∈ B X , N i=1 y k,i = m.This is a slightly more restrictive form of the ZTL.

FIG. 2 .
FIG.2.Example of modular architecture for constructing universal cluster states from three-photon GHZ states with linear optics.Each module, shown in a, arranges three-photon GHZ state generators and fusion gates such that a single d-dimensional qudit is linked to the four neighbouring qudits in the universal square lattice in b.

Appendix 1 :
Calculation of output amplitudes and success probabilities in the N -photon d-dimensional heralded GHZ generator
the ZTL implies that the output configurations with non-zero amplitudes are only d, and are of the following form: m − N photons are the in the 0-th mode, and the remaining photons are in a configuration b k = (y k,1 , y k,2 , . . ., y k,N ), with y k,i ∈ {1, . . ., m − 1} and m i=1 y k,i = 0 mod m.Property 3 for the set χ also implies configurations satisfying m i=1 y k,i = m.Each of the d possible configurations b k is mapped directly into the logical state of the N d-dimensional qudits: b k → |k, k, . . ., k .In order to show that the states generated are in fact GHZ states, we will prove that the amplitude for each configuration b k is the same.For a given k ∈ {0, . . ., d − 1}, the amplitude of the associated state is given by k, k, . . ., k|Ψ = 1

Algorithm 1
Different-Sums Iterative (DSI) algorithmRequire: Number of qudits in the GHZ N , dimensionality for each qudit d.1:M ← N − 1 2: Y ← {1, . . ., M } Initialize the set Y. 3: for k ∈ 1 → M do 4: S k ← {y1 + ... + y k |yi ∈ Y}Calculate set of sums S k of all choices of k elements from Y (with repetitions).5:for n ∈ M → M d − 1 doIteratively adds elements to the set Y such that any subset of M elements, including repetitions, always has different sum.