Multi-dimensional super- and subradiance in waveguide quantum electrodynamics

We study the collective decay rates of multi-dimensional quantum networks in which one-dimensional waveguides form an intersecting hyper-rectangular lattice, with qubits located at the lattice points. We introduce and motivate the \emph{dimensional reduction of poles} (DRoP) conjecture, which identifies all collective decay rates of such networks via a connection to waveguides with a one-dimensional topology (e.g. a linear chain of qubits). Using DRoP, we consider many-body effects such as superradiance, subradiance, and bound-states in continuum in multi-dimensional quantum networks. We find that, unlike one-dimensional linear chains, multi-dimensional quantum networks have superradiance in distinct levels, which we call multi-dimensional superradiance. Furthermore, we generalize the $N^{-3}$ scaling of subradiance in a linear chain to $d$-dimensional networks. This work represents the first systematic study of collective behavior in waveguide quantum electrodynamics beyond one-dimensional topologies.

We study the collective decay rates of multi-dimensional quantum networks in which onedimensional waveguides form an intersecting hyper-rectangular lattice, with qubits located at the lattice points. We introduce and motivate the dimensional reduction of poles (DRoP) conjecture, which identifies all collective decay rates of such networks via a connection to waveguides with a onedimensional topology (e.g. a linear chain of qubits). Using DRoP, we consider many-body effects such as superradiance, subradiance, and bound-states in continuum in multi-dimensional quantum networks. We find that, unlike one-dimensional linear chains, multi-dimensional quantum networks have superradiance in distinct levels, which we call multi-dimensional superradiance. Furthermore, we generalize the N −3 scaling of subradiance in a linear chain to d-dimensional networks. This work represents the first systematic study of collective behavior in waveguide quantum electrodynamics beyond one-dimensional topologies.
Given this potential for realizing quantum networks in waveguide QED systems, it is perhaps surprising that the theoretical study of many-body effects in waveguide QED has been mostly confined to single [19][20][21][22][23][24][25][26] (see, for example, Fig. 1a) or coupled [7,[27][28][29][30][31] waveguides with linear topologies. Nonetheless, a multi-dimensional quantum network, as idealized in ), provides compactness and higher connectivity within the network with its many nodes and channels. Yet, such a fruitful concept has been untouched in the waveguide QED literature so far, perhaps due to the inefficiency of current analytical and computational techniques for generalizing to larger dimensions. On the other hand, a systematic theoretical study of many-body effects in large multi-dimensional networks will pave the way for developing complex quantum networks, and ultimately, the quantum internet [9]. With quantum computing becoming more of a reality with each passing day [4], it is now a crucial time to address this problem.
In this paper, we introduce a strategy for investigating previously unexplored multi-dimensional * fdinc@stanford.edu collective phenomena in waveguide QED. Specifically, we highlight a connection between multidimensional lattices and linear chains, and show that this connection can be used to compute the collective decay rates of multi-dimensional networks. We study collective phenomena such as super-and subradiance [11,32,33] and bound-states in continuum (BIC) [34,35] in these multi-dimensional systems. In our investigations, we discover the concept of multi-dimensional superradiance, where, unlike the well-known phenomenon discussed by Dicke [36], superradiance becomes segregated for these multidimensional networks. As a consequence, we show FIG. 1. Schematic of quantum networks in (a) d = 1, (b) d = 2, and (c) d = 3. Quantum emitters, located at waveguide intersections, function as nodes, which are two-level systems with energy gap Ω that generate, store and process quantum information. Loss-less 1D waveguides form quantum channels that transport information, in the form of quantum light, between the nodes. In (b), blue arrows denote coupling to horizontal waveguides and green arrows denote coupling to vertical waveguides.
that the dimensionality of BIC becomes smaller with increasing lattice dimension, presenting a trade-off between compact and efficient network design versus quantum memory capability. Moreover, we discover that the most subradiant decay rate in a ddimensional hyper-cubic lattices scales as N −3/d for large N , generalizing the known N −3 scaling in linear chains [26,32].
We consider a d-dimensional lattice of intersecting loss-less waveguides with linear topologies, 1 such as the ones illustrated in Fig. 1. Each waveguide along the nth direction contains N n identical qubits, and each qubit couples to a waveguide in the nth dimension with the decay rate γ n . The Hamiltonian for the entire system is H = H 0 + H I , where H 0 includes the self-energies of the light and qubits, and H I contains point-like interaction terms 2 located at the positions of the nodes [22].
Of particular interest in waveguide QED systems are the reflection and transmission parameters, r and t respectively, and excitation amplitudes e corresponding to a plane wave with momentum k. Collectively, these quantities are known as scattering parameters, and can be found from the Hamiltonian H by following the procedure outlined in Ref. [37].
Let us label each qubit using a d-dimensional lattice coordinate σ = (σ 1 , . . . , σ d ). If we concentrate on a particular direction n ≤ d with corresponding unit vectorn, we can define an adjacent qubit by the coordinate σ + an = (σ 1 , . . . , σ n + a, . . . , σ d ), where a is the lattice constant. The scattering parameters then satisfy the equations of motion (EoMs) which we derive in App. A. Within these equations, ∆ k = E k − Ω is the photon detuning energy, E k is the energy of the system, k is the momentum of the photonic degree of freedom and a is the lattice constant. Here we linearize the phase picked up by light propagating between two adjacent qubits such that ka Ωa = θ, which is accurate as long as time retardation effects inside the network are negligible. This assumption is valid in the Markovian regime, where the qubits are separated microscopically [37,38]. We note that, with our definition of scattering parameters, Eq. (1) is a high-dimensional generalization of the linear chain of qubits discussed in Eq. (6) in Ref [26]. Now, we turn our attention to the collective decay rates Γ, which contain information about the system's many-body structure. The collective decay rates are complex-valued, with their real components dictating the exponential decay of observables in time and their imaginary components capturing both the oscillatory behavior of the system as well as a characteristic frequency shift in energy levels [26].
The collective decay rates can, in principle, be found by solving the EoMs and examining the poles of the scattering parameters, which yields a polynomial characteristic equation for Γ [37]. In a network of N = d n=1 N n qubits, the behaviour of the entire system is thus governed by a total of (2d+1)N linear equations. For the special case of d = 1, the boundary of the quantum network does not scale with N as the linear chain has only one input and one output port. Hence, the equations of motion can be solved by eliminating all but the two boundary scattering parameters that are defined by the initial conditions. The rest of the scattering parameters can be found via back-propagation using transfer matrices. For this very special case, the scattering problem can be solved analytically and fairly efficiently for up to N = 500 qubits using the transfer matrix method [37]. For d ≥ 2, however, finding the scattering parameters by solving the EoMs becomes intractable for large N . In this case, even if the system's internal degrees of freedom (corresponding to scattering within the system) can be eliminated, as is usually done in d = 1 through transfer matrix methods [37], the number of external parameters (and, correspondingly, the number of equations to solve) scales with the size of the quantum network's boundary.
To overcome some of these limitations, we introduce an idea that we call the dimensional reduction of poles (DRoP), whereby we conjecture that there exists an effective mapping between the collective decay rates (which are obtained from the poles of the scattering parameters [37]) for the multi-dimensional network and waveguides with linear topologies. The DRoP conjecture makes it possible to find the decay rates of previously intractable higher-dimensional quantum networks: first, one applies the DRoP conjecture to divide the multi-dimensional network into a subset of linear chains, and then one finds the collective decay rates of these linear chains using efficient transfer matrix methods [37], obtaining the decay rates of the quantum network in the process. Through this method, one avoids performing calculations for d-dimensional, N -qubit networks directly, and instead works with a chain of size ∼ O(N 1/d ), for which it is possible to eliminate internal degrees of freedom.
To motivate the conjecture, we begin with a simple example using decay rates obtained analytically by solving the EoMs directly. Consider a d = 2 network with N n = 2 nodes and single emitter decay rates γ n along each direction n = 1, 2. In this case, there are four collective decay rates, given by (2) One can express these decay rates in terms of those corresponding to a d = 1 waveguide along the direction n with N = 2 nodes, each with single emitter decay rate γ n . For such a one-dimensional set-up, the decay rates are Here, the superscript in Γ corresponds to the direction n, whereas the subscript distinguishes between distinct decay rates. Comparing Eqs.
(2) and (3) reveals that the decay rates for the d = 2 system can be written in terms of the dimensionless decay rates z i /γ n (with i = 1, 2 denoting the decay rates along n = 1, 2 direction) for the linear system as Remarkably, as we show throughout this paper, an analogous construction appears to hold for arbitrary d and set of N n 's. Based on this observation, we propose the following conjecture. i /γ n denote the dimensionless collective decay rates along the direction n, where γ n is the single emitter decay rate corresponding to the same direction, such that 1 ≤ n ≤ d and 1 ≤ i ≤ N n . Then the complete set of collective decay rates belonging to the d-dimensional quantum network is with s = σ/a denoting the set of indices of the decay rates (with | s| = d) and |Γ| = N equal to the total number of qubits inside the network. Note that there are N unique sets s, each corresponding to a single collective decay rate Γ s .
We first demonstrate that decay rates obtained using DRoP match those found directly using the EoMs and later we use DRoP to discover new physics by probing regions inaccessible via EoM-motivated methods. We restrict our analysis to physically relevant dimensions d = 2 and d = 3, and make three different kinds of comparisons depending on the complexity of the system.
In the first comparison, we considered systems where N 1 , N 2 , N 3 ≤ 2 for d = 3, or N 1 , N 2 ≤ 3 for d = 2. In these cases, analytical results for collective decay rates can be found by directly solving Eq. (1). We compare these results to those obtained using DRoP and find that they agree exactly. The case for N 1 = N 2 = 3 is illustrated explicitly in App. B.
Next, for N ∼ 10, finding the decay rates analytically by solving the EoMs is intractable. DRoP, on the other hand, yields analytical results. As an example, we consider a 4 × 4 two-dimensional waveguide lattice. In Fig. 2, we show analytical results using DRoP compared with results obtained numerically from the EoMs (App. C contains details on the numerical approach used). We check all possible cases with N n ≤ 3 in d = 3 and N n ≤ 4 in d = 2, and find that the direct numerical results agree with our analytical DRoP results within machine precision.
Analytical results using DRoP are possible even for large N , but analysis becomes cumbersome due to the fact that the analytical expressions for the decay rates end up having many branch cuts in the complex plane. We therefore turn to using DRoP numerically, and make comparisons for particular values of θ. We check for various cases and again find that the decay rates found by numerically solving the EoMs directly and those found with DRoP (in combination with transfer matrix methods, see App. C) agree within machine precision. An example case for a three dimensional 5 × 3 × 4 waveguide lattice is shown in Fig. 3.
In all cases that we have checked, we find that the collective decay rates from DRoP match those found via EoM-motivated methods either exactly analytically or, when analytical comparison is not possible, to within machine precision. We also find that DRoP is quite robust to random noise introduced to the single-qubit decay rates (see App. D).
Guided by an intuition developed through studying the linear case [37], we notice in our explorations that for θ ≈ mπ (with m = 0, 1, . . .), the decay rates tend to cluster around certain regions of the complex plane. To understand the nature of this behavior, we now focus on physical phenomena such as superradiance and subradiance (BIC for when θ = mπ) for waveguide QED systems [37]. For the results that follow, we use DRoP to discover new physics by probing regions previously inaccessible.
Superradiance (subradiance) occurs when constructive (destructive) interference enhances (suppresses) spontaneous emission. Both physical phenomena are known to occur in a linear chain of qubits when θ = mπ [24,26,37]. When subradiance occurs, N −1 decay rates converge to the origin, whereas superradiance is when one of the decay rates converges to N γ, where γ is the single qubit decay rate. Out of N possible first excited states, N − 1 are dark states, i.e. states that have zero decay rate, owing to subradiance. Thus, one can construct the subspace containing the first excited states in terms of N − 1 dark states, which do not couple to electro-FIG. 2. Collective decay rates (in units of γ1) for a d = 2, 4 × 4 waveguide lattice (d = 2 and N1 = N2 = 4, with γ2/γ1 = 0.4). Each panel shows how particular decay rates behave in the complex plane for θ ∈ [0, 2π]. The coloured solid lines correspond to results obtained analytically using DRoP while the black crosses correspond to direct results obtained from the EoMs using the numerical method discussed in App. C).
FIG. 3. Collective decay rates Γi (in units of γ1) for a d = 3, 5×3×4 waveguide lattice (d = 3, N1 = 5, N2 = 3, and N3 = 4 for θ = π/2 with γ2/γ1 = 4, γ3/γ1 = 2). The system has 60 collective decay rates, all shown here. The circles correspond to results obtained numerically using DRoP while the black crosses correspond to direct results obtained numerically from the equations of motion using a numerical method known as the condition number method (CNM), as discussed in App. C. magnetic radiation, and one superradiant state that does. The collective system behaves as a two-level system between the one superradiant state and the ground state, which explains the Lorentzian shape of the transmission and reflection amplitudes, as discussed in [24]. We find that this behaviour no longer holds true for multi-dimensional networks.
We find that in a d-dimensional quantum network, the dimension of the superradiant subspace is larger than one. As a result, the collective system no longer behaves as a qubit and is hence no longer described by Lorentzian transmission and reflection amplitudes. Consequently, the superradiant and subradiant states emerge differently than for the linear case. In a d-dimensional system, the DRoP conjecture predicts n (N n − 1) subradiant states, with the rest showing superradiant features. We additionally find that, unlike in d = 1, for d > 1 superradiance is also segregated. We define n-dimensional superradiance as the case where the decay rates of qubits along n different directions are summed constructively. Then, there exist states with 1-, 2-, . . ., and d-dimensional superradiance, a previously unobserved phenomena which we call multi-dimensional superradiance. This segregation of the superradiant behavior is illustrated for a small network in Fig. 4 FIG. 4. Collective decay rates (in units of γ1) for a d = 3 system. Notice the difference in scale between the xand y-axes. We see that the decay rates cluster around the superradiant and subradiant values with vanishing imaginary part. Here, we used the parameters: N1 × N2 × N3 = 2 × 3 × 4, γ1 = γ2 = γ3 and θ = 0.9999π. As in Figs. 2 and 3, the black crosses correspond to results obtained from the numerical condition number method (CNM) discussed in App. C.
with θ = 0.9999π. We pick θ = π but close to π to show the dimensionality of each cluster. We use this small network to validate our DRoP results with EoM-motivated methods. While the segregated nature of superradiance may be observable without the application of DRoP, the multi-dimensionality aspect cannot. The grouping into 1D, 2D and 3D superradiance is only possible through the DRoP conjecture and is consistent among large structures that cannot be accessed through the numerical algorithm discussed in App. C. With DRoP, we find the origin of multi-dimensional superradiance, the dimension of each superradiant subspace as well as the strength of the corresponding d-dimensional superradiance.
Additionally, DRoP allows us to determine how a system's most-subradiant decay rate scales with its dimension. As pointed out in Refs. [26,32], for a linear chain of qubits, the most subradiant decay rate scales as N −3 . Here, we find that for a d-dimensional quantum network, the corresponding scaling is N −3 min with N min ≡ min n [N n ]. For a hyper-cubic lattice, N min = N 1/d , leading to a subradiant decay rate that scales as N −3/d . As expected, with increased dimensionality, subradiant behavior is washed out more gradually and the network couples to the con-tinuum more strongly. Since subradiance has been considered crucial for quantum memory applications [25,37], stronger design restrictions are expected to apply for higher dimensional quantum networks in comparison to linear chains.
For θ = mπ and in the absence of non-radiative decay (which is the case considered in this paper), the subradiant states become BIC (or dark states). Consequently, the dimensionality of BICs is equal to the dimensionality of subradiant states, i.e. n (N n − 1). The condition of BIC for a linear chain has been considered in Refs. [35,37]. As shown in App. E, the overall condition for BIC can be generalized to higher dimensions as Here, p c = 1 for even number of qubits along the chain and ±1 (alternating along the chain) for odd number of qubits along the chain. Since n (N n − 1) ≤ N − 1, the dimensionality of the BIC subspace in a d > 1 dimensional quantum network is smaller than the one belonging to N qubits in a linear chain. Therefore, it may be beneficial to use a lower-dimensional quantum network for memory applications. Here, again, design restrictions are expected to dictate the dimensionality of the quantum circuit being designed. Within this work, we have introduced the DRoP conjecture and illustrated its accuracy for a variety of examples via analytical and numerical methods. We have used DRoP to probe superradiance, subradiance and BIC in multi-dimensional quantum networks. We emphasize that, while EoM-motivated methods such as those discussed in App. C can provide some numerical information on the collective decay rates of small networks, our main results are derived from the analyticity that accompanies the DRoP conjecture and apply to multi-dimensional networks of arbitrarily large sizes. Previous research in waveguide QED mainly focused on linear structures and scattering parameters. While finding the scattering parameters efficiently in a large quantum network is still an open and important question, the DRoP conjecture opens the door for investigating multi-dimensional networks via their collective decay rates.
We have focused our studies on cases where time retardation effects resulting from the inter-system photon propagation are neglected. In future work, it will be interesting to consider whether the DRoP conjecture holds in regimes where non-Markovian and time-delayed quantum coherent feedback ef-fects become dominant [25,39,40] and, in particular, to check whether recently-discovered supersuperradiance effects for linear chains [41,42] persist for higher-dimensional networks. We believe that the proof of the DRoP conjecture lies in deriving the matrix equation for the spontaneous emission dynamics (discussed in Ref. [37] for a linear chain), which would allow for studies of time evolution in multi-dimensional quantum networks. As collective decay rates are linked to time evolution in waveguide QED [37], it is possible that the dimensional reduction produced by DRoP could be present for the time evolution as well. Due to their potential to aid in the investigation of multi-dimensional net-works and their time evolution, we expect DRoP to be useful for designing complex quantum networks for future quantum technologies.
Acknowledgments -FD acknowledges discussions with Diego Garcia and Jairo Rojas on Mathematica coding language. We also thank Gabriela Secara for help with preparing Figure 4 In this section, we describe the Hamiltonian and the corresponding stationary states for a d-dimensional quantum network. The free and interaction Hamiltonians are given by Here, a † σ is the excitation operator for the qubit whose position is given by the set σ = {σ 1 , ..., σ d }. γ n is the single qubit decay rate along the nth direction, ψ † n,m,L/R (x) is the bosonic creation operator for the left/right moving photons at position x in the mth waveguide along the nth direction, and m = m(n, σ). Ω is the energy separation of the qubit, N n is the number of atoms along direction n and |e σ = a † σ |0 is the excited state for the σth qubit with |0 being the vacuum state. γ n is a constant and not a function of frequency, which is inline with the assumption that we are interested in energies E k ∼ Ω ± O(γ n ), where γ n Ω [22]. Throughout the paper, we use natural units such that = v g = 1, where v g is the group velocity of photons inside the waveguides. We can construct a Bethe Ansatz as are the piece-wise field amplitudes and a is the lattice constant. At the boundary of the network, the field amplitudes include only one Heaviside function rather than two, meaning that the photon can radiate out of the system. This expression is a generalization of Eqs. (3-4) in Ref. [26]. t (n) σ and r (n) σ are, respectively, the transmission and reflection coefficients along the nth dimension belonging to the qubit σ, whereas e σ is the corresponding excitation coefficient. B.C. refers to the boundary terms of the photonic field, which can be hand-picked depending on the type of solution sought, due to the degeneracy of scattering eigenstates. Here we omit discussion the set of boundary conditions, as they do not have any effect on the collective decay rates. Applying the condition H |E k = E k |E k , we show below that we obtain the equations of motion (EoMs) given in Eq. (1) following the usual position space approach [19,26].
Applying the free Hamiltonian to the energy eigenstate, we find that where ∆ k = E k − Ω and E k = |k|. For now, we do not put much emphasis on the boundary terms, although their shape will emerge at the end of our calculations. Applying the interaction Hamiltonian gives Now, shifting the indices, re-arranging some terms and using field continuity at the atomic positions, we find that Here, σ + an is defined as (σ 1 , . . . , σ n + a, . . . , σ d ) and we note that such terms originate from the index shifting in summations. Within this equation, the shape of the boundary terms arise analogous to the 1D case [26,37], and the boundary terms include incoming and out-radiating photonic components. For |E k to be an energy eigenstate, all the other terms in Eq. (A5) should be zero, which leads to the EoMs given in Eq.
(1) of the main text.
There is another more straightforward and elegant proof for deriving local EoMs in waveguide QED systems with delta-function point interactions, for which we describe the strategy here. For a given system with many waveguides and qubits, one can divide the Hamiltonian into smaller pieces, with each piece containing a qubit and portion of all the waveguides that interact with it. These portions can be picked such that each is halved between the two adjacent qubits that are coupled to the same waveguide. Then, the Hamiltonian divides into sub-pieces such that where Q sums over all the qubits inside the system. Without loss of generality, H Q can be defined as Here, W stands for waveguides that interact with the qubit Q and the upper and lower bounds of the integrals are not relevant, as they depend on the sub-division of waveguides into H Q . As long as different sub-pieces are patched such that the photonic components are continuous at the patch points, no further equations of motion arise from the boundaries. Now, one can use the Bethe Ansatz approach to find the equation of motion around this single qubit as done in [19,26], there are 2N What makes this alternative proof more elegant is the fact that it does not use the specific geometry of the problem at hand. In fact, the local equations of motion are all the same for any type of waveguide QED system. The geometric properties of the system become important at the patching stage, where the inputs and outputs of patches should be properly defined to be continuous at the patch points (hence the irrelevance of space integral bounds). In a way, it is not the fundamental physics behind the EoMs that result in different emerging properties, such as the ones we discover in this paper, but rather the different phase relations that arise from patching in different geometries that lead to the different emergent phenomena.

Γ
(2) 4 = 1 2 2 + e 2iθ + e iθ 8 + e 2iθ γ 1 + 1 2 2 + e 2iθ − e iθ 8 + e 2iθ γ 2 , Re-writing these decay rates in terms of z (1) , we obtain the effective mapping predicted by the DRoP conjecture such that Similar analytical correspondence can be shown for a 3 × 2 × 2, 4 × 3 systems, both of which seem to be the boundary cases where Mathematica gives analytical results within less than a few hours on a standard personal computer.
Appendix C: Methods used in the main text for finding the collective decay rates Finding the scattering parameters by solving Eq. (1) for d = 1 and a general N is straightforward via the transfer matrix method [26]. Once the scattering parameters are known, the collective decay rates can be read off from the poles ∆ (0) = {∆ (0) k } of the scattering parameters. However, for d ≥ 2, solving for the scattering parameters become computationally intractable for even N ∼ 10. As an example, for {d, N 1 , N 2 , N 3 } = {3, 2, 3, 4}, one needs to solve a set of 168 coupled equations, which is not solvable on a standard personal computer within a time-span of few hours. As a result, in contrast to the 1-D case, solving for the scattering parameters is not a viable method to investigate the collective decay rates of a large multi-dimensional quantum network. For a 2-D quantum network, solving Eq. (1) becomes analytically intractable for N ∼ O(10) and numerically intractable for N ∼ O(100) on a standard personal computer.
Fortunately, if one is only interested in collective decay rates, instead of solving the linear system of equations multiple times, one can simply consider the matrix A, which contains the left hand side of Eq. (1), and find the set of poles ∆ (0) , for which A is singular. The decay rates can be found by rotation the poles in the complex plane via Γ = 2i∆ (0) [37].
Claim. The complete set of the poles of the scattering parameters is given by the values of ∆ k for which the matrix A is singular.
Proof. First, let us denote ∆ (0) as the ∆ k values for which A is singular. Moreover, let us introduce the notation A = A(∆ k ). For this sketch, we use proof by contra-positive, meaning we shall show that if A is non-singular, then ∆ k / ∈ ∆ (0) . If A is non-singular, then it is invertible. Let us call denote the inverse matrix by A −1 . Then, the solution to the matrix equation can be given as If A is invertible, then all entries of A −1 are finite. Similarly, all entries of b are finite by construction. Therefore, the scattering parameters are finite, since multiplication of two finite-dimensional matrices with finite entries results in a matrix with finite entries. By definition, if ∆ k ∈ ∆ (0) , then the scattering parameters should diverge by the existence of a pole, leading to a contradiction. Hence ∆ k / ∈ ∆ (0) . This result simplifies our search for the poles, as we no longer need to solve the system of equations. Moreover, this method gives us important information regarding the maximum number of poles. Specifically, if the matrix A is singular, then its determinant is zero. Now, the determinant of matrix A(∆ k ) is a polynomial with a degree of N , and hence there are at most N poles. This result is in line with the findings of the literature so far [26,37], since the number of poles is expected to be bounded by the number of atoms inside the system.
To find the poles of A, one can solve the condition det(A) = 0. While the determinant algorithm provides useful insight, its implementation is cumbersome as the determinant of the A matrix is a highly oscillating function of ∆k. Fortunately, one can probe the singularity of a matrix by its eigenvalues since a matrix A is singular if it has a zero eigenvalue. Thus, any pole ∆ (0) satisfies the property We shall denote the algorithm using this approach as the "eigenvalue method". While this method provides accurate results, it is slow. One can speed up the process by simply considering the condition number of the matrix A instead of its whole eigenvalue spectrum. Then, the pole can be given by the so called "condition number method" (CNM) as This equation gives a pole depending on the seeding of the algorithm. By seeding the algorithm many times, one can find all N distinct poles. The eigenvalue and condition number methods give the same results for two distinct cases: (a) 5×3×4, γ2/γ1 = 4, γ3/γ1 = 2 and θ = 0.5π, (b) 3 × 2 × 6 with max = 0, γ2/γ1 = 3, γ3/γ1 = 2 and θ = 0.65π with max = 0.05 (see App. D for the definition of max) . Note that these numerical plots correspond to configurations given in Fig. 3 and Fig. 6, respectively.
within numerical precision. Throughout this paper, we use the CNM to find the collective decay rates of a high-dimensional quantum network, since it is fastest.
In this work, we use the predictions made by the DRoP method to seed the condition number method, which is then used to find the poles of the quantum network. One might be concerned that the search algorithm finds a local minimum rather than a minimum corresponding to a zero. However, this concern can be addressed using the minimum modulus principle of complex analysis: Since f (∆ k ) = det(A(∆ k )) is an analytical function of ∆ k (specifically here, a polynomial), |f | can only have local minima at the position of its zeros.
To check whether the results from the CNM are valid, we compare three different methods. In the first method, we consider the log-absolute determinant value of the matrix A and show that there are indeed N distinct minima, each corresponding to a pole, as plotted in Fig. 5. Second, we use the eigenvalue method to verify the results found by the CNM. Finally, we use Mathematica's NRoots function when possible to provide another check for our results. In all cases, the values found agree with the ones predicted by DRoP conjecture.
Having shown how the numerical approaches mentioned in the text work, let us now focus our attention on finding the collective decay rates of N qubits in a linear chain efficiently. First, we present a compact two-equation system that includes all information about such decay rates. We then propose another method to find a polynomial characteristic equation of degree N .
When d = 1, the scattering problem has been solved via the transfer matrix method in Ref. [26]. According to calculations performed in this reference, the set of equations describing the poles are given by Here, 0 ≤ Re[λ] ≤ π is a complex parameter, γ is the single qubit decay rate and ∆ k are the poles of the scattering parameters. This result shows that one can define dimensionless poles such that z p = ∆ (0) k /γ is the same for any γ, which means that the collective decay rates of a linear chain of atoms depend on γ only linearly. Thus, we can define the dimensionless decay rates that describe N qubit in a chain, regardless of the specific value of γ.
In order to obtain analytical expressions for the decay rates, below we find a polynomial characteristic equation. To do so, we use the transfer matrix method, but employ a different approach than in Ref. [26] for higher computational efficiency. The transfer matrix for a single unit cell, which includes a qubit and a propagation phase θ = Ωa, is where χ k = γ/(2∆ k ). Here, t j and r j are the transmission and reflection coefficients for the jth atom. By construction, t 0 = 1 and r N = 0. Then, one can relate the output field amplitudes as From this relation, one can find the final transmission coefficient as The characteristic polynomial describing the poles of the system is where (S N ) 11 (χ k can be easily found.