  We propose a new methodology, called numerical canonical quantization, to
solve quantum Maxwell's equations useful for mathematical modeling of quantum
optics physics, and numerical experiments on arbitrary passive and lossless
quantum-optical systems. It is based on: (1) the macroscopic (phenomenological)
electromagnetic theory on quantum electrodynamics (QED), and (2) concepts
borrowed from computational electromagnetics. It was shown that canonical
quantization in inhomogeneous dielectric media required definite and proper
normal modes. Here, instead of ad-hoc analytic normal modes, we numerically
construct complete and time-reversible normal modes in the form of traveling
waves to diagonalize the Hamiltonian. Specifically, we directly solve the
Helmholtz wave equations for a general linear, reciprocal, isotropic,
non-dispersive, and inhomogeneous dielectric media by using either
finite-element or finite-difference methods. To convert a scattering problem
with infinite number of modes into one with a finite number of modes, we impose
Bloch-periodic boundary conditions. This will sparsely sample the normal modes
with numerical Bloch-Floquet-like normal modes. Subsequent procedure of
numerical canonical quantization is straightforward using linear algebra. We
provide relevant numerical recipes in detail and show an important numerical
example of indistinguishable two-photon interference in quantum beam splitters,
exhibiting Hong-Ou-Mandel effect, which is purely a quantum effect. Also, the
present methodology provides a way of numerically investigating existing or new
macroscopic QED theories. It will eventually allow quantum-optical numerical
experiments of high fidelity to replace many real experiments as in classical
electromagnetics.
