  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.
