A multistate displaced oscillator system strongly coupled to a heat bath is considered a model of an electron transfer (ET) reaction system. By performing canonical transformation, the model can be reduced to the multistate system coupled to the Brownian heat bath defined by a non-ohmic spectral distribution. For this system, we have derived the hierarchy equations of motion for a reduced density operator that can deal with any strength of the system bath coupling at any temperature. The present formalism is an extension of the hierarchy formalism for a two-state ET system introduced by Tanimura and Mukamel into a low temperature and multistate system. Its ability to handle a multistate system allows us to study a variety of problems in ET and nonlinear optical spectroscopy. To demonstrate the formalism, the time-dependent ET reaction rates for a three-state system are calculated for different energy gaps. Electron transfer (ET) processes play an important role in many fields in physics, chemistry, and biology.1) Most ET processes occur in condensed phases where the surrounding molecules provide the energetic fluctuations and dissipations needed in the reactions. In the case of ET in a polar solvent, the interaction energy that consists of the solvent-solute interaction energy and solvent dipole-dipole energy is used as the ET reaction coordinate.2-4) The ET reaction is then characterized by parabolic free-energy surfaces for the reactant and product states as a function of the reaction coordinate. Although ET rates are commonly studied using free-energy surfaces, the ultrafast dynamics of ET processes is investigated using the potential energy surfaces as a function of molecular coordinates.5,6) This potential energy description is similar to that used in nonlinear ultrafast spectroscopy. 7,8) By adopting this description, one may study ET processes by nonlinear optical measurements. 9,10) A variety of approaches that have been developed to study not only ET processes but also nonlinear optical responses are used to investigate ET dynamics. [5][6][7][8][9][10][11][12][13][14] Unlike most of the above-mentioned approaches based on the perturbative treatment of nonadiabatic transitions, the reduced equation of motion approach that describes the dynamics of density matrix of an ET system coupled to the environment can handle any strength of nonadiabatic coupling. This approach is successful in a classical case characterized by free-energy surfaces, 5,6,15) but it has to employ crucial assumptions such as rotating wave approximation and perturbative system-bath interaction in a quantum case. This strongly limits its applicability. Meanwhile, it was found that the hierarchy equation approach is a powerful means of describing a system strongly coupled to a bath without using rotating wave approximation and white-noise (Born or Van Hove) approximation. [16][17][18][19][20][21] This approach was introduced to investigate the connection between the phenomenological stochastic Liouville equation theory and dynamical Hami...