Fast simulation for multi-photon, atomic-ensemble quantum model of linear optical systems addressing the curse of dimensionality

TL;DR

This paper introduces a novel simulation method for multi-photon quantum systems in linear optical setups, significantly reducing computational complexity and enabling detailed analysis of quantum phenomena like interference and entanglement.

Contribution

The authors develop a decomposition technique for multi-photon evolution operators, addressing the curse of dimensionality in quantum optical simulations.

Findings

Accurately simulates Hong-Ou-Mandel interference

Reproduces Bell-CHSH inequality violations

Provides visual insights into photon propagation

Abstract

Photons are elementary particles of light in quantum mechanics, whose dynamics can be difficult to gain detailed insights, especially in complex systems. Simulation is a promising tool to resolve this issue, but it must address the curse of dimensionality, namely, that the number of bases increases exponentially in the number of photons. Here we mitigate this dimensionality scaling by focusing on optical systems composed of linear optical objects, modeled as an ensemble of two-level atoms. We decompose the time evolutionary operator on multiple photons into a group of time evolution operators acting on a single photon. Since the dimension of a single-photon time evolution operator is exponentially smaller than that of a multi-photon one in the number of photons, the decomposition enables the multi-photon simulations to be performed at a much lower computational cost. We apply this method to basic single- and multi-photon phenomena, such as Hong-Ou-Mandel interference and violation of the Bell-CHSH inequality, and confirm that the calculated properties are quantitatively comparable to the experimental results. Furthermore, our method visualizes the spatial propagation of photons hence provides insights that aid experiment designs for quantum-enabled technologies.

Tables

Table 1. Table 1: Parameter list of the numerical calculation. The row and column names indicate corresponding figure names and parameters, respectively. The parameters 𝐫 ¯ ¯ 𝐫 \bar{{\bf r}} and 𝐤 ¯ ¯ 𝐤 \bar{{\bf k}} are shown in Eq. 22 . The sixth column show the parameter of optical objects, such as mirror (M), beamsplitter (BS), phase shifter (PS), scatterer, polarization beamsplitter (PBS), and polarization rotator (PR). The parameters D j subscript 𝐷 𝑗 D_{j} , ω j subscript 𝜔 𝑗 \omega_{j} , and D j , p subscript 𝐷 𝑗 𝑝 D_{j,p} are shown in Eqs. 12 and 17 .

	Space size	Grid	$\overline{𝐫}$	$\overline{𝐤}$	Optical components
Fig. 1 (MZ)	($10\pi$, $10\pi$)	(256, 256)	(1.7, 7.7)	(10.0, 0)	M: 88 atoms $\times$ 8 layers, $D_j=0.56$, ${\omega }_j=5.0$ BS: 88 atoms $\times$ 1 layers, $D_j=2.8$, ${\omega }_j=0.31$ PS: 120 atoms $\times$ 0$\sim$20 layers, $D_j=0.56$, ${\omega }_j=1.0$
Fig. 2 (Scatterer)	($20\pi$, $10\pi$)	(512, 256)	(2.0, $5\pi$)	(10.0, 0)	Scatterer: 4 atoms $\times$ 4 layers, $D_j=1.0$, ${\omega }_j=5.0$
Fig. 3 (HOM)	($15\pi$, $15\pi$)	(384, 384)	$\xi$: (5.0, $7.5\pi$) $\eta$: ($7.5\pi$, 5.0)	$\xi$: (5.0, 0) $\eta$: (0, 5.0)	BS: 128 atoms $\times$ 1 layers, $D_j=2.0$, ${\omega }_j=0.34$
Fig. 4 (Bell-CHSH)	($20\pi$, $10\pi$)	(512, 256)	$\xi$: ($10\pi$, $5\pi$) $\eta$: ($10\pi$, $5\pi$)	$\xi$: (-10.0, 0) $\eta$: (0, 10.0)	PBS: 148 atoms $\times$ 8 layers, $D_{j,V}=0.56$, ${\omega }_j=5.0$ PR: 128 atoms $\times$ 16 layers, $D_{j,S}=0.56$, ${\omega }_j=1.2$

Equations

$$|1_{\xi },1_{\eta }⟩={\hat{a}}_{\xi }^{†}{\hat{a}}_{\eta }^{†}|0⟩,$$

$$V_{\text{BS}}{\hat{a}}_{\xi }^{†}{V}_{\text{BS}}^{†}=\frac{1}{\sqrt{2}}({\hat{a}}_{\xi }^{†}+\mathrm{i}{\hat{a}}_{\eta }^{†})$$

$$V_{\text{BS}}{\hat{a}}_{\eta }^{†}{V}_{\text{BS}}^{†}=\frac{1}{\sqrt{2}}({\hat{a}}_{\eta }^{†}+\mathrm{i}{\hat{a}}_{\xi }^{†}).$$

$$V_{\text{BS}}\otimes V_{\text{BS}}|1_{\xi },1_{\eta }⟩=\frac{\mathrm{i}}{\sqrt{2}}(|2_{\xi },0_{\eta }⟩+|0_{\xi },2_{\eta }⟩).$$

$$p(\Delta x)=|⟨{\Phi }_0(T)⟩\Phi (\Delta x,T){|}^2,$$

$$\displaystyle p(\Delta x)=\frac{1}{2}\Bigl{(}1-\exp\Bigl{(}-\frac{\Delta x^{2}}{2\sigma^{2}}\Bigr{)}\Bigr{)},$$

$$S(\theta )=E({\theta }_{\mathbf{a}},{\theta }_{\mathbf{b}})+E({\theta }_{{\mathbf{a}}^{\prime }},{\theta }_{\mathbf{b}})+E({\theta }_{{\mathbf{a}}^{\prime }},{\theta }_{{\mathbf{b}}^{\prime }})-E({\theta }_{\mathbf{a}},{\theta }_{{\mathbf{b}}^{\prime }}),$$

$$P_{p{p}^{\prime }}({\theta }_{\mathbf{a}},{\theta }_{\mathbf{b}})=\frac{1}{2}|{\delta }_{p,p^{\prime }}-{\cos }^2{\theta }_{\mathbf{ab}}|,$$

$$S(\theta )=3\cos 2\theta -\cos 6\theta ,$$

$$\hat{\rho }(c)=|\Phi ⟩⟨\Phi |.$$

$$\displaystyle\hat{\rho}^{\prime}(c)=\sum_{p}\bra{p}\hat{\rho}(c)\ket{p}=\frac{1}{1+c^{2}}\Bigl{(}\ket{H}\bra{H}+c^{2}\ket{V}\bra{V}\Bigr{)}.$$

$$P_{p{p}^{\prime }}({\theta }_{\mathbf{a}},{\theta }_{\mathbf{b}},\mathbf{r})=\frac{1}{2}\sum _{{\mathbf{r}}^{\prime }}|\left(⟨1_{\mathbf{r},p}|⟨1_{{\mathbf{r}}^{\prime },p^{\prime }}|+⟨1_{{\mathbf{r}}^{\prime },p}|⟨1_{\mathbf{r},p^{\prime }}|\right)|\Phi (t)⟩{|}^2,$$

$$S(\theta ,\mathbf{r})=E({\theta }_{\mathbf{a}},{\theta }_{\mathbf{b}},\mathbf{r})+E({\theta }_{{\mathbf{a}}^{\prime }},{\theta }_{\mathbf{b}},\mathbf{r})+E({\theta }_{{\mathbf{a}}^{\prime }},{\theta }_{{\mathbf{b}}^{\prime }},\mathbf{r})-E({\theta }_{\mathbf{a}},{\theta }_{{\mathbf{b}}^{\prime }},\mathbf{r}).$$

$$\hat{h}$$

$${\hat{h}}_0$$

$${\hat{h}}_I$$

$$g(j,\mathbf{k})$$

$$|\varphi (t)⟩$$

$${\hat{a}}_{\mathbf{k}}|1_{\mathbf{k}}⟩=|0⟩,{\hat{a}}_{\mathbf{k}}^{†}|0⟩=|1_{\mathbf{k}}⟩,$$

$${\hat{a}}_j^{†}{\hat{a}}_{\mathbf{k}}^{†}|0⟩={\hat{a}}_{\mathbf{k}}^{†}{\hat{a}}_j^{†}|0⟩=0,$$

$$\mathrm{i}\frac{\partial |\varphi (t)⟩}{\partial t}$$

$$|\varphi (t)⟩$$

$$|{\varphi }_p(t)⟩$$

$${\hat{h}}_I$$

$$g_p(j,\mathbf{k})$$

$$|1_{\mathbf{k},F}⟩$$

$$|1_{j,F}⟩$$

$$\sum _{\mathbf{k},p^{\prime }\in \{F,S\}}|1_{\mathbf{k},p^{\prime }}⟩⟨1_{\mathbf{k},p^{\prime }}⟩\varphi (t)+\sum _{j,p^{\prime }\in \{F,S\}}|1_{j,p^{\prime }}⟩⟨1_{j,p^{\prime }}⟩\varphi (t).$$

$$e^{-\mathrm{i}\hat{h}\delta t}$$

$${\hat{h}}_I=\sum _{p,j}(W_{j,p}{\hat{a}}_{j,p}^{†}{\hat{a}}_{{\mathbf{r}}_j,p}+W_{j,p}^{\ast }{\hat{a}}_{j,p}{\hat{a}}_{{\mathbf{r}}_j,p}^{†}),$$

$$\sum _{\mathbf{k}}{g}_p(j,\mathbf{k}){\hat{a}}_{j,p}^{†}{\hat{a}}_{\mathbf{k},p}$$

$${\hat{a}}_{\mathbf{r},p}|1_{\mathbf{r},p}⟩=|0⟩,{\hat{a}}_{\mathbf{r},p}^{†}|0⟩=|1_{\mathbf{r},p}⟩.$$

Full text

††preprint: APS/123-QED

I Introduction

Quantum electrodynamics (QED) is a subset of quantum field theory [1] that describes the interaction between matter and the electromagnetic field. This theory treats light as particles called photons. Arguably, QED is the most relevant theory considering the length, time, and energy scale required to observe the fundamental forces. If the matter appears in the form of atoms and molecules, then the combined matter-light system becomes a quantum optical system [2]. The dynamics of quantum optical systems, which we call quantum optical dynamics (QOD), is also governed by the basic equations of QED. However, the equations are often impossible to solve analytically and even numerically. One of the most critical difficulties is so-called the curse of dimensionality, which comes from the dimension (the size of a basis) of the Hilbert space corresponding to a quantum optical system. A basis can be characterized by “modes”, which roughly corresponds to the different ways in which particles (both matter and photons) can be excited. The number of possible modes is infinite. Even if the number of modes is restricted, the size of the basis increases exponentially in the number of particles present in QOD. A common simplification is to limit the number of modes and/or the total number of particles of each mode to just a few. This approximation returns reasonable results if most contributing degrees of freedom can be confidently identified from previous experiments. Should that fail, we have no other means but to numerically solve high-dimensional equations to accurately analyze QOD.

Known QOD simulation techniques [3, 4, 5, 6, 7] fall into two categories: cavity QED [8] and wave packet dynamics [9]. In a typical cavity QED simulation, the photonic modes are restricted to a single standing wave mode in the cavity, and the atom-light interaction modeled by the Jaynes–Cummings Hamiltonian [10, 11]. This suffices to explain several physical phenomena such as vacuum Rabi oscillation. On the other hand, the wave-packet-dynamics approach treats photons as a localized wave packet. Previous works treat Gaussian-shaped photons propagating through one-dimensional waveguides and interacting with two-level atoms [9, 12]. However, the restriction to the one-dimensional waveguide hinders their application to spatially two or three-dimensional systems. Havukainen et al. [13] conducted numerical simulations of wave packet dynamics where a single photon propagates through an ensemble of atoms laid out in a two-dimensional space. The simulation divides the space into smaller grids and retain all the spatial modes up to the resolution determined by the size of the grids. The atoms are located at the grid points. The dimension of the Hilbert space handled by this simulation exceeds that of the others by several orders of magnitude. The high dimensionality allows this simulation to reproduce a variety of QOD such as reflection and interference of a single photon by mirrors, beamsplitters, and double slits.

All these previous methods address the curse of dimensionality by either limiting the spatial modes and/or the total number of particles. Nevertheless, certain characteristic quantum phenomena require multi-photon and spatially multi-dimensional QOD to simulate (e.g., entanglement). These phenomena are not only important for basic science, but also are the underpinnings of quantum-enabled applications in cryptography, sensing, and imaging [14, 15, 16, 17]. Polarization is a degree of freedom of photons that is used in many quantum applications due to the availability of means to manipulate this degree of freedom at quantum precision. It is desirable that QOD simulations incorporate polarization.

In this paper, we introduce a numerical method to analyze multi-photon, spatially multi-dimensional QOD. Our method includes the polarization degree of freedom, based on the Hamiltonian in the quantum optical systems proposed by Havukainen et al. [13]. We construct the time evolution operator of the multi-photon system by composing a group of time evolution operators, each of which acting on a particular photon. This treatment exponentially reduces the dimension of the time evolution operator that need to be computed, with respect to the number of photons (see Methods section for the details). The numerical stability of simulations is improved by implementing a symplectic integration of the QED equations based on the Suzuki-Trotter decomposition. All combined, we succeed to simulate basic one- and two-photon phenomena, namely, the Mach–Zehnder (MZ) interference, Hong-Ou-Mandel (HOM) interference and violation of Bell-CHSH inequality. We use this simulation method to visualize the photon propagation dynamics of these phenomena and to understand the physical origins of the computed results. We also simulate a photon directed toward a scattering object (scatterer) and discover a counter-intuitive behavior when detecting this photon after the scatterer. The simulation results are presented in Sec. II, followed by discussions in Sec. III. The details of the present methods and parameters used in the numerical simulations are summarized in the Methods section.

II Results

II.1 Mach–Zehnder interference

The MZ interference is used for a variety of applications in optics, including optical switches [18], modulators [19], sensors [20], and quantum computing [21]. A typical MZ interferometer uses a set of representative linear optical objects, namely, two beamsplitters (BS1, BS2), two mirrors (M1, M2) and one phase shifter (PS) deployed as shown in Fig. 1(a). The BS1 splits an incident beam of light into two beams, one which runs through M1 followed by PS and the other beam through M2. The two beams are then combined by BS2 to produce an interference pattern.

We simulate the MZ interference of a single photon. The optical objects are each implemented by an ensemble of atoms, indicated by a gray filled rectangle in Fig. 1(a). The parameters of the atoms are tuned so as to serve as the desired linear optical object (see Table 1 for specific parameters). At $t=0$, a single photon is generated just left of BS1 directed toward BS1. The probability $P_{\textrm{right}}$ that the photon is ejected from the right side of BS2 is determined by the phase shift $\varphi$ imposed by PS. A simplified theory predicts that $P_{\textrm{right}}(\varphi)=\cos^{2}(\varphi/2)$, assuming that a photon propagates in at most two modes at all times with no photons absorbed by the atoms.

Figure 1(a) shows snapshots of the time evolution in the case of $\varphi=\pi$. Because the photon is ejected only in the upward direction, the result confirms the solution $P_{\textrm{right}}(\pi)=0$. Figure 1(b) shows a comparison of $P_{\textrm{right}}(\varphi)$ simulated by the Suzuki-Trotter decomposition and the Runge-Kutta method used in the present and previous works, respectively (see Methods section for the details). The result of using the Suzuki-Trotter decomposition is in a good agreement with the theoretical prediction, despite the fact that more than two modes are involved in the process and photons are absorbed by the atoms at the intermediate time steps. On the other hand, the result of the Runge-Kutta method shows that the probabilities significantly deviate from the theoretical prediction, and some of the values exceed one. Furthermore, even using a shorter time step than that employed in the Suzuki-Trotter method, a numerical instability appears in QOD as shown in the inset of Fig. 1(b). We conclude that the present method offers more reliable simulation than previously possible.

II.2 Photon detection in the presence of scatterer

Given a practical use of a quantum sensing and imaging [14, 15], obstacles in a space may scatter photons and affect the detection accuracy. Such a system cannot be described by an idealized theory, where the photonic modes are limited to one or just a few. Here, we demonstrate scattering of a single photon by an obstacle (scatterer), as shown in Fig. 2(a), exploiting the simulation capacity of our method for high-dimensional QOD. The scatterer, which is composed of $4\times 4$ atoms, is located at the center of the simulation space and are shifted by $\Delta x$ and $\Delta y$ in the $x$ and $y$ direction, respectively. A detector is placed at the far right edge of the space, indicated by the gray rectangle. The probability $P$ of the propagated photon being detected is given by the photon number density within the detector region. The error rate is calculated by 1- $P$. We computed the error rate for various values of the widths of the detector region. In general, we expect that the scatterer should have less effect on the photon as the scatterer moves further away from the optical path. We shall see that this general intuition does not hold.

Figure 2(b) shows the dependence of the error rates on the widths of the detector window. We change $\Delta x$ but fix $\Delta y=0.0$. As a general trend, the error rates monotonously decrease as the width increases. This is because the wider the width of the detector region, the more likely it is to find the photon inside the detector region. Yet, we observe a plateau profile of the error rate in $\Delta x=11.0$ around the width indicated by the vertical gray solid line. The number densities at the final states are visualized for $\Delta x=-7.4$ and $11.0$ in the lower panels of Fig. 2(b). When $\Delta x=11.0$, the number density has a particular structure, that is caused by an interference of a scattered photon coming from the upper and lower side of the scatterer. This fringe structure gradually disappears by diffusion of the photon number density as the distance to the scatterer increases after passing it, as in the case of $\Delta x=-7.4$. The plateau profile appears when an edge of the detector region is in the sparse region of the interference fringes. In a sense, the detector “fails” to capture more photons despite increasing its size.

Figure 2(c) compares the profiles of the error rates by changing $\Delta y$ values and fixing $\Delta x=7.4$. As in Fig. 2(b), the error rates generally decrease as the width increases. When $\Delta y$ increases, the error rates decrease because the scatterer moves further away from the optical path. We observe that the order of the error rate curves changes as the detector region width changes (cf. the black arrow). The result shows that there is a range of the widths where the error rate increases when the obstacle moves away from the center of the optical path (also see the inset of Fig. 2(c)). To understand this counter-intuitive behavior, we visualize the final states of $\Delta y=0.0,0.6$, and $1.2$ as shown in the lower panels of Fig. 2(c). In the case of $\Delta y=0$, the interference fringes almost completely land within the detector region. On the other hand, a shift $\Delta y$ destroys the fringe structure. Due to the cancellation of the interference, a part of the wave-packet distributes outside the detector region (indicated by the white arrows in the lower panels of Fig. 2(c)), resulting in the higher error rate than the $\Delta y=0$ case.

II.3 Hong-Ou-Mandel interference

The Hong-Ou-Mandel (HOM) interference is a quintessential quantum effect of two indistinguishable photons [22, 23], which cannot be analyzed by simulations based on classical electrodynamics or QOD simulation limited to a single-photon states. Figure 3(a) shows a minimal model for the HOM interference, where two photons characterized by $\xi$ and $\eta$ indices are simultaneously injected to the ports of a beamsplitter. The beamsplitter divides a injected photon in half into transmitted and reflected parts. In QED, the injected state is expressed by

[TABLE]

where $\hat{a}^{\dagger}$ is the creation operator of a photon. The operator of the beamsplitter $V_{\textrm{BS}}$ performs as,

[TABLE]

Therefore, the outgoing state from the beamsplitter becomes

[TABLE]

This result indicates that the two photons are always ejected together from either right ($\xi$) or upper ($\eta$) outlet port. In other words, the probability that each photon is emitted in a separate port disappears by quantum interference of the two photons.

This HOM interference is simulated by the model shown in Fig. 3(a). The calculated quantum state is denoted by $\ket{\Phi(\Delta x,t)}$, where $\Delta x$ is shift in the x direction of the initial position of the $\xi$ photon. First, we visualize the spatial distribution of the probability of finding two photons at a given location at different times of the evolution. This bunching probability can be calculated by $\rho(\Delta x,t,{\bf r})=|\bra{1_{\bf r}}\braket{1_{\bf r}}{\Phi(\Delta x,t)}|^{2}$, where $\ket{1_{\bf r}}\ket{1_{\bf r}}$ is a state in which two photons are in the same position. Figure 3(b) shows the dynamics of $\rho(0,t,{\bf r})$. Before the two photons reach at the beamsplitter ($t<17$), there is almost no value of the bunching probability because they are far from each other. After the photons pass through the beamsplitter, we observe that the bunching appears and then the distribution is separated into two portions that travel in the right and upper direction. Then, we perform this numerical experiment for different values of $\Delta x$. We remove the beamsplitter in Fig. 3(a) and calculate the two-photons system. This free propagating state is denoted by $\ket{\Phi_{0}(t)}$ in which the two photons are observed separately in the different outlets, which we call the coincidence probability. The coincidence probability can be measured by

[TABLE]

where time $T=45$ is chosen to assure that the wave-packets have separated enough after passing through the beamsplitter. Figure 3(c) shows the profile of $p(\Delta x)$. Around $\Delta x=0$, the coincidence probability is almost zero which corresponds to the theoretical result in Eq. 2. Increasing $\Delta x$, we observe a dip in the coincidence probability. This HOM dip was observed in experiments as an evidence of the quantum nature of light [22, 23]. The shape of the dip is characterized by an overlap of the two photons. In fact, a theory of the HOM interference offers an analytical solution [22].

[TABLE]

where $\sigma$ is the Gaussian width of the photon. This analytical solution, which is plotted in Fig. 3(c), shows a good agreement with the numerical one.

II.4 Violation of Bell-CHSH inequality

In addition to the multi-photon states, the photon polarization degree of freedom is included in the present method. The polarization is commonly chosen in quantum experiments, including the photonic experiments confirming the violation of Bell-CHSH inequality [24, 25]. The CHSH inequality provides a limit on a particular type of correlations of two separated systems, should their behavior be determined by classical mechanics. Quantum mechanics, however, violates the inequality under specific conditions. This violation was successfully confirmed by photonic experiments [26, 27, 28] and various physical systems including atoms [29] and superconducting qubits [30].

Here, we model the experiment performed by Aspect et al. [31], as shown in Figs. 4(a) and (b). We deploy optical objects that rotate the polarization of the photon passing through them. One of the polarization rotators shifts by $\theta_{\mathbf{a}}$ or $\theta_{\mathbf{a}^{\prime}}$, while the other by $\theta_{\mathbf{b}}$ or $\theta_{\mathbf{b}^{\prime}}$. Two photons are directed to the two rotators individually, and then the horizontal and vertical components are spatially separated by the polarization beamsplitters. Given the rotation angles, we observe a correlation of the polarizations of the two photons as defined by

[TABLE]

where $E(\theta_{\mathbf{a}},\theta_{\mathbf{b}})=P_{HH}(\theta_{\mathbf{a}},\theta_{\mathbf{b}})+P_{VV}(\theta_{\mathbf{a}},\theta_{\mathbf{b}})-P_{HV}(\theta_{\mathbf{a}},\theta_{\mathbf{b}})-P_{VH}(\theta_{\mathbf{a}},\theta_{\mathbf{b}})$. The probability $P_{pp^{\prime}}(\theta_{\mathbf{a}},\theta_{\mathbf{b}})$ ($p,p^{\prime}\in\{H,V\}$) is the probability that a photon is detected in the $p$-polarization outlet in the left and, simultaneously, another photon detected in the $p^{\prime}$-polarization outlet, when the rotation angles are set at $\theta_{\mathbf{a}}$ and $\theta_{\mathbf{b}}$. According to the Bell-CHSH inequality, $S(\theta)$ must range $-2\leq S(\theta)\leq 2$ if the local realism is correct. However, if we prepare two photons that are in a maximally entangled state, quantum mechanics asserts that

[TABLE]

where $\theta_{\mathbf{ab}}=|\theta_{\mathbf{a}}-\theta_{\mathbf{b}}|$ and $\delta$ is the Kronecker delta [32]. Equation 6 leads to $E(\theta_{\mathbf{a}},\theta_{\mathbf{b}})=\cos 2\theta_{\mathbf{ab}}$. When setting the angles at $\theta_{\mathbf{ab}}=\theta_{\mathbf{a^{\prime}b}}=\theta_{\mathbf{a^{\prime}b^{\prime}}}=\theta$ (see also the inset of Fig. 4(c) ), and hence $\theta_{\mathbf{ab^{\prime}}}=3\theta$, Eq. 5 becomes

[TABLE]

in the case of the maximally entangled state. Equation 7 breaks the inequality $-2\leq S(\theta)\leq 2$ as shown by the blue solid line of Fig. 4(c).

Figure 4(a) shows the dynamics of the number density of two photons prepared in a maximally entangled state. Each photon propagates to the left and right direction, as indicated by the white arrows. The angles of the two polarization rotators are set at $\theta_{\bf a}=0$ and $\theta_{\bf b}=\pi/2$. The initial state is given by $\ket{\Phi(t=0)}=(\ket{H}\ket{H}+\ket{V}\ket{V})/\sqrt{2}$, where $\ket{H}\ket{H}$ ($\ket{V}\ket{V}$) denotes the two photons polarized in the horizontal (vertical) direction. Figure 4(b) shows the case of a product state $\ket{\Phi(t=0)}=\ket{H}\ket{H}$. The photons of the product change their polarization angle by passing through the polarization rotator, while not for those of the maximally-entangled state. This invariance can be explained by computing the state of one of the photons when the initail state of the two photons is $\ket{\Phi}=\frac{\ket{H}\ket{H}+c\ket{V}\ket{V}}{\sqrt{1+c^{2}}}$, where $0\leq c\leq 1$. The density matrix corresponding to the two photons is

[TABLE]

A partial trace of the density matrix on the right photon yields the reduced density matrix

[TABLE]

If $c=1$ which corresponds to the maximally-entangled state, $\hat{\rho}^{\prime}(1)=\hat{I}/2$ where $\hat{I}$ is the identity operator. Considering a polarization rotation $\hat{U}$ that operates $\hat{U}\rho^{\prime}(c)\hat{U}^{\dagger}$, it is obvious that the polarization rotator does not change the partial state of this photon. On the other hand, we have $\hat{\rho}^{\prime}(0)=\ket{H}\bra{H}$, thus the state of the photon of a product state is changed by the operation of $\hat{U}$. These results are visualized in Figs.  4(a) and (b).

Figure 4(c) shows the simulated $S(\theta)$. The results of the maximally entangled state (the blue plots) show a good agreement with the theoretical prediction (the blue solid line) of Eq. 7. The violations of the Bell-CHSH inequality appear around $\theta=\pi/8,\ 3\pi/8$. In addition, we calculate $S(\theta)$ by changing the parameter $c$ of the initial state as in Eq. 8. The violations occur if $c\geq 0.25$.

To obtain deeper insights into the difference between the maximally entangled and product states, we visualize $S(\theta)$ in Fig 4(c). Namely, we first define

[TABLE]

and $E(\theta_{\mathbf{a}},\theta_{\mathbf{b}},{\bf r})=P_{HH}(\theta_{\mathbf{a}},\theta_{\mathbf{b}},{\bf r})+P_{VV}(\theta_{\mathbf{a}},\theta_{\mathbf{b}},{\bf r})-P_{HV}(\theta_{\mathbf{a}},\theta_{\mathbf{b}},{\bf r})-P_{VH}(\theta_{\mathbf{a}},\theta_{\mathbf{b}},{\bf r})$. The basis $\ket{1_{{\bf r},p}}$ is a state in which a photon with the polarization $p$ is at the position r. Accordingly, a correlation density is defined by

[TABLE]

Note that we can retrieve $S(\theta)$ by $S(\theta)=\sum_{\bf r}S(\theta,{\bf r})$. The distribution of $S(\theta,{\bf r})$ for the maximally entangled state is found to be invariant with respect to $\theta$ except for its magnitude. This result is relevant to Fig. 4(a) in which the number density apparently does not change by the polarization rotators due to the fact that the reduced density matrix becomes the identity operator. Therefore, the correlations between the sets of $\theta_{\mathbf{a}}$ and $\theta_{\mathbf{b}}$ tend to be magnified, resulting in clear violations of the inequality at particular angles. On the other hand, since the product state is changed by the polarization rotator as shown in Fig. 4(b), terms composed of the correlations tend to cancel each other. This fact prevents $S(\theta)$ to exceed the limit of the local realism. In short, the present method reproduces the theory and experimental results, and its visualization facilitates understanding of the mechanism.

III Discussion

The curse of dimensionality does not exist if the system is modeled and solved by classical electrodynamics. Many aspects of optics can be analyzed by classical methods of electromagnetic fields, such as method of moments [33], finite-difference time-domain methods [34, 35, 36], finite-difference frequency-domain methods [37], and finite element method [38]. These techniques are based on the Maxwell equations that approximate light-matter interactions as a dissipation of the electromagnetic field or as spatially and spectrally nonuniform permittivity and/or permeability (so-called macroscopic treatment). However, the classical framework cannot capture multi-photon phenomena which are demonstrated in this study, as well as nonlinear optical effects such as the spontaneous parametric down conversion. Our method simulates an optical system as a high-dimensional quantum system, solving the quantum version of electrodynamical equations, and hence is inherently compatible with these extensions.

While any QOD simulation has to contend with some omission of details when used to understand experiments, visualization of quantum dynamics of complex systems, as treated in this paper, is a unique advantage of numerical simulation, not present in experiments. A quantum state changes itself in principle when observed. It is, therefore, significantly resource intensive to experimentally trace the time-evolution process of photons propagating through complex optical setups. The presented high-dimensional QOD would be useful to reveal the intermediate process in detail. For instance, the visualizations of Sec. II.2 helped us to clarify the cause of the unexpected behavior of the error rate. This flexible and accurate simulation of the high-dimensional QOD should also aid in creating a viable minimal model that guides designing experiments and applications that exploit quantum states of light.

IV Methods

We explain the present method in a two-dimensional space that spans over $0\leq x\leq L_{x}$ and $0\leq y\leq L_{y}$. The space is discretized by a grid, and its boundaries are periodic. We set $\hbar=c=1$ in the following formulations for a notational simplicity.

IV.1 Single-photon system

We begin with a theoretical framework of the single-photon simulation as introduced by [13] (also see the Supplemental Information for the details). The system consists of $N_{A}$ atoms and one photon that has a wave-number vector ${\bf k}$. The total Hamiltonian is

[TABLE]

where

[TABLE]

and $\hat{a}$ is an annihilation operator. The frequencies $\omega_{{\bf k}}$ and $\omega_{j}$ indicate eigenenergies of the corresponding photon mode and atom, respectively. The Hamiltonian $\hat{h}_{I}$ represents a dipole-dipole interaction between the photon and atoms,

[TABLE]

where $L=\sqrt{L_{x}L_{y}}$. The parameter $D_{j}$ sets the strength of the dipole-dipole interaction, and the vector ${\bf r}_{j}$ the position of the $j$th atom. This interaction makes the constituent atoms play a role of linear optical objects such as mirror, beamsplitter, and scatterer by controlling the parameters $\omega_{j}$ and $D_{j}$.

A quantum state of this system is restricted so that the total number of excitations is one, namely,

[TABLE]

where $c(t,{\bf k})$ and $c_{j}(t)$ are the complex amplitudes satisfying $\sum_{{\bf k}}|c(t,{\bf k})|^{2}+\sum_{j}|c_{j}(t)|^{2}=1$. The atoms are deployed at the grid lattice points to form necessary optical objects. The indices $1_{{\bf k}}$ and $1_{j}$ in the kets denote one photon of ${\bf k}$ mode and one excitation of a $j$-th atom, respectively. Namely, operations of $\hat{a}_{\bf k}$ to the bases are

[TABLE]

and those of $\hat{a}_{j}$ are the same manner. In addition,

[TABLE]

due to the restriction of the total number of excitations to one.

The Schr$\ddot{\textrm{o}}$dinger equation and the solution of the time evolution are written by

[TABLE]

IV.2 Polarization degrees of freedom

The method is extended simply by adding a polarization index $p$ that specifies the horizontal and vertical polarization by $H$ and $V$, respectively. The state of the single-photon system is expressed by

[TABLE]

The Hamiltonians $\hat{h}_{0}$ and $\hat{h}_{I}$ in Eqs. 10 and 11 can be extended by just replacing $\hat{a}_{{\bf k}}\rightarrow\hat{a}_{{\bf k},p},\hat{a}_{j}\rightarrow\hat{a}_{j,p}$. Note that these Hamiltonians do not influence the polarization degrees of freedom. The interaction Hamiltonian $\hat{h}_{I}$ becomes

[TABLE]

where $D_{j,p}$ denotes the strength of the dipole interaction between the $j$-th atom and the photon mode $p$. Tuning $D_{j,p}$ creates optical objects that change the polarization. For example, a polarization beamsplitter reflects only the vertical polarization but transmits the horizontal one. It can be created by setting $\omega_{j}$ and $D_{j,V}$ at appropriate values for a role of mirrors, but $D_{j,H}=0$. A polarization rotator can be constructed by a rotated half wave plate which introduces $\pi$ phase shift on the polarization directed toward the slow axis of the wave plate. This can be achieved by using the basis of the fast and slow axes of the wave plate $\{F,S\}$, as

[TABLE]

The basis transformation of the single-photon system is

[TABLE]

Then, $\Theta$ is set as $\Theta=\theta_{\textrm{rot}}/2+\theta_{\textrm{pol}}$ depending on a rotation angle $\theta_{\textrm{rot}}$ given by the polarization rotator, where $\theta_{\textrm{pol}}$ denotes an angle of the linear polarization of the input photon against a direction of the horizontal polarization. $\omega_{j}$ and $D_{j,S}$ are set at appropriate values to add $\pi$ phase shift on the slow axis component, while $D_{j,F}=0$.

This scheme for the polarization is examined by a single-photon simulation. The time evolution is solved by the Suzuki-Trotter decomposition which will be explained in the Sec. IV.3. Figure 5(a) shows the system where the initial polarization of a photon is $H$, and the photon is directed toward a polarization rotator. In this case, the polarization is rotated by an angle $\theta_{\textrm{rot}}=\pi/4$. Then, only the $V$ component is reflected to the upward by the polarization beamsplitter, while the $H$ component passes through the object. We observe the probabilities of the $V$ components as a function of $\theta_{\textrm{rot}}$, as shown in Fig. 5(b). The calculated probabilities are well-fitted by the behavior of an ideal polarization rotator $\sin^{2}\theta_{\textrm{rot}}$. These results confirm that the quantum state including the polarization and the relevant optical objects used in Sec.II.4 work as inteded.

IV.3 Time evolution by Suzuki-Trotter decomposition

The previous study [13] rewrites the Schr$\ddot{\textrm{o}}$dinger equation in Eq. 13 in the interaction picture and numerically solves it by the fourth-order Runge-Kutta method. We have decided to formulate the Hamiltonian in such a way to apply the Suzuki-Trotter decomposition as it is known that the Suzuki-Trotter decomposition has numerically advantages in computing the time evolution of a Hamiltonian system [39, 40, 41, 42, 43]. The time evolution operator in Eq. 14 can be approximated by the Suzuki-Trotter decomposition, as

[TABLE]

where $\delta t$ is a time step of the simulation. Each Hamiltonian is diagonalized to compute the time evolution operator. The $\hat{h}_{0}$ is already diagonal in the wave-number basis ${\bf k}$. To diagonalize $\hat{h}_{I}$, it is numerically convenient to transform the interaction Hamitonian $\hat{h}_{I}$ from the basis ${\bf k}$ into the position basis ${\bf r}$

[TABLE]

because by using a Fourier transformation the term in Eq. 11 becomes

[TABLE]

where $W_{j,p}=-\mathrm{i}D_{j,p}\sqrt{N\omega_{j}}/\sqrt{2}L$. The operation of $\hat{a}_{{\bf r},p}$ is

[TABLE]

As seen in Eq. 20, $\hat{h}_{I}$ exchanges the photon and the excitation of the atom located at the same position. This compact representation thanks to the ${\bf r}$ basis is advantageous in the numerical diagonalization. In fact, multiplying $\ket{\phi}$ in Eq. 15 gives

[TABLE]

By using matrix representation, this is equivalent to

[TABLE]

where $V=2^{-1/2}\begin{pmatrix}1&-\mathrm{i}\\ 1&\mathrm{i}\end{pmatrix}$ is a unitary matrix and we used the fact that $W_{j,p}^{*}=-W_{j,p}$. Therefore, the time evolution can be operated, as

[TABLE]

We assessed a numerical performance of the Suzuki-Trotter decomposition by a simple single-photon simulation without the polarization degree of freedom. As shown in Fig. 6(a), the injected photon collides with the tilted mirror at about $t=10$ and is reflected upwards. Then, due to the periodic boundary condition, the photon emerges from the lower region and is reflected by the mirror again at about $t=35$ to go in the right direction. Figure 6(b) shows the trajectories of the probability $\braket{\phi(t)}{\phi(t)}$ that ideally should remain $1$ all the time. The result obtained by the previous scheme based on the Runge-Kutta method increases significantly above $1$. On the other hand, the present scheme keeps the probability at $1$, even with a larger time step. This is a benefit of the Suzuki-Trotter decomposition in Eq. 19 that uses only unitary operators. Figure 6(c) shows the total energy of the quantum state. By using the previous scheme, the total energy increases and violates the law of energy conservation. The energy trajectory by the present scheme shows small deviations from the initial energy at the times when the photon interacts with the mirror. Afterwards, however, the energy recovers to the initial energy, and the energy conservation tends to be maintained in the dynamics. Such a favourable feature is known as the result of a symplectic condition the Suzuki-Trotter decomposition has [42, 43]. The time evolution by the Suzuki-Trotter decomposition significantly improves stability of the calculation, which gives a reliable quantitative evaluation by the QOD simulation.

IV.4 Multi-photon system

For notation simplicity, this subsection does not write the polarization index $p$ explicitly, because the procedure with $p$ is almost the same, just replacing ${\bf r}\rightarrow{\bf r},p$ and $j\rightarrow j,p$. When naively written, an $N$-photon state without the polarization is

[TABLE]

A critical issue is that handling the coefficient $c(t,{\bf r}_{1},{\bf r}_{2},\cdots{\bf r}_{N})$ requires at least $O(M^{N})$ of numerical operations and memory size, where $M$ is the number of grids.

To expedite the calculation, we exploit the fact that $N$-partite states such as $\ket{\Phi(t)}$ can be expressed with a fewer number of terms with a suitable choice of local states [44]. In addition our simulation only use the atoms to mimic linear optical objects such as mirrors and beamsplitters. These optical objects only need to reproduce their appropriate output electromagnetic wave for a given incoming wave. The atoms can effectively induce interactions between photons, but linear optical elements do not exhibit such photon-photon interaction. Based on this intuition, we simulate the time evolution of each photon separately neglecting the presence of the other photons in doing so. More precisely, we consider $N$ atoms virtually located at the same position ${\bf r}_{j}$, which we call virtual atoms, and let the $N$ atoms only interact with their respective partner photon to simulate $N$-photon states. This treatment limits interaction $\hat{h}_{I}$ within a pair of a photon and the corresponding virtual atom. Our $N$-photon state becomes

[TABLE]

where $\xi_{n}$ is a sets of the properties of $n$-th photon. Deviation of this treatment from the true $N$-photon evolution occurs when two photons are present at the same location and an atom located at that point absorbs the photons. In our simulation, photon densities are kept sufficiently low to suppress this error. These tricks render a time evolution operator of the $N$-photon system to be a direct product form as

[TABLE]

Note that the dimension of each $\exp(-\mathrm{i}\hat{h}t)$ is $M^{2}$, whereas the naive approach results in treating a time evolution operator of size $M^{2N}$. This reduction of the dimensionality dramatically improves speed of the time evolution and makes the multi-photon simulation feasible.

IV.5 Initial states

All the initial states of the atoms at $t=0$ are set in their ground states. For example, $c_{j}(0)=0$ in the case of the single-photon system, and a Gaussian-shaped photon is injected in the space as

[TABLE]

where $\sigma$ is a width of the Gaussian. The vectors $\bar{{\bf k}}$ and $\bar{{\bf r}}$ are parameters to decide the initial velocity and position, respectively. For the simulations with the polarization, the initial state is given by

[TABLE]

where $\theta$ is an initial polarization angle and $\delta$ is the Kronecker delta.

As in Sec. II.3 and Sec. II.4, the two Gaussian-shaped photons are injected at $t=0$. The modes of the photons are characterized by $\xi$ and $\eta$, where $\xi$ and $\eta$ indicate sets of the photon properties that are $\sigma$, $\bar{{\bf k}}$ and $\bar{{\bf r}}$. For the HOM interference, we used an initial state as

[TABLE]

where $(\xi\leftrightarrow\eta)$ refers to exchange of the indices in the previous term, which originates from the commutative relation of bosons. For the test of the Bell-CHSH inequality, the initial state entangles their polarization $p$.

[TABLE]

IV.6 Computational details

Table 1 summarizes the parameters list that are used in the calculation shown in the main text. All the simulations are performed by the time step $\delta t=0.1$ and Gaussian width $\sigma=2.0$.

Supplementary information

Hamiltonian of single-photon system

In this section, we prepossess a Hamiltonian of the one-photon system proposed in a previous study [13] in a form suitable for the main text. According to Ref. [13], the one-photon Hamiltonian that involves $N_{A}$ atoms is defined by

[TABLE]

The Hamiltonian of the free propagation of the photon, $\hat{h}_{F}$, consists of an annihilation operator $\hat{a}_{{\bf k}}$ for one photon that has a wave-number vector ${\bf k}$. That of the atoms, $\hat{h}_{A}$, is composed of the Pauli z operator $\hat{\sigma}_{j}$ of the $j$-th atom. The parameters $\omega_{{\bf k}}$ and $\omega_{j}$ indicate eigen energies of the photon and atom, respectively. The $\hat{h}_{I}$ represents interaction between the photon and atoms.

A quantum state of the corresponding system is written by

[TABLE]

where $c({\bf k})$ and $c_{j}$ are coefficients of the probability amplitudes. When operating $\hat{h}_{A}$ to the state, we obtain

[TABLE]

where $\omega=\sum_{j}\omega_{j}/N_{A}$, and we defined $\hat{h}_{a}$ as

[TABLE]

Therefore, the Hamiltonian becomes

[TABLE]

where $\hat{h}_{0}=\hat{h}_{F}+\hat{h}_{a}$. Because the constant energy $-N_{A}\omega$ does not influence the time evolution, we omit this term when displaying the Hamiltonian in the main text.

The Hamiltonian $\hat{h}_{I}$ that describes the dipole-dipole interaction between the atoms and photons as

[TABLE]

We set the frequencies of the atoms and photon to be in resonance in simulation, thus $g(j,{\bf k})$ is approximated as

[TABLE]