  The recent experimental discoveries about excitation energy transfer (EET) in
light harvesting antenna (LHA) attract a lot of interest. As an open
non-equilibrium quantum system, the EET demands more rigorous theoretical
framework to understand the interaction between system and environment and
therein the evolution of reduced density matrix. A phonon is often used to
model the fluctuating environment and convolutes the reduced quantum system
temporarily. In this paper, we propose a novel way to construct complex-valued
Gaussian processes to describe thermal quantum phonon bath exactly by
converting the convolution of influence functional into the time correlation of
complex Gaussian random field. Based on the construction, we propose a rigorous
and efficient computational method, the covariance decomposition (CD) and
conditional propagation scheme, to simulate the temporarily entangled reduced
system.
  The new method allows us to study the non-Markovian effect without
perturbation under the influence of different spectral densities of the linear
system-phonon coupling coefficients. Its application in the study of EET in the
Fenna-Matthews-Olson (FMO) model Hamiltonian under four different spectral
densities is discussed. Since the scaling of our algorithm is linear due to its
Monte Carlo nature, the future application of the method for large LHA systems
is attractive. In addition, this method can be used to study the effect of
correlated initial condition on the reduced dynamics in the future.
