$\displaystyle=-i(J_{0}+i\Omega)\ket{\psi_{R}(t)}+\sqrt{J_{0}}C_{R}^{\dagger}(t)\ket{0},$
(64b) $\displaystyle i\partial_{t}\ket{\psi_{L}(t)}$
$\displaystyle=-i(J_{0}+i\Omega)\ket{\psi_{R}(t)}+\sqrt{J_{0}}C_{L}^{\dagger}(-t)\ket{0}.$
(64c)

Next, let us find the action of the free Hamiltonian ($H_{0}$) on each part of
$\ket{\psi(t)}$

$\displaystyle H_{0}\ket{\psi_{q}(t)}$
$\displaystyle=\Omega\ket{\psi_{q}(t)},$ (65a) $\displaystyle
H_{0}\ket{\psi_{R}(t)}$
$\displaystyle=-i(J_{0}+i\Omega)\ket{\psi_{R}(t)}-\sqrt{J_{0}}e^{-(J_{0}+i\Omega)t}C_{R}^{\dagger}(0)\ket{0}+\sqrt{J_{0}}C_{R}^{\dagger}(t)\ket{0},$
(65b) $\displaystyle H_{0}\ket{\psi_{L}(t)}$
$\displaystyle=-i(J_{0}+i\Omega)\ket{\psi_{L}(t)}-\sqrt{J_{0}}e^{-(J_{0}+i\Omega)t}C_{L}^{\dagger}(0)\ket{0}+\sqrt{J_{0}}C_{L}^{\dagger}(-t)\ket{0}.$
(65c)

Finally, we consider the action of the interaction Hamiltonian ($H_{I}$) on
each part of $\ket{\psi(t)}$

$\displaystyle H_{I}\ket{\psi_{q}(t)}$
$\displaystyle=\sqrt{J_{0}}e^{-(J_{0}+i\Omega)t}[C_{R}^{\dagger}(0)+C_{L}^{\dagger}(0)]\ket{0},$
(66a) $\displaystyle H_{I}\ket{\psi_{R}(t)}$
$\displaystyle=-iJ_{0}\Theta(0)e^{-(J_{0}+i\Omega)t}\ket{e},$ (66b)
$\displaystyle H_{I}\ket{\psi_{L}(t)}$
$\displaystyle=-iJ_{0}e^{-(J_{0}+i\Omega)t}(1-\Theta(0))\ket{e}.$ (66c)

From our calculations, we see that
$i\partial_{t}\ket{\psi_{q}(t)}=H_{0}\ket{\psi_{q}(t)}+H_{I}\ket{\phi_{R}(t)}+H_{I}\ket{\phi_{L}(t)}$.
Similarly,
$i\partial_{t}(\ket{\psi_{R}(r)}+\ket{\psi_{L}(t)})=H_{I}\ket{\psi_{q}(t)}+H_{0}(\ket{\psi_{R}(t)}+\ket{\psi_{L}(t)})$.
Consequently, the time evolved state given in (60-61) satisfies the
Schrödinger equation for $t>0$ and describes the spontaneous decay of an
initially excited qubit. As can be seen from this expression, the excited
qubit component of the state decays with a decay rate $J_{0}$ and the position
component consists of a superposition of two decaying exponentials travelling
away from the qubit, as we argued at the beginning of this section. This state
has an interesting feature: it is exactly a solution to the Schrödinger
equation, whereas in evaluation of (9b), approximations are performed. This is
due to the fact that all these approximations are exact in the limit
$J_{0}/\Omega\to 0$, and the time evolution of single qubit states has similar
behavior, i.e. follows a similar physical pattern, for all $J_{0}/\Omega$
values within the accuracy of the rotating-wave approximation. This property
will be discussed further in future work, where we plan to show that the cure
of causality violations in the RWA by the extension of energy integrals to
negative energies (as explained in [64] and had been employed by Fermi [62])
is not a trick, but rather has a physical basis.

##  Appendix B Waveguide QED with networks of quantum emitters

In this appendix, we extend the analysis beyond three qubits in a 1D
waveguide. We re-derive the results for quantum networks with various emitters
that can be modelled as identical two-level systems. We show that, in the end,
the generalization comes down to simply summing over all possible stationary
states, in analogy to the case of three qubits. This demonstrates the power of
the real-space approach. At the end of this section, we argue that this
approach can also be used for non-identical quantum emitters.

For the general network, the Hamiltonian can be written as
$H=H_{0}+H_{\textsc{i}}$, where $H_{0}$ contains the self energies and
$H_{\textsc{i}}$ is the interaction Hamiltonian that includes point-like
interactions. In real space, the Ansatz for such a system can be written as

$\begin{split}|E_{k}^{(i)}\rangle&=\sum_{\Xi\in\\{\rm
WG\\}}\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}x\phi_{\Xi}^{(i)}(x)C_{\Xi}^{\dagger}(x)\ket{0}+\sum_{s\in\\{\rm
qubits\\}}e_{s}^{(i)}\ket{e_{s}},\end{split}$ (67)

where, $\Xi$ represents all possible modes of all 1D waveguides inside the
network, $\phi_{\Xi}^{(i)}(x)$ includes the field amplitudes with the initial
condition $i$, $C_{\Xi}^{\dagger}(x)$ is the creation operator corresponding
to the mode $\Xi$ and $e_{s}^{(i)}$ is the excitation coefficient for the
qubit denoted by $s$. If the interactions are point interactions, then the
field amplitudes $\phi_{\Xi}^{(i)}(x)$ change only at the atomic positions
while behaving as free fields inside the waveguide.

Without loss of generality, we assume that the energy eigenstates are
normalized as $\langle
E_{k}^{(i)}|E_{p}^{(j)}\rangle=2\pi\delta(k-p)\delta_{ij}$. Here, there is one
aspect that requires further attention. When finding the energy eigenbasis, we
usually consider the scattering of light that is initially incident from one
side of the system. It is not always clear whether the (degenerate)
eigenstates corresponding to the same momentum value found via this method are
orthogonal. [35] doesn’t mention this for the linear chain of qubits, but one
can prove this using the transfer matrix properties. For the general case, we
keep in mind that this normalization is only valid in the limit
$J_{0}/\Omega\to 0$ such that the left and right moving particle subspace does
not mix. This ensures the orthogonality condition for the degenerate
eigenstates found by considering one-sided excitations.

###  B.1 Number of Markovian collective decay rates

In a network with $N$ qubits, the maximum number of Markovian collective decay
rates can be at most $N$. To illustrate this fact, we first consider a
generalization of Equation (6c) from [35]:

$\sum_{\Xi}\sqrt{J_{\Xi}}\phi_{\Xi}(x_{s})-\Delta_{k}e_{s}^{(i)}=0,$ (68)

where $x_{s}$ is the position of the qubit denoted by $s$ and $J_{\Xi}$ is the
coupling energy of the $\Xi$ mode to the qubits. This equation needs to be
satisfied at each atomic position, such that $|E_{k}^{(i)}\rangle$ is an
energy eigenstate.

Now, let us construct a coefficient matrix $A$, such that $Ax=b$, with $x$
including the scattering parameters and $b$ including the initial conditions.
The coefficient matrix $A$ includes the coefficient of the scattering
parameters in the generalized versions of the equations (6a-c) from [35]. The
scattering parameters, $x$, diverge for $\Delta_{k}$ values, for which $A$ is
singular. As a result, the characteristic equation for the poles of scattering
parameters can be given as $\det(A(\Delta_{k}))=0$. Since $A$ has only $N$
distinct $\Delta_{k}$ values, each at different rows according to (68), and
the phase is linearized in the Markovian limit; the characteristic polynomial
for the poles is of $N$th order for $\Delta_{k}$ and can have only $N$
distinct solutions $p^{\prime}$. Then, let us define $p$ equal to $p^{\prime}$
if $p^{\prime}$ is in the lower-half plane or on the real axis; or as
$p^{{}^{\prime}*}$ if $p^{\prime}$ is in the higher-half plane. In analogy
with the three qubits case, these solutions relate to the decay rates via a
$\Gamma=2ip$. In a future work, we will prove that the poles $p^{\prime}$ of
the scattering parameters are indeed in the LHP or on the real axis,
consequently $p=p^{\prime}$. This theorem requires a long discussion of the
causality principle in waveguide QED and is therefore left as a future work.

###  B.2 Time evolution

As before, we start by assuming that there are no BICs present in the system.
Then, the time evolution operator is

$U(t)=\sum_{i\in\\{\rm
I.C.\\}}\int_{0}^{\infty}\frac{\mathop{}\\!\mathrm{d}k}{2\pi}|E_{k}^{(i)}\rangle\langle
E_{k}^{(i)}|e^{-iE_{k}t}.$ (69)

Here, $i$ represents the initial conditions (I.C.), two for each 1D waveguide
present in the system. Then, the time evolution of any regular state is

$\displaystyle\ket{\psi(t)}$ $\displaystyle=\sum_{i\in\\{\rm
I.C.\\}}\int_{0}^{\infty}\frac{\mathop{}\\!\mathrm{d}k}{2\pi}|E_{k}^{(i)}\rangle\langle
E_{k}^{(i)}|\psi(0)\rangle e^{-iE_{k}t},$ (70a)
$\displaystyle=\int_{-\Omega}^{\infty}\frac{\mathop{}\\!\mathrm{d}\Delta_{k}}{2\pi}\ket{g(k)}\exp{-i\Delta_{k}t},$
(70b)
$\displaystyle\simeq\int_{-\infty}^{\infty}\frac{\mathop{}\\!\mathrm{d}\Delta_{k}}{2\pi}\ket{g(k)}\exp{-i\Delta_{k}t},$
(70c)
$\displaystyle=\sum_{p}\underset{\Delta_{k}=p}{\text{Res}}\left[\ket{g(k)}\exp{-i\Delta_{k}t}\right],$
(70d)

where in the second line, we substitute $k\to\Delta_{k}$ and define
$\ket{g(k)}=\sum_{i\in\\{\rm I.C.\\}}\langle
E_{k}^{(i)}|\psi(0)\rangle|E_{k}^{(i)}\rangle$. The rest is analogous to the
case of three qubits considered earlier, with the final step being meaningful
only in the Markovian limit as we have discussed in Section 3. Here, $p$ are
the lower half plane poles of $|E_{k}^{(i)}\rangle$ and $\langle
E_{k}^{(i)}|$, and consequently of scattering parameters. As in the case of
three qubits, the poles contain complete information of the collective decay
rates of the system.

As a result, the collective decay rates in a quantum network consisting of 1D
waveguides can be read-off from the scattering parameters. The collective
decay rates govern the time evolution of the overall system in the Markovian
limit. This intuitive approach to the time evolution in waveguide QED can be
employed as long as the coupling energy between the qubits and light is
smaller than the resonance frequency and the distance between the emitters is
$L\sim O(\Omega^{-1})$ such that any exponential in the characteristic
equation can be linearized:

$kL_{ij}=(\Delta_{k}+\Omega)L_{ij}\simeq\Omega L_{ij}=\theta_{ij},$ (71)

where $L_{ij}$ is the distance between the two emitter and $\theta_{ij}$ is
the corresponding phase acquired by light.

Networks with non-identical emitters and/or with non-radiative decay can also
be described in the Markovian limit, as long as the resonance energies of the
emitters differ only $\sim O(J_{0})$, with $J_{0}=\min\\{J_{\Xi}\\}$, so that
the phase linearization is valid. We will do this for a linear chain in
Appendix C. The non-Markovian limit is trivial to describe by calculating the
final step of time evolution as a numerical integral. Overall, adding more
qubits or waveguides to the system changes only the scattering eigenstates
$|E_{k}^{(i)}\rangle$ and not the theory itself. Since the scattering
parameters can be found via solving a linear set of equations, the real space
approach is easily scale-able to large systems.

###  B.3 Single-photon pulse scattering

In section 5.2, we discussed pulse scattering in the Markovian limit for a
linear chain of three qubits and found that the transmission and reflection
amplitudes are determined by the $k$-mode scattering eigenstates, inline with
the findings of [30]. The calculations in Appendix H present a direct proof
that this property is more general than the specific example of a linear
quantum emitter chain and are not specific to the Markovian limit. The
scattering by any “black box” system can be described by the external degrees
of freedom, such that the scattered light amplitudes vary via the scattering
parameters of the stationary states. In such a case, the asymptotic field
amplitudes for each $k$-mode can be written as

$|\psi_{s}(k,t\to\infty)|^{2}=|\psi(k,0)t_{s}(k)|^{2}\,,$ (72)

where $|\psi_{s}(k,t\to\infty)|^{2}$ is the output field amplitude
corresponding to the (normalized) scattering parameter $t_{s}(k)$ and
$|\psi(k,0)|^{2}$ is the input field amplitude.

We thus see that the real-space approach can be used to find collective decay
rates in complicated systems with multiple quantum emitters and multiple
waveguides. In doing so, the approach can be used to study subradiance and
superradiance, BICs, single-photon pulse scattering and time evolution of
system observable such as excitation probabilities.

Figure 16: The Bethe ansatz for a linear chain of non-identical qubits with
various energy separation $\Omega_{j}$ and coupling energies $J_{j}$.

##  Appendix C Waveguide QED with non-identical emitters: a general theory

In this appendix, we set up the theory for studying non-identical emitters. We
do this for a single waveguide coupled to $N$ qubits, but the derivations are
more general and can be employed for multi-waveguide systems.

An illustration of the linear chain consisting of $N$ non-identical qubits is
given in Fig. 16. The real-space Hamiltonian for this system is
$H=H_{\textsc{q}}+H_{\textsc{f}}+H_{\textsc{i}}$, where

$\displaystyle H_{\textsc{q}}={}$
$\displaystyle\sum_{j=1}^{N}\Omega_{j}\ket{e_{j}}\bra{e_{j}}$ (73)

is the free Hamiltonian of the qubits, where $\ket{e_{j}}$ is the excited
state of the $j^{th}$ qubit and $\Omega_{j}$ is the energy separation for
qubit $j$,

$\displaystyle\begin{split}H_{\textsc{f}}&=i\hbar
v_{g}\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}x\left(C_{L}^{\dagger}(x)\frac{\partial}{\partial
x}C_{L}(x)-C_{R}^{\dagger}(x)\frac{\partial}{\partial
x}C_{R}(x)\right)\end{split}$ (74)

is the free Hamiltonian of the field, where $v_{g}$ is the group velocity of
photons inside the waveguide and $C_{R/L}(x)$ are annihilation operators for
right/left moving photons, and

$\displaystyle\begin{split}H_{\textsc{i}}={}&\sum_{j=1}^{N}\sqrt{J_{j}}\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}x\delta(x-(j-1)L)\left((C_{R}^{\dagger}(x)+C_{L}^{\dagger}(x))\sigma_{j}+\text{H.c.}\right)\,,\end{split}$
(75)

is the interaction Hamiltonian between the systems, where $J_{j}$ is the
coupling energy between qubit $j$ and light. From now on, we set
$\hbar=v_{g}=1$ for algebraic simplicity, as usual.

Similar to the identical qubit case, we start by writing a Bethe Ansatz for a
general eigenstate

$\begin{split}\ket{E_{k}}&=\sum_{\Xi\in\\{WG\\}}\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}x\phi_{\Xi}(x)C_{\Xi}^{\dagger}(x)\ket{0}+\sum_{j=1}^{N}e_{j}\ket{e_{j}},\end{split}$
(76)

where $e_{j}$ is the excitation coefficient for the $j$th qubit and
$\phi_{R/L}(x)$ are field amplitudes for left/right moving photons. Note that
$e_{j}$ and $\phi_{R/L}(x)$ depend on $\theta$ implicitly. Owing to delta-
function interactions, the field becomes free for $x\neq(j-1)L$. Hence, we can
write an Ansatz for the field amplitudes $\phi_{R/L}(x)$ as:

$\displaystyle\phi_{R}(x)$ $\displaystyle=\begin{cases}t_{1}e^{ikx}&x<0,\\\
t_{j+1}e^{ik[x-jL]}&(j-1)L<x<jL,\\\ t_{N+1}e^{ik[x-NL]}&x>(N-1)L,\end{cases}$
(77a) $\displaystyle\phi_{L}(x)$
$\displaystyle=\begin{cases}r_{1}e^{-ikx}&x<0,\\\
r_{j+1}e^{-ik[x-jL]}&(j-1)L<x<jL,\\\
r_{N+1}e^{-ik[x-NL]}&x>(N-1)L,\end{cases}$ (77b)

where $t_{1}$, $r_{1}$, $t_{N+1}$ and $r_{N+1}$ are picked for different
initial conditions. For a photon incident from the far left ($\ket{E_{k}}$):
$t_{1}=0$, $r_{N+1}=0$ and for non-radiating conditions (for BIC,
$\ket{D_{i}}$): $t_{1}=r_{1}=t_{N+1}=r_{N+1}=0$. For $\ket{E_{-k}}$, one can
mirror the state $\ket{E_{k}}$ w.r.t. the center of the linear chain, as in
Section 3. For now, we suppress the implicit $k$-dependence for brevity and
discuss the Markovian limit and corresponding additional conditions. The non-
Markovian regime calculations can be performed analogous to Section 6.

Having set an Ansatz for the energy eigenstates, we apply the condition
$H\ket{E_{k}}=E_{k}\ket{E_{k}}$, with $E_{k}=|k|$, to obtain the equation of
motion for the scattering parameters:

$\displaystyle t_{j+1}e^{-i\theta}-t_{j}+i\sqrt{J_{j}}e_{j}$
$\displaystyle=0,$ (78a) $\displaystyle
r_{j+1}e^{i\theta}-r_{j}-i\sqrt{J_{j}}e_{j}$ $\displaystyle=0,$ (78b)
$\displaystyle\sqrt{J_{j}}(t_{j}+r_{j})-(\Delta_{k}+\delta_{j})e_{j}$
$\displaystyle=0,$ (78c)

where we linearize the propagation phase as usual $\theta\simeq\Omega_{1}L$,
and define $\Delta_{k}=E_{k}-\Omega_{1}$ and $j=1,...,N$. Moreover, for $j\neq
1$, we define $\delta_{j}=\Omega_{1}-\Omega_{j}$ such that
$E_{k}-\Omega_{j}=\Delta_{k}+\delta_{j}$. For the case of non-identical
emitters, the Markovian limit has an additional constraint, apart from
$J_{j}/\Omega_{j}\ll 1$. So that the interacting frequencies are confined
within $|\Delta_{k}|\leq O(J_{i})$, the detuning of the qubit separation
frequencies should satisfy $|\delta_{j}|\sim O(J_{0})$. Consequently, the
linearization of $\theta$ around $\Omega_{j}$ are all equivalent. If this is
not satisfied, the linearization assumption of the phase is no longer valid.

There are two possible methods that could be used to solve (78):

  1. 1.

By writing the equations as a matrix equation with $3N$ unknowns and solving
this linear system. We used this method in Section 3, and it is reasonable to
do so here as long as the system is small. But for large systems, this method
requires unnecessarily high computation and is therefore inefficient.

  2. 2.

By using the transfer matrix method and obtaining a recursive algorithm to
find each scattering parameter. This method preserves the polynomial shape of
the scattering parameters and leads to the fully analytical results that we
desire. We will use this method to find the scattering parameters.

There is also an additional method considered in [35]. We could use the
transfer matrix method and draw parallels to one-dimensional photonic crystals
[74]. But apart from the fact that this method can only be used for identical
emitters, even in the identical qubit case, the numerator and denominator of
scattering parameters are no longer polynomial in $\Delta_{k}$, hence the
poles of the system cannot be obtained easily. As we have seen so far, the
complex analysis plays an important role in time evolution of single-photon
states. Therefore, this method is also not preferred for obtaining fully
analytical results.

To use the transfer matrix method, we shall first reshape (78) and eliminate
the $e_{j}$ degree of freedom:

$\displaystyle t_{j+1}e^{-i\theta}$
$\displaystyle=\left(1-i\frac{J_{j}}{\Delta_{k}+\delta_{j}}\right)t_{j}-i\frac{J_{j}}{\Delta_{k}+\delta_{j}}r_{j},$
(79a) $\displaystyle r_{j+1}e^{i\theta}$
$\displaystyle=i\frac{J_{j}}{\Delta_{k}+\delta_{j}}t_{j}+\left(1+i\frac{J_{j}}{\Delta_{k}+\delta_{j}}\right)t_{j}.$
(79b)

Writing this in matrix form and re-arranging terms, we obtain the matrix
equation

$\displaystyle\begin{split}\begin{pmatrix}t_{j}\\\
r_{j}\end{pmatrix}={}&\begin{pmatrix}1+i\frac{J_{j}}{\Delta_{k}+\delta_{j}}&i\frac{J_{j}}{\Delta_{k}+\delta_{j}}\\\
-i\frac{J_{j}}{\Delta_{k}+\delta_{j}}&1-i\frac{J_{j}}{\Delta_{k}+\delta_{j}}\end{pmatrix}\begin{pmatrix}e^{-i\theta}&0\\\
0&e^{i\theta}\end{pmatrix}\begin{pmatrix}t_{j+1}\\\
r_{j+1}\end{pmatrix}.\end{split}$ (80)

Then, we find the transfer matrix for a unit cell, which contains the
scattering from the qubit $j$ and a propagation phase $\theta$, as

$\begin{split}T_{j}&=\begin{pmatrix}1+i\frac{J_{j}}{\Delta_{k}+\delta_{j}}&i\frac{J_{j}}{\Delta_{k}+\delta_{j}}\\\
-i\frac{J_{j}}{\Delta_{k}+\delta_{j}}&1-i\frac{J_{j}}{\Delta_{k}+\delta_{j}}\end{pmatrix}\begin{pmatrix}e^{-i\theta}&0\\\
0&e^{i\theta}\end{pmatrix}\implies\begin{pmatrix}t_{j}\\\
r_{j}\end{pmatrix}=T_{j}\begin{pmatrix}t_{j+1}\\\
r_{j+1}\end{pmatrix}.\end{split}$ (81)

By recursively relating them, we can eliminate the internal scattering
parameters and obtain the relation between the external fields for the initial
condition, where the photon is incident from far left, as

$\begin{pmatrix}1\\\ r_{1}\end{pmatrix}=T\begin{pmatrix}t_{N+1}\\\
0\end{pmatrix},$ (82)

where $T=T_{1}T_{2}...T_{N}=\prod_{j=1}^{N}T_{j}$. From this equation, we can
find the external transmission and reflection coefficients as

$t_{N+1}=1/(T)_{11},\quad r_{1}=(T)_{21}/(T)_{11}.$ (83)

Here, $(A)_{ij}$ represents the $ij$th element of the matrix $A$. Then, the
rest of the transmission and reflection coefficients can be found via

$\begin{pmatrix}t_{j}\\\ r_{j}\end{pmatrix}=S_{j}\begin{pmatrix}1\\\
r_{1}\end{pmatrix},$ (84)

where we define $S_{j}=T_{j-1}^{-1}...T_{2}^{-1}T_{1}^{-1}$. Finally, the
excitation coefficients can be obtained from (78) as

$e_{j}=\sqrt{J_{j}}\frac{(t_{j}+r_{j})}{\Delta_{k}+\delta_{j}}.$ (85)

As an example, let us consider the two qubit case and find the collective
decay rates. A perturbative result for the decay rates has been obtained in
[53], here we find the complete expressions. Since the collective decay rates
can be read-off from any scattering parameter, let us find the reflection
coefficient $r_{1}$ by employing the transfer matrix approach:

$r_{1}=-\frac{i\left(J_{1}(\delta_{2}+\Delta_{k}+iJ_{2})+J_{2}e^{2i\theta}(\Delta_{k}-iJ_{1})\right)}{(\Delta_{k}+iJ_{1})(\delta_{2}+\Delta_{k}+iJ_{2})+J_{1}J_{2}e^{2i\theta}}.$
(86)

Then, the two collective decay rates are

$\displaystyle\Gamma_{1/2}$
$\displaystyle=\left(J_{1}+J_{2}-i\delta_{2}\pm\Delta\right),$ (87)

where $\Delta=\sqrt{(J_{1}+i\delta_{2}-J_{2})^{2}+4J_{1}J_{2}e^{2i\theta}}$.
For $J_{1}=J_{2}$ and $\theta=\pi$, expanding around small $\delta$ leads to
the perturbative results obtained in [53]. Furthermore, we emphasize that for
a photon incident from far right, exchanging $J_{1}\iff J_{2}$, setting
$\delta_{2}\to-\delta_{2}$ results in the same expressions for the collective
decay rates up to a global imaginary shift of $i\delta_{2}$, which is simply
due to the redefinition of $\Delta_{k}$ and is just a re-normalization of
energy levels. This shows clearly that collective decay rates are indeed
independent of initial conditions, as expected.

Moreover, the two-qubit system with identical emitters becomes transparent for
a Fano minimum ([75]) such that $r_{1}=0$, where the transmission becomes
unity [35]. By the same logic, for the general two-qubit system, we can find
the $\Delta_{k}$ value corresponding to the zero reflection as

$\Delta_{k}=-\frac{J_{1}\left(\delta_{2}-iJ_{2}e^{2i\theta}+iJ_{2}\right)}{J_{1}+J_{2}e^{2i\theta}}.$
(88)

However, this frequency detuning, for which $r_{1}=0$, is complex, whereas
only real $\Delta_{k}$ values are physical. This means that the expression in
(88) is not the Fano dip and the reflected pulse intensity may not be zero for
the general case. Consequently, the two-qubit system is only transparent if
the qubits are identical. In fact, the minimum value of the reflection
intensity takes a very complicated form, which is different than the real part
of (88). The transmission and reflection intensities are illustrated in Fig.
17 for a specific example. As apparent from the figure, the intensities show
Fano type line-shapes and the system no longer becomes completely transparent
for the Fano minimum. We emphasize that (88) reduces to the result found in
[35] for $J_{1}=J_{2}=J_{0}$ and $\delta_{2}=0$, which is
$\Delta_{k}=-J_{0}\tan(\theta)$. The two-qubit system becomes completely
transparent for this special case.

Figure 17: The transmission and reflection intensities, $|t_{2}|^{2}$ and
$|r_{1}|^{2}$ respectively, for a photon incident from far left to a two qubit
system with $J_{2}/J_{1}=2$, $\delta_{2}=0.3J_{1}$ and $\theta=0.85\pi$. The
intensities show Fano type line-shapes and the system no longer becomes
completely transparent for the Fano minimum.

We note that while we have not considered the non-radiative decay, the non-
radiative modes can be modelled as waveguides (as long as the dispersion
relation can be linearized) and interactions between the qubits can be
described via this formalism by using multiple waveguides, whose only effect
is to increase the dimensions of the transfer matrices.

##  Appendix D Parity (mirror) symmetry

We define the general parity operator, $\hat{P}$, that acts on the combined
Hilbert space of the single-excitation states. This operator simply mirrors
the state w.r.t. the center of the multi-qubit system. When the multi-qubit
system is centered at $x=0$, its action on creation operators is defined as

$\displaystyle\hat{P}\sigma_{j}^{\dagger}\hat{P}^{\dagger}=\sigma_{N-j+1}^{\dagger}$
$\displaystyle\implies\hat{P}\ket{e_{j}}=\ket{e_{N-j+1}},$ (89a)
$\displaystyle\hat{P}C_{R/L}^{\dagger}(x)\hat{P}^{\dagger}=C_{L/R}^{\dagger}(-x)$
$\displaystyle\implies\hat{P}\ket{x}=\ket{-x}.$ (89b)

Then, the symmetry states are defined as the eigenvectors of $\hat{P}$ that
have $+1$ eigenvalues, whereas anti-symmetric states correspond to $-1$
eigenvalues. For example, for odd $N$, the state where the initial qubit is
excited is a symmetric state since
$\hat{P}\ket{e_{(N+1)/2}}=\ket{e_{(N+1)/2}}$. Note that the general parity
operator leads to the natural definition of even (symmetric) and odd (anti-
symmetric) basis operators, as discussed in the literature (for example [6])

$\displaystyle C_{e}^{\dagger}(x)$
$\displaystyle=\frac{1}{\sqrt{2}}\left(C_{R}^{\dagger}(x)+C_{L}^{\dagger}(-x)\right),$
(90a) $\displaystyle C_{o}^{\dagger}(x)$
$\displaystyle=\frac{1}{\sqrt{2}}\left(C_{R}^{\dagger}(x)-C_{L}^{\dagger}(-x)\right),$
(90b)

such that $\hat{P}C_{e/o}^{\dagger}(x)\hat{P}^{\dagger}=\pm
C_{e/o}^{\dagger}(x)$. Since a single qubit couples to light symmetrically,
odd (e.g. anti-symmetric) states cannot excite single qubits, which will be
discussed in details in Section 5.2.

Let us now use the parity operator to find the energy eigenstates. First, it
is straightforward to show that the parity operator commutes with the
Hamiltonian $[H,\hat{P}]=0$. Then, let $\ket{E_{k}}$ be an eigenstate of the
Hamiltonian but not an eigenstate of the parity operator. (Any scattering
state, where the photon is initially incident from one side, fits this
definition.) Then, $\hat{P}\ket{E_{k}}$ is also an eigenstate of the
Hamiltonian since:

$H(\hat{P}\ket{E_{k}})=\hat{P}H\ket{E_{k}}=\hat{P}E_{k}\ket{E_{k}}=E_{k}(\hat{P}\ket{E_{k}}).$
(91)

Consequently, $\ket{E_{-k}}=\hat{P}\ket{E_{k}}$. Throughout the paper, we
refer to this procedure simply by stating that $\ket{E_{-k}}$ can be found via
symmetry considerations.

##  Appendix E An example for how residues contain information on BIC

For the bright state $\ket{B}$, the system completely couples with the light
and the residues corresponding to subradiant poles are zero. The reason why
the subradiant residues are non-zero for cases, where the initial state has
non-zero overlap with subradiant states, can be explained as follows. Consider
$|\theta-n\pi|\simeq\delta$. In this case, the subradiant state couples to the
light slightly and decay is described solely by scattering eigenstates, since
there are no BICs. Any regular state $\ket{\psi}$ can be written as

$\ket{\psi}=\sum_{i}\alpha_{i}\ket{D_{i}}+\beta\ket{B}.$ (92)

For any $\delta>0$, no matter how small, the residue has a finite value for
$t=0$, which is the complex coefficients corresponding to states $\ket{D_{i}}$
and $\ket{B}$ (which is, for example, $\alpha_{1}$ for the coefficient of
$\ket{D_{1}}$). Taking the limit $\delta\to 0$ then results in $\alpha_{i}$
for $\ket{D_{i}}$, since $\alpha_{i}$ is obtained for each value of
$\delta>0$. The exponential term corresponding to dark states become unity,
since the subradiant poles become zero.

To illustrate how a state approaches the dark state as $\theta\rightarrow\pi$,
let us compute the overlap $|\innerproduct{D_{2}}{\psi(t)}|$, where
$\ket{\psi(0)}=\ket{e_{0}}$. This overlap is important, because it is a
continuous function of $\theta\neq\pi$ for any $t$ for the time evolution
operator presented in (8), whereas it becomes discontinuous when
$\theta\neq\pi$. In this case, the time evolution given in (10) should be
used. Now, of course, the immediate question is whether this discontinuity is
a physical one, that is whether we can remove it by using (8) asymptotically
at $\theta\rightarrow\pi$ and omit the definition of (10). The overlap is
discontinuous at $\theta=pi$, since the scattering part, which constitutes
(8), does not contribute to this overlap exactly at $\theta=\pi$, whereas they
do for any other value $\theta\neq\pi$.

Defining $\delta=\pi-\theta$, we compare the overlap as a function of time $t$
(with units $J_{0}^{-1}$) for various $\delta$. As $\delta$ smaller, this
illustrates taking the limit $\theta\rightarrow\pi$. Since $\ket{D_{2}}$ does
not decay as $\delta\to 0$, the coefficient $|\innerproduct{D_{2}}{\psi(t)}|$
is also expected to be non-zero and non-decaying in this limit. For
$\delta=0$, this term can only be obtained from the BIC contribution in (10),
where the residue obtained from the scattering states would not contribute to
this coefficient since $\ket{D_{2}}$ is a dark state. Fig. 18 shows that the
residue obtained from the scattering parameter is non-zero in the limit
$\delta\rightarrow 0$. Thus, the time evolution in the limit $\delta\to 0$ can
be described by taking the limit of the residues found using only the
scattering eigenstates.

Figure 18: The overlap $|\innerproduct{D_{2}}{\psi(t)}|$ where
$\ket{\psi(0)}=\ket{e_{0}}$ for various $\delta=\pi-\theta$. As $\delta$ gets
smaller, the system approaches the regime where BICs appear. Taking the limit
$\delta\to 0$, the information on BIC is stored in the residues of the poles.

This shows that the proposed strategy yields the correct results, because the
limit of the residue converges to the coefficient that would have been
obtained via the complete time evolution operator with the BIC contribution.

##  Appendix F An alternative method for finding collective decay rates

In this appendix, we show an alternative method for finding collective decay
rates for the special case of a linear chain of $N$ qubits. This method is
less efficient than simply finding the poles of the scattering parameters, but
it provides further insight into the physics of the system. Specifically, this
method shows that each collective decay rate corresponds to a specific
coherent excitation of the system. This means that one can always prepare
initial states that excited only one of the $N$ collective decay rates.

Let us start with the following Ansatz for an initially excited system
following the approach by [26]

$\ket{\psi(t)}=\sum_{j=1}^{N}\alpha_{j}(t)e^{-i\Omega
t}\ket{e_{j}}+\ket{\chi(t)},$ (93)

where $\ket{\chi(t)}$ includes the photon contribution with the initial
condition $\ket{\chi(0)}=0$. Then, the relation between $\alpha_{j}(t)$ is
found in [26] and can be given in terms of the variables defined in this paper
as

$\dot{\alpha}_{j}(t)=-\sum_{k=1}^{N}J_{0}e^{i\theta|j-k|}\alpha_{k}(t-L|j-k|),$
(94)

where $\dot{\alpha}_{j}(t)$ denotes the time derivative of $\alpha_{j}(t)$.
[26] uses a more general form of this equation to calculate observables (upon
pulse incidence). Here, we will take another approach and show that this
equation can yield the collective decay rates. We also note that here we are
not interested in the time evolution of $\ket{\chi(t)}$.

Since we are interested in the collective decay rates, we take the Markovian
limit such that $\alpha_{k}(t-L|j-k|)\simeq\alpha_{k}(t)$ since $L\sim
O(\Omega^{-1})$ and $t\sim O(J_{0}^{-1})$. Consequently, (94) can be written
as a matrix equation

$\dot{x}(t)=-Jx(t),$ (95)

where $x(t)=[\alpha_{1}(t),\ldots,\alpha_{N}(t)]$ includes the qubit
excitation coefficients and $(J)_{jk}=J_{0}e^{i\theta|j-k|}$ is the collective
coupling matrix. The solution to this equation is trivial (as long as $J$ is
non-singular, which it is for $\theta\neq n\pi$; for $\theta=n\pi$, the
existence of BIC changes the time evolution as discussed in Section 3.2) and
can be given as

$x(t)=\sum_{l=1}^{N}\beta_{l}e^{-0.5\Gamma_{l}t}\xi_{l},$ (96)

where $\xi_{l}$ and $\Gamma_{l}/2$ are eigenvectors and eigenvalues of the
matrix $J$, and $\beta_{l}$ are some complex constants that can be determined
by the initial conditions 888Here, we assumed, for simplicity, that $J$ has a
non-degenerate spectrum, which does not affect any of the arguments that
follows.. This expression ties back to (12c), where $\Gamma_{l}$ are indeed
the collective decay rates.

This derivation reveals an important property of collective decay rates: they
have one-to-one correspondence with states ($\ket{\xi_{l}}$) in the qubit
subspace. It is important to note that these states are not guaranteed to be
orthogonal, although they are distinct and span the qubit subspace. While [26]
discusses the eigenvalues and eigenvectors of the coupling matrix and
associates, in passing, eigenvalues with collective decay rates when
considering a pulse scattering problem, here we proved this relationship,
which is only accurate for the Markovian limit.

The eigenstates of the matrix are the basis states of the $N$ distinct decay
modes. Consequently, each collective decay rate corresponds to a physical
decay mode. This leads to the following phase space picture: the initial
coherent excitation of the qubits can be written in terms of a linear
combination of decay mode basis states. Consequently, a decay mode can only be
accessed if the overlap is nonzero999It is important to clarify that the
eigenstates are not necessarily orthogonal, so the overlap is not taken as
orthogonal projections, but according to the angles between eigenstates. This
leads to the fact that a system prepared in a certain eigenmode has non-zero
probability to be observed at another mode. Nonetheless, this is not an
interference effect of different modes, since any observed quantity decays
only with the decay rate of the initially prepared mode.. Hence, there is
always a specific coherent excitation (i.e. one that overlaps perfectly with
the corresponding eigenvector) of qubits that can excite a single decay mode
only. We discuss this property in the next section in the context of pulse
shaping by engineering collective decay rates. Moreover, this property also
explains how decay rates can signal the existence of BIC and why the
dimensionality ($N-1$) of BIC is linked to the number ($N-1$) of zero
collective decay rates. Since for $\theta=n\pi$, the $J$ matrix has $N-1$ zero
eigenvalues, the corresponding subspace has dimensionality $N-1$ and can be
constructed with orthogonal basis states. We also infer that subradiant states
become BIC in a continuous manner as $\theta$ approaches $n\pi$.

On another note, combining the one-to-one correspondence of decay modes and
$\ket{\xi_{l}}$ with the symmetric and anti-symmetric collective decay rates
conjecture, we realize that the states $\ket{\xi_{l}}$ have either even or odd
parity such that $P\ket{\xi_{l}}=\pm\ket{\xi_{l}}$. The even (odd) parity
states correspond to symmetric (anti-symmetric) states. The implication is
easy to prove via proof by contradiction. Assume $\ket{\xi_{1}}$ has both
symmetric and anti-symmetric parts. Then, it can be decomposed into both parts
and hence can excite certain symmetric and anti-symmetric modes (that have
non-zero overlap with the symmetric/anti-symmetric part of $\ket{\xi_{1}}$),
which contradicts the conjecture. The fact that the decay mode states,
$\ket{\xi_{l}}$, are either symmetric or anti-symmetric, and not a mixture of
two, is intriguing. For now, we do not have a conclusive proof for this,
although we believe that the highly special shape of the $J$ matrix might be
the first step towards understanding this phenomenon.

We emphasize that finding the collective decay rates via this method is
inefficient, since it requires diagonalization of a $N\times N$ matrix. Using
the transfer matrix method is efficient, since it eliminates the internal
degrees of freedom and deals with only $2\times 2$ matrices, as shown in
Appendix B. For example, using the transfer matrix method, the collective
decay rates for $N=30$ can be found almost instantly, whereas it is nearly
impossible to diagonalize the coupling matrix. This phenomenon illustrates
clearly why real-space approach outperforms the existing methods [26] by a
large margin. We also emphasize that this approach works only in the Markovian
limit. To find non-Markovian collective decay rates, one needs to use the
machinery of the real-space approach.

##  Appendix G Pulse scattering examples for the Markovian regime

Here we consider two more distinct pulse shapes: 1) Decaying Exponential, i.e.
the shape of a photon emitted from a two-level quantum emitter [72] and 2)
Rising Exponential, which is known to give maximum excitation for a single
qubit [31] (see Appendix A).

####  G.0.1 Decaying exponential

A decaying exponential pulse corresponds to a regular state, therefore we can
obtain analytical expressions for each atomic excitation probability, as well
as the transmitted and reflected pulse shapes. For a (resonant) decaying
exponential incident from the left, $f(x)$ in (29) takes the form

$f(x)=\sqrt{2\xi}\mathrm{e}^{\xi x}\Theta(-x).$ (97)

Then, the state of the pulse at time $t$ is

$\begin{split}&\ket{S(t)}=\int_{-\infty}^{\infty}\frac{\mathop{}\\!\mathrm{d}\Delta_{k}}{2\pi}\frac{\sqrt{2\xi}e^{-i\Delta_{k}(t-x_{0})}}{\xi-i\Delta_{k}}\ket{E_{k}}=\sum_{p^{\prime}}\underset{\Delta_{k}=p^{\prime}}{\text{Res}}\left[\frac{\sqrt{2\xi}e^{-i\Delta_{k}(t-x_{0})}}{\xi-i\Delta_{k}}\ket{E_{k}}\right],\end{split}$
(98)

where $p^{\prime}$ includes, in addition to three poles of the system, an
additional pole ($-i\xi$) introduced by the decaying exponential. It is
important to note that this expression has no UHP poles, hence the interaction
between the photon and the system occurs only for $t\leq x_{0}$, which is a
consequence that photons travel with a speed $v_{g}=1$. Consequently, the
causality principle manifests itself in this scattering problem by the absence
of poles in the upper-half plane. The derivation for the single photon case
([31]) can be found in Appendix A. The emitted photon probability density
$\mathcal{P}(x,t)$ can also be found analytically, but we don’t show it here
as it doesn’t offer any more insight into the formalism.

Figure 19: Single-photon scattering in a system of three qubits in the
Markovian limit. The excitation probability of each atom during interaction
with a single photon pulse shaped as a decaying exponential centred at
$x_{0}=10/J_{0}$. For each value of $\theta$, the parameter $\xi$ was
optimized to maximise the excitation probability for the left atom; the
optimal values are a) $\xi=0.94J_{0}$, b) $\xi=1.12J_{0}$, c) $\xi=0.73J_{0}$
and d) $\xi=3J_{0}$. In d), the linearization assumption causes all three
excitation probabilities to be the same, which is accurate as long as $L\sim
O(\Omega^{-1})$ and $J_{0}/\Omega\ll 1$.

We do, however, study the atom excitation probabilities. In particular, we
optimize the decaying exponential and $\theta$ to maximize or minimize the
(maximum) excitation probabilities.

The maximum excitation probability for the left atom is
$P^{\mathrm{max}}_{-1}(t)\simeq 0.454$, and occurs when the system is
parameterized by $\theta=\pi/2$, and the pulse is parameterized by
$\xi=0.73J_{0}$ and $t=x_{0}+1.43/J_{0}$.

The fact that the maximum excitation is obtained for $\theta=\pi/2$ can be
explained as follows. The imaginary parts of the poles correspond to decaying
exponential. If the decay is large, then the coupling between the light and
the atom is strong and the atom can be excited by the incoming pulse more
easily. Nonetheless, an excited atom can decay faster through high decay rate
modes, leading to a trade-off. Therefore, from an heuristic point of view, it
is expected that the highest probability of excitation is achieved for the
case where the relevant decay modes have similar decay rates, which is when
$\theta=\pi/2$.

In the other extreme, the minimum (maximum) excitation probability is
$P^{\mathrm{min}}_{-1}(t)=\frac{2}{3e^{2}}\simeq 0.09$, and occurs when the
system is parameterized by $\theta=\pi$, and the pulse is parameterized by
$\xi=3J_{0}$ and $t=x_{0}+(3J_{0})^{-1}$. This can be understood as follows.
When $\theta=\pi$, out of three decay rates, two of them become zero, leading
to dark states that do not couple to the incoming pulse. Therefore, the
collective system behaves as a two-level system between the ground state and
the bright state $\ket{B}$. This explains the excitation probability
$\frac{2}{3e^{2}}$ of individual atoms, since the overall qubit excitation of
$\frac{2}{e^{2}}$ (See Appendix A) is divided equally between them.

We can also optimize the system when a subset of the parameters are fixed.
Fig. 19, show excitation probabilities for which $P_{-1}(t)$ is maximized for
different $\theta$.

####  G.0.2 Rising exponential

A rising exponential photon corresponds to a regular state. Therefore, the
discussion on rising exponential pulse is identical to the decaying
exponential. For brevity, we report directly the result. The rising
exponential pulse results in a maximum excitation of $\simeq 0.6808$ for the
first atom for $\xi=0.97J_{0}$, $\tau=0$ (with $x_{0}\to\infty$) and
$\theta=\pi/2$. This shows that the rising exponential is a more suitable one-
sided pulse to achieve a maximum excitation probability in comparison to
decaying exponential and Gaussian pulses.

##  Appendix H Derivation of scattering pulse shape

In this appendix, we carry out the derivations of $S(x<0,2x_{0})$ and
$S(x>0,2x_{0})$ from the general formula of $S(x,t)$ given in (31). Here, the
time is chosen as $t=2x_{0}$ so that the qubits are practically in the ground
state. We shall start by dividing the integral into two parts

$\displaystyle
S(x,t)=\frac{1}{2\pi}\left(\int_{0}^{\infty}\mathop{}\\!\mathrm{d}k\innerproduct{x}{E_{k}}\innerproduct{E_{k}}{S(0)}e^{-iE_{k}t}+\int_{0}^{\infty}\mathop{}\\!\mathrm{d}k\innerproduct{x}{E_{-k}}\innerproduct{E_{-k}}{S(0)}e^{-iE_{k}t}\right),$
(99)

where by $E_{-k}$ we mean the scattering energy-eigenstate for a photon far
from right. For the next step, we find

$\displaystyle\innerproduct{E_{k}}{S(0)}$
$\displaystyle\simeq\int_{-\infty}^{\infty}f(x+x_{0})e^{-i(k-k_{0})x}\mathop{}\\!\mathrm{d}x=\sqrt{2\pi}\tilde{f}(k-k_{0})e^{i(k-k_{0})x_{0}},$
(100a) $\displaystyle\innerproduct{E_{-k}}{S(0)}$ $\displaystyle\simeq 0,$
(100b)

where we make use of the fact that $f(x+x_{0})\simeq 0$ for $x>0$ and
$[C_{R}^{\dagger}(x),C_{L}(x)]=0$.

Now, we set $t=2x_{0}$ and assume $x>0$ to find the transmitted signal, then

$\displaystyle S(x,2x_{0})$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}k\left(t_{3}e^{ikx}\right)\tilde{f}(k-k_{0})e^{i(k-k_{0})x_{0}}e^{-2iE_{k}x_{0}},$
(101a)
$\displaystyle=\frac{e^{-ik_{0}x_{0}}}{\sqrt{2\pi}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}k\left(t_{3}e^{ikx}\right)\tilde{f}(k-k_{0})e^{-ikx_{0}}.$
(101b)

Here, we recall that $E_{k}=k$ for $k>0$ and $\tilde{f}(k-k_{0})\simeq 0$ for
$k<0$. Then, we have

$\displaystyle S(x,2x_{0})$
$\displaystyle=\frac{e^{-ik_{0}x_{0}}}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}k\,t_{3}\tilde{f}(k-k_{0})e^{ik(x-x_{0})},$
(102)

where we realize that
$\tilde{S}_{+}(k,2x_{0})=e^{-i(k_{0}+k)x_{0}}\,t_{3}\tilde{f}(k-k_{0})$. Now,
as $|\tilde{S}(k,0)|^{2}=|\tilde{f}(k-k_{0})|^{2}$, we have that
$|\tilde{S}_{+}(k,2x_{0})|^{2}=|t_{3}\tilde{f}(k-k_{0})|^{2}=|t_{3}\tilde{S}(k,0)|^{2}$.

For the reflected signal, we assume $x<0$. Then, we find

$\displaystyle S(x,2x_{0})$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}k\left(e^{ikx}+r_{1}e^{-ikx}\right)\tilde{f}(k-k_{0})e^{i(k-k_{0})x_{0}}e^{-2iE_{k}x_{0}},$
(103a)
$\displaystyle=\frac{e^{-ik_{0}x_{0}}}{\sqrt{2\pi}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}k\left(r_{1}e^{-ikx}\right)\tilde{f}(k-k_{0})e^{-ikx_{0}},$
(103b)

as the first term of the integral (which includes $e^{ik(x-x_{0})}$)
corresponds to a function localized at $x=x_{0}$ and is therefore $\simeq 0$
for $x<0$ and $E_{k}=k$ for $k>0$. We first realize that
$\tilde{f}(k-k_{0})\simeq 0$ for $k<0$ and expand the order of integration to
$-\infty$. Then, we change $k\to-k$

$\displaystyle S(x,2x_{0})$
$\displaystyle=\frac{e^{-ik_{0}x_{0}}}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}k\,r_{1}\tilde{f}(-k-k_{0})e^{ik(x+x_{0})},$
(104)

where we realize that
$\tilde{S}_{-}(k,2x_{0})=e^{i(k-k_{0})x_{0}}r_{1}\tilde{f}(-k-k_{0})$. Then,
$|\tilde{S}_{-}(-k,2x_{0})|^{2}=|r_{1}\tilde{f}(k-k_{0})|^{2}=|r_{1}\tilde{S}(k,0)|^{2}$.
The narrow band pulses scatter indeed corresponding to the stationary
eigenstate scattering coefficients $t_{3}$ and $r_{1}$.

##  Appendix I The time-evolved state $\ket{\psi(t)}$ satisfies the
Schrödinger equation ($N=2$ Qubits)

In this appendix, we show that the time evolved state given in (45) satisfies
the Schrödinger equation

$i\frac{\partial}{\partial t}\ket{\psi(t)}=(H_{0}+H_{I})\ket{\psi(t)},$ (105)

where the Hamiltonian for the $2$ qubit system is

$\displaystyle H_{0}$
$\displaystyle=\Omega\sum_{j=\\{-1,1\\}}\ket{e_{i}}\bra{e_{i}}+i\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}x\left(C_{L}^{\dagger}(x)\frac{\partial}{\partial
x}C_{L}(x)-C_{R}^{\dagger}(x)\frac{\partial}{\partial x}C_{R}(x)\right),$
(106a) $\displaystyle H_{I}$
$\displaystyle=\sqrt{J_{0}}\sum_{j=\\{-1,1\\}}\left(\sigma_{j}^{\dagger}[C_{L}(jL/2)+C_{R}(jL/2)]+\sigma_{j}[C_{L}^{\dagger}(jL/2)+C_{R}^{\dagger}(jL/2)]\right),$
(106b)

where we pick $x=0$ as the center of the $2$ qubit system.

We re-write (45) once more while scaling $\Gamma_{1}\to 2\Gamma_{1}$ for
convenience:

$\begin{split}\ket{\psi(t)}=\frac{1}{\sqrt{2}}\Bigg{(}e^{-\Gamma_{1}t-i\Omega
t}(\ket{e_{1}}+\ket{e_{2}})&-i\int_{0}^{t}\mathop{}\\!\mathrm{d}x\frac{\Gamma_{1}}{\sqrt{J_{0}}}e^{-(\Gamma_{1}+i\Omega)(t-x+L/2)}C_{R}^{\dagger}(x)\ket{0}\\\
&-i\int_{-t}^{0}\mathop{}\\!\mathrm{d}x\frac{\Gamma_{1}}{\sqrt{J_{0}}}e^{-(\Gamma_{1}+i\Omega)(t+x+L/2)}C_{L}^{\dagger}(x)\ket{0}\Bigg{)},\end{split}$
(107)

Our calculations are performed in the Markovian limit, where $\Omega
L=\theta$, $J_{0}L\simeq 0$ and $C_{R/L}^{\dagger}(\pm L/2)\simeq e^{\mp
i\theta/2}C_{R/L}^{\dagger}(0)$. The latter can be proven using the definition
of $C_{R/L}^{\dagger}(x)$ as a Fourier Transform of momentum mode creation
operator $a_{R}^{\dagger}(k)$. This is consistent with the time evolution in
(9b), where the Markovian limit is assumed implicitly to obtain a polynomial
characteristic equation for the collective decay rates.

First, we divide
$\ket{\psi(t)}=\ket{\psi_{q}(t)}+\ket{\psi_{R}(t)}+\ket{\psi_{L}(t)}$ into
three parts as

$\displaystyle\ket{\psi_{q}(t)}$
$\displaystyle=\frac{1}{\sqrt{2}}e^{-\Gamma_{1}t-i\Omega
t}(\ket{e_{1}}+\ket{e_{2}}),$ (108a) $\displaystyle\ket{\psi_{R}(t)}$
$\displaystyle=\frac{-i}{\sqrt{2}}\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}x\frac{\Gamma_{1}}{\sqrt{J_{0}}}e^{-(\Gamma_{1}+i\Omega)(t-x+L/2)}[\Theta(x)-\Theta(x-t)]C_{R}^{\dagger}(x)\ket{0},$
(108b) $\displaystyle\ket{\psi_{L}(t)}$
$\displaystyle=\frac{-i}{\sqrt{2}}\int_{-\infty}^{\infty}\mathop{}\\!\mathrm{d}x\frac{\Gamma_{1}}{\sqrt{J_{0}}}e^{-(\Gamma_{1}+i\Omega)(t+x+L/2)}[\Theta(t+x)-\Theta(x)]C_{L}^{\dagger}(x)\ket{0}.$
(108c)

Let us start with the time derivatives

$\displaystyle i\partial_{t}\ket{\psi_{q}(t)}$
$\displaystyle=-i(\Gamma_{1}+i\Omega)\ket{\psi_{q}(t)},$ (109a) $\displaystyle
i\partial_{t}\ket{\psi_{R}(t)}$
$\displaystyle=-i(\Gamma_{1}+i\Omega)\ket{\psi_{R}(t)}+\frac{\Gamma_{1}}{\sqrt{2J_{0}}}e^{-i(\Gamma_{1}+i\Omega)L/2}C_{R}^{\dagger}(t)\ket{0},$
(109b) $\displaystyle i\partial_{t}\ket{\psi_{L}(t)}$
$\displaystyle=-i(\Gamma_{1}+i\Omega)\ket{\psi_{R}(t)}+\frac{\Gamma_{1}}{\sqrt{2J_{0}}}e^{-i(\Gamma_{1}+i\Omega)L/2}C_{L}^{\dagger}(-t)\ket{0}.$
(109c)

Now, let us find the action of the free Hamiltonian ($H_{0}$) on each part of
$\ket{\psi(t)}$

$\displaystyle H_{0}\ket{\psi_{q}(t)}={}$
$\displaystyle\Omega\ket{\psi_{q}(t)},$ (110a)
$\displaystyle\begin{split}H_{0}\ket{\psi_{R}(t)}={}&-i(\Gamma_{1}+i\Omega)\ket{\psi_{R}(t)}-\frac{\Gamma_{1}}{\sqrt{2J_{0}}}e^{-(\Gamma_{1}+i\Omega)t}e^{-i\theta/2}C_{R}^{\dagger}(0)\ket{0}\\\
&+\frac{\Gamma_{1}}{\sqrt{2J_{0}}}e^{-(\Gamma_{1}+i\Omega)L/2}C_{R}^{\dagger}(t)\ket{0},\end{split}$
(110b)
$\displaystyle\begin{split}H_{0}\ket{\psi_{L}(t)}={}&-i(\Gamma_{1}+i\Omega)\ket{\psi_{L}(t)}-\frac{\Gamma_{1}}{\sqrt{2J_{0}}}e^{-(\Gamma_{1}+i\Omega)t}e^{-i\theta/2}C_{L}^{\dagger}(0)\ket{0}\\\
&+\frac{\Gamma_{1}}{\sqrt{2J_{0}}}e^{-(\Gamma_{1}+i\Omega)L/2}C_{L}^{\dagger}(-t)\ket{0},\end{split}$
(110c)

where we define $e^{-(\Gamma_{1}+i\Omega)L/2}\simeq e^{-i\Omega
L/2}=e^{-i\theta/2}$.

Finally, we consider the action of the interaction Hamiltonian ($H_{I}$) on
each part of $\ket{\psi(t)}$

$\displaystyle H_{I}\ket{\psi_{q}(t)}$
$\displaystyle=\frac{\sqrt{J_{0}}}{\sqrt{2}}e^{-\Gamma_{1}t-i\Omega
t}[C_{R}^{\dagger}(L/2)+C_{L}^{\dagger}(L/2)+C_{R}^{\dagger}(-L/2)+C_{L}^{\dagger}(-L/2))]\ket{0},$
(111a) $\displaystyle H_{I}\ket{\psi_{R}(t)}$
$\displaystyle=\frac{-i}{\sqrt{2}}\Gamma_{1}e^{-(\Gamma_{1}+i\Omega)t}\ket{e_{2}},$
(111b) $\displaystyle H_{I}\ket{\psi_{L}(t)}$
$\displaystyle=\frac{-i}{\sqrt{2}}\Gamma_{1}e^{-(\Gamma_{1}+i\Omega)t}\ket{e_{1}}.$
(111c)

From our calculations, we see that
$i\partial_{t}\ket{\psi_{q}(t)}=H_{0}\ket{\psi_{q}(t)}+H_{I}\ket{\phi_{R}(t)}+H_{I}\ket{\phi_{L}(t)}$.
Eliminating other components in a similar fashion, we obtain

$\displaystyle\frac{\Gamma_{1}}{\sqrt{2J_{0}}}e^{-(\Gamma_{1}+i\Omega)t}e^{-i\theta/2}\left[C_{R}^{\dagger}(0)+C_{L}^{\dagger}(0)\right]\ket{0}$
(112a) $\displaystyle=\frac{\sqrt{J_{0}}}{\sqrt{2}}e^{-\Gamma_{1}t-i\Omega
t}[C_{R}^{\dagger}(L/2)+C_{L}^{\dagger}(L/2)+C_{R}^{\dagger}(-L/2)+C_{L}^{\dagger}(-L/2))]\ket{0}.$
(112b)

Now, we use the condition $C_{R}^{\dagger}(\pm L/2)\simeq e^{\mp
i\theta/2}C_{R}^{\dagger}(0)$, $C_{L}^{\dagger}(\pm L/2)\simeq e^{\pm
i\theta/2}C_{L}^{\dagger}(0)$ (which are consequences of the Markovian
approximation) and $\Gamma_{1}=J_{0}(1+e^{i\theta})$ (with the re-scaling
$\Gamma_{1}\to 2\Gamma_{1}$ at the beginning taken into account) to obtain

$\displaystyle
2J_{0}\cos(\theta)\left[C_{R}^{\dagger}(0)+C_{L}^{\dagger}(0)\right]\ket{0}=2J_{0}\cos(\theta)[C_{R}^{\dagger}(0)+C_{L}^{\dagger}(0))]\ket{0},$
(113)

which shows that (45) is a solution to the Schrödinger equation in the
Markovian limit.

##  Appendix J New insights gained

In this Appendix, we discuss various areas where we developed new insights as
a result of applying the real-space formalism to the problems considered here.

###  J.1 Collective decay rates

The most important intuition that we gained is the idea of collective decay
rates and how they dictate the collective behavior. When there is an effective
coupling between two or more isolated systems, the coupling shifts the energy
levels as well as the individual decay rates. In waveguide QED, this effective
coupling is mediated by bosonic fields (photons). The real part of the
collective decay rate corresponds to the coupling energy of the corresponding
interaction basis state, whereas the imaginary part (divided by 2) corresponds
to the shift in energies. Here, the interaction basis states, let’s call them
$\ket{S_{i}}$, are coherent excitations of single qubit excitation states such
that $\ket{S_{i}}=\sum_{j}\alpha_{j}^{(i)}\ket{e_{j}}$, where
$\alpha_{j}^{(i)}$ can be found by diagonalizing the coupling matrix as
described in Appendix F. Consequently, by finding poles of the scattering
parameters (i.e. highly mathematical objects), we gain information about the
shifts in energy levels and individual decay rates.

###  J.2 Causality

Another important aspect we saw is how the causality principle is related to
the position of scattering parameter poles (and correspondingly the collective
decay rates) in the complex plane. As we saw in the example of the decaying
exponential, the causality principle dictates that the poles should be in the
LHP for the linear chain of $N=3$ qubits. In fact, this is a more general
phenomenon and the scattering parameter poles should be in the LHP for _any_
waveguide QED system. We plan to prove this in future work. Here, we emphasize
another important feature of this mathematical connection. The imaginary part
of the poles gives the decay rates, and so the poles should be in the LHP for
the decay rate of the multi-qubit system to be positive. If the system had a
pole in the UHP, then one of the decay rates would be negative, that is the
system would get excited with more than unity probability after some time.
This is clearly not physical. What is striking is the fact that both this
unphysical behavior and the causality principle are linked by the same
mathematical identity, which can be found simply by considering _the steady-
state solution_ without any time dynamics calculations. We have not seen this
intuition discussed in the literature.

###  J.3 The Markovian limit

In [35], the authors linearize the phase of photons propagating between two
adjacent qubits as $kL\simeq\Omega L$. A similar substitution is later
referred to as a Markovian approximation in [33]. In the Markovian limit, the
inter-system propagation time of photons is neglected and only the phase
acquired by the propagating photon is accounted for. The latter is trivial to
see in the linearization process $e^{ikL}\simeq e^{i\Omega L}$ as the acquired
phase can be approximated by the phase of the resonant photon. Neglecting the
time delay portion, however, is not as straightforward to see in this
linearization. As we have seen in Section 3, the $k$\- dependence of the phase
$e^{ikL}$ shifts the time $t$ by $t-L$ (See, for example, Eq. (9b)). To see
this clearly, let us calculate the excitation probability of a qubit upon
pulse scattering:

$P_{m}(t)=\left|\int_{-\infty}^{\infty}\frac{\mathop{}\\!\mathrm{d}k}{2\pi}e_{m}S(k,0)\exp(-i\Delta_{k}t)\right|^{2}.$
(114)

Here, let us assume that the qubit is situated at position $x=L$ instead of
$x=0$. Then, $e_{m}$ would gain a phase of $e^{ikL}$, which would later be
picked up by $\exp(-i\Delta_{k}t)$ term to shift the time $t\to t-L$. This
leads to a delay of $L$ in time. If, however, $e^{ikL}$ was to be linearized
by replacing $k\simeq\Omega$, the time delay would not be accessible since
$e^{i\Omega L}$ does not depend on $k$ and can be taken out of the integral.
By linearizing $kL\simeq\Omega L$, we effectively disregard any time delay
that is of order $L$. Consequently, the real-space approach shows how this
linearization is equivalent to applying the Markovian approximation.
Consequently, the non-Markovian behavior of multi-qubit systems is extremely
easy to obtain using the real-space approach, where we simply omit the
linearization.

###  J.4 Black-box behaviour

The real-space approach also explains the black box behavior we see in
waveguide QED systems. When a pulse is incident on a system, the modulation of
the transmitted/reflected pulse depends only on the external scattering
parameters. Each momentum mode of the initial pulse gets modulated via the
external parameters and the asymptotic shape of the pulse at $t\to\infty$ can
be found without any time dynamics consideration. We proved this using the
real-space approach in Appendix H for a specific case, however the general
proof is analogous.

###  J.5 Non-Markovian dynamics

Finally, the real-space approach makes it possible to calculate exact and
fully analytical non-Markovian dynamics. The fact that non-Markovian behavior
for a single qubit system was explored only recently [32, 46], means that
probing fully analytical multi-qubit non-Markovian dynamics would be a big
step forward. We will discuss an approach to this problem in upcoming work.
For now, we note that this process is analogous to scattering from a finite
well in introductory quantum mechanics. The elementary nature of the real-
space approach guides our intuition about simple problems toward solving more
complicated ones.