  We extend a standard stochastic theory to study open quantum systems coupled
to generic quantum environments including the three fundamental classes of
noninteracting particles: bosons, fermions and spins. In this unified
stochastic approach, the generalized stochastic Liouville equation (SLE)
formally captures the exact quantum dissipations when noise variables with
appropriate statistics for different bath models are applied. Anharmonic
effects of a non-Gaussian bath are precisely encoded in the bath multi-time
correlation functions that noise variables have to satisfy. Staring from the
SLE, we devise a family of generalized hierarchical equations by averaging out
the noise variables and expand bath multi-time correlation functions in a
complete basis of orthonormal functions. The general hiearchical equations
constitute systems of linear equations that provide numerically exact
simulations of quantum dynamics. For bosonic bath models, our general
hierarchical equation of motion reduces exactly to an extended version of
hierarchical equation of motion which allows efficient simulation for arbitrary
spectral densities and temperature regimes. Similar efficiency and exibility
can be achieved for the fermionic bath models within our formalism. The spin
bath models can be simulated with two complementary approaches in the presetn
formalism. (I) They can be viewed as an example of non-Gaussian bath models and
be directly handled with the general hierarchical equation approach given their
multi-time correlation functions. (II) Alterantively, each bath spin can be
first mapped onto a pair of fermions and be treated as fermionic environments
within the present formalism.
