Quantum Electromagnetic Finite-Difference Time-Domain Solver

: We employ another approach to quantize electromagnetic ﬁelds in the coordinate space, instead of the mode (or Fourier) space, such that local features of photons can be efﬁciently, physically, and more intuitively described. To do this, coordinate-ladder operators are deﬁned from mode-ladder operators via the unitary transformation of systems involved in arbitrary inhomogeneous dielectric media. Then, one can expand electromagnetic ﬁeld operators through the coordinate-ladder operators weighted by non-orthogonal and spatially-localized bases, which are propagators of initial quantum electromagnetic (complex-valued) ﬁeld operators. Here, we call them QEM-CV-propagators . However, there are no general closed form solutions available for them. This inspires us to develop a quantum ﬁnite-difference time-domain (Q-FDTD) scheme to numerically time evolve QEM-CV-propagators. In order to check the validity of the proposed Q-FDTD scheme, we perform computer simulations to observe the Hong-Ou-Mandel effect resulting from the destructive interference of two photons in a 50/50 quantum beam splitter.

The fundamental assumption in the above is that photon is treated as the smallest energy lump, which is carried in the form of EM fields, having a definite value ofhω while being spatially indeterministic (for example, the expectation value of the energy density for a monochromatic photon in the vacuum is uniform over all space due to the characteristics of a plane wave). As a consequence, these formulations correctly account for the anomalous observations including black-body radiation and photoelectric effects; however, it is possible, though, inefficient and, more importantly, physically less intuitive to characterize local features of photons observed in many quantum optics experiments. Basically, this difficulty comes from the wave-particle duality principle, which has been one of the bizarre properties of quantum physics. In a word, canonical quantization is not a universal formulation to describe the behaviors of photons; but one way to characterize photons from the viewpoint of waves, that is, deterministic momentum and indeterministic position.
In this report, we (1) employ another approach to quantize EM fields in the coordinate space [1,26], instead of the conventional mode space, and (2) develop a quantum finite-difference time-domain scheme to numerically time evolve EM field operators such that it is best suited for capturing the local features of photons. To do this, we define coordinate-ladder (CL) operators from mode-ladder (ML) operators via the unitary transformation of the systems involved in arbitrary inhomogeneous dielectric media. Then, EM field operators are expanded by the CL operators weighted by non-orthogonal and spatially-localized bases, which can be interpreted as propagators of initial quantum electromagnetic (complex-valued) field operators (QEM-CV-propagators); however, unlike the classical propagator [27,28] there are no closed form solutions available for general QEM-CV-propagators. This inspires us to develop a quantum finite-difference time-domain (Q-FDTD) scheme which can numerically time evolve QEM-CV-propagators. In order to test the validity of the proposed Q-FDTD scheme, we perform computer simulations to observe the Hong-Ou-Mandel effect [29], which is widely utilized in many quantum optics experiments [30,31], resulting from the destructive interference of two photons in a 50/50 quantum beam splitter.
It is to be mentioned that there are some previous works adopting the coordinate space to describe coupled-resonator optical waveguide based on the tight-binding model or nearest hopping terms [32,33]. But, the present work provides the specific computational framework by fully taking into account the long-range characteristics of the hopping terms to accurately describe the quantum information propagation in arbitrary non-medium-dispersive and lossless EM environments.
It should be also emphasized that the classical FDTD method [34,35] is one of the most popular and powerful numerical solvers widely used in a variety of research areas. This is due to the simplicity w.r.t. both formulation and implementation with a high reliability on simulation results. Hence, the proposed Q-FDTD scheme would be a useful and accessible tool for theoretical/experimental scientists/laboratories in quantum optics.

Numerical Canonical Quantization
Many previous theoretical works [2,17,18] showed that, in principle, the concept of the standard canonical quantization in the vacuum can be extended to that of inhomogeneous and anisotropic lossless media once normal modes of the systems are properly found even though their closed form solutions are often restricted due to the geometrical complexity. As an effective solution to such difficulty, for the first time, we recently proposed in Reference [36] the so-called numerical canonical quantization in which normal modes are numerically obtained by solving the Helmholtz wave equations for arbitrary inhomogeneous dielectric media through computational electromagnetic (CEM) tools such as finite-difference or -element methods.
Suppose that an arbitrary inhomogeneous dielectric object is present in the vacuum. We assume that a finite-sized periodic vacuum box includes the inhomogeneous dielectric object, as illustrated in Figure 1. To extract traveling-wave normal modes, one should employ Bloch-periodic boundary conditions instead of periodic boundary conditions that only permit the existence of standing-wave normal modes. Thus, the local dynamics of photons can be correctly captured within a proper time window. Consider a given N x number of grid points evenly-spaced by ∆x, that is, X = {x m : x m = m∆x for m ∈ [1, 2, · · · , N x ]} that discretizes the primitive cell of the periodic media [37], as depicted in Figure 1. Using a countably-finite set (N k = N x number) of numerical normal modes, one can quantize a vector potential operator [36] aŝ where k p , ω p , andâ k p are Bloch wavenumber, eigenfrequency, and mode-annihilation operator for m-th numerical normal mode, denoted by Φ p (x), respectively. The corresponding electric field operator becomes (based on Φ = 0 gauge or generalized radiation or transverse gauge [38]) For convenience, we employ the linear algebra notation to represent the positive frequency component of (2) asÊ where in what follows, a boldface letter stands for a column vector and a boldface letter with a bar represents a matrix and a dimensionless quantity is measured in the unit of [one]. The numerical normal modes satisfy the orthonormal condition [2,39] as where I is the identity matrix, M is a mass matrix, incorporating medium and metric properties, and C is the Cholesky decomposition of M, that is, M = C † · C; hence, C · Φ becomes a unitary matrix. It should be mentioned that both M and C are diagonal matrices when using a finite-difference method and their elements are given by M m,m = m ∆x and C m,m = √ m ∆x where m = (x m ). Consequently, the Hamiltonian operator can be diagonalized aŝ where D (e) p,p =hω p and E 0 denotes the zero-point energy. One can think of numerical normal modes as uncoupled harmonic oscillators. Therefore, there is no energy coupling between different uncoupled harmonic oscillators; consequently, the Hamiltonian operator in the mode space becomes diagonalizable.

Relation between Mode-and Coordinate-Ladder Operators
One can define coordinate-ladder (CL) operators, denoted byb orb † , from mode-ladder (ML) operators via the unitary transformation or matrix shown in (4) where b n =b x n for n-th grid point. Inversely, one can also writê Note that CL operators still hold the commutator relation as and the lowering and raising process aŝ Note that, in what follows, the index n is reserved to indicate n-th grid point (i.e., x n = n∆x) for CL operators.

Hamiltonian Operator in the Coordinate Space
By substituting (8) and (9) into (5), one finds that the Hamiltonian operator in the coordinate space becomes non-diagonalizable aŝ where J n,n is called hopping terms taking the form of The hopping terms basically dictate the extent of energy coupling between different coupled harmonic oscillators. Unlike to the nearest hopping terms often used to model electron tunneling, Ising model, or coupled-resonator optical waveguides [32,33,40,41], the present hopping terms have the relatively long-range characteristics. These long-range hopping terms of the Hamiltonian operator play a key role in having the equivalence of solutions from Schrödinger-like quantum state equations (parabolic) and quantum Maxwell's equations (hyperbolic). In other words, over-truncation of hopping terms for arbitrary EM environments may cause the significant amount of artificial dispersion errors to the quantum information propagation.

Electric Field Operator in the Coordinate Space
Substituting (8) into (3), one can rewrite (3) w.r.t. CL operators aŝ The new basis, denoted by K (+) n (x m , t), can be thought as a propagator for n-th CL operator (QEM-CV-propagator). For more intuitive understanding, consider n-th QEM-CV-propagator for the vacuum which can be simplified as for m = 1, 2, · · · , N x . Except for the existence of hω p /2N k , QEM-CV-propagators look quite similar to the classical counterpart (CEM-CV-propagator) [27,28] of which continuum expression can be defined by In other words, the resultant classical (real-valued) fields are given by the sum of two convolutions which are equivalent to (1) the delta function sifting property and (2) the Hilbert transform w.r.t. x as where Ω denotes the 1-D problem domain and b x is a complex-valued intial field scalar value at x = x . Note that the Hilbert transform plays a key role to model one-way wave propagations, provided that a set of b x describes the phasor form of one-way traveling waves. This formulation is particularly feasible to deal with a wave equation (for a single-species unknown variable) rather than Maxwell's curl equations. It should be noted that previous classical 1-D Maxwell-curl-propagators [27,28,42] consist of the time-derivative of retard and advanced time-domain Green's functions and are supposed to carry real-valued initial electric and magnetic field values as where E x and H x are real-valued initial electric and magnetic field values at x = x and θ(t) is a unit time-step function. As seen in (22), the underlying principle of one-way wave propagations is to cancel one of the two delta functions by controlling the initial electric or magnetic field values.
In the quantum realm, positive-and negative-frequency components of QEM-CV-propagators become bases of annihilation and creation operators of which actions on quantum states have their own physical meaning such as annihilation and creation of photons, respectively. Thus, one has to treat positive-and negative-frequency propagators individually for quantum electromagnetics. The discrete representation of (20) can be written as for m = 1, 2, · · · , N x . One can think of QEM-CV-propagators as the convolution between the inverse Fourier transform of (band-limited) hω p /2N k and CEM-CV-propagators (see Figure 2). Furthermore, employing the concept of fractional derivatives, one can associate CEM-CV-propagator with QEM-CV-propagator as where D 0.5 t denotes the half derivative w.r.t. time. As implied in (24), the vector space of the expectation value of the energy in the quantum regime becomes the subset of that of the classical EM energy. This results from the restriction to field amplitudes of quantum EM fields so as to enforce the energy of N number of photons in a monochromatic wave to be Nhω whereas there is no restriction in classical field amplitudes, hence, classical EM energy can take any positive real value. Figure 2 depicts an example of the temporal behaviors of a CEM-/QEM-CV-propagators. It is worth noting that the factor 1/ω p in (23) prevents to arrive at the correspondence principle between classical and quantum propagators even for the vacuum, that is, (23) and (20). This is a consequence of non-covariant normalization of the plane wave basis [43]. The invariant normalization approach [43,44] can be used to eliminate the factor. We will investigate the invariant normalization approach for our future work.
It should be emphasized that the concept of propagators for classical EM fields has been already discussed in References [27,28]. Their propagators are designed to carry real-valued initial field (scalar) values. However, propagators shown here can carry arbitrary complex-valued initial fields. Thus, it is suited to analyze the quantum nature of EM fields wherein probability amplitudes of quantum state can be arbitrary complex values as long as the normalization condition is satisfied. Then, let us initialize a single photon spatially-localized at x n whose initial quantum state is given by |ψ =b † x n |0 . The expectation value of the electric energy operator, denoted byĤ e , for the single photon can be evaluated as which corresponds to the half of the average value of all possible eigenenergies of the systems. Hence, we can say that QEM-CV-propagators describe how EM field operators are time evolving when a single photon is initialized at a certain grid, having its energy expectation value as 1 N k ∑ N k p=1h ω p . In other words, a single photon described in the coordinate space becomes more (and less) deterministic w.r.t. position (and energy). Note that the opposite situation happens in the mode space.

Quantum Finite-Difference Time-Domain Scheme
We have seen in the previous section that electric field operators could be expanded by CL operators weighted by QEM-CV-propagators which are given by the summation of numerical normal modes. However, there are no closed form solutions for general QEM-CV-propagators due to the energy scaling factor, viz., hω p /2. Hence, we develop a quantum finite-difference time-domain scheme to numerically time evolve QEM-CV-propagators.
Here, we discretize the continuum of (16) and (23), that is, x m is the variable to be discretized in the present Q-FDTD scheme whereas x n denotes the position of initializing n−th QEM-CV propagator.
In other words, one can think of x m → x and x n → x . The Hamilton's equations of motion (or wave equations) for the positive frequency component of electric field operators can be written as Substituting (16) into (26), one can have Since coordinate-annihilation operators do not rely on space and time, the equation (27) can be further simplified by It is to be noted that there should be no coupling between different QEM-CV-propagators; hence, one can independently find the time evolution of QEM-CV-propagators as for n = 1, 2, · · · , N x . Applying a finite difference method to (29) w.r.t. space and time, one can derive the discrete counterpart of (29), called a quantum finite-difference time-domain (Q-FDTD) scheme, for n-th QEM-CV-propagator as Note that the above Q-FDTD scheme is embarrassingly-parallel due to the independence of index n while marching-on in time. Finally, numerical solutions of the electric field operator can be found aŝ

Initial Quantum States for Few Photons
In the mode space, an initial quantum state for a single photon spatially localized in the form of Gaussian wavepacket is given by wherew p is a probability amplitude forâ † k p |0 = |1 p , which corresponds to a spectral amplitude of the Gaussian wavepacket. Substituting (9) into (33), one can rewrite the initial quantum state in the coordinate space as 2 e ik w (x n −x w ) and x w , σ w , and k w are localization center position, localization degree, and carrier wavenumber of the Gaussian wavepacket, respectively, while satisfying the normalization condition w † · w = 1. Note that one can also model wavepackets using Lorentzian or sinc functions.
An initial composite quantum state for non-entangled two (indistinguishable) photons can be given by the tensor product of each photon's quantum state as It is to be noted that (35) is invariant under the swapping process, viz., b , since photons are assumed to be indistinguishable.

Initial Conditions of Quantum Finite-Difference Time-Domain Scheme
To run the Q-FDTD scheme, one needs initial values of QEM-CV-propagators and their time derivatives. Again, due to the non-existence of the closed form solutions for general QEM-CV-propagators, in order to get the initial values, one has to know numerical normal modes by solving the Helmholtz wave equations before running the Q-FDTD scheme, but, this becomes redundant.
Here, we show an alternative way of the initialization of the Q-FDTD scheme devoid of solving for the normal modes. First, assume that an initialized photon is spatially localized around x w in the vacuum. This allows one to only take into account some of QEM-CV-propagators that significantly affect the quantum information propagation. For example, one can select QEM-CV-propagators defined over grid point index n ∈ Y where Y = x n ∈ X : |[w] n /max (w)| > 1/e 2 . Second, assume that the photon initialized in the vacuum is spatially far away from dielectric objects. In other words, the photon cannot recognize the presence of dielectric objects before touching them. As a result, QEM-CV-propagators can be approximated by analytic plane wave solutions as for n ∈ Y. Thus, the Q-FDTD scheme can be initialized by letting n (x m , t) ∂t

Numerical Simulations of Quantum Beam Splitter
The Hong-Ou-Mandel effect [29] refers to a two-photon destructive interference effect in a 50/50 quantum beam splitter [45] and widely used in quantum optics experiments. Specifically, when a pair of photons is impinging on the beam splitter with the perfect temporal overlap, that is, τ = 0, the two photons will always exit through a same output port while being bunched, having 50:50 chance of exiting from either one or the other output port. This is often measured through the degree of intensity coherence, called second order correlation, denoted as g (2) , originally related to the Hanbury Brown and Twiss (HBT) effect [46]. The quantum version of the second order correlation [47] is defined by (2) . (40) In our computer simulation, we set x 1 = x l and x 2 = x r , as illustrated in Figure 3, and τ = δx 0 /c [s]. We used same simulation parameters used in Reference [36]. Also, One can evaluate (40) in the same way explained in Reference [36], which is the crucial part of creating quantum effects.   (36) to (39). There are excellent agreement in all cases, thus, the Q-FDTD scheme with approximate initialization is correctly validated. Furthermore, as expected, there are the HOM dips in all cases when τ = 0 since the destructive interference of the pair of photons occurs.

Conclusions
We have employed another approach to quantize electromagnetic (EM) fields in the coordinate space, instead of the mode space. In the coordinate space, EM field operators were expanded by coordinate-ladder operators weighted by non-orthogonal and spatially-localized bases, which were propagators of quantum electromagnetic initial field operators (QEM-CV-propagator). Due to the property of QEM-CV-propagators, the present formulation is suited to describe local features of photons. Since there were no closed form solutions available for them in arbitrary EM environments, we developed a quantum finite-difference time-domain method to numerically time evolve QEM-CV-propagators. In order to check the validity of the proposed quantum finite-difference time-domain scheme, we performed computer simulations to observe the two-photon indistinguishability in quantum beam splitters, known as the Hong-Ou-Mandel effect.