  We use the "generalized hierarchical equation of motion" proposed in Paper I
to study decoherence in a system coupled to a spin bath. The present
methodology allows a systematic incorporation of higher order anharmonic
effects of the bath in dynamical calculations. We investigate the leading order
corrections to the linear response approximations for spin bath models. Two
types of spin-based environments are considered: (1) a bath of spins
discretized from a continuous spectral density and (2) a bath of physical spins
such as nuclear or electron spins. The main difference resides with how the
bath frequency and the system-bath coupling parameters are chosen to represent
an environment. When discretized from a continuous spectral density, the
system-bath coupling typically scales as $\sim 1/\sqrt{N_B}$ where $N_B$ is the
number of bath spins. This scaling suppresses the non-Gaussian characteristics
of the spin bath and justify the linear response approximations in the
thermodynamic limit. For the physical spin bath models, system-bath couplings
are directly deduced from spin-spin interactions with no reason to obey the
$1/\sqrt{N_B}$ scaling. It is not always possible to justify the linear
response approximations. Furthermore, if the spin-spin Hamiltonian and/or the
bath parameters are highly symmetrical, these additional constraints generate
non-Markovian and persistent dynamics that is beyond the linear response
treatments.
