  An extremely useful evolution equation that allows systematically calculating
the two-time correlation functions (CF's) of system operators for non-Markovian
open (dissipative) quantum systems is derived. The derivation is based on
perturbative quantum master equation approach, so non-Markovian open quantum
system models that are not exactly solvable can use our derived evolution
equation to easily obtain their two-time CF's of system operators, valid to
second order in the system-environment interaction. Since the form and nature
of the Hamiltonian are not specified in our derived evolution equation, our
evolution equation is applicable for bosonic and/or fermionic environments and
can be applied to a wide range of system-environment models with any factorized
(separable) system-environment initial states (pure or mixed). When applied to
a general model of a system coupled to a finite-temperature bosonic environment
with a system coupling operator L in the system-environment interaction
Hamiltonian, the resultant evolution equation is valid for both L = L^+ and L
\neq L^+ cases, in contrast to those evolution equations valid only for L = L^+
case in the literature. The derived equation that generalizes the quantum
regression theorem (QRT) to the non-Markovian case will have broad applications
in many different branches of physics. We then give conditions on which the QRT
holds in the weak system-environment coupling case, and apply the derived
evolution equation to a problem of a two-level system (atom) coupled to a
finite-temperature bosonic environment (electromagnetic fields) with L \neq
L^+.
