Optical precursors in waveguide quantum electrodynamics

When a broadband signal propagates through a dispersive medium, some frequency components move faster than the center of the pulse. This leads to the appearance of precursors, transient signals that emerge from the medium earlier than the main part of the pulse and seem to propagate superluminally. Here, we investigate the microscopic origin of precursors in a minimal setup: an array of qubits coupled to a waveguide. The linear transmission function only converges to that of a continuous medium for large qubit numbers. Nevertheless, the dispersion produced by only two qubits is enough to produce oscillatory transients. Precursors are best observed under conditions of electromagnetically-induced transparency, as the center of the pulse is significantly delayed. Under these conditions, just a single qutrit is enough to generate a precursor. Our results pave the way towards dispersion engineering of light with just a few qubits, and can be realized with superconducting qubits coupled to transmission lines or atoms coupled to optical waveguides.


I. INTRODUCTION
The propagation of light through a continuous dispersive medium is a canonical problem in electrodynamics [1,2]. One of the most fascinating aspects of the transmitted radiation is the formation of transients that precede and follow the main pulse when the input signal has sharp edges, compared to the response time of the system. These transients are known as precursors and were theoretically predicted by Sommerfeld and Brillouin in the early 20th century [3,4]. They occur due to the different group velocities of the high-and low-frequency components in the spectrum of the original pulse. Since their prediction, precursors have been extensively studied theoretically [5][6][7][8]. They have also been observed in a plethora of systems and frequency ranges, such as in microwaves propagating in transmission lines [9], optical photons traversing atomic clouds [10][11][12], gamma rays in Mössbauer spectroscopy [13][14][15][16][17] and mechanical waves in fluids [18,19].
In waveguide quantum electrodynamics (wQED), where a collection of qubits are coupled to each other via a one-dimensional (1D) photonic reservoir, dispersion arises due to the narrow spectral response of the qubits. These two-level systems (which represent neutral atoms, molecules, or superconducting qubits, to name a few examples) only interact with photons of frequency resonant with their ground-to excited-state transition. Photon-mediated interactions between qubits give rise to the emergence of collective states that can either decay rapidly (superradiant) or be protected from dissipation (subradiant) [20]. In the last decade, wQED has attracted significant interest due to the possibility of exploiting these states for quantum information processing and storage (for instance, to produce quantum states * ana.asenjo@columbia.edu of light [21] and to compute via decoherence-free subspaces [22]) as well as for exploring many-body physics in open quantum systems (many-body localization [23], spin dimerization [24], and fermionization [20], among other examples). A lot of work has been devoted to single- [25], few- [26][27][28][29][30][31], and many-photon [32] transport in wQED, although most of it (except for a few exceptions [33][34][35]) is focused on the steady-state regime or on propagation of quasi-monochromatic light. In parallel, experimental realizations of wQED systems have multiplied, with setups ranging from neutral atoms coupled to single-mode optical fibers [36][37][38][39] and photonic-crystal waveguides [40][41][42], to superconducting qubits coupled to microwave transmission lines [43,44].
Here, we investigate the transport of broadband photon pulses in wQED, for a system consisting of N qubits coupled to a waveguide. Employing a transmission coefficient in terms of collective frequency shifts and decay rates, we demonstrate that the temporal response under a short pulse coincides with that of a continuous medium for N 1. The macroscopic description of the medium breaks down for a small qubit number. Nevertheless, just two qubits generate enough dispersion to produce an intensity profile that oscillates rapidly in time. The delay between the main signal and its precursor is evident under conditions of electromagnetically induced transparency, where a single qutrit is enough to generate a precursor.

II. CONTINUOUS MEDIUM
The amplitude and phase modulation acquired by a monochromatic field after propagation through a (either classical or quantum) linear dispersive medium is encoded in the complex transmission function t(ω). For non-monochromatic input pulses, the transmitted field is simply where E 0 (ω) is the Fourier transform of the temporal profile of the input pulse. In standard electromagnetism, the complex relative permittivity (ω) determines the transmission coefficient, and is usually postulated phenomenologically to match the optical response of a continuous medium. A conventional model is that of a Lorentz oscillator with a single resonance of frequency ω 0 and damping coefficient Γ ω 0 . Then, the transmission coefficient reads [45] where b quantifies the strength of the light-matter interaction. Throughout this paper, we consider an input field of central frequency ω p with a square tempo- Assuming that the duration of the input pulse is larger than the time it takes the system to reach the steady state, we approximate the input signal as a step function to calculate the transients in the transmitted intensity right after switching off the input field. Similarly, to calculate the transmitted field immediately after switching on the input field, we approximate the pulse as E 0 (t) Setting t f = 0 (t i = 0), the Fourier transform for the rising (falling) edge is where +(−) corresponds to the raising (falling) edge and P stands for Cauchy principal value. Plugging E 0,{R,F } (ω) and t cont (ω) into the transmitted field expression in Eq. (1) results in where z ≡ ω − ω 0 + iΓ /2 and ∆ = ω p − ω 0 is the detuning between the central and resonance frequencies. After solving the integral (see Appendix A), the final form for the transients reads where J n is a Bessel function of the first kind. Here, t 0 = t i corresponds to the rising edge and t 0 = t f to the falling one. The sharp edges of the input signal translate into a broad spectrum in Fourier space, and the interference between different frequency components propagating at different group velocities gives rise to temporal oscillations in the transmitted intensity. The transmitted intensity consists only of the precursor and the final transient, since the main pulse has been absorbed and scattered to free space.

A. Model
We demonstrate that temporal oscillations in the transmitted field intensity (so-called "dynamical beats" in the Mössbauer literature [13][14][15][16][17]) are not unique to continuous classical media, but also occur in "granular" quantum systems, such as a chain of N > 1 qubits coupled to a 1D waveguide, as shown in Fig. 1. In this case, the optical response can be obtained by first tracing out the field and solving for the dynamics of the qubits [46], which interact with each other as they share a common electromagnetic environment. Then, one recovers the field evolution via an input-output formalism. In the linear (or single excitation) regime, the qubits' evolution is governed by the effective Hamiltonian H = H 1D + H + H drive [47][48][49], where Here, H 1D describes the qubit-qubit interaction, which occurs at a rate Γ 1D , depends on the lattice constant d, and is mediated by photons of wave-vector k 1D . The qubits may decay (independently from each other) into other non-guided modes at a rate Γ , and the presence of the waveguide imparts a Lamb shift J on their resonance frequency, as described by H . The qubit ensemble is being driven by a propagating pulse of Rabi frequency Ω(t), whose central frequency is detuned from the qubit resonance frequency by ∆ = ω p − ω 0 . In the above equations,σ i eg = |e i g i | is the coherence operator between the ith-qubit excited and ground states,σ i ee = |e i e i |, and H.c. stands for Hermitian conjugate. We note that the rotating wave approximation is justified as counterrotating terms produce rapidly oscillating contributions (at frequency ∼ 2ω 0 ) that average out in the timescales relevant for the transients (proportional to the inverse of Temporal evolution of the intensity at the detector, I d with (red) and without (blue) qubits, normalized to the maximum input intensity I0. The central part of the pulse is mostly absorbed, and the transmitted intensity is low except at the edges, where a precursor can be observed at the first edge. The qubits' response immediately after switch-on and off is transient, taking some time for the induced field to build up (and cancel the forward-propagating input field) and decay, respectively. The optical depth is N Γ 1D /Γ = 5 (N = 20, Γ 1D /Γ = 0.25), the detuning between the resonance and central frequency of the pulse is ∆ = 0.37Γ , the lattice constant is k 1D d = π/2, and t i(f ) represents the time at which the pulse is switched on (off).
the qubit linewidth). The transmitted light intensity is found via the input-output equation [48,49] whereÊ + is the positive-frequency component of the right-propagating field (normalized to have units of Rabi frequency), and the field is measured by a detector at a position z that lies beyond the last qubit. In this manuscript, the dispersion is solely due to the qubits, and the waveguide is considered to be dispersionless.

B. Transmission coefficient
The steady-state transmission is mostly determined by the optical depth OD ≡ N Γ 1D /Γ , as shown in Fig. 2(a), where systems with different number of qubits and decay rates but fixed optical depth display almost identical transmittance spectra. To calculate the transmitted light for a continuous wave drive (i.e., E 0 (t) = Ω 0 e iωpt ), we solve for the expectation value of the steady-state coherences (such that σ i eg = 0) and plug the result into the above input-output equation. The transmission coefficient is defined as where z left is a point immediately to the left of the qubit 1 and z right is a point immediately to the right of the qubit N . E + (z) is the expectation value of the positivefrequency component of the total field operatorÊ(z) in the steady state and E + p (z) = Ω 0 e ik 1D z is the input field. The transmission coefficient can be expressed in terms of collective shifts and decay rates, as we now derive (see [50] for full details).
The input-output equation states that the total field is the sum of the input field and the field radiated by the qubits, i.e., where σ n ge ≡ σ n ge is the expectation value for the coherences in the steady state and g(z, z ) = −i(Γ 1D /2)e ik 1D |z−z | . Defining g nm = g(z n , z m ), and E + p a N -dimensional vector whose entries are E + p,n ≡ E + p (z n ), the evolution equations for the expectation value of the coherences arė The steady state solutions of these equations (for whicḣ σ n ge = 0) are We express this in terms of collective modes, since the eigenvectors of g satisfy Using this identity, we find where {λ ξ } are the eigenvalues of g (and of H 1D defined in Eq. (6)). Plugging the steady state solution (12) into the expression for the field we obtain In the last expression we have adopted the shorthand notation (g(z)) j = g(z, z j ). Here, g(z, z ) is the propagator of the guided field that has been projected in the direction of the qubits' dipole transition. Physically, g(z, z ) describes the field at z that is generated by a dipole at z . By means of the trace-determinant lemma [50], this expression can be written in terms of eigenvalues only, yielding whereω 0 = ω 0 + J . The real (J ξ = Re{λ ξ }) and imaginary (Γ ξ = −2Im{λ ξ }) parts of these eigenvalues correspond to the frequency shifts and decay rates of the collective modes. The decay rates can be either superradiant (with Γ ξ > Γ 1D , and the largest one scaling as Γ ξ ∼ N Γ 1D ) or subradiant (with Γ ξ < Γ 1D , and the smallest one scaling as Γ ξ ∼ Γ 1D /N 3 ), and their actual values depend on the specific lattice constant [20,23]. For lattice constants such that k 1D d = nπ, with n being an integer, there is only one superradiant eigenvalue of decay N Γ 1D . In this so-called 'mirror configuration', there is only one collective mode coupled to the waveguide, and the array of qubits behaves effectively as a single qubit with a large decay rate [50].

C. Transients
The temporal dynamics of the transmitted intensity for a broadband input pulse is not determined solely by the optical depth, but is instead sensitive to the specific values of N and Γ 1D /Γ separately, as shown in Fig. 2(b). From Eqs. (1), (3) and (14), the transmitted field at the beginning and end of the pulse is We solve the remaining integral using the residue theorem and a semicircle that closes in the lower-half plane. There are two contributions: one from the simple pole at ω = ω p , and one due to the singularities in t(ω) in the lower-half plane. They read = −πi t N (ω p )e −iωpt , where {ω ξ } are the singularities of t N (ω).
Reintroducing the times t i and t f yields the final expression for the transmitted field where t 0 = t i at the rising edge and t 0 = t f at the falling edge.
The term originated by the simple pole ξ in the transmission coefficient for a discrete chain, t N (ω), describes the light emitted by the corresponding collective mode. The measured intensity at the rising and falling edges is thus a coherent sum over all the contributions of the collective modes and reads where∆ = ω p −ω 0 . The contributions from different modes give rise to a time-dependent slope in the decay. The most superradiant modes play an important role for shorter times. As these modes become depopulated, the most significant contributions originate from less superradiant states, yielding a smaller decay rate [51]. The transients cannot be faithfully reproduced by just including a few modes in Eq. (19).
Oscillations only occur for N ≥ 2, as they arise from interference between different collective modes. Two qubits is the minimum required number to produce oscillations that -only in this case -are periodic with a frequency equal to half of the difference between the two collective modes. For N 1, the intensity calculations agree with those obtained for a continuous medium [i.e., as described by Eq. (5)]. The agreement can be readily understood by noting that for |λ ξ | |ω−ω 0 +iΓ /2|, with ξ = {1, ..., N }, the transmission coefficient for a finite array reduces to which is precisely the transmission coefficient of a continuous Lorentz medium (Beer-Lambert law), as captured by Eq. (2), with resonant frequency ω 0 + J , damping coefficient Γ , and coupling strength b = | ξ λ ξ | = N Γ 1D /2. For a fixed optical depth, the region where the series expansion is valid increases with qubit number [see Fig. 2(c)], thus the approximation works better for N 1. As exemplified in Fig. 3, in this limit the temporal evolution of the intensity is independent of the qubit spatial configuration, is robust against imperfect filling of the array, and is dictated only by the optical depth. This occurs for any lattice constant different from that of the mirror configuration.
Moreover, Eq. (5) captures the temporal response even for large optical depths (as long as N 1, N Γ 1D /Γ can take any value), as shown in Fig. 2(b) for N = 200 qubits and N Γ 1D /Γ = 5, even if t N cannot be approximated by t cont at resonance. This occurs because the temporal response for a broadband pulse involves an integral over frequencies, which is less sensitive to the specific details of the response function at resonance, compared to the steady-state transmission.
The breakdown of the continuous approximation for a few qubits can be understood by analyzing the pole structure of the two transmission coefficients in the complex plane, as shown in Fig. 2(c). For a finite array, each qubit (or more specifically, each collective mode) contributes with one simple pole. As N → ∞, the poles cluster around ω * =ω 0 − iΓ /2, which is an essential singularity of t cont (ω), and the response coincides to that of a continuous medium. As the number of qubits decreases, the poles do not densely cover the region of the essential singularity. Nevertheless, due to the frequency splitting between the collective modes, the poles have finite real parts, producing oscillations in the time domain. Lastly, N = 1 has a single pole at a purely imaginary frequency, −i(Γ 1D + Γ )/2, giving rise to exponential, nonoscillatory decay. The approximation in Eq. (20) is still valid for a single atom as long as Γ 1D Γ . However, in this limit, t N 1 (ω) roughly predicts a purely exponential decay, since the oscillations would occur in a timescale (∼ Γ −1 1D ) that is much larger than the decay (∼ Γ −1 ). Eq. (19) describes the transients only when the Markov approximation is valid, i.e., when retardation is negligible, and when the bandwidth of the reservoir is much larger than the linewidth of the qubits [thus making Γ 1D and k 1D approximately constants in the frequency interval with the most important contributions to the integral in Eq. (1)]. If this approximation is not valid [for instance, for total chain lengths comparable to or larger than c/(Γ + Γ 1D )], precursors are expected to appear even at the mirror configuration. The output field in this regime is wheret(ω) is the transmission coefficient in the non-Markovian limit. The transmission coefficient takes the form [50], Here, g(z, z , ω) = −[iΓ 1D (ω)/2]e ik 1D (ω)|z−z | with (g(z, ω)) j = g(z, z j , ω). In the non-Markovian regime, both the eigenvalues {λ ξ (ω)} and eigenvectors {v ξ (ω)} of g ij = g(z i , z j , ω) depend on frequency. Even if the qubits are in the mirror configuration at resonance (i.e., k 1D (ω 0 )d = nπ),t(ω) will generically have more than one simple pole, since k 1D (ω)d = nπ for most frequencies, so g ij will have more than one nonzero eigenvalue. Hence the contributions from the singularities oft(ω) to the output field will interfere and produce oscillations.

D. Role of spatial disorder and inhomogeneous broadening
The transients for high qubit numbers depend exclusively on the optical depth. For low qubit numbers, the transmitted intensity is not solely dictated by the OD. As shown in Fig. 4, transients change for different spatial configurations. The transients produced by the perfect array can be approximately recovered by averaging over many realizations with imperfect arrays with the same qubit number, as shown in Fig. 4(a). Even if the qubits have random positions, as in Fig. 4(b), the average over many realizations is qualitatively similar to the signal of a perfect array.
Next, we analyze the effect of inhomogeneous broadening on the transients. We consider that atom i is detuned from the central frequency of the input pulse by ∆ i , which is chosen randomly from a Gaussian distribution of mean µ = 0 and standard deviation σ. As shown in Fig. 5, the transients described by Eqs. (5) and (19) are robust against typical disorder levels found in experimental realizations for both large ( σ ∼ Γ [42]) and low qubit numbers (σ ∼ 0.01Γ 1D [44]).

IV. N QUTRITS
To observe the delay between the main pulse and its precursor, one can employ qutrits (three-level systems) under electromagnetic-induced transparency (EIT) conditions [7]. EIT has been used to observe precursors in both coherent [10] and single-photon [52] pulses propagating through dilute atomic clouds in free space. By coupling the excited state to a metastable level |s via a control field of Rabi frequency Ω c , a transparency window of width ∼ Ω 2 c / √ N Γ 1D Γ opens up and a pulse that is spectrally narrower than the window propagates without being absorbed or reflected, at a reduced group velocity v g = 2Ω 2 c d/Γ 1D . As shown in Fig. 6(a), a square pulse is also delayed and the precursor is measured before the main signal arrives at the detector. The system is described by the Hamiltonian where ∆ s = ω p −ω c −ω s is the two-photon detuning (with ω c being the frequency of the control field) andσ i es = |e i s i | is the coherence operator between the excited and the metastable states.
The transmission coefficient describing light propagation through this system is derived in an analogous manner to t N (ω) and reads [50] This response function has 2N poles at (complex) frequencies whereδ ξ ≡ λ ξ − iΓ /2. Proceeding in a similar way as in section III C, the transmitted intensity is Plugging in the poles and simplifying results in with Ω ξ = 4Ω 2 c +δ 2 ξ . The transmitted field after a large number of qutrits consists of an initial precursor, the main pulse, and a final transient. For large enough optical depth, the precursor is clearly separated from the (delayed) main pulse, as shown in Fig. 6(a).
At resonance, the intensity of the precursor is always the intensity of the main pulse (I d (t = t i ) = I 0 ), as can be seen from Eq. (27). Furthermore, the delay time of the main pulse can be inferred by estimating the timescale of decay of the second term. For simplicity, we consider ∆ = 0. The dependence in time of the contribution of mode ν to the second term of Eq. (27) is The most important contributions to Eq. (27) come from the most superradiant modes, for which we can neglect the shift and approximate λ ν = J ν − iΓ ν /2 ∼ −iΓ ν /2. Furthermore, we can assume Γ ν Γ , Ω c , which yields Plugging this expression into Eq. (28), we find where the first term has been neglected as Γ ν Ω c . The term with the slowest decay thus belongs to the most superradiant mode, for which Γ ν N Γ 1D . Hence, the time of arrival of the main pulse (i.e., the time at which the output signal stabilizes to I 0 ) scales as ∼ N Γ 1D Ω 2 c , which agrees with the delay of a monochromatic wave propagating through a dilute atomic cloud [10,53].
Under conditions of EIT, a single qutrit is enough to produce both a precursor and oscillations in the intensity after switching-off the input field, as shown in Fig. 6(b). Due to the coupling of the excited and metastable state, each qutrit contributes to two poles (of finite and opposite frequency) to the transmission coefficient. Right after switching off the field, there is no radiation as the excited state is unpopulated. The emission of light follows a slower trend, as light can only be emitted when the qutrits oscillate into the excited state. In the limit where Ω c Γ 1D , the oscillations have a period that scales as ∼ 1/ 4Ω 2 c − (Γ + Γ 1D ) 2 /4. The dips in the intensity correspond to times when the population of |e is minimal. For more qutrits, the oscillations are almost periodic if Ω c N Γ 1D , or resemble those of two-level systems (Fig. 2) in the opposite case.

V. CONCLUSIONS
In summary, waveguide QED constitutes a versatile platform for dispersion engineering, which can be employed to tailor the temporal shape of propagating photons. The transport of spectrally-broad photon pulses can be understood from the location of the poles of transmission coefficient in the complex plane, which correspond to the collective frequency shifts and decay rates arising from photon-mediated qubit-qubit interactions. For large qubit number, the response of the discrete system approaches that of a continuous medium, where the temporal oscillations in the intensity (arising from interference between frequency components of the original pulse) are fully determined by the optical depth N Γ 1D /Γ . In contrast, for low qubit number, there is a breakdown of the macroscopic response, and a single qutrit is enough to give rise to precursors (under EIT conditions). In this regime, the optical depth is no longer a good figure of merit, and the Purcell factor and the number of qubits, separately, play a significant role in the dynamics.
Precursors provide information about the number of qubits and their coupling separately, in contrast to the Beer-Lambert law that is obeyed by systems under continuous-wave illumination. For low atom numbers, this feature allows one to count how many atoms are coupled to a nanostructure, which is difficult to do with continuous-wave measurements [42]. Both the oscillation timescale for the transients and the delay between precursor and main pulse in the EIT regime increase with the ratio Γ 1D /Γ . Efficient coupling to the waveguide mode can be achieved in state-of-the-art experimental platforms, such as in superconducting qubits coupled to microwave transmission lines [44] (Γ 1D /Γ > 100) and quantum dots coupled to photonic crystal waveguides (Γ 1D /Γ 10 [54,55]). Acknowledgments -The authors acknowledge D. H. Phong, B. Macke, and S. Du for useful discussions. We gratefully acknowledge support from the Air Force Office of Scientific Research through their Young Investigator Prize (grant No. 21RT0751), the A. P. Sloan foundation, and the David and Lucile Packard foundation. S. C.-L. acknowledges support from the Chien-Shiung Wu Family Foundation. This work was supported in part by CONICYT-PAI grant 77190033, FONDECYT grant N • 11200192 from Chile.