For a detailed description of the method, we refer readers to Ref. 54 and the literature therein. Here, we only briefly introduce the main ideas of the method in twin space. Every component in Eq. ( ) can be interpreted as a dissipative bosonic mode with a complex-valued frequency γ α p and coupling strength to the system η α p . A specific configuration n α = (n α0 ,n α1 ,n α2 ,⋯n αP ) specifies the state with n α p phonons in the pth bosonic mode (p runs from 0 to P) of bath α. For a super index n = (⋯,n α ,⋯), one can introduce an auxiliary density operator (ADO), ρ n (t). This set of auxiliary density operators is closed with respect to the time-derivative operation and yields a hierarchical equation of motion