Point-coupling Hamiltonian for frequency-independent linear optical devices

TL;DR

This paper introduces the point-coupling Hamiltonian as a versatile model for frequency-independent linear optical devices, enabling formal analysis and simulation of complex optical systems with feedback.

Contribution

It provides a method to construct the Hamiltonian from classical scattering matrices and demonstrates its application in simulating feedback systems using matrix-product states.

Findings

Derived the quantum scattering matrix from the Hamiltonian.

Showed equivalence between quantum scattering and classical inverse scattering.

Applied the model to simulate time-delayed feedback systems.

Abstract

We present the point-coupling Hamiltonian as a model for frequency-independent linear optical devices acting on propagating optical modes described as a continua of harmonic oscillators. We formally integrate the Heisenberg equations of motion for this Hamiltonian, calculate its quantum scattering matrix, and show that an application of the quantum scattering matrix on an input state is equivalent to applying the inverse of classical scattering matrix on the annihilation operators describing the optical modes. We show how to construct the point-coupling Hamiltonian corresponding to a general linear optical device described by a classical scattering matrix, and provide examples of Hamiltonians for some commonly used linear optical devices. Finally, in order to demonstrate the practical utility of the point-coupling Hamiltonian, we use it to rigorously formulate a matrix-product-state based simulation for time-delayed feedback systems wherein the feedback is provided by a linear optical device described by a scattering matrix as opposed to a hard boundary condition (e.g. a mirror with less than unity reflectivity).

Equations

$$H=\sum _{n=1}^N{\int }_{-\infty }^{\infty }\omega a_n^{†}(\omega )a_n(\omega )\text{d}\omega +\sum _{n=1}^N\sum _{m=1}^N{\text{V}}_{n,m}{a}_n^{†}(x_n^0)a_m(x_m^0),$$

$$a_n(x_n)={\int }_{-\infty }^{\infty }{a}_n(\omega )e^{\text{i}\omega x_n}\frac{\text{d}\omega }{\sqrt{2\pi }}.$$

$$\text{i}\frac{\text{d}{a}_n(\omega ;t)}{\text{d}t}=\omega a_n(\omega ;t)+\sum _{m=1}^N\frac{{\text{V}}_{n,m}}{\sqrt{2\pi }}{a}_m(x_m^0;t)e^{-\text{i}\omega x_n^0}\forall n\in \{1,2\dots N\}.$$

$$\displaystyle\begin{bmatrix}a_{1}(x_{1}^{0}+x;t_{0}+\tau)\\ a_{2}(x_{2}^{0}+x;t_{0}+\tau)\\ \vdots\\ a_{N}(x_{N}^{0}+x;t_{0}+\tau)\end{bmatrix}=\bigg{[}\textbf{I}-\textrm{i}\textbf{V}\bigg{(}\textbf{I}+\frac{\textrm{i}\textbf{V}}{2}\bigg{)}^{-1}\Theta(0\leq y\leq\tau)\bigg{]}\begin{bmatrix}a_{1}(x_{1}^{0}+x-\tau;t_{0})\\ a_{2}(x_{2}^{0}+x-\tau;t_{0})\\ \vdots\\ a_{N}(x_{N}^{0}+x-\tau;t_{0})\end{bmatrix},$$

$$\Theta (x_1\le x\le x_2)=\{\begin{array}{ll}1& \text{if }x\in (x_1,x_2)\\ 0& \text{if }x\in (-\infty ,x_1)\cup (x_2,\infty )\\ \frac{1}{2}& \text{if }x\in \{x_1,x_2\}\end{array}.$$

$$\left[\begin{array}{c}{a}_1(x_1^0+x;t_0+\tau )\\ a_2(x_2^0+x;t_0+\tau )\\ ⋮\\ a_N(x_N^0+x;t_0+\tau )\end{array}\right]=\text{S}\left[\begin{array}{c}{a}_1(x_1^0+x-\tau ;t_0)\\ a_2(x_2^0+x-\tau ;t_0)\\ ⋮\\ a_N(x_N^0+x-\tau ;t_0)\end{array}\right],$$

$$\displaystyle\textbf{S}=\textbf{I}-\textrm{i}\textbf{V}\bigg{(}\textbf{I}+\frac{\textrm{i}\textbf{V}}{2}\bigg{)}^{-1}=\bigg{(}\textbf{I}-\frac{\textrm{i}\textbf{V}}{2}\bigg{)}\bigg{(}\textbf{I}+\frac{\textrm{i}\textbf{V}}{2}\bigg{)}^{-1}.$$

$$\displaystyle\textbf{V}=-\textbf{U}\ \text{diag}\bigg{[}2\tan\bigg{(}\frac{\bm{\phi}}{2}\bigg{)}\bigg{]}\ \textbf{U}^{\dagger}.$$

$$\displaystyle H_{\text{phase-shifter}}=\int_{-\infty}^{\infty}\omega a^{\dagger}(\omega)a(\omega)\textrm{d}\omega-2\tan\bigg{(}\frac{\varphi}{2}\bigg{)}a^{\dagger}(x=0)a(x=0),$$

$$\text{S}=\left[\begin{array}{cc}\cos \theta & \sin \theta e^{\text{i}\phi }\\ -\sin \theta e^{-\text{i}\phi }& \cos \theta \end{array}\right].$$

$$\text{V}=\left[\begin{array}{cc}0& 2\text{i}\tan (\theta /2)e^{\text{i}\phi }\\ -2\text{i}\tan (\theta /2)e^{-\text{i}\phi }& 0\end{array}\right],$$

$$\displaystyle H_{\text{beam-splitter}}=\sum_{k\in\{1,2\}}\int_{-\infty}^{\infty}\omega a_{k}^{\dagger}(\omega)a_{k}(\omega)\textrm{d}\omega+\bigg{[}2\textrm{i}\tan\bigg{(}\frac{\theta}{2}\bigg{)}e^{\textrm{i}\varphi}a_{1}^{\dagger}(x_{1}=0)a_{2}(x_{2}=0)+\textrm{h.c.}\bigg{]},$$

$$\text{S}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right],$$

$$\text{V}=\left[\begin{array}{ccc}0& 2\text{i}& -2\text{i}\\ -2\text{i}& 0& 2\text{i}\\ 2\text{i}& -2\text{i}& 0\end{array}\right],$$

$$\displaystyle H_{\text{circulator}}=\sum_{k\in\{1,2,3\}}\int_{-\infty}^{\infty}\omega a_{k}^{\dagger}(\omega)a_{k}(\omega)\textrm{d}\omega+\bigg{[}\sum_{k\in\{1,2,3\}}2\textrm{i}a_{k}^{\dagger}(x_{k}=0)a_{k+1}(x_{k+1}=0)+\textrm{h.c.}\bigg{]},$$

$$\mathcal{S}=\underset{\begin{array}{c}{t}_{+}\to \infty \\ t_{-}\to -\infty \end{array}}{\lim }{e}^{\text{i}{H}_0{t}_{+}}{e}^{-\text{i}H(t_{+}-t_{-})}{e}^{-\text{i}{H}_0{t}_{-}},$$

$$H_0=\sum _{k=1}^N{\int }_{-\infty }^{\infty }\omega a_k^{†}(\omega )a_k(\omega )\text{d}\omega .$$

$$\displaystyle\mathcal{S}(\textbf{x},\bm{\mu};\textbf{x}^{\prime},\bm{\mu}^{\prime})=\bra{\text{vac}}\bigg{[}\prod_{i=1}^{K}a_{\mu_{i}}(x_{i})\bigg{]}\mathcal{S}\bigg{[}\prod_{i=1}^{K}a_{\mu_{i}^{\prime}}^{\dagger}(x_{i}^{\prime})\bigg{]}\ket{\text{vac}}.$$

$$\displaystyle\mathcal{S}(\textbf{x},\bm{\mu};\textbf{x}^{\prime},\bm{\mu}^{\prime})=\lim_{\begin{subarray}{c}t_{+}\to\infty\\ t_{-}\to-\infty\end{subarray}}\bra{\text{vac}}\bigg{[}\prod_{i=1}^{K}a_{\mu_{i}}(x_{i}+t_{+};t_{+})\bigg{]}\bigg{[}\prod_{i=1}^{K}a_{\mu_{i}^{\prime}}^{\dagger}(x_{i}^{\prime}+t_{-};t_{-})\bigg{]}\ket{\text{vac}}.$$

$$a_{{\mu }_i}(x_i+t_{+};t_{+})=\sum _{{\nu }_i=1}^N{\text{S}}_{{\mu }_i,{\nu }_i}{a}_{{\nu }_i}(x_i-x_{{\mu }_i}^0+x_{{\nu }_i}^0+t_{-};t_{-}),$$

$$\displaystyle\mathcal{S}(\textbf{x},\bm{\mu};\textbf{x}^{\prime},\bm{\mu}^{\prime})=\sum_{\mathcal{P}_{K}}\prod_{l=1}^{K}\text{S}_{\mu_{l},\mathcal{P}_{K}\mu^{\prime}_{l}}\delta\big{(}(x_{l}-x_{\mu_{l}}^{0})-(\mathcal{P}_{K}x_{l}^{\prime}-x^{0}_{\mathcal{P}_{K}\mu^{\prime}_{l}})\big{)},$$

$$\displaystyle\mathcal{S}(\bm{\omega},\bm{\mu};\bm{\omega}^{\prime},\bm{\mu}^{\prime})=\int_{\mathbb{R}^{K}}\int_{\mathbb{R}^{K}}S(\textbf{x},\bm{\mu};\textbf{x}^{\prime},\bm{\mu}^{\prime})\bigg{[}\prod_{l=1}^{K}\frac{e^{\textrm{i}(\omega_{l}^{\prime}x_{l}^{\prime}-\omega_{l}x_{l})}}{2\pi}\bigg{]}\textrm{d}^{K}\textbf{x}\ \textrm{d}^{K}\textbf{x}^{\prime}=\sum_{\mathcal{P}_{K}}\prod_{l=1}^{K}\textrm{S}_{\mathcal{P}_{K}\mu_{l},\mu_{l}^{\prime}}\ e^{-\textrm{i}\omega_{l}\big{(}x^{0}_{\mathcal{P}_{K}\mu_{l}}-x^{0}_{\mu_{l}^{\prime}}\big{)}}\delta(\mathcal{P}_{K}\omega_{l}-\omega_{l}^{\prime}).$$

$$\displaystyle\ket{\psi_{\text{in}}}=\sum_{\bm{\mu}}\int_{\mathbb{R}^{K}}\psi_{\text{in}}(\bm{\omega},\bm{\mu})\bigg{[}\prod_{l=1}^{K}a_{\mu_{l}}^{\dagger}(\omega_{l})\bigg{]}\ket{\text{vac}}\textrm{d}^{K}\bm{\omega},$$

$$\displaystyle\ket{\psi_{\text{out}}}=\sum_{\bm{\mu}}\int_{\mathbb{R}^{K}}\psi_{\text{out}}(\bm{\omega},\bm{\mu})\bigg{[}\prod_{l=1}^{K}a_{\mu_{l}}^{\dagger}(\omega_{l})\bigg{]}\ket{\text{vac}}\textrm{d}^{K}\bm{\omega},$$

$${\psi }_{\text{out}}(\mathbf{\omega },\mathbf{\mu })=\frac{1}{N!}\sum _{{\mathbf{\mu }}^{\prime }}{\int }_{{\mathbb{R}}^K}S(\mathbf{\omega },\mathbf{\mu };{\mathbf{\omega }}^{\prime },{\mathbf{\mu }}^{\prime }){\psi }_{\text{in}}({\mathbf{\omega }}^{\prime },{\mathbf{\mu }}^{\prime }){\text{d}}^K{\mathbf{\omega }}^{\prime }.$$

$${\psi }_{\text{out}}(\mathbf{\omega },\mathbf{\mu })$$

$$|{\psi }_{\text{out}}⟩$$

$$[b_n(\omega ),H]=\omega b_n(\omega ),$$

$$[b_n(\omega ),b_m^{†}({\omega }^{\prime })]={\delta }_{n,m}\delta (\omega -{\omega }^{\prime }),$$

$$H=\sum _{n=1}^N{\int }_{-\infty }^{\infty }\omega b_n^{†}(\omega )b_n(\omega )\text{d}\omega ,$$

$$b_n(\omega )={\int }_{-\infty }^{x_n^0}{e}^{-\text{i}\omega (x-x_n^0)}{a}_n(x_n)\frac{\text{d}{x}_n}{\sqrt{2\pi }}+\sum _{m=1}^N{\text{S}}_{m,n}^{\ast }{\int }_{x_m^0}^{\infty }{e}^{-\text{i}\omega (x-x_m^0)}{a}_m(x_m)\frac{\text{d}{x}_m}{\sqrt{2\pi }},$$

Full text

††preprint: APS/123-QED

I Introduction

Linear-optical elements are key components in any quantum information processing system Blais et al. (2007); Imamoglu et al. (1999); Duan et al. (2001). They can be used for applying phase shifts, interfering, and splitting quantum light propagating in waveguide modes or collimated optical beams. Consequently, they are traditionally described by scattering matrices Fan et al. (2003); Ruan and Fan (2009). The response of the linear optical element to an incoming state is often calculated by applying the inverse of the classical scattering matrix on the annihilation operators Mandel and Wolf (1995); Kok et al. (2007); Skaar et al. (2004). The full time evolution of a quantum state incident on linear-optical elements has also been captured within the framework of quantum stochastic differential equations Gough and James (2009); Petersen (2010); Parthasarathy (2012). Several attempts have been made to calculate a Hamiltonian that captures the dynamics of linear optical elements Prasad et al. (1987); Fearn and Loudon (1987); Garcia-Escartin et al. (2016); Leonhardt and Neumaier (2003). In particular, Refs. Leonhardt and Neumaier (2003) and Garcia-Escartin et al. (2016) show the existence and explicit construction of an anti-hermitian matrix whose exponential reproduces the classical scattering matrix of a linear-optical device. A Hamiltonian describing the quantum physics of the linear-optical device can then be constructed from this matrix. However, this approach needs to introduce a fictitious time corresponding to the duration for which the incident photons interact with the linear-optical device. While the introduction of this fictitious time does not impact the relationship between the input and output states, it makes the model unphysical if the full dynamics of the quantum state is desired. Alternatively, as shown in Ref. Dalton et al. (1999), a physically accurate Hamiltonian for the linear-optical device can be written down if the normal modes of the linear-optical device are known.

In this paper, we propose a Hamiltonian that can model an arbitrary frequency-independent linear optical device acting on propagating optical modes. The Hamiltonian assumes that the optical modes couple to each other at a single-point in space — we therefore call it the point-coupling Hamiltonian. For a given frequency-independent classical scattering matrix implemented by the linear optical device, we provide a recipe to construct a point-coupling Hamiltonian describing the device. We formally integrate the Heisenberg equations of motion for the point-coupling Hamiltonian, and use the resulting solution to calculate its quantum scattering matrix Fan et al. (2010); Xu and Fan (2015); Trivedi et al. (2018). It is shown that an application of the quantum scattering matrix on an incoming quantum state is equivalent to applying the inverse of the classical scattering matrix on the annihilation operators of the optical modes in the incoming quantum state, thereby reproducing the commonly used procedure for analyzing the quantum physics of linear-optical devices. We also ‘diagonalize’ the point-coupling Hamiltonian to calculate its normal modes — this provides a connection between the Hamiltonian presented in this paper and the quantization approach described in Ref. Dalton et al. (1999)

Finally, in order to demonstrate the utility of the point-coupling Hamiltonian proposed in this paper, we use it to rigorously derive and implement a matrix-product-state (MPS) Schollwöck (2011); Orús (2014) based update to simulate time-delayed feedback systems with linear optical devices providing feedback Pichler and Zoller (2016); Grimsmo (2015). Time-delayed feedback systems typically have a low-dimensional quantum system, such as a two-level system, coupling to a well-defined optical mode and the emission from the quantum system is used to re-excite the quantum system via a feedback path. The feedback path can be constructed using a linear-optical device such as a mirror, or using another quantum system. Time-delayed feedback systems possibly provide a platform for generation of highly entangled states for quantum computation Pichler et al. (2017); Lubasch et al. (2018), quantum simulation Dhand et al. (2018) as well as implementation of quantum memories using long-lived bound states Calajó et al. (2019). MPS based simulations of such systems were proposed in Ref. Pichler and Zoller (2016), although the formulation relied on hard boundary conditions e.g. it assumed that the feedback to the quantum system was provided by an ideal mirror. Using the point-coupling Hamiltonian, we outline a method to extend the MPS based simulation of this system to account for the situation wherein the mirror is described by a scattering matrix, and thus enable an analysis of the impact of non-ideality in the mirror on the dynamics of the feedback-system. We expect the point-coupling Hamiltonian proposed in this paper to be of utility in analyzing quantum systems with complicated linear optical devices that can only be described by a full scattering matrix, and cannot be well-approximated by hard boundary conditions.

This paper is organized as follows — Section II introduces the point-coupling Hamiltonian and analyzes it in the Heisenberg picture. Using the solution of the Heisenberg equations of motion, we explicitly calculate the quantum scattering matrix for the linear-optical device. It is shown that an application of the quantum scattering matrix on an incoming state is equivalent to applying the inverse of the classical scattering matrix on the annihilation operators in the state. We also provide the point-coupling Hamiltonians corresponding to some commonly used linear optical devices such as phase shifters, beam splitters and optical circulators. Finally, we diagonalize the point-coupling Hamiltonian and show that the classical modes corresponding to the annihilation operators that diagonalize the Hamiltonian can be interpreted as the normal modes of the linear-optical device. In section III, we use the proposed point-coupling Hamiltonian to analyze feedback into a two-level system from a partially transmitting mirror using MPS update.

II Point-coupling Hamiltonian

II.1 Dynamics in the Heisenberg picture

Consider $N$ propagating optical modes (which can physically be waveguide modes, or collimated optical beams) interacting with each other through a linear-optical device (Fig. 1). Note that optical modes that are identical but for the direction of propagation are counted as separate modes. Labelling by $a_{n}(\omega)$ the annihilation operator of the $n^{\text{th}}$ optical mode at frequency $\omega$, we propose the following Hamiltonian for describing the dynamics of the system:

[TABLE]

where $a_{n}(x_{n})$ is the position-domain annihilation operator for the $n^{\text{th}}$ optical mode at coordinate $x_{n}$ along its direction of propagation:

[TABLE]

Here the position $x_{n}$ is expressed in units of time such that the group velocity, assumed to be frequency independent and uniform across all the optical modes, is unity. The coordinate of the linear-optical device along the direction of propagation of the $n^{\text{th}}$ optical mode is $x_{n}^{0}$. Note that we allow for the coordinate $x_{n}$ for different optical modes to be expressed in different coordinate systems — in general, each optical mode can have its own coordinate system with the $x$ axis of the coordinate system being along its direction of propagation. For e.g. the $x$ axis for a forward propagating mode will be in a direction opposite to the $x$ axis for a backward propagating mode, or modes waveguides oriented in physically different directions will have differently oriented coordinate systems. We note that $a_{n}(\omega)$, and consequently $a_{n}(x_{n})$, satisfy the usual bosonic commutation relations: $[a_{n}(\omega),a_{m}^{\dagger}(\omega^{\prime})]=\delta_{n,m}\delta(\omega-\omega^{\prime})$, $[a_{n}(\omega),a_{m}(\omega^{\prime})]=0$, $[a_{n}(x_{n}),a_{m}^{\dagger}(x^{\prime}_{m})]=\delta_{n,m}\delta(x_{n}-x^{\prime}_{n})$ and $[a_{n}(x_{n}),a_{m}(x^{\prime}_{m})]=0$. Moreover, as a consequence of the Hermiticity of $H$, the coupling coefficients $\text{V}_{k,m}$ satisfy $\text{V}_{n,m}=\text{V}_{m,n}^{*}$.

 

The dynamics of this Hamiltonian can be easily analyzed using the Heisenberg’s equations of motion, which are given by:

[TABLE]

As is shown in appendix A, these equations can easily be integrated to relate the position-domain annihilation operators at time $t=t_{0}+\tau$ and displaced by a distance $x\in(-\infty,\infty)$ from the linear optical device to the position-domain annihilation operators at time $t=t_{0}$:

[TABLE]

where V is a $N\times N$ Hermitian matrix formed by $\text{V}_{m,n}$ as its elements, I is the identity matrix of size $N$ and $\Theta(x_{1}\leq x\leq x_{2})$ is defined by:

[TABLE]

To intuitively interpret the result in Eq. 4, note that if $y<0$, then the position-domain annihilation operator at $t=t_{0}+\tau$ is simply a propagated version of itself at $t=t_{0}$. This is a simply a consequence of the physical points described by $x<0$ lying before the linear optical device along the propagation direction, and consequently being unaffected by scattering from the linear-optical device. This situtation is the same for $x>\tau$, since the scattered light from the linear-optical device has not had sufficient time to propagate to the points in question from the location of the linear optical device (i.e. from $x=0$). For $0<x<\tau$, the optical mode annihilation operator at $t=t_{0}+\tau$ is a sum of a propagated version of itself and contributions from other optical modes as scattered by the linear optical device, and Eq. 4 can be simplified to:

[TABLE]

where

[TABLE]

Since S relates the optical mode annihilation operators before and after scattering from the linear optical-device has occurred, it can be interpreted as the classical scattering matrix corresponding to the linear optical-device. It can readily be verified that a consequence of V being Hermitian is that the matrix S is unitary. Moreover, if S has the diagonalization \textbf{S}=\textbf{U}\ \text{diag}\big{[}e^{\textrm{i}\bm{\phi}}\big{]}\ \textbf{U}^{\dagger}, it follows from Eq. 7 that:

[TABLE]

Therefore, given the classical matrix S of a linear-optical device, Eq. 8 allows us to construct the matrix V and by extension the point-coupling Hamiltonian that can model the device. As examples, we consider some commonly used linear optical devices (Fig. 2) and construct the point-coupling Hamiltonians that describe their dynamics:

1. (a)

Phase shifter: The classical scattering matrix S of the phase-shifter [Fig. 2(a)] is a single-element matrix: $\textbf{S}=[e^{\textrm{i}\varphi}]$. From Eq. 8, we obtain $\textbf{V}=[-2\tan(\varphi/2)]$. The Hamiltonian for the phase shifter can thus be written as:

[TABLE]

where $a(\omega)$ [$a(x)$] is the frequency-domain (position-domain) annihilation operator for the optical mode that the phase-shifter is acting on. 2. (b)

Beam-splitter: The beam-splitter [Fig. 2(b)] is described by the following classical scattering matrix:

[TABLE]

Again, using Eq. 8, we obtain

[TABLE]

from which we can construct the beam splitter Hamiltonian:

[TABLE]

where $a_{1,2}(\omega)$ [$a_{1,2}(x_{1,2})$] are the frequency-domain (position-domain) annihilation operators for the optical modes that the beam-splitter is acting on. 3. (c)

Optical circulator: The optical circulator [Fig. 2(c)] routes an excitation in optical mode 1 to optical mode 2, optical mode 2 to optical mode 3 and optical mode 3 to optical mode 1. It can thus be described by the following $3\times 3$ classical scattering matrix:

[TABLE]

and therefore

[TABLE]

from which we can construct the optical circulator Hamiltonian:

[TABLE]

where $a_{1,2,3}(\omega)$ [$a_{1,2,3}(x_{1,2,3})$] are the frequency-domain (position-domain) annihilation operators for the optical modes that the optical circulator acts on, and $a_{4}(\omega)/a_{4}(x_{4})$ is to be interpreted as $a_{1}(\omega)/a_{1}(x_{1})$.

II.2 Quantum scattering matrix of the point-coupling Hamiltonian

While the previous subsection analyzed the point-coupling Hamiltonian in the Heisenberg picture, in this subsection we analyze its dynamics in the Schroedinger picture. In particular, we consider the problem of exciting the linear-optical device with an input state, and attempt to calculate the output state produced by the device.

 

The key object relating the input state (assumed to be the asymptote Taylor (2006) to the system at $t\to-\infty$) of the system to its output state (assumed to be the asymptote Taylor (2006) to the state of the system at $t\to\infty$) is the quantum Taylor (2006):

[TABLE]

where $H_{0}$ is the Hamiltonian of the optical modes without accounting for the linear-optical device:

[TABLE]

Consider now the computation of the following $K$ photon matrix element of the scattering matrix ($\textbf{x}\equiv\{x_{1},x_{2}\dots x_{K}\}$, $\textbf{x}^{\prime}\equiv\{x_{1}^{\prime},x_{2}^{\prime}\dots x_{K}^{\prime}\}$, $\bm{\mu}=\{\mu_{1},\mu_{2}\dots\mu_{K}\}$ and $\bm{\mu}^{\prime}=\{\mu_{1}^{\prime},\mu_{2}^{\prime}\dots\mu_{K}^{\prime}\}$):

[TABLE]

Note that the relations $\exp(-\textrm{i}H_{0}t)a_{k}(x_{k})\exp(\textrm{i}H_{0}t)=a_{k}(x_{k}+t)$ and $\exp(\textrm{i}Ht)a_{k}(x_{k})\exp(-\textrm{i}Ht)=a_{k}(x_{k};t)$ together with $H_{0}\ket{\text{vac}}=H\ket{\text{vac}}=0$ immediately imply the following relationship between the scattering matrix element in Eq. 18 and the Heisenberg picture optical mode position-domain annihilation operators:

[TABLE]

Using Eq. 4, and in the limit of $t_{+}\to\infty$ and $t_{-}\to-\infty$, it follows that :

[TABLE]

where $S_{i,j}$ are the elements of the classical scattering matrix S. With this, the following explicit expression for $\mathcal{S}(\textbf{x},\bm{\mu};\textbf{x}^{\prime},\bm{\mu}^{\prime})$ as given by Eq. 19 can be obtained:

[TABLE]

where $\mathcal{P}_{K}$ is a $K$ element permutation. It can immediately be noticed that the quantum scattering matrix elements are completely determined in terms of the classical scattering matrix elements $\text{S}_{i,j}$. From Eq. 21, we can also evaluate the frequency domain scattering matrix elements $S(\bm{\omega},\bm{\mu};\bm{\omega}^{\prime},\bm{\mu}^{\prime})$ ($\bm{\omega}=\{\omega_{1},\omega_{2}\dots\omega_{K}\}$, $\bm{\omega}^{\prime}=\{\omega_{1}^{\prime},\omega_{2}^{\prime}\dots\omega_{K}^{\prime}\}$, $\bm{\mu}=\{\mu_{1},\mu_{2}\dots\mu_{K}\}$ and $\bm{\mu}^{\prime}=\{\mu_{1}^{\prime},\mu_{2}^{\prime}\dots\mu_{K}^{\prime}\}$):

[TABLE]

Note that the frequency domain scattering matrix doesn’t have any connected parts Xu and Fan (2015) i.e. scattering of a $K$ photon wave-packet from the linear optical device conserves the individual input frequencies. This is a direct consequence of the ‘linearity’ of the optical device. To gain more insight into the form of the scattering matrix in Eqs. 21 and 22, consider the calculation of the output state corresponding to a $K$-photon input state $\ket{\psi_{\text{in}}}$:

[TABLE]

where the amplitude $\psi_{\text{in}}(\bm{\omega},\bm{\mu})$ can be chosen, without any loss of generality, to be symmetric with respect to a simultaneous permutation of the indices $\bm{\mu}$ and $\bm{\omega}$: $\psi_{\text{in}}(\bm{\omega},\bm{\mu})=\psi_{\text{in}}(\mathcal{P}_{K}\bm{\omega},\mathcal{P}_{K}\bm{\mu})\ \forall\ K$-element permutations $\mathcal{P}_{K}$. The output state is then given by:

[TABLE]

where $\psi_{\text{out}}(\bm{\omega},\bm{\mu})$ is given by:

[TABLE]

Assuming that the coordinate systems of the optical modes are chosen such that $x_{k}^{0}=0\ \forall\ k\in\{1,2\dots N\}$, from Eqs. 22 and 25 we obtain:

[TABLE]

Therefore,

[TABLE]

where $\tilde{a}_{i}(\omega)=\sum_{i=1}^{N}\text{S}_{j,i}^{*}a_{j}(\omega)$. The application of the quantum scattering matrix of a linear optical device to an input state is equivalent to replacing the annihilation operators in the input state with $\textbf{S}^{\dagger}$ times the annihilation operators. This is the usual procedure used to analyze the impact of a linear optical element on an input state Mandel and Wolf (1995); Kok et al. (2007), and our analysis derives it from a Hamiltonian based description of the linear optical element and thus lends rigor to this procedure.

II.3 Normal modes of the point-coupling Hamiltonian

In this subsection, we diagonalize the point-coupling Hamiltonian and explicitly calculate the normal modes of the linear-optical device. Given that the frequency-domain annihilation operator of the $n^{\text{th}}$ normal mode is $b_{n}(\omega)$, it should satisfy the following commutation relations:

[TABLE]

where $H$ is the point-coupling Hamiltonian. Eq. 28a is a consequence of the Hamiltonian $H$ being expressible as a sum of independent continua of harmonic oscillators in the normal mode basis:

[TABLE]

and Eq. 28b is a result of the different normal modes being physically independent modes. As is shown in appendix B, using Eqs. 28a and 28b, the normal mode annihilation operators $b_{n}(\omega)$ can be chosen to be:

[TABLE]

where $\text{S}_{n,m}$ are the elements of the classical scattering matrix. To obtain a physical interpretation of this result, consider creating a photon in the normal mode described by $b_{n}(\omega)$ and calculating its projection on the modes described by $a_{m}(x_{m})$ — this projection is given by the expectation $\bra{\text{vac}}a_{m}(x_{m})b_{n}^{\dagger}(\omega)\ket{\text{vac}}$ which can be readily evaluated using Eq. 30:

[TABLE]

This is exactly the same field profile that would be obtained on classically exciting the linear optical device through the $n^{\text{th}}$ input port at frequency $\omega$, and calculating the fields scattered in all the output ports — the $n^{\text{th}}$ normal mode is simply the continuum of quantum harmonic oscillators associated with this field profile.

 

Finally, Eq. 30 can be inverted to relate the annihilation operators $a_{n}(x)$ to the normal mode annihilation operators $b_{n}(\omega)$ (refer to appendix B for details):

[TABLE]

where, again, we see that the modes described by $a_{n}(x)$ are identical to the normal modes at points before the linear-optical device (i.e. $x_{n}<x_{n}^{0}$ in Eq. 32) and their linear combination at points after the linear-optical device. Eqs. 30 and 32 thus provides a connection between the point-coupling Hamiltonian for linear-optical device and a direct quantization of the normal modes of the device Dalton et al. (1999).

III Matrix-product-state based simulations of time-delayed feedback systems

In this section, we use the point-coupling Hamiltonian proposed in the section II to develop an MPS update for a two-level system coupled to a waveguide with a partially transmitting mirror. Using the formulated MPS update, we study the impact of the less than unity reflectivity of the mirror on the dynamics of the time-delayed feedback system.

The system under consideration is a two-level system with a time-delayed feedback shown in Fig. 3. The Hamiltonian for this system can be expressed as a sum of a two-level system Hamiltonian, waveguide Hamiltonian and the mirror Hamiltonian (which is modeled as a point interaction between the forward and backward propagating modes):

[TABLE]

where $H_{\text{TLS}}(t)$ is the Hamiltonian of the two-level system including a coherent drive, $H_{\text{wg}}$ is the Hamiltonian describing the forward and backward propagating waveguide modes, $H_{\text{mirror}}$ is the Hamiltonian of the mirror providing feedback to the two-level system and $H_{\text{wg-TLS}}$ is the interaction Hamiltonian between the two-level system and the two waveguide modes:

[TABLE]

Here $\sigma$ is the de-excitation operator for a two-level system with resonant frequency $\omega_{e}$, $a_{+}(\omega)$ and $a_{-}(\omega)$ are the frequency-domain annihilation operators for the forward and backward propagating waveguide modes respectively and $a_{+}(x_{+})$ and $a_{-}(x_{-})$ are the position-domain annihilation operators for the forward and backward propagating waveguide modes respectively. Note that the $x$-axis for the coordinate systems for the forward and backward propagating waveguide modes are chosen to be in their direction of propagation with the same origin. A mirror, modeled with a beam splitter Hamiltonian between the forward and backward propagating waveguide modes with parameters $(\theta,\varphi)$, is located at $x_{\pm}=0$. The two-level system interacts with the waveguide (with a decay rate $\gamma_{+}$ into the forward propagating waveguide mode and $\gamma_{-}$ into the backward propagating waveguide mode) at a distance of $t_{d}$ from the mirror — this corresponds to the point with coordinates $x_{+}=t_{d}$ and $x_{-}=-t_{d}$ in the coordinate systems of the forward and backward propagating waveguide modes. Furthermore, the two-level system is driven with a laser pulse at frequency $\omega_{0}$ and pulse-shape $\Omega(t)$.

 

We first go into a frame rotating as per the Hamiltonian $H_{0}=\omega_{0}\sigma^{\dagger}\sigma+H_{\text{wg}}+H_{\text{mirror}}$. In this frame, the state of the system evolves as per the Hamiltonian $\tilde{H}(t)$:

[TABLE]

where $\delta_{e}=\omega_{e}-\omega_{0}$ and $a_{\pm}(x_{\pm}=\pm t_{d};t)$ are the Heisenberg picture operators corresponding to $a_{\pm}(x_{\pm}=\pm t_{d})$ with respect to $H_{0}$ at time $t$, subject to the initial condition $a_{\pm}(x_{\pm}=\pm t_{d};t=0)=a_{\pm}(x_{\pm}=\pm t_{d})$. From Eq. 4 and for $t>0$, we obtain:

[TABLE]

where we have defined the operators $A_{\pm}(t)$ via:

[TABLE]

An MPS update for the state of the system in the rotating frame with respect to $H_{0}$ can now be framed using the procedure outlined in Ref. Pichler and Zoller (2016) — the first step is to discretize the waveguide Hilbert space. Using a discretization step $\Delta t$, we can define the waveguide bin operators $A_{\pm}[k]$ in terms of the operators $A_{\pm}(t)$ via:

[TABLE]

which satisfy the commutation relations $[A_{+}[i],A_{+}^{\dagger}[j]]=\delta_{i,j}$ and $[A_{-}[i],A_{-}^{\dagger}[j]]=\delta_{i,j}$. The state of the entire system, including the waveguide and the two-level system, can be represented by a matrix-product state with the two-level system corresponding to the first site in the matrix product state, and the subsequent sites corresponding to the waveguide bins. The Hilbert space of the $k^{\text{th}}$ waveguide bin is spanned by the tensor product of the Hilbert spaces of the harmonic oscillators whose annihilation operators are $A_{+}[k]$ and $A_{-}[k]$. To compute the the state at $t=k\Delta t$ from the state at $t=(k+1)\Delta t$, we act on the matrix product state with the unitary operator $U[k+1,k]$ defined by:

[TABLE]

where

[TABLE]

where $n_{d}=\left\lfloor t_{d}/\Delta t\right\rfloor$ and $\Omega_{k}=\Omega(k\Delta t)\Delta t$. The application of $U[k+1,k]$ on the matrix product state at time step $k$ requires the implementation of a long-range gate, since it acts on the site corresponding to the two-level system, the $k^{\text{th}}$ waveguide bin and the $(k-2n_{d})^{\text{th}}$ waveguide bin. Following the approach introduced in Ref. Pichler and Zoller (2016), we implement this long range gate using a sequence of swap operations Shi et al. (2006) followed by a short-range gate corresponding to $U[k+1,k]$ Vidal (2003). The update is implemented using the tensor network state python library tncontract Darmawan and Grimsmo (2016) along with qutip Johansson et al. (2013)

 

We first validate our MPS update implementation git (2019) against the implementation introduced in Ref. Pichler and Zoller (2016) for an ideal mirror (i.e. $\theta=\pi/2$). Fig. 4 shows a comparison between the two implementation for two distinct settings — Fig. 4(a) in which the emitter is initially prepared in its excited state and allowed to decay into the waveguide without any external driving, and Fig. 4(b) in which the emitter is initially in its ground state and then driven by an exponentially decaying pulse ($\Omega(t)=\Omega_{0}e^{-\alpha t}$). We simulate both of these settings for different mirror phase $\varphi$ — as is known in such feedback systems with ideal mirrors, a properly chosen mirror phase can result in the emitter not decaying completely into the waveguide mode, rather exciting the bound state that exists between the emitter and the waveguide mode. We observe such incomplete decay for $\varphi=0$, and a complete decay of the emitter into the waveguide mode for other mirror phases. Moreover, the MPS update implementation presented in this section agrees perfectly with the MPS update implementation introduced in Ref. Pichler and Zoller (2016). This perfect agreement can be analytically explained by considering the normal modes of the point-coupling Hamiltonian describing the mirror (as described in section II.3). Using Eq. 30 for a perfect mirror ($\theta=\pi/2$), the two normal modes $b_{L,R}(\omega)$ can be expressed in terms of $a_{\pm}(x)$ via:

[TABLE]

Clearly, $b_{L}^{\dagger}(\omega)$ only creates excitations to the left of the mirror (which do not interact with the two-level system) while $b_{R}^{\dagger}(\omega)$ only creates excitations to the right of the mirror (which interact with the two-level system). Indeed, using Eq. 32, it can easily be seen that the interaction Hamiltonian between the waveguide and the two-level system ($H_{\text{wg-TLS}}$ defined in Eq. 34d) can be entirely expressed in terms of $b_{R}(\omega)$:

[TABLE]

which is identical to the interaction Hamiltonian used in Ref. Pichler and Zoller (2016).

Next, we study the impact of the non-ideality (i.e. the mirror reflection being less than 1) in the mirror on the dynamics of the feedback system. Changing the mirror reflection is equivalent to changing the parameter $\theta$ describing the mirror. Fig. 5 shows the impact of $\theta$ on the dynamics of the feedback system. Again, we simulate two different settings — Fig. 5(a) in which the emitter is initially prepared in its excited state and allowed to decay into the waveguide without any external driving, and Fig. 5(b) in which the emitter is initially in its ground state and then driven by an exponentially decaying pulse ($\Omega(t)=\Omega_{0}e^{-\alpha t}$). We note that unlike the case of an ideal mirror, having a less than unity reflectivity implies that there is no bound state that the emitter can decay into. Consequently, the emitter always decays to the ground state — however, as is seen in Fig. 5, the decay rate can be controlled by controlling the reflectivity of the mirror. While this seems reminiscent of the Purcell effect, we note that even with a less than unity reflectivity, the emitter does not decay exponentially into its ground state — this is a consequence of the non-Markovian nature of the feedback system.

 

Finally, we point out that an alternative to using the point-coupling Hamiltonian for analyzing the feedback system considered in this section is to directly formulate the Hamiltonian using the normal modes of the mirror-waveguide system. Such a formulation would also be able to capture the impact of the mirror reflectivity on the dynamics of the two-level system if it was accounted for in the construction of the normal modes. The point-coupling Hamiltonian provides an alternative, and equivalent (as shown in section II.3), approach for modelling linear-optical devices. Apart from the didactic importance of this result, we expect it to be of utility in understanding complicated feedback systems with multiple linear optical devices providing feedback.

IV Conclulsion

This paper resolves the problem of calculating the Hamiltonian for a frequency-independent linear-optical device from its classical scattering matrix. It is shown that an application of the quantum scattering matrix corresponding to the proposed Hamiltonian on an input state is equivalent to applying the inverse of the classical scattering matrix of the linear-optical device on the annihilation operators in the input state. We also diagonalize the point-coupling Hamiltonian and provide a connection between the point-coupling Hamiltonian and the quantization of normal modes of the linear-optical device Dalton et al. (1999). Finally, we demonstrate the practical utility of the proposed Hamiltonian by using it to rigorously formulate an MPS-based update for a time-delayed feedback system wherein the linear-optical device providing feedback is described by a full scattering matrix as opposed to a hard boundary condition.

V Acknowledgements

R.T. acknowledges funding from Kailath Stanford Graduate Fellowship.

Appendix A Solving Heisenberg equations of motion for point coupling Hamiltonian

The Heisenberg equation of motion (Eq. 3) can be easily integrated from $t_{0}$ to $t_{0}+\tau$ (where $\tau>0$) to obtain:

[TABLE]

Therefore, it follows from Eq. 2 that $a_{n}(x_{n};t_{0}+\tau)$ is given by:

[TABLE]

where the function $\Theta(x_{1}\leq x\leq x_{2})$ is defined in Eq. 5. Imposing Eq. A at $x_{n}=x_{n}^{0}\ \forall\ n\in\{1,2\dots N\}$, we obtain:

[TABLE]

from which we can obtain $a_{n}(x_{n}^{0};t)$ in terms of operators at $t=t_{0}$:

[TABLE]

where V is a $N\times N$ Hermitian matrix formed by $\text{V}_{m,n}$ as its elements and I is the identity matrix of size $N$. Substituting Eq. 46 into Eq. A together with the substitution $x_{n}=x_{n}^{0}+y$, we obtain Eq. 4

Appendix B Calculating normal modes of the point coupling Hamiltonian

Since the normal mode annihilation operator $b_{n}(\omega)$ is expressible as a linear combination of the annihilation operators $a_{n}(x_{n})$, we assume the following ansatz for $b_{n}(\omega)$:

[TABLE]

where $\mathcal{F}_{n,m}(\omega,x)$ is a function that is to be determined. Using Eq. 28a, we obtain the following differential equation for $\mathcal{F}_{n,m}(\omega,x)$:

[TABLE]

Here $\bm{\mathcal{F}}(\omega,x)$ is a $N\times N$ complex matrix with the functions $\mathcal{F}_{n,m}(\omega,x)$ as its elements. The solution for Eq. 48 is of the form:

[TABLE]

To evaluate the matrices $\mathcal{F}_{+}(\omega)$ and $\mathcal{F}_{-}(\omega)$, we use the boundary condition obtained on integrating Eq. 48 across a small interval around $x=0$:

[TABLE]

from which it immediately follows that $\bm{\mathcal{F}}_{+}(\omega)=\bm{\mathcal{F}}_{-}(\omega)\textbf{S}^{\dagger}$, where S is the classical scattering matrix of the linear-optical device. Moreover, using Eq. 48 along with Eq. 28b, we obtain $\bm{\mathcal{F}}_{-}^{\dagger}(\omega)\bm{\mathcal{F}}_{-}(\omega)=\textbf{I}$ – any choice for $\bm{\mathcal{F}}(\omega)$ is valid as long as it satisfies this constraint. Choosing $\bm{\mathcal{F}}_{-}(\omega)=\textbf{I}$ and $\bm{\mathcal{F}}_{+}(\omega)=\textbf{S}^{\dagger}$, we obtain Eq. 30.

 

To derive Eq. 32, we assume the following anstaz for $a_{n}(x_{n})$ in terms of $b_{n}(\omega)$:

[TABLE]

where $\mathcal{G}_{n,m}(x,\omega)$ is a function that is to be determined. Note that $\mathcal{G}_{n,m}(x,\omega)=\sqrt{2\pi}[a_{n}(x+x_{n}^{0}),b_{m}^{\dagger}(\omega)]$ and therefore from Eq. 47 it follows that $\mathcal{G}_{n,m}(x,\omega)=\mathcal{F}_{m,n}^{*}(\omega,x)$. We thus immediately obtain:

[TABLE]

Substituting this into Eq. 51, we immediately obtain Eq. 32.

Appendix C Convergence of Matrix-product-state simulations

Simulating time-delayed feedback with MPS, as described in section III, requires discretizing the waveguide modes which are described as continuas of harmonic oscillators into waveguide bins (Eq. 38). This introduces a simulation parameter $\delta t$ which controls the coarseness of this discretization. Moreover, application of a unitary impacting multiple bins in the MPS requires us to perform Schmidt decompositions on the state, and neglect Schmidt vectors with coefficients smaller than a specified tolerance. For the MPS-based simulation to be accurate, it is necessary to ensure that it has converged with respect to these two parameters. In this appendix, we present convergence studies with respect to these parameters to show that the choice of $\Delta t$ ($=0.05/\gamma$) and tolerance (0.01) assumed in section III indeed results in accurate simulations.

 

For the convergence study, we study the problem of an emitter decaying into the waveguide with a time-delayed feedback provided by the mirror and calculate the probability amplitude $|\varepsilon(t)|$ of the two-level system being in its excited state at time $t$ calculated using $|\varepsilon(t)|^{2}=\langle\sigma^{\dagger}\sigma\rangle$. We also assume that the mirror is perfect with phase $\varphi=0$ (which, as can be seen from Fig. 4(a) results in the two-level system not completely decaying into the ground state). In this case, the time-dependence of the complex amplitude of the excited state $\varepsilon(t)$ is given by the following ordinary differential equation (ODE) Tufarelli et al. (2013):

[TABLE]

where $\gamma=\gamma_{+}+\gamma_{-}$ is the total decay rate of the two-level system into the forward and backward propagating waveguide modes and this equation is solved subject to the initial condition $\varepsilon(0)=1$. $\varepsilon(t)$ obtained on solving this ODE thus provides a benchmark simulation that the MPS simulation can be compared against to gauge its convergence.

 

Fig. 6 shows the results of the convergence studies. The dependence of $|\varepsilon(t)|$ obtained from the MPS update on the discretization time-step $\Delta t$ is shown in Fig. 6(a) — it can be seen that choosing $\Delta t<0.15/\gamma$ is sufficient to ensure that the MPS simulation has converged and agrees with the ODE simulation. A similar study performed with respect to the tolerance used in the Schmidt decomposition is shown in Fig. 6(b) from which it can be seen that a tolerance smaller than 0.1 is sufficient to ensure that the MPS simulation has converged and agrees with the ODE simulation.