Volume 33 Issue 6
Jan.  2021
Turn off MathJax
Article Contents

Yu-Chen Wang, Yi Zhao. The Hierarchical Stochastic Schrödinger Equations: Theory and Applications[J]. Chinese Journal of Chemical Physics , 2020, 33(6): 653-667. doi: 10.1063/1674-0068/cjcp2009165
Citation: Yu-Chen Wang, Yi Zhao. The Hierarchical Stochastic Schrödinger Equations: Theory and Applications[J]. Chinese Journal of Chemical Physics , 2020, 33(6): 653-667. doi: 10.1063/1674-0068/cjcp2009165

The Hierarchical Stochastic Schrödinger Equations: Theory and Applications

doi: 10.1063/1674-0068/cjcp2009165
More Information
  • Corresponding author: Yi Zhao, E-mail:yizhao@xmu.edu.cn
  • Received Date: 2020-09-18
  • Accepted Date: 2020-10-20
  • Publish Date: 2020-12-27
  • The hierarchical stochastic Schrödinger equations (HSSE) are a kind of numerically exact wavefunction-based approaches suitable for the quantum dynamics simulations in a relatively large system coupled to a bosonic bath. Starting from the influence-functional description of open quantum systems, this review outlines the general theoretical framework of HSSEs and their concrete forms in different situations. The applicability and efficiency of HSSEs are exemplified by the simulations of ultrafast excitation energy transfer processes in large-scale systems.
  • 加载中
  • [1] R. P. Feynman and F. L. Vernon, Ann. Phys. 24, 118 (1963).
    [2] Y. Tanimura, J. Chem. Phys. 153, 020901 (2020).
    [3] Y. Tanimura and R. Kubo, J. Phys. Soc. Jpn. 58, 101 (1989).
    [4] Y. A. Yan, F. Yang, Y. Liu, and J. Shao, Chem. Phys. Lett. 395, 216 (2004).
    [5] R. X. Xu, P. Cui, X. Q. Li, Y. Mo, and Y. Yan, J. Chem. Phys. 122, 041103 (2005).
    [6] A. Ishizaki and Y. Tanimura, J. Phys. Soc. Jpn. 74, 3131 (2005).
    [7] J. Hu, R. X. Xu, and Y. Yan, J. Chem. Phys. 133, 101106 (2010).
    [8] C. Meier and D. J. Tannor, J. Chem. Phys. 111, 3365 (1999).
    [9] Z. Tang, X. Ouyang, Z. Gong, H. Wang, and J. Wu, J. Chem. Phys. 143, 224112 (2015).
    [10] C. Duan, Z. Tang, J. Cao, and J. Wu, Phys. Rev. B 95, 214308 (2017).
    [11] Y. Tanimura, J. Phys. Soc. Jpn. 75, 082001 (2006).
    [12] M. Schröder, M. Schreiber, and U. Kleinekathöfer, J. Chem. Phys. 126, 114102 (2007).
    [13] Q. Shi, L. Chen, G. Nan, R. X. Xu, and Y. Yan, J. Chem. Phys. 130, 084105 (2009).
    [14] H. Liu, L. Zhu, S. Bai, and Q. Shi, J. Chem. Phys. 140, 134106 (2014).
    [15] C. Kreisbeck, T. Kramer, M. Rodriguez, and B. Hein, J. Chem. Theory Comput. 7, 2166 (2011).
    [16] M. Tsuchimoto and Y. Tanimura, J. Chem. Theory Comput. 11, 3859 (2015).
    [17] J. Strümpfer and K. Schulten, J. Chem. Theory Comput. 8, 2808 (2012).
    [18] C. Kreisbeck, T. Kramer, and A. Aspuru-Guzik, J. Chem. Theory Comput. 10, 4045 (2014).
    [19] J. Cao, L. W. Ungar, and G. A. Voth, J. Chem. Phys. 104, 4189 (1996).
    [20] J. T. Stockburger and H. Grabert, Phys. Rev. Lett. 88, 170407 (2002).
    [21] J. T. Stockburger, Chem. Phys. 296, 159 (2004).
    [22] J. Shao, J. Chem. Phys. 120, 5053 (2004).
    [23] X. Chen, J. Cao, and R. J. Silbey, J. Chem. Phys. 138, 224104 (2013).
    [24] Y. A. Yan and J. Shao, Front. Phys. 11, 110309 (2016).
    [25] Y. A. Yan and J. Shao, Phys. Rev. A 97, 042126 (2018).
    [26] Y. Zhou, Y. Yan, and J. Shao, Europhys. Lett. 72, 334 (2005).
    [27] Y. Zhou and J. Shao, J. Chem. Phys. 128, 034106 (2008).
    [28] J. M. Moix and J. Cao, J. Chem. Phys. 139, 134106 (2013).
    [29] L. Zhu, H. Liu, and Q. Shi, New J. Phys. 15, 095020 (2013).
    [30] D. Suess, A. Eisfeld, and W. Strunz, Phys. Rev. Lett. 113, 150403 (2014).
    [31] K. Song, L. Song, and Q. Shi, J. Chem. Phys. 144, 224105 (2016).
    [32] Y. Ke and Y. Zhao, J. Chem. Phys. 145, 024101 (2016).
    [33] R. Hartmann and W. T. Strunz, J. Chem. Theory Comput. 13, 5834 (2017).
    [34] R. Hartmann, M. Werther, F. Grossmann, and W. T. Strunz, J. Chem. Phys. 150, 234105 (2019).
    [35] M. Lian, Y. C. Wang, Y. Ke, and Y. Zhao, J. Chem. Phys. 151, 044115 (2019).
    [36] J. Hubbard, Phys. Rev. Lett. 3, 77 (1959).
    [37] R. Stratonovich, Sov. Phys. Dokl. 2, 416 (1957).
    [38] P. P. Zhang, C. D. B. Bentley, and A. Eisfeld, J. Chem. Phys. 148, 134103 (2018).
    [39] A. T. A. Wood and G. Chan, J. Comput. Graph. Stat. 3, 409 (1994).
    [40] J. M. Moix, J. Ma, and J. Cao, J. Chem. Phys. 142, 094108 (2015).
    [41] X. Zhong and Y. Zhao, J. Chem. Phys. 138, 014111 (2013).
    [42] Y. Ke and Y. Zhao, J. Chem. Phys. 147, 184103 (2017).
    [43] Y. C. Wang, Y. Ke, and Y. Zhao, WIREs Comput. Mol. Sci. 9, e1375 (2019).
    [44] L. Han, V. Chernyak, Y. A. Yan, X. Zheng, and Y. Yan, Phys. Rev. Lett. 123, 050601 (2019).
    [45] L. Han, A. Ullah, Y. A. Yan, X. Zheng, Y. Yan, and V. Chernyak, J. Chem. Phys. 152, 204105 (2020).
    [46] A. Ullah, L. Han, Y. A. Yan, X. Zheng, Y. Yan, and V. Chernyak, J. Chem. Phys. 152, 204106 (2020).
    [47] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. 42, 1698 (1979).
    [48] T. Holstein, Ann. Phys. 8, 325 (1959).
    [49] R. E. Peierls, Quantum Theory of Solids, New York: Oxford University Press, (1996).
    [50] R. W. Munn and R. Silbey, J. Chem. Phys. 83, 1843 (1985).
    [51] H. Fröhlich, H. Pelzer, and S. Zienau, Philos. Mag. 41, 221 (1950).
    [52] H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems, New York: Oxford University Press, (2002).
    [53] Y. Ke and Y. Zhao, Mol. Phys. 116, 813 (2018).
    [54] Y. C. Wang and Y. Zhao, J. Comput. Chem. 40, 1097 (2019).
    [55] Y. Ke and Y. Zhao, J. Chem. Phys. 146, 174105 (2017).
    [56] P. P. Zhang and A. Eisfeld, J. Phys. Chem. Lett. 7, 4488 (2016).
    [57] Y. Ke and Y. Zhao, J. Chem. Phys. 149, 014104 (2018).
    [58] M. Mohseni, Y. Omar, G. S. Engel, and M. B. Plenio, Quantum Effects in Biology, New York: Cambridge University Press, (2014).
    [59] R. E. Blankenship, Molecular Mechanisms of Photosynthesis, John Wiley & Sons, (2014).
    [60] R. Fenna, B. Matthews, J. Olson, and E. Shaw, J. Mol. Biol. 84, 231 (1974).
    [61] R. E. Fenna and B. W. Matthews, Nature 258, 573 (1975).
    [62] G. S. Engel, T. R. Calhoun, E. L. Read, T. K. Ahn, T. Mančal, Y. C. Cheng, R. E. Blankenship, and G. R. Fleming, Nature 446, 782 (2007).
    [63] A. Ishizaki and G. R. Fleming, Proc. Natl. Acad. Sci. USA 106, 17255 (2009).
    [64] G. Panitchayangkoon, D. Hayes, K. A. Fransted, J. R. Caram, E. Harel, J. Wen, R. E. Blankenship, and G. S. Engel, Proc. Natl. Acad. Sci. USA 107, 12766 (2010).
    [65] G. D. Scholes, G. R. Fleming, A. Olaya-Castro, and R. van Grondelle, Nat. Chem. 3, 763 (2011).
    [66] C. Kreisbeck and T. Kramer, J. Phys. Chem. Lett. 3, 2828 (2012).
    [67] V. Tiwari, W. K. Peters, and D. M. Jonas, Proc. Natl. Acad. Sci. USA 110, 1203 (2013).
    [68] A. W. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. F. Huelga, and M. B. Plenio, Nat. Phys. 9, 113 (2013).
    [69] N. Lambert, Y. N. Chen, Y. C. Cheng, C. M. Li, G. Y. Chen, and F. Nori, Nat. Phys. 9, 10 (2013).
    [70] D. M. Wilkins and N. S. Dattani, J. Chem. Theory Comput. 11, 3411 (2015).
    [71] H. G. Duan, V. I. Prokhorenko, R. J. Cogdell, K. Ashraf, A. L. Stevens, M. Thorwart, and R. J. D. Miller, Proc. Natl. Acad. Sci. USA 114, 8493 (2017).
    [72] A. Ben-Shem, F. Frolow, and N. Nelson, FEBS Lett. 564, 274 (2004),
    [73] D. E. Tronrud, J. Wen, L. Gay, and R. E. Blankenship, Photosynth. Res. 100, 79 (2009).
    [74] M. Schmidt am Busch, F. Müh, M. El-Amine Madjet, and T. Renger, J. Phys. Chem. Lett. 2, 93 (2011).
    [75] C. Olbrich, T. L. Jansen, J. Liebers, M. Aghtar, J. Strümpfer, K. Schulten, J. Knoester, and U. Kleinekath öfer, J. Phys. Chem. B 115, 8609 (2011).
    [76] X. Jia, Y. Mei, J. Zhang, and Y. Mo, Sci. Rep. 5, 17096 (2015).
    [77] N. Christensson, H. F. Kauffmann, T. Pullerits, and T. Mančal, J. Phys. Chem. B 116, 7449 (2012).
    [78] A. Chenu, N. Christensson, H. F. Kauffmann, and T. Mančal, Sci. Rep. 3, 2029 (2013).
    [79] L. Z. Sharp and D. Egorova, J. Chem. Phys. 139, 144304 (2013).
    [80] V. Butkus, D. Zigmantas, D. Abramavicius, and L. Valkunas, Chem. Phys. Lett. 587, 93 (2013).
    [81] M. B. Plenio, J. Almeida, and S. F. Huelga, J. Chem. Phys. 139, 235102 (2013).
    [82] V. Perlík, C. Lincoln, F. Šanda, and J. Hauer, J. Phys. Chem. Lett. 5, 404 (2014).
    [83] C. Kreisbeck, T. Kramer, and A. Aspuru-Guzik, J. Phys. Chem. B 117, 9380 (2013).
    [84] H. Lee, Y. C. Cheng, and G. R. Fleming, Science 316, 1462 (2007).
    [85] E. Collini, C. Y. Wong, K. E. Wilk, P. M. G. Curmi, P. Brumer, and G. D. Scholes, Nature 463, 644 (2010).
    [86] F. D. Fuller, J. Pan, A. Gelzinis, V. Butkus, S. S. Senlik, D. E. Wilcox, C. F. Yocum, L. Valkunas, D. Abramavicius, and J. P. Ogilvie, Nat. Chem. 6, 706 (2014).
    [87] P. Nalbach, C. A. Mujica-Martinez, and M. Thorwart, Phys. Rev. E 91, 022706 (2015).
    [88] Y. Sato and B. Doolittle, J. Chem. Phys. 141, 185102 (2014).
    [89] J. M. Womick and A. M. Moran, J. Phys. Chem. B 115, 1347 (2011).
    [90] A. Kolli, E. J. OReilly, G. D. Scholes, and A. Olaya-Castro, J. Chem. Phys. 137, 174109 (2012).
    [91] M. d. Rey, A. W. Chin, S. F. Huelga, and M. B. Plenio, J. Phys. Chem. Lett. 4, 903 (2013).
    [92] A. G. Dijkstra, C. Wang, J. Cao, and G. R. Fleming, J. Phys. Chem. Lett. 6, 627 (2015).
    [93] V. Perlík, J. Seibt, L. J. Cranston, R. J. Cogdell, C. N. Lincoln, J. Savolainen, F. Šanda, T. Mančal, and J. Hauer, J. Chem. Phys. 142, 212434 (2015).
    [94] J. C. Dean, T. Mirkovic, Z. S. D. Toa, D. G. Oblinsky, and G. D. Scholes, Chem 1, 858 (2016).
    [95] E. K. Irish, R. Gómez-Bombarelli, and B. W. Lovett, Phys. Rev. A 90, 012510 (2014).
    [96] C. A. Mujica-Martinez and P. Nalbach, Ann. Phys. (Berlin) 527, 592 (2015).
    [97] V. I. Novoderezhkin, E. Romero, and R. van Grondelle, Phys. Chem. Chem. Phys. 17, 30828 (2015).
    [98] P. Malý, O. J. G. Somsen, V. I. Novoderezhkin, T. Mančal, and R. Van Grondelle, ChemPhysChem 17, 1356 (2016).
    [99] D. Abramavicius and L. Valkunas, Photosynth. Res. 127, 33 (2016).
    [100] J. Iles-Smith, A. G. Dijkstra, N. Lambert, and A. Nazir, J. Chem. Phys 144, 044110 (2016).
    [101] M. H. Lee and A. Troisi, J. Chem. Phys. 146, 075101 (2017).
    [102] Y. Fujihashi, L. Chen, A. Ishizaki, J. Wang, and Y. Zhao, J. Chem. Phys. 146, 044101 (2017).
    [103] Z. Huang, Y. Fujihashi, and Y. Zhao, J. Phys. Chem. Lett. 8, 3306 (2017).
    [104] D. M. Monahan, L. Whaley-Mayda, A. Ishizaki, and G. R. Fleming, J. Chem. Phys. 143, 065101 (2015).
    [105] Y. Fujihashi, G. R. Fleming, and A. Ishizaki, J. Chem. Phys. 142, 212403 (2015).
    [106] D. I. G. Bennett, P. Malý, C. Kreisbeck, R. van Grondelle, and A. Aspuru-Guzik, J. Phys. Chem. Lett. 9, 2665 (2018).
    [107] S. M. Blau, D. I. G. Bennett, C. Kreisbeck, G. D. Scholes, and A. Aspuru-Guzik, Proc. Natl. Acad. Sci. USA 115, E3342 (2018).
    [108] Y. Fujihashi, M. Higashi, and A. Ishizaki, J. Phys. Chem. Lett. 9, 4921 (2018).
    [109] J. P. Paz and W. H. Zurek, Phys. Rev. Lett. 82, 5181 (1999).
    [110] T. Holstein, Ann. Phys. 8, 343 (1959).
  • 加载中
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Figures(4)  / Tables(1)

Article Metrics

Article views(28) PDF downloads(4) Cited by()

Proportional views
Related

The Hierarchical Stochastic Schrödinger Equations: Theory and Applications

doi: 10.1063/1674-0068/cjcp2009165

Abstract: The hierarchical stochastic Schrödinger equations (HSSE) are a kind of numerically exact wavefunction-based approaches suitable for the quantum dynamics simulations in a relatively large system coupled to a bosonic bath. Starting from the influence-functional description of open quantum systems, this review outlines the general theoretical framework of HSSEs and their concrete forms in different situations. The applicability and efficiency of HSSEs are exemplified by the simulations of ultrafast excitation energy transfer processes in large-scale systems.

Yu-Chen Wang, Yi Zhao. The Hierarchical Stochastic Schrödinger Equations: Theory and Applications[J]. Chinese Journal of Chemical Physics , 2020, 33(6): 653-667. doi: 10.1063/1674-0068/cjcp2009165
Citation: Yu-Chen Wang, Yi Zhao. The Hierarchical Stochastic Schrödinger Equations: Theory and Applications[J]. Chinese Journal of Chemical Physics , 2020, 33(6): 653-667. doi: 10.1063/1674-0068/cjcp2009165
  • Over the years, exciting phenomena observed in femosecond spectroscopic experiments have been stimulating growing research interests in the quantum effects emerging in a variety of ultrafast photochemical and photophysical processes. These effects include but not limited to non-Markovianity, entanglement, quantum fluctuation and dissipation, environment-assisted quantum transport, nuclear tunneling, etc. The requirement of accurate theoretical descriptions for deeply understanding them has significantly boosted the development of various numerically exact quantum dynamics methods.

    The path integral is an invaluable tool for theoretical researches in the field of open quantum systems, since within this formalism it is straightforward to take partial trace over the environmental (the so-called bath) degrees of freedom and absorb all of the bath influence into the well-known Feynman-Vernon influence functional [1]. One powerful approach that was originally established within the path-integral formalism is the non-perturbative and non-Markovian hierarchical equations of motion (HEOM) [2]. In this approach, the hierarchical expansion technique is applied to the influence functional and leads to a closed set of differential equations of density matrix with a hierarchical structure. These equations are simultaneously solved, and the zeroth-order equation gives the exact reduced density matrix of the system. Pioneered by Tanimura and Kubo in their seminal work [3], where the high-temperature approximation was adopted, great efforts from different contributors have been devoted to deriving a numerically exact veresion of HEOM [4-6]. Various powerful techniques have also been proposed to optimize its numerical performance, for instance, new decomposition schemes for spectral densities [7, 8] and bath correlation function [9, 10], better truncation schemes [6, 11, 12], highly efficient filtering algorithm [13, 14], high-performance computational platform and efficient parallel algorithm [15-18], etc. For a comprehensive introduction of HEOM, the readers can refer to Ref.[2]. With all of these techniques, HEOM has become a standard approach in the field of open quantum systems nowadays. Nevertheless, despite the great achievement, HEOM is still partially limited by the form of the spectral density function, and the computational efforts grow factorially in terms of the system size.

    Another disparate strategy is to perform complete stochastic unraveling to the influence functional by the introduction of complex Gaussian stochastic fields, which results in the stochastic Liouville-von Neumann equation [19-25]. This approach is not limited to the form of the spectral density function, and can be directly decomposed into a forward and a backward stochastic Schrödinger equations in the case of a pure initial system state. However, it may suffer from severe stochastic convergence problem in long-time simulations or when the system-bath interaction is not weak. Recently, several density-matrix approaches [11, 26-29] that combine the hierarchical expansion and stochastic unraveling techniques have also been pro posed, which inherit both merits of HEOM and the stochastic Liouville-von Neumann equation, in the sense that the resulting hierarchical structure is simpler and the stochastic convergence property is better.

    One less satisfactory feature of density-matrix approaches is that the number of elements that needs to be calculated scales quadratically with the system size. As a counterpart of density-matrix approaches, wavefunction-based methods have also been playing a prominent role in the dynamics of open quantum systems due to the beneficial linear scaling of Hilbert space with respect to the system size. Inspired by the spirits of mixed deterministic-stochastic scheme, Suess et al. [30] applied the hierarchical expansion technique to the functional derivative term of non-Markovian quantum state diffusion and proposed a numerically exact approach named hierarchy of pure states (HOPS), which is a closed set of stochastic Schrödinger equations with hierarchical structure. Whereafter, at about the same time, Song et al. [31] and Ke et al. [32] independently derived two similar approaches. In the new approaches, temperature effects are mostly or entirely accounted for by the stochastic fields such that the resulting equations of motion are more stable at finite temperatures as compared with HOPS. Soon afterwards, with a similar spirit, Hartmann, Strunz and coworkers [33, 34] improved the original HOPS by including temperature effects via a stochastic contribution to the system Hamiltonian. In this review, all of the above approaches refer to the hierarchical stochastic Schrödinger equations (HSSE), since they share very similar equations of motion.

    This review aims at giving a comprehensive overview of HSSEs. A general theoretical framework of HSSEs is given and their concrete forms in different situations are discussed. Some recent applications of HSSEs, including the excitation energy transfer (EET) in the Fenna-Matthews-Olson (FMO) complex and the resonant vibration-assisted EET mechanism, are also introduced.

  • We consider the following Hamiltonian that describes an open quantum system interacting with a bosonic bath,

    Here, $ \hat{H}_ {\rm{S}} $, $ \hat{H}_ {\rm{B}} $, and $ \hat{H}_ {\rm{I}} $ are the system Hamiltonian, the bath Hamiltonian, and the interaction between them. Typically, the system refers to a few degrees of freedom that are primarily interested, for instance, free carriers in the charge transport process in photovoltaic devices and excitons in EET process in photosynthetic systems. The bath Hamiltonian is described by a set of independent vibrational modes, the total number of which can be either finite or infinite (hereafter, we set $ \hbar $ = $ k_ {\rm{B}} $ = 1 for the sake of simplicity),

    where $ \hat{b}_\mu^\dagger $ and $ \hat{b}_\mu $ are the creation and annihilation operators for the $ \mu $-th mode with frequency $ \omega_\mu $. These normal modes can represent local intramolecular vibrations, nonlocal collective nuclear motions in solids, slow solvent motions, etc., depending on the situation being considered. The linear system-bath interaction is given by

    where $ \hat{f}_\mu $ is a system operator that characterizes the form and strength of the interaction between the system and the $ \mu $-th mode.

    For many important issues in the field of open quantum systems, it often happens that what is under primary consideration is the information of the system dynamics but not the motion of the bath degrees of freedom. Hence, it is sufficient to focus on the reduced density matrix of the system defined as $ \hat{\rho}_ {\rm{S}}(t) $ = Tr$ _ {\rm{B}}\{ {\rm{e}}^{-i\hat{H}t}\hat{\rho}(0) {\rm{e}}^{i\hat{H}t}\} $, where Tr$ _ {\rm{B}}\{\cdot\} $ represents to take partial trace over the bath degrees of freedom. Assuming that the initial total density matrix has a factorized form, $ \hat{\rho}(0) $ = $ \hat{\rho}_ {\rm{S}}(0) $$ \otimes $$ \hat{\rho}_ {\rm{B}}(0) $, where $ \hat{\rho}_ {\rm{B}}(0) $ = e$ ^{-\beta\hat{H}_ {\rm{B}}} $/Tr$ \{ {\rm{e}}^{-\beta\hat{H}_ {\rm{B}}}\} $ and $ \beta $ is the inverse temperature, one valuable advantage of the Hamiltonian Eq.(1) is that the formal expression of $ \hat{\rho}_ {\rm{S}}(t) $ can be elegantly expressed via the well-known Feynman-Vernon influence functional [1],

    where $ \hat{U}_0(t) $$ \equiv $e$ ^{-i\hat{H}_ {\rm{S}}t} $ is the bare system propagator, $ \mathcal{T}_+ $ denotes the chronological time-ordering operation, and $ \mathcal{F}(t) $ is the influence functional in the interaction picture, the expression of which is given by

    with [35]

    Here, $ n_\mu $ = $ 1/( {\rm{e}}^{\beta\omega_\mu}-1) $ is the thermal average occupation number of the $ \mu $-th phonon. $ \hat{A}(\tau) $ = $ \hat{U}_0(\tau)\hat{A}\hat{U}_0^\dagger(\tau) $ represents the operator $ \hat{A} $ at time $ \tau $ in the interaction picture. $ \hat{A}^{(+)} $ and $ \hat{A}^{(-)} $ are superoperators of $ \hat{A} $ that obey the rules $ \hat{A}^{(+)}\hat{O} $$ \equiv $$ \hat{A}\hat{O} $ and $ \hat{A}^{(-)}\hat{O} $$ \equiv $$ \hat{O}\hat{A} $, where $ \hat{O} $ is an arbitrary operator. $ + $ and $ - $ in the superoperators correspond to forward and backward paths, respectively. Note that the multiplication order of the operators in Eq.(6) and in the following derivations does not matter since we have utilized the benefit from the chronological time-ordering operation introduced in Eq.(4).

    It is easy to see that $ \hat{f}_\mu $ and $ \hat{f}_\mu^\dagger $ are not symmetric in $ \Phi_\mu(\tau, \tau') $. This is because they play different roles in the system-bath interaction. $ \hat{f}_\mu $ corresponds to the process that one phonon with index $ \mu $ is created as the system-bath scattering occurs, whereas $ \hat{f}_\mu^\dagger $ corresponds to the inverse process. As such, the former process dominates at low temperatures, and the latter one gradually takes effects as the temperature increases. For the further derivation, it is helpful to divide $ \Phi_\mu(\tau, \tau') $ into symmetric and antisymmetric parts in terms of $ \hat{f}_\mu $ and $ \hat{f}_\mu^\dagger $. Note that $ \Phi_\mu(\hat{f}_\mu, \hat{f}_\mu^\dagger) $$ \equiv $$ \Phi_\mu(\tau, \tau') $, the symmetric and antisymmetric parts are given by $ \Phi_\mu^{\mathrm{sym}}(\tau, \tau') $ = $ [\Phi_\mu(\hat{f}_\mu, \hat{f}_\mu^\dagger) $+$ \Phi_\mu(\hat{f}_\mu^\dagger, \hat{f}_\mu)]/2 $ and $ \Phi_\mu^{\mathrm{asym}}(\tau, \tau') $ = $ [\Phi_\mu(\hat{f}_\mu, \hat{f}_\mu^\dagger) $$ - $$ \Phi_\mu(\hat{f}_\mu^\dagger, \hat{f}_\mu)]/2 $, respectively. Therefore, we have

    and

    where

    which are the correlation functions. One reason for adopting this division scheme is that in many situations the asymmetric part is equal to zero and does not play a role in the dynamics. For instance, this occurs when $ \hat{f}_\mu $ are Hermitian operators or when Eq.(1) possesses translational symmetry. As such, one can just omit the antisymmetric part and simplify the problem to a large extent. In what follows, we will retain the antisymmetric part to make sure that the theory being derived is as general as possible.

    Although Eq.(4) offers a formally exact expression for $ \hat{\rho}_ {\rm{S}}(t) $, it cannot directly be used in calculations. Owing to the system-bath interaction, the operators of the system at different times are highly correlated with each other through the correlation functions Eq.(9). This leads to the time-nonlocality of the influence functional and the entanglement between forward and backward paths, which impedes a simple numerical solution to the dynamics. In the standard derivation of HSSEs, two powerful techniques, namely the stochastic unraveling [20-23], and the hierarchical expansion [3-6], are cleverly combined to deal with this issue.

  • The main idea of stochastic unraveling is to make use of the fact that the correlation between system operators at different times is of Gaussian type (i.e. only the second-order cummulant is nonzero) so that one can introduce correlated Gaussian stochastic fields via the Hubbard-Stratonovich transformation [36, Stratonovich–1957-p416] to decouple this correlation. It is known that a stochastic unraveling on the entire influence functional will completely remove this correlation and result in a simple equation of motion called the stochastic Liouville-von Neumann equation [19-25]. However, this equation may suffer severe stochastic convergence problems at long times or when the system-bath interaction is not weak. Instead of full stochastic unraveling, one often applies partial stochastic unraveling on the primary part of the influence functional and manages the residual part by other methods. The partition scheme plays a decisive role in the stochastic convergence properties of the resulting HSSE. Since the influence functional is completely characterized by the correlation functions Eq.(9), it is convenient to propose the partition scheme by dividing the correlation functions into a primary and a residue part, $ \alpha_\mu(t) $ = $ \alpha_\mu^{\mathrm{pri}}(t)+\alpha_\mu^{\mathrm{res}}(t) $ and $ \widetilde{\alpha}_\mu(t) $ = $ \widetilde{\alpha}_\mu^{\mathrm{pri}}(t) $+$ \widetilde{\alpha}_\mu^{\mathrm{res}}(t) $, where

    Here, we have introduced two parameters, $ \gamma_\mu $ and $ \widetilde{\gamma}_\mu $, to determine the partition scheme, and we requires that 0$ \le $$ \gamma_\mu, \widetilde{\gamma}_\mu $$ \le $$ n_\mu $+1/2. The other four parameters are $ \bar{n}_{\mu, 1} $ = $ (n_\mu $$ - $$ \gamma_\mu $)/2, $ \bar{n}_{\mu, 2} $ = ($ n_\mu $+1$ - $$ \gamma_\mu)/2 $, $ \bar{n}_{\mu, 3} $ = $ (n_\mu $$ - $$ \widetilde{\gamma}_\mu)/2 $, and $ \bar{n}_{\mu, 4} $ = $ (n_\mu $+1$ - $$ \widetilde{\gamma}_\mu)/2 $. It is obvious that the sum of primary and residue parts indeed recovers the original correlation functions. Then, based on Eq.(10), we separate $ \Phi_\mu(\tau, \tau') $ as $ \Phi_\mu(\tau, \tau') $ = $ \Phi_{\mu, \mathrm{pri}}(\tau, \tau') $+$ \Phi_{\mu, \mathrm{res}, +}(\tau, \tau') $+$ \Phi_{\mu, \mathrm{res}, -}^\dagger(\tau, \tau') $, where

    and

    It is worth noting that all of the cross terms between forward ($ + $) and backward ($ - $) paths are incorporated into $ \Phi_\mu^{\mathrm{pri}}(\tau, \tau') $. Therefore, the entanglement between the two paths can be removed if we apply the partial stochastic unraveling to the corresponding part of the influence functional. To this end, we introduce two sets of Gaussian stochastic fields, $ \{\xi_\mu^{\pm}(t)\} $ and $ \{\chi_\mu^{\pm}(t)\} $, which satisfy $ \langle\xi_\mu^{\pm}(t)\rangle $ = $ \langle\chi_\mu^{\pm}(t)\rangle $ = 0 and

    It is easy to see that $ \{\xi_\mu^{\pm}(t)\} $ and $ \{\chi_\mu^{\pm}(t)\} $ will be independent of each other if $ \Phi_\mu^{\mathrm{asym}}(\tau, \tau') $ is zero. With the aid of the stochastic fields, the influence functional can now be reformulated as

    where $ \mathcal{M}\{\cdot\} $ denotes to take the average over the stochastic fields, and

    One can find that the time-nonlocality of the influence functional is partially removed. To go further, assuming that the system is in a pure state at the beginning, $ \hat{\rho}_ {\rm{S}}(0) $ = $ |\psi_0\rangle\langle\psi_0| $, and substituting Eqs. (14)$ - $(15) into Eq.(4), one arrives at

    and the stochastic propagator is given by

    where we have explicitly written out the residual term with $ \hat{A}_\mu(t) $ = $ \int_0^t\mathrm{d}\tau {\rm{e}}^{i\omega_\mu(t-\tau)}\hat{f}_\mu(\tau) $ and $ \hat{B}_\mu(t) $ = $ \int_0^t\mathrm{d}\tau {\rm{e}}^{-i\omega_\mu(t-\tau)}\hat{f}_\mu(\tau) $.

    So far, via the introduction of stochastic fields, we have decomposed the reduced density matrix into the product of forward and backward wavefunctions. Nevertheless, the propagator of the wavefunctions is still not manageable due to the existence of the residual influence functional. This difficulty can be overcome by the hierarchical expansion technique [3-6].

  • Because of the double integral in the stochastic propagator, it is impossible to derive a self-consistent equation of motion of $ |\psi_\pm(t)\rangle $ itself without any approximation. Therefore, we must expand the forward and backward wavefunctions to a complete set which consists of infinite number of auxiliary wavefunctions so that the joint equation of motion of all the wavefunctions is closed. Concretely speaking, we define

    where $ k_\mu $, $ l_\mu $, $ m_\mu $ and $ n_\mu $ are non-negative integer indexes that characterize the auxiliary wavefunction, $ {\bf{k}} $ = $ (\cdots, k_\mu, \cdots) $ is a short-hand vector notation of $ \{k_\mu\} $, and so do the other three vectors. Note that the original forward and backward wavefunctions correspond to the one with all the indexes being zero, $ |\psi_\pm(t)\rangle $ = $ |\Psi_\pm^{[{\bf{0}}, {\bf{0}};{\bf{0}}, {\bf{0}}]}(t)\rangle $. Taking the time derivative of Eq.(19), one readily arrives at the final result, which is a closed set of equations of motion with a hierarchical structure,

    where $ {\bf{k}}_\mu^\pm = (\cdots, k_\mu\pm1, \cdots) $, and so do the other three vectors. Eq.(20) is the HSSE in a general form. Setting $ t $ = 0 in Eq.(19), one readily obtains the initial conditions for Eq.(20),

    Theoretically speaking, an infinite number of auxiliary wavefunctions are coupled with each other through Eq.(20), and one has to solve all of them simultaneously, which is obviously an impossible task. In practical implementation, Eq.(20) is truncated to a certain hierarchical level $ L_\mathrm{max} $, which is a non-negative integer. The simplest way of truncation is to maintain auxiliary wavefunctions with $ \sum\limits_\mu(k_\mu+l_\mu+m_\mu+n_\mu) $$ \le $$ L_\mathrm{max} $, and omit all of the rest ones in higher levels. $ L_\mathrm{max} $ is required to be large enough to obtain a convergent result. One can also use a more flexible and efficient truncation scheme like the $ n $-particle or $ n $-mode schemes [38].

    The numerical performance of HSSEs is reflected in two aspects. One is the stochastic convergence property of the stochastic fields, the other is the minimal hierarchical level that is required to obtain an accurate result. It should be noted that Eq.(20) has yet to be completely specified because all of the parameters $ \{\gamma_\mu, \widetilde{\gamma}_\mu\} $ introduced in Eq.(10) have not been assigned so far. As a result, Eq.(20) actually represents an infinite number of HSSEs, and only a small part of them bear satisfactory numerical performance. Furthermore, in many situations, it is possible that the stochastic fields obeying Eq.(13) can be generated under the assumptions that $ \xi_\mu^{+}(t) $ = $ \xi_\mu^{-}(t) $ and $ \chi_\mu^{+}(t) $ = $ \chi_\mu^{-}(t) $. As such, the forward and backward paths exactly coincide with each other, and we only need to propagate one of them in calculations. We point out that the HSSEs suggested in the literatures usually consider the situation where the antisymmetric part of the influence functional do not take part in the dynamics, and so do the parameters $ \{\widetilde{\gamma}_\mu\} $. More specifically, $ \gamma_\mu $ = 0, $ \gamma_\mu $ = $ \sqrt{\bar{n}_\mu(\bar{n}_\mu+1)} $, $ \gamma_\mu $ = $ \bar{n}_\mu $+1/2, and $ \gamma_\mu $ = $ \bar{n}_\mu $ have been suggested in Refs.[30-33], respectively. In addition, the HSSE proposed in Ref.[32] propagates both forward and backward paths, as opposed to other ones which only propagate a single path.

    Stochastic fields that obey Eq.(13) can be generated by various generation schemes that have been suggested ever. The most universal one is the discrete Fourier transformation scheme [27, 28, 39, 40], where white noises are translated into colored ones by a filtering kernel originating from the Cholesky decomposition on the discrete correlation function. This scheme is generally applicable for Gaussian stochastic fields regardless of their concrete correlation properties. Another popular scheme [28, 32, 41, 42] is to make a clever ansatz for the stochastic fields so that they can be expressed by simple analytical formulas in terms of trivial Gaussian noises. Past experience shows that ansatz methods usually provide stochastic fields with better convergence properties, but they are only suitable for specific forms of correlation properties. One can also flexibly separate the stochastic fields into different uncorrelated parts and generate each part via different schemes. It would be preferable that the primary and residual parts are generated by the ansatz and the discrete Fourier transformation schemes, respectively, so as to produce well-behaved stochastic fields for arbitrary correlation properties.

    Another factor of practical aspect is that Eq.(20) is a linear stochastic differential equation that usually does not conserve the norm of the wavefunction in each random trajectory, which is unfavorable for a good stochastic convergence property. One way to circumvent this problem is to resort to a good generation scheme for the stochastic fields, and the other way is to adopt a nonlinear counterpart of Eq.(20). Various numerical calculations [30, 33, 34] have proven that a non-linear version of HOPS can effectively preserve the norm and significantly improve the stochastic convergence. The effectiveness of nonlinearization for other HSSEs remains further demonstrations.

    In deriving Eq.(20), we do not make any approximations, neither do we assume a special property for $ \hat{f}_\mu $ or for the total Hamiltonian Eq.(1). Therefore, Eq.(20) is suitable for the investigation of a large variety of issues in the field of open quantum systems, and its form can often be further simplified to some extent, depending on the specific situation being considered. In Section 2.3, we will discuss the concrete form of Eq.(20) in different situations.

    It is worth mentioning that another strategy to deal with the residual influence functional is to perform perturbation expansion in the system-bath interaction strength and truncate at the second order, which leads to several approximate but highly efficient perturbative stochastic Schrödinger equations (PSSE) [41-43]. Due to the relatively cheap numerical costs, PSSEs have been widely used in the studies of singlet fission, exciton relaxation and dissociation, carrier transport, absorption spectra of organic aggregates, etc. For details, the readers can refer to Ref.[43].

    Finally, we point out that Eq.(20) is derived on the basis of a bosonic environment. Nevertheless, the key idea may also apply to an open quantum system coupled with a fermionic environment, for instance, the quantum impurity system. In doing so, the most difficult part is the numerical realization of the fermionic bath-induced stochastic fields, which are no longer c-number noises but Grassmann fields. Only very recently, the first numerically feasible scheme has been proposed for the realization of stochastic Grassmann fields in an approximate manner [44-46]. Extending the HSSE methodology to the fermionic case may be a promising way to arrive at an exact stochastic description of fermionic open systems.

  • When the interaction operators are Hermitian, $ \hat{f}_\mu $ = $ \hat{f}_\mu^\dagger $, we have $ \hat{A}_\mu(t) $ = $ \hat{B}_\mu^\dagger(t) $. In this case, half of the indexes in Eq.(19) are redundant, and we can recast the expression of the auxiliary wavefunctions into

    where $ \bf{r} $ = $ [\cdots, r_\mu, \cdots] $ and $ \bf{s} $ = $ [\cdots, s_\mu, \cdots] $ with $ r_\mu $ = $ k_\mu+n_\mu $ and $ s_\mu $ = $ m_\mu $+$ l_\mu $. Furthermore, it is easy to prove that

    and

    Substituting Eqs. (22)$ - $(24) into Eq.(20), we obtain a simplified version of HSSE,

    It is readily seen that Eq.(25) depends only on $ \gamma_\mu $ but not on $ \widetilde{\gamma}_\mu $, and the statistical properties of the stochastic fields degenerate to

    because there is only one set of the fields in Eq.(25). This is not surprising, since in this situation the antisymmetric part of the influence functional is equal to unity and does not have an impact on the dynamics.

  • For ultrafast quantum dynamics processes occurring in solids, one often considers a system with translational symmetry and the periodic boundary condition. Typical models of this kind are the Su-Schrieffer-Heeger model [47], the Holstein-Peierls model [48, Peierls–1996-p, 50], the Fröhlich polaron model [51], etc. As such, the system and the bath are both well-described by the crystal momentum, and their interaction conserves the total crystal momentum. Therefore, we can rewrite Eqs. (2)$ - $(3) as

    and

    Here, the index $ \mu $ has been replaced by the crystal momentum defined as $ q $ = $ 2\pi m/N $, where $ N $ is the total number of $ k $ points in the first Brillouin zone, $ m $$ = 0, 1, \ \cdots, \ $$ \pm \frac{N-1}{2} $ for an odd $ N $ and $ m $ = $ - \frac{N}{2}+1, -\frac{N}{2}+2, \ \cdots, \frac{N}{2} $ for an even $ N $. Furthermore, due to the symmetry of the Brillouin zone and the Hermiticity of the Hamiltonian, we have $ \hat{f}_q^\dagger $ = $ \hat{f}_{-q} $ and $ \omega_q $ = $ \omega_{-q} $, which lead to $ \hat{A}_q(t) $ = $ \hat{B}_{-q}^\dagger(t) $ and $ \hat{B}_q(t) $ = $ \hat{A}_{-q}^\dagger(t) $. Therefore, again, half of the indexes in Eq.(19) are redundant, and we can re-express Eq.(19) as

    where $ \bf{u} $ = $ [\cdots, u_q, \cdots] $ and $ \bf{v} $ = $ [\cdots, v_q, \cdots] $ with $ u_q $ = $ k_q $+$ n_{-q} $ and $ v_q $ = $ m_q $+$ l_{-q} $. Likewise, one can directly prove that

    and

    To impose the symmetry of the Hamiltonian to Eq.(20), we further set $ \gamma_q $ = $ \gamma_{-q} $, which leads to $ \bar{n}_{q, 1} $ = $ \bar{n}_{-q, 1} $ and $ \bar{n}_{q, 2} $ = $ \bar{n}_{-q, 2} $. In the end, with all of the above properties, we can transform Eq.(20) into a simpler form. This is accomplished by invoking the following relationship,

    where $ Y $ is a symmetric function of $ q $. For example, using this relationship, it is easy to show that

    Another example is that

    where $ \widetilde{\xi}_q^\pm(t) $ = $ \left(\xi_q^\pm(t) + \xi_{-q}^\pm(t)\right))/2 $ and $ \widetilde{\chi}_q^\pm(t) $ = $ \left(\chi_q^\pm(t)-\chi_{-q}^\pm(t)\right)/2 $ are two newly-defined stochastic fields. Note that $ \widetilde{\chi}_0^\pm(t) = \widetilde{\chi}_\pi^\pm(t) $ = 0. Continuing this procedure, we finally obtain

    The statistical properties of the newly-defined stochastic fields are

    Again, we can see that Eq.(35) is irrelevant to $ \widetilde{\gamma}_q $ and $ \widetilde{\alpha}_q(t) $ because the antisymmetric part of the influence functional reduces to unity in this situation.

  • Without any assumptions of the Hamiltonian, the hierarchical structure of Eq.(20) can also be naturally reduced if we properly choose the values of $ \gamma_\mu $ and $ \widetilde{\gamma}_\mu $. Setting $ \gamma_\mu $ = $ \widetilde{\gamma}_\mu $ = $ n_\mu $, which corresponds to $ \bar{n}_{\mu, 1} $ = $ \bar{n}_{\mu, 3} $ = 0 and $ \bar{n}_{\mu, 2} $ = $ \bar{n}_{\mu, 4} $ = 1/2, it is not hard to find that the zeroth order auxiliary wavefunction is completely decoupled from higher-order auxiliary wavefunctions with nonzero $ {\bf{k}} $, $ {\bf{l}} $ or $ {\bf{n}} $. Therefore, it is no longer necessary to solve all of the wavefunctions simultaneously. In fact, we only need to solve those with $ {\bf{k}} $ = $ {\bf{l}} $ = $ {\bf{n}} $ = $ {\bf{0}} $. Defining

    and carefully inspecting Eq.(20), one obtains

    In the cases that $ \hat{f}_\mu $ is Hermitian or the Hamiltonian has translational symmetry, Eq.(38) can be further simplified. Therefore, Eq.(38) has the great advantage that the total number of auxiliary wavefunctions are tremendously reduced so that it can be applied to larger systems as compared with other HSSEs. For practical purposes, careful numerical tests are still required for the stochastic convergence property of Eq.(38).

  • In the condensed phase, the environment such as the solvent or the protein scaffold often exhibits local and overdamped features. Consequently, the interaction operators $ \{\hat{f}_\mu\} $ are very local in space, and each of them couples to their own local environment. Therefore, it is valid to invoke the independent bath approximation and assume that each local environment consists of a continuum of harmonic oscillators. In the end, Eq.(2) and Eq.(3) are reformulated as $ \hat{H}_ {\rm{B}} $ = $ \sum\limits_\mu\sum\limits_j\omega_{\mu j}\left(\hat{b}_{\mu j}^\dagger\hat{b}_{\mu j}+\frac{1}{2}\right) $ and $ \hat{H}_ {\rm{I}} $ = $ \sum\limits_\mu\hat{f}_\mu\otimes\sum\limits_{j}c_{\mu j}(\hat{b}_{\mu j}^\dagger+\hat{b}_{\mu j}) $, respectively. Here, the spectrum of each bath are characterized by the frequency $ \omega_{\mu j} $ and the system-bath interaction strength constant $ c_{\mu j} $, which can be elegantly described by the spectral density function $ J_{\mu} $ = $ \pi\sum\limits_{j}c_{\mu j}^2\delta(\omega-\omega_{\mu j}) $. Under this situation, the bath correlation function can often be written as the sum of exponential functions, $ \alpha_\mu(t) $ = $ \sum\limits_{m}d_{\mu m} {\rm{e}}^{-\Omega_{\mu m}t} $. Here, $ d_{\mu m} $ and $ \Omega_{\mu m} $ represent no longer the properties of a single oscillator but the overall properties of a substantial part of the bath. In general, $ d_{\mu m} $ and $ \Omega_{\mu m} $ can be either real or complex numbers. For instance, one of frequently adopted models in the condensed phase is the Debye-Drude spectral density function [52]

    where $ \lambda $ and $ \omega_c $ denote the reorganization energy and the bath cutoff frequency, respectively. For this model, the corresponding correlation function reads

    where $ \omega_\nu $ = $ 2\pi\nu/\beta $ is the Matsubara frequency. We can still use the same techniques introduced previously to deal with this model. Concretely speaking, adopting Eq.(39) for all $ \hat{f}_\mu $ with the same $ \lambda $ and $ \omega_c $ and defining the residual correlation function as $ \alpha_\mathrm{res}(t) $ = $ -i\lambda\omega_c {\rm{e}}^{-\omega_ct} $, we straightforwardly obtain the following HSSE [32]:

    and the stochastic fields obey

    A good generation scheme for such stochastic fields has been proposed in Ref.[32, 42]. One great merit of Eq.(41) is that all of the temperature effects are included in the stochastic fields so that the equations possess the simplest hierarchical structure. The truncation level of Eq.(41) is usually less than that of HEOM, since the primary part of the system-bath interaction has already been incorporated into the stochastic fields. Furthermore, various numerical tests [32, 43, 53] have also demonstrated that Eq.(41) has a satisfied stochastic convergence property over a large parameter regimes if the generation scheme for stochastic fields is properly chosen. In the next section, we will present the applications of Eq.(41) in the EET processes in large-scale systems to demonstrate its power for the efficient dynamics simulations of open quantum systems.

  • As compared with other numerically exact methods, HSSEs are more suitable for the dynamics simulations in relatively large systems owing to the preferable scaling relationship between the numerical costs of HSSEs and the system size. So far, HSSEs have been used to investigate EET in the FMO trimer [32], resonant vibration-assisted EET [54], the absorption and circular dichroism spectra of light-harvesting complex Ⅱ [55], the non-Markovianity of a system coupled with an ultraslow bath [53], and the calculation of two-dimensional electronic spectroscopy [56, 57]. In this section, we review two recent applications of HSSEs, that is, the EET process in the FMO complex [32, 43] and the resonant vibration-assisted EET mechanism [54], to demonstrate the applicability and efficiency of HSSEs in simulating relatively large-scale systems.

  • EET in the photosynthetic process are known to have an incredible quantum efficiency of nearly unity [Mohseni–2014-p, Blankenship–2014-p], which has aroused great research interest during the decades. One prototype for the studies of such a process is the FMO [60, 61], a homotrimer found in green sulfur bacteria and funnelling excitation energy from light-harvesting chlorosome toward the reaction center. Over ten years ago, Engel et al. [62] carried out elaborate two-dimensional electronic spectroscopy (2DES) experiments on the FMO complex at cryogenic temperature and found unexpected long-lasting quantum beating phenomenon. This rapidly fuels intense debates about the role played by quantum coherence in EET [58, 63-71].

    From the theoretical perspective, it is a key point to give theoretical descriptions with enough accuracy so as to ravel out the ambiguity of experimental results. Using HEOM under the high-temperature approximation, Ishizaki and Fleming [63] performed the first non-perturbative and non-Markovian quantum dynamics calculations of EET in the FMO subunit. Their work is based on a primal model of apo-FMO that consists of seven bacteriocholoriphyll (BChl). They assumed that each BChl site is coupled to an independent thermal bath. These baths share the same Debye-Drude spectral density function, Eq.(39), which is well-justified and parameterized experimentally. It is adequate to consider EET within one subunit because the inter-subunit couplings are very weak in this primal model. Nevertheless, it has been confirmed recently [72, 73] that in each FMO subunit there exists an eighth BChl (holo-FMO), which is suggested to link two neighbouring subunits [74, 75]. To ascertain the function of the eighth BChl in EET, it is necessary to consider EET in a full 24-site FMO trimer, which is much more demanding for theoretical simulations. In the following, we demonstrate the applicability and efficiency of HSSE, Eq.(41), by calculating EET in both 7-site FMO subunit and 24-site FMO trimer and comparing the results with those from the numerically exact HEOM.

    Population evolutions of the first four BChls in a 7-site FMO subunit calculated by HEOM and HSSE are both shown in FIG. 1, where the same model used by Ishizaki and Fleming is adopted, and four different sets of $ T $ and $ \omega_c $ with a fixed bath reorganization energy of $ \lambda $ = 35 cm$ ^{-1} $ are considered. The local excitation at BChl 1 is set as the initial condition. Not surprisingly, the results of HEOM and HSSE coincide with each other since they are both numerically exact methods. It can be seen that the truncation level of hierarchy for HEOM, $ L_\mathrm{HEOM} $, increases as the temperature increases or $ \omega_c $ decreases, while the number of Matsubara frequencies $ K $ increases in the opposite situation. By contrast, the truncation levels of hierarchy for HSSE, $ L_\mathrm{HSSE} $, are nearly independent of the two parameters and are all smaller than those for HEOM. In addition, as mentioned previously, since all of the temperature effects are absorbed in the stochastic fields, Eq.(41) does not require Matsubara or Padé frequency decomposition schemes. As a consequence, the memory space needed for solving HSSE is far smaller than that for solving HEOM. The stochastic nature of HSSE also makes it straightforward to include static disorder in the system Hamiltonian without extra numerical efforts.

    Figure 1.  Population evolution of the first four BChl sites for a 7-site FMO complex. The reorganization energy $ \lambda $ is set as 35 cm$ ^{-1} $. The truncation levels of hierarchy, $ L_\mathrm{HEOM} $, and the number of Matsubara frequencies, $ K $, that are needed in HEOM are (a) (4, 3), (b) (5, 1), (c) (6, 2), (d) (10, 1). The truncation levels of hierarchy for HSSE, $ L_\mathrm{HSSE} $, are 3 in all four cases, while the random trajectory numbers are (a) 2$ \times $10$ ^4 $, (b) $ 10^4 $, (c) $ 10^4 $, (d) 5$ \times $10$ ^3 $. Reprinted with permission from Ref.[43] ©2018 Wiley Periodicals, Inc.

    The advantages of HSSE are more evident for a larger system. FIG. 2 presents the results of EET in a full 24-site FMO trimer calculated by HSSE. The intra-subunit excitonic Hamiltonian and the inter-subunit couplings are adopted from Refs. [76] and [75], respectively (for the details, see Ref.[32]). It can be seen that $ L_\mathrm{HSSE} $ is unchanged as compared with that of a 7-site result, and the random trajectory numbers slightly decrease. The calculations of HEOM for the population evolutions in a FMO trimer at the bath conditions of $ \omega_c^{-1} $ = 50 fs and $ T $ = 300 K have been done by Wilkins and Dattani [70] with the aid of a more efficient propagation scheme. Table Ⅰ lists the ratio of computational costs in a 24-site simulation to those in a 7-site simulation [43], where for HEOM in the 24-site case we have assumed that $ L_\mathrm{HEOM} $ and $ K $ are the same as the 7-site case. One can see that the computational costs of both HEOM and HSSE increase rapidly, but the increasing rate of the former one is far greater than that of the latter one, especially in the case $ \omega_c^{-1} $ = 50 fs and $ T $ = 300 K. FIG. 2 also presents the approximate results calculated by PSSE. The remarkable accuracy of PSSE in the FMO complex is well justified.

    Figure 2.  Population evolution of the most populated eight BChls for a full 24-site FMO complex. The reorganization energy $ \lambda $ is set as $ 35 $ cm$ ^{-1} $. $ L_\mathrm{HSSE} $ are 3 in all four cases, while the random trajectory numbers of HSSE are (a) $ 10^4 $, (b)5$ \times $10$ ^3 $, (c) $ 10^4 $, (d) 5$ \times $10$ ^3 $. Also shown are the approximate results calculated by PSSE. Reprinted with permission from Ref.[43] ©2018 Wiley Periodicals, Inc.

    Table Ⅰ.  The ratio of computational costs in a 24-site simulation to those in a 7-site simulation.

  • Recently, there exists a conjecture [67, 77-83] that the long-lived beating signals in the cross-peak of 2DES arising in the FMO complex [62, 64] and other systems [84-86] originate from the interaction between the exciton and an underdamped vibration. Among the relevant studies, a strand of researches focus on the phenomenon of resonant vibration-assisted EET, which is a EET process accelerated by an underdamped vibration as the vibrational frequency is resonant with the excitonic energy gap. Theoretically, various efforts have been devoted to confirming the existence of such a phenomenon [87-101]. It has been found [88, 102, 103] that the resonant frequency is relevant to the total reorganization energy of EET. Synergistic effects of multiple modes have also been reported [102].

    Nevertheless, several studies have found that resonant underdamped vibrations have little effect on the EET processes in the FMO complex [104, 105] and light-harvesting complex Ⅱ [18]. Researches on the phycobilliproteins and photosystem Ⅱ reaction center [106-108] have also pointed out that underdamped vibrations may even hamper the EET and charge separation processes. In addition, most of the theoretical studies have assumed a local exciton-phonon coupling, whereas in realistic systems nonlocal contributions (non-Condon effect) may also be substantial. To clarify the conflicts and give a definite condition required for the occurrence of resonant vibration-assisted EET, it is necessary to consider both local and nonlocal exciton-phonon couplings for a broad parameter regimes.

    We have used Eq.(41) to investigate the resonant vibration-assisted EET in a heterodimer system and assumed that the donor and the acceptor are respectively coupled to two independent continuous baths which share the same Debye-Drude spectral density function, Eq.(39). In this situation, the two independent baths can be exactly mapped into one anti-correlated bath, and this trick can be used to reduce the computational cost [54]. The Hamiltonian of the excitonic states coupled with an underdamped mode is given by

    Here, $ |D\rangle $ or $ |A\rangle $ represent a local excitation at the donor or the acceptor, respectively. $ \Delta G $ is the site-energy difference between $ |D\rangle $ and $ |A\rangle $. $ J_\mathrm{DA} $ is the excitonic coupling. $ \hat{P} $ and $ \hat{Q} $ are the mass-weighted momentum and coordinate operators of the mode with frequency $ \Omega $. $ \kappa $ determines the coupling strength between excitonic states and the mode. The variable $ \phi $, which we refer to as the coupling angle hereafter, is introduced to specify the coupling form of the underdamped mode to the exciton. In general, $ \phi $ could vary from zero to $ \pi $. In the special case $ \phi $ = 0 or $ \phi $ = $ \pi $, the mode only diagonally coupled to the excitonic states, which corresponds to a pure local exciton-phonon coupling. On the other hand, $ \phi $ = $ \pi/2 $ corresponds to a pure nonlocal exciton-phonon coupling. The last term in Eq.(43) accounts for the reorganization energy from the local component of the exciton-phonon coupling. In all calculations, the temperature is 277 K, the bath cutoff frequency is $ \omega_c^{-1} $ = 50 cm$ ^{-1} $, and we fix $ \kappa^2/2\Omega^2 $ = 1 cm$ ^{-1} $ to ensure a very small coupling between excitonic states and the mode. The initial excitation is located at the donor.

    For the convenience of discussions, we further introduce three quantities. The first one is the mixing angle defined as $ \theta $ = $ \tan^{-1}(2J_\mathrm{DA}/\Delta G) $, which describes the mixing degree of $ |D\rangle $ and $ |A\rangle $ in the electronic eigenstates. The second one is the energy gap between two electronic eigenstates, calculated as $ \Delta E $ = $ \sqrt{\Delta G^2+4J_\mathrm{DA}^2} $. The last one is the population ratio $ R $ that measures the promotion effect of the underdamped mode on EET. The population ratio is defined as the ratio of average transferred population in the presence of the mode to that without the presence of the mode. The larger the population ratio, the more significant the promotion effect. The readers can refer to Ref.[54] for a detailed expression of $ R $.

    We first investigate the resonant vibration-assisted EET in a downhill energy landscape with weak thermal fluctuations. FIG. 3 presents the population ratios $ R $ for different combinations of $ \Omega $ and $ \phi $ at $ \Delta G $ = 250 cm$ ^{-1} $ and $ \lambda $ = 5 cm$ ^{-1} $. The contour maps exhibit colorful features. It is immediately seen that a striking peak clearly appears as the vibrational frequency $ \Omega $ is resonant with $ \Delta E $. One can also find that the promotion effect is strongly dependent on the coupling angle $ \phi $. For instance, at the resonant frequency, the population ratio reaches its maximum at $ \phi $ = $ \theta $+$ \pi/2 $, while it rapidly declines to unity (i.e., no promotion effect) as $ \phi $ approaches $ \theta $. The results demonstrate that it is not enough to only concentrate on the vibrational frequency in the resonant vibration-assisted EET.

    Figure 3.  Contour maps of the population ratio $ R $ that measures the magnitude of vibration-assisted promotion in EET. The value of $ J_\mathrm{DA} $ is given in each subgraph, and other parameters are $ \Delta G $ = 250 cm$ ^{-1} $, $ \lambda $ = 5 cm$ ^{-1} $. The black and white thick dashed lines mark $ \Omega $ = $ \Delta E $ and $ \phi $ = $ \theta $, respectively. All of the areas with $ R $ smaller than 1 are painted black and looped by thin white dashed curves. Adapted with permission from Ref.[54] ©2018 Wiley Periodicals, Inc.

    These results can be explained by inspecting the downhill flow of the population between excitonic eigenstates, which form a natural basis set for the excitonic subspace at weak thermal fluctuations [109]. Transforming the excitonic degree of freedom into the eigenstate representation, Eq.(43) can be rewritten as

    where $ |E_1\rangle $ = cos$ \frac{\theta}{2}|D\rangle $+sin$ \frac{\theta}{2}|A\rangle $ and $ |E_2\rangle $ = $ - $sin$ \frac{\theta}{2}|D\rangle $+ cos$ \frac{\theta}{2}|A\rangle $ are high- and low-energy eigenstates, respectively. In the case that $ \Delta G $$ \gg $$ J_\mathrm{DA} $, the mixing angle $ \theta $ is quite small, and consequently $ |D\rangle $ and $ |A\rangle $ mostly comprise $ |E_1\rangle $ and $ |E_2\rangle $, respectively. As such, the population transfer between $ |E_1\rangle $ and $ |E_2\rangle $ is very in line with that between $ |D\rangle $ and $ |A\rangle $. Inspecting Eq.(44), one can find that the mode-generated $ |E_1\rangle $$ \rightarrow $$ |E_2\rangle $ channel is controlled by the factor $ \sin(\phi-\theta) $. In the case $ \phi $ = $ \theta $, the channel is completely shut off, and there is no promotion effect no matter what the value of $ \Omega $ is. As $ \phi $ gradually deviates from $ \theta $, transition channels $ |E_1, \nu\rangle $$ \rightarrow $$ |E_2, \nu+1\rangle $ are turned on, where $ |\psi, \nu\rangle $ denotes that the excitonic state is $ |\psi\rangle $ and the underdamped mode is in the $ \nu $-th energy level. Hence, the EET process is accelerated when $ \Omega $ is resonant with $ \Delta E $. Finally, at $ \phi $ = $ \theta $+$ \pi/2 $, the coupling between $ |E_1, \nu\rangle $ and $ |E_2, \nu+1\rangle $ reaches its maximum, which leads to a most favorable $ \phi $ for EET. In addition, multiphonon processes that occur at $ \Omega $$ \approx $$ \Delta E/n $ ($ n = 2, 3, \cdots $) are also allowed, which are confirmed by the appearance of secondary maximum at $ \Omega $ = $ \Delta E/2 $ in FIG. 3.

    Nevertheless, the above picture gradually loses its validity as thermal fluctuations become more and more strong. It is known that strong background phonon fluctuations will dress the exciton to behave like a quasi-particle called the polaron [110], and the excitonic coupling $ J_\mathrm{DA} $ is reduced to $ \widetilde{J}_\mathrm{DA} $ = $ DJ_\mathrm{DA} $, where $ D $$ < $1 is the Debye-Waller factor. The stronger $ \lambda $ is, the smaller the effective excitonic coupling. As such, the effective mixing angle $ \widetilde{\theta} $ will gradually vary to zero as $ \lambda $ becomes very large, and the optimal coupling angle changes from $ \phi $ = $ \theta $+$ \pi/2 $ to $ \phi $ = $ \pi/2 $. FIG. 4 presents the results at very strong thermal fluctuations. Indeed, the optimal $ \phi $ always appears at $ \pi/2 $ no matter what the value of the mixing angle is. Furthermore, since the strong thermal fluctuations break down the energy-level description of the exciton, the promotion effect arises over a broad range of $ \Omega $ from tens to hundreds of wavenumbers.

    Figure 4.  Contour maps of the population ratio $ R $ that measures the magnitude of vibration-assisted promotion in EET. The value of $ J_\mathrm{DA} $ is given in each subgraph, and other parameters are $ \Delta G $ = 40 cm$ ^{-1} $, $ \lambda $ = 150 cm$ ^{-1} $. The black and white thick dashed lines mark $ \Omega $ = $ \Delta E $ and $ \phi $ = $ \theta $, respectively. All of the areas with $ R $ smaller than 1 are painted black and looped by thin white dashed curves. Reprinted with permission from Ref.[54] ©2018 Wiley Periodicals, Inc.

    Finally, it should be pointed out that in the case $ J_\mathrm{DA} $$ \gg $$ \Delta G $, the promotion effect completely disappears regardless of other parameters. This is numerically proved elsewhere [54] and can be inferred from both FIG. 3 and FIG. 4. In the two situations, the population ratio gradually approaches unity as $ J_\mathrm{DA} $ becomes larger and larger.

  • The accurate and efficient dynamics simulation of open quantum systems is a long-standing subject with a profound history. Here, we have reviewed a class of numerically exact approaches of high efficiency termed as HSSEs, which are wavefunction-based approaches that benefit from the linear scaling property of Hilbert space with respect to the system size. Starting from the influence-functional description of an open quantum system coupled to a bosonic bath, we have presented in great detail the general theoretical framework of HSSEs, discussed their concrete forms in different situations, and demonstrated their practicality by applying them to the simulations of EET processes in a large-scale photosynthetic complex and to the investigation of resonant vibration-assisted EET. Due to their unique merits, it is expected that HSSEs will continue to play an active role in the fields of open quantum systems in the future.

  • This work is supported by the National Natural Science Foundation of China (No.22033006, No.21833006, and No.21773191).

Reference (110)

Catalog

    /

    DownLoad:  Full-Size Img  PowerPoint
    Return
    Return