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

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.


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.One of the most critical difficulties in solving the equations of QOD is 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, independently.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 matter-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][13][14][15].However, the restriction to the onedimensional waveguide hinders their application to spatially two or three-dimensional systems.Havukainen et al. [16] conducted numerical simulations of wave packet dynamics where a single photon propagates through an ensemble of two-level atoms laid out in a two-dimensional space.The simulation divides the space into smaller grids and retains all the spatial modes up to the resolution determined by the size of the grids.The two-level atoms are located at the grid points (see the top left panel of Fig. 1(a)).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.For example, Havukainen et al. [16] considered a large number of possible spatial modes, but they treated the dynamics of one photon.The multi-photon system gives rise to characteristic quantum phenomena that involve quantum entanglement.These phenomena are not only important for basic science, but also are the underpinnings of quantum-enabled applications in cryptography, sensing, and imaging [17][18][19][20].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. [16].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 needs 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 present the detailed interplay between the detection pattern and the single-photon interference caused by 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.

A. Mach-Zehnder interference
The MZ interference is used for a variety of applications in optics, including optical switches [21], modulators [22], sensors [23], and quantum computing [24].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 two-level atoms, indicated by a gray filled rectangle in Fig. 1(a).The parameters of the two-level atoms are tuned so as to serve as the desired linear optical object (see Table I for specific parameters).At t = 0, a single photon is generated just left of BS1 directed toward BS1.The probability P right that the photon is ejected from the right side of BS2 is determined by the phase shift φ imposed by PS.A simplified theory predicts that P right (φ) = cos 2 (φ/2), assuming that a photon propagates in at most two modes at all times with no photons absorbed by the two-level atoms.
Figure 1(a) shows snapshots of the time evolution in the case of φ = π.Because the photon is ejected only in the upward direction, the result confirms the solution P right (π) = 0. Figure 1(b) shows a comparison of P right (φ) 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 two-level 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.

B. Photon detection in the presence of scatterer
Given a practical use of a quantum sensing and imaging [17,18], 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 × 4 two-level atoms, is located at the center of the simulation space and is shifted by ∆x and ∆y The contour shows the photon number density.The gray lines and filled rectangles correspond to the beamsplitters (BS1, BS2) and mirrors (M1, M2), respectively.The unfilled rectangle corresponds to the phase shifter (PS), which adds the phase-shift exp(iφ) to the photon.(b) Plot of the probability P right (φ) obtained by the Suzuki-Trotter (blue) and the Runge-Kutta (orange) method.dt is the time step used in each method.The black dotted curve shows the standard theoretical prediction.The horizontal gray dotted line is drawn at P right = 1 as a guide for the eye.
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 ∆x but fix ∆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 ∆x = 11.0 around the width indicated by the vertical gray solid line.
The number densities at the final states are visualized for ∆x = −7.4 and 11.0 in the lower panels of Fig. 2(b).When ∆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 ∆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 ∆y values and fixing ∆x = 7.4.As in Fig. 2(b), the error rates generally decrease as the width increases.When ∆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 scatterer moves away from the center of the optical path (also see the inset of Fig. 2(c)).To understand this behavior, we visualize the final states of ∆y = 0.0, 0.6, and 1.2 as shown in the lower panels of Fig. 2(c).In the case of ∆y = 0, the interference fringes almost completely land within the detector region.On the other hand, a shift ∆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 ∆y = 0 case.

C. Hong-Ou-Mandel interference
The Hong-Ou-Mandel (HOM) interference is a quintessential quantum effect of two indistinguishable photons [25,26], 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 ξ and η 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 where â † is the creation operator of a photon.The operator of the beamsplitter V BS performs as, Therefore, the outgoing state from the beamsplitter becomes This result indicates that the two photons are always ejected together from either right (ξ) or upper (η) 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 |Φ(∆x, t)⟩, where ∆x is a shift in the x direction of the initial position of the ξ 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 ρ(∆x, t, r) , where |1 r ⟩ |1 r ⟩ is a state in which two photons are in the same position.Figure 3(b) shows the dynamics of ρ(0, t, 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 ∆x.We remove the beamsplitter in Fig. 3(a) and calculate the two-photons system.This free propagating state is denoted by |Φ 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 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(∆x).Around ∆x = 0, the coincidence probability is almost zero which corresponds to the theoretical result in Eq. 2. Increasing ∆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 [25,26].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 [25].
where σ 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.

D. 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 [27,28].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 [29][30][31] and various physical systems including atoms [32] and superconducting qubits [33].
Here, we model the experiment performed by Aspect et al. [34], 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 θ a or θ a ′ , while the other by θ b or θ b ′ .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 where is the probability that a photon is detected in the p-polarization outlet in the left and, simultaneously, another photon detected in the p ′ -polarization outlet, when the rotation angles are set at θ a and θ b .According to the Bell-CHSH inequality, S(θ) must range −2 ≤ S(θ) ≤ 2 if the local realism is correct.However, if we prepare two photons that are in a maximally entangled state, quantum mechanics asserts that where θ ab = |θ a − θ b | and δ is the Kronecker delta [35].Equation 6 , where 0 ≤ c ≤ 1.The density matrix corresponding to the two photons is A partial trace of the density matrix on the right photon yields the reduced density matrix If c = 1 which corresponds to the maximally entangled state, ρ′ (1) = Î/2 where Î is the identity operator.Considering a polarization rotation Û that operates Û ρ ′ (c) Û † , it is obvious that the polarization rotator does not change the partial state of this photon.On the other hand, we have ρ′ (0) = |H⟩ ⟨H|, thus the state of the photon of the product state is changed by the operation of Û .These results are visualized in Figs.4(a) and (b).
Figure 4(c) shows the simulated S(θ).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 θ = π/8, 3π/8.In addition, we calculate S(θ) by changing the parameter c of the initial state as in Eq. 8.The violations occur if c ≥ 0.25.
To obtain deeper insights into the difference between the maximally entangled and product states, we visualize S(θ) in Fig 4(c).Namely, we first define and E(θ a , θ b , r) = P HH (θ a , θ b , r) + P V V (θ a , θ b , r) − P HV (θ a , θ b , r) − P V H (θ a , θ b , r).The basis |1 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 Note that we can retrieve S(θ) by S(θ) = r S(θ, r).The distribution of S(θ, r) for the maximally entangled state is found to be invariant with respect to θ 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 θ a and θ 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(θ) 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 [36], finite-difference time-domain methods [37][38][39], finite-difference frequency-domain methods [40], and finite element method [41].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 (see the recent review [42]), such as quantum interference and quantum entanglement, which are demonstrated in this study.By modeling a three-level atomic system, our numerical scheme can simulate nonlinear optical effects such as spontaneous parametric down conversion.However, the curse of dimensionality would reappear because a direct product form of the time evolution operator is no longer applicable.We will tackle this in a future work.
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 B helped us to clarify the cause of the 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 ≤ x ≤ L x and 0 ≤ y ≤ L y .The space is discretized by numbers of grids (M x , M y ), and its boundaries are periodic.We set ℏ = c = 1 in the following formulations for a notational simplicity.

A. Single-photon system
We begin with a theoretical framework of the single-photon simulation as introduced by [16] (also see the Supplemental Information for the details).The system consists of N A two-level atoms (we write "atom" in the following for the simplicity) and one photon that has a wave-number vector k.The total Hamiltonian is ĥ = ĥ0 + ĥI , and â is an annihilation operator.The frequencies ω k and ω j indicate eigenenergies of the corresponding photon mode and atom, respectively.The Hamiltonian ĥI represents a dipole-dipole interaction between the photon and atoms, where L = L x L y .The parameter D j sets the strength of the dipole-dipole interaction, and the vector r j the position of the jth atom.This interaction makes the constituent atoms play a role of linear optical objects such as mirror, beamsplitter, and scatterer by controlling the parameters ω j and D j .
A quantum state of this system is restricted so that the total number of excitations is one, namely, where c(t, k) and c j (t) are the complex amplitudes satisfying k |c(t, k)| 2 + j |c j (t)| 2 = 1.The atoms are deployed at the grid lattice points to form necessary optical objects.The indices 1 k and 1 j in the kets denote one photon of k mode and one excitation of a j-th atom, respectively.Namely, operations of âk to the bases are and those of âj are the same manner.In addition, due to the restriction of the total number of excitations to one.The Schrödinger equation and the solution of the time evolution are written by i ∂ |ϕ(t)⟩ ∂t = ĥ |ϕ(t)⟩ ( 13)

B. 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 The Hamiltonians ĥ0 and ĥI in Eqs. 10 and 11 can be extended by just replacing âk → âk,p , âj → âj,p .Note that these Hamiltonians do not influence the polarization degrees of freedom.The interaction Hamiltonian ĥI becomes ĥI = p,j,k 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 ω 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 π 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 The basis transformation of the single-photon system is Then, Θ is set as Θ = θ rot /2 + θ pol depending on a rotation angle θ rot given by the polarization rotator, where θ pol denotes an angle of the linear polarization of the input photon against a direction of the horizontal polarization.ω j and D j,S are set at appropriate values to add π 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 C. 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 θ rot = π/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 θ rot , as shown in Fig. 5(b).The calculated probabilities are well-fitted by the behavior of an ideal polarization rotator sin 2 θ rot .These results confirm that the quantum state including the polarization and the relevant optical objects used in Sec.II D work as intended.

C. Time evolution by Suzuki-Trotter decomposition
The previous study [16] rewrites the Schrö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 [43][44][45][46][47][48][49].The time evolution operator in Eq. 14 can be approximated by the Suzuki-Trotter decomposition, as where δt is a time step of the simulation.Each Hamiltonian is diagonalized to compute the time evolution operator.The ĥ0 is already diagonal in the wave-number basis k.To diagonalize ĥI , it is numerically convenient to transform the interaction Hamitonian ĥI from the basis k into the position basis r ĥI = p,j (W j,p â † j,p ârj,p + W * j,p âj,p â † rj ,p ), (20) because by using a Fourier transformation the term in Eq. 11 becomes where W j,p = −iD j,p M ω j / √ 2L.M = M x M y is the number of grids.The operation of âr,p is As seen in Eq. 20, ĥI exchanges the photon and the excitation of the atom located at the same position.This compact representation thanks to the r basis is advantageous in the numerical diagonalization.In fact, multiplying |ϕ⟩ in Eq. 15 gives .
By using matrix representation, this is equivalent to where V = 2 −1/2 1 −i 1 i is a unitary matrix and we used the fact that W * j,p = −W j,p .Therefore, the time evolution can be operated, as .
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 ⟨ϕ(t)|ϕ(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 [46,47].The above observations suggest that the time step δt = 0.1 is sufficiently small to justify the approximation of Eq. 19.Although an error of a time step in the second-order Suzuki-Trotter decomposition is an order of δt 3 which is larger than that of the fourth-order Runge-Kutta method (δt 5 ), its property of the energy conservation hinders the accumulation of the error in the long-time dynamics [50], making the method preferable.To keep the error of the Runge-Kutta method comparable to that of the Suzuki-Trotter method, a significantly smaller time step than that used in the Suzuki-Trotter method is required (see Fig. 6(c)).Such a small time step leads to the longer computational time for the Runge-Kutta method.The time evolution by the Suzuki-Trotter decomposition significantly improves stability of the calculation, which gives a reliable quantitative evaluation by the QOD simulation.

D. Multi-photon system
For notational simplicity, we omit the polarization index p, since the inclusion of p can be done by replacing r → r, p and j → j, p.When naively written, an N -photon state without the polarization is A straightforward updating of the coefficients c(t, r 1 , r 2 , • • • r N ) would require to represent the time evolution operator as an O(M N × M N ) matrix, as the dimension of the Hilbert space of the N -photon states in M grids is O(M N ).This can be reduced if the time evolution operator can be diagonalized but still requires manipulation of O(M N ) elements.
To expedite the calculation, we exploit the fact that N -partite states such as |Φ(t)⟩ can be expressed with a fewer number of terms with a suitable choice of local states [51].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 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 ĥI within each single-photon system.Our N -photon state becomes where ξ 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 and the time evolution of the multi-photon system is written as, Note that the dimension of each exp(−i ĥt) that needs to be compute is M , whereas the naive approach treats the time evolution operator Û (t) of size M N .This reduction of the dimensionality dramatically improves speed of the time evolution and makes the multi-photon simulation feasible.

E. 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 where σ is a width of the Gaussian.The vectors k and r are parameters to decide the initial velocity and position, respectively.For the simulations with the polarization, the initial state is given by

F. Computational details
Table I summarizes the parameters list that are used in the calculation shown in the main text.All the simulations are performed by the time step δt = 0.1 and Gaussian width σ = 2.0.Because we used the arbitrary unit by using ℏ = c = 1, here we evaluate the spatio temporal scale in the real units.To recover the units, at a typical optical wavelength of 500 nm, the system size of the Mach-Zehnder interferometer (see Fig. 1) is ∼ 0.025 mm and simulation time t = 46 is corresponding to 0.12 ps.Naively, for the lab scale simulation with a photon at an optical wavelength, we need to set a 10 5 times larger space size.To reduce the computational cost, Krylov subspace techniques [52,53] and quantics tensor trains [54] might be applied.

FIG. 1 .
FIG. 1. (a)Snapshots of the simulation of the MZ interference when imposed the phase shift φ = π.The contour shows the photon number density.The gray lines and filled rectangles correspond to the beamsplitters (BS1, BS2) and mirrors (M1, M2), respectively.The unfilled rectangle corresponds to the phase shifter (PS), which adds the phase-shift exp(iφ) to the photon.(b) Plot of the probability P right (φ) obtained by the Suzuki-Trotter (blue) and the Runge-Kutta (orange) method.dt is the time step used in each method.The black dotted curve shows the standard theoretical prediction.The horizontal gray dotted line is drawn at P right = 1 as a guide for the eye.

2 FIG. 2 .
FIG. 2. (a) The time evolution of the photon number density with the scattering obstacle located at the center of space.Two typical widths of the detector region are illustrated by the dashed and solid lines.The width between the two dashed lines is narrower and that between the two solid lines is wider.(b) The upper panel shows the error rate as a function of the width of the detector window.The results are obtained by changing ∆x but fixing ∆y = 0.The gray dotted line shows the result without a scatterer.The vertical lines correspond to the typical widths of the detector region.The lower panels show the final state of ∆x = −7.4 and 11.0.(c) The same as (b) but the plots denote the different values of ∆y.We set ∆x = 7.4.The inset shows detailed plots around the region indicated by the black arrow.The lower panels show the final states for ∆y = 0.0, 0.6, and 1.2.
FIG. 3. (a) Schematic of a setup for HOM interference.The beamsplitter is located at the center of the space tilted at 45 degrees.The initial position of the ξ photon in the x direction is shifted by ∆x.(b) Time evolution of the probability distribution of the two photons being found at the same position, i.e., bunching probability distribution.(c) Simulation result and analytical solution in Eq. 4 of the HOM interference.The vertical axis indicates the coincidence probability, namely, the probability of detecting a photon simultaneously in the two separate output directions (up and right).The inset panels show the bunching probability distribution within the white rectangle in (b).
leads to E(θ a , θ b ) = cos 2θ ab .When setting the angles at θ ab = θ a ′ b = θ a ′ b ′ = θ (see also the inset of Fig. 4(c) ), and hence θ ab ′ = 3θ, Eq. 5 becomes S(θ) = 3 cos 2θ − cos 6θ, (7) in the case of the maximally entangled state.Equation 7 breaks the inequality −2 ≤ S(θ) ≤ 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 θ a = 0 and θ b = π/2.The initial state is given by |Φ(t = 0)⟩ = (|H⟩ |H⟩ + |V ⟩ |V ⟩)/ √ 2, where |H⟩ |H⟩ (|V ⟩ |V ⟩) denotes the two photons polarized in the horizontal (vertical) direction.Figure 4(b) shows the case of a product state |Φ(t = 0)⟩ = |H⟩ |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 initial state of the two photons is |Φ⟩ = |H⟩|H⟩+c|V ⟩|V ⟩ √ 1+c 2 FIG. 4. (a) Time evolution of two photons in the Bell-CHSH experiment, in the case of the maximally entangled state.The left and right panel show the photon number density of the horizontal and the vertical polarization component, respectively.The filled and unfilled rectangles represent the polarization beamsplitter (PBS) and the polarization rotator, respectively.The angles of the two polarization rotators are θa = 0, θ b = π/2.(b) The same as (a) but the input photon is a product state.(c) The correlation value S(θ) in Eq. 5.The horizontal dashed gray lines indicate the limit of the local realism.The inset explains a configuration of the angles set in this simulation: θa = 0 and θ ab = θ a ′ b = θ a ′ b ′ = θ.The upper panels show the correlation densities S(θ, r) of the maximally entangled and product states at θ = 0, π/8, π/4.

FIG. 5 .
FIG. 5. (a)The result of the simulation with the polarization rotator and the polarization beamsplitter (PBS).We set the rotation angle of the polarization rotator as θrot = π/4.The polarization rotator consists of 2048 atoms composed with 16 atomic layers, and its parameters are set as Dj,S = 0.56 and ωj = 1.2.PBS consists of 1185 atoms composed with 8 atomic layers, and its parameters are set as Dj,V = 0.56 and ωj = 5.0.The cavity lengths are Lx = 20π and Ly = 10π, and the grids are cut by 512 × 256.The Gaussian width of the photon is σ = 2.0.The initial position and wave number of the photon are r = (2.0,Ly/2) and k = (10.0,0).(b) Plot of the probability PV (θrot) which is the probability of the photon polarized in the vertical direction at the end of simulation.The solid line is a theoretical prediction assuming an idealized PBS.

FIG. 6 .
FIG. 6.(a) Schematics of a test system.The space lengths are Lx = Ly = 10π and the grids are 256 × 256.The boundary is periodic.The mirror consists of 1,584 atoms (Dj = 0.5, ωj = 2.5), and the width of mirror is composed of 8 atomic layers.The mirror denoted by "M" is located at the center of space with tilted at 45 degrees.The initial position and wave number of the injected photon are r = (5.0,Ly/2) and k = (5.0,0).The width is σ = 1.0.Trajectories of (b) probability of the quantum state defined by ⟨ϕ(t)|ϕ(t)⟩ and (c) total energy ⟨ϕ(t)| ĥ |ϕ(t)⟩.The present method using Suzuki-Trotter decomposition is compared with the method used in[16], namely, the fourth-order Runge-Kutta in the interaction picture.