AS and MZ are funded by the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No 714235). MZ is also supported by the Austrian Science Fund FWF within the DK-ALM (W1259-N27). MLJ is funded by MaQsens project (European Union’s Horizon 2020 research and innovation program, grant agreement No 736943 and European Research Council under the European Union’s Horizon 2020 research and innovation program, grant agreement No 649008). J.J.G.R. acknowledges support from project PGC2018-094792-B-I00 (MCIU/AEI/FEDER, UE), CSIC Research Platform PTI-001, and CAM/FEDER project No. S2018/TCS-4342 (QUITEMAD-CM). Facilities use was supported by the KIT Nanostructure Service Laboratory (NSL).

APPENDICES
Wigner-Weisskopf model in a rectangular cross section waveguide
In this section we discuss to what extend the Wigner-Weisskopf approach, developed originally for an atomic decay in 3D vacuum is relevant to our experimental situation of a rectangular cross section waveguide.
We follow [1]} and start from the very general equation for the state vector of the atom-photon system:
\({\Psi (t)}=c_{e}(t){e,0}+\sum _{k}c_{g,k}(t){g,1_{k}}\) 
with initial condition \(c_{e}(0)=1\)  and \(c_{g,k}(0)=0\) . We are going to use the standard Hamiltonian in the interaction picture and within rotating wave approximation:
\(H=\hbar \sum _{k}(g_{k}^{*}\sigma _{+}a_{k}e^{i(\omega _{q}-\omega _{k})t}+h.c.)\) 
where \(a_{k}\)  is the annihilation operator for the field mode with wave vector \(k\) , \(g_{k}=g_{k}(z_{0})=g_{k}(0)e^{-ikz_{0}}\)  is the coupling between the field mode \(k\)  and the atom positioned in the waveguide at the point \(z_{0}\) , and h.c. refers to hermitian conjugate. In the further discussion we will assume that the atom is located at the point \(z_{0}=0\) . The time dependent Schrödinger equation for these \({\Psi }\)  and \(H\)  leads to a set of equations for \(c_{e}(t)\)  and \(c_{g,k}(t)\) :
\(\dot{c}_{e}(t)=-i\sum _{k}g_{k}^{*}e^{i(\omega _{q}-\omega _{k})t}c_{g,k}(t)\) 
\(\dot{c}_{g,k}(t)=-ig_{k}e^{-i(\omega _{q}-\omega _{k})t}c_{e}(t)\) 
Reducing this set to only one equation for \(c_{e}(t)\) , one can write:
\(\dot{c}_{e}(t)=-\sum _{k}|g_{k}|^{2}\int _{0}^{t}dt^{\prime }e^{i(\omega _{q}-\omega _{k})(t-t^{\prime })}c_{e}(t^{\prime })\) 
Up to this moment the environmental properties have not played their role yet. The next step - changing the sum with the integral - has to take the 3D waveguide into account:
\(\sum _{k}\rightarrow 2\frac{L}{2\pi }\int _{0}^{\infty }dk.\) 
The factor 2 here comes from two possible directions of \(k\) , while there is only one possible polarisation direction, as we consider only decay into the fundamental mode of the waveguide. \(L\)  is a "quantization length" playing a similar role in the waveguide as the quantization volume in an open 3D space. With respect to REF , it is important to notice that we understand \(k=|k|\)  as an effective wavenumber in the waveguide, i.e. as \(\sqrt{{k}_{r}^{2}-{k}_{c}^{2}}\) , where \({k}_{r}\)  is the wavenumber in free space, and \(k_{c}=\sqrt{(\frac{m\pi }{a})^2+(\frac{n\pi }{b})^2}\)  is the cutoff wavenumber of the rectangular waveguide with cross section dimensions \(a\)  and \(b\) . As we are interested only in the first TE spatial mode (\(m\) =1 and \(n\) =0, providing \(a>b\) ), in our case \(k_{c}=\frac{\pi }{a}\) .
<FIGURE>To move further we have to assume a specific model of interaction between the field and the emitter; for our experimental situation an electric dipole approximation is reasonable. Therefore we define atom-field coupling as \(g_{k}=\frac{(E_{k}d)}{\hbar }\)  and use expression for an electric field amplitude created by a photon with frequency \(\omega _{k}\)  in a specific location \(x_{0}\)  inside the waveguide (see figure REF ) \(E_{k}=(\frac{\hbar \omega _{k}}{\epsilon V})^{\frac{1}{2}}\cdot sin(\frac{x_{0}\pi }{a})\) , where \(V=abL\)  is the quantization volume. This brings us to:
\(g_{k}=\sqrt{\frac{\omega _{k}}{\hbar \epsilon V}}|d_{eg}|cos(\vartheta )sin(\frac{x_{0}\pi }{a}).\) 
with \(d_{eg}\)  being the transitional dipole moment of the atom and \(\theta \)  is the angle between \(E_{k}\)  and \(d_{eg}\) .
Using REF  and REF  we rewrite REF :
\(\dot{c}_{e}(t)=-2\frac{L}{2\pi }\int _{0}^{\infty }|d_{eg}|^2cos^{2}(\vartheta )sin^{2}(\frac{x_{0}\pi }{a})\frac{\omega _{k}}{\hbar \epsilon V}dk\cdot \int _{0}^{t}dt^{\prime }e^{i(\omega _{q}-\omega _{k})(t-t^{\prime })}c_{e}t^{\prime }.\) 
Using the dispersion relation in the waveguide \(k=\sqrt{\epsilon \mu \omega ^{2}_{k}-k^{2}_{c}}\)  we switch to frequencies instead of wavenumbers:
\(\dot{c}_{e}(t)=-\frac{\mu |d_{eg}|^{2}cos^{2}(\vartheta )sin^{2}(\frac{x_{0}\pi }{a})}{\pi \hbar ab}\cdot \int _{\omega _{c}}^{\infty }\frac{\omega ^{2}_{k}d\omega _{k}}{\sqrt{\epsilon \mu \omega ^{2}_{k}-k_{c}^{2}}}e^{i(\omega _{q}-\omega _{k})(t-t^{\prime })}\cdot \int _{0}^{t}c_{e}(t^{\prime })dt^{\prime },\) 
where \(\omega _{c}=k_{c}/\sqrt{\epsilon \mu }\)  is the cutoff frequency of the waveguide.
At this moment we make two other assumptions which are often referred to as "Wigner-Weisskopf approximation". Both assumptions are based on the idea that at any time during the decay only frequencies close to \(\omega _{q}\)  can be populated. With this in mind, we exchange in REF  \(\omega ^{2}_{k}\)  with \(\omega ^{2}_q\)  and extend the lower integration limit to \(-\infty \) . This means that \(\omega _{q}\)  is sufficiently higher than \(\omega _{c}\) , so that in the frequency range we are interested in \(\frac{\omega ^{2}_{k}}{\sqrt{\epsilon \mu \omega ^{2}_{k}-k_{c}^{2}}}\)  can be assumed to be constant. In this case REF  becomes:
\(\dot{c}_{e}(t)=-\frac{\mu |d_{eg}|^{2}cos^{2}(\vartheta )sin^{2}(\frac{x_{0}\pi }{a})}{\pi \hbar ab}\frac{\omega ^{2}_{q}}{\sqrt{\epsilon \mu \omega ^{2}_{q}-k_{c}^{2}}}\cdot \int _{-\infty }^{\infty }d\omega _{k}e^{i(\omega _{q}-\omega _{k})(t-t^{\prime })}\int _{0}^{t}c_{e}(t^{\prime })dt^{\prime }.\) 
Using one of the Dirac delta function definitions \(2\pi \delta (t^\prime -t)=\int _{-\infty }^{\infty }d\omega _{k}e^{i(\omega _{q}-\omega _{k})(t-t^{\prime })}\)  and after integration over \(t^\prime \)  we finally arrive at the usual exponential decay:
\(\dot{c}_{e}(t)=-\frac{\gamma }{2}c_{e}(t)\) 
\(\gamma =\frac{2\mu |d_{eg}|^{2}cos^{2}(\vartheta )sin^{2}(\frac{x_{0}\pi }{a})}{\hbar ab}\frac{\omega ^{2}_{q}}{\sqrt{\epsilon \mu \omega ^{2}_{q}-k_{c}^{2}}}=\frac{2|d_{eg}|^{2}cos^{2}(\vartheta )sin^{2}(\frac{x_{0}\pi }{a})}{\hbar ab}\frac{\omega ^{2}_{q}\cdot Z_{r}}{\sqrt{\omega ^{2}_{q}-\omega _{c}^{2}}}.\) 
\(Z_{r}=\sqrt{\frac{\mu }{\epsilon }}\)  in the last expression is the impedance of the waveguide medium. Plugging the usual exponential solution of REF  into REF , one can arrive at the Wigner-Weisskopf picture from the equations REF -REF  from the main text. Therefore, providing all the aforementioned approximations are correct, one should expect the typical Wigner-Weisskopf picture of exponential decay to be accurate in the 3D waveguide as well, apart of expression for \(\gamma \) . Note, that for the atomic frequencies \(\omega _{q}\)  close to the cutoff \(\omega _{c}\)  the assumptions of the Wigner-Weisskopf approximation are not correct and the decay is not exponential. For the experiment, described in the main text we are detecting frequency components of the decay within a 60 MHz frequency window around the qubit frequency 7.3 GHz. Within this window the factor \(\frac{\omega ^{2}_{k}}{\sqrt{\epsilon \mu \omega ^{2}_{k}-k_{c}^{2}}}\)  changes around 1.7% which is small enough not to affect the decay.

Methods
In the experiment we use a standard 12cm long WR90 copper rectangular waveguide with cross section dimensions 10.16 mm and 22.86 mm and cutoff frequency 6.557 GHz. On both ends it is connected to 50\(\Omega \)  coaxial cables with cable-to-waveguide adapters. Superconducting coils placed on the external surface of the waveguide were used to tune the qubits frequencies. The DC lines of the coils were heavily filtered to avoid additional flux noise at the qubit locations, including home-made steel powder filters [2]} installed on the top of the shielding can. The waveguide was thermally anchored to the baseplate of a dilution refrigerator and operated at \(\sim \) 20 mK. Qubits were inserted into the waveguide through specially designed slits, and mechanically fixed with copper clamps, allowing for appropriate thermalization. Appropriate RF and DC shielding of the sample was provided by a double layer cryoperm shield, and a niobium shield.
For qubit fabrication we used a standard bridge-free Al/AlO\(_{x}\) /Al technique, described in details elsewhere [3]}. In short, electron beam lithography was used to pattern junctions and antennas on a 2-inch 330 \(\mu \) m thick commercial sapphire substrate. An electron-gun evaporator was used to deposit 2 layers of Al on the substrate with a thickness of 20 and 30 nm respectively. Josephson junctions were formed between the 2 layers by 2.5 minutes oxidation under 15 mbar pressure, resulting in a critical current density \(\sim \) 50 A/cm\(^2\) . To achieve appropriate balance between flux tunability and susceptibility to the flux noise, both qubits had nonidentical junctions. The sapphire wafer was subsequently diced with a diamond saw, to 2.5 mm-wide and 2.5 cm-long chips.

Experimental setup
<FIGURE>The experimental setup is shown in figure REF . We are generating our pulses at low (\(\sim \) 200 MHz) frequency from 2 independent channels of an Arbitrary Waveform Generator and up-convert them to the qubit frequency with an IQ-mixer. The pulses are attenuated by 50 dB before entering the waveguide installed at the baseplate of the fridge. After interaction with the qubits, the signal is first amplified by 20 dB within a narrow frequency band (\(\sim \) 5 MHz) at the base plate of the refrigerator using a JPC and consequently at the 4K stage by a High Electron Mobility Transistor amplifier (+40 dB) and at room temperature by another low-noise amplifier (also +40 dB). The input line of the fridge was equipped with a home-made Eccosorb infrared filter as well as a 12 GHz cutoff low-pass filter (both are not shown). We use a conventional mixer for down converting the signal to 50 MHz and an analog-to-digital converter for digitizing and averaging the signal amplitude. The effective noise temperature without and with JPC was estimated to be about 10.5 K and 80 mK respectively.
Superconducting coils, used to tune the qubit frequencies, were powered from current sources through a set of low-pass filters: 1.9 MHz cutoff, homemade 500 kHz and 3 kHz cutoff resistive Pi-filters and steel powder filters with a \(\sim \) 30 kHz cutoff. The filters were installed inside an aluminum housing on the base plate of the refrigerator outside the shielding-can of the sample.

Theoretical model: "Infinitely long" pulse
The next two appendices describe the model, based on input-output approach, used to fit the experimental data from the main text. This subsection is devoted to the situation in the main text referred to as "infinitely long" pulse, where everything can be derived analytically. In the next subsection we describe the model we solve numerically for arbitrary length pulses.
Let us consider a driven 2-level system in the framework of optical Bloch equations and input-output theory:
\(\begin{split}&\frac{\mathrm {d}}{\mathrm {dt}} \left< \sigma ^z(t) \right> =- 2\mathrm {i}\sqrt{\frac{\gamma \beta }{2}}\Bigl [ \alpha e^{-\mathrm {i}\omega t}\left<\sigma ^+(t)\right>- \alpha ^{*} e^{\mathrm {i}\omega t}\left<\sigma ^-(t)\right>\Bigr ]- \gamma \left( \left< \sigma ^z(t) \right> + 1\right)\\&\frac{\mathrm {d}}{\mathrm {dt}} \left< \sigma ^-(t) \right> = -\left( \mathrm {i}\omega _{q} + \frac{\gamma }{2}\right)\left<\sigma ^-(t)\right> + \mathrm {i}\alpha e^{-\mathrm {i}\omega t}\sqrt{\frac{\gamma \beta }{2}} \left< \sigma ^z(t) \right>\\&\left<a_{out}\right> = \left<a_{in}\right> - \mathrm {i}\sqrt{\frac{\gamma \beta }{2}} \left< \sigma ^-(t) \right>,\end{split}\) 
where \(\omega _{q}\)  is the qubit transition frequency, \(\gamma \)  is the qubit linewidth, related to both, spontaneous emission into the waveguide and internal losses with \(\gamma =\gamma _{waveguide}+\gamma _{internal}\) ), \(\alpha \)  is the drive amplitude and \(\beta =\gamma _{waveguide}/(\gamma _{waveguide}+\gamma _{internal})\)  is a parameter characterising the portion of the qubit energy emitted into the waveguide.
Switching into the rotating (with the drive) frame and changing notation for \(\left< \sigma ^z(t)\right>\) :
\(\left< \sigma ^-(t) \right> = e^{-\mathrm {i}\omega t}s\textnormal {, }\left< \sigma ^+(t) \right> = e^{\mathrm {i}\omega t}s^{*}\textnormal {, }\left< \sigma ^z(t) \right> = 2\rho - 1,\) 
where \(s\)  represents the coherence and \(\rho \)  is the excited population of the qubit. Both are changing slowly in comparison with \(\omega \) . This yields:
\(\begin{split}&\frac{\mathrm {d}}{\mathrm {dt}} \rho =- \mathrm {i}\sqrt{\frac{\gamma \beta }{2}}\left[ \alpha s^{*} - \alpha ^{*}s \right] - \gamma \rho \\&\frac{\mathrm {d}}{\mathrm {dt}} s = -\left[ \mathrm {i}\left( \omega _{q} -\omega \right) + \frac{\gamma }{2}\right]s + \mathrm {i}\alpha \sqrt{\frac{\gamma \beta }{2}} \left(2\rho -1\right)\end{split}\) 
which for steady state reduces to:
\(s_{0} = \frac{\mathrm {i} \alpha \sqrt{\frac{\gamma \beta }{2}}}{\mathrm {i}\left( \omega _{q} -\omega \right)+ \frac{\gamma }{2}}\left( 2\rho _{0} - 1\right),\) 
and:
\(\rho _{0}= -\mathrm {i}\sqrt{\frac{\beta }{2\gamma }} \left( \alpha s_{0}^{*} - \alpha ^{*} s_{0} \right)\) 
Choosing the initial phase so \(\alpha \)  is real, we finally obtain:
\(\begin{split}& s_{0} = -\frac{\alpha \sqrt{\gamma \beta /2}\left( \omega _{q} - \omega +\mathrm {i}\gamma /2\right)}{\left( \omega _{q} - \omega \right)^2 + \gamma ^2/4 +\alpha ^2\beta \gamma }\\&\rho _{0} = \frac{\alpha ^2\beta \gamma /2}{\left( \omega _{q} - \omega \right)^2 + \gamma ^2/4 +\alpha ^2\beta \gamma }\end{split}\) 
We are going to use the obtained steady state solution as an initial condition to fit the qubit decay without drive. The drive is abruptly switched off at a given time \(t=0\) . The system then evolves according to:
\(\begin{split}&\frac{\mathrm {d}}{\mathrm {dt}} \left< \sigma ^z(t) \right> = - \gamma \left( \left< \sigma ^z(t) \right> + 1\right)\\&\frac{\mathrm {d}}{\mathrm {dt}} \left< \sigma ^-(t) \right> = -\left( \mathrm {i}\omega _{q} + \frac{\gamma }{2}\right)\left<\sigma ^-(t)\right>\\&\left<a_{out}\right> = - \mathrm {i}\sqrt{\frac{\gamma \beta }{2}} \left< \sigma ^-(t) \right>,\end{split}\) 
This leads to the usual exponential decay with initial condition set by \(s_{0}\)  and \(\rho _{0}\) :
\(\begin{split}&\left< \sigma ^z(t) \right> = 2e^{-\gamma t}\rho _{0} - 1\textnormal {, }\left< \sigma ^-(t) \right> = e^{-\mathrm {i}\omega _{q} t - \gamma t /2}s_{0}\\&\left< a_{out}(t) \right> = -\sqrt{\frac{\gamma \beta }{2}} e^{-\mathrm {i}\omega _{q} t - \gamma t /2}s_{0}\end{split}\) 
Now we have to add the amplifier to the picture:
\(a_{amp}(t) = \sqrt{G} \int _{-\infty }^{t} \left<a_{out}(\tau )\right> e^{-\left[\mathrm {i}\omega _{amp}+\kappa \right]\left(t-\tau \right)}\mathrm {d}\tau +h,\) 
where \(\omega _{amp}\)  is the center frequency of the amplifier band \(\kappa \) , \(G\)  is the gain, while \(h\)  is a noise operator.
For the transmitted field after the amplifier, we then obtain:
\(\left<a_{amp}(t)\right>=\sqrt{G} \int _{-\infty }^t e^{-\left[\mathrm {i}\omega _{amp}+\kappa \right]\left(t-\tau \right)} \Bigl [ e^{-\mathrm {i}\omega t}\left( \alpha -\mathrm {i} s_{0} \sqrt{\frac{\gamma \beta }{2}}\right)\Theta (-\tau )-s_{0}\sqrt{\frac{\gamma \beta }{2}} e^{-\mathrm {i}\omega _{q} t - \gamma t /2} \Theta (+\tau ) \Bigr ] \mathrm {d}\tau \) 
where \(\Theta (t)\)  is the Heaviside step function. This can be evaluated analytically and after some algebra one arrives at the final result:
\(\begin{split}\left< a_{amp}(t)\right> = & \frac{e^{-t (\kappa + \mathrm {i} \omega _{amp}) + (\kappa - i (\omega - \omega _{amp})) \mathrm {Min}[0,t]} \left( \alpha - \frac{\alpha \beta \gamma /2 (\gamma /2 - \mathrm {i} (\omega _{q} - \omega )) }{ \alpha ^2 \beta \gamma + \gamma ^2/4 +(\omega _{q} - \omega )^2}\right)}{\kappa -\mathrm {i} (\omega - \omega _{amp})}+\\& + \frac{e^{-t (\kappa + \mathrm {i} \omega _{amp})} \left( -1 +e^{(-(\gamma /2) - \mathrm {i} \omega _{q} + \kappa +i \omega _{amp}) \mathrm {Max}[0,t]} \right) \alpha \beta \gamma /2 (\gamma /2 - \mathrm {i} (\omega _{q} - \omega ))}{(\alpha ^2 \beta \gamma + \gamma ^2/4 + (\omega _{q} - \omega )^2) (\gamma /2 + \mathrm {i} (\omega _{q} + \mathrm {i} \kappa - \omega _{amp}))}.\end{split}\) 
In this expression the first term describes the system behaviour at the times \(t<0\)  when the drive is on and the second term is 0. After the drive stops (at \(t=0\) ), the second term describes the qubit decay (with the rate \(\gamma /2\)  instead of \(\gamma \)  as it corresponds to the amplitude reduction rate), while the first one describes the amplifier decay (with the rate \(\kappa \) ). Formula REF  was used to fit our experimental data for the from the figures 4(a) and (c) in the main text. Figure REF  shows how well the fit works for the state of one qubit and the collective state of two qubits for different values of the drive-JPC detuning. In the experiment we intend to keep the drive frequency in resonance with the qubit. However, qubit frequency stability was noticeably affected by magnetic field noise in the cryostat. Since the theory does not require us to fix \(\omega =\omega _{q}\)  the qubit-drive detuning can be extracted as a result of the fit. Fitted experimental data for different detunings \(\delta \omega \)  for the one- and two-qubits cases, similar to figures 4(a) and (c) from the main text, are shown in figures REF (a) and REF (c). The qubit drive detuning extracted from the fit for every frequency is shown in figure REF (b), and mostly does not exceed 200 kHz. Comparison of the figure REF (a) and 4(a) from the main text leads us to the conclusion that this detuning is the main source of noise for our experimental data.
The theory only describes the system dynamics in the case of the "infinitely long" pulse, after the dynamical equilibrium between the drive and qubit decay has been established. The case of short (\(\pi /2\) , \(\pi \) , \(3\pi /2\) ...) pulses has been treated numerically and is described in the next subsection.
<FIGURE><FIGURE>
Theoretical model: Arbitrary length pulses
In order to simulate the time-dynamics during the onset of the drive pulse, the system was also modeled numerically using QuTiP [4]}, [5]}. The model, accounting for multiple qubits, is a modification of the model described in [6]} to add the direct 'dipole-dipole' coupling term between the qubits \(\mathcal {G}_{ij}\) :
\(H =H_{qubits}+H_{qubit-drive}+H_{qubit-qubit}\) 
\(H_{qubits} = \sum _{i=1}^2\omega _{i} \left|e_i\right> \left<e_i\right|\) 
\(H_{qubit-drive} = \sum _{i=1}^2\left(\epsilon _i \sigma _i^{+}+h.c.\right)\) 
\(H_{\textnormal {qubit-qubit}} = \sum _{j> i=1}^2\left( J_{ij} + \mathcal {G}_{ij}\right) \left( \sigma _i^- \sigma _j^+ + \sigma _i^+\sigma _j^-\right)\) 
\(\dot{\rho }= -\frac{i}{\hbar }\left[H,\rho \right] + \sum _{i,j=1}^2 \gamma _{ij}\left[ \sigma _i^- \rho \sigma _j^+ -\frac{1}{2}\lbrace \sigma _i^+\sigma _j^-, \rho \rbrace \right].\) 
Here \(\omega _{i}\)  is the frequency of the i-th qubit, \(g_{i}\)  is its dimensionless coupling to the waveguide, \(\epsilon _i\)  is the drive amplitude at the position of the i-th qubit, \(J_{ij}\)  is the exchange interaction amplitude and \(\gamma _{ij}\)  are the elements of the radiative decay matrix. In the case when two qubits are tuned in resonance and only driven from one side, the expressions for \(J_{ij}\) , \(\gamma _{ij}\)  and \(\epsilon _{ij}\)  are:
\(\epsilon _i = -i\sqrt{ \frac{\gamma ^0_{ii}\omega }{2\omega _i}}\left<a_{in}\right>\textnormal {e}^{-i \omega t_i}\) 
\(J_{i j} = 2\pi g_i g_j \omega _i \textnormal {sin}\left( \omega _i t_{i j}\right)\) 
\(\gamma _{i j} = 4\pi g_i g_j \omega _i \textnormal {cos}\left( \omega _i t_{i j}\right) + \delta _{ij}\gamma _{internal}.\) 
Here \(t_i=x_i/\nu \)  defines the phase of the input signal \(\left<a_{in}\right>\)  propagating through the waveguide with a speed \(\nu \)  seen by i-th qubit located at point \(x_i\)  along the waveguide. Similarly, \(t_{ij}=|x_i-x_j|/\nu \)  defines the time delay of a signal propagating from the i-th to the j-th qubit. As before, \(\omega \)  is the drive frequency and \(\gamma _{internal}\)  is the non-radiative decay rate. In turn, \(\gamma _{ii}^0\)  is the radiative decay rate of the i-th qubit, i.e. \(\gamma _{ii}\)  in equation REF  without \(\gamma _{internal}\) . In our configuration the qubits are located in a 3D waveguide effectively forming a closely packed system. In this case, the coupling term obtained form the waveguide input-output formalism gives \(J_{ij}\) =0 since the near-field components of the electromagnetic field are not considered. Such near-field term is instead included by introducing a direct qubit-qubit coupling term, \(\mathcal {G}_{ij}\) , which is related to the real part of the free space Green function. The output field is defined in the usual way, which is:
\(\left<a_{out}\right> = \left<a_{in}\right> + \sum _{i=1}^2 \textnormal {e}^{i \omega _i t_i} \sqrt{\frac{\gamma ^0_{ii}}{2}}\left<\sigma _i^-\right>.\) 
Solving these equations numerically allows us to reconstruct the field propagating along the waveguide regardless of the length of the drive pulse. This method was used for state the tomography presented in figure 3(a,b) of the main text. Figure REF  shows \(\left<a_{out}\right>\)  during the entire 1\(\mu \) s pulse: (a) with qubit, and (b) without qubit.
<FIGURE>
Two-"atom" hybridization
Two resonantly interacting transmon qubits in the rotating wave approximation can be described by the interaction Hamiltonian \(H_{int}=g(\sigma ^{+}_1\sigma ^{-}_{2}+h.c.)\) , which in the single excitation manifold has two eigenstates \(\frac{1}{\sqrt{2}}({01}\pm {10})\) . When two identical transmons are placed symmetrically in the waveguide, within the electric dipole approximation it is easy to show, that the dashed transitions in figure REF (a) are prohibited. Indeed, the transition probability is proportional to the square of the corresponding matrix element of the electric dipole operator \(\sim |d_{fi}|^{2}\) . This element can be expressed in terms of the wave functions of both transmons, given in charge representation. For example for the \({11}\rightarrow {01}\pm {10}\)  transition it reads:
\(\begin{split}&d_{{11}\rightarrow {01}\pm {10}}=\\&=\frac{l}{\sqrt{2}}\int \left(\psi _{{1}}(Q_{1})\psi _{{0}}(Q_{2})\pm \psi _{{0}}(Q_{1})\psi _{{1}}(Q_{2})\right) ^{*}(Q_{1}+Q_{2})\left(\psi _{{1}}(Q_{1})\psi _{{1}}(Q_{2})\right)dQ_{1}dQ_{2}=\\&\hspace{199.16928pt}=\frac{1}{\sqrt{2}}(d_{{1}\rightarrow {0}}\pm d_{{1}\rightarrow {0}})=\sqrt{2}d_{{1}\rightarrow {0}}\;or\;0\end{split}\) 
where \(l\)  is proportional to the transmon length and \(Q_{1}(Q_{2})\)  is the charge operator of the first (second) transmon. According to equation REF  the probability of allowed transition is, in turn, increased twofold with respect to the single transmon case. States \(\frac{1}{\sqrt{2}}({01}\pm {10})\)  are therefore referred as dark and bright respectively. Note, that their corresponding eigenvalues depend on the sign of the interaction constant \(g\) . The experiment with a photon emitted from a collective 2-transmon state has been done using the transition marked in figure REF (a) with a red arrow. Figure REF (b) shows transmission measurements through the waveguide with 2 transmons. The color code on the picture corresponds to the transmission coefficient \(|S_{21}|\) . The current in a coil coupled to both qubits is swept along the X-axis. To emphasise the difference between dark and bright states around the avoided crossing for this measurements, the qubits were installed in a way to have a larger coupling to the waveguide.
<FIGURE>