The concept of synthetic dimensions has recently opened the door to novel perspectives for expanding the dimensionality of well-understood physical systems [1–5]. One strategy to explore synthetic dimensions consists in driving the associated dynamical systems in order to activate the coupling between different internal modes, which under normal conditions remain uncoupled [6]. By doing so, the resulting coupled modes exhibit lattice-like structures that exist in an abstract space, which is nonetheless physical. The importance of synthetic lattices lies in the fact that they allow us to explore a variety of effects that are not available in spatial or temporal domains.

To illustrate the basic idea of activating synthetic dimensions, and to set the stage for the present work, we begin by elucidating how a 1D quantum harmonic oscillator generates a lattice in Fock space. The oscillator’s Hamiltonian is given as H^=ω(a^†a^+12), and its dynamics is governed by the Schrödinger equation i∂t|Ψ(t)⟩=H^|Ψ(t)⟩. Here, ω is the angular frequency of the oscillator, and a^ and a^† denote, respectively, the annihilation and creation operators [7]. Note that we have set the reduced Planck constant and the oscillator mass to unity, i.e., ℏ=1 and mo=1. When the oscillator is initially prepared in the eigenstate |Ψ(0)⟩=|n⟩, it will remain in this state, only acquiring a time-dependent phase factor during evolution, i.e., |Ψ(t)⟩=e−i(n+12)ωt|n⟩. No transitions to other eigenstates occur. However, by subjecting the oscillator to a time-dependent displacement, x^(t)=f(t)(a^†+a^), the Hamiltonian acquires the form H^(t)=ω(a^†a^+12)+f(t)·(a^†+a^). Substituting the general state vector |Ψ(t)⟩=∑m=0∞cm(t)|m⟩, where cm(t)=⟨m|U^(t)|Ψ(0)⟩ are the transition amplitudes from the initial state |Ψ(0)⟩ to the final state |m⟩ and U^(t) is the time evolution operator, into the Schrödinger equation, we find that the amplitudes cm(t) obey the semi-infinite set of coupled differential equations These equations clearly illustrate that the time-dependent displacement x^(t) activates transitions among the amplitudes cm−1(t), cm(t), and cm+1(t). This implies that, in Fock space, the oscillator generates a lattice, where it can “hop” from eigenstate |m⟩ to the adjacent eigenstates |m−1⟩ and |m+1⟩ with hopping rates f(t)m and f(t)m+1, respectively [8–14].

In general, applying dynamic modulations to the potentials associated with physical systems induces coupling among the supported eigenstates. Using this technique, a photonic topological insulator in synthetic dimensions has been recently implemented via modulated waveguide lattices [15,16]. Synthetic dimensions have also been explored in harmonic traps [17], optical lattices [18], cavities [19], and even in room-temperature Rydberg atoms [20].

Within the realm of optics and photonics, synthetic dimensions can be created by exploiting the spatial, temporal, polarization, and frequency degrees of freedom of light [6]. For instance, large-scale parity-time symmetric lattices have been implemented in the temporal domain using optical fiber loops endowed with gain and loss [21,22] and a driven-dissipative analogon of the 4D quantum Hall effect has been observed in a spatially 3D resonator lattice [4].

In this work, we show that high-dimensional lattices emerge in photon-number space when a photonic lattice of M ports [23,24] is excited by N indistinguishable photons (see Fig. 1). More precisely, the Fock-representation of N-photon states in systems composed of M evanescently coupled single-mode waveguides yields to a new layer of abstraction, where the associated states can be visualized as the energy levels of a synthetic atom, which features a number of allowed and disallowed transitions between its energy levels.

In photonic waveguide lattices, where all the waveguides are coupled to each other, the quantum optical Hamiltonian in paraxial approximation is given as H^=∑j=1Mβja^j†a^j+∑i≠jMκija^i†a^j [25], where a^j† and a^j, respectively, are bosonic creation and annihilation operators for photons in the jth waveguide. Further, βj denotes the propagation constant of the jth waveguide, and κij is the coupling coefficient between the ith and jth waveguides.

For simplicity, we restrict our subsequent analysis to the simplest scenario of (in real space) essentially 1D waveguide arrays with nearest-neighbor couplings

Under these premises, the propagation of a single-photon along the waveguides can be described using the Heisenberg equations of motion for the bosonic creation operators [26,27] where m=1,…,M. Accordingly, the single-photon response is computed through the input–output transformation a^n†(0)→∑m=1MUn,m(z)a^m†(z), where Un,m(z) denotes the (n,m) matrix element of the evolution operator U^(z)=exp(−izH^) [28]. Using this formalism, it is straightforward to show that an initial N-photon state |n1,n2,…,nM⟩, with N=∑m=1Mnm, will transform into the output state In the context of waveguide lattices, the input–output formalism is by far the most common approach used to compute the output states [29]. Nonetheless, as we will demonstrate in the remainder of the paper, the input–output scheme fails to expose the intrinsic coupling interactions between the emerging states.

In what follows, we use the equivalent Schrödinger-picture formalism to unveil the high-dimensional lattice structures arising from the propagation of multiple photons through multiport waveguide systems. To do so, we first notice that N indistinguishable photons exciting M coupled waveguides, give rise to a total of NF=(N+M−1)!/(N!(M−1)!) states which are given by all permutations of the integer partitions of N among the M sites.

For the trivial case of N=1 photon, we simply obtain a set of M states with m=1,…,M. By computing the matrix elements of the Hamiltonian given in Eq. (3) for N=1, Hn,m=⟨1n|H^|1m⟩=βnδn,m+κn,m−1δn,m−1+κn,m+1δn,m+1, one can readily see that the single-photon states are coupled to each other as displayed by the equations in agreement with Eq. (4).

We now consider the more interesting scenario of N photons propagating through a waveguide beam splitter, M=2, with propagation constants β1 and β2 and symmetric coupling, i.e., κ1,2=κ2,1≡κ. In this case, there exists a total of (N+1) states, namely, (|0,N⟩,|1,N−1⟩,…,|N−1,1⟩,|N,0⟩), and the Hamiltonian given in Eq. (3) acquires the form

Computing the matrix elements H^(m,n),(p,q)=⟨m,n|H^|p,q⟩ reveals that the states obey the (N+1) equations of motion with Cm=κm(n+1) and n=N−m [30]. This indicates that, inside a waveguide beam splitter, the amplitudes of two-mode N-photon states evolve coupled to each other with hopping rates Cm, and the corresponding phases depend on both propagation constants.

For the case of two identical waveguides, we have β1=β2=β so that the first term on the r.h.s. of Eq. (9) becomes βN|m,N−m⟩, which indicates that all the states will exhibit the same effective propagation constant. Interestingly, it has been recently shown that waveguide beam splitters produce the discrete fractional Fourier transform (DFrFT) of N-photon states [30] as well as exceptional points of arbitrary order, provided that losses are introduced in one of the waveguides [31].

On the other hand, when considering two nonidentical waveguides, β1≠β2, the first term on the r.h.s. of Eq. (9) acquires the form ((β1−β2)m+β2N)|m,N−m⟩. Remarkably, the term ((β1−β2)m) indicates that the state evolution will be influenced by an effective ramping potential in the same fashion as in the case of classical waves in Bloch oscillator systems [12,32,33]. Consequently, we can tailor the dynamics of N-photon states by simply adjusting the Bloch slope (β1−β2) in order to suppress and/or create certain output states. As an illustration, we depict in Fig. 2 the probability evolution for the initial state |5,5⟩ in a waveguide beam splitter with coupling coefficient κ=1 for β1=β2=1 [Fig. 2(a)] and  for β1=0,β2=4 [Fig. 2(b)]. While the case of Fig. 2(a) corresponds to discrete “diffraction” of the initial state in state space, the case of Fig. 2(b) corresponds to “Bloch oscillations” in state space. Note that, throughout this work, we present all simulations using the normalized propagation coordinate z=κZ, where Z is the actual propagation distance, and κ stands for the nearest-neighbor coupling coefficient. After the above introductory examples, we now proceed to consider the most interesting case where multiple photons N>1 excite more than two waveguides M>2. In order to motivate the concept of pseudo-energy, we first examine the simplest case of a waveguide trimer, M=3, that is excited by N=2 photons and then move on to the general case.

For a waveguide trimer and two identical photons, the Hamiltonian takes the form In this scenario, we have a total of six photon-number states obeying the following coupled set of equations of motion:

As in the earlier examples, here we also have the possibility of molding the state dynamics via tuning the propagation constants and coupling coefficients. For instance, for equal coupling coefficients κ1=κ2=1 and identical waveguides β1=β2=β3=0, we observe periodic spreading and contraction of the two-photon wave function, as illustrated in Fig. 3(a). In contrast, choosing a different propagation constant for the central waveguide, β2=2, leads to a quasi-periodic evolution [Fig. 3(b)]. Indeed, this quasi-periodic evolution occurs because the ratios between the eigenvalues of the coupling matrix are irrational numbers. We would like to emphasize that, at the propagation distance indicated by the dashed line in Fig. 3(b), the input state |101⟩ evolves into a quasi-two-photon NOON state in state space, which is reminiscent of the Hong–Ou–Mandel effect [34].

To describe the photon dynamics in the waveguide trimer, we have obtained an even number of equations. At this point, the way in which the states should be arranged into a synthetic lattice is not at all clear. To be precise, the six states representing the sites of the synthetic lattice can be sorted into at least two distinct natural sequences, as shown in Table 1.

Clearly, arranging the states into a lattice (i.e., sorting) and analyzing the corresponding equations of motion becomes rather cumbersome when considering higher photon numbers in multiple coupled waveguides. In the following section, we therefore introduce a concise and universal method that facilitates studying the general case of N>1 photons propagating in arrays formed by M>2 waveguides. The resulting structures follow from physical and mathematical considerations that eventually allow us to describe multiphoton processes in waveguide arrays in a surprising and remarkable way that resembles the quantum-mechanical description of multilevel atoms.

We now introduce a concept analogous to the concept of energy, which we refer to as pseudo-energy. As we show below, the concept of pseudo-energy is rather useful since it facilitates a unique sorting of multiphoton Fock states in a physically meaningful way and allows for establishing correspondence between Fock states and the energy levels of a synthetic atom. Concurrently, we identify pseudo-energy ladder operators along with pseudo-exchange energies in order to define the corresponding selection rules in Fock space for transitions between the pseudo-energy levels of the synthetic atom.

We consider N indistinguishable photons propagating in an array of M lossless evanescently coupled waveguides, which give rise to NF=(N+M−1)!/(N!(M−1)!) Fock states |n1,…,nM⟩, fulfilling the condition ∑m=1Mnm=N. The first issue to be addressed is to determine a way to sort the multiphoton states in Fock space in a meaningful way. To do so, we associate a unique numerical value to every state |n1,…,nM⟩ as follows: Here, the subscript N+1 indicates that the numbers in the square brackets have to be expressed in base N+1, and the least-significant digit is the left-most number n1. Observing that [n1,…,nM]N+1=∑m=1M(N+1)m−1nm allows us to define the pseudo-energy operator such that its action on the N-photon–M-mode Fock states |n1,…,nM⟩ yields with eigenspectrum K(n1,…,nM)=∑m=1M(N+1)m−1nm. From Eq. (19), we readily infer the smallest and largest eigenvalues Kmin=K(N,0,…,0,0)=[N.0.….0.0]N+1=N and Kmax=K(0,0,…,0,N)=[0.0.….0.N]N+1=N(N+1)M−1, respectively. Accordingly, the eigenvalues are bounded by Kmin≤Kν≤Kmax.

As a result, in order to sort the associated Fock states, we have to compute the corresponding Kν and arrange them in ascending order. The resulting ladder of Kν then defines the synthetic lattice formed by the states. We refer to this ordering as the pseudo-energy representation of the N-photon–M-mode Fock states.

For illustration, we revisit the above case of N=2 photons propagating in an array of M=3 waveguides. Accordingly, there are NF=6 states, and the spectrum of the pseudo-energy operator K^(2,3) comprises six integers: Using these numbers, we readily obtain the pseudo-energy representation of the two-photon–three-mode Fock space: Consequently, we designate Kν as the pseudo-energy of the νth Fock state in the N-photon–M-mode Fock space with ν=1,…,NF. In general, for any given N, M, and pseudo-energy Kν, the inverse mapping onto the mode-occupation numbers is where the symbol / corresponds to integer division and # is the modulo operator.

We now proceed to show how the pseudo-energy representation of Fock states allows us to express the equations of motion of N photons in M waveguides in a concise way. To do so, we take a closer look at the action of the operator a^i†a^j on a Fock state: If the state |n1,…,nM⟩ corresponds to the pseudo-energy Kν, then the resulting state on the r.h.s. of Eq. (24) must have the pseudo-energy Therefore, the action of a^i†a^j changes the pseudo-energy of Fock states by the amount which we denote as the pseudo-exchange energy associated with the tunneling process taking place between waveguides i and j. In this sense, the operators a^i†a^j can be thought of as pseudo-energy ladder operators, which raise or lower the pseudo-energy of Fock states. Consequently, we can write The physical significance of Eq. (27) is that a direct transition between the states |Kμ⟩ and |Kν⟩ is only possible if there exists a pseudo-exchange energy ΔKij such that Obviously, Eq. (28) defines the selection rules in Fock space. Together with the action of the photon number operators n^m, the full system of coupled equations governing the propagation of N photons through M coupled waveguides in the pseudo-energy representation is given by For the case of nearest-neighbour coupled, identical waveguides, where all the propagation constants are the same, the relevant pseudo-exchange energies are ΔKi=ΔKi+1,i=N(N+1)i−1, and the set of coupled equations reduces to To further illustrate the resulting coupling system in Fock space, we revisit the case of a single photon N=1 propagating in M=3 waveguides. The effective coupling behavior, of allowed and forbidden transitions in Fock space, can now be visualized within a pseudo-energy term diagram, as illustrated in Fig. 4(a). In this particular case, the nearest-neighbour coupling of the waveguides is retained in Fock space, and any given Fock state |Kν⟩ only couples to its nearest neighbors |Kν±1⟩.

A similar picture arises in the case of two waveguides M=2 and N=2 photons, as depicted in Fig. 4(b). Here, we obtain a term diagram that is essentially isomorphic to Fig. 4(a), where, again, only nearest-neighbor Fock states are coupled to each other.

The nearest-neighbor picture radically changes when applying the pseudo-energy approach to the case of N=2 photons and M=3 waveguides, as displayed in the corresponding term diagram in Fig. 4(c). Importantly, even when the waveguides are, in real space, only coupled to their nearest neighbors, in photon number space certain states become coupled to next-nearest neighbor states. For instance, in Fig. 4(c), we observe that the state |K2⟩=|4⟩=|1,1,0⟩ not only couples to its neighbors |K1⟩=|2⟩=|2,0,0⟩ and |K3⟩=|6⟩=|0,2,0⟩ but also to the next-nearest neighbor state |K4⟩=|10⟩=|1,0,1⟩. For illustrative purposes, we present in Fig. 5 the coupling matrix for this particular set of states when the three-waveguide system is formed by identical waveguides, β1=β2=β3=0, and balanced coupling coefficientsκ1=κ2=1.

At this point, it is rather evident that the richness and complexity of the emerging synthetic configurations will become more prominent when a larger number of photons and waveguides are considered. Moreover, it is worth stressing that, in order to generate the present synthetic structures, we did not require any modulation of the system parameters, as the states naturally couple due to the system’s internal dynamics.

In this section, we introduce a more convenient way of representing the Hamiltonian matrix of N-photons exciting M-waveguides. To do so, we interpret the states as vertices of a graph (Fock graph) where the allowed interstate transitions represent the edges. A practical representation of finite graphs is the so-called adjacency matrix whose entries indicate whether or not pairs of vertices are adjacent. In the present context, the effective Hamiltonian Hμν=⟨Kμ|H^|Kν⟩ in the N-photon–M-mode pseudo-energy representation determines such an adjacency matrix where Θ is the step function, and Aμν(N,M)=1(or0) indicates a connection (or no connection) between the vertices μ and ν. In what follows, we assume identical waveguides with β1=⋯=βM=0 in order to omit self-loops in the graph representation. As an example, in Fig. 6(a), we depict the Fock graph arising from the effective Hamiltonian of Fig. 5, which we have already discussed in the previous section. In Fig. 7(a), we depict further examples for photon numbers up to N=5 and up to M=6 waveguides. The first row, which corresponds to single-photon graphs, simply reflects the 1D spatial configuration of the waveguides. By introducing a second photon, we observe that the Fock graphs become 2D [Fig. 7(b)], except for the case M=2. The inclusion of more photons leads to nonplanar graphs, i.e., graphs that cannot be drawn in 2D without intersecting edges, which exhibit a layered structure in three dimensions as indicated by the different coloring of the nodes in different layers.

A prominent feature to highlight is the symmetry observed among graphs emerging for the combinations (M,N) and (M−l,N+l) and for (M,N) and (M+l,N−l), where l is an integer. In other words, every Fock graph has an isomorphic partner graph with an identical adjacency matrix, up to a trivial permutation of the node labels. In Fig. 7(b), we depict the smallest nontrivial pair of Fock graphs and the corresponding adjacency matrices that are induced by the pseudo-energy representation. If we were to start from Aμν(3,3) and permute its rows and columns according to (1,…,10)→(1,2,4,7,3,5,8,6,9,10), we will exactly obtain Aμν(2,4).

Indeed, this underlying symmetry in the space of possible Fock graphs has interesting implications. For instance, in Ref. [30] we have shown that it is possible to implement the number-resolved (N+1)-dimensional discrete fractional Fourier transform (DFrFT) with a single waveguide beam splitter by launching N indistinguishable photons. Furthermore, using the same photon-number-resolved mapping in Ref. [31], we have shown how to attain so-called exceptional points of N+1 order, by way of exciting a semi-lossy waveguide beam splitter with high photon number states. In fact, it is now clear that these results emerge as special cases of Eq. (32), which pertains to the identity of the first row and column in Fig. 7(a). Thus, by following similar ideas, it is possible, in principle, to find the corresponding effects for waveguide systems with M≥3 excited by N≥2 photons.

Additionally, by exploiting the graph symmetry, it becomes apparent that a specific transformation, which requires N photons and M waveguides could likewise be implemented with M−1 photons and N+1 waveguides. Of course, such alternative pathways of implementing a transformation are not always guaranteed because of the different dimensions of the experimentally accessible parameter spaces. Nonetheless, this may serve as a useful Ansatz to overcome concrete experimental difficulties.

Quite interestingly, Wang et al. explored synthetic Fock lattices in the context of QED circuits [19]. In such a study, the joint excitation states of an atom coupled to the N-photon three-cavity Fock space form a 2D, hexagonal Haldane-like synthetic lattice, which facilitates the generation of high-photon-number NOON states. Crucially, the realization of this scheme demands the judicious implementation of the coupling between atom and cavity as well as the precise modulation of the cavity resonance frequencies. In contrast, the multiphoton synthetic dimensions explored in the present work are intrinsically active by virtue of the indistinguishability of the photons; as such, they do not require any external driving of the system’s parameters.

The Fock graphs offer a rich variety of synthetic coupled structures. This variety can be further enhanced by considering different spatial arrangements of the waveguides, for instance, ring- or star-shaped structures instead of the simple planar configuration studied here. Importantly, the evolution of multiphoton states in synthetic lattices and graphs can be dynamically reconfigured by using programmable photonic chips [35], that is, integrated optical devices where the waveguides’ refractive index and coupling coefficients can be modified externally. Nevertheless, even with this simple 1D arrangement comprising a few waveguides, small photon numbers, and a time-independent Hamiltonian, one encounters interesting effects that are only possible due to the multidimensionality of the corresponding Fock graphs.

To show possible applications of the pseudo-energy synthetic lattices, we discuss the generation of all-optical dark states [36] and parallel multiphoton quantum random walks. The simplest dark states are encountered in three-level atomic or molecular systems, where radiative transitions between, e.g., |1⟩↔|2⟩↔|3⟩, are allowed, but the transition |1⟩↔|3⟩ is forbidden. In this simple scenario, a dark state is a superposition of the uncoupled states |D⟩=cos(θ)|1⟩−sin(θ)|3⟩, where θ is given in terms of the Rabi frequencies of the allowed transitions [36]. Once the system is in such a state, adiabatic changes in the Rabi frequencies allow for the tuning of the populations of the states |1⟩ and |3⟩, while the probability of |2⟩ remains 0. This interesting behavior, which seemingly evades the radiative selection rules, forms the basis of the method of stimulated Raman adiabatic passage (STIRAP) and has been applied in a variety of physical contexts, such as atoms, molecules, electrons, photons, magnons, and phonons [37]. In what follows, we demonstrate that combining the principles of STIRAP with the pseudo-energy representation of Fock states allows the generation of complex all-optical dark states, i.e., classes of physical multiphoton states in waveguide arrays that are immune to discrete diffraction effects. These states obviously lend themselves to numerous applications in quantum information science [36].

To do so, we revisit one more time the case of M=3 waveguides, with equal propagation constants β1=β2=β3=0 and balanced coupling coefficients κ1=κ2=12, excited by N=2 photons. The pseudo-energy representation of the effective Hamiltonian takes the form With these parameters, the spectrum of Hμν is integer-valued which indicates that the third and fourth eigenstates are degenerate with eigenvalues λ3=λ4=0. We now consider the evolution of a coherent superposition |ψ⟩ of the eigenstates with corresponding eigenvalues λ3=0 and λ5=1, specifically In the standard Fock representation, |ψ⟩ reads as

The probability evolution for this state is shown in Fig. 8. As one can see, this state displays the characteristic behavior of a dark state. That is, the initial state evolves exhibiting oscillating transitions between the states |200⟩ and |002⟩ with period 2πλ5−λ3=2π. These transitions occur in spite of the fact that the direct transition |200⟩↔|002⟩ is forbidden (⟨200|H^|002⟩=0), and those states have the maximum possible distance within the graph, that is, at least four single-photon tunneling processes are required to transform one state into the other. All probabilities of the intermediate states remain constant and, in a way, assist the simultaneous tunneling of two photons between the outermost waveguides. We stress that this six-level dark state is induced by a time-independent Hamiltonian; further, it occurs naturally without the need of adiabatic fine-tuning of the external parameters. We would also like to note that the state |020⟩ exhibits zero probability for all z, further attesting a multiphoton tunneling (in this case co-tunneling) effect taking place between the two waveguides. Geometrically speaking, this effect arises due to destructive interference taking place in the two-way branching of the Fock graph shown in Fig. 6(a). This branching effectively allows for the flow of the amplitudes to take a “detour” around the |020⟩ node.

As an alternative, one may attempt to implement a real space structure in one or two dimensions consisting of six coupled waveguides in order to emulate an equivalent Hamiltonian for just a single photon; in Fig. 6(b), we depict one such attempt. However, this would be topologically impossible since there always exists additional crosstalk between the waveguides representing the nodes at the center of the graph. Furthermore, from a mathematical point of view, the dark-state behavior of |ψ⟩ is based on a critical relationship between the components un(3) and un(5) (n indicating the row of the element) of the nondegenerate eigenvectors |ϕ3⟩ and |ϕ5⟩. Specifically, as one can see in Eq. (35), the interference terms In=un(3)*un(5)+un(5)*un(3)=0 vanish for all n, except for n=1 and n=6. Consequently, the probabilities of the pseudo-energy states |K2⟩,…,|K5⟩ remain constant during propagation. From this, we conclude that a coupled structure can only support all-optical dark states if there exist pairs of nondegenerate eigenstates, whose mutual interference terms vanish except for two components. Thus far and despite intense efforts, we have only been able to observe this specific property in Fock graphs with N≥2, M≥3. In other words, our Fock graph analysis of multiphoton propagation in waveguide arrays allows the realization of functionalities beyond those that can be realized with planar (single-photon) networks.

Quite interestingly, exciting waveguide lattices with multiphoton states comprising infinite coherent superpositions, e.g., coherent states |α⟩=exp(−|α|2/2)∑n=0∞(αn/n!)|n⟩ or two-mode squeezed vacuum states |ξ⟩=1−|ξ|2·∑n=0∞ξn|n,n⟩, opens a route to generating, in principle, an infinite number of lattices or graphs with different numbers of lattice sites and many dimensions simultaneously. This possibility is appealing for realizing parallel quantum random walks where the corresponding walkers can perform different numbers of steps on different Fock graphs that depend on the number of photons involved in each process.

We stress that the observation of the above effects is possible by utilizing on-chip direct laser-written waveguides [38], which are fed by bright parametric downconversion light in combination with photon-number-resolving detectors [39]. This type of waveguide can be designed to work at the typical wavelengths of 815 nm as well as 1500 nm; further, waveguide-coupling strengths can be engineered to range from κ=0.5cm−1 up to κ=2.5cm−1. Moreover, such photonic chips can be easily implemented to host a large number of waveguides with effective propagation distances ranging from several micrometers up to 15 cm [24]. As multiphoton sources, nonlinear periodically poled potassium titanyl phosphate (ppKTP) waveguides are available. Using such systems, the generation of two-mode quantum light with a mean photon number of 50 and a maximum number of 80 photons for each of the two modes has been demonstrated recently [40]. Finally, in order to translate the pseudo-energy states into real-space measurements, we require to resolve the number of photons in each waveguide. This may be achieved via so-called photon-number-resolving transition edge sensors, i.e., superconducting devices capable of detecting visible and near-infrared light at the single-photon level [41].

In this final section, we obtain an analytical expression for the eigensystem of an M-waveguide system (or tight-binding network) with arbitrary coupling coefficients κm excited by N indistinguishable photons. With the help of the pseudo-energy representation, we will be able to find a concise expression, which also introduces a natural ordering of the N-photon–M-waveguide eigenstates. As we have seen, in the case of a single-photon N=1, the Hamiltonian takes on a bi-diagonal form in the pseudo-energy representation. In some cases, it is possible to find an analytical closed-form expression for the eigensystem, as for example in the case of the DFrFT [42]. Even if no analytical solution is available, numerical algorithms are known [43] that deal with bi-diagonal matrices efficiently. Therefore, without loss of generality, we assume that we know the complete eigensystem of the single-photon–M-waveguide Hamiltonian, which we denote as where n=1,…,M. In the above equation, um(n) is the mth component of the nth eigenvector of the matrix H^m,n=⟨Km|H^|Kn⟩, and it defines the single-particle eigenstates When the same waveguide system is excited by N>1 photons, it is clear that the many-particle eigenstates arise from the tensor products of the single-particle eigenstates. Formally, we may write the resulting states as but now the occupation numbers n˜m pertain to the number of photons occupying the mth single-particle eigenmode. Consequently, we can apply the pseudo-energy ordering to the N-particle eigenstates by defining K˜ν=[n˜1(ν),…,n˜M(ν)]N+1. The νth eigenstate of the N-photon system is then given by Note that, in most cases, it is necessary to normalize the resulting expression on the r.h.s. of Eq. (42). By requiring |K˜ν⟩=∑μ=1NFcμ(ν)|Kμ⟩, where |Kμ⟩ denotes N-photon–M-waveguide Fock states, we find for the components cμ(ν)It is now rather straightforward to show that the N-particle eigenvalues are given as the sum of the eigenvalues of the involved single-particle eigenstates Using Eqs. (42) and (44), it is straightforward to find the N-photon–M-waveguide time-evolution operator U^(t)=∑ν=1NFe−iλ˜νt|K˜ν⟩⟨K˜ν|. We would like to emphasize that the numerical evaluation of Eq. (42) is far more efficient than the direct diagonalization of the full matrix representation of H^ in N-photon–M-waveguide Fock space. Due to the size and highly nontrivial structure of the resulting matrices, general eigensystem solvers produce a significant amount of overhead, which we avoid in our approach. Essentially, we do not even require a calculation of the full matrix representation Hμν. Instead, knowledge of the single-particle eigensystem and the bosonic nature of photons suffices.

In summary, we have shown that the propagation of multiphoton states through multiport waveguide systems (tight-binding networks) gives rise to multiple synthetic lattices and multidimensional Fock graphs that allow for transparent analyses of the relevant physical processes and the design of novel functionalities, such as multilevel all-optical dark states, beyond the linear (single-photon) realm. Since such synthetic structures emerge in the photon-number space, we have been able to associate coherent multiphoton processes to parallelized multidimensional quantum random walks. This parallelization brings about novel opportunities for the implementation of random walks, where the randomness is not only present in the dynamics of the walkers but also in the simultaneous occurrence of different walks.

Acknowledgment. We acknowledge support by the Deutsche Forschungsgemeinschaft (DFG) within the framework of the DFG priority program 1839 Tailored Disorder. We thank Friedemann Motzkus for helpful discussions on the encoding of Fock states. R.J.L.M. thankfully acknowledges financial support by CONACYT under the project CB-2016-01/284372 and by DGAPA-UNAM under the project UNAM-PAPIIT IN102920.