Waveguide QED with Quadratic Light-Matter Interactions

Quadratic light-matter interactions are nonlinear couplings such that quantum emitters interact with photonic or phononic modes exclusively via the exchange of excitation pairs. Implementable with atomic and solid-state systems, these couplings lead to a plethora of phenomena that have been characterized in the context of cavity QED, where quantum emitters interact with localized bosonic modes. Here, we explore quadratic interactions in a waveguide QED setting, where quantum emitters interact with propagating fields confined in a one-dimensional environment. We develop a general scattering theory under the Markov approximation and discuss paradigmatic examples for spontaneous emission and scattering of biphoton states. Our analytical and semi-analytical results unveil fundamental differences with respect to conventional waveguide QED systems, such as the spontaneous emission of frequency-entangled photon pairs or the full transparency of the emitter to single-photon inputs. This unlocks new opportunities in quantum information processing with propagating photons. As a striking example, we show that a single quadratically-coupled emitter can implement a two-photon logic gate with unit fidelity, circumventing a no-go theorem derived for conventional waveguide-QED interactions.


I. INTRODUCTION
The study of light-matter interactions is one of the main research pillars of quantum science. The implementation of systems where quantum emitters interact strongly with confined modes of the electromagnetic fields has allowed us to achieve an unprecedented level of control over quantum degrees of freedom [1,2]. In the framework of cavity quantum electrodynamics (QED), the confinement of the electromagnetic field makes it possible to observe a coherent exchange of excitations between localized photonic modes and single quantum emitters [3]. Similarly, the coupling of quantum emitters to propagating fields can be strongly enhanced using waveguide structures, which confine photons to a onedimensional environment [4,5]. This setup, known as waveguide QED, has been realized in a variety of platforms such as atoms [6][7][8][9] or quantum dots [10,11] coupled to photonic waveguides, as well as superconducting qubits [12][13][14][15][16][17][18][19][20] coupled to microwave transmission lines. Waveguide QED structures have a great potential to implement building blocks of quantum networks [21,22], since propagating photons are ideal to transport flying qubits over long distances while emitters can provide the strong quantum non-linearity necessary for quantum information processing. Therefore, there has been intense theoretical [23][24][25][26][27][28][29][30] and experimental [31][32][33][34][35][36] research in the control and characterization of few-photon correlations generated by single quantum emitters, as well as in the implementation of photonic devices working at the few-photon level [37][38][39][40][41][42][43].
In the vast majority of platforms, quantum emitters couple linearly, via dipolar interactions, to photonic or phononic modes. Such interactions support only the exchange of individual quanta-e.g. an atom decays emit-ting a single photon-. Transitions involving multiple quanta appear as higher-order processes and are therefore strongly suppressed. Only recently, it has been shown how to implement nonlinear couplings where a quantum emitter interacts with localized bosonic modes via the direct exchange of two excitation quanta (e.g. an atom decays emitting a photon pair). Often dubbed two-photon interactions, these quadratic light-matter couplings have been proposed in cavity QED settings such as superconducting qubits non-linearly coupled to quantum microwave resonators [44][45][46] or nanomechanical oscillators [47]. Alternatively, quadratic couplings can also be effectively induced via parametric pumping and other simulation schemes intrinsic to superconducting circuits [48], hybrid spin-nanomechanical oscillators [49,50], trapped ions [51][52][53][54], and ultracold atoms [55,56]. Notice that non-dipolar [57] couplings have already been observed using superconducting artificial atoms.
In this work, we develop a quantum optics theory to describe a single quantum emitter interacting quadratically with the photons that propagate along a one-dimensional waveguide. We study this problem with a two-photon scattering theory based on a Wigner-Weisskopf approach and the Born-Markov approximation. We derive the general form of the scattering matrix, including semianalytical solutions for arbitrary photonic input states and full analytical solutions for Gaussian inputs. Applying this theory, we unveil observable features of the emitter's response that are fundamentally different with respect to conventional waveguide QED setups (see Fig. 1). These include (i) the spontaneous emission of correlated biphoton states, (ii) the strong interaction with spectrally narrow two-photon pulses, and (iii) full transparency to single-photon inputs. Finally, we show that these effects can be exploited in quantum information applications, designing a deterministic controlled-phase gate that acts on pairs of propagating photons with perfect fidelity. This result seems to contradict a famous no-go theorem for photonic gates [40,79,80]. However, it is the use of quadratic interactions that allows us to eliminate the wavepacket distortions that are intrinsic to the localized non-linearities induced by linear light-matter coupling [40]. The proposed gate is based on dual-rail encoding, which is of increasing relevance for superconducting quantum-computing applications [81]. All results discussed in this work can be implemented using state-of-the-art superconducting quantum emitters [44][45][46]48] interacting with propagating microwave fields or nanomechanical oscillators [47] and we expect that similar regimes will soon be achievable with solid-state [49,50] or atom-based [51][52][53][54][55][56] nanophotonic devices.
The work is structured as follows. In Sec. II, we develop a general theory for the scattering of pairs of photons by a single quantum emitter interacting quadratically with a one-dimensional waveguide. In Sec. III, we derive consequences from this theory, such as the twophoton spontaneous emission law and the two-photon scattering probabilities. We derive general bounds for the scattering cross-section valid for arbitrary photon-pair (biphoton) input states. We then consider specific examples for spontaneous emission and for the scattering of biphoton states with Gaussian and Lorentzian frequency distributions. In Sec. IV, we provide a striking example of application in quantum-information tasks, showing how to design a controlled-phase gate with unit fidelity. In Sec. V, we discuss possible physical implementations of the considered model. In Sec. VI, we provide conclusions and the outlook for future directions.

II. SCATTERING THEORY FOR TWO-PHOTON INTERACTIONS
This section introduces a new scattering theory for two-photon states impinging on a quadratically-coupled FIG. 1. Photon scattering on a quadratically-coupled twolevel emitter. The nonlinearity of the light-matter interaction is such that the emitter is completely transparent to one-photon pulses, while it strongly interacts with the multiphoton components of the input state. When a two-photon state is sent as input there are three allowed output channels: reflection, splitting, and transmission. The frequency distribution of the input photons plays a key role in determining each channel probability. The frequency of reflected and split photons are strongly anti-correlated, while the transmitted photons are positively correlated. quantum emitter. It describes the model's Hamiltonian (II A), the derivation of the scattering equations using the Wigner-Weisskopf formalism (II B), and the computation of the two-photon scattering matrix under the Markov approximation (II C).

A. Model
We consider a two-level quantum emitter with general quadratic coupling to a continuum of bosonic modes in a one-dimensional waveguide. Under the rotating-wave approximation, the system Hamiltonian can be written as (ℏ = 1), Here, ω 0 andσ + (σ − ) are the resonance frequency and the raising (lowering) operator of the emitter, respectively. The operatorâ µ ω (â µ † ω ) annihilates (creates) a waveguide photon with frequency ω and propagating direction µ ∈ {±}. We denote with g µµ ′ ωω ′ the function that gives the quadratic coupling strength between the emitter and the field frequency components. Note that Eq. (1) is over-parametrized. This is because the relation [â µ ω , (â µ ′ ω ′ ) † ] = δ µ,µ ′ δ ω,ω ′ implies that all vector elements appear twice in the summation and integration. Therefore, we can arbitrarily set the symmetric condition g µµ ′ ωω ′ = g µ ′ µ ω ′ ω , holding for any µ, µ ′ , ω and ω ′ . We will keep the coupling function as general as possible in the derivation, specifying it only in the examples, in order to identify the intrinsic properties of the quadratic coupling. The scattering theory that we develop is applicable to different experimental implementations of the model.

B. Wigner-Weisskopf ansatz
The Hamiltonian in Eq. (1) commutes with the oper-atorN = µ â µ † ωâ µ ω dω + 2σ +σ− . This is a continuoussymmetry generator associated with the conservation of the weighted number of excitations [51]. This conservation law implies that the whole dynamics can be described using a Wigner-Weisskopf ansatz wavefunction with a fixed number of (weighted) excitations as defined byN . In the vacuum sector or in the sector with individual photons, this dynamics is trivial: in both cases, the photonic states are decoupled from the emitter, which is transparent to individual photon states.
Consequently, we may focus on solving the scattering problem for input states with either two incoming photons and an emitter in the ground state, or an excited emitter in the vacuum. The most general such state, is constructed by creating excitations on top of |0⟩, the tensor product between the bosonic vacuum and the emitter's ground state. The two-photon Wigner-Weisskopf state is normalized as and it holds that C µµ ′ ωω ′ (t) = C µ ′ µ ω ′ ω (t) for any µ, µ ′ , ω, ω ′ . The Wigner-Weisskopf ansatz (2) describes a solution to the Schrödinger equation, i∂ t |Φ(t)⟩ =Ĥ |Φ(t)⟩, provided the coefficients satisfy a set of coupled, linear ordinary differential equations: It is convenient to reparameterize the wavefunction's element using the sum and difference of photon frequencies, ω = ω ′ + ω and ∆ = ω ′ − ω, as C µµ ′ ω∆ (t) and g µµ ′ ω∆ . At this stage, let us introduce the Markovian approximation, assuming that the emitter will spontaneously decay, releasing photons to the waveguide at a rate that is much slower than its intrinsic frequency, Γ ≪ ω 0 . In this limit one can formally solve Eq. (5) and replace the solution into Eq. (4) to obtain (see App. A 1): The total spontaneous emission rate is self-consistently defined as In the standard Markovian approach with linear interactions, the decay rate depends on the frequency of the emitted photons, although this dependency is usually weak or ignored [2]. In the case of quadratic lightmatter interactions, the decay rate depends on the coupling strength evaluated at the two-photon resonance, i.e. at ω = ω 1 + ω 2 = ω 0 , and integrated over frequency differences ∆. We can assume that the coupling strength does not change significantly with ω, at least in a band of frequencies around ω 0 , but we cannot fully eliminate the dependency on ∆, thus g µµ ′ ωω ′ ≡ g µµ ′ ∆ . This dependency may be factored into a product of emission rates γ µµ ′ ∈ R and an envelope u(∆) = u(−∆) ∈ C, i.e.
We consider ω 0 to be the dominant energy scale, that is we take ω 0 → ∞ for the upper integration limit in Eq. (7). This corresponds to assuming that the coupling function quickly decays with ∆. Then, we take u(∆) to be a bell-shaped function normalized to 1 in the l 2 -norm, i.e. ∥u∥ l 2 / √ 2 = 1 [82]. This convenient normalization implies that the function u(∆) encodes the differencefrequency dependence of the coupling strength, but it does not affect the total emission rate Γ. As we will see in the following, the explicit dependence on ∆ has highly non-trivial consequences on the scattering phenomenology and on potential applications.

C. Two-photon scattering matrix
The scattering matrix connects linearly the amplitude coefficient of the input photons C µµ ′ ω∆ (t 0 ) with the output field C µµ ′ ω∆ (t 1 ) at asymptotic times t 1 ≫ t 0 , when all the interaction with the quantum emitter has ceased. As explained in App. A 2, it is derived by integrating Eq. (5) using the initial and the final times as boundary conditions and equating the results. We obtain: where we definedC e (ω) = 1 √ 2π t1 t0 C e (τ )e iωτ . This quantity may be obtained from the Fourier transform of Eq. (6), i.e.
assuming that the emitter is initially in the ground state C e (t 0 ) = 0 and that there is a large temporal separation between the input and the output events. More precisely, we require that (t 1 − t 0 ) ≫ 1/Γ, so that the emitter is completely decayed at the time in which the output state is formally defined C e (t 1 ) = 0. For simplicity, in the following we consider the narrow-bandwidth limit and extend the integration region to ∞ in Eq. (10) (see App. A 2). This assumption is not strictly needed to obtain semi-analytical results.
For the sake of clarity, let us first write the scattering relations for a fixed propagation direction of the input field along the waveguide. In other words, we set C µµ ′ ω∆ (t 0 ) = C ω∆ (t 0 )δ µ λ δ µ ′ λ ′ for some fixed λ and λ ′ , where δ µ λ is a Kronecker delta. In this way, we can remove the summation in the direction index in Eq. (10). This assumption will be dropped in the final expression of the scattering matrix. After all these manipulations, the output and input states are related through As in the standard case of dipolar interactions, the output field is the interference between the input biphoton state (first line) and the two photons re-emitted by the atom (second line). Also as in standard waveguide QED, the Markovian approximation manifests itself in the Lorentzian dependency with respect to the total energy of the input field ω. However, differently from the dipolar case, for quadratic interactions the re-emitted field depends on the relative frequency ∆ through the interaction envelope u(∆). In the next sections, we will explore the consequences of this relationship. Let us now write a general relation for the two-photon scattering matrix.

Proposition 1. [Scattering matrix]
Consider a quadratic interaction between bosonic modes and a two-level emitter, defined by the Hamiltonian in Eq. (1), under Markovian approximation Γ ≪ ω 0 . Given a two-photon input state |Φ(t 0 )⟩ as in Eq. (2), the probability amplitudes at a time t 1 such that ( where a summation is implicit over the repeated indexes α and α ′ . The term S µµ ′ αα ′ (ω, ∆, ∆ ′ ) is a generic element of the scattering matrix, which can be expressed as where we defined In Eqs. (12) and (13), we identify that the integration with u(∆)u(∆ ′ ) plays the role of a projection onto a specific biphoton distribution u(∆). The action of the scattering matrix onto the input state is better understood if we explicitly separate the parallel and orthogonal components to u(∆) as with The scattering relation can then be rewritten as which is analogous of the scattering phase for a single photon in waveguide QED with linear Markovian interaction [23].

III. PARADIGMATIC CASES
In this section, we study the phenomenology of spontaneous emission (III A) and two-photon scattering (III B), analyzing the cross-section and the spectral features of the output channels. In the case of scattering of biphoton states, we provide general bounds on the scattering probabilities which hold for arbitrary frequency distributions of the input. We then focus on the Gaussian case, where we are able to give explicit analytical solutions for the scattering matrix and gather physical intuition on the features of the scattering problem. We also reproduce the analysis in the semi-analytical case of Lorentzian coupling function, which confirms that the results of the Gaussian case are valid for biphoton states of high experimental relevance.

A. Spontaneous emission
When we excite the quantum emitter in a waveguide without photons and both interact via Eq. (1), we expect the emitter to relax to the ground state emitting two photons. We can use the Wigner-Weisskopf theory [see Eq. (6)] to derive the frequency distribution of the emitted photons at asymptotic times. In the Markovian regime, the amplitude of the emitter's excited state decays exponentially with the rate Γ introduced in Eq. (7): The emitted photon wavepacket in frequency space at long times t 1 ≫ t 0 is recovered by inserting Eq. (18) into Eq. (9). This wavefunction, depends on the sum ω and on the difference of frequencies ∆. This profile resembles the Lorentzian lineshape that is characteristic of Markovian dynamics, where the wavepacket is centered around the total energy of ω of the incoming photons. However, the actual wavepacket is now shaped by the dependence of the interaction on the relative frequency of the photons g µµ ′ ∆ . This is a general result that can be applied to very different experimental scenarios.
Example: Gaussian and Lorentzian coupling functions-As an example, let us compute the spectrum of the spontaneously emitted photon pairs in Eq. (19), assuming two different frequency distributions for the coupling function. The first coupling function is Gaussian with zero mean and variance β 2 , i.e.
For comparison, we consider also a Lorentzian distributed coupling This gives us an expression for |C µµ ′ ω∆ (t 1 )| 2 via Eq. (19). In Fig. 2, we plot the emission rate for three different values of β, in the isotropic case where γ µµ ′ = Γ/4. Notice that, unless a non-isotropic emitter or a chiral waveguide is considered [83], there are three possible output channels, such that the photon pairs are emitted on the same direction or they can split. The frequency distribution is shown in Fig. 2 and it is the same for the three possible channels.
Summarizing, similarly to standard linear couplings the decay of the emitter excitation is exponential. Strong elements of novelty are: (i) Spontaneous decay of a single emitter results in the generation of a photon pair. (ii) The emitted biphoton state has a non-trivial frequency distribution, which is Lorentzian in the sum-frequency variable (main diagonal in Fig. 2, width set by the total decay rate Γ), while it reproduces the coupling function in the difference-frequency variable [anti-diagonal in Fig. 2, width set by the full width at half maximum, FWHM, of u(∆)]. (3) Interestingly, it is possible to generate frequency-correlated or anticorrelated photon pairs by changing the ratio between the total emission rate Γ and the FWHM of the coupling function u(∆).

B. Scattering of two-photon inputs
Let us discuss the phenomenology that appears in the scattering dynamics computed in Sec. II C. As sketched in Fig. 1, for two copropagating input photons λλ ′ = ++, there are three possible outcomes classified by the output directions: (i) reflection (µµ ′ = −−), (ii) transmission (µµ ′ = ++) and (iii) splitting (µ ̸ = µ ′ ), depending on whether both directions are preserved or reversed, or become opposite to each other. The respective scattering probabilities integrated over the complete frequency range are These values depend on the particular input field, but they exhibit general symmetries P −+ = P +− and constraints that can be derived from Prop. 1. Let us find an upper bound to the reflection probability under reasonable assumptions on the input state. This corresponds to maximizing Eq. (22) for µµ ′ = −−. Let us assume the atom to be initially in its ground state and the input field probability amplitude to be separable in ∆ and ω, i.e.
where f (ω) and h(∆) are complex functions normalized in the l 2 -norm, i.e. ∥f ∥ l 2 = 1 = ∥h∥ l 2 / √ 2. By using the Cauchy-Schwarz inequality in Eq. (12), we obtain The inequality in Eq. (24) is saturated if and only if h(∆) = e iφ u(∆) * , that is when the frequency distribution of the input field matches the coupling function g µµ ′ ∆ = g µµ ′ u(∆). Inserting Eq. (24) into Eq. (22) and using Holder's inequality, we obtain that the reflection probability R ≡ P −− is bounded by This inequality is saturated when the bandwidth δ of the input field f (ω)-measured as the FWHM or by some other means-is much narrower than the emitter's natural linewidth Γ ≫ δ.
We can also derive exact expressions for the maximal achievable rates of splitting S ≡ 2P −+ and transmission T ≡ P ++ processes, defined as While in waveguide QED with linear interactions a single emitter can perfectly reflect a single photon [37], the perfect two-photon reflection by a quadratically-coupled emitter is only possible when the splitting process is suppressed, i.e. γ −+ = 0 and γ ++ = γ −− . In the case of isotropic interactions γ µµ ′ = Γ/4, this is thus not possible. Here, even for optimal frequency-matched input h(∆) = e iφ u(∆) * , one obtains R max = 1/4, S max = 1/2 and T min = 1/4. Example: Gaussian wavepackets-Let us discuss a scattering experiment with Gaussian couplings and Gaussian input wavepackets that has analytical solutions and which provides insight that is extensible to other situations, e.g., the Lorentzian case shown in App. B. Without loss of generality, we use the isotropic Gaussian coupling function with variance β 2 from Eq. (20). The input state is formed by two Gaussian wavepackets with frequencies ω 1,2 , propagating along the same direction (++) with identical widths α 2 . The state is separable both on the individual photons' frequencies and the canonical variables {ω, ∆} We can now use the scattering matrix formalism to retrieve the expression of the output amplitude probabilities (see Prop. 1). Indeed, using Eq. (12) and Eq. (20), we find the frequency distribution of the field for the different output channels Let us now analyze the phenomenology described by this analytical result, focusing on the spectral properties of the output field and the scattering probabilities. The input field has been chosen as the product of two independent Gaussian FDs in ω and ∆, centered at ω0 (corresponding to the two-photon resonance condition), variance α 2 = (0.02ω0) 2 and FWHM= 2α √ 2 ln 2 = 0.047ω0. We consider a Gaussian distributed coupling [Eq. (20)] with β = α and we assume an isotropic spontaneous emission rate, i.e. γ µµ ′ = Γ/4. Focusing on the scattered field (rightmost plots), we see that the width in the sum-frequency variable ω (main diagonal) is proportional to Γ, while the width in the difference-frequency variable ∆ (anti-diagonal) is given by the input field variance α 2 (which has been set equal to that of the coupling function). Interestingly, we see that, when both photons are transmitted, their frequency distribution is correlated, while it is strongly anti-correlated for the reflection and splitting processes.
FIG. 4. Plot of the reflection probability, i.e. P −− , as a function of the ratio β/α. To obtain the plot, we considered a Gaussian input [Eq. (28)] centered in ω0 with variance α 2 and a Gaussian coupling with variance β 2 [Eq. (20)]. In addition, we assumed an isotropic spontaneous emission rate. We compare the reflection probability for different values of Γ and we show that the maximum value of P −− is reached when α = β and Γ ≫ α.
In Fig. 3, we plot the spectral distribution of the reflected (µ, µ ′ = −, −), split (µ, µ ′ = ±, ∓) and transmitted (µ, µ ′ = +, +) states, derived from Eq. (29) for two values of Γ. The distribution of the scattered/split photons is strongly anti-correlated in frequency, even though the input state was a separable one. A careful observation of the formulas reveals how the output field can be tailored by choosing the input state and by engineering the coupling function. In the sum-frequency variable ω [main diagonal in Fig. 3], the distribution has a Lorentzian shape whose width is established by the intensity of the light-matter coupling strength Γ. A larger value of the parameter Γ translates into a larger FWHM for the reflected frequency distribution. In the differencefrequency variable ∆ [anti-diagonal in Fig. 3], the distribution has a Gaussian shape whose width depends on the widths of the input distribution α and of the coupling function β.
In Fig. 4, we analyze how the total reflection rate depends on the input frequency distribution and the emitter properties. We plot P −− as a function of the ratio between the standard deviation β of the coupling function [see Eq. (20)] and the standard deviation α of the sumfrequency distribution of the input field [see Eq. (28)]. As expected, the scattering rate is maximal when the input frequency distribution matches the coupling function, a specific occurrence of the general results derived in Sec. II C. By probing different emission rates Γ, Fig. 4 also confirms that the upper bound on the reflection probability P −− = 1/4 is attainable when the width of the input frequency distribution is small compared to the emitter line width, Γ ≫ α.

C. Entanglement generation
From the phenomenological analysis of Sec. III A and Sec. III B, one can easily infer that two-photon states scattered or spontaneously emitted by a quadraticallycoupled emitter have non-trivial entanglement properties. In general, the scattered/emitted states are hyperentangled, i.e. non-classical correlations are present both in the position and in the frequency degrees of freedom. For the sake of brevity, we will focus here on entanglement generated in the (continuous) frequency variable. We consider a specific measurement setting relevant for practical experiments and applications such as quantum key distribution. Let us take the case of a spontaneously emitted two-photon state, which can be obtained by plugging the result of Eq. (19) in the definition of the state Eq. (2). We consider the setup shown in Fig. 5(a), that is we take the post-selected state |ψ PS ⟩ in which the photons are emitted in opposite directions, i.e.
where N is a renormalization factor and where we have re-expressed the output frequency distribution C(ω, ω ′ ) = C +− ω,∆ (t 1 ) as a function of the frequencies of the output photons ω and ω ′ . Bipartite entanglement in a continuous variable such as the frequency can be fully characterized with a Schmidt decomposition [84,85]. However, here we consider a simpler approach that is of direct experimental relevance. We identify two discrete frequency modes by applying frequency filters, defined by the functions f (ω). The effect of the filters is modeled by applying the following transformation to the state distribution: We choose the filters to have two narrow transmission windows, so that they are well approximated by the sum of two Dirac deltas centered around frequencies ω a and After the filters, the (post-selected state) can be written as where the first(second) quantum number in the ket notation encodes the frequency of left(right) propagating photons, that is |ω a , ω b ⟩ ≡ (â + ωaâ − ω b ) † |0⟩. We can now characterize the entanglement of the two-qubit state defined by Eq. (32). As entanglement measure we use the entanglement entropy defined as S(|ψ PS ⟩) = − Tr [ρ 1 log 2 ρ 1 ], where we defined the reduced density operator ρ 1 = Tr sys2 {|ψ PS ⟩ ⟨ψ PS |}, obtained by tracing out system 2 to the global state of Eq. (32). In Fig. 5(b) and (c) we analyze the behavior of the entanglement for different system parameters, considering the frequency distribution of Eq. (19) with the Lorentzian coupling function of Eq. (21). We take the filters to be centered around frequencies symmetric with respect to the atomic resonance, that is we fix ω a = ω 0 −δ and ω b = ω 0 +δ in Eq. (32). In Fig. 5(b) we show that the entanglement entropy grows with the frequency-offset δ, until it saturates to a plateau. This is consistent with the intuition gathered by looking at Fig. 2. Indeed, the output frequency distribution is separable close to the resonance condition, while positive or negative correlations can be observed far from resonance. To further characterize correlations, in Fig. 5(c) we show the entanglement entropy S(|ψ PS ⟩) as a function of the ratio between the width β of the coupling function and the total emission rate Γ. For β/Γ = 1 we observe a minimum of the entanglement entropy, as the frequency distribution is approximately separable [see Fig. 2]. For small values of β/Γ [ Fig. 5(d)] the state is strongly correlated, and in the limit β/Γ → 0 the state asymptotically approaches a Bell state lim β/Γ→0 |ψ PS ⟩ = |ψ − ⟩ ≡ 1 2 (|ω a , ω a ⟩ − |ω b , ω b ⟩). On the contrary, for large values of β/Γ [see Fig. 5(e)], the output photons are anticorrelated in frequency and the state asymptotically approaches the following Bell state lim β/Γ→∞ |ψ PS ⟩ = |ϕ + ⟩ ≡ 1 2 (|ω a , ω b ⟩ + |ω b , ω a ⟩). Therefore, quadratically-coupled atoms represent a compelling tool to generate frequency-entangled states, which are relevant for a variety of quantum-information applications [85][86][87][88]. In the following section, we provide another striking example of potential applications.

IV. IDEAL CONTROLLED-PHASE GATE FOR PROPAGATING PHOTONS
Let us now provide a compelling application of this phenomenology in quantum computing tasks. We consider a controlled-phase gate implemented on qubits defined with a dual-rail encoding, which is of increasing relevance for quantum computing with superconducting circuits [81]. In particular, we consider an encoding where each qubit is defined as a single photon distributed along two propagating modesâ andb, such that the logical states are encoded as |0⟩ L =â † |0, 0⟩ = |1, 0⟩ and |1⟩ L =b † |0, 0⟩ = |0, 1⟩. Regardless of the physical encoding, a c-phase (or control-Z) gate is defined by the transformation where c and s denote the control and signal logical qubits, respectively. It has been shown in Ref. [40] that, by first principles, it is impossible to implement this gate with unit fidelity using a single (point-like) two-level quantum emitter as the nonlinear element. The reason is that a distortion of the wavefunction frequency distribution is unavoidably introduced, as conflicting requirements on the input state should be imposed for the single-and two-photon scattering events. This no-go theorem has been derived assuming standard linear interactions. There are alternative proposals to circumvent such an intrinsic limitation, but they come at the cost of using multi-level artificial atoms [89], an array of many emitters [90], or an active time-dependent modulation of the coupling strength [91].
In the following, we propose a straightforward scheme that can overcome this limitation and implement a passive protocol for a c-phase gate with (in principle) unit fidelity, by using a single quadratically coupled two-level emitter per path. The scheme (sketched in Fig. 6) is composed of four semi-open waveguides (two for the control and two for the signal photons), each encoding a state of one of the two logical qubits. A balanced beamsplitter (BS), represented as a directional coupler in the sketch, connects the |1 c ⟩ and |1 s ⟩ paths. As the BS is traversed back and forth, it does not modify the logical input states. As a result, since quadratically coupled emitters are transparent to single-photon signals, it can be shown by direct inspection that the first three operations of Eq. (33) are correctly implemented. When the input state is |1 c , 1 s ⟩, after the first passage in the BS the state will be 1 √ 2 (|2 c , 0 s ⟩ + |0 c , 2 s ⟩), i.e. the photons will come out from the same path. The role of the BS is then to send both photons toward the same emitter if and only if the input state is |1 c , 1 s ⟩.
Let us analyze the process in which two photons |2 c , 0 s ⟩ (or equivalently |0 c , 2 s ⟩) scatter from the quadratically coupled emitters. We will compute which wavefunction of the two input photons maximizes the gate fidelity, studying the problem in frequency space. When the emitter is placed at the end of a semi-infinite waveguide, all photons are reflected, and the probability amplitudes of the output field are given by Prop. 1 with γ −− = γ −+ = 0, and the only non-zero coupling term is (7), it follows that Γ = γ ++ . As shown in Sec. III B, the optimal interaction is achieved when the input state wavefunction in frequency space matches the distribution of the coupling. Consequently, we assume a biphoton input state given by C ++ ω∆ (t 0 ) = f (ω)u(∆) * , where ∥f ∥ l 2 = 1. Following Prop 1, the output field is where we dropped the propagation-direction indexes for the sake of brevity. Note how the output field equals the input field C ++ ω∆ (t 0 ) = f (ω)u(∆) * times a Lorentzian (in square brackets) in ω, which in the limit of ω ≃ ω 0 approaches −1. We may use this limit to implement a c-phase gate, by imposing that the input state f (ω) is a Bell-shaped function, with FWHM δ, concentrated around the resonance frequency ω 0 . In the limit of narrow photons, Γ ≫ δ, the output field acquires the desired nonlinear phase C ω∆ (t 1 ) ≈ −e −iω0(t1−t0) C ω∆ (t 0 ), up to a phase e −iω0(t1−t0) common to all logical input states.
Let us now analyze the performance of the controlledphase gate, estimating the worst-case fidelity over all possible logical qubit states |ψ⟩ = a |0 c 0 s ⟩ + b |0 c 1 s ⟩ + c |1 c 0 s ⟩ + d |1 c 1 s ⟩. In a loss-less waveguide, the worstcase fidelity is defined as [40] whereÛ is the ideal c-phase gate andÊ are the actual operation experienced by the single and biphoton states. In our setup, only the two-photon state experiences a FIG. 6. Sketch of the circuit scheme that implements a controlled-phase gate using two-level quantum emitters. The scheme has been proposed in [40] to show intrinsic limitations to the achievable fidelity. Here, we show that in the case of quadratic interactions, there are no fundamental limits and the fidelity can in principle be arbitrarily close to 1. The origin of this advantage is that, for quadratic couplings, only the input state |1c, 1s⟩ interacts with the emitters. For all other input states, only single photons impinge on the emitters, which are fully transparent in the one-photon subspace.
FIG. 7. Worst-case gate logarithmic infidelity as a function of Γ, expressed in units of the input FWHM, in the case of a chiral waveguide with coupling g ++ ∆ = g ++ u(∆). The input amplitude probability is chosen as C ++ ω∆ = f (ω)u(∆) * . We compare two different relevant inputs: f (ω) real Normal distributed (blue) and f (ω) real Lorentzian distributed (orange). Both the input frequency distributions f (ω) are centered in ω0 and have the same FWHM. We see how in both cases the worst-case gate fidelity approaches unity. In particular, we observe a steeper slope for the Normal distributed f (ω) with respect to the Lorentzian one. nontrivial scattering. We may thus writê with |a| 2 + |b| 2 + |c| 2 + |d| 2 = 1 and 1 c1s given by Eq. (34). Inserting Eq. (36) into Eq. (35), we obtain a formula that only depends on the overlap 1 c 1 s 1 c1s In Fig. 7, we show the worst-case gate logarithmic infidelity as a function of Γ expressed in units of the input FWHM. We compared two different relevant input cases: f (ω) real Normal distributed and f (ω) real Lorentzian distributed. Both distributions are centered on the resonance frequency ω 0 and have the same FWHM. We observe a steeper slope for the Normally distributed f (ω) with respect to the Lorentzian one, due to the higher variance of the latter. Most importantly, in the limit of narrow inputs, FWHM≪ Γ, the worst-case gate fidelity approaches 1, thus overcoming the intrinsic limitation [40] found when considering linear interactions. This unanticipated result is of high conceptual value, as part of a prototypical analysis based of an effective model. An assessment of the fidelity achievable in practice should be done considering a specific device and keeping into account the corresponding main sources of dissipation and decoherence.

V. DISCUSSION ON PHYSICAL IMPLEMENTATIONS
Let us now provide comments on possible physical implementations of the considered model. There are two main approaches to implement quadratic interactions in the context of cavity QED. The first one consists in designing effective implementations, or analog quantum simulations, where an external driving is used to selectively activate pump-induced interactions. This approach can be applied to a broad variety of quantum platforms and it was considered to design quadratic interactions first for atomic systems [51][52][53][54][55][56] and then in solid state devices [48][49][50]. The second approach consists in engineering genuine quadratic couplings, such that the nonlinearity of the interaction is intrinsic and the exchange of excitations is not mediated by external drivings. Genuine implementations have been designed for superconducting quantum circuits [44][45][46] and nanomechanical resonators [47]. In particular, the most promising platform to implement the proposed controlled-gate is circuit QED [92], where decoherence and noise sources can be typically controlled within the 1%. In this framework, it was shown that it is possible to engineer a genuine quadratic qubit-cavity coupling [44] while fully inhibiting the linear coupling. It was also shown that a similar implementation can even be pushed into the ultrastrong coupling regime [45]. Such proposals are based on a flux qubit inductively coupled to a SQUID device, where the latter is operated in a weakly-nonlinear regime and is well approximated by a harmonic mode. A quadratic atom-cavity system linearly coupled to external waveguides. In the bad-cavity limit, the cavity mode can be adiabatically eliminated, obtaining an effective quadratic interaction between qubit and propagating modes.
Currently, it is therefore well understood how to design genuine implementations of quadratic interactions between a quantum emitter and localized harmonic modes. Going a step further, here we envision two ways in which these cavity QED settings could be directly generalized to implement quadratic couplings to propagating modes: (i) The most straightforward generalization is shown in Fig. 8(a) and consists in implementing an array of coupled cavities which supports propagating modes [93,94]. A quadratically-coupled emitter could be added using known designs into one of the resonators of the array. In this case, the quadratic coupling function to propagating modes would depend on the collective mode structure of the resonator array, and in the weak coupling limit we expect the effective dynamics to be Markovian [93] as assumed in this work. (ii) An alternative solution is shown in Fig. 8(b) and consists in considering a quadratic atomresonator system in the bad cavity limit, where the resonator dissipation rate is large compared to the coupling so that its modes can be adiabatically eliminated [95,96]. In this case, we expect an effective quadratic coupling with Lorentzian-like distribution between the artificial atom and the propagating modes of waveguides coupled to the leaky resonator.
These two solutions allow the implementation of the considered model with known results. A possible solution that would require further analyses consists in generalizing the superconducting circuit designs [44,45] to a native waveguide-QED setting. In the simplest case, such a configuration could consist in a superconducting transmission-line waveguide that is galvanically coupled to a flux qubit through a SQUID. More complex configurations could be designed considering SQUID-based metamaterials. A detailed analysis of any of these implementations should be performed considering the microscopic modeling of the specific platform and this is outside the scope of this paper.

VI. CONCLUSIONS AND OUTLOOK
This work introduces the study of quadratic lightmatter couplings in the field of waveguide QED. We have developed a general scattering theory with a broad range of applicability to solid-state and atomic quantum technologies. The predicted phenomenology is fundamentally different from what is observable with standard dipolar couplings and it bears great potential for quantum information applications. Remarkably, we show that a single quadratically coupled emitter can overcome a well-known no-go theorem and implement a controlled-phase gate on flying qubits with unit fidelity. Let us discuss the future perspectives in this novel research line, focusing first on the observation of an unconventional quantum phenomenology and then on the potential for quantum information applications.
In the framework of waveguide QED, quadraticallycoupled emitters can be embedded in very different setups and configurations. For instance, the placement of artificial atoms in arrays with sub-wavelength spacing can be used to enhance/inhibit their individual emission rate and scattering properties, both in standard and chiral waveguides. The emergence of collective effects arising from many quadratically-coupled emitters has so far been studied only in the context of cavity QED [71,72]. Nonlinear waveguide QED devices can also undergo interaction-induced phase transitions [97,98]. Alternatively, quadratically coupled emitters can be used as nonlinear mirrors that confine multiphoton states, by analogy with similar studies in the dipolar coupling regime [99]. Furthermore, the quadratic coupling in waveguide-QED setups can be pushed into the non-Markovian or even in the ultrastrong-coupling regime, requiring the development of new theoretical and numerical tools [13] and opening the door to new correlated phases and quantum phase transitions. Finally, nonlinear couplings can also be combined with the framework of giant atoms [100].
As demonstrated in this work, quadratic light-matter couplings are a compelling tool to process quantum information encoded in propagating photons. This possibility should be thoroughly addressed considering each specific experimental setting, with a detailed microscopic model of the system and realistic noise sources. The most interesting case is provided by superconducting quantum circuits, for which the dual-rail encoding is of high relevance for quantum computing [81]. Finally, besides gates in the dual-rail encoding, quadratic light-matter couplings can be used to deterministically generate propagating biphoton states or to implement highly non-trivial photon subtraction processes. The possibility of tailoring the spectral features of emitted/scattered biphoton states paves also the way to the generation of entanglement in the time-frequency domain, which is relevant for quantum computation [85] and sensing [86][87][88]  In this section, we provide a detailed derivation of the scattering matrix (Prop. 1) starting from the system of coupled differential equations for the amplitude probability coefficients in Eq. (4) and Eq. (5).

Markovian approximation
We start by solving formally the differential equation for the field amplitude probability coefficient C µµ ′ ωω ′ (t): Inserting this result in Eq. (4) allows us to decouple the system of equations and find a linear self-consistent differential equation for C e (t): The last term in Eq. (A2) can be further simplified. First, we change the integration variables defining ω = ω ′ + ω, ∆ = ω ′ − ω and dωdω ′ = dωd∆ 2 , so that Then we rewrite the integral as where we have used g µµ ′ ω∆ = g µµ ′ ω−∆ and we have defined the kernel So far the treatment is exact. Now we perform the Markovian approximation on the basis of two assumptions. First, we assume the coupling strength to be weak with respect to the emitter frequency ω 0 . First-order perturbation theory allows us to rewrite the solution of Eq. (A2) as C e (t) ≈ e −iω0t S e (t), where S e (t) is a slowlyvarying function of time. Second, we point out that, if the coupling strength has a weak dependence on ω, the Kernel K(t) is a rapidly-decaying function. Note that indeed, in the limit in which the coupling is flat in frequency, the Kernel is proportional to a Dirac delta in time. So, if the coupling strength is perturbative and sufficiently smooth with respect to ω, we can approximate S e (t) as a constant in the short time interval in which K(t) is non-vanishing. Accordingly, we can set S e (t − τ ) ∼ S e (t) in the time integral in Eq. (A4) and recombine the terms so that Note that for any specific form of the coupling strength, the integral in the last equation is now explicit and it can be evaluated by analytical or numerical means. However, to have a general analytical expression, we can take a further step and extend to infinity the upper limit of integration in Eq. (A6). This corresponds to assuming that the Kernel correlation length is zero (which is standard in the Markov approximation, as it is exactly true when the coupling strength is natively constant with respect to ω). Under this limit, performing the time-integration in Eq. (A4), we obtain FIG. 9. Two photon scattering, Lorentzian case. We compare the transmitted and reflected/split FDs for two different values of the total emission rate, Γ = 0.004ω0 for the first and Γ = 0.012ω0 for the second line. The FD for the reflection and splitting processes are identical. The FDs are shown in function of the two frequency variables ω, ω ′ . The input field has been chosen as the product of two independent Gaussian FDs in ω and ∆, centered at ω0 (corresponding to the two-photon resonance condition), variance α 2 = (0.02ω0) 2 and FWHM= 2α √ 2 ln 2 = 0.047ω0. We consider a Lorentzian distributed coupling as in [Eq. (21)] with β = 2α √ 2 ln 2, in order to have the same FWHM for the input field and for the coupling. We assume an isotropic spontaneous emission rate, i.e. γ µµ ′ = Γ/4.
where we have defined the total emission rate as We have also defined the Lamb shift which is a small shift in the emitter resonance frequency. Here and in the main text, we absorb δ e in the definition of ω 0 , as in most cases the experimental characterization of the emitter frequency already includes the Lamb shift. Now that all the terms in Eq. (A2) have been simplified, we are able to write a compact expression for the atom amplitude probability differential equation [Eq. (6) in the main text], iĊ e (t) = ω 0 C e (t) − i Γ 2 C e (t) where for convenience we have used the relation and g µµ ′ ω∆ = g µµ ′ ω−∆ , C µµ ′ ω∆ (t) = C µµ ′ ω−∆ (t).

Scattering states
As done in Eq. (A1), we can formally solve for C µµ ′ ωω ′ (t) setting the initial conditions for a time t 1 ≫ t 0 : Subtracting Eq. (A1) from Eq. (A13), we find an inputoutput relation for the field before (t = t 0 ) and after (t = t 1 ) the scattering event: where we switched to the frequency-sum and difference variables and we have used that g µµ ′ ω∆ = γ µµ ′ 2π u(∆). We now take the Fourier transform of Eq. (A11) [or Eq. (6)] and integrate by parts the left side, assuming that C e (t 0/1 ) = 0. We obtaiñ where we have used that, for t 1 ≫ t 0 , we can approximate t1 t0 e iτ (ω−ω) = 2πδ(ω − ω). Finally, for the output state we find In the main text, we consider the narrow-bandwidth limit and extend to infinity the range of the integration over ∆ in the previous expression. This corresponds to assuming that the input distribution is non-vanishing only in a region around a central frequency which is of the order of ω 0 (for two-photon interactions the input photons are resonant when the sum of their frequencies matches the emitter frequency). We also remind that u(∆) is bell-shaped and non-zero only for ∆ ≪ ω 0 .

Appendix B: Lorentzian case
In Sec. III B, we have presented results for the scattering of two-photon states in the Gaussian case, which can be analytically solved and allows one to understand the fundamental features of the two-photon scattering problem on a quadratically-coupled emitter. Let us now show that those features are not modified in the case of Lorentzian coupling functions [see Eq. (21)], which are more frequently encountered in realistic physical systems. This case can be solved semi-analytically using the results of Prop. 1, by numerically integrating Eq. (12). The results are presented in Fig. 9 and good qualitative agreement is found with the Gaussian case presented in Fig. 3.