Multiphoton Scattering Tomography with Coherent States

TL;DR

This paper presents a new experimental method for characterizing the scattering properties of quantum systems using coherent states and homodyne detection, enabling detailed analysis of elastic and inelastic scattering processes.

Contribution

It introduces a robust tomography technique that uses pulsed lasers and coherent states to efficiently probe both linear and nonlinear scattering matrices of unknown quantum systems.

Findings

Method is resilient to detector noise

Errors can be minimized by varying laser power

Effective for both elastic and inelastic scattering analysis

Abstract

In this work we develop an experimental procedure to interrogate the single- and multiphoton scattering matrices of an unknown quantum system interacting with propagating photons. Our proposal requires coherent state laser or microwave inputs and homodyne detection at the scatterer's output, and provides simultaneous information about multiple ---elastic and inelastic--- segments of the scattering matrix. The method is resilient to detector noise and its errors can be made arbitrarily small by combining experiments at various laser powers. Finally, we show that the tomography of scattering has to be performed using pulsed lasers to efficiently gather information about the nonlinear processes in the scatterer.

Figures

Equations

$$S_{p_1\dots p_n{k}_1\dots k_m}=⟨0⟩A_{p_1}\dots A_{p_n}U{A}_{k_1}^{†}\dots A_{k_m}^{†}0,$$

$$\ket{\Psi_{\rm in}}=|\alpha_{\vec{k}}\rangle=e^{-\frac{1}{2}|\alpha|^{2}}\mathrm{exp}\Big{(}\sum_{j=1}^{M}\alpha_{k_{j}}A^{\dagger}_{k_{j}}\Big{)}\ket{0},$$

$$⟨B_{p_1}\dots B_{p_n}⟩=N^{-n/2}⟨{\alpha }_{\vec{k}}|U^{†}{A}_{p_1}\dots A_{p_n}U|{\alpha }_{\vec{k}}⟩,$$

$$\alpha_{k_{j}}^{l,\vec{s}}=s_{j}e^{i\phi_{l}}|\alpha_{k_{j}}|,\;\;\;\;\left\{\begin{array}[]{l}s_{1}=1\\ s_{m\geq 2}=\pm 1\end{array}\right.,\;\;\;\;\phi_{l}=\frac{\pi}{M}l.$$

$$S_{p_1\dots p_n{k}_1\dots k_m}=\frac{N^{n/2}{e}^{|\alpha {|}^2}}{2^{M}M}\sum _{l=1}^{2M}\sum _{\vec{s}}\frac{F_n(l,\vec{s})}{{\prod }_{j=1}^m{\alpha }_{k_j}^{l,\vec{s}}}+{\epsilon }_m^{(1)},$$

$$\displaystyle\alpha_{k_{j}}^{\vec{s}}=s_{j}|\alpha_{k_{j}}|,\;\;\;\;\;\left\{\begin{array}[]{l}s_{1}=1\\ s_{j\geq 2}=\pm 1\end{array}\right.,$$

$$F_n(\vec{s})=⟨{\Psi }_{\mathrm{in}}(\vec{s})⟩U^{†}{B}_{p_1}\dots B_{p_n}U{\Psi }_{\mathrm{in}}(\vec{s}),$$

$$S_{p_1\dots p_m{k}_1\dots k_m}=\frac{N^{m/2}{e}^{|\alpha {|}^2}}{2^{M-1}}\sum _{\vec{s}}\frac{F_m(\vec{s})}{{\prod }_{j=1}^m{\alpha }_{k_j}^{\vec{s}}}+{\epsilon }_m^{(1)},$$

$$S_{p_1{k}_1}=⟨0⟩A_{p_1}U{A}_{k_1}^{†}0=e^{|\alpha {|}^2}\frac{⟨B_{p_1}⟩}{{\alpha }_{k_1}}+{\epsilon }_1^{(1)}.$$

$$S_{p_1{p}_2{k}_1{k}_2}=e^{|\alpha {|}^2}\frac{\left[F_2(1,1)-F_2(1,-1)\right]}{|{\alpha }_{k_1}||{\alpha }_{k_2}|}+{\epsilon }_2^{(1)},$$

$$S=\frac{bE(|\alpha {|}^2)-E(b|\alpha {|}^2)}{b-1}+{\epsilon }_m^{(2)},$$

$$A_{k_j}^{†}=\int \mathrm{d}{k}^{\prime }{\psi }_{k_j}(k^{\prime })a_{k^{\prime }}^{†},$$

$${\bar{\mathcal{S}}}_{p_1\dots p_n{k}_1\dots k_m}=⟨0⟩a_{p_1}\dots a_{p_n}U{a}_{k_1}^{†}\dots a_{k_m}^{†}0.$$

$${\bar{\mathcal{S}}}_{p_1{p}_2{k}_1{k}_2}=$$

$$+$$

$$S_{p_1\dots p_n,k_1\dots k_m}=\int \dots \int$$

$$\times \prod _{j=1}^m{\psi }_{k_j}(k_j^{\prime })\prod _{r=1}^n{\psi }_{p_r}^{\ast }(p_r^{\prime }).$$

$$G_{\sigma }(k^{\prime })$$

$$S_{p_1{k}_1}=e^{-\frac{(p_1-k_1{)}^2}{4{\sigma }^2}}\int \mathrm{d}{k}_1^{\prime }{G}_{\sigma }(k_1^{\prime }-\frac{[k_1+p_1]}{2})t_{k_1^{\prime }},$$

$$T_{p_1{p}_2{k}_1{k}_2}=2\sigma \sqrt{\pi }{e}^{-(p_1+p_2-k_1-k_2{)}^2/(8{\sigma }^2)}\times$$

$$\times \iiint \mathrm{d}{p}_1^{\prime }\mathrm{d}{k}_1^{\prime }\mathrm{d}{k}_2^{\prime }\mathcal{W}(p_1^{\prime },k_1^{\prime },k_2^{\prime }){\bar{\mathcal{T}}}_{p_1^{\prime }{p}_2^{\prime }{k}_1^{\prime }{k}_2^{\prime }}.$$

$$S_{p_1{k}_1}=\mathcal{O}(1),T_{p_1{p}_2{k}_1{k}_2}=\mathcal{O}({\sigma }^1).$$

$$t_{k_1}=\int \mathrm{d}{k}_1^{\prime }{\mathcal{K}}_{\sigma }(k_1^{\prime }-k_1)S_{k_1^{\prime }{k}_1^{\prime }},\text{mbox}where$$

$${\mathcal{K}}_{\sigma }(k)=G_{\sigma }(k)\sum _{q=0}^{\infty }\frac{(-1{)}^q}{2^{q}q!}{H}_{2q}\left(\frac{k}{\sigma }\right).$$

$${\bar{\mathcal{T}}}_{p_1{p}_2{k}_1{k}_2}=\frac{1}{\sqrt{\pi }\sigma }\iiint \mathrm{d}{\hat{k}}^{\prime }\mathrm{d}{\Delta }_p^{\prime }\mathrm{d}{\Delta }_k^{\prime }{T}_{{\hat{k}}^{\prime }-{\Delta }_p^{\prime },{\hat{k}}^{\prime }+{\Delta }_p^{\prime },{\hat{k}}^{\prime }-{\Delta }_k^{\prime },{\hat{k}}^{\prime }+{\Delta }_k^{\prime }}$$

$$\times {\mathcal{K}}_{\sigma }(\sqrt{2}[{\hat{k}}^{\prime }-\hat{k}]){\mathcal{K}}_{\sigma }({\Delta }_p^{\prime }-{\Delta }_p){\mathcal{K}}_{\sigma }({\Delta }_k^{\prime }-{\Delta }_k).$$

$$S_{nm}=⟨0|A^{n}U(A^{†}{)}^m|0⟩.$$

$$S_{p_1\dots p_n{k}_1\dots k_m}=⟨0|A_{p_1}\dots A_{p_n}U{A}_{k_1}^{†}\dots A_{k_m}^{†}|0⟩.$$

$$|{\Psi }_{\mathrm{in}}⟩=|\alpha ⟩=e^{-\frac{1}{2}|\alpha {|}^2}{e}^{\alpha A^{†}}|0⟩,$$

$$⟨A^n⟩=⟨{\Psi }_{\mathrm{out}}|A^n|{\Psi }_{\mathrm{out}}⟩=⟨\alpha |U^{†}{A}^{n}U|\alpha ⟩,$$

$$⟨Q⟩=e^{-|\alpha {|}^2}\sum _{t=0}^{\infty }\sum _{r=\lambda }^{\infty }|\alpha {|}^{r+t}{e}^{i\varphi (r-t)}{C}_{tr}.$$

$$C_{tr}=\frac{⟨0|A^t{U}^{†}QU(A^{†}{)}^r|0⟩}{t!r!}.$$

$$⟨Q⟩(l)=e^{-|\alpha {|}^2}\sum _{t=0}^{\infty }\sum _{r=\lambda }^{\infty }|\alpha {|}^{r+t}{e}^{i{\varphi }_l(r-t)}{C}_{tr}.$$

$$\frac{e^{|\alpha {|}^2}}{R}\sum _{l=1}^R\frac{⟨Q⟩(l)}{e^{im{\varphi }_l}}=\sum _{u=\lambda }^{\infty }\underset{d=-u+2\lambda }{\overset{u}{\text{scalerel}*\mathrm{\Sigma }\sum }}|\alpha {|}^u{T}_{ud}\sum _{l=1}^R\frac{e^{i{\varphi }_l(d-m)}}{R},$$

Full text

It is now possible to achieve strong and ultrastrong coupling between quantum emitters and propagating photons using superconducting qubits Astafiev et al. (2010); Hoi et al. (2011, 2013); Forn-Díaz et al. (2017), atoms Tiecke et al. (2014); Goban et al. (2015); Solano et al. (2017) or quantum dots Arcari et al. (2014) in photonic circuits, or even molecules in free space Hwang et al. (2009). This has motivated a stunning progress in the theory of single- and multiphoton scattering using wave functions Shen and Fan (2005) and the Bethe ansatz Shen and Fan (2007), as well as input-output theory Fan et al. (2010); Caneva et al. (2015), diagrammatic calculations Pletyukhov and Gritsev (2012); Laakso and Pletyukhov (2014); Hurst and Kok (2017), and path integral formalism Shi and Sun (2009); Shi et al. (2015). Very recently, the theory has even covered the ultrastrong coupling regime Shi et al. (2017). Experiments, however, cannot yet recover all the scattering information predicted by those studies, and are limited to comparing low-power coherent state transmission coefficients Astafiev et al. (2010); Forn-Díaz et al. (2017); Pechal et al. (2016), cross-Kerr phases Hoi et al. (2013), and antibunching Hoi et al. (2011). We therefore need an ambitious framework for reconstructing the complete one-, two-, or ideally any multiphoton scattering matrix. Such framework would allow studying the elastic Peropadre et al. (2013), and inelastic Goldstein et al. (2013) properties of quantum impurities in waveguides, quasi-particle spectroscopy Kurcz et al. (2014), interactions Gorshkov et al. (2010) in quantum simulators, and even characterizing all-optical quantum processors.

In this Letter we present a theoretical and experimental framework for estimating the scattering matrix

[TABLE]

which describes the transition amplitude from an input state of $m$ photons with generic quantum numbers $k_{1},\ldots,k_{m}$, to an asymptotic output state of $n$ photons with labels $p_{1},\ldots,p_{n}$ [cf. Fig. 1a]. Operator $U$ represents the evolution in the limit of infinitely long time, and $A^{\dagger}_{k}$ are generic input and output bosonic operators acting on the vacuum state $|0\rangle$. Our proposal assumes an experimental setup that injects coherent states and performs homodyne detection at the output of a generalized multiport beam splitter [cf. Fig. 1b]. Combining measurements with different input phases and amplitudes, we can approximate Eq. (1) with arbitrarily small error. The scheme is ideally suited for superconducting circuit and nanophotonic experiments, because all noise from amplifiers or detectors is canceled without previous calibration, similar to the dual-path method Menzel et al. (2010); Candia et al. (2014); da Silva et al. (2010).

In this Letter we also prove that the nonlinear contributions in the scattering are efficiently activated only for finite length input wave packets $A_{k}^{\dagger}$. However, standard deconvolution techniques applied to such experiments accurately and robustly reconstruct the scattering matrix in the monochromatic limit. We exemplify all these ideas for a two-level scatterer, whose scattering matrix is known Fan et al. (2010) and has only been probed at the single-photon level Astafiev et al. (2010); Hoi et al. (2011). This Letter is closed with a discussion on the generality of the protocol, and possible experimental implementations.

Scattering tomography protocol.–

Our presentation begins with the tomography architecture, continues with a discrete set of input states and a corresponding set of measurements, and closes with a general reconstruction formula for the scattering matrix. As shown in Fig. 1b, the setup consists of the quantum scatterer to be analyzed —or any active or passive optical medium—, a photonic channel that couples the light in and out of the scatterer, a signal generator to prepare the input states, and a multiport beam splitter that divides the output signal into $N$ independent measurement ports.

Our protocol requires a specific set of input states to probe the scatterer: coherent state wave packet inputs $|\alpha_{\vec{k}}\rangle$, created on top of the vacuum 111Throughout this work we assume that the scatterer has unique initial and final states, which are implicit in the definition of $\ket{0}$, but this restriction can be easily lifted. through a superposition of $M$ bosonic modes $A_{k_{j}}^{\dagger}$

[TABLE]

where $\alpha_{k_{j}}$ are complex weights and $|\alpha|^{2}$ the mean photon number 222The normalization of the coherent state (2) implies that the mean photon number is given by $|\alpha|^{2}=\sum_{j,j^{\prime}=1}^{M}\alpha_{k_{j}}^{\ast}\alpha_{k_{j^{\prime}}}f_{jj^{\prime}}$, where the commutators $[A_{k_{j}},A_{k_{j^{\prime}}}^{\dagger}]=f_{jj^{\prime}}\mathbb{I}$ are not orthonormal for wave packets.. These multimode coherent states can be prepared using a signal generator [cf. Fig. 1b], or combining laser pulses through beam-splitters [cf. Fig. 2]. While later examples assume wave packets $A_{k_{j}}^{\dagger}$ centered around momenta $k_{j}$, the formalism allows $k_{j}$ to label any set of quantum numbers: frequency, polarization, path, etc.

The output of the scatterer $|\Psi_{\rm out}\rangle=U|\Psi_{\rm in}\rangle$ is led through a balanced multiport beam splitter into $N$ homodyne detectors. At each output port, we filter outgoing photons with quantum numbers $\{p_{r}\}_{r=1}^{N}$, and measure the quadratures $\{X_{p_{r}},P_{p_{r}}\}$ to reconstruct the Fock operators $B_{p_{r}}=X_{p_{r}}+iP_{p_{r}}$. The nature of the beam splitter transformation is irrelevant, but all detectors should get a similar fraction of the scattered output, typically $B_{p_{r}}=N^{-1/2}A_{p_{r}}+\text{(N-1 vacuum inputs)}.$ Combining the homodyne measurements we estimate any correlation function of the form

[TABLE]

where the filtered labels $p_{1},\ldots,p_{n}$ are to the outgoing indices in the scattering matrix sector to be estimated (1). Note that any set size $n\leq N$ is possible.

The expectation value (3) is an analytic function $F_{n}[\vec{\alpha}]$ of the complex amplitudes $\alpha_{\vec{k}}$. Each order of the Taylor expansion of $F_{n}$ is determined by a different sector of the scattering matrix. As explained in the Supplemental Material Sup , the values $F_{n}[\vec{\alpha}]$ for different complex amplitudes —i.e. the correlation functions for different input states (2)—, can be added together, so as to cancel all terms except a desired sector of the scattering matrix and a small error. This is the basis for our protocols.

[TABLE]

For each input state, we measure the $2N$ amplitudes $B_{p_{r}}$ repeatedly, gathering statistics to reconstruct all correlations $F_{n}(l,\vec{s})=\braket{\Psi_{\rm in}(l,\vec{s})}{U^{\dagger}B_{p_{1}}\ldots B_{p_{n}}U}{\Psi_{\rm in}(l,\vec{s})}$, for $n=1\ldots N$. The scattering matrix is approximated as

[TABLE]

with a controlled error scaling as $\varepsilon_{m}^{(1)}=\mathcal{O}(|\alpha|^{2})$.

Note how the same set of measurements outcomes provides a simultaneous reconstruction of all scattering matrices from sizes $1\times 1$ up to $N\times M$.

Elastic scatterers.–

The reconstruction protocol simplifies when scattering conserves the total number of photons. This happens for emitters with one ground state and Jaynes-Cummings type interactions with $U(1)$ symmetry (no cyclic transitions). The scattering matrix (1) is exactly zero for $n\neq m$, quadrature measurements do not depend on global input phases and the total number of measurement setups reduces to $2^{M-1}$.

[TABLE]

with an error bounded by $|\varepsilon_{m}^{(1)}|\leq\mathcal{O}(e^{|\alpha|^{2}}-1)$.

Examples.–

Experiments with superconducting or optical qubits at low power are well described by a RWA Hamiltonian Fan et al. (2010), and we can apply Protocol 2. For a single photon we require only one input with arbitrary, but small, complex amplitude $\alpha_{k_{1}}$, and the measurement of one quadrature $B_{p_{1}}$, obtaining

[TABLE]

This formula includes the limit of state-of-the-art experiments Astafiev et al. (2010); Hoi et al. (2011); Forn-Díaz et al. (2017), where the transmission and reflection coefficients of single photons with momentum $k$ and fixed polarization, $S_{k,k}$ and $S_{-k,k}$, are recovered from the ratio between the input amplitude $\alpha_{k}$ of a monochromatic coherent beam, and the scattered amplitude $\braket{B_{\pm k}}$.

The reconstruction of the two-photon scattering matrix demands at least two measurement ports, $B_{p_{1}}$ and $B_{p_{2}}$, and a set of two input modes, $A_{k_{1}}^{\dagger}$ and $A_{k_{2}}^{\dagger}$. As shown in Fig. 2, an experiment could combine two independent pulses through a beam splitter, and then direct the scattering output to two homodyne measurement devices for estimating the correlations $\braket{B_{p_{1}}B_{p_{2}}}$. For a RWA model, we reconstruct

[TABLE]

using only two different input phases. If we cannot ensure $U(1)$ symmetry because of inelastic channels Goldstein et al. (2013), ultrastrong coupling Shi et al. (2017), external driving on the scatterer Bishop et al. (2009), etc., we need 8 input states with varying global phase $\phi_{l}=(\pi/2)l$, and the general reconstruction formula (5). However, the same measurements provide us with estimates for $S_{p_{1}k_{1}}$, $S_{p_{1}k_{1}k_{2}}$, $S_{p_{1}p_{2}k_{1}}$, and $S_{p_{1}p_{2}k_{1}k_{2}}$.

Arbitrary reconstruction error.–

There are four sources of error in our reconstruction protocol: (i) quantum fluctuations in quadrature measurements, (ii) detector noise, (iii) imperfect input preparation, and (iv) approximation error. The first source of error scales as $\mathcal{O}(\mathcal{N}^{-1/2})$ and can be decreased arbitrarily by increasing the number of repetitions $\mathcal{N}$ of the experiment. By design, our protocol is intrinsically resilient to the second source of errors, because detector noise averages out when combining odd powers of quadratures from different detectors —similar to Refs. Menzel et al. (2010); Candia et al. (2014); da Silva et al. (2010).

Imperfections in the preparation of relative phases $s_{j}=\pm 1+\delta s_{j}^{(\pm)}$, laser power $|\alpha|^{2}+\delta n$, and global phases $\phi_{l}=\pi l/M+\delta\phi_{l}$, add linear contributions to the error, $\varepsilon_{\rm sign}\sim|\delta s_{j}^{(\pm)}|$, $\varepsilon_{\alpha}\sim|\delta n|/|\alpha|^{2}$ and $\varepsilon_{\phi}\sim|\delta\phi_{l}|/|\alpha|^{m-1}$Sup . These errors can be controlled using standard calibration and phase stability techniques, and ensuring to work at moderate powers, $|\alpha|^{2}\sim 1$.

Indeed, a feature of our method is that we are not restricted to working at infinitesimal $|\alpha|^{2}$. Instead, we can combine different estimates of the scattering matrix, $E(|\alpha|^{2})=S-\varepsilon_{m}^{(1)}$, reconstructed from Eqs. (5)-(8) at different laser powers $|\alpha|^{2}$, to create a refined estimate with a higher order truncation error $\varepsilon_{m}^{(Z)}=\mathcal{O}(|\alpha|^{2Z})$. The simplest instance of this idea requires one extra estimate

[TABLE]

at a larger power $b>1$, to give $\varepsilon_{m}^{(2)}=\mathcal{O}(|\alpha|^{4})$. Higher order formulas can be derived analytically Sup , with error estimates $\varepsilon^{(Z)}_{m}=\mathcal{O}(|\alpha|^{2Z})$ that are strongly suppressed and allow working at $|\alpha|^{2}\gtrsim 1$. This is illustrated by Fig. 3a, where we plot an upper bound for $\varepsilon^{(Z)}_{m}(|\alpha|^{2})$ for elastic two-photon scattering, as a function of the smallest laser power $|\alpha|^{2}$, for different approximation orders $Z$. For $|\alpha|^{2}\sim 1$, we just need $Z\sim 10$ and $b=1.05$, to get error bounds $\varepsilon_{\rm b}^{(Z)}\sim 10^{-4}$, showing the potential of this method to estimate multiphoton scattering matrices.

The need of wave packets.–

We will now discuss the case in which $k_{j}$ describes the momentum degree of freedom. We will argue that our reconstruction protocol requires input states that are finite length wave packets

[TABLE]

built from normalized superpositions of plane waves $a_{k^{\prime}}^{\dagger}$, with $[a_{k},a_{k^{\prime}}^{\dagger}]=\delta(k-k^{\prime})$ and $\int|\psi_{k_{j}}(k^{\prime})|^{2}\mathrm{d}k^{\prime}=1$. The discussion below does not include other discrete degrees of freedom, which can be straightforwardly added Sup .

The use of pulsed light contrasts with existing theory, which computes the scattering matrix elements for well defined momentum modes Fan et al. (2010); Caneva et al. (2015); Shi et al. (2015), as in

[TABLE]

The reasons for studying the monochromatic $\bar{\mathcal{S}}$ are (i) the possibility of analytical calculations and that (ii) it reveals the underlying nonlinearity of the scatterer. Take, for instance, the two-photon scattering matrix for a two-level system, which can be decomposed as  Fan et al. (2010); Xu et al. (2013)

[TABLE]

with $\omega_{k}$ the photon dispersion relation and $c$ the velocity of light. The first two terms in Eq. (The need of wave packets.–) connect independent single-photon events $\bar{\mathcal{S}}_{pk}$, while the last one is the nonlinear contribution $\bar{\mathcal{T}}_{p_{1}p_{2}k_{1}k_{2}}$ that describes photon-photon interaction mediated by simultaneous interaction with the scatterer, such as the two-photon Kerr effect Hoi et al. (2013).

Interestingly, $S$ and $\bar{\mathcal{S}}$ are related through an integral,

[TABLE]

We will evaluate this for the two-photon scattering matrix in a tomography experiment using Gaussian pulses $\psi_{k}(k^{\prime})=G_{\sigma}(k^{\prime}-k)^{1/2}$, with

[TABLE]

We are specially interested in analyzing how the nonlinearity $\bar{\mathcal{T}}$ manifests in the monochromatic limit of negligible bandwidth $\sigma\to 0$. To do so, we focus on forward scattering $(k_{j},p_{r}>0)$ in a waveguide with linear dispersion relation $\omega_{k}=c|k|$, but this is easily extended Sup . The main result is that the measured two-photon scattering matrix $S_{p_{1}p_{2},k_{1}k_{2}}=S_{p_{1}k_{1}}S_{p_{2}k_{2}}+S_{p_{1}k_{2}}S_{p_{2}k_{1}}+iT_{p_{1}p_{2}k_{1}k_{2}}$ also splits into single- and two-photon contributions, which are given by the convolutions,

[TABLE]

The transmission coefficient $t_{k_{1}^{\prime}}$ appears in Eq. (17) as a consequence of energy conservation, $\bar{\mathcal{S}}_{pk}=t_{p}\delta(p-k)$, as well as the relation $p_{2}^{\prime}=k_{1}^{\prime}+k_{2}^{\prime}-p_{1}^{\prime}$ for the integration momenta in Eq. (18). Because the kernel $\mathcal{W}$ in Eq. (18) is a product of three Gaussians Sup and $\bar{\mathcal{T}}$ is typically a smooth and bounded function, any nonlinear contribution to the scattering experiment vanishes as we make the wave packet width $\sigma$ tend to zero,

[TABLE]

This argument can be extended to higher order processes, $T_{p_{1}\ldots p_{n}k_{1}\ldots k_{m}}\sim\sigma^{(m+n-2)/2}$ 333This derives from Eq. (15) by a purely dimensional analysis, assuming $\psi_{k}(k^{\prime})=\tilde{\psi}_{k}(k^{\prime}/\sigma)/\sigma^{1/2}$, and a monochromatic nonlinear contribution $\bar{\mathcal{S}}\sim\bar{\mathcal{T}}\delta(\sum_{r}\omega_{p_{r}}-\sum_{j}\omega_{k_{j}})$ that satisfies energy conservation as in Eq. (The need of wave packets.–)., illustrating the fact that nonlinear terms can only be activated when photons coexist in the scatterer, and the probability of this overlap tends to zero as the wave packet length $1/\sigma$ tends to infinity. We, therefore, conclude that an efficient reconstruction of the full scattering matrix for two or more photons requires working with finite duration wave packets.

Deconvolution formulas.–

Even if we need wave packets to get an experimentally measurable signal, we can still reconstruct the monochromatic properties from such experiments using standard deconvolution techniques Ulmer and Kaissl (2003). We illustrate this by deriving the single- and two-photon forward scattering coefficients $t_{k_{1}}$ and $\bar{\mathcal{T}}_{p_{1}p_{2}k_{1}k_{2}}$ from the measured $S_{k_{1}k_{1}}$ and $T_{p_{1}p_{2}k_{1}k_{2}}$. This requires inverting Eqs. (17)-(18), which can be done analytically for Gaussian wave packets Ulmer and Kaissl (2003); Ulmer (2010); Fang et al. (1994). For the single-photon transmission we obtain,

[TABLE]

The inverse kernel $\mathcal{K}_{\sigma}(k)$ contains Hermite polynomials $H_{q}(x)=(-1)^{q}e^{x^{2}}\partial_{x}^{q}(e^{-x^{2}})$ and produces a convergent series provided the function to deconvolve is $L_{1}$ integrable Ulmer and Kaissl (2003); Ulmer (2010). As discussed above, the single-photon reconstruction still works with monochromatic beams, recovering the state-of-the-art experimental formula $t_{k_{1}}=\lim_{\sigma\to 0}S_{k_{1}k_{1}}$.

The reconstruction of the two-photon scattering strength $\bar{\mathcal{T}}$ from the measured values of $T$ involves a three-dimensional deconvolution using a product of the same inverse kernels $\cal{K}_{\sigma}$ Ulmer (2010). Energy conservation imposes that the only nonzero elements of $\bar{\mathcal{T}}_{p_{1}p_{2}k_{1}k_{2}}$ have to be functions of the conserved average momentum $\hat{k}=(k_{1}+k_{2})/2=(p_{1}+p_{2})/2$, and the relative differences, $\Delta_{p}=(p_{2}-p_{1})/2$ and $\Delta_{k}=(k_{2}-k_{1})/2$. For these elements we get,

[TABLE]

As an illustration, we evaluate the measured scattering matrix $T$ and the reconstructed monochromatic version $\bar{\mathcal{T}}$, for a gedanken experiment with Gaussian pulses and a two-level scatterer. This problem admits an analytical solution Shen and Fan (2007); Fan et al. (2010); Sup with which we can test the reconstruction formulas (18) and (Deconvolution formulas.–). As shown in Figs. 3b and 3d, the measured matrix $|T|^{2}$ is broader than the monochromatic $|\bar{\mathcal{T}}|^{2}$, due to the convolution with the Gaussians. However, we have found that provided the wave packet size remains on order of the scatterer linewidth $\sigma\sim\gamma/c$, we can efficiently reconstruct $\bar{\mathcal{T}}$ from $T$. We exemplify this in Fig. 3c, where we show how cross sections of $|T|^{2}$ obtained for two different widths $\sigma=0.4\gamma/c$ (brown) and $\sigma=0.8\gamma/c$ (blue) both reconstruct the exact result for $|\bar{\mathcal{T}}|^{2}$ after the deconvolution (brown and blue dots). The fact that $\sigma\sim\gamma/c$ is a good compromise should not be a surprise, as this is the regime which maximizes the nonlinear effects and the coexistence of photons. In general, the gaussian deconvolution with kernel (20) allows us to efficiently reconstruct any smooth sector of the scattering matrix, provided the measured function to deconvolve is $L_{1}$ integrable Ulmer and Kaissl (2003); Ulmer (2010) and the order of the width $\sigma$ is chosen according to the ‘bandwidth’ of the specific scatterer.

Summary and outlook.–

This Letter introduced a tomography protocol for reconstructing the scattering matrix of a photonic field interacting with a quantum scatterer, using coherent states and correlated homodyne measurements. We have demonstrated that pulsed spectroscopy is needed to gather information about the nonlinear processes in scattering. This could remind the reader of two-dimensional pulsed spectroscopy methods in the optical and NMR realms Keusters et al. (1999); Cho (2008), but those constitute a time-resolved interrogation of the scatterer, whereas our protocol studies the asymptotic transformation (1) imparted by an optical medium in a propagating field.

While our protocol is inspired by recent progress in the fields of waveguide QED and nanophotonics, the idea, setup, and formulas can be used to probe any system that is in contact with a linear bosonic field. This includes not only superconducting qubits in strong-coupling Astafiev et al. (2010); Hoi et al. (2013) or ultrastrong-coupling 1D setups Forn-Díaz et al. (2017), but also studying single molecule emitters in three dimensions Hwang et al. (2009), or other extended optical media. The reconstruction protocol is so general that it does not require any a-priori knowledge of the quantum emitter, and can be applied in the presence of decoherence and dissipation. We believe that under such circumstances our protocol is optimal, but particular symmetries or a better understanding of the models can lead to substantial simplifications to be considered in future work.

Contents

• •

I.—Derivation of the general scattering tomography relations.

• •

II.—Arbitrary reconstruction error by combining multiple estimates of the scattering matrix.

• •

III.—Deconvolution with gaussian wave packets in dispersive channels.

• •

IV.—Imperfections in the input state preparation.

I I. Scattering tomography relations

In this section, we derive the central relations underlying our multiphoton scattering tomography protocol. We relate the components of an unknown unitary matrix $U$ with the measurement of certain photonic quadrature correlations at the output, when probing the system with coherent state inputs (see Fig. 4). In a first order approximation, a small reconstruction error requires attenuated coherent state inputs, but this restriction is lifted in Sec. II by combining various estimates at different input powers.

The derivation is split into two parts. First, in Sec. I.A, we show how to reconstruct the scattering matrix elements in the total photon number basis, namely

[TABLE]

These matrix elements describe processes involving $m$ incoming and $n$ outgoing photons, without resolving other photonic degrees of freedom. Then, in Sec. I.B, we extend the protocol to identify photons by a generic quantum number $k$, which accounts for any combination of external and/or internal photonic properties such as momentum, frequency, polarization, etc. In particular, we can reconstruct general scattering matrix elements of the form,

[TABLE]

describing processes from $m$ incoming photons with general quantum numbers $k_{1},\ldots,k_{m}$ to $n$ outgoing with $p_{1},\ldots,p_{n}$. As explained below, the injection and detection of photons with specific quantum numbers requires the use of a signal generator at the laser input and a multiport beam splitter at the measurement output (see also Fig. 4b).

I.1 A. Tomography of a scattering matrix in the photon number basis

We start by considering the situation sketched in Fig. 4a, where photons propagate along a one-dimensional (1D) channel and interact with an unknown quantum medium according to a scattering matrix $U$. As a result, any photonic input state $|\Psi_{\rm in}\rangle$ is transformed into the output $|\Psi_{\rm out}\rangle=U|\Psi_{\rm in}\rangle$, with the only restriction that the vacuum state $|0\rangle$ remains invariant, i. e. $U|0\rangle=|0\rangle$. Besides this last condition, the unitary $U$ is completely arbitrary, and we show in the following how to reconstruct its elements using homodyne detection.

First, our scattering tomography protocol requires the preparation of a coherent state input,

[TABLE]

where $A^{\dagger}$ is a bosonic creation operator of a photon in the channel, satisfying $[A,A^{\dagger}]=\mathbb{I}$, and $\alpha=|\alpha|e^{i\phi}$ is the complex vacuum displacement with module $|\alpha|$ and phase $\phi$.

Secondly, we require the measurement of photon output correlations of the form

[TABLE]

which can be determined to any order $n$ by measuring the output quadratures $X$ and $P$ via homodyne detection, with $A=X+iP$.

To derive the relation between the scattering matrix elements (22) and the correlations (25), let us first consider the measurement of a general operator $Q$ at the output, namely $\langle Q\rangle=\langle\alpha|U^{\dagger}QU|\alpha\rangle$. If we replace here the definition of the coherent state (24) and decompose $\alpha$ in terms of its module and phase, we obtain

[TABLE]

where the coefficients $C_{tr}$ read,

[TABLE]

In addition, $\lambda$ can only take two values: $\lambda=1$ in the special case $Q|0\rangle=0$ (since then $C_{t0}=0$), and $\lambda=0$ in any other case. In particular, for the choice $Q=A^{n}$ we can simply take $\lambda=1$, but we keep the derivation general.

Our aim now is to solve for the coefficients of the form $C_{0m}=\langle 0|U^{\dagger}QU(A^{\dagger})^{m}|0\rangle/m!$ from Eq. (26), since they correspond to the scattering matrix elements (22) when $Q=A^{n}$. As we show in the following, we can isolate the coefficients $C_{0m}$, for $m=\lambda,\ldots,M$, by preparing $R=2(M+1-\lambda)$ different coherent input states $|\alpha(l)\rangle$ of the form (24), and measuring $\langle Q\rangle(l)=\langle\alpha(l)|U^{\dagger}QU|\alpha(l)\rangle$ for each of them ($l=1,\ldots,R$). Importantly, we choose the different coherent state inputs $|\alpha(l)\rangle$ to have the same module $|\alpha|$, but different global phase $\phi_{l}$ as $\alpha(l)=|\alpha|e^{i\phi_{l}}$. As a result, we can replace $\langle Q\rangle\rightarrow\langle Q\rangle(l)$ and $\phi\rightarrow\phi_{l}$ in Eq. (26), obtaining a set of $l=1,\ldots,R$ independent equations for the unknowns $C_{tr}$, which read

[TABLE]

To see how to choose the phases $\phi_{l}$ in order to solve for the coefficients $C_{0m}$ in Eq. (28), let us multiply Eq. (28) by $e^{|\alpha|^{2}}e^{-im\phi_{l}}/R$ and sum over all values $l=1,\ldots,R$, obtaining the equation,

[TABLE]

after the convenient change of variables, $u=r+t$ and $d=r-t$. We have defined the new coefficients,

[TABLE]

and the sum over $d$, denoted with a bold sum symbol $\operatorname*{\scalerel*{\Sigma}{\sum}}_{y=a}^{b}$, corresponds to a sum where $y$ goes from $a$ to $b$ in steps of $2$. Choosing the phases as

[TABLE]

we see that the last factor on the right hand side of Eq. (29) becomes the discrete Fourier transform of a periodic Kronecker delta,

[TABLE]

with $q$ taking any integer value due to the function periodicity of size $R$. When replacing (32) into Eq. (29), we see that all non-vanishing terms must satisfy $d=m+qR$, with $q\in\mathbb{Z}$ and $d$ of the same parity as $u$. This allows us to select specific terms in the sum. In particular, if we set

[TABLE]

then all non-zero terms with $u\leq M+1$ in Eq. (29) can only take $d=m$ or equivalently $q=0$, as they all lie within the first period of the delta function. In addition, since $d$ is bounded by $-u+2\lambda\leq d\leq u$, then all terms with $u<m\leq M$ vanish, and our target coefficients $C_{0m}$ appear to lowest non-zero order $\sim|\alpha|^{m}$, for $m=\lambda,\ldots,M$, namely

[TABLE]

Here, the bounds of $q$ are given by $q_{\rm min}=-I[(u+m-2\lambda)/R]$ and $q_{\rm max}=I[(u-m)/R]$ where the function $I[x]$ rounds the number $x$ to its nearest smaller integer. From Eq. (34) it is clear that in the case of attenuated coherent state inputs, $|\alpha|\ll 1$, all terms on order $\sim|\alpha|^{m+2}$ or higher are further suppressed and we can solve for $C_{0m}$ with a small relative error. By formally solving for the coefficient $C_{0m}$ in Eq. (34), we obtain the main result of this subsection,

[TABLE]

which relates the desired matrix elements $\langle 0|QU(A^{\dagger})^{m}|0\rangle$, for powers $m=\lambda,\ldots,M$, with the measurements $\langle Q\rangle(l)=\langle\alpha(l)|U^{\dagger}QU|\alpha(l)\rangle$ at the output. In particular, when replacing $Q=A^{n}$ in Eq. (35), we access all the scattering matrix elements $S_{nm}=\langle 0|A^{n}U(A^{\dagger})^{m}|0\rangle$ in Eq. (22), describing $m=\lambda,\ldots,M$ incident photons to any number $n$ of outgoing photons in the channel $A$. As $Q|0\rangle=A^{n}|0\rangle=0$, we can take $\lambda=1$ and the number of measurements needed is $R=2M$.

Regarding the error $\epsilon_{m}^{(1)}$, we find a precise expression for it in terms of a series with even powers of the coherent state amplitude $|\alpha|$ as,

[TABLE]

From here it is clear that $\epsilon_{m}^{(1)}$ scales as $|\epsilon_{m}^{(1)}|\sim|\alpha|^{2}\ll 1$ in the case of attenuated coherent states. The coefficients $h_{vm}$ are independent of the coherent input power, and are given in general by,

[TABLE]

In addition, using Eqs. (27), (36) and (37), as well as the triangular inequality, we can show that the error (36) is strictly bounded by $|\epsilon_{m}^{(1)}|\leq\!\sum_{v=1}^{\infty}|\alpha|^{2v}|h_{vm}|$, where the bound for the coefficients reads

[TABLE]

and $|v\rangle=(A^{\dagger})^{v}/\sqrt{v!}$ are the photon number states.

It is worth highlighting the special case where we a priori know that the number of photons is conserved during the scattering, i. e. $[U,A^{\dagger}A]=0$, since then the tomography protocol is strongly simplified. In particular, only the diagonal components $S_{mm}$ are nonzero and a single measurement, $R=1$, is enough to determine them. To see this, notice that the conservation of photon number and $Q=A^{m}$ implies that $T_{ud}$ is only nonzero for $d=m$. Using this in Eq. (29) allows us to derive the same Eqs. (34) and (35), but now valid for $R=1$ and all $m=1,\ldots,\infty$, since in this case all coefficients appearing at lower order than $\sim|\alpha|^{m}$ are zero regardless the value of $R$. Another peculiarity of this elastic scattering case is that the bound for the error in Eq. (38) takes a very simple closed form. Again, using $Q=A^{m}$ and the conservation of photon number, we can set $U|v\rangle=e^{i\varphi_{v}}|v\rangle$ and $q=0$ in Eq. (38), and additionally using $A^{m}|v\rangle=\sqrt{v!/(v-m)!}|v-m\rangle$, we find a simple bound for the coefficients in this elastic case,

[TABLE]

As a result, the error $\epsilon_{m}$ in Eq. (36) is strictly bounded by a displaced exponential for all $m$, namely

[TABLE]

I.2 B. Tomography of a scattering matrix with resolution on photon degrees of freedom

In this subsection, we extend the protocol to measure scattering matrix elements with resolution on a generic quantum number $k$, describing any combination of external or internal photon degrees of freedom as in Eq. (23). First, in Sec. I.B.1. , we generalize Eq. (35) to determine matrix elements of the form $\langle 0|QUA^{\dagger}_{k_{1}}\ldots A^{\dagger}_{k_{m}}|0\rangle$, which describe input states of $m$ photons with different generic quantum numbers $k_{1},\ldots,k_{m}$. Then, in Sec. I.B.2. , we give details on a beam splitter setup to measure the output quadrature correlations $\langle B_{p_{1}}\ldots B_{p_{n}}\rangle$ (see Fig. 4b), and we show how to reconstruct the general scattering matrix elements (23) from these correlations.

I.2.1 1. Multimode coherent state input

First, we need a mechanism to inject photons, with various values of a combined quantum number $k$, at the input of the scatterer (see Fig. 4b). To this end we prepare a coherent state input $|\alpha_{\vec{k}}\rangle$ as in Eq. (24), but now for a multimode superposition $A^{\dagger}_{\vec{k}}$, defined for $M$ different components $\vec{k}=(k_{1},\ldots,k_{M})$ as

[TABLE]

Here, $A_{k_{j}}^{\dagger}$ for $j=1,\ldots,M$, are the creation operators of a photon with generic quantum number $k_{j}$, and $\xi_{j}$ is a projection vector normalized as $\sum_{j,v=1}^{M}\xi_{j}^{\ast}\xi_{v}[A_{k_{j}},A^{\dagger}_{k_{v}}]=\mathbb{I}$, such that $[A_{\vec{k}},A_{\vec{k}}^{\dagger}]=\mathbb{I}$. Notice that the modes $A_{k_{j}}$ do not satisfy standard commutation relations in general, allowing us to describe photon wave packets as inputs. For instance, for Gaussian wave packet modes as defined in Eqs. (12) and (16) of the main text, we obtain $[A_{k_{j}},A^{\dagger}_{k_{v}}]=e^{-(k_{j}-k_{v})^{2}/(2\sigma)^{2}}\mathbb{I}$.

The coherent state input corresponding to the multimode superposition (41) can be expressed more explicitly as

[TABLE]

where the weights on each component are given by $\alpha_{k_{j}}=\alpha\xi_{j}$, with $\alpha=|\alpha|e^{i\phi}$, and the mean photon number reads $|\alpha|^{2}=\sum_{j,v=1}^{M}\alpha_{k_{j}}^{\ast}\alpha_{k_{v}}[A_{k_{j}},A^{\dagger}_{k_{v}}]$. This multimode coherent state can be prepared using a signal generator (see Fig. 4b), which is readily implemented in microwave and optical photonic setups.

When we vary the global phase of the displacement weights as $\alpha_{k_{j}}(l)=|\alpha|e^{i\phi_{l}}\xi_{j}$, with $l=1,\ldots,R$ and $\phi_{l}$ defined in the previous subsection, we can apply the result (35) to the multimode superposition operator in Eq. (41), and thereby determine the matrix elements $\langle 0|QU(A_{\vec{k}}^{\dagger})^{m}|0\rangle$, with the associated error $\epsilon_{m}^{(1)}=|\alpha|^{2}$. However, this is not enough the determine matrix elements of the form $\langle 0|QUA^{\dagger}_{k_{1}}\ldots A^{\dagger}_{k_{m}}|0\rangle$, since when expanding the multinomial in $\langle 0|QU(A_{\vec{k}}^{\dagger})^{m}|0\rangle$, we see that it contains many unwanted terms, namely

[TABLE]

Notice that the sums over $n_{1},\ldots,n_{M}$ are constrained to values that satisfy $\sum_{j=1}^{M}n_{j}=m$, as demanded by the multinomial theorem.

Our aim in this subsection is therefore to extract from the above relation the term $\langle 0|QUA^{\dagger}_{k_{1}}\ldots A^{\dagger}_{k_{m}}|0\rangle$ and relate it to $\langle 0|QU(A_{\vec{k}}^{\dagger})^{m}|0\rangle$ by canceling all unwanted terms. We achieve this by performing additional measurements as in the previous subsection, but now we vary the relative phases in the multimode superposition $A_{\vec{k}}^{\dagger}(\vec{s})=\sum_{j=1}^{M}\xi_{j}^{(s_{j})}A^{\dagger}_{k_{j}}$ as $\xi_{j}^{(s_{j})}=e^{i\Theta_{j}^{(s_{j})}}|\xi_{j}|$. Here, $\Theta_{j}^{(s_{j})}$ for $s_{j}=1,\ldots R_{j}$ denote the $R_{j}$ different values that the phase of each component $j=1,\ldots,M$ can take, and we keep all the moduli $|\xi_{j}|$ fixed. We also use the shorthand notation, $\vec{s}=(s_{1},\ldots,s_{M})$, to arrange the values of the $M$ different indices. To see how to choose the phases $\Theta_{j}^{(s_{j})}$ in order to cancel the unwanted terms in Eq. (44), we divide this equation by $m!\ \xi_{1}^{(s_{1})}\ldots\xi_{m}^{(s_{m})}$ and sum over all independent values of $\vec{s}$, for $s_{j}=1,\ldots R_{j}$, and $j=1,\ldots,M$, obtaining

[TABLE]

Choosing $\Theta_{j}^{(s_{j})}=2\pi(s_{j}-1)/R_{j}$, for all $j=1,\ldots,M$, we can use the same property (32) as in the previous subsection and form periodic Kronecker deltas in Eq. (45),

[TABLE]

with $q_{j}$ and $t_{j}$ arbitrary integers. Importantly, due to the constrains $n_{j}\geq 0$ and $\sum_{j}n_{j}=m$, it is enough to choose $R_{j}$ as

[TABLE]

so that the deltas in Eq. (45) cancel all terms, except for the one with $n_{1}=\ldots=n_{m}=1$ and $n_{m+1}=\ldots=n_{R}=0$, and we obtain the desired relation,

[TABLE]

With the choice (47), we see that the relative phases take the simple values $\Theta_{j\geq 2}^{(1,2)}=(0,\pi)$ and $\Theta_{1}^{(1)}=0$, resulting only in sign changes of the vector components as, $\xi_{j\geq 2}^{(1,2)}=\pm|\xi_{j}|$ and $\xi_{1}^{(1)}=|\xi_{1}|$. Therefore, from now on (and also in the main text), we conveniently redefine the indices $\vec{s}=(s_{1},\ldots,s_{M})$ as

[TABLE]

such that the different values of vector components are simply given by

[TABLE]

Finally, replacing the result (35) into Eq. (48), we obtain a closed relation for the matrix elements $\langle 0|QUA^{\dagger}_{k_{1}}\ldots A^{\dagger}_{k_{m}}|0\rangle$, with $m=\lambda,\ldots,M$, and the various measurements $\langle Q\rangle(l,\vec{s})=\langle\alpha_{\vec{k}}(l,\vec{s})|U^{\dagger}QU|\alpha_{\vec{k}}(l,\vec{s})\rangle$, given by

[TABLE]

Here, we have used the new notation (49) for $\vec{s}$, and defined the weights $\alpha_{k_{j}}^{l,\vec{s}}$ of the coherent states $|\alpha_{\vec{k}}(l,\vec{s})\rangle$ in Eq. (43) as

[TABLE]

with $l=1,\ldots,R$, and $\phi_{l}$ as previously defined in Eqs. (31) and (33).

The error in Eq. (51) is also of the form,

[TABLE]

where the modified error coefficients $f_{vm}^{(1)}$ are given by

[TABLE]

with $h_{vm}(\vec{s})$, defined in Eq. (37), now depends on $\vec{s}$ via $A_{\vec{k}}(\vec{s})$ in the coefficients (27).

I.2.2 2. Homodyne detection at beam splitter outputs

Notice that Eq. (51) directly gives the desired scattering matrix elements,

[TABLE]

in the case we replace $Q=A_{p_{1}}\ldots A_{p_{n}}$. Using the resulting relation, however, demands the measurement of specific correlation functions at the scatterer’s output on channel $A$, namely

[TABLE]

To distinguish the contribution of photons with different generic quantum numbers $p_{1},\ldots,p_{n}$ in the above correlations (56), a filtering mechanism at the scatterer’s output is needed. Therefore, we connect the output of the scatterer to a multiport beam splitter, which divides the scattered signal into $N$ independent channels (see Fig. 4b). At each independent output port of the beam splitter, labeled by $r=1,\ldots,N$, we measure the quadratures $X_{p_{r}}$ and $P_{p_{r}}$, and reconstruct the field operator $B_{p_{r}}=X_{p_{r}}+iP_{p_{r}}$ for outgoing photons with quantum number $p_{r}$. When gathering enough statistics via homodyne detection, we can determine the correlations,

[TABLE]

which can be related to the correlations in Eq. (56) via the beam splitter transformation $U_{\rm BS}$.

Notice that the specific form of the beam splitter transformation $U_{\rm BS}$ is irrelevant, as long as each output port $r$ gets a similar fraction of the scattered signal $B_{p_{r}}\sim A_{p_{r}}$ and that all other $N-1$ input ports are in vacuum (see Fig. 4). As a practical example, let us consider the transformation for a balanced multiport beam splitter, for which the photonic amplitude operators at each independent output, $r=1,\ldots,N$, read

[TABLE]

Here, $\Upsilon_{p}^{r^{\prime}}$ with $r^{\prime}=1,\ldots,N-1$ denote the annihilation operators on the $N-1$ vacuum input ports different than channel $A$. As these extra channels are independent between each other and with channel $A$, they satisfy the commutation relations, $[\Upsilon_{k}^{r^{\prime}},\Upsilon_{p}^{r}{}^{\dagger}]=\delta_{rr^{\prime}}\delta_{pk}$, and $[A_{k},\Upsilon_{p}^{r}{}^{\dagger}]=[U,\Upsilon_{p}^{r}]=0$, implying $[B_{p_{r}},B_{p_{r^{\prime}}}^{\dagger}]=\delta_{rr^{\prime}}$.

Using expression (58) in Eq. (57), we see that the required correlations (56) can be accessed by measuring at the outputs of the beam splitter. In fact, both correlations are just proportional to each other:

[TABLE]

Finally, we just replace $Q=A_{p_{1}}\ldots A_{p_{n}}$ into Eq. (51) and use the relation (59) to obtain the general scattering matrix elements $S_{p_{1}\ldots p_{n}k_{1}\ldots k_{m}}$ in terms of the measurable correlations $F_{n}(l,\vec{s})=\langle\alpha_{\vec{k}}(l,\vec{s})|U^{\dagger}B_{p_{1}}\ldots B_{p_{n}}U|\alpha_{\vec{k}}(l,\vec{s})\rangle$, as

[TABLE]

This is the main result of the paper and is shown in Eqs. (5) and (8) of the main text, for the general ($R=2M$) and elastic ($R=1$) cases, respectively.

The truncation error $\varepsilon_{m}^{(1)}$ appearing in the final result (60) has a closed analytical upper bound in the case that the scattering matrix $U$ conserves the number of photons (Protocol 2 of the main text). From Eqs. (53)-(54) we see that this problem is equivalent to find a closed bound for the coefficients $|h_{vm}(\vec{s})|$. Therefore, we replace $Q=A_{p_{1}}\ldots A_{p_{m}}$ in Eq. (38) and use Eq. (58), to re-express the bound in terms of the $B_{p_{r}}$ operators as

[TABLE]

Since $U$ conserves the total number of photons, we can set $q=0$ and $U|v\rangle=e^{i\varphi_{v}}|v\rangle$ in Eq. (61), and additionally using the Cauchy-Schwarz inequality, we obtain

[TABLE]

The next step is to decompose the total number state $|v+m\rangle$ into components along all $N$ channels of the interferometer, and using the fact the the sum of photons in all channels is fixed, $N=\sum_{r=1}^{N}B_{p_{r}}^{\dagger}B_{p_{r}}$, we can bound the numerator in Eq. (62) by $\sqrt{(v+m)!/v!}$. Using this result in Eq. (62) we obtain a simple bound for the $h_{vm}$ coefficients similar to Eq. (39) as $|h_{vm}(\vec{s})|\leq N^{m/2}/v!$, and replacing this into Eq. (54), we find the bound for the $f_{vm}$ coefficients as well,

[TABLE]

For simplicity, let us consider the case that we prepare the multimode coherent state inputs with equal projections, $|\xi_{j}|=|\xi|$. The normalization condition, stated below Eq. (43), implies $1/|\xi|=\sqrt{\sum_{j,v=1}^{M}e^{-(k_{j}-k_{v})^{2}/(2\sigma)^{2}}}\leq M$, in the case of gaussian wave-packets as in the main text. Using this and $N=M$, we find a strict bound for the coefficients that is easy to evaluate,

[TABLE]

and replacing this into Eq. (53), the strict bound for the error reads,

[TABLE]

II II. Arbitrary reconstruction error by combining multiple estimates

An intrinsic error that limits the precision of our tomography protocol is the reconstruction error $\varepsilon_{m}^{(1)}=\sum_{v=1}^{\infty}|\alpha|^{2v}f_{vm}$. In this section, we show how to reduce this error by performing the scattering tomography protocol in Eq. (60) at different laser powers and combining the outcomes in a clever way. Importantly, this allows us to perform the tomography at non-perturbative mean photon numbers $|\alpha|^{2}\sim 1$, lifting the requirement of attenuated coherent state inputs.

A first order estimate $E$ of the scattering matrix $S$ is given by the outcome of our protocol in Eq. (60) as,

[TABLE]

with an error $\varepsilon_{m}^{(1)}=S-E(|\alpha|^{2})=\sum_{v=1}^{\infty}|\alpha|^{2v}f_{vm}$, scaling as $\varepsilon_{m}^{(1)}=\mathcal{O}(|\alpha|^{2})$. Interestingly, if we consider another first order estimate measured at a larger power $E(b|\alpha|^{2})$, with $b>1$, we can combine it with $E(|\alpha|^{2})$, and cancel the error term of order $\sim|\alpha|^{2}$. As a result, we obtain a second order estimate,

[TABLE]

with a higher order approximation error,

[TABLE]

that scales as $\varepsilon_{m}^{(2)}=\mathcal{O}(|\alpha|^{4})$. Here, the modified error coefficients read $f^{(2)}_{vm}=f_{vm}[(b-b^{v})/(b-1)]$. This process can be iterated up to arbitrary order using the recursion relation,

[TABLE]

which gives the estimate of order $\mu$ from two estimates of a lower order: $E^{(\mu-1)}(|\alpha|^{2})$ and $E^{(\mu-1)}(b|\alpha|^{2})$. In general, a $Z$-order estimate $E^{(Z)}(|\alpha|^{2})$ has an error $\varepsilon_{m}^{(Z)}=\mathcal{O}(|\alpha|^{2Z})$ given by

[TABLE]

where the $Z$-order coefficients can expressed in terms of $f_{vm}$ in Eq. (54) as

[TABLE]

The final step is to express a general $Z$-order estimate $E^{(Z)}(|\alpha|^{2})$ in terms of first order estimates only, as these are the ones that we measure in practice via (66). Using $Z-1$ times the recursion relation (69), it is straightforward to check that $E^{(Z)}(|\alpha|^{2})$ can always be expressed as a linear combination of $Z$ different first order estimates $E(x_{q})$, measured at increasing laser powers $x_{q}=b^{q-1}|\alpha|^{2}$, with $q=1,\ldots,Z$ and $b>1$, as

[TABLE]

Here, the $Z$ coefficients $w_{q}^{(Z)}$ are determined recursively from Eqs. (67) and (69). As an example, for $Z=1,\ldots,4$, the coefficients are explicitly given by

[TABLE]

In summary, the elements of the scattering matrix $S$ can be estimated with a $Z$-order error $\varepsilon_{m}^{(Z)}=\mathcal{O}(|\alpha|^{2Z})$ given in Eqs. (70)-(71) by combining $Z$ first order estimates (66), obtained at incresing laser powers as

[TABLE]

Importantly, the $Z$-order error $\varepsilon^{(Z)}_{m}=\mathcal{O}(|\alpha|^{2Z})$ is strongly suppressed by its prefactor, which allows one to perform the scattering tomography protocol at non-perturbative laser powers $|\alpha|^{2}\gtrsim 1$. We show this explicitly by numerically evaluating the upper bound of the $Z$-order error $\varepsilon_{\rm b}^{(Z)}$ in the case of elastic scattering $m=n$, which is obtained from Eqs. (65), (70), and (71) as

[TABLE]

In Fig. 5a-b, we plot $\varepsilon^{(Z)}_{\rm b}/M^{3m/2}$ in logarithmic scale, as a function of the minimal mean photon number $|\alpha|^{2}$ and the number of estimates $Z$, for for two values of the factor $b$. For $b=1.05$, we indeed get a good error bound of $\varepsilon^{(Z)}_{\rm b}/M^{3m/2}\sim 10^{-5}$ for $|\alpha|^{2}\sim 1$ and a moderate number of estimates $Z\sim 10$. Figure 3a of the main text shows the same bound $\varepsilon^{(Z)}_{\rm b}$ and the same parameters, but evaluated in the specific case $M=m=2$.

In general, the error suppression is stronger for smaller $b$, as it becomes clear by comparing Figs. 5a and 5b. Experimentally, the value of this factor $b$ is restricted by the precision to increase the laser power as $x_{q}=b^{q-1}|\alpha|^{2}$. Finally, notice that for a given starting mean photon number $|\alpha|^{2}$ and a fixed value of the $b$ factor, there will be an optimal number of estimates $Z$ that gives the maximum supression of the error. For $Z$ higher than the optimal, the error $\varepsilon^{(Z)}_{m}$ will actually increase due to its prefactors (see Fig. 5b).

III III. Deconvolution with gaussian wave packets in dispersive channels

To be sensitive to nonlinear multiphoton scattering events, our scattering protocol requires wave packet input modes in the case the quantum number $k_{j}$ describes a continuous degree of freedom such as momentum or frequency. Therefore, in the following analysis we want to distinguish between continuous and discrete photon degrees of freedom, which we label by $k_{j}$ and $\lambda_{j}$, respectively. For instance, photon polarization and/or path can be included in the $\lambda_{j}$ quantum numbers. Using this separation of labels, we assume that the generic bosonic modes in Eq. (43) are prepared as

[TABLE]

where $\psi_{k_{j}}(k^{\prime})$ is the square-normalized wave packet profile for the continuous variable $k^{\prime}$, centered at $k_{j}$, and $a_{k^{\prime}}^{\lambda_{j}}$ describe monochromatic photons with quantum numbers $k^{\prime}$ and $\lambda_{j}$, satisfying standard communation relations $[a_{k}^{\lambda_{j}},(a_{k^{\prime}}^{\lambda_{j}^{\prime}})^{\dagger}]=\delta_{\lambda_{j}\lambda_{j^{\prime}}}\delta(k-k^{\prime})$.

The choice of input modes in Eq. (79) implies that our scattering protocol gives us direct access to the ‘wave packet’ scattering matrix elements only, given in general by

[TABLE]

As explained in the main text, the scattering matrix elements in the monochromatic Fock basis, $\bar{\mathcal{S}}_{p_{1}\ldots p_{n}k_{1}\ldots k_{m}}^{\lambda_{1}^{\prime}\ldots\lambda_{n}^{\prime},\lambda_{1}\ldots\lambda_{m}}=\langle 0|a_{p_{1}}^{\lambda_{1}^{\prime}}\ldots a_{p_{n}}^{\lambda_{n}^{\prime}}U(a_{k_{1}}^{\lambda_{1}})^{\dagger}\ldots(a_{k_{m}}^{\lambda_{m}})^{\dagger}|0\rangle$, can still be accessed by using Eq. (79) and inverting the integral relation,

[TABLE]

Therefore, one needs to perform a deconvolution on the continuous variables $k_{1},\ldots,k_{m},p_{1},\ldots,p_{n}$ for each combination of the discrete variables $\lambda_{1},\ldots,\lambda_{m},\lambda_{1}^{\prime},\ldots,\lambda_{n}^{\prime}$, which in turn requires the measurement of the wave packet matrix elements for all convinations of quantum numbers $S_{p_{1}\ldots p_{n},k_{1}\ldots k_{m}}^{\lambda_{1}^{\prime}\ldots\lambda_{n}^{\prime},\lambda_{1}\ldots\lambda_{m}}$. As we are mainly interested in illustrating how to perform a single deconvolution, in the remainder of the section, as well as in the main text, we assume that all discrete quantum numbers are equal and we absorb them in the definition of the photon operators $A_{k_{j}}$ and $a_{k_{j}}$. This simplifies the equations considerably, and we just have to keep in mind that to describe additional internal degrees of freedom of photons, we have to repeat the same deconvolution procedure for each combination of of the discrete labels $\lambda_{1},\ldots,\lambda_{m},\lambda_{1}^{\prime},\ldots,\lambda_{n}^{\prime}$.

In this section, we show how to analytically solve the deconvolution integral for the monochromatic single- and two-photon scattering matrices in Eq. (The need of wave packets.–) of the main text, which corresponds to the more general Eq. (III) in the case the continuous degree of freedom $k_{j}$ corresponds to momentum or frequency, and all photons have the same discrete degrees of freedom. To do so, we assume wave-packets with gaussian shape and a possibly dispersive photonic channel, as shown below.

In Sec. III.C, we show the analytical expressions for the single- and two-photon scattering matrices used in Figs. 3(b)-(d) of the main text to test the gaussian deconvolution formulas.

III.1 A. Gaussian deconvolution in a dispersive channel

As we want to discuss photons in a dispersive channel, the deconvolution integrals turn out to be simpler when working with frequency modes instead of momentum. Therefore, in this subsection we derive all deconvolution formulas in terms of frequency-direction modes, and in the next subsection we map the results to momentum, connecting to the specific formulas shown in the main text.

In the special case that the dispersion relation is symmetric, $\omega_{k}=\omega_{-k}$, and the group velocity is antisymmetric, $v_{k}=-v_{-k}$, the propagating momentum modes $a_{k}$ can be unambiguously related to frequency-direction modes $a^{\rho}_{\omega}$, as

[TABLE]

Here, the mapping assumes that photons with $k>0$ propagate to the right ($v_{k}>0$), while photons with $k<0$ propagate to the left ($v_{k}<0$), and we identify these directions with the quantum number $\rho={\rm sign}(k)=\pm$. In addition, the frequency $\omega$ of the photons has a one-to-one correspondence with the wavenumber $|k|$ via the invertible function $\omega=\omega_{|k|}$. Note that this relation (82) is not valid for $k=0$, since the group velocity necessarily vanishes $v_{0}=0$. Nevertheless, we will always consider propagating wave packets (79) with negligible components on this $k=0$ mode.

The wave packet operators in Eq. (79) can be expressed in terms of frequency-direction modes as

[TABLE]

where the corresponding wave packet profile is related to its momentum counterpart by

[TABLE]

In particular, we choose wave packets with a gaussian profile defined as

[TABLE]

such that the wave packet mode $A_{\omega}^{\rho}{}^{\dagger}=\int d\omega^{\prime}G_{\tilde{\sigma}}(\omega^{\prime}-\omega)^{1/2}a^{\rho}_{\omega^{\prime}}{}^{\dagger}$ is built exclusively from monochromatic frequency modes $a^{\rho}_{\omega^{\prime}}$ propagating in the same direction $\rho$. We also assume that the central frequency of the wave packet $\omega$ is much larger than its width, $\tilde{\sigma}\ll\omega$, to ensure a vanishing component on the non-propagating mode $\omega=0$.

Using all these considerations, we can re-express the general integral equation (15) of the main text, in terms of frequency-direction modes, as

[TABLE]

where $\bar{\mathcal{S}}_{\omega_{1}\ldots\omega_{n}\nu_{1}\ldots\nu_{m}}^{\rho_{1}\ldots\rho_{n}\eta_{1}\ldots\eta_{m}}=\langle 0|a_{\omega_{1}}^{\rho_{1}}\ldots a_{\omega_{n}}^{\rho_{n}}Ua_{\nu_{1}}^{\eta_{1}}{}^{\dagger}\ldots a_{\nu_{m}}^{\eta_{m}}{}^{\dagger}|0\rangle$ denote the monochromatic scattering elements from $m$ incoming photons with frequencies $\nu_{1},\ldots,\nu_{m}$, and directions $\eta_{1},\ldots,\eta_{m}$, to $n$ outgoing ones with frequencies $\omega_{1},\ldots,\omega_{n}$, and directions $\rho_{1},\ldots,\rho_{n}$, respectively.

In the following, we specialize the discussion to the single- and two-photon matrix elements of a scatterer with a single ground state. In this case, the conservation of energy allows us to decompose the monochromatic scattering components as in Eq. (The need of wave packets.–) of the main text:

[TABLE]

Here, $\chi_{\omega}^{\rho\eta}$ contain the reflection $r_{\omega}=\chi_{\omega}^{-\rho,\rho}$ and transmission $t_{\omega}=\chi_{\omega}^{\rho,\rho}$ coefficients, and the nonlinear contribution $\bar{\mathcal{T}}_{\omega_{1}\omega_{2}\nu_{1}\nu_{2}}^{\rho_{1}\rho_{2}\eta_{1}\eta_{2}}$ describes photon-photon interactions mediated by the scatterer.

Using Eqs. (88)-(89) in Eq. (87), simplifies the integral relations by reducing their dimensionality. In particular, the single-photon scattering matrix satisfies,

[TABLE]

On the other hand, the two-photon ‘wave packet’ scattering matrix can be also decomposed into linear and nonlinear contributions,

[TABLE]

where the non-linear part reads,

[TABLE]

Here, the kernel $\mathcal{W}$ is the product of three Gaussians,

[TABLE]

and we conveniently defined the new variables $\bar{\omega}=(\omega_{1}+\omega_{2})/2$, and $\bar{\nu}=(\nu_{1}+\nu_{2})/2$, $\Delta_{\nu}=(\nu_{2}-\nu_{1})/2$, and $\Delta_{\omega}=(\omega_{2}-\omega_{1})/2$.

Interestingly, both Eqs. (90) and (92) can be analytically inverted as they involve a product of Gaussian kernels for independent integration variables. Therefore, following Refs. Ulmer and Kaissl (2003); Ulmer (2010); Fang et al. (1994) of the main text, we obtain

[TABLE]

where the inverse Gaussian kernel $\mathcal{K}_{\tilde{\sigma}}(\omega)$ is given by

[TABLE]

Here, $H_{q}(x)=(-1)^{q}e^{x^{2}}\partial_{x}^{q}(e^{-x^{2}})$ are Hermite polynomials and the series is convergent because $S_{\omega_{1}\omega_{1}}^{\rho_{1}\eta_{1}}$ and $T_{\omega_{1}\omega_{2}\nu_{1}\nu_{2}}^{\rho_{1}\rho_{2}\eta_{1}\eta_{2}}$ decay exponentially fast.

These last Eqs. (94)-(95) are direct analytical relations for the single- and two-photon scattering coefficients in terms of the ‘wave packet’ counterparts, measurable with our tomography protocol. Moreover, they are valid for any dispersive photonic channel, satisfying $\omega_{k}=\omega_{-k}$ and $v_{k}=-v_{-k}$. In the next subsection, we connect to the expressions stated in the main text, by transforming to momentum variables, and specializing to a linear dispersion relation $\omega_{k}=c|k|$.

III.2 B. Gaussian deconvolution in momentum modes for a non-dispersive channel

The monochromatic scattering elements in momentum and frequency basis can be related using Eq. (82) as

[TABLE]

where we used the identifications, $\omega_{j}=\omega_{p_{j}}$, $\nu_{j}=\omega_{k_{j}}$, $\rho_{j}={\rm sign}(p_{j})$, and $\eta_{j}={\rm sign}(k_{j})$. In addition, using the explicit formulas (88)-(89), we obtain the same formulas for the momentum monochromatic elements in the main text,

[TABLE]

but now with the coefficients given for dispersive media as,

[TABLE]

The wave packet profile used in Eq. (86) to perform the deconvolutions, can be recast in momentum as

[TABLE]

where the $\sigma=\tilde{\sigma}/c$ is the momentum width used in the main text.

Finally, when evaluating the above expressions for a linear dispersion relation $\omega_{k}=c|k|$ and forward scattering $k_{j},p_{j}>0$, as in the main text. In particular, the wave packet profile reduces simply to $\psi_{k}(k^{\prime})=G_{\sigma}\left(k^{\prime}-k\right)^{1/2}$ and the scattering coefficients read,

[TABLE]

where $\hat{k}=(k_{1}+k_{2})/2=(p_{1}+p_{2})/2$, $\Delta_{p}=(p_{2}-p_{1})/2$, and $\Delta_{k}=(k_{2}-k_{1})/2$.

III.3 C. One- and two-photon scattering matrix for a qubit weakly coupled to a 1D waveguide

In this subsection, we explicitly state the results shown in Refs. Shen and Fan (2007); Fan et al. (2010) of the main text, for the one- and two-photon scattering matrices of a two-level system weakly coupled to a 1D photonic waveguide. We use these analytical results in Figs. 3(b)-(d) of the main text to test the deconvolution formulas shown above. In particular, the monochromatic single-photon scattering matrix has the form of Eq. (88) with the reflection $r_{\omega}=\chi_{\omega}^{-\rho,\rho}$ and transmission $t_{\omega}=\chi_{\omega}^{\rho,\rho}$ coefficients given by

[TABLE]

Here $\omega_{0}$ is the transition frequency of the two-level system and $\gamma$ its decay rate. The two-photon scattering matrix is of the form of Eq. (89), with the two-photon interaction strength given by

[TABLE]

This is the nonlinearity we replace in Eq. (18) of the main text (with the proper change of variables as explained in the previous subsections) to evaluate the prediction for a scattering experiment with gaussian wave packets and a weakly coupled qubit. We then perform a deconvolution of the result using Eq. (Deconvolution formulas.–) of the main text, and check that the result agrees with the initial Eq. (107), which allows us to test the deconvolution formulas.

IV IV. Imperfections in the input state preparation

In this section we discuss errors in the reconstruction of the scattering matrix elements, which are caused by imperfections in the preparation of the coherent state inputs. The analysis includes deviations in the relative phases $s_{j}=\pm 1+\delta s_{j}^{(\pm)}$ [cf. Sec. IV.A], laser power $|\alpha|^{2}+\delta n$ [cf. Sec. IV.B], and global phases $\phi_{l}=\pi l/M+\delta\phi_{l}$ [cf. Sec. IV.C].

IV.1 A. Deviations in relative phases

In the case the scatterer conserves the photon number, we can apply Protocol 2 of the main text, and therefore the preparation of the various input states requires only the control of relative phases $s_{j}=\pm 1$, for $j=1,\ldots,M$. If we go back to Eq. (48) in the derivation of the protocol and assume a small deviation, $s_{j}=\pm 1+\delta s_{j}^{(\pm)}$, in the ideal relative phases, we obtain the relation,

[TABLE]

where the ‘sign’ error $\varepsilon_{\rm sign}$ is given to first order in the small deviations $|\delta s_{j}^{(\pm)}|\ll 1$ as

[TABLE]

with $\delta\bar{s}_{j}$ he average deviation $\delta\bar{s}_{j}=(\delta s_{j}^{(+)}+\delta s_{j}^{(-)})/2$.

The rest of the derivation for the reconstruction formula is the same as in Sec. I.B. In particular, when replacing Eq. (35) into the imperfect relation (108), we obtain the usual estimate of the scattering matrix in Eqs. (51) and (66), but now with the extra ’sign’ error given above,

[TABLE]

Notice that this small error in the relative phases is linear in the deviations, and thus can be controlled provided $\varepsilon_{\rm sign}\sim|\delta s_{j}^{(\pm)}|\ll 1$. Since $\varepsilon_{\rm sign}$ does not depend on the laser power $|\alpha|^{2}$, the formula for the $Z$-order estimate is straightforwardly derived by reapeating the procedure in Sec. II., obtaining simply

[TABLE]

with $E^{(Z)}(|\alpha|^{2})$ and $\varepsilon_{m}^{(Z)}$ as given in Eqs. (72) and (70).

IV.2 B. Imperfections in laser power input

In the case there are imperfections in setting the mean photon number of the input states, $|\alpha|^{2}+\delta n$, the procedure in Sec. II for combining $Z$ estimates at well defined laser powers, $x_{q}=b^{q-1}|\alpha|^{2}$, will be affected by an error as well. For a first order estimate, this error appears as

[TABLE]

Using the above expression in the procedure of Sec. II to combine $Z$ estimates, we generalize Eq. (111),

[TABLE]

with the extra error $\varepsilon_{\alpha}^{(Z)}$ due to imperfections in the laser power, given by

[TABLE]

Notice that for $\delta n/|\alpha|^{2}\ll 1$ this error is also linear in the deviations,

[TABLE]

and it is small $|\varepsilon^{(Z)}_{\alpha}|\sim\delta n/|\alpha|^{2}\ll 1$ since high order terms in the sum (with $v\gg 1$) are strongly suppressed by the coefficients $|f_{vm}^{(Z)}|\sim 1/v!$ in Eq. (71).

IV.3 C. Deviations in global phases

In this last section we discuss errors due to deviations in the global phase of the input states, $\phi_{l}=\pi l/M+\delta\phi_{l}$, for $l=1,\ldots,2M$. This is needed in the case that we do not know a priori if the scatterer conserves the photon number, and we have to apply the general Protocol 1 of the main text. If we go back to Eq. (35) in the derivation of the protocol (with $\lambda\!=\!1$), the deviations $\delta\phi_{l}$ will add an extra error to this expression, obtaining

[TABLE]

In addition to the usual approximation error $\epsilon_{m}^{(1)}$, the error due to imperfect global phases $\delta\phi_{l}$ reads,

[TABLE]

Including the error due to imperfect relative phases, we replace Eq. (117) into Eq. (108), and following the usual derivation, we obtain a more general expression for the first order estimate of $S$,

[TABLE]

Here, $\varepsilon_{m}^{(1)}$ and $\varepsilon_{\rm sign}$ are the approximation and sign errors already discussed, and the first order error due to imperfections in the global phases reads

[TABLE]

Additionally, if we now include the error due to imperfections in the laser power as in Eq. (114), with $\delta n/|\alpha|^{2}\ll 1$, and combine $Z$ estimates as in Sec. II, we derive the most general $Z$-order estimate for $S$,

[TABLE]

with the $Z$-order error due to imperfect global phases given by

[TABLE]

Here, only the few first terms in the expansion contribute to the error as $|T_{ud}|\sim 1/u!$ in $t_{um}$ strongly suppresses higher order terms with $u\gg 1$. Finally, notice that if we use moderate laser powers $|\alpha|\gtrsim 1$ and assume small deviations in the global phases $|\delta\phi_{l}|\ll 1$, the resulting error is small and linear in the deviations, scaling approximately as $|\varepsilon_{m}^{\phi,(Z)}|\sim|\delta\phi_{l}|/|\alpha|^{m-1}\sim|\delta\phi_{l}|\ll 1$.