††thanks: Source code for the numerical simulations in this paper is available
here

# Few-particle scattering from localized quantum systems in spatially
structured bosonic baths

Rahul Trivedi  Max Planck Institut für Quantenoptik, Garching bei München,
85748, Germany<EMAIL_ADDRESS>Kevin Fischer  Stanford University,
Stanford, CA, 94305, USA.  Shanhui Fan  Stanford University, Stanford, CA,
94305, USA.  Jelena Vuckovic  Stanford University, Stanford, CA, 94305, USA.

###### Abstract

Understanding the dynamics of localized quantum systems embedded in engineered
bosonic environments is a central problem in quantum optics and open quantum
system theory. We present a formalism for studying few-particle scattering in
a localized quantum system interacting with an bosonic bath described by an
inhomogeneous wave-equation. In particular, we provide exact relationships
between the quantum scattering matrix of this interacting system and frequency
domain solutions of the inhomogeneous wave-equation thus providing access to
the spatial distribution of the scattered few-particle wave-packet. The
formalism developed in this paper paves the way to computationally
understanding the impact of structured media on the scattering properties of
localized quantum systems embedded in them without simplifying assumptions on
the physics of the structured media.

##  I Introduction

Non-classical light sources lodahl2015interfacing ; ding2016demand ; michler ;
senellart2017high ; zhang2018strongly ; kok2007linear ; o2009photonic ;
roy2017colloquium ; reiserer2015cavity ; duan2010colloquium ;
sangouard2011quantum ; nemoto2014photonic  are an essential resource in all
optical quantum information processing and communication systems.
Traditionally, such light sources have been designed by interfacing localized
quantum systems (quantum dots, color centers etc.) with a bosonic bath (an
electromagnetic structure such as an optical fiber, waveguide or optical
cavity). The localized quantum system can modify the statistics and the
entanglement structure of an incoming state propagating in the bosonic bath
englund2005controlling ; daveau2017efficient ; pelton2002efficient ;
he2013demand ; fischer2017signatures ; hanschke2018quantum ;
metcalfe2010resolved ; miao2019electrically ; lukin2020spectrally . This has
opened up the possibility of generating high-fidelity single photon sources
he2013demand ; wang2019towards , photon pair sources li2004all ;
harder2013optimized  as well as high-dimensional entangled photon states
pichler2016photonic ; xu2018generate  by engineering both the bath as well as
the Hamiltonian of the localized quantum systems.

From a theoretical standpoint, a quantitative theory of the dynamics and
scattering properties of this coupled system is required for designing any of
the above-mentioned applications. If only the dynamics of the localized system
is of interest, the usual approach is to trace out the Hilbert space of the
bath to obtain a master equation for the density operator of the localized
system carmichael2009open ; gardiner2004quantum . While traditionally the
master equation is formulated for a one-dimensional homogeneous bath
gardiner2004quantum , it is easily extended to cases where the bath is
described by an inhomogeneous 3D wave-equation asenjo2017atom ;
dung2002resonant . However, a significant drawback of the master equation
formalism is that it is not amenable to computing the state of the bath,
although recent works have attempted to resolve this problem for 1D
homogeneous baths kiilerich2019input .

An alternative set of computational and analytical tools which allow us to
calculate the bath’s quantum state is provided by scattering theory
taylor2006scattering . Here, the localized system is treated as a scatterer
for the quantum states that otherwise propagate in the bath. Significant
progress has been made on scattering theory for systems where a modal
decomposition of the bosonic bath is easily accessible. In particular, there
has been great success in developing techniques for analyzing the scattering
matrices of the localized quantum system when the electromagnetic structure is
a photonic waveguide fan2010input ; xu2013analytic ; xu2015input ;
trivedi2018few ; trivedi2021optimal ; fischer2017scattering ;
shi2015multiphoton ; caneva2015quantum , a cavity haroche1993cavity ;
trivedi2019photon ; shi2013two  or a network of these components. There has
also been progress towards understanding the interaction of localized quantum
systems with 3D homogeneous baths liu2016quantum  and periodically structured
baths orenstein2020two . Computing the scattering of few-photon states in
systems where the bosonic bath is describable as arbitrary spatially varying
permittivity distributions remains an open problem. Such a framework, if it
exists, can enable both analysis and design of more complex quantum photonic
devices similar to how computational electromagnetic methods capable of
simulating arbitrary permittivity distributions revolutionized the design of
classical photonic devices miller1988selective ; molesky2018inverse .

In this paper, we present a framework for exactly calculating the scattering
matrix of a localized quantum system interacting with a bosonic bath whose
dynamics are described by an inhomogeneous wave equation. There are two key
ingredients that go towards making this framework possible —

  1. 1.

A second quantized description of the inhomogeneous bosonic bath in terms of
both Langevin noise operators and plane-wave annihilation operators (Section
II). Here, we build upon the formalism presented in Refs. dung1998three ;
gruner1996green  for quantizing lossy frequency-dependent and inhomogeneous
electromagnetic environments.

  2. 2.

Calculation of few-particle scattering matrices that builds on this second
quantized description (Section III). We express the spatial and spectral
dependence of the few-particle scattering matrices in terms of two quantities
— _first_ , the frequency-domain solutions of the inhomogeneous wave-equation
describing the bath dynamics and _second_ , the time-ordered expectations of
the localized system operators (Section III). While our general approach is
applicable to arbitrary localized systems, we provide exact expressions for
the single and two-particle scattering matrices for the paradigmatic example
of a two-level localized system.

##  II Quantization of an inhomogeneous bosonic bath

Figure 1: Schematic. A localized two-level system interacting with an
inhomogeneous bosonic bath described by a permittivity distribution
$\varepsilon(\mathbf{x})$ which is asymptotically homogeneous
($\varepsilon(\mathbf{x})\to\varepsilon_{0}$ as $|\mathbf{x}|\to\infty$). An
incident few-particle state $\ket{\psi_{\text{in}}}$ propagates through the
inhomogeneous bath and scatters off the localized system to produce an
outgoing few-particle state $\ket{\psi_{\text{out}}}$.

Figure 1 schematically depicts the systems under consideration comprising of a
localized quantum system interacting with an inhomogeneous bosonic bath. We
assume that the bosonic bath is classically described by the inhomogenous
wave-equation

$\nabla^{2}\psi(\mathbf{x},t)=\varepsilon(\mathbf{x})\frac{\partial^{2}}{\partial
t^{2}}\psi(\mathbf{x},t),$ (1)

where $\psi(\mathbf{x},t)$ is the time-dependent field variable and
$\varepsilon(\mathbf{x})$ is the permittivity distribution describing the
inhomogeneity in the bath. While we restrict ourselves to a simple scalar
wave-equation model in this paper, an extension to a full vector-wave equation
model is also possible. Furthermore, we will assume that the bosonic bath in
isolation is lossless and consequently the permittivity distribution
$\varepsilon(\mathbf{x})$ is real and positive. Furthermore, we assume that
the inhomogeneous media is embedded in a background homogeneous media with
permittivity $\varepsilon_{0}$ i.e.
$\varepsilon(\mathbf{x})\to\varepsilon_{0}$ as $|\mathbf{x}|\to\infty$.

###  II.1 First quantization

A quantum theory for such a bosonic bath can easily be developed by applying
Dirac’s quantization weinberg1995quantum  procedure on its classical
Lagrangian $\mathcal{L}_{\text{bath}}[\psi]$

$\mathcal{L}_{\text{bath}}[\psi]=\frac{1}{2}\int_{\mathbb{R}^{3}}\bigg{(}\varepsilon(\mathbf{x})\bigg{|}\frac{\partial\psi(\mathbf{x},t)}{\partial
t}\bigg{|}^{2}-|\nabla\psi(\mathbf{x},t)|^{2}\bigg{)}d^{3}\mathbf{x}.$ (2)

From this lagrangian, the canonically conjugate field to $\psi(\mathbf{x},t)$
is easily identified to be

$\displaystyle\pi(\mathbf{x},t)=\varepsilon(\mathbf{x})\frac{\partial\psi(\mathbf{x},t)}{\partial
t}.$ (3)

Dirac’s quantization then promotes both of these fields to operators,
$\psi(\mathbf{x})$ and $\pi(\mathbf{x})$, with commutator
$[\psi(\mathbf{x}),\pi(\mathbf{x}^{\prime})]=i\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime})$.
Furthermore, the dynamics of the quantum state of the bath is governed by the
Hamiltonian

$\displaystyle
H_{\text{bath}}=\int_{\mathbb{R}^{3}}\bigg{(}\frac{\pi^{2}(\mathbf{x})}{2\varepsilon(\mathbf{x})}+\frac{1}{2}|\nabla\psi(\mathbf{x})|^{2}\bigg{)}d^{3}\mathbf{x}$
(4)

It can be readily verified that the Heisenberg equations of motion for the
operators $\psi(\mathbf{x})$ and $\pi(\mathbf{x})$ reproduce the equations of
motion for the classical fields as given by Eqs. 1 and 3.

###  II.2 Second quantization

While the first quantized description provides a complete quantum theory for
the bosonic bath, it does not allow for an easy and physically meaningful
description of incident and scattered quantum states. It is, therefore, of
interest to perform second quantization and introduce particle annihilation
and creation operators. We accomplish a second quantization of an
inhomogeneous bosonic bath in two steps. _First_ , we consider a more general
case of lossy inhomogeneous media and introduce a second quantization based on
Langevin noise formalism. _Second_ , we treating lossless inhomogeneous media
as a limiting case of lossy inhomogeneous media, we construct a second
quantization in which an annihilation operator is associated with each
possible plane-wave incident on the inhomogeneous media.

####  II.2.1 Langevin noise formalism

In this subsection, we build upon the formalism introduced in Refs.
dung1998three ; gruner1996green  which quantize the wave-equation in lossy
inhomogeneous Kramer-Kronig dielectrics using a Langevin noise-operator
approach. However, since we are fundamentally interested in lossless
inhomogeneous baths, to make use of this formalism we fictitiously introduce
loss in the bosonic bath, and consider a frequency-dependent permittivity
distribution $\varepsilon_{\alpha}(\mathbf{x},\omega)$ given by

$\displaystyle\varepsilon_{\alpha}(\mathbf{x},\omega)=\varepsilon(\mathbf{x})+i\alpha(\omega),$
(5)

where $\alpha(\omega)>0$ is the loss in the permittivity distribution. For
simplicity, we assume the loss to be position independent — the permittivity
distribution $\varepsilon_{\alpha}(\mathbf{x},\omega)$ is thus that of a lossy
inhomogeneous permittivity distribution in a homogeneous background with
permittivity $\varepsilon_{0}+i\alpha(\omega)$. While we could work with a
position dependent loss as well, since the objective is to study the limiting
case of a lossless permittivity distribution, the position dependence of the
loss does not matter.

With the introduction of this loss, following Ref. dung1998three  we introduce
a bosonic noise operator $\phi_{\alpha}(\mathbf{x},\omega)$ for every position
$\mathbf{x}$ and for every frequency $\omega>0$ which satisfy
$[\phi_{\alpha}(\mathbf{x},\omega),\phi^{\dagger}_{\alpha}(\mathbf{x}^{\prime},\omega^{\prime})]=\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime})\delta(\omega-\omega^{\prime})$.
Furthermore, the bath’s Hamiltonian is expressible as111We use the notation
that $\mathbb{R}^{+}:=[0,\infty)$,

$H_{\text{bath}}=\int_{\begin{subarray}{c}\mathbf{x}\in\mathbb{R}^{3}\\\
\omega\in\mathbb{R}^{+}\end{subarray}}\omega\phi^{\dagger}_{\alpha}(\mathbf{x},\omega)\phi_{\alpha}(\mathbf{x},\omega)\
d^{3}\mathbf{x}\ d\omega.$ (6)

We define the single-frequency field operator in the presence of loss,
$\psi_{\alpha}(\mathbf{x},\omega)$, in terms of the noise operators
$\phi_{\alpha}(\mathbf{x},\omega)$ by

$\displaystyle\psi_{\alpha}(\mathbf{x},\omega)=\bigg{(}\frac{\omega^{2}\alpha(\omega)}{\pi}\bigg{)}^{1/2}\times$
$\displaystyle\qquad\int_{\mathbf{x}^{\prime}\in\mathbb{R}^{3}}\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)\phi_{\alpha}(\mathbf{x}^{\prime},\omega)d^{3}\mathbf{x}^{\prime}d\omega,$
(7)

where $\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)$ is the
Green’s function corresponding to the lossy inhomogeneous bath, which is
obtained by solving the partial differential equation

$\big{[}\nabla_{\mathbf{x}}^{2}+\omega^{2}\varepsilon_{\alpha}(\mathbf{x},\omega)\big{]}\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)=\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime}).$
(8)

under the boundary condition that
$\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)$ is an outgoing wave
as $|\mathbf{x}|\to\infty$. We point out that since
$\text{Im}[\varepsilon_{\alpha}(\mathbf{x},\omega)]=\alpha(\omega)>0$, this
boundary condition also implies that
$\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)$ decays
exponentially with $|\mathbf{x}|$ as $|\mathbf{x}|\to\infty$. The field
operators in the presence of loss, $\psi_{\alpha}(\mathbf{x})$ and
$\pi_{\alpha}(\mathbf{x})$, are related to the single-frequency field operator
via

$\displaystyle\psi_{\alpha}(\mathbf{x})=i\int_{\omega\in\mathbb{R}^{+}}\big{(}\psi_{\alpha}(\mathbf{x},\omega)-\psi_{\alpha}^{\dagger}(\mathbf{x},\omega)\big{)}d\omega,$
(9a)
$\displaystyle\pi_{\alpha}(\mathbf{x})=\varepsilon(\mathbf{x})\int_{\omega\in\mathbb{R}^{+}}\omega\big{(}\psi_{\alpha}(\mathbf{x},\omega)-\psi_{\alpha}^{\dagger}(\mathbf{x},\omega)\big{)}d\omega.$
(9b)

It is shown in appendix A that the single-frequency field operator
$\psi_{\alpha}(\mathbf{x},\omega)$ has the commutator
$[\psi_{\alpha}(\mathbf{x},\omega),\psi_{\alpha}^{\dagger}(\mathbf{x}^{\prime},\omega^{\prime})]=-{\delta(\omega-\omega^{\prime})}\text{Im}\big{[}\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)\big{]}/{2\pi}$.
Furthermore, this implies that the operators $\psi_{\alpha}(\mathbf{x})$ and
$\pi_{\alpha}(\mathbf{x})$ satisfy the canonical commutation relations
$[\psi_{\alpha}(\mathbf{x}),\pi_{\alpha}(\mathbf{y})]=i\delta^{3}(\mathbf{x}-\mathbf{y})$.
Finally, we point out the field operators for a lossless bosonic bath,
$\psi(\mathbf{x},\omega),\psi(\mathbf{x}),\pi(\mathbf{x})$, can be constructed
from the operators
$\psi_{\alpha}(\mathbf{x},\omega),\psi_{\alpha}(\mathbf{x}),\pi_{\alpha}(\mathbf{x})$
by taking $\alpha(\omega)\to 0$. Since $\lim_{\alpha(\omega)\to
0}\text{Im}\big{[}\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)]$
exists and is well defined, this immediately implies the following commutator
for the single-frequency field operator

$\displaystyle[\psi(\mathbf{x},\omega),\psi^{\dagger}(\mathbf{x}^{\prime},\omega^{\prime})]=-\frac{\delta(\omega-\omega^{\prime})}{2\pi}\text{Im}\big{[}\text{G}(\mathbf{x},\mathbf{x}^{\prime};\omega)\big{]},$
(10)

and the canonical commutation relation,
$[\psi(\mathbf{x}),\pi(\mathbf{x}^{\prime})]=i\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime})$,
holds between the full field operators $\psi(\mathbf{x})$ and
$\pi(\mathbf{x})$.

####  II.2.2 Second quantization with plane-wave operators

Since the physics of an inhomogeneous bath can fundamentally be captured by
the fields scattered from the inhomogeneity when excited by incident plane-
waves, we seek its second quantized description in terms of incident plane-
wave annihilation operators. This can be derived from the noise-operator
description by taking the limit of $\alpha(\omega)\to 0$. In order to take
this limit, we express the Green’s function
$\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)$ as

$\displaystyle\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)=$
$\displaystyle\quad\text{G}_{\alpha}^{0}(\mathbf{x},\mathbf{x}^{\prime};\omega)+\int_{\mathbb{R}^{3}}\text{S}_{\alpha}(\mathbf{x},\mathbf{y};\omega)\text{G}_{\alpha}^{0}(\mathbf{y},\mathbf{x}^{\prime};\omega)d^{3}\mathbf{y},$
(11)

where $\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)$ is the
Green’s function corresponding to the (lossy) background i.e. it is the
solution of the partial differential equation

$\displaystyle\big{[}\nabla_{\mathbf{x}}^{2}+\omega^{2}(\varepsilon_{0}+i\alpha(\omega))\big{]}\text{G}_{\alpha}^{0}(\mathbf{x},\mathbf{x}^{\prime};\omega)=\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime}),$
(12)

with the boundary condition that
$\text{G}_{\alpha}^{0}(\mathbf{x},\mathbf{x}^{\prime};\omega)$ is an outgoing
spherical wave when $|\mathbf{x}|\to\infty$ and
$\text{S}_{\alpha}(\mathbf{x},\mathbf{y};\omega)$ is the position-domain
representation of the classical scattering operator of the (lossy)
inhomogeneous bath. The scattering operator is the operator that maps an
incident field to the corresponding fields scattered from the inhomogeneous
permittivity distribution. Substituting this expansion into Eq. II.2.1, we
obtain

$\psi_{\alpha}(\mathbf{x},\omega)=\psi^{0}_{\alpha}(\mathbf{x},\omega)+\int_{\mathbb{R}^{3}}\text{S}_{\alpha}(\mathbf{x},\mathbf{y};\omega)\psi^{0}_{\alpha}(\mathbf{y},\omega)d^{3}\mathbf{y},$
(13a) where
$\psi^{0}_{\alpha}(\mathbf{x},\omega)=\bigg{(}\frac{\omega^{2}\alpha(\omega)}{\pi}\bigg{)}^{1/2}\int_{\mathbb{R}^{3}}\text{G}_{\alpha}^{0}(\mathbf{x},\mathbf{x}^{\prime};\omega)\phi_{\alpha}(\mathbf{x}^{\prime};\omega)d^{3}\mathbf{x}^{\prime}.$
(13b)

Physically, Eq. 13a expresses the fields inside the inhomogeneous bath as the
sum of an incident field and a scattered field. Furthermore, since the
operator $\psi^{0}(\mathbf{x};\omega)$ describes the excitations incident on
the scattering region, it is expressible as a sum of plane-wave modes
propagating in the background media. Formally, this is accomplished by
introducing the plane-wave decomposition of
$\text{G}_{\alpha}^{0}(\mathbf{x},\mathbf{x}^{\prime};\omega)$,

$\text{G}_{\alpha}^{0}(\mathbf{x},\mathbf{x}^{\prime};\omega)=\frac{1}{(2\pi)^{3}}\int_{\begin{subarray}{c}\Omega\in\mathbb{S}^{2}\\\
k\in\mathbb{R}^{+}\end{subarray}}\frac{e^{ik\Omega\cdot(\mathbf{x}-\mathbf{x}^{\prime})}\
k^{2}d^{2}\Omega\
dk}{\omega^{2}\varepsilon_{0}-k^{2}+i\omega^{2}\alpha(\omega)},$ (14)

into Eq. 13b. This allows us to express $\psi^{0}_{\alpha}(\mathbf{x};\omega)$
as

$\psi^{0}_{\alpha}(\mathbf{x},\omega)=N_{0}(\omega)\int_{\Omega\in\mathbb{S}^{2}}a_{\alpha}(\mathbf{x},\Omega;\omega)e^{ik_{0}(\omega)\Omega\cdot\mathbf{x}}\
d^{2}\Omega,$ (15)

where $N_{0}(\omega)=({\omega\sqrt{\varepsilon_{0}}}/{16\pi^{3}})^{1/2}$ and
$k_{0}(\omega)=\omega\sqrt{\varepsilon_{0}}$ and

$\displaystyle
a_{\alpha}(\mathbf{x},\Omega,\omega)=\bigg{(}\frac{16\pi^{2}\omega\alpha(\omega)}{\sqrt{\varepsilon_{0}}}\bigg{)}^{1/2}e^{-ik_{0}(\omega)\Omega\cdot\mathbf{x}}\times$
$\displaystyle\qquad\int_{\begin{subarray}{c}\mathbf{x}^{\prime}\in\mathbb{R}^{3}\\\
{k\in\mathbb{R}^{+}}\end{subarray}}\frac{e^{ik\Omega\cdot(\mathbf{x}-\mathbf{x}^{\prime})}}{\omega^{2}\varepsilon_{0}-k^{2}+i\omega^{2}\alpha(\omega)}\phi_{\alpha}(\mathbf{x}^{\prime},\omega)\
d^{3}\mathbf{x}^{\prime}\ dk.$ (16)

It is shown in appendix B that taking the limit of $\alpha(\omega)\to 0$ makes
$a_{\alpha}(\mathbf{x},\Omega,\omega)$ independent of $\mathbf{x}$, enabling
us to define $a(\Omega,\omega)=\lim_{\alpha(\omega)\to
0}a_{\alpha}(\mathbf{x},\Omega,\omega)$. Furthermore, $a(\Omega,\omega)$
satisfies the commutator

$\displaystyle[a(\Omega,\omega),a^{\dagger}(\Omega^{\prime},\omega^{\prime})]=\delta^{2}(\Omega-\Omega^{\prime})\delta(\omega-\omega^{\prime}),$
(17)

and Eq. 15 in the limit of $\alpha(\omega)\to 0$ reduces to

$\psi^{0}(\mathbf{x},\omega)=N_{0}(\omega)\int_{\Omega\in\mathbb{S}^{2}}a(\Omega,\omega)e^{ik_{0}(\omega)\Omega\cdot\mathbf{x}}d^{2}\Omega.$
(18)

The operators $a(\Omega,\omega)$ can thus be interpreted as the annihilation
operator for a plane-wave mode incident on the inhomogeneous bath from the
surrounding background and the full field operator $\psi(\mathbf{x},\omega)$
is then given by the sum of the incident and scattered field operators

$\psi(\mathbf{x},\omega)=N_{0}(\omega)\int_{\Omega\in\mathbb{S}^{2}}\mathcal{E}(\mathbf{x},\Omega;\omega)a(\Omega,\omega)d^{2}\Omega,$
(19)

where

$\displaystyle\mathcal{E}(\mathbf{x},\Omega;\omega)=e^{ik_{0}(\omega)\Omega\cdot\mathbf{x}}+\mathcal{E}_{s}(\mathbf{x},\Omega;\omega)\text{
with }$ (20a)
$\displaystyle\mathcal{E}_{s}(\mathbf{x},\Omega;\omega)=\int_{\mathbb{R}^{3}}\text{S}(\mathbf{x},\mathbf{y};\omega)e^{ik_{0}(\omega)\Omega\cdot\mathbf{y}}d^{3}\mathbf{y},$
(20b)

with $\text{S}(\mathbf{x},\mathbf{y};\omega)=\lim_{\alpha(\omega)\to
0}\text{S}_{\alpha}(\mathbf{x},\mathbf{y};\omega)$. Finally, noting from Eqs.
6 and II.2.2 that $[a(\Omega,\omega),H_{\text{bath}}]=\omega
a(\Omega,\omega)$, and consequently

$H_{\text{bath}}=\int_{\begin{subarray}{c}\Omega\in\mathbb{S}^{2}\\\
\omega\in\mathbb{R}^{+}\end{subarray}}\omega
a^{\dagger}(\Omega,\omega)a(\Omega,\omega)\ d^{2}\Omega\ d\omega.$ (21)

##  III Scattering from a localized two-level system

In this section, we consider the problem of embedding a localized quantum
system, which for simplicity is assumed to be a two-level system (TLS) with
ground state $\ket{g}$ and excited state $\ket{e}$, in the inhomogeneous bath.
We make the point-dipole approximation and assume that the localized quantum
system interacts with the bath’s field at only $\mathbf{x}=\mathbf{x}_{d}$.
Furthermore, we assume a rotating-wave approximation for the interaction
Hamiltonian for this two-level system to obtain

$\displaystyle
H_{\text{int}}=iV_{0}\int_{\omega\in\mathbb{R}^{+}}\psi(\mathbf{x}_{d},\omega)\sigma^{\dagger}d\omega+\text{h.c.},$
(22)

where $\sigma^{\dagger}=\ket{e}\bra{g}$, $\sigma=\ket{g}\bra{e}$ and $V_{0}$
is the coupling constant between the two-level system and the inhomogeneous
bath. This Hamiltonian can be expressed in terms of the operators
$a(\Omega,\omega)$ as

$\displaystyle H_{\text{int}}=$ $\displaystyle
iV_{0}\int_{\begin{subarray}{c}\Omega\in\mathbb{S}^{2}\\\
\omega\in\mathbb{R}^{+}\end{subarray}}N_{0}(\omega)\mathcal{E}(\mathbf{x}_{d},\Omega,\omega)a(\Omega,\omega)\sigma^{\dagger}\
d^{2}\Omega\ d\omega+\text{h.c.}$ (23)

It can be pointed out that we do not make the Markovian approximation
carmichael2009open ; breuer2002theory  in this interaction Hamiltonian — we
intend our analysis to account for any strong frequency dependent interactions
that can be induced by a spatial structuring of the bath. The total
Hamiltonian $H$ for the system can be expressed as a sum of this interaction
Hamiltonian, the Hamiltonian of the bath and the Hamiltonian of the two-level
system: $H=H_{0}+H_{\text{int}}$, where
$H_{0}=H_{\text{bath}}+\omega_{0}\sigma^{\dagger}\sigma$.

We now consider the problem of exciting the two-level system with $N$ incident
particles, and compute the spectral and spatial distribution of the scattered
state. The input state that we consider is that of $N$ plane waves at
frequencies $\vec{\nu}=(\nu_{1},\nu_{2}\dots\nu_{N})$222Throughout this paper,
we will adopt that notation that the symbol $\vec{f}$ represents the vector
$(f_{1},f_{2}\dots f_{N})$, where the length of the vector will be clear from
the context. and propagating in directions
$\vec{\Omega}=(\Omega_{1},\Omega_{2}\dots\Omega_{N})$

$\displaystyle\ket{\psi_{\text{in}}^{N}(\vec{\Omega},\vec{\nu})}=\prod_{i=1}^{N}a^{\dagger}(\Omega_{i},\nu_{i})\ket{\text{vac},g},$
(24)

and we measure the overlap of the scattered state with
$\ket{\psi_{\text{out}}(\vec{\mathbf{x}},\vec{\omega})}$, where

$\displaystyle\ket{\psi_{\text{out}}^{N}(\vec{\mathbf{x}},\vec{\omega})}=\prod_{i=1}^{N}\psi^{\dagger}(\mathbf{x}_{i},\omega_{i})\ket{\text{vac},g}.$
(25)

This corresponds to monitoring the scattered $N-$particle quantum field at
positions
$\vec{\mathbf{x}}=(\mathbf{x}_{1},\mathbf{x}_{2}\dots\mathbf{x}_{N})$ and at
frequencies $\vec{\omega}=(\omega_{1},\omega_{2}\dots\omega_{N})$. The
$N-$particle scattering matrix element corresponding to this measurement is
thus given by

$\displaystyle
S^{N}\big{(}\vec{\mathbf{x}},\vec{\omega};\vec{\Omega},\vec{\nu}\big{)}=\bra{\psi_{\text{out}}(\vec{\mathbf{x}},\vec{\omega})}\hat{\text{S}}\ket{\psi_{\text{in}}(\vec{\Omega},\vec{\nu})}$
(26)

where $\hat{\text{S}}$ is the scattering matrix defined via

$\displaystyle\hat{\text{S}}=\lim_{\begin{subarray}{c}t_{f}\to\infty\\\
t_{i}\to-\infty\end{subarray}}e^{iH_{0}t_{f}}e^{-iH(t_{f}-t_{i})}e^{-iH_{0}t_{i}}.$
(27)

Using this definition of $\hat{\text{S}}$, the scattering matrix-element can
be expressed in terms of the Heisenberg picture field operators
$\psi(\mathbf{x},\omega;t)=e^{iHt}\psi(\mathbf{x},\omega)e^{-iHt}$ and
$a(\Omega,\omega;t)=e^{iHt}a(\Omega,\omega)e^{-iHt}$

$\displaystyle
S^{N}\big{(}\vec{\mathbf{x}},\vec{\omega};\vec{\Omega},\vec{\nu}\big{)}=\lim_{\begin{subarray}{c}t_{f}\to\infty\\\
t_{i}\to-\infty\end{subarray}}e^{i\sum_{n=1}^{N}(\omega_{n}t_{f}-\nu_{n}t_{i})}\times$
$\displaystyle\qquad\langle
g|\mathcal{T}\bigg{[}\prod_{n=1}^{N}\psi(\mathbf{x}_{n},\omega_{n};t_{f})\prod_{n=1}^{N}a^{\dagger}(\Omega_{n},\nu_{n};t_{i})\bigg{]}|g\rangle,$
(28)

where $\mathcal{T}[\cdot]$ is the chronological time-ordering operator and we
have used
$e^{iH_{0}t}\psi(\mathbf{x},\omega)e^{-iH_{0}t}=\psi(\mathbf{x},\omega)e^{-i\omega
t}$ and $e^{iH_{0}t}a(\Omega,\omega)e^{-iH_{0}t}=a(\Omega,\omega)e^{-i\omega
t}$. Using the commutators in Eqs. 10 and 17, the Heisenberg equations of
motion for $\psi(\mathbf{x},\omega;t)$ and $a(\Omega,\omega;t)$ can be setup
and integrated to obtain

$\displaystyle a(\Omega,\omega;t_{f})=e^{-i\omega
t_{f}}\bigg{[}a(\Omega,\omega;t_{i})e^{i\omega t_{i}}-$
$\displaystyle\qquad\qquad
V_{0}N_{0}(\omega)\mathcal{E}^{*}(\mathbf{x}_{d},\Omega;\omega)\int_{t_{i}}^{t_{f}}\sigma(\tau)e^{i\omega\tau}d\tau\bigg{]},$
(29a) $\displaystyle\psi(\mathbf{x},\omega;t_{f})=e^{-i\omega
t_{f}}\bigg{[}\psi(\mathbf{x},\omega;t_{i})e^{i\omega t_{i}}+$
$\displaystyle\qquad\qquad\frac{V_{0}}{\pi}\text{Im}\big{[}\text{G}(\mathbf{x},\mathbf{x}_{d};\omega)\big{]}\int_{t_{i}}^{t_{f}}\sigma(\tau)e^{i\omega\tau}d\tau\bigg{]},$
(29b)

where $\sigma(t)=e^{iHt}\sigma e^{-iHt}$. As is shown in appendix C, Eqs. 29
and III can be used to obtain

$\displaystyle
S^{N}(\vec{\mathbf{x}},\vec{\omega};\vec{\Omega},\vec{\nu})=\bigg{(}\prod_{n=1}^{N}N_{0}(\nu_{n})\bigg{)}\sum_{k=1}^{N}\bigg{[}\sum_{\mathcal{B}_{k},\mathcal{D}_{k}}\bigg{(}\sum_{\mathcal{P}_{k}}\prod_{n=1}^{k}$
$\displaystyle\quad\mathcal{E}(\mathbf{x}_{\mathcal{P}_{k}\mathcal{B}_{k}^{N}(n)},\Omega_{\mathcal{D}_{k}^{N}(n)},\nu_{\mathcal{D}_{k}(n)})\delta\big{(}\omega_{\mathcal{P}_{k}\mathcal{B}_{k}^{N}(n)}-\nu_{\mathcal{D}_{k}(n)}\big{)}\bigg{)}\times$
$\displaystyle\quad\bigg{(}\prod_{n=1}^{N-k}\frac{V_{0}^{2}}{\pi}\text{Im}\big{[}\text{G}(\mathbf{x}_{\bar{\mathcal{B}}_{k}^{N}(n)},\mathbf{x}_{d};\omega_{\bar{\mathcal{B}}_{k}^{N}(n)})\big{]}\times$
$\displaystyle\quad\mathcal{E}(\mathbf{x}_{d},\Omega_{\bar{\mathcal{D}}_{k}^{N}(n)},\nu_{\bar{\mathcal{D}}^{N}_{k}(n)})\bigg{)}\mathcal{G}^{N-k}\big{(}\vec{\omega}_{\bar{\mathcal{B}}^{N}_{k}};\vec{\nu}_{\bar{\mathcal{D}}^{N}_{k}}\big{)}\bigg{]},$
(30)

where the notation used in Eq. III is

  1. 1.

$\mathcal{B}_{k}^{N}$ and $\mathcal{D}_{k}^{N}$ are two (independent)
unordered $k-$element subsets of $\\{1,2\dots N\\}$,

  2. 2.

$\bar{\mathcal{B}}_{k}^{N}$ and $\bar{\mathcal{D}}_{k}^{N}$ are complements of
$\mathcal{B}_{k}^{N}$ and $\mathcal{D}_{k}^{N}$ respectively,

  3. 3.

$\mathcal{P}_{k}$ is a $k-$element permutation,

  4. 4.

$\mathcal{P}_{k}\mathcal{B}_{k}^{N}$ ($\mathcal{P}_{k}\mathcal{D}_{k}^{N}$) is
a permutation of the $k-$element set $\mathcal{B}_{k}^{N}$
($\mathcal{D}_{k}^{N}$),

  5. 5.

For any set $\mathcal{A}$, $\mathcal{A}(i)$ is its $i^{\text{th}}$ element and

  6. 6.

$\vec{\omega}_{\bar{B}_{k}^{N}}=(\omega_{\bar{\mathcal{B}}_{k}^{N}(1)},\omega_{\bar{\mathcal{B}}_{k}^{N}(2)}\dots\omega_{\bar{\mathcal{B}}_{k}^{N}(N-k)})$
(similarly for $\bar{D}_{k}^{N}\vec{\nu}$).

Furthermore, $\mathcal{G}^{N}(\vec{\omega};\vec{\nu})$ is the $N-$excitation
Green’s function for the two-level system and is given by

$\displaystyle\mathcal{G}^{N}(\vec{\omega};\vec{\nu})=$
$\displaystyle\int_{\begin{subarray}{c}t_{1},t_{2}\dots t_{N}\in\mathbb{R}\\\
s_{1},s_{2}\dots
s_{N}\in\mathbb{R}\end{subarray}}\bra{\text{vac},g}\mathcal{T}\bigg{[}\prod_{i=1}^{N}\sigma(t_{i})\sigma^{\dagger}(s_{i})\bigg{]}\ket{\text{vac},g}\times$
$\displaystyle\qquad\prod_{j=1}^{N}e^{i(\omega_{j}t_{j}-\nu_{j}s_{j})}dt_{j}ds_{j}.$
(31)

We note that the spatial dependence of the scattering matrix-element in Eq.
III depends on $\mathcal{E}(\mathbf{x},\Omega,\omega)$, which are the fields
produced by the permittvity distribution on excitation by an incident plane-
wave, as well $\text{Im}\big{[}\text{G}(\mathbf{x},\mathbf{x}_{d};\omega)]$
that captures the fields radiated by the two-level system. However, while
$\text{Im}[\text{G}(\mathbf{x},\mathbf{x}_{d};\omega)]$ has contributions from
both incoming and outgoing waves, when considering the action of the
scattering-matrix on a normalizable $N-$particle wave-packet, only the
outgoing waves have a non-zero contribution. Mathematically, this follows from
the fact that $\text{G}^{*}(\mathbf{x},\mathbf{x}^{\prime};\omega)$ is
analytic in the lower-half of the complex plane and thus

$\displaystyle\lim_{t\to\infty}\int_{0}^{\infty}f(\omega)\text{G}^{*}(\mathbf{x},\mathbf{x}^{\prime};\omega)e^{-i\omega
t}d\omega=0,$ (32)

for any smooth and square-integrable function $f(\omega)$. Therefore, we can
only retain the outgoing part of
$\text{Im}\big{[}\text{G}(\mathbf{x},\mathbf{x}_{d};\omega)]$ in Eq. III and
obtain

$\displaystyle
S^{N}(\vec{\mathbf{x}},\vec{\omega};\vec{\Omega},\vec{\nu})=\bigg{(}\prod_{n=1}^{N}N_{0}(\nu_{n})\bigg{)}\sum_{k=1}^{N}\bigg{[}\sum_{\mathcal{B}_{k},\mathcal{D}_{k}}\bigg{(}\sum_{\mathcal{P}_{k}}\prod_{n=1}^{k}$
$\displaystyle\quad\mathcal{E}(\mathbf{x}_{\mathcal{P}_{k}\mathcal{B}_{k}^{N}(n)},\Omega_{\mathcal{D}_{k}^{N}(n)},\nu_{\mathcal{D}_{k}(n)})\delta\big{(}\omega_{\mathcal{P}_{k}\mathcal{B}_{k}^{N}(n)}-\nu_{\mathcal{D}_{k}(n)}\big{)}\bigg{)}\times$
$\displaystyle\quad\bigg{(}\prod_{n=1}^{N-k}\frac{V_{0}^{2}}{2\pi
i}\text{G}(\mathbf{x}_{\bar{\mathcal{B}}_{k}^{N}(n)},\mathbf{x}_{d};\omega_{\bar{\mathcal{B}}_{k}^{N}(n)})\times$
$\displaystyle\quad\mathcal{E}(\mathbf{x}_{d},\Omega_{\bar{\mathcal{D}}_{k}^{N}(n)},\nu_{\bar{\mathcal{D}}^{N}_{k}(n)})\bigg{)}\mathcal{G}^{N-k}\big{(}\vec{\omega}_{\bar{\mathcal{B}}^{N}_{k}};\vec{\nu}_{\bar{\mathcal{D}}^{N}_{k}}\big{)}\bigg{]},$
(33)

Equation III is one of the key results of this paper — to compute the quantum
scattering matrix, we need to compute $\mathcal{E}(\mathbf{x},\Omega;\omega)$
and $\text{G}(\mathbf{x},\mathbf{x}^{\prime};\omega)$, which can be obtained
by solving the (classical) inhomogeneous wave equation, and
$\mathcal{G}^{N}(\vec{\omega};\vec{\nu})$, the Green’s function of the quantum
emitter. Furthermore, the $N-$particle scattering matrix in Eq. III,
consistent with the cluster decomposition principleweinberg1995quantum ,
comprises of terms where $k<N$ particles are not scattered by the two-level
system (which results in the $\delta(\cdot)$ functions conserving the
individual frequencies) while the remaining $N-k$ particles interact with the
quantum emitter (which is captured by the $N-k$ excitation Green’s function
$\mathcal{G}^{N-k}(\vec{\omega};\vec{\nu})$).

In the following subsections, we explicitly calculate and study the properties
of single- and two-particle scattering matrices. We provide analytical results
for single and two-particle scattering when the two-level system is coupled to
a homogeneous bath (i.e. $\varepsilon(\mathbf{x})=\varepsilon_{0}$ for all
$\mathbf{x}$). We also numerically study properties of the single and two-
particle scattering when the inhomogeneous bath is a 2D particleic crystal
with a point defect. While we restrict our analysis to two-level emitters
throughout this paper, in appendix E we extend the formalism of this paper to
an emitter with a general level structure but under the Markovian
approximation.

###  III.1 Single-particle scattering

We first consider scattering of a single particle propagating in the
inhomogeneous bath from the localized two-level system. The spectral
properties of the single-particle scattering matrix are governed by the
single-excitation Green’s function

$\displaystyle\mathcal{G}^{1}(\omega;\nu)=$
$\displaystyle\int_{t,s\in\mathbb{R}}\bra{g,\text{vac}}\mathcal{T}\big{[}\sigma(t)\sigma^{\dagger}(s)\big{]}\ket{g,\text{vac}}e^{i(\omega
t-\nu s)}dtds.$ (34)

Since $\sigma(t)\ket{\text{vac},g}=0$, it follows that
$\bra{g,\text{vac}}\mathcal{T}\big{[}\sigma(t)\sigma^{\dagger}(s)\big{]}\ket{g,\text{vac}}=A_{e}(t-s)\Theta(t\geq
s)$ where $A_{e}(t)=\bra{e,\text{vac}}e^{-iHt}\ket{e,\text{vac}}$. It then
follows that

$\displaystyle\mathcal{G}^{1}(\omega;\nu)=2\pi\delta(\omega-\nu)\underbrace{\int_{0}^{\infty}A_{e}(t)e^{i\nu
t}dt}_{\mathcal{G}_{0}(\nu)}.$ (35)

As is shown in appendix D.1, $A_{e}(t)$ can be computed by solving the
following (non-Markovian) ODE

$\displaystyle\frac{dA_{e}(t)}{dt}=-i\omega_{0}A_{e}(t)+\frac{V_{0}^{2}}{\pi}\times$
$\displaystyle\qquad\int_{\begin{subarray}{c}t^{\prime}\in[0,t]\\\
\omega\in\mathbb{R}^{+}\end{subarray}}\text{Im}\big{[}\text{G}(\mathbf{x}_{d},\mathbf{x}_{d};\omega)\big{]}e^{-i\omega(t-t^{\prime})}A_{e}(t^{\prime})dt^{\prime}\
d\omega.$ (36)

subject to the boundary condition $A_{e}(t)=0$ for $t\leq 0$. Specializing Eq.
III to the case of single-particle scattering ($N=1$), we obtain

$\displaystyle
S^{1}\big{(}\mathbf{x},\omega;\Omega,\nu\big{)}=N_{0}(\nu)\bigg{(}e^{ik_{0}(\nu)\Omega\cdot\mathbf{x}}+\mathcal{E}_{s}(\mathbf{x},\Omega,\nu)-$
$\displaystyle\qquad
iV_{0}^{2}\mathcal{E}(\mathbf{x}_{d},\Omega,\nu)\text{G}(\mathbf{x},\mathbf{x}_{d};\nu)\mathcal{G}_{0}(\nu)\bigg{)}\delta(\omega-\nu).$
(37)

This form of the single-particle scattering matrix is immediately physically
interpretable — we notice that as a consequence of energy conservation, the
frequencies of the incoming and outgoing single-particles are constrained to
be equal. The spatial dependence of the scattering matrix can be seen to be
the sum of three terms — first is the incident plane-wave
$e^{ik_{0}(\nu)\Omega\cdot\mathbf{x}}$, second is the fields produced by
direct scattering of the incident plane-wave from the inhomogeneous bath,
$\mathcal{E}_{s}(\mathbf{x},\Omega,\nu)$, and third are dipolar fields
radiated by the two-level system due to excitation by the incident plane-wave.
We note that the strength of the dipolar field depends on the coupling
strength $V_{0}$, as well as the total fields at the position of the dipole
$\mathcal{E}(\mathbf{x}_{d},\Omega,\nu)$.

Figure 2: Single-particle scattering from an inhomogeneous bosonic bath. a)
Schematic depiction of the permittivity distribution — a 2D square lattice
photonic crystal cavity structure which is translationally invariant and
infinite along the $z-$axis. The radii of the central (defect) cylinder is
$0.65a$, and of the surrounding cylinders is $0.2a$. The permittivity of the
cylinder is assumed to be 8.9 and the two-level system is assumed to be
positioned at $\mathbf{x}_{d}=(0,0.5a)$ with respect to the center of the
central cylinder. b) Spectral dependence of the imaginary part of the Green’s
function $\text{Im}[G(\mathbf{x}_{d},\mathbf{x}_{d};\omega)]$ at the TLS
position and the total electric field at the TLS position on excitation with a
plane-wave propagating along the $x$ direction. The photonic crystal structure
supports two (lossy) localized modes that are shown as an inset. c) The
single-excitation Green’s function (normalized to its maximum value) as a
function of frequency for different coupling strengths $V_{0}$. The resonant
frequency of the two-level system is assumed to be $\omega_{0}=0.35(2\pi/a)$.
d) The scattering cross-section corresponding to the fields emitted by the
two-level system for different coupling strengths. e) The total scattering
cross-section seen by the incident plane-wave field. f) Far field angular
distribution of the fields scattered by the bath, the two-level system and the
total scattered fields in the $xy$-plane as a function of the angle with
respect to the $x-$axis.

We first consider a homogeneous bath, for which
$\mathcal{E}_{s}(\mathbf{x},\Omega,\nu)=0$, and assume that $\mathbf{x}_{d}=0$
and $\mathcal{E}(\mathbf{x}_{d},\Omega,\nu)=1$. The spatial dependence of the
single-particle scattering matrix is thus only determined by interference
between the incident plane-wave and an outgoing spherical waves emitted by the
two-level system at $\textbf{x}=\textbf{x}_{d}$

$\displaystyle S^{1}\big{(}\mathbf{x},\omega;\Omega,\nu\big{)}=$
$\displaystyle\qquad
N_{0}(\nu)\bigg{(}e^{ik_{0}(\nu)\Omega\cdot\mathbf{x}}+iV_{0}^{2}\frac{e^{ik_{0}(\nu)r}}{4\pi
r}\mathcal{G}_{0}(\nu)\bigg{)}\delta(\omega-\nu).$ (38)

Furthermore, the Weisskopf-Wigner theory wang1974wigner  can be used to
approximate the single-excitation Green’s function by a complex Lorentzian
(refer to appendix D.1)

$\displaystyle\mathcal{G}_{0}(\nu)=\frac{1}{i(\omega_{0}-\nu)+\gamma/2},$ (39)

where
$\gamma=2V_{0}^{2}\text{G}(\mathbf{x}=0,\mathbf{x}^{\prime}=0;\omega_{0})=V_{0}^{2}\omega_{0}\sqrt{\varepsilon_{0}}$.
Our result for the single-particle scattering matrix Eq. III.1 can be used to
calculate the single-particle scattering cross-section by integrating the
square of the amplitude of the scattered field over a spherical surface to
obtain

$\displaystyle\sigma_{\text{3D}}(\nu)=\frac{V_{0}^{4}}{4\pi^{2}}|\mathcal{G}_{0}(\nu)|^{2}=\frac{1}{\pi^{2}\omega_{0}^{2}\varepsilon_{0}}\bigg{[}\frac{\gamma^{2}/4}{(\nu-\omega_{0})^{2}+\gamma^{2}/4}\bigg{]}.$
(40)

This result agrees with the classical scattering cross-section for a dipole
scatterer, and the single-particle scattering cross-section for a dipole two-
level system coupled to a homogeneous bath derived in Ref. liu2016quantum . A
similar analysis yields the following result for the scattering cross-section
for a 2D problem

$\displaystyle\sigma_{\text{2D}}(\nu)=\frac{4}{\omega_{0}\sqrt{\varepsilon_{0}}}\bigg{[}\frac{\gamma^{2}/4}{(\nu-\omega_{0})^{2}+\gamma^{2}/4}\bigg{]},$
(41)

where $\gamma=V_{0}^{2}/2$. We note that at resonance ($\omega=\omega_{0}$),
both $\sigma_{\text{3D}}(\nu)$ and $\sigma_{\text{2D}}(\nu)$ are independent
of the coupling strength $V_{0}$, and are completely determined by the
resonant frequency of the two-level system.

As an example of an inhomogeneous bath, we consider a 2D photonic crystal
structure with a point defect (Fig. 2a). This structure supports two resonant
modes that interact with the two-level system (Fig. 2b). The single-excitation
Green’s function of the two-level system interacting with the photonic crystal
structure computed from Eqs. 35 and III.1 is shown in Fig. 2c (refer to
appendix D.1 for computational details). We observe that for low coupling
strengths $V_{0}$, the single-excitation Green’s function is approximately a
Lorentzian in frequency and shows multiple resonant features at higher
coupling strengths — this is consistent with polaritonic splitting seen in
strongly coupled cavity QED systems. To gain more insight into the
modification of the scattering properties of the two-level system due to its
interaction with the photonic crystal structure, we calculate the total
scattering cross-section corresponding to fields scattered at frequency
$\omega$ by the two-level system, $\sigma_{\text{TLS}}(\omega;\Omega)$, the
bath, $\sigma_{\text{bath}}(\omega;\Omega)$, and both the two-level system and
bath, $\sigma_{\text{full}}(\omega;\Omega)$, when excited with an incident
plane-wave propagating along $\Omega$,

$\displaystyle\sigma_{\text{TLS}}(\omega;\Omega)=V_{0}^{4}|\mathcal{E}(\mathbf{x}_{d},\Omega;\omega)|^{2}|\mathcal{G}_{0}(\omega)|^{2}\times$
$\displaystyle\qquad\qquad\bigg{[}\lim_{R\to\infty}\int_{\mathbb{S}^{1}}R\big{|}\text{G}(R\Omega^{\prime},\mathbf{x}_{d};\omega)\big{|}^{2}d\Omega^{\prime}\bigg{]},$
(42a)
$\displaystyle\sigma_{\text{bath}}(\omega;\Omega)=\lim_{R\to\infty}\int_{\mathbb{S}^{1}}R\big{|}\mathcal{E}_{s}(R\Omega^{\prime},\mathbf{x}_{d};\omega)\big{|}^{2}d\Omega^{\prime},$
(42b)
$\displaystyle\sigma_{\text{full}}(\omega;\Omega)=\lim_{R\to\infty}\int_{\mathbb{S}^{1}}R\bigg{|}\mathcal{E}_{s}(R\Omega^{\prime},\mathbf{x}_{d};\omega)-$
$\displaystyle\qquad\qquad
iV_{0}^{2}\mathcal{E}(\mathbf{x}_{d},\Omega;\omega)\mathcal{G}_{0}(\omega)\text{G}(R\Omega^{\prime},\mathbf{x}_{d};\omega)\bigg{|}^{2}d\Omega^{\prime}.$
(42c)

Figures 2d shows $\sigma_{\text{TLS}}(\omega;\hat{x})$ as a function of
frequency $\omega$ — for small coupling strength, this scattering cross
section resembles a single Lorentzian similar to the result in the homogeneous
bath Eq. 41 (the maximum homogeneous scattering cross-section,
$4/\omega_{0}\sqrt{\varepsilon_{0}}$, is indicated with the black dashed
line). For larger coupling strengths, this scattering cross-section is no
longer a Lorentzian in frequency due to the strong frequency dependence of the
inhomogeneous bath and can become significantly larger than the scattering
cross-section in homogeneous bath. Figure 2e shows
$\sigma_{\text{full}}(\omega;\hat{x})$ as a function of frequency $\omega$ —
for small coupling strengths, this is dominated by scattering from the
inhomogeneous bath (i.e.
$\sigma_{\text{full}}(\omega;\hat{x})\approx\sigma_{\text{bath}}(\omega;\hat{x})$
which is shown in the dashed black line) while for larger coupling strengths,
the contribution of the two-level system is observable. The angular
distribution of the scattered single particle in the far field is shown in
Fig. 2f.

###  III.2 Two-particle scattering

The two-particle spectrum of the scattered wave-packet is governed by the two-
excitation Green’s function $\mathcal{G}^{2}(\vec{\omega};\vec{\nu})$, which
is given by

$\displaystyle\mathcal{G}^{2}(\vec{\omega};\vec{\nu})=$
$\displaystyle\int_{\begin{subarray}{c}t_{1},t_{2}\in\mathbb{R}\\\
s_{1},s_{2}\in\mathbb{R}\end{subarray}}\bra{g,\text{vac}}\mathcal{T}\big{[}\sigma(t_{1})\sigma(t_{2})\sigma^{\dagger}(s_{1})\sigma^{\dagger}(s_{2})\big{]}\ket{g,\text{vac}}\times$
$\displaystyle\qquad\qquad\prod_{i=1}^{2}e^{i(\omega_{i}t_{i}-\nu_{i}s_{i})}dt_{i}ds_{i}.$
(43)

As is shown in Appendix D.2, following the procedure introduced in Ref.
shi2015multiphoton , this Green’s function can be obtained entirely from the
single-excitation Green’s function

$\displaystyle\mathcal{G}^{2}(\vec{\omega};\vec{\nu})=\sum_{\mathcal{P}_{2}}\mathcal{G}^{1}(\omega_{1};\nu_{\mathcal{P}_{2}(1)})\mathcal{G}^{1}(\omega_{2};\nu_{\mathcal{P}_{2}(2)})$
$\displaystyle\quad\underbrace{-\pi\frac{\mathcal{G}_{0}(\nu_{1})\mathcal{G}_{0}(\nu_{2})\mathcal{G}_{0}(\omega_{1})\mathcal{G}_{0}(\omega_{2})}{\Gamma(\omega_{1}+\omega_{2})}\delta\bigg{(}\sum_{j=1}^{2}\big{(}\omega_{j}-\nu_{j}\big{)}\bigg{)}}_{\mathcal{G}^{2}_{\text{C}}(\vec{\omega};\vec{\nu})},$
(44a) where $\mathcal{P}_{2}$ is a two element permutation and
$\displaystyle\Gamma(E)=\int_{\mathbb{R}}\mathcal{G}_{0}(\omega)\mathcal{G}_{0}(E-\omega)\frac{d\omega}{2\pi}.$
(44b)

Specializing Eq. III provides the following expression for the two-particle
scattering matrix

$\displaystyle
S^{2}\big{(}\vec{\mathbf{x}},\vec{\omega};\vec{\Omega},\vec{\nu}\big{)}=\sum_{\mathcal{P}_{2}}\bigg{(}S\big{(}\mathbf{x}_{1},\omega_{1};\Omega_{\mathcal{P}(1)},\nu_{\mathcal{P}_{2}(1)}\big{)}\times$
$\displaystyle\qquad
S\big{(}\mathbf{x}_{2},\omega_{2};\Omega_{\mathcal{P}(2)},\nu_{\mathcal{P}_{2}(2)}\big{)}\bigg{)}+S_{\text{C}}^{2}\big{(}\vec{\mathbf{x}},\vec{\omega};\vec{\Omega},\vec{\nu}\big{)},$
(45a) where $\mathcal{P}$ is a permutation of the set $\\{1,2\\}$,
$S^{1}\big{(}\mathbf{x},\omega;\Omega,\nu)$ is the single-particle scattering
matrix and
$S_{\text{C}}^{2}(\vec{\mathbf{x}},\vec{\omega};\vec{\Omega},\vec{\nu})$ is
the connected part of the two-particle scattering matrix that is given by
$\displaystyle
S_{\text{C}}^{2}\big{(}\vec{\mathbf{x}},\vec{\omega};\vec{\Omega},\vec{\nu}\big{)}=-\frac{V_{0}^{4}}{4\pi^{2}}\bigg{(}\prod_{i=1}^{2}N_{0}(\nu_{i})\mathcal{E}(\mathbf{x}_{d},\Omega_{i},\nu_{i})\times$
$\displaystyle\qquad\qquad\text{G}(\mathbf{x}_{i},\mathbf{x}_{d};\omega_{i})\bigg{)}\mathcal{G}_{\text{C}}^{2}(\vec{\omega};\vec{\nu}).$
(45b)

We point out that Eq. 45 is consistent with the form of the two-particle
scattering matrix as expected by the cluster-decomposition principle
weinberg1995quantum , with the connected part of the scattering matrix
capturing the particle-particle interactions induced by the nonlinearity of
the two-level system. As can be seen from Eq. 45, the spectral properties of
the entanglement generated in the scattered two-particles due to such
particle-particle interactions is captured by the connected part of the two-
excitation Green’s function, while the spatial distribution of the
entanglement is simply given by the Green’s function corresponding to the
permittivity distribution.

For a homogeneous bath ($\varepsilon(\mathbf{x})=\varepsilon_{0}$) and
assuming $\mathbf{x}_{d}=0$, the connected part of the scattering matrix
reduces to a product of two outgoing radially symmetric spherical waves

$\displaystyle
S_{\text{C}}^{2}\big{(}\vec{\mathbf{x}},\vec{\omega};\vec{\Omega},\vec{\nu}\big{)}=-\frac{V_{0}^{4}}{4\pi^{2}}\bigg{(}\prod_{i=1}^{2}N_{0}(\nu_{i})\frac{e^{ik_{0}(\omega_{i})r_{i}}}{4\pi
r_{i}}\bigg{)}\mathcal{G}_{\text{C}}^{2}(\vec{\omega};\vec{\nu}),$ (46)

where, in the Weisskopf-Wigner approximation,
$\mathcal{G}_{\text{C}}(\vec{\omega};\vec{\nu})$ is given by (refer to
appendix D.2 for details)

$\displaystyle\mathcal{G}_{\text{C}}^{2}(\vec{\omega};\vec{\nu})=-\frac{i\pi(\sum_{j=1}^{2}\omega_{j}-2\omega_{0}-i\gamma)\delta(\sum_{j=1}^{2}(\omega_{j}-\nu_{j}))}{\prod_{j=1}^{2}(\nu_{j}-\omega_{0}-i\gamma/2)(\omega_{j}-\omega_{0}-i\gamma/2)},$
(47)

where $\gamma=V_{0}^{2}\omega_{0}\sqrt{\varepsilon_{0}}$. We point out that
this result for the connected part of the two-excitation Green’s function is
identical to the one obtained in Waveguide QED under a Markovian approximation
to the waveguide two-level system interaction fan2010input .

Figure 3: Two-particle scattering from a localized two-level system embedded
in the inhomogeneous bath shown in Fig. 2a. We assume both the incident
particles to be propagating along the $x$ axis (i.e.
$\Omega_{1}=\Omega_{2}=\hat{x}$) for our simulations. a) The connected part of
the two-excitation Green’s function as a function of the frequency difference
between the two scattered particles assuming that the incident particles are
at frequency $\nu_{0}$ resonant with the two-level system
($\nu_{0}=\omega_{0}$). b) The frequency spectrum and c) the far field angular
distribution of the first 5 Schmidt vectors of the output state corresponding
to the connected part of the two-particle scattering matrix when both the
input particles are Gaussian wave-packets in frequency (Eq. 48) with
$\tau=220/\omega_{0}$.

Next, we consider two-particle scattering from the two-level system when it is
embedded in the photonic crystal structure shown in Fig. 2a. Figure 3b shows
the dependence of the connected part of the two-excitation Green’s function on
the frequency difference of the output particles computed using Eq. 47. To
study the spectral and spatial entanglement properties of the outgoing two-
particles, we consider exciting the two-level system with a two-particle Fock
state with both particles propagating along the $\hat{{x}}$ axis but with a
Gaussian spectrum:

$\displaystyle\ket{\psi_{\text{in}}}\propto\bigg{(}\int_{\mathbb{R}^{+}}e^{-(\nu-\omega_{0})^{2}\tau^{2}/2}a^{\dagger}(\hat{x},\nu)d\nu\bigg{)}^{2}\ket{\text{vac},g},$
(48)

where $\tau$ governs the spectral bandwidth of the incident wave-packet. For
large $\omega_{0}\tau\gg 1$, application of the connected part of the two-
particle scattering matrix on this wave-packet yields the following two-
particle wave-packet $\psi_{\text{out,C}}(\vec{\mathbf{x}},\vec{\omega})$

$\displaystyle\psi_{\text{out,C}}(\vec{\mathbf{x}},\vec{\omega})\approx-\frac{V_{0}^{4}N_{0}(\omega_{0})\mathcal{E}^{2}(\mathbf{x}_{d},\hat{x},\omega_{0})\mathcal{G}_{0}^{2}(\omega_{0})}{4\pi^{2}}\times$
$\displaystyle\qquad\bigg{(}\prod_{i=1}^{2}\mathcal{G}_{0}(\omega_{i})\text{G}(\mathbf{x}_{i},\mathbf{x}_{d};\omega_{i})\bigg{)}\frac{e^{-(\omega_{1}+\omega_{2}-2\omega_{0})^{2}\tau^{2}}}{\Gamma(\omega_{1}+\omega_{2})}.$
(49)

The Schmidt decomposition of
$\psi_{\text{out,C}}(\vec{\mathbf{x}},\vec{\omega})$ can then express it as a
superposition of two-particle fock states

$\displaystyle\psi_{\text{out,C}}(\vec{\mathbf{x}},\vec{\omega})=\sum_{k=1}^{\infty}\lambda_{i}\psi_{i}(\mathbf{x}_{1},\omega_{1})\psi_{i}(\mathbf{x}_{2},\omega_{2}),$
(50)

with $\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}\dots$. Figures 3b and c show
the spectral and spatial dependence of the first five Schmidt components of
$\psi_{\text{out,C}}$ for a two-level system strongly coupled to the photonic
crystal bath. The frequency spectrum of each Schmidt vector (shown in Fig. 3b)
is computed by integrating the magnitude square of the Schmidt vector over a
sufficiently large circle surrounding the photonic crystal structure and the
angular distribution of the far field (shown in Fig. 3c) is computed by
integrating the magnitude square of the Schmidt vector with respect to
frequency. The multi-Lorentzian nature of the spectrum of the emitted
particles can be attributed to a strong hybridization between the two-level
system and the electromagnetic modes of the photonic crystal cavity structure,
while the far field pattern is almost completely determined by the field
profile of the photonic crystal cavity mode.

##  IV Conclusion

In conclusion, we presented an approach to exactly analyze spatial and
spectral properties of few-particle scattering from a localized system
embedded in a an inhomogeneous bath. We develop exact solutions for the single
and two-particle scattering matrices when the localized system is a two-level
system. The formalism presented in this paper paves the way for
computationally understanding the impact of structured media on quantum
scattering from embedded localized systems. While the formalism in this paper
assumes that the bosonic bath can be described by an inhomogeneous scalar wave
equation, it can be extended to quantum optical systems where the bosonic bath
(i.e. the photonic field) interacting with the localized system is described
by a vector wave-equation.

##  V Acknowledgements

We thank David A. B. Miller and Geun Ho Ahn for providing feedback on the
manuscript. R. T. acknowledges support from Kailath Graduate fellowship. This
work is funded by the Air Force Office of Scientific Research under AFOSR MURI
programme (award no. FA9550-17-1-0002).

## References

  * (1) P. Lodahl, S. Mahmoodian, and S. Stobbe, “Interfacing single photons and single quantum dots with photonic nanostructures,” _Reviews of Modern Physics_ , vol. 87, no. 2, p. 347, 2015. [Online]. Available: https://doi.org/10.1103/RevModPhys.87.347
  * (2) X. Ding, Y. He, Z.-C. Duan, N. Gregersen, M.-C. Chen, S. Unsleber, S. Maier, C. Schneider, M. Kamp, S. Höfling _et al._ , “On-demand single photons with high extraction efficiency and near-unity indistinguishability from a resonantly driven quantum dot in a micropillar,” _Physical Review Letters_ , vol. 116, no. 2, p. 020401, 2016. [Online]. Available: https://doi.org/10.1103/PhysRevLett.116.020401
  * (3) P. Michler, Ed., _Quantum Dots for Quantum Information Technologies_. Springer, 2017. [Online]. Available: https://doi.org/10.1007/978-3-319-56378-7
  * (4) P. Senellart, G. Solomon, and A. White, “High-performance semiconductor quantum-dot single-photon sources,” _Nature Nanotechnology_ , vol. 12, no. 11, pp. nnano–2017, 2017. [Online]. Available: https://doi.org/10.1038/nnano.2017.218
  * (5) J. L. Zhang, S. Sun, M. J. Burek, C. Dory, Y.-K. Tzeng, K. A. Fischer, Y. Kelaita, K. G. Lagoudakis, M. Radulaski, Z.-X. Shen _et al._ , “Strongly cavity-enhanced spontaneous emission from silicon-vacancy centers in diamond,” _Nano letters_ , vol. 18, no. 2, pp. 1360–1365, 2018. [Online]. Available: https://doi.org/10.1021/acs.nanolett.7b05075
  * (6) P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowling, and G. J. Milburn, “Linear optical quantum computing with photonic qubits,” _Reviews of Modern Physics_ , vol. 79, no. 1, p. 135, 2007. [Online]. Available: https://doi.org/10.1103/RevModPhys.79.135
  * (7) J. L. O’brien, A. Furusawa, and J. Vučković, “Photonic quantum technologies,” _Nature Photonics_ , vol. 3, no. 12, p. 687, 2009. [Online]. Available: https://doi.org/10.1038/nphoton.2009.229
  * (8) D. Roy, C. M. Wilson, and O. Firstenberg, “Colloquium: Strongly interacting photons in one-dimensional continuum,” _Reviews of Modern Physics_ , vol. 89, no. 2, p. 021001, 2017. [Online]. Available: https://doi.org/10.1103/RevModPhys.89.021001
  * (9) A. Reiserer and G. Rempe, “Cavity-based quantum networks with single atoms and optical photons,” _Reviews of Modern Physics_ , vol. 87, no. 4, p. 1379, 2015\. [Online]. Available: https://doi.org/10.1103/RevModPhys.87.1379
  * (10) L.-M. Duan and C. Monroe, “Colloquium: Quantum networks with trapped ions,” _Reviews of Modern Physics_ , vol. 82, no. 2, p. 1209, 2010. [Online]. Available: https://doi.org/10.1103/RevModPhys.82.1209
  * (11) N. Sangouard, C. Simon, H. De Riedmatten, and N. Gisin, “Quantum repeaters based on atomic ensembles and linear optics,” _Reviews of Modern Physics_ , vol. 83, no. 1, p. 33, 2011. [Online]. Available: https://doi.org/10.1103/RevModPhys.83.33
  * (12) K. Nemoto, M. Trupke, S. J. Devitt, A. M. Stephens, B. Scharfenberger, K. Buczak, T. Nöbauer, M. S. Everitt, J. Schmiedmayer, and W. J. Munro, “Photonic architecture for scalable quantum information processing in diamond,” _Physical Review X_ , vol. 4, no. 3, p. 031022, 2014. [Online]. Available: https://doi.org/10.1103/PhysRevX.4.031022
  * (13) D. Englund, D. Fattal, E. Waks, G. Solomon, B. Zhang, T. Nakaoka, Y. Arakawa, Y. Yamamoto, and J. Vučković, “Controlling the spontaneous emission rate of single quantum dots in a two-dimensional photonic crystal,” _Physical review letters_ , vol. 95, no. 1, p. 013904, 2005. [Online]. Available: https://doi.org/10.1103/PhysRevLett.95.013904
  * (14) R. S. Daveau, K. C. Balram, T. Pregnolato, J. Liu, E. H. Lee, J. D. Song, V. Verma, R. Mirin, S. W. Nam, L. Midolo _et al._ , “Efficient fiber-coupled single-photon source based on quantum dots in a photonic-crystal waveguide,” _Optica_ , vol. 4, no. 2, pp. 178–184, 2017\. [Online]. Available: https://doi.org/10.1364/OPTICA.4.000178
  * (15) M. Pelton, C. Santori, J. Vucković, B. Zhang, G. S. Solomon, J. Plant, and Y. Yamamoto, “Efficient source of single photons: a single quantum dot in a micropost microcavity,” _Physical review letters_ , vol. 89, no. 23, p. 233602, 2002. [Online]. Available: https://doi.org/10.1103/PhysRevLett.89.233602
  * (16) Y.-M. He, Y. He, Y.-J. Wei, D. Wu, M. Atatüre, C. Schneider, S. Höfling, M. Kamp, C.-Y. Lu, and J.-W. Pan, “On-demand semiconductor single-photon source with near-unity indistinguishability,” _Nature Nanotechnology_ , vol. 8, no. 3, pp. 213–217, 2013. [Online]. Available: https://doi.org/10.1038/nnano.2012.262
  * (17) K. A. Fischer, L. Hanschke, J. Wierzbowski, T. Simmet, C. Dory, J. J. Finley, J. Vučković, and K. Müller, “Signatures of two-photon pulses from a quantum two-level system,” _Nature Physics_ , vol. 13, no. 7, pp. 649–654, 2017. [Online]. Available: https://doi.org/10.1038/nphys4052
  * (18) L. Hanschke, K. A. Fischer, S. Appel, D. Lukin, J. Wierzbowski, S. Sun, R. Trivedi, J. Vučković, J. J. Finley, and K. Müller, “Quantum dot single-photon sources with ultra-low multi-photon probability,” _npj Quantum Information_ , vol. 4, no. 1, pp. 1–6, 2018. [Online]. Available: https://doi.org/10.1038/s41534-018-0092-0
  * (19) M. Metcalfe, S. M. Carr, A. Muller, G. S. Solomon, and J. Lawall, “Resolved sideband emission of inas/gaas quantum dots strained by surface acoustic waves,” _Physical review letters_ , vol. 105, no. 3, p. 037401, 2010. [Online]. Available: https://doi.org/10.1103/PhysRevLett.105.037401
  * (20) K. C. Miao, A. Bourassa, C. P. Anderson, S. J. Whiteley, A. L. Crook, S. L. Bayliss, G. Wolfowicz, G. Thiering, P. Udvarhelyi, V. Ivády _et al._ , “Electrically driven optical interferometry with spins in silicon carbide,” _Science Advances_ , vol. 5, no. 11, p. eaay0527, 2019\. [Online]. Available: https://doi.org/10.1126/sciadv.aay0527
  * (21) D. M. Lukin, A. D. White, R. Trivedi, M. A. Guidry, N. Morioka, C. Babin, Ö. O. Soykal, J. Ul-Hassan, N. T. Son, T. Ohshima _et al._ , “Spectrally reconfigurable quantum emitters enabled by optimized fast modulation,” _npj Quantum Information_ , vol. 6, no. 1, pp. 1–9, 2020. [Online]. Available: https://doi.org/10.1038/s41534-020-00310-0
  * (22) H. Wang, Y.-M. He, T.-H. Chung, H. Hu, Y. Yu, S. Chen, X. Ding, M.-C. Chen, J. Qin, X. Yang _et al._ , “Towards optimal single-photon sources from polarized microcavities,” _Nature Photonics_ , vol. 13, no. 11, pp. 770–775, 2019. [Online]. Available: https://doi.org/10.1038/s41566-019-0494-3
  * (23) X. Li, J. Chen, P. Voss, J. Sharping, and P. Kumar, “All-fiber photon-pair source for quantum communications: Improved generation of correlated photons,” _Optics express_ , vol. 12, no. 16, pp. 3737–3744, 2004. [Online]. Available: https://doi.org/10.1364/OPEX.12.003737
  * (24) G. Harder, V. Ansari, B. Brecht, T. Dirmeier, C. Marquardt, and C. Silberhorn, “An optimized photon pair source for quantum circuits,” _Optics express_ , vol. 21, no. 12, pp. 13 975–13 985, 2013. [Online]. Available: https://doi.org/10.1364/OE.21.013975
  * (25) H. Pichler and P. Zoller, “Photonic circuits with time delays and quantum feedback,” _Physical review letters_ , vol. 116, no. 9, p. 093601, 2016. [Online]. Available: https://doi.org/10.1103/PhysRevLett.116.093601
  * (26) S. Xu and S. Fan, “Generate tensor network state by sequential single-photon scattering in waveguide qed systems,” _APL Photonics_ , vol. 3, no. 11, p. 116102, 2018. [Online]. Available: https://doi.org/10.1063/1.5044248
  * (27) H. Carmichael, _An open systems approach to quantum optics: lectures presented at the Université Libre de Bruxelles, October 28 to November 4, 1991_. Springer Science & Business Media, 2009, vol. 18. [Online]. Available: https://doi.org/10.1007/978-3-540-47620-7
  * (28) C. Gardiner and P. Zoller, _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_. Springer Science & Business Media, 2004, vol. 56. [Online]. Available: https://doi.org/10.1007/9783540223016
  * (29) A. Asenjo-Garcia, J. Hood, D. Chang, and H. Kimble, “Atom-light interactions in quasi-one-dimensional nanostructures: A green’s-function perspective,” _Physical Review A_ , vol. 95, no. 3, p. 033818, 2017. [Online]. Available: https://doi.org/10.1103/PhysRevA.95.033818
  * (30) H. T. Dung, L. Knöll, and D.-G. Welsch, “Resonant dipole-dipole interaction in the presence of dispersing and absorbing surroundings,” _Physical Review A_ , vol. 66, no. 6, p. 063810, 2002. [Online]. Available: https://doi.org/10.1103/PhysRevA.66.063810
  * (31) A. H. Kiilerich and K. Mølmer, “Input-output theory with quantum pulses,” _Physical review letters_ , vol. 123, no. 12, p. 123604, 2019. [Online]. Available: https://doi.org/10.1103/PhysRevLett.123.123604
  * (32) J. R. Taylor, _Scattering theory: the quantum theory of nonrelativistic collisions_. Courier Corporation, 2006\.
  * (33) S. Fan, Ş. E. Kocabaş, and J.-T. Shen, “Input-output formalism for few-photon transport in one-dimensional nanophotonic waveguides coupled to a qubit,” _Physical Review A_ , vol. 82, no. 6, p. 063821, 2010. [Online]. Available: https://doi.org/10.1103/PhysRevA.82.063821
  * (34) S. Xu, E. Rephaeli, and S. Fan, “Analytic properties of two-photon scattering matrix in integrated quantum systems determined by the cluster decomposition principle,” _Physical review letters_ , vol. 111, no. 22, p. 223602, 2013\. [Online]. Available: https://doi.org/10.1103/PhysRevLett.111.223602
  * (35) S. Xu and S. Fan, “Input-output formalism for few-photon transport: A systematic treatment beyond two photons,” _Physical Review A_ , vol. 91, no. 4, p. 043845, 2015. [Online]. Available: https://doi.org/10.1103/PhysRevA.91.043845
  * (36) R. Trivedi, K. Fischer, S. Xu, S. Fan, and J. Vuckovic, “Few-photon scattering and emission from low-dimensional quantum systems,” _Physical Review B_ , vol. 98, no. 14, p. 144112, 2018. [Online]. Available: https://doi.org/10.1103/PhysRevB.98.144112
  * (37) R. Trivedi, D. Malz, S. Sun, S. Fan, and J. Vučković, “Optimal two-photon excitation of bound states in non-markovian waveguide qed,” _Physical Review A_ , vol. 104, no. 1, p. 013705, 2021.
  * (38) K. A. Fischer, R. Trivedi, V. Ramasesh, I. Siddiqi, and J. Vučković, “Scattering into one-dimensional waveguides from a coherently-driven quantum-optical system,” _Quantum_ , vol. 2, p. 69, 2018. [Online]. Available: https://doi.org/10.22331/q-2018-05-28-69
  * (39) T. Shi, D. E. Chang, and J. I. Cirac, “Multiphoton-scattering theory and generalized master equations,” _Physical Review A_ , vol. 92, no. 5, p. 053834, 2015. [Online]. Available: https://doi.org/10.1103/PhysRevA.92.053834
  * (40) T. Caneva, M. T. Manzoni, T. Shi, J. S. Douglas, J. I. Cirac, and D. E. Chang, “Quantum dynamics of propagating photons with strong interactions: a generalized input–output formalism,” _New Journal of Physics_ , vol. 17, no. 11, p. 113001, 2015. [Online]. Available: http://dx.doi.org/10.1088/1367-2630/17/11/113001
  * (41) S. Haroche and D. Kleppner, “Cavity quantum electrodynamics,” _Phys. Today_ , vol. 42, no. 1, pp. 24–30, 1989. [Online]. Available: http://dx.doi.org/10.1063/1.881201
  * (42) R. Trivedi, M. Radulaski, K. A. Fischer, S. Fan, and J. Vučković, “Photon blockade in weakly driven cavity quantum electrodynamics systems with many emitters,” _Physical review letters_ , vol. 122, no. 24, p. 243602, 2019. [Online]. Available: https://doi.org/10.1103/PhysRevLett.122.243602
  * (43) T. Shi, S. Fan _et al._ , “Two-photon transport through a waveguide coupling to a whispering-gallery resonator containing an atom and photon-blockade effect,” _Physical Review A_ , vol. 87, no. 6, p. 063818, 2013. [Online]. Available: https://doi.org/10.1103/PhysRevA.87.063818
  * (44) J. Liu, M. Zhou, and Z. Yu, “Quantum scattering theory of a single-photon fock state in three-dimensional spaces,” _Optics letters_ , vol. 41, no. 18, pp. 4166–4169, 2016. [Online]. Available: https://doi.org/10.1364/OL.41.004166
  * (45) M. Orenstein, S. Fan _et al._ , “Two-level quantum system as a macroscopic scatterer for ultraconfined two-dimensional photonic modes,” _Physical Review A_ , vol. 102, no. 1, p. 013717, 2020. [Online]. Available: https://doi.org/10.1103/PhysRevA.102.013717
  * (46) E. K. Miller, “A selective survey of computational electromagnetics,” _IEEE Transactions on Antennas and Propagation_ , vol. 36, no. 9, pp. 1281–1305, 1988. [Online]. Available: https://doi.org/10.1109/8.8607
  * (47) S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, “Inverse design in nanophotonics,” _Nature Photonics_ , vol. 12, no. 11, pp. 659–670, 2018. [Online]. Available: https://doi.org/10.1038/s41566-018-0246-9
  * (48) H. T. Dung, L. Knöll, and D.-G. Welsch, “Three-dimensional quantization of the electromagnetic field in dispersive and absorbing inhomogeneous dielectrics,” _Physical Review A_ , vol. 57, no. 5, p. 3931, 1998. [Online]. Available: https://doi.org/10.1103/PhysRevA.57.3931
  * (49) T. Gruner and D.-G. Welsch, “Green-function approach to the radiation-field quantization for homogeneous and inhomogeneous kramers-kronig dielectrics,” _Physical Review A_ , vol. 53, no. 3, p. 1818, 1996. [Online]. Available: https://doi.org/10.1103/PhysRevA.53.1818
  * (50) S. Weinberg, _The quantum theory of fields_. Cambridge university press, 1995, vol. 2. [Online]. Available: https://doi.org/10.1017/CBO9781139644167
  * (51) H.-P. Breuer and F. Petruccione, _The theory of open quantum systems_. Oxford University Press on Demand, 2002. [Online]. Available: https://doi.org/10.1093/acprof:oso/9780199213900.001.0001
  * (52) Y. Wang and I. Khoo, “On the wigner-weisskopf approximation in quantum optics,” _Optics Communications_ , vol. 11, no. 4, pp. 323–326, 1974. [Online]. Available: https://doi.org/10.1016/0030-4018(74)90227-2

##  Appendix A Second quantization with noise operators

We begin by outlining some useful properties of the scalar Green’s function
corresponding to the inhomogeneous wave-equation, followed by showing that the
noise-operator based quantum theory presented in the main text satisfies the
canonical commutation relations stemming from Dirac’s quantization.

###  A.1 Properties of the scalar inhomogeneous Green’s function

Throughout this section,
$\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)$ will refer to the
scalar’s green’s function which satisfies

$\displaystyle\big{[}\nabla^{2}_{\mathbf{x}}+\omega^{2}\varepsilon_{\alpha}(\mathbf{x},\omega)\big{]}\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)=\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime}),$
(51)

where $\varepsilon_{\alpha}(\omega)=\varepsilon(\mathbf{x})+i\alpha(\omega)$.
Here $\varepsilon(\mathbf{x})$ is the inhomogeneous permittivity distribution
corresponding to the inhomogenuous media and $\alpha(\omega)$ governs the loss
in the media. We point out that while physically $\omega>0$, it is convenient
to consider the solution of Eq. 51 for $\omega<0$ as well. We assume that for
$\omega>0,\alpha(-\omega)=-\alpha(\omega)$ — this ensures that for
$\omega<0,\alpha(\omega)<0$ and hence the outgoing solution to Eq. 51 for
$\omega<0$ decays exponentially as $|\mathbf{x}|\to\infty$. Furthermore, we
also make a physically reasonable assumption of $\alpha(\omega)\to 0$ as
$|\omega|\to\infty$.

Property 1:

$\displaystyle\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)=\text{G}_{\alpha}^{*}(\mathbf{x},\mathbf{x}^{\prime};-\omega).$
(52)

_Proof_ : This immediately follows by conjugating Eq. 51 and by noting that
$\varepsilon_{\alpha}(\mathbf{x},-\omega)=\varepsilon_{\alpha}^{*}(\mathbf{x},\omega)$.

Property 2:

$\displaystyle
2i\omega^{2}\alpha\int_{\mathbb{R}^{3}}\text{G}_{\alpha}^{*}(\mathbf{x}^{\prime},\mathbf{y};\omega)\text{G}_{\alpha}(\mathbf{x},\mathbf{y};\omega)d^{3}\mathbf{y}=\text{G}_{\alpha}^{*}(\mathbf{x},\mathbf{x}^{\prime};\omega)-\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega).$
(53)

_Proof_ : We note that the scalar Green’s function is necessarily reciprocal
i.e.
$\text{G}_{\alpha}(\mathbf{x},\mathbf{y};\omega)=\text{G}_{\alpha}(\mathbf{y},\mathbf{x};\omega)$.
Therefore, it follows from Eq. 51 that

$\displaystyle\big{[}\nabla^{2}_{\mathbf{y}}+\omega^{2}\varepsilon_{\alpha}(\mathbf{y};\omega)\big{]}\text{G}_{\alpha}(\mathbf{x},\mathbf{y};\omega)=\delta^{3}(\mathbf{x}-\mathbf{y}).$
(54)

Multiplying this equation with
$\text{G}_{\alpha}^{*}(\mathbf{x}^{\prime},\mathbf{y};\omega)$ and integrating
with respect to $\mathbf{y}$ we obtain

$\displaystyle\text{G}_{\alpha}^{*}(\mathbf{x},\mathbf{x}^{\prime};\omega)=-\int_{\mathbb{R}^{3}}\nabla_{\mathbf{y}}\text{G}_{\alpha}^{*}(\mathbf{x}^{\prime},\mathbf{y};\omega)\cdot\nabla_{\mathbf{y}}\text{G}_{\alpha}(\mathbf{x},\mathbf{y};\omega)d^{3}\mathbf{y}+\omega^{2}\int_{\mathbb{R}^{3}}\varepsilon_{\alpha}(\mathbf{y};\omega)\text{G}_{\alpha}^{*}(\mathbf{x}^{\prime},\mathbf{y};\omega)\text{G}_{\alpha}(\mathbf{x},\mathbf{y};\omega)d^{3}\mathbf{y}.$
(55)

Conjugating Eq. 55 followed by swapping $\mathbf{x}$ and
$\mathbf{x}^{\prime}$, we obtain

$\displaystyle\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)=-\int_{\mathbb{R}^{3}}\nabla_{\mathbf{y}}\text{G}_{\alpha}^{*}(\mathbf{x}^{\prime},\mathbf{y};\omega)\cdot\nabla_{\mathbf{y}}\text{G}_{\alpha}(\mathbf{x},\mathbf{y};\omega)d^{3}\mathbf{y}+\omega^{2}\int_{\mathbb{R}^{3}}\varepsilon_{\alpha}^{*}(\mathbf{y};\omega)\text{G}_{\alpha}^{*}(\mathbf{x}^{\prime},\mathbf{y};\omega)\text{G}_{\alpha}(\mathbf{x},\mathbf{y};\omega)d^{3}\mathbf{y}.$
(56)

Subtracting Eq. 56 from Eq. 55, we obtain Eq. 53.

Property 3: Assuming that $\exists\ \varepsilon_{\text{ub}}$ such that
$\varepsilon(\mathbf{x})<\varepsilon_{\text{UB}}\ \forall\ \mathbf{x}$ and
$\varepsilon(\mathbf{x})>0\ \forall\ \mathbf{x}$,

$\displaystyle\int_{-\infty}^{\infty}\omega\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)\
d\omega=-\frac{i\pi}{\varepsilon(\mathbf{x})}\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime}).$
(57)

_Proof_ : Defining
$\theta(\mathbf{x})=\varepsilon_{\text{ub}}-\varepsilon(\mathbf{x})$, Eq. 51
can be rewritten as

$\displaystyle\big{[}\nabla^{2}_{\mathbf{x}}+\omega^{2}(\varepsilon_{\text{ub}}+i\alpha(\omega))\big{]}\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)=\omega^{2}\theta(\mathbf{x})\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)+\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime}).$
(58)

A Dyson series for $\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)$
can be deduced from Eq. 58 to obtain

$\displaystyle\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)=\text{G}_{\alpha}^{\text{ub}}(\mathbf{x},\mathbf{x}^{\prime};\omega)+\sum_{n=1}^{\infty}\omega^{2k}\int_{\mathbf{y}_{1},\mathbf{y}_{2}\dots\mathbf{y}_{k}\in\mathbb{R}^{3}}\bigg{[}\prod_{k=0}^{n}\text{G}_{\alpha}^{\text{ub}}(\mathbf{y}_{k},\mathbf{y}_{k+1};\omega)\bigg{]}_{\begin{subarray}{c}\mathbf{y}_{0}=\mathbf{x}\\\
\mathbf{y}_{n+1}=\mathbf{x}^{\prime}\end{subarray}}\bigg{[}\prod_{k=1}^{n}\theta(\mathbf{y}_{k})d^{3}\mathbf{y}_{k}\bigg{]},$
(59)

where $\text{G}_{\alpha}^{\text{ub}}(\mathbf{x},\mathbf{x}^{\prime};\omega)$
satisfies

$\displaystyle\big{[}\nabla_{\mathbf{x}}^{2}+\omega^{2}(\varepsilon_{\text{ub}}+i\alpha(\omega))\big{]}\text{G}_{\alpha}^{\text{ub}}(\mathbf{x},\mathbf{x}^{\prime};\omega)=\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime}).$
(60)

It can easily be seen from Eq. 60 that

$\displaystyle\text{G}_{\alpha}^{\text{ub}}(\mathbf{x},\mathbf{x}^{\prime};\omega)=\int_{\mathbb{R}^{3}}\frac{e^{i\mathbf{k}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}}{\omega^{2}\varepsilon_{\text{ub}}-|\mathbf{k}|^{2}+i\omega^{2}\alpha(\omega)}\frac{d^{3}\mathbf{k}}{(2\pi)^{3}}.$
(61)

We next evaluate the integral in Eq. 57 using the expansion for
$\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)$ in Eq. 59. Consider
the following integral

$\displaystyle\int_{-\infty}^{\infty}\omega^{2k+1}$
$\displaystyle\bigg{[}\prod_{k=0}^{n}\text{G}_{\alpha}^{\text{ub}}(\mathbf{y}_{k},\mathbf{y}_{k+1};\omega)\bigg{]}d\omega=\int_{\omega\in\mathbb{R}}\int_{\mathbf{k}_{0},\mathbf{k}_{1}\dots\mathbf{k}_{n}\in\mathbb{R}^{3}}\omega^{2k+1}\bigg{[}\prod_{k=0}^{n}\frac{e^{i\mathbf{k}_{k}\cdot(\mathbf{y}_{k}-\mathbf{y}_{k+1})}}{\omega^{2}\varepsilon_{\text{ub}}-|\mathbf{k}_{k}|^{2}+i\omega^{2}\alpha(\omega)}\frac{d^{3}\mathbf{k}_{k}}{(2\pi)^{3}}\bigg{]}d\omega.$
(62)

Noting that the integrand is analytic as a function of $\omega$ in the upper
complex half-plane, it follows that the integral with respect to $\omega$ can
be evaluated along $\omega=Re^{i\theta}$, with $R\to\infty$ with $\theta$
going from $\pi$ to $0$. This immediately yields

$\displaystyle\int_{-\infty}^{\infty}\omega^{2k+1}\bigg{[}\prod_{k=0}^{n}\text{G}_{\alpha}^{\text{ub}}(\mathbf{y}_{k},\mathbf{y}_{k+1};\omega)\bigg{]}d\omega=-\frac{i\pi}{\varepsilon_{\text{ub}}}\prod_{k=0}^{n}\delta^{3}(\mathbf{y}_{k}-\mathbf{y}_{k+1}).$
(63)

It thus follows from Eq. 59 that

$\displaystyle\int_{-\infty}^{\infty}\omega\text{G}_{\alpha}^{\text{ub}}(\mathbf{x},\mathbf{x}^{\prime};\omega)=-i\pi\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime})\sum_{k=0}^{\infty}\frac{\theta^{k}(\mathbf{x})}{\varepsilon_{\text{ub}}^{k+1}}=-\frac{i\pi}{\varepsilon(\mathbf{x})}\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime}).$
(64)

###  A.2 Verifying canonical commutation relations between $\psi(\mathbf{x})$
and $\pi(\mathbf{x})$

We begin by deriving the commutator
$[\psi(\mathbf{x},\omega),\psi^{\dagger}(\mathbf{x}^{\prime},\omega^{\prime})]$.
Using Eq. II.2.1 and the commutator for the noise operator
$\psi(\mathbf{x},\omega)$, we obtain

$\displaystyle[\psi(\mathbf{x},\omega),\psi^{\dagger}(\mathbf{x}^{\prime};\omega^{\prime})]=\frac{1}{\pi}\omega^{2}\alpha(\omega)\delta(\omega-\omega^{\prime})\int_{\mathbb{R}^{3}}\text{G}_{\alpha}^{*}(\mathbf{x}^{\prime},\mathbf{y};\omega)\text{G}_{\alpha}(\mathbf{x},\mathbf{y};\omega)d^{3}\mathbf{y}.$
(65)

Using Eq. 53, we obtain

$\displaystyle[\psi(\mathbf{x},\omega),\psi^{\dagger}(\mathbf{x}^{\prime};\omega^{\prime})]$
$\displaystyle=\frac{\delta(\omega-\omega^{\prime})}{2\pi
i}\big{[}\text{G}_{\alpha}^{*}(\mathbf{x},\mathbf{x}^{\prime};\omega)-\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)\big{]},$
$\displaystyle=-\frac{\delta(\omega-\omega^{\prime})}{2\pi}\text{Im}\big{[}\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)\big{]}.$
(66)

Consider now the commutator $[\psi(\mathbf{x}),\pi(\mathbf{x}^{\prime})]$:
Using Eqs. 9a and A.2, we obtain

$\displaystyle[\psi(\mathbf{x}),\pi(\mathbf{x}^{\prime})]=\frac{\varepsilon(\mathbf{x}^{\prime})}{\pi}\int_{\omega\in\mathbb{R}^{+}}\omega\big{[}\text{G}_{\alpha}^{*}(\mathbf{x},\mathbf{x}^{\prime};\omega)-\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)\big{]}d\omega.$
(67)

Using Eqs. 52 and 57, it immediately follows that
$[\psi(\mathbf{x}),\pi(\mathbf{x}^{\prime})]=i\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime})$
— the noise operator based quantization thus reproduces the canonical
commutation relations for inhomogeneous bosonic baths.

##  Appendix B Evaluating $a_{\alpha}(\mathbf{x},\Omega,\omega)$ as
$\alpha(\omega)\to 0$

In this appendix, we detail the evaluation of the $\lim_{\alpha(\omega)\to
0}a_{\alpha}(\mathbf{x},\Omega;\omega)$, where
$a_{\alpha}(\mathbf{x},\Omega;\omega)$ is defined in Eq. II.2.2. In
particular, we show that in the limit of $\alpha(\omega)\to 0$,
$a_{\alpha}(\mathbf{x},\Omega;\omega)$ becomes independent of $\mathbf{x}$.
Since $a_{\alpha}(\mathbf{x},\Omega;\omega)$ are, by construction, bosonic
annihilation operators, in order to show that $\lim_{\alpha(\omega)\to
0}a_{\alpha}(\mathbf{x},\Omega;\omega)$ is independent of $\mathbf{x}$, it is
sufficient to show that $\lim_{\alpha(\omega)\to
0}[a_{\alpha}(\mathbf{x},\Omega;\omega),a_{\alpha}^{\dagger}(\mathbf{x}^{\prime},\Omega^{\prime};\omega^{\prime})]$
is independent of $\mathbf{x}$ and $\mathbf{x}^{\prime}$. Using Eq. II.2.2 and
$[\phi_{\alpha}(\mathbf{x};\omega),\phi^{\dagger}_{\alpha}(\mathbf{x}^{\prime};\omega^{\prime})]=\delta^{3}(\mathbf{x}-\mathbf{x}^{\prime})\delta(\omega-\omega^{\prime})$,
we can explicitly evaluate this commutator to obtain

$\displaystyle[a_{\alpha}(\mathbf{x},$
$\displaystyle\Omega;\omega),a_{\alpha}^{\dagger}(\mathbf{x}^{\prime},\Omega^{\prime};\omega^{\prime})]=$
$\displaystyle\frac{\alpha(\omega)\omega}{4\pi^{4}\sqrt{\varepsilon_{0}}}\delta(\omega-\omega^{\prime})\int_{\begin{subarray}{c}\mathbf{y}\in\mathbb{R}^{3}\\\
k,k^{\prime}\in\mathbb{R}^{+}\end{subarray}}\frac{k^{2}k^{\prime
2}e^{i(k\Omega-k^{\prime}\Omega^{\prime})\cdot\mathbf{y}}e^{i(k-k_{0}(\omega))\Omega\cdot\mathbf{x}}e^{-i(k^{\prime}-k_{0}(\omega))\Omega^{\prime}\cdot\mathbf{x}^{\prime}}}{(\omega^{2}\varepsilon_{0}-k^{2}+i\omega^{2}\alpha(\omega))(\omega^{2}\varepsilon_{0}-k^{\prime
2}-i\omega^{2}\alpha(\omega))}dkdk^{\prime}d^{3}\mathbf{y}.$ (68)

Noting that

$\displaystyle\int_{\mathbb{R}^{3}}e^{i(k\Omega-k^{\prime}\Omega^{\prime})\cdot\mathbf{y}}d^{3}\mathbf{y}=\frac{8\pi^{3}\delta(k-k^{\prime})\delta^{2}(\Omega-\Omega^{\prime})}{k^{2}},$
(69)

we obtain

$\displaystyle[a_{\alpha}(\mathbf{x},\Omega;\omega),a_{\alpha}^{\dagger}(\mathbf{x}^{\prime},\Omega^{\prime};\omega^{\prime})]=\frac{2\alpha(\omega)\omega}{\pi\sqrt{\varepsilon_{0}}}\delta(\omega-\omega^{\prime})\delta^{2}(\Omega-\Omega^{\prime})\int_{\mathbb{R}^{+}}\frac{k^{2}e^{i(k-k_{0}(\omega))\Omega\cdot(\mathbf{x}-\mathbf{x}^{\prime})}}{(\omega^{2}\varepsilon_{0}-k^{2})^{2}+\omega^{4}\sigma^{2}(\omega)}dk.$
(70)

We note that in the limit of $\alpha(\omega)\to 0$, the integral in the above
equation can be approximated as follows

$\displaystyle\int_{\mathbb{R}^{+}}\frac{k^{2}e^{i(k-k_{0}(\omega))\Omega\cdot(\mathbf{x}-\mathbf{x}^{\prime})}}{(\omega^{2}\varepsilon_{0}-k^{2})^{2}+\omega^{2}\sigma^{4}(\omega)}dk\approx\frac{1}{4}\int_{\mathbb{R}^{+}}\frac{1}{(k-\omega\sqrt{\varepsilon_{0}})^{2}+\omega^{2}\sigma^{2}(\omega)/4\varepsilon_{0}}\approx\frac{\pi\sqrt{\varepsilon_{0}}}{2\omega\alpha(\omega)}.$
(71)

Consequently, we obtain

$\displaystyle\lim_{\alpha(\omega)\to
0}[a_{\alpha}(\mathbf{x},\Omega;\omega),a_{\alpha}^{\dagger}(\mathbf{x}^{\prime},\Omega^{\prime};\omega^{\prime})]=\delta(\omega-\omega^{\prime})\delta^{2}(\Omega-\Omega^{\prime}).$
(72)

Thus, we conclude that the operators $a_{\sigma}(\mathbf{x},\Omega;\omega)$
are independent of $\mathbf{x}$ in the limit of $\alpha(\omega)\to 0$.
Defining $a(\Omega;\omega)$ via $a(\Omega;\omega)=\lim_{\alpha(\omega)\to
0}a_{\alpha}(\mathbf{x},\Omega;\omega)$, we thus obtain the following
commutator

$\displaystyle[a(\Omega;\omega),a^{\dagger}(\Omega^{\prime};\omega^{\prime})]=\delta(\omega-\omega^{\prime})\delta^{2}(\Omega-\Omega^{\prime}).$
(73)

##  Appendix C Relating scattering matrix to the localized system’s Green’s
functions

Our starting point is the Heisenberg picture representation of the
$N-$particle scattering matrix element in Eq. III. We first substitute for
$\psi(\mathbf{x},\omega;t_{f})$ from Eq. 29 to obtain

$\displaystyle
e^{i\sum_{n=1}^{N}\omega_{n}t_{f}}\bra{\text{vac},g}\mathcal{T}\bigg{[}\prod_{i=1}^{N}\psi(\mathbf{x}_{i},\omega_{i};t_{f})\prod_{i=1}^{N}a^{\dagger}(\Omega_{i},\nu_{i};t_{i})\bigg{]}\ket{\text{vac},g}$
$\displaystyle=\bra{\text{vac},g}\sum_{k=0}^{N}\sum_{\mathcal{B}_{k}^{N}}\mathcal{T}\bigg{[}\bigg{(}\prod_{n\in\bar{\mathcal{B}}_{k}^{N}}\frac{V_{0}}{\pi}\text{Im}\big{[}\text{G}(\mathbf{x}_{n},\mathbf{x}_{d};\omega_{n})\big{]}\int_{t_{i}}^{t_{f}}\sigma(\tau)e^{i\omega_{n}\tau}d\tau\bigg{)}\bigg{(}\prod_{n\in\mathcal{B}_{k}^{N}}\psi(\mathbf{x}_{n},\omega_{n};t_{i})e^{i\omega_{n}t_{i}}\bigg{)}\times$
$\displaystyle\bigg{(}\prod_{n=1}^{N}a^{\dagger}(\Omega_{n},\nu_{n};t_{i})\bigg{)}\bigg{]}\ket{\text{vac},g},$
(74)

where $\mathcal{B}_{k}^{N}$ is an unordered $k-$element subset of
$\\{1,2,\dots N\\}$, and $\bar{\mathcal{B}}_{k}^{N}$ is its complement. Using
the commutator
$[\psi(\mathbf{x},\omega),a^{\dagger}(\Omega,\nu)]=N_{0}(\nu)\mathcal{E}(\mathbf{x},\Omega,\nu)\delta(\omega-\nu)$
and the relation $\psi(\mathbf{x},\omega;t_{i})\ket{\text{vac},g}=0$, we
obtain

$\displaystyle\bigg{(}\prod_{n\in\mathcal{B}_{k}^{N}}\psi(\mathbf{x}_{n},\omega_{n};t_{i})\bigg{)}\prod_{n=1}^{N}a^{\dagger}(\Omega_{n},\nu_{n};t_{i})\ket{\text{vac},g}=$
$\displaystyle\sum_{\mathcal{D}_{k}^{N},\mathcal{P}_{k}}\bigg{(}\prod_{n=1}^{N-k}a^{\dagger}(\Omega_{\bar{\mathcal{D}}_{k}^{N}(n)},\nu_{\bar{\mathcal{D}}_{k}^{N}(n)};t_{i})\bigg{)}\bigg{(}\prod_{n=1}^{k}N_{0}\big{(}\nu_{\mathcal{D}_{k}^{N}(n)}\big{)}\mathcal{E}\big{(}\mathbf{x}_{\mathcal{P}_{k}\mathcal{B}_{k}^{N}(n)},\Omega_{\mathcal{D}_{k}^{N}(n)};\nu_{\mathcal{D}_{k}^{N}(n)}\big{)}\delta\big{(}\omega_{\mathcal{P}_{k}\mathcal{B}_{k}^{N}(n)}-\nu_{\mathcal{D}_{k}^{N}(n)}\big{)}\bigg{)},$
(75)

where $\mathcal{D}_{k}^{N}$ is another unordered $k-$element subset of
$\\{1,2\dots N\\}$ with $\bar{\mathcal{D}}^{N}_{k}$ being its complement,
$\mathcal{P}_{k}$ is a $k-$ element permutation and
$\mathcal{P}_{k}\mathcal{D}_{k}^{N}$ is a permutation of the $k-$element
subset $\mathcal{D}_{k}^{N}$. Next, we substitute Eq. C into Eq. C and use Eq.
III to substitute for $a^{\dagger}(\Omega,\nu;t_{i})$ in terms of
$a^{\dagger}(\Omega,\nu;t_{f})$. We note that due to the time-ordering
operator, $\hat{a}^{\dagger}(\Omega,\nu;t_{f})$ annihilate
$\bra{\text{vac},g}$. With this simplification, we obtain Eq. III in the main
text.

##  Appendix D Computing single and two-excitation Green’s functions

In this section, we outline the numerical procedures used for computing the
single and two-excitation Green’s function that are needed for calculating the
single and two-excitation scattering matrices.

###  D.1 Single-excitation Green’s function

The single-excitation Green’s function in time-domain,
$\mathcal{G}(t;t^{\prime})$, is given by

$\displaystyle\mathcal{G}(t,s)=\bra{g,\text{vac}}\mathcal{T}\big{[}\sigma(t)\sigma^{\dagger}(s)\big{]}\ket{g,\text{vac}},$
(76)

where $\sigma(t)$ and $\sigma^{\dagger}(t)$ are the Heisenberg picture
operators corresponding to the two-level system’s lowering and raising
operators. Since $\sigma(t)\ket{g,\text{vac}}=0$ and
$\bra{g,\text{vac}}\sigma^{\dagger}(t)=0$ for all $t$, we obtain an
alternative representation of the single-excitation Green’s function

$\displaystyle\mathcal{G}(t,s)=\bra{e,\text{vac}}e^{-iH(t-s)}\ket{e,\text{vac}}\Theta(t\geq
t^{\prime}),$ (77)

i.e. it can be obtained by solving for the dynamics of an initially excited
emitter emitting into the bosonic bath. To this end, we consider the system to
initially be in the state $\ket{\psi(0)}=\sigma^{\dagger}\ket{G}$ — using the
noise operator representation of the bath, the following ansatz can be assumed
for the state of the system

$\ket{\psi(t)}=\bigg{[}A_{e}(t)\sigma^{\dagger}+\int_{\mathbf{x}\in\mathbb{R}^{3}}\int_{\omega\in\mathbb{R}^{+}}A_{g}(\mathbf{x},\omega;t)\phi^{\dagger}_{\alpha}(\mathbf{x},\omega)d^{3}\mathbf{x}\
d\omega\bigg{]}\ket{G},$ (78)

with $A_{e}(0)=1$ and $A_{g}(\mathbf{x},\omega;0)=0$. From the Schroedinger’s
equation, it then immediately follows that

$\displaystyle i\frac{dA_{e}(t)}{dt}=\omega_{0}$ $\displaystyle
A_{e}(t)+iV_{0}\int_{\mathbf{x}\in\mathbb{R}^{3}}\int_{\omega\in\mathbb{R}^{+}}\bigg{(}\frac{\omega^{2}\alpha(\omega)}{\pi}\bigg{)}^{1/2}\text{G}_{\alpha}(\mathbf{x}_{d},\mathbf{x};\omega)A_{g}(\mathbf{x},\omega;t)d^{3}\mathbf{x}d\omega,$
(79) $\displaystyle i\frac{dA_{g}(\mathbf{x},\omega;t)}{dt}$
$\displaystyle=\omega
A_{g}(\mathbf{x},\omega;t)-iV_{0}\bigg{(}\frac{\omega^{2}\alpha(\omega)}{\pi}\bigg{)}^{1/2}\text{G}_{\alpha}^{*}(\mathbf{x}_{d},\mathbf{x};\omega)A_{e}(t).$
(80)

Integrating the second equation, we obtain

$A_{g}(\mathbf{x},\omega;t)=-V_{0}\bigg{(}\frac{\omega^{2}\alpha(\omega)}{\pi}\bigg{)}^{1/2}\text{G}_{\alpha}^{*}(\mathbf{x}_{d},\mathbf{x};\omega)\int_{0}^{t}A_{e}(t^{\prime})e^{-i\omega(t-t^{\prime})}dt^{\prime}.$
(81)

Substituting this back into Eq. 79 and using Eq. 53, we thus obtain

$\frac{dA_{e}(t)}{dt}=-i\omega_{0}A_{e}(t)+\frac{V_{0}^{2}}{\pi}\int_{t^{\prime}=0}^{t}\int_{\omega=0}^{\infty}\text{Im}\big{[}\text{G}_{\alpha}(\mathbf{x}_{d},\mathbf{x}_{d};\omega)\big{]}e^{-i\omega(t-t^{\prime})}A_{e}(t^{\prime})dt^{\prime}\
d\omega.$ (82)

We note that while the above analysis assumes a loss $\alpha(\omega)$ in the
dielectric distribution under consideration, it can be taken to $0$ in Eq. 82
since $\lim_{\alpha(\omega)\to
0}\text{Im}\big{[}\text{G}_{\alpha}(\mathbf{x}_{d},\mathbf{x}_{d};\omega)]$
exists.

_Numerical solution_ : For a given permittivity distribution
$\varepsilon(\mathbf{x})$,
Im$\big{[}G(\mathbf{x}_{d},\mathbf{x}_{d};\omega)\big{]}$ can be computed
using various numerical methods (such as FDTD or FDFD). To compute $A_{e}(t)$,
we approximate Im$\big{[}G(\mathbf{x}_{d},\mathbf{x}_{d};\omega)\big{]}$ as a
sum of Lorentzians within a sufficiently large bandwidth around the emitter
frequency

$\text{Im}\big{[}\text{G}(\mathbf{x}_{d},\mathbf{x}_{d};\omega)\big{]}=-\sum_{n}\frac{p_{n}}{(\omega-\omega_{n})^{2}+\gamma_{n}^{2}}.$
(83)

Eq. 82 can the be expressed as

$\displaystyle\frac{dA_{e}(t)}{dt}=-i\omega_{0}A_{e}(t)-\frac{V_{0}^{2}}{\pi}\sum_{n}\xi_{n}(t),$
(84)

where

$\displaystyle\xi_{n}(t)=\int_{t^{\prime}=0}^{t}\int_{0}^{\infty}\frac{p_{n}e^{-i\omega(t-t^{\prime})}}{(\omega-\omega_{n})^{2}+\gamma_{n}^{2}}A_{n}(t^{\prime})\
dt^{\prime}d\omega.$ (85)

Assuming that $\omega_{n}/\gamma_{n}\gg 1$,

$\displaystyle\int_{0}^{\infty}\frac{p_{n}e^{-i\omega(t-t^{\prime})}}{(\omega-\omega_{n})^{2}+\gamma_{n}^{2}}d\omega\approx\int_{-\infty}^{\infty}\frac{p_{n}e^{-i\omega(t-t^{\prime})}}{(\omega-\omega_{n})^{2}+\gamma_{n}^{2}}d\omega=\frac{\pi
p_{n}}{\gamma_{n}}e^{-i\omega_{n}(t-t^{\prime})}e^{-\gamma_{n}|t-t^{\prime}|},$
(86)

which immediately yields

$\displaystyle\xi_{n}(t)=\frac{\pi
p_{n}}{\gamma_{n}}\int_{0}^{t}e^{-i\omega_{n}(t-t^{\prime})}e^{-\gamma_{n}(t-t^{\prime})}A_{n}(t^{\prime})dt^{\prime}.$
(87)

A differential equation for $\xi_{n}(t)$ can be derived from this expression
by differentiating it with respect to $t$

$\displaystyle\frac{d\xi_{n}(t)}{dt}=-(i\omega_{n}+\gamma_{n})\xi_{n}(t)+\frac{\pi
p_{n}}{\gamma_{n}}A_{e}(t).$ (88)

In frequency-domain scattering matrix calculations, the quantity of interest
is the one-sided fourier transform of $A_{e}(t)$

$\displaystyle A_{e}(\omega)=\int_{0}^{\infty}A_{e}(t)e^{i\omega t}dt.$ (89)

Taking the one-sided fourier transform of Eqs. 84 and 88, we obtain an
expression for $A_{e}(\omega)$

$\displaystyle
A_{e}(\omega)=\bigg{[}i(\omega_{0}-\omega)+V_{0}^{2}\sum_{n}\frac{p_{n}/\gamma_{n}}{i(\omega_{n}-\omega)+\gamma_{n}}\bigg{]}^{-1}.$
(90)

_Weisskopf-Wigner approximation_ : For permittivity distribution
$\varepsilon(\mathbf{x})$ with a frequency-response much broader than the
line-width of the emitter, Eq. 82 can be approximated using the Weisskopf-
Wigner formalism

$\displaystyle\frac{dA_{e}(t)}{dt}=-\bigg{(}i\omega_{0}+\frac{\gamma}{2}\bigg{)}A_{e}(t),$
(91)

where
$\gamma=-2V_{0}^{2}\text{Im}\big{[}G(\mathbf{x}_{d},\mathbf{x}_{d};\omega_{0})]$
is the spontaneous emission decay rate of the emitter into the inhomogeneous
bath. This yields an exponentially decaying solution for $A_{e}(t)$ and a
complex Lorentzian for $A_{e}(\omega)$

$\displaystyle A_{e}(t)=e^{-i\omega_{0}t}e^{-\gamma t/2}\text{ and
}A_{e}(\omega)=\frac{1}{i(\omega_{0}-\omega)+\gamma/2}.$ (92)

###  D.2 Two-excitation Green’s function

The evaluation of the two-excitation Green’s function for the two-level system
can be done by following the procedure outlined in Ref. shi2015multiphoton .
As an initial step, the two-level system is replaced with a boson with an
onsite repulsion as described by the Hamiltonian

$\displaystyle
H_{\text{0}}=\omega_{0}\sigma^{\dagger}\sigma+\frac{U_{0}}{2}(\sigma^{\dagger})^{2}\sigma^{2},$
(93)

instead of $\omega_{0}\sigma^{\dagger}\sigma$ and the commutation relation
$[\sigma,\sigma^{\dagger}]=\sigma_{z}$ replaced by the bosonic commutation
$[\sigma,\sigma^{\dagger}]=1$. The two-level system is recovered from this
model by taking $U_{0}\to\infty$. Next, we decompose the full hamiltonian of
the TLS coupled to the inhomogeneous bath, $H$, as

$\displaystyle H=H_{\text{L}}+H_{\text{NL}},$ (94)

where
$H_{\text{L}}=H_{\text{bath}}+H_{\text{int}}+\omega_{0}\sigma^{\dagger}\sigma$
and $H_{\text{NL}}=U_{0}(\sigma^{\dagger})^{2}\sigma^{2}/2$. We point out that
$H_{\text{L}}$ is a quadratic hamiltonian when expressed in terms of the
bosonic annihilation operators of the bath and $\sigma$, and consequently its
dynamics are completely captured by two-point correlations between the
different bosonic operators. We now consider the time-ordered expectation
$\bra{\text{vac},g}\mathcal{T}\big{[}\sigma(t_{1})\sigma(t_{2})\sigma^{\dagger}(s_{1})\sigma^{\dagger}(s_{2})\big{]}\ket{\text{vac},g}$
and rewrite it as

$\displaystyle\bra{\text{vac},g}\mathcal{T}\big{[}\sigma(t_{1})\sigma(t_{2})$
$\displaystyle\sigma^{\dagger}(s_{1})\sigma^{\dagger}(s_{2})\big{]}\ket{\text{vac},g}=$
$\displaystyle\bra{\text{vac},g}\mathcal{T}\bigg{[}\tilde{\sigma}(t_{1})\tilde{\sigma}(t_{2})\tilde{\sigma}^{\dagger}(s_{1})\tilde{\sigma}^{\dagger}(s_{2})\exp\bigg{(}{-i\int_{-\infty}^{\infty}\tilde{H}_{\text{NL}}(t)dt}\bigg{)}\bigg{]}\ket{\text{vac},g},$
(95)

where, for an operator $O$,
$\tilde{O}(t)=e^{iH_{\text{L}}t}Oe^{-iH_{\text{L}}t}$. The exponential in Eq.
D.2 can be expanded to obtain

$\displaystyle\bra{\text{vac},g}\mathcal{T}\big{[}$
$\displaystyle\sigma(t_{1})\sigma(t_{2})\sigma^{\dagger}(s_{1})\sigma^{\dagger}(s_{2})\big{]}\ket{\text{vac},g}=\bra{\text{vac},g}\mathcal{T}\big{[}\tilde{\sigma}(t_{1})\tilde{\sigma}(t_{2})\tilde{\sigma}^{\dagger}(s_{1})\tilde{\sigma}^{\dagger}(s_{2})\big{]}\ket{\text{vac},g}+$
$\displaystyle\sum_{k=1}^{\infty}\int_{\tau_{1},\tau_{2}\dots\tau_{k}\in\mathbb{R}}\frac{1}{k!}\bigg{(}-\frac{iU_{0}}{2}\bigg{)}^{k}\bra{\text{vac},g}\mathcal{T}\bigg{[}\tilde{\sigma}(t_{1})\tilde{\sigma}(t_{2})\tilde{\sigma}^{\dagger}(s_{1})\tilde{\sigma}^{\dagger}(s_{2})\prod_{i=1}^{k}\big{(}\tilde{\sigma}^{\dagger}(\tau_{k}))^{2}\tilde{\sigma}^{2}(\tau_{k})\bigg{]}\ket{\text{vac},g}\prod_{i=1}^{k}d\tau_{i}.$
(96)

Each of the time-ordered expectations in the summation in the above equation
can be evaluated through an application of the Wick’s theorem since the
operators involved in the expectation are Heisenberg operators with respect to
a quadratic Hamiltonian $H_{\text{L}}$

$\displaystyle\bra{\text{vac},g}\mathcal{T}\big{[}\tilde{\sigma}(t_{1})\tilde{\sigma}(t_{2})\tilde{\sigma}^{\dagger}(s_{1})\tilde{\sigma}^{\dagger}(s_{2})\big{]}\ket{\text{vac},g}=g(t_{1};s_{1})g(t_{2};s_{2})+g(t_{1};s_{2})g(t_{2};s_{1}),$
(97a)
$\displaystyle\bra{\text{vac},g}\mathcal{T}\bigg{[}\tilde{\sigma}(t_{1})\tilde{\sigma}(t_{2})\tilde{\sigma}^{\dagger}(s_{1})\tilde{\sigma}^{\dagger}(s_{2})\prod_{i=1}^{k}\big{(}\tilde{\sigma}^{\dagger}(\tau_{k}))^{2}\tilde{\sigma}^{2}(\tau_{k})\bigg{]}\ket{\text{vac},g}=$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \
2^{k-1}\sum_{\mathcal{P}_{k}}g(t_{1};\tau_{\mathcal{P}_{k}(1)})g(t_{2};\tau_{\mathcal{P}_{k}(1)})g(s_{1};\tau_{\mathcal{P}_{k}(k)})g(s_{2};\tau_{\mathcal{P}_{k}(k)})\prod_{i=1}^{k-1}g^{2}(\tau_{\mathcal{P}_{k}(i)};\tau_{\mathcal{P}_{k}(i+1)}).$
(97b)

where $\mathcal{P}_{k}$ is a $k-$element permutation of $\\{1,2\dots k\\}$ and
$g(t;s)=\bra{\text{vac},g}\mathcal{T}[\tilde{\sigma}(t)\tilde{\sigma}^{\dagger}(s)]\ket{\text{vac},g}$.
We note that $H_{\text{NL}}$ does not impact the dynamics of the system within
the single-excitation subspace, and consequently
$g(t;s)=\bra{\text{vac},g}\mathcal{T}[\sigma(t)\sigma^{\dagger}(s)]\ket{\text{vac},g}=A_{e}(t-s)\Theta(t\geq
s)$ where $A_{e}(t)$ is defined in appendix D.1. Finally, we compute the two-
excitation Green’s function
$\mathcal{G}(\omega_{1},\omega_{2};\nu_{1},\nu_{2})$ defined in Eq. III.2 by
Fourier transforming the time-ordered expectation in Eq. D.2 and using Eqs.
D.2 and 97. We thus obtain

$\displaystyle\mathcal{G}^{2}(\vec{\omega};\vec{\nu})=\mathcal{G}^{1}(\omega_{1};\nu_{1})\mathcal{G}^{2}(\omega_{2};\nu_{2})+\mathcal{G}^{1}(\omega_{1};\nu_{2})\mathcal{G}^{2}(\omega_{2};\nu_{1})+\mathcal{G}_{C}^{2}(\vec{\omega};\vec{\nu}),$
(98)

where $\mathcal{G}^{1}(\omega;\nu)$ is the single-excitation Green’s function
and $\mathcal{G}^{2}_{C}(\vec{\omega};\vec{\nu})$ is the connected part of the
two-excitation Green’s function

$\displaystyle\mathcal{G}_{C}^{2}(\vec{\omega};\vec{\nu})=-i\pi\frac{U_{0}\mathcal{G}_{0}(\nu_{1})\mathcal{G}_{0}(\nu_{2})\mathcal{G}_{0}(\omega_{1})\mathcal{G}_{0}(\omega_{2})}{1+iU_{0}\Gamma(\omega_{1}+\omega_{2})}\delta(\omega_{1}+\omega_{2}-\nu_{1}-\nu_{2}),$
(99)

where $\mathcal{G}_{0}(\omega)$ is defined in Eq. 35 and

$\displaystyle\Gamma(E)=\int_{0}^{\infty}A_{e}^{2}(t)e^{iEt}dt=\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathcal{G}_{0}(E-\omega)\mathcal{G}_{0}(\omega)d\omega.$
(100)

Finally, in order to obtain the result for a two-level system, we take the
limit of $U_{0}\to\infty$, which yields the following expression for
$\mathcal{G}_{C}^{2}(\omega_{1},\omega_{2};\nu_{1},\nu_{2})$

$\displaystyle\mathcal{G}_{C}^{2}(\omega_{1},\omega_{2};\nu_{1},\nu_{2})=-\pi\frac{\mathcal{G}_{0}(\nu_{1})\mathcal{G}_{0}(\nu_{2})\mathcal{G}_{0}(\omega_{1})\mathcal{G}_{0}(\omega_{2})}{\Gamma(\omega_{1}+\omega_{2})}\delta(\omega_{1}+\omega_{2}-\nu_{1}-\nu_{2}).$
(101)

_Weisskopf-Wigner approximation_ : Within the Weisskopf-Wigner approximation,

$\displaystyle\mathcal{G}_{0}(\omega)=\frac{1}{i(\omega_{0}-\omega)+\gamma/2}.$
(102)

Therefore, it follows from Eq. 100 that

$\displaystyle\Gamma(E)=\frac{1}{i(2\omega_{0}-E)+\gamma},$ (103)

with which we obtain

$\displaystyle\mathcal{G}_{C}^{2}(\vec{\omega};\vec{\nu})=-\frac{i\pi(\omega_{1}+\omega_{2}-2\omega_{0}-i\gamma)\delta(\omega_{1}+\omega_{2}-\nu_{1}-\nu_{2})}{(\nu_{1}-\omega_{0}-i\gamma/2)(\nu_{2}-\omega_{0}-i\gamma/2)(\omega_{1}-\omega_{0}-i\gamma/2)(\omega_{2}-\omega_{0}-i\gamma/2)}.$
(104)

##  Appendix E Extension to Markovian localized systems with general level
structures

The results presented in Section III can easily be generalized to localized
quantum systems with complex level structures (e.g. $\Lambda$-systems, V-level
systems etc.) under the Markovian approximation. We consider a localized
system with Hamiltonian $H_{s}$ interacting with the inhomogeneous bosonic
bath with an interaction Hamiltonian identical to that in Eq. III with
$\sigma$ being the operator through which the low-dimensional system couples
to the bosonic bath. Since the Heisenberg equations of motion for the bath
operators (Eq. 29) are independent of the level structure of the localized
system, the $N-$particle scattering matrix element defined in Eq. 26 can still
be related to the Green’s function of the low-dimensional system using Eq.
III. In this section, we show that under the Markovian approximation, the low
dimensional system’s Green’s functions are completely determined by the
propagator corresponding to a dissipative (non-Hermitian) Hamiltonian defined
entirely within the system’s Hilbert space. This is an extension of a similar
well known result in waveguide QED xu2015input ; trivedi2018few ;
caneva2015quantum  to the case of an inhomogeneous high-dimensional bath.

For notational convenience and generality, we restrict our attention to
evaluating the following time-ordered expectation value

$\displaystyle\mathcal{E}(t_{1},t_{2}\dots
t_{M})=\bra{\text{vac},g}\mathcal{T}\bigg{[}\prod_{i=1}^{M}c_{i}(t_{i})\bigg{]}\ket{\text{vac},g},$
(105)

where $c_{i}$ are operators defined on the localized system’s Hilbert space
and $\ket{g}$ is the ground state of the localized system. We note that the
$N-$excitation Green’s function defined in Eq. III is simply the Fourier
transform of such a time-ordered expectation with the system operators $c_{i}$
being either $\sigma$ or $\sigma^{\dagger}$. Choosing to describe the Hilbert
space of the bath via the noise operators $\phi_{\alpha}(\mathbf{x},\omega)$,
this expectation can be expressed as the following path integral
weinberg1995quantum

$\displaystyle\mathcal{E}(t_{1},t_{2}\dots
t_{M})=\int_{\begin{subarray}{c}\phi_{\alpha}(\mathbf{x},\omega;t_{i})=0\\\
\phi_{\alpha}(\mathbf{x},\omega;t_{f})=0\end{subarray}}\mathcal{D}[\phi_{\alpha}(\mathbf{x},\omega;t),\phi^{*}(\mathbf{x},\omega;t)]\mathcal{D}s\
\bigg{[}\prod_{i=1}^{M}c_{i}(t_{i})\bigg{]}e^{-i(\mathcal{S}_{\text{sys}}+\mathcal{S}_{\text{bath}}+\mathcal{S}_{\text{int}})},$
(106)

where $t_{i}\to-\infty$, $t_{f}\to\infty$,
$\mathcal{D}[\phi_{\alpha}(\mathbf{x},\omega;t),\phi^{*}(\mathbf{x},\omega;t)]$
is the path integral measure over the Hilbert space of the bath,
$\mathcal{D}s$ is the path integral measure over the Hilbert space of the
localized system, $\mathcal{S}_{\text{sys}}$ is the classical action for the
localized system from $t=t_{i}$ to $t=t_{f}$ and

$\displaystyle\mathcal{S}_{\text{bath}}=\int_{t=t_{i}}^{t_{f}}\int_{\mathbf{x}\in\mathbb{R}^{3}}\int_{\omega\in\mathbb{R}^{+}}\bigg{(}\phi^{*}(\mathbf{x},\omega;t)\frac{\partial}{\partial
t}\phi_{\alpha}(\mathbf{x},\omega;t)-\omega\big{|}\phi_{\alpha}(\mathbf{x},\omega;t)\big{|}^{2}\bigg{)}dt\
d^{3}\mathbf{x}\ d\omega,$ (107a)
$\displaystyle\mathcal{S}_{\text{int}}=-iV_{0}\int_{t=t_{i}}^{t_{f}}\int_{\omega\in\mathbb{R}^{+}}\bigg{(}\psi_{\alpha}(\mathbf{x}_{d},\omega;t)\sigma^{*}(t)-\psi_{\alpha}^{*}(\mathbf{x}_{d},\omega;t)\sigma(t)\bigg{)}dt\
d\omega,$ (107b) where
$\displaystyle\psi_{\alpha}(\mathbf{x},\omega;t)=\bigg{(}\frac{\omega^{2}\alpha(\omega)}{\pi}\bigg{)}^{1/2}\int\text{G}_{\alpha}(\mathbf{x},\mathbf{x}^{\prime};\omega)\phi_{\alpha}(\mathbf{x}^{\prime},\omega;t)d^{3}\mathbf{x}^{\prime}.$
(107c)

We note that the integrand in Eq. 106 is gaussian in
$\phi_{\alpha}(\mathbf{x},\omega;t)$ and consequently the integral with
respect to $\phi_{\alpha}(\mathbf{x},\omega;t)$ can be evaluated exactly using
the stationary phase method. The stationary point of the integrand in Eq. 106
with respect to $\phi_{\alpha}(\mathbf{x},\omega;t)$ is given by the classical
path $\phi_{\text{cl}}(\mathbf{x},\omega;t)$ that minimizes the action
$\mathcal{S}_{\text{bath}}+\mathcal{S}_{\text{int}}$, which is governed by

$\displaystyle i\frac{\partial}{\partial
t}{\phi}_{\text{cl}}(\mathbf{x},\omega;t)=\omega\phi_{\text{cl}}(\mathbf{x},\omega;t)+iV_{0}\bigg{(}\frac{\omega^{2}\alpha(\omega)}{\pi}\bigg{)}^{1/2}\text{G}_{\alpha}^{*}(\mathbf{x},\mathbf{x}_{d};\omega)\sigma(t).$
(108)

Since we are interested in evaluating expectations in which the bath is in the
vacuum state, $\phi_{\text{cl}}(\mathbf{x},\omega;t_{i})=0$ — the above
equation of motion can then be integrated to obtain

$\displaystyle\phi_{\text{cl}}(\mathbf{x},\omega;t)=V_{0}\bigg{(}\frac{\omega^{2}\alpha(\omega)}{\pi}\bigg{)}^{1/2}\text{G}_{\alpha}^{*}(\mathbf{x},\mathbf{x}_{d};\omega)\int_{t_{i}}^{t}\sigma(t^{\prime})e^{-i\omega(t-t^{\prime})}dt^{\prime}.$
(109)

Next, we perform a change of variables in the path integral from
$\phi_{\alpha}(\mathbf{x},\omega;t)$ to
$\delta\phi_{\alpha}(\mathbf{x},\omega;t)=\phi_{\text{cl}}(\mathbf{x},\omega;t)+\delta\phi_{\alpha}(\mathbf{x},\omega;t)$.
Since $\mathcal{S}_{\text{bath}}+\mathcal{S}_{\text{int}}$ is quadratic in
$\phi_{\alpha}(\mathbf{x},\omega;t)$ and
$\phi_{\text{cl}}(\mathbf{x},\omega;t)$ is a stationary point of this action,
it is expressible as a quadratic in $\delta\phi_{\alpha}(\mathbf{x},\omega;t)$
without any linear terms in $\delta\phi_{\alpha}(\mathbf{x},\omega;t)$

$\displaystyle\mathcal{S}_{\text{bath}}+\mathcal{S}_{\text{int}}=\mathcal{S}_{0}+\delta\mathcal{S},$
(110)

where

$\displaystyle\delta\mathcal{S}=\int_{t=t_{i}}^{t_{f}}\int_{\mathbf{x}\in\mathbb{R}^{3}}\int_{\omega\in\mathbb{R}^{+}}\bigg{(}\delta\phi^{*}(\mathbf{x},\omega;t)\frac{\partial}{\partial
t}\delta\phi_{\alpha}(\mathbf{x},\omega;t)-\omega\big{|}\delta\phi_{\alpha}(\mathbf{x},\omega;t)\big{|}^{2}\bigg{)}dt\
d\mathbf{x}\ d\omega,$ (111)

and $\mathcal{S}_{0}$ is the action
$\mathcal{S}_{\text{bath}}+\mathcal{S}_{\text{int}}$ evaluated at
$\phi_{\text{cl}}(\mathbf{x},\omega;t)$

$\displaystyle\mathcal{S}_{0}$
$\displaystyle=-iV_{0}\int_{t=t_{i}}^{t_{f}}\int_{\mathbf{x}\in\mathbb{R}^{3}}\int_{\omega\in\mathbb{R}^{+}}\bigg{(}\frac{\omega^{2}\alpha(\omega)}{\pi}\bigg{)}^{1/2}\text{G}_{\alpha}(\mathbf{x}_{d},\mathbf{x};\omega)\phi_{\text{cl}}(\mathbf{x},\omega;t)\sigma^{*}(t)dt\
d^{3}\mathbf{x}\ d\omega,$
$\displaystyle=\frac{iV_{0}^{2}}{\pi}\int_{t=t_{i}}^{t_{f}}\int_{t^{\prime}=t_{i}}^{t}\int_{\omega\in\mathbb{R}^{+}}\text{Im}\big{[}\text{G}(\mathbf{x}_{d},\mathbf{x}_{d};\omega)\big{]}\sigma^{*}(t)\sigma(t^{\prime})e^{-i\omega(t-t^{\prime})}dt\
dt^{\prime}\ d\omega,$ (112)

wherein in the second step we have used the fact that $\lim_{\alpha(\omega)\to
0}\text{Im}[\text{G}_{\alpha}(\mathbf{x}_{d},\mathbf{x}_{d};\omega)]$ exists
and have explicitly taken this limit. The expectation
$\mathcal{E}(t_{1},t_{2}\dots t_{N})$ then evaluates to

$\displaystyle\mathcal{E}(t_{1},t_{2}\dots t_{M})=\mathcal{A}\int\mathcal{D}s\
\bigg{[}\prod_{i=1}^{M}c_{i}(t_{i})\bigg{]}e^{-i\mathcal{S}_{\text{eff}}}$
(113)

where

$\displaystyle\mathcal{S}_{\text{eff}}=\mathcal{S}_{\text{sys}}+\frac{iV_{0}^{2}}{\pi}\int_{t=t_{i}}^{t_{f}}\int_{t^{\prime}=t_{i}}^{t}\int_{\omega\in\mathbb{R}^{+}}\text{Im}\big{[}\text{G}(\mathbf{x}_{d},\mathbf{x}_{d};\omega)\big{]}\sigma^{*}(t)\sigma(t^{\prime})e^{-i\omega(t-t^{\prime})}dt\
dt^{\prime}\ d\omega,$ (114)

and

$\displaystyle\mathcal{A}=\int_{\begin{subarray}{c}\delta\phi_{\alpha}(\mathbf{x},\omega;t_{i})=0\\\
\delta\phi_{\alpha}(\mathbf{x},\omega;t_{f})=0\end{subarray}}\mathcal{D}[\delta\phi_{\alpha}(\mathbf{x},\omega;t),\delta\phi^{*}(\mathbf{x},\omega;t)]e^{-i\delta\mathcal{S}}.$
(115)

We immediately notice that $\mathcal{A}$ is just the expectation
$\bra{\text{vac}}e^{-iH_{\text{bath}}(t_{f}-t_{i})}\ket{\text{vac}}$ and thus
evaluates to 1. This yields

$\displaystyle\mathcal{E}(t_{1},t_{2}\dots t_{M})=\int\mathcal{D}s\
\bigg{[}\prod_{i=1}^{M}c_{i}(t_{i})\bigg{]}e^{-i\mathcal{S}_{\text{eff}}}.$
(116)

We note that the effective action $\mathcal{S}_{\text{eff}}$ is, in general, a
nonlocal action consequently making the integral in Eq. 116 difficult to
evaluate directly. However, within the Markovian approximation
carmichael2009open , this effective action becomes local. This immediately
follows from Eq. 114 by noting that the Markovian approximation around a
(system resonance) frequency $\omega_{0}$ approximates
$\text{Im}[\text{G}(\mathbf{x}_{d},\mathbf{x}_{d};\omega)]$ by
$\text{Im}[\text{G}(\mathbf{x}_{d},\mathbf{x}_{d};\omega_{0})]$ and extends
the integral with respect to $\omega$ from $-\infty$ to $\infty$ to obtain:

$\displaystyle\mathcal{S}_{\text{eff}}=\mathcal{S}_{\text{sys}}+\frac{i\gamma}{2}\int_{t=t_{i}}^{t_{f}}|\sigma(t)|^{2}dt$
(117)

where
$\gamma=V_{0}^{2}\text{Im}[\text{G}(\mathbf{x}_{d},\mathbf{x}_{d};\omega_{0})]$.
From this local lagrangian, we immediately obtain that the expectation in Eq.
116 can be evaluated by evolving the localized system with respect to a
dissipative (non-Hermitian) Hamiltonian $H_{\text{eff}}$ defined by:

$\displaystyle
H_{\text{eff}}=H_{\text{sys}}-\frac{i\gamma}{2}\sigma^{\dagger}\sigma$ (118)